Velocity Jumps for Molecular Dynamics

Nicolaï Gouraud, Louis Lagardère, Olivier Adjoua, Thomas Plé, Pierre

Monmarché, Jean-Philip Piquemal

► To cite this version:

Nicolaï Gouraud, Louis Lagardère, Olivier Adjoua, Thomas Plé, Pierre Monmarché, et al. Velocity Jumps for Molecular Dynamics. Journal of Chemical Theory and Computation, 2025, 21 (6), pp.2854-2866. 10.1021/acs.jctc.5c00023 . hal-04852362

HAL Id: hal-04852362

https://hal.science/hal-04852362v1

Submitted on 10 Dec 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Velocity Jumps for Molecular Dynamics

Nicolaï Gouraud<sup>1,2,3</sup>, Louis Lagardère<sup>1,3</sup>, Olivier Adjoua<sup>1</sup>, Thomas Plé<sup>1</sup>, Pierre Monmarché<sup>1,2,4,*</sup>, and Jean-Philip Piquemal<sup>1,3,*</sup>

<sup>1</sup>Sorbonne Université, CNRS, LCT UMR 7616, 75005 Paris, France

<sup>2</sup>Sorbonne Université, Université Paris Cité, CNRS, INRIA, LJLL UMR 7598, 75005 Paris, France

<sup>3</sup>Qubit Pharmaceuticals, Advanced Research Department, 75014 Paris, France

<sup>4</sup>Institut Universitaire de France, 75005 Paris, France

*pierre.monmarche@sorbonne-universite.fr, jean-philip.piquemal@sorbonne-universite.fr

Abstract

We introduce the Velocity Jumps approach, denoted as JUMP, a new class of Molecular dynamics integrators, replacing the Langevin dynamics by a hybrid model combining a classical Langevin diffusion and a piecewise deterministic Markov process, where the expensive computation of long-range pairwise interactions is replaced by a resampling of the velocities at random times. This framework allows for an acceleration in the simulation speed while preserving sampling and dynamical properties such as the diffusion constant. It can also be integrated in classical multi-timestep methods, pushing further the computational speedup, while avoiding some of the resonance issues of the latter thanks to the random nature of jumps. The JUMP, JUMP-RESPA and JUMP-RESPA1 integrators have been implemented in the GPU-accelerated version of the Tinker-HP package and are shown to provide significantly enhanced performances compared to their BAOAB, BAOAB-RESPA and BAOAB-RESPA1 counterparts respectively.

# 1 Introduction

Molecular dynamics (MD) [1, 2] is a popular numerical tool to infer macroscopic properties of matter from simulations at a microscopic level, with applications ranging from material sciences to biology [3, 4, 5]. Because of the discrepancy between the microscopic and macroscopic time and space scales, in order to obtain relevant results, one needs to simulate complex molecular systems with a high precision and for a very long time [6, 7, 8, 9, 10, 11]. In recent years, this has been allowed by the development of empirical force fields (such as CHARMM [12], AMBER [13], OPLS [14], GROMOS [15], AMOEBA [16], SIBFA [17] and others [18, 19, 20]), as well as the rise of high-performance computing with GPUs [21, 22, 23, 24], allowing to run massively parallel codes. Nowadays, many molecular simulation softwares are available (LAMMPS [25], NAMD [26], GROMACS [27], AMBER [28], Tinker-HP [29] and others [9, 30, 31, 32]) and a lot of attention is drawn to improving the efficiency of the simulation algorithms.

In the framework of Langevin dynamics, multi-timestep methods [33, 34, 35, 36] are now a standard way to make faster simulations than classical integrators such as BAOAB [37, 38, 39, 40], but they are limited by resonance effects that bound the maximum usable time step at a given accuracy [41, 42]. Various alternative strategies have been proposed such as the Generalized Langevin Equation (GLE) [43] or the stochastic isokinetic extended phase-space algorithm [44, 45, 46]. They yield significant acceleration but rely either on some empirical fitting [43] or have an important negative impact on the dynamical properties such as the diffusion coefficient [44], which is an indication of a limitation in the effective sampling rate. In the last case, the computational gain is reduced by the fact that a longer trajectory is necessary to keep the same amount of sampling. To speed up molecular dynamics without resorting to some fitting and while not affecting too much the dynamics [43, 47], one possibility is to combine well-chosen integrators within a multi-split approach like for BAOAB-RESPA1 to push forward the stability limit [34].

Here we will focus on an alternative approach that involves replacing the Langevin dynamics by a hybrid model combining a classical Langevin diffusion and a piecewise deterministic Markov process, in which some long-range, computationally expensive force are treated in an adaptive way at random times. The initial idea
of this strategy was first introduced in [48] in a particular case, and further studied theoretically in [49]. This hybrid model still samples exactly from $\mu$ but can be simulated using a numerical splitting scheme that requires fewer gradient computations per time step, with a precision of the same order (in the time step) as the classical splitting schemes of the Langevin diffusion such as BAOAB. Moreover, it can be tuned to be arbitrarily close to the Langevin dynamics in terms of stochastic trajectories (see [50, Theorem 3.6]), which makes it suitable to estimate the dynamical properties of the process (with, of course, a trade-off between the accuracy of these dynamical properties and the numerical cost of the simulation, as would be with any numerical approximation of the Langevin equation). Finally, this versatile framework is parallelizable (allowing GPU implementations) and can be combined with the multi-timestep methods, pushing further the computational speedup while avoiding some of the resonance issues of the latest [41, 42]. This new integrator has been implemented in the Tinker-HP software [29], both on CPU and GPU versions [21], the latter allowing simulations on much larger systems.

The paper is organized as follows. Section 2 is devoted to the description of the numerical method: we define in Sections 2.1 and 2.2 the continuous-time hybrid process, Sections 2.3 and 2.4 describe the simulation techniques and how they yield a computational advantage. Then, in Section 3, we see how our approach can be applied to Coulomb and van der Waals interactions in a classical force field. Finally, Section 4 is devoted to numerical experiments, both on CPU and GPU versions of the Tinker-HP software [29]. The supporting information gives the proof of some of the results and other technical details.

## 1.1 Notations

In all the following, $|\cdot|$ denotes the Euclidean norm in $\mathbb{R}^d$ (or the Frobenius norm when it is applied to matrices), and $\cdot$ the standard dot product. For all $x \in \mathbb{R}$, $(x)_+ = \max(0, x)$ denotes the positive part of $x$.

# 2 Numerical method

## 2.1 The velocity jump Langevin process

Let us consider a system of $N$ interacting atoms. We denote $(x, v) \in \mathbb{R}^{6N}$ their positions and velocities, $M = \text{diag}(m_1 I_3, \dots, m_N I_3)$ their mass matrix, $\beta = 1/(k_B T)$ the inverse temperature of the system (where $k_B$ is the Boltzmann constant), $U: \mathbb{R}^{3N} \to \mathbb{R}$ the potential energy function encoding the interactions between the particles and finally $H(x, v) = U(x) + \frac{1}{2}v^T M v$ the Hamiltonian of the system, corresponding to its total energy.

In the canonical ensemble, the system is described by a probability measure that gives, for any given state $(x, v) \in \mathbb{R}^{6N}$, the probability that the system is in the configuration $(x, v)$. This probability measure is called the Boltzmann-Gibbs measure, defined by

$$d\mu(x, v) = \frac{1}{Z_\mu} \exp(-\beta H(x, v))dxdv, \qquad (1)$$

where $Z_\mu = \int_{\mathbb{R}^{6N}} \exp(-\beta H(x, v))dxdv$. Macroscopic quantities are then described as expectations of an observable with respect to this Gibbs measure. In many cases, computing them analytically is impossible, since they are high-dimensional integrals involving an unknown constant $Z_\mu$. However, it is possible to approximate $E_\mu[\varphi]$ by simulating a long trajectory of a Markov process $(X_t, V_t)_{t \ge 0}$ that is ergodic with respect to $\mu$, which means that for any $\varphi$ (in some suitable class of functions), almost surely:

$$\frac{1}{t} \int_{0}^{t} \varphi(X_s, V_s)ds \xrightarrow{t \to \infty} \int_{\mathbb{R}^{6N}} \varphi(x, v)d\mu(x, v) = E_{\mu}[\varphi].$$

A popular process that has this property (under mild conditions on the potential $U$) is the Langevin dynamics, defined as the solution of the following SDE:

$$\left\{ \begin{aligned} dX_t &= V_t dt \\ dV_t &= -M^{-1}\nabla U(X_t)dt - \gamma V_t dt - M^{-1/2}\sqrt{2\gamma\beta^{-1}}dW_t, \end{aligned} \right. \qquad (2)$$

where $(W_t)_{t \ge 0}$ is a standard Brownian motion in $\mathbb{R}^{3N}$ and $\gamma > 0$ is a friction parameter.

In practice, in many molecular simulations, periodic boundary conditions are enforced (typically to ensure bounded solutions in systems that do not have a confining potential), which means that the position lies in the periodic flat torus $\mathbb{T}^{3N} = \mathbb{R}^{3N}/\mathbb{Z}^{3N}$. The measure (1) and the SDE (2) are then understood in $\mathbb{T}^{3N} \times \mathbb{R}^{3N}$.
The following general procedure was first described in [48]. The first step is to decompose the forces as a sum of $K \ge 1$ vector fields $F_i$:

$$ \nabla U(x) = \sum_{i=0}^{K} F_i(x), $$

such that, typically, $F_0$ gathers the computationally inexpensive components (e.g. short-range interactions exhibiting fast variation), while the $F_i$ terms, for $i \ge 1$, represent longer-range forces, which are more numerically intensive than $F_0$ (as each atom interacts with all others through these forces, unlike the short-range ones).

We now introduce the velocity jump Langevin process $(X_t, V_t)_{t \ge 0}$ as the following. The dynamics follows the Langevin diffusion associated with the force $F_0$:

$$ \left\{ \begin{aligned} dX_t &= V_t dt \\ dV_t &= -M^{-1}F_0(X_t)dt - \gamma V_t dt - M^{-1/2}\sqrt{2\gamma\beta^{-1}}dW_t, \end{aligned} \right. \quad (3) $$

with the velocity $V_t$ undergoing additional random jumps at rate $\lambda_i(x, v)$ following a jump kernel $q_i(x, v, dv')$ (that describes the probability distribution of the velocity after a jump), both defined below in Section 2.2, as also explained in [48] in a particular case. The jump mechanism $(\lambda_i, q_i)$ depends on the force $F_i$ in a way that ensures that the equilibrium measure of the process in indeed the canonical measure (1). Informally, the process follows (3) but, between times $t$ and $t + \delta$ for a small $\delta$, with probability $\lambda_i(X_t, V_t)\delta + o(\delta)$, its velocity $v$ is re-sampled according to the kernel $q_i(X_t, V_t, dv')$. It can be seen as a hybrid between a diffusion and a piecewise deterministic Markov process (PDMP) such as considered in [51, 52, 53], which have recently drawn much attention in numerical probabilities and Bayesian statistics, see [52, 54, 55, 56, 57, 58]. As we will explain, simulating this process in a certain way yields a computational advantage over the simulation of the Langevin diffusion, because the computation of the force $F_i$ will only be required when a jump of type $i$ is proposed, which will not happen at each time step for every $i$. We refer to the SI, Section 1, for a more rigorous definition.

## 2.2 The jump mechanism

The velocity jump Langevin process defined in Section 2.1 was described in [48] in the particular case where the jump mechanism corresponds to the Bouncy Particle Sampler (BPS) [59], i.e. with

$$ \lambda_i(x, v) = \beta(v \cdot F_i(x))_+, \quad (4) $$

and deterministic jumps $q_i(x, v, dv') = \delta_{R_i(x,v)}(dv')$, where

$$ R_i(x, v) = v - 2 \frac{v \cdot F_i(x)}{F_i(x) \cdot M^{-1} F_i(x)} M^{-1} F_i(x). \quad (5) $$

Although this jump mechanism gives good results in a certain framework [48], it also suffers some limitations, that we will explain below, in Section 3. We now define a more general jump mechanism based on [50], that can be seen as an interpolation between the BPS and Hamiltonian dynamics. Let $\varepsilon(x) > 0$ be a positive function of the position. For each $1 \le i \le K$, We define the following jump rate:

$$ \lambda_i(x, v) = \frac{\sqrt{\beta}}{\varepsilon(x)} |M^{-1/2} F_i(x)| \Theta(\eta), \quad (6) $$

where

$$ \eta = \varepsilon(x) \sqrt{\beta} \frac{v \cdot F_i(x)}{|M^{-1/2} F_i(x)|}, $$

and $\Theta(\eta) = \mathbb{E}[(\eta + G)_+]$ with $G \sim \mathcal{N}(0, 1)$ a standard Gaussian variable. The transition kernel $q_i(x, v, dv')$ is such that at the moment of a jump, the velocity is updated as

$$ v \leftarrow v - \frac{2\varepsilon(x)}{\sqrt{\beta(1 + \varepsilon^2(x))}} \left(\eta + \tilde{G}\right) \frac{M^{-1} F_i(x)}{|M^{-1/2} F_i(x)|}, \quad (7) $$

where $\tilde{G}$ is a one-dimensional random variable with density

$$ f_{\eta}(y) = \frac{1}{\Theta(\eta)\sqrt{2\pi}} (\eta + y)_+ e^{-y^2/2}. $$
In practice, depending on the context, the parameter $\varepsilon$ will be chosen to be either a positive constant or proportional to the norm of the force: $\varepsilon(x) = \varepsilon_0|F_i(x)|$ with $\varepsilon_0$ a positive constant. It is shown in [50] that in particular, $\varepsilon = \infty$ corresponds to the BPS (i.e. the jump process described in (4) and (5)), and that when $\varepsilon \to 0$, the corresponding velocity jump process (with a free transport on positions) converges to the Hamiltonian dynamics associated with the force $F_i$, i.e the process solving:

$$ \left\{ \begin{array}{l} dX_t = V_t dt \\ dV_t = -M^{-1}F_i(X_t)dt. \end{array} \right. $$

Moreover, as explained in SI, section 4, the definitions (6) and (7) ensure that the Gibbs measure $\mu$ is left invariant by the velocity jump Langevin process defined in Section 2.1, and the ergodicity of the process in a theoretical setting is proven in [49]. In other words, when $\varepsilon$ is small enough, the velocity jump Langevin process is close, in terms of stochastic trajectories, to the Langevin dynamics, which yields (as we will see in Section 4), a good preservation of dynamical properties of the system on top of a computational speedup without loss of accuracy in sampling.

## 2.3 Splitting schemes and integration to multi-time-step methods

### 2.3.1 JUMP integrator

In practice, since the solution of the SDE (3) cannot be simulated exactly, time discretization methods need to be used to approximate the process. Similarly to the BAOAB [38, 39, 40] scheme, the continuous-time dynamics of the velocity jump Langevin process can be approximated following a Trotter/Strang splitting scheme. We will describe three of them, and refer to the SI, section 2, for more technical details.

The first splitting consists in decomposing the dynamics into a free transport (A), an acceleration due to $F_0$ (B), a jump (J) and a friction-dissipation part (O). For a given time-step $\delta$, one step of the numerical scheme is given by the succession of steps BJAOAJB, where:

$$ (B) \quad v \leftarrow v - \frac{\delta}{2} M^{-1} F_0(x) \quad (\text{force } F_0) $$

$$ (J) \quad v \leftarrow W_{\delta/2} \quad (\text{jumps}) $$

$$ (A) \quad x \leftarrow x + \frac{\delta}{2} v \quad (\text{free transport}) $$

$$ (O) \quad v \leftarrow e^{-\gamma \delta} v + \sqrt{\beta^{-1}(1 - e^{-2\gamma \delta})} M^{-1/2} G \quad \text{with } G \sim \mathcal{N}(0, I_d) \quad (\text{friction/dissipation}), $$

and where $(W_s)_{s \ge 0}$ is the piecewise-constant process initialized at $W_0 = v$, that jumps at rate $\lambda_i$ according to $q_i$ for each $i \in [1, K]$. Such a scheme is similar to BAOAB, but where two half-time jump steps are added between the forces and the transport part. As explained and proved in [49], this palindromic form gives rise to a second-order scheme in the time-step, that is, the discretization bias on the invariant measure is of order $\delta^2$. This was the procedure used in [48] for the particular jump process defined in equations (4) and (5). For numerical stability reasons that we will detail in Section 3, it will be sometimes useful to merge the (B) and (J) parts into one (C) step, and run a splitting scheme of the form CAOAC, where the (C) corresponds to

$$ (C) \quad v \leftarrow \widetilde{W}_{\delta/2}, $$

with $(\widetilde{W}_s)_{s \ge 0}$ is defined as $(W_s)_{s \ge 0}$, but such as between jumps, it undergoes the constant acceleration $-M^{-1}F_0(x)$ instead of being constant.

In fact, the BJAOAJB and CAOAC schemes can be combined: we can for instance treat some of the forces $F_i$ in (C), and the others separately in (J), that is, doing a scheme of the form C'J'AOAJ'C', where (J') treats, through a piecewise-constant velocity jump process, the forces $F_i$ for $i \in [1, K']$ (with $K' < K$) and (C') treats the $F_i$ for $i \in [K' + 1, K]$ and $F_0$ in the same way as the step (C) described above. This situation will also be encountered in the following applications, see Section 3. We will now refer to those splitting schemes (BJAOAJB, CAOAC or C'J'AOAJ'C') as the JUMP integrator.

### 2.3.2 JUMP-RESPA and JUMP-RESPA1 integrators

All these steps can be integrated in a classical multi-timestep framework. Let us say that $F_0$ is itself a sum of two distinct terms $F_0 = F_{0,0} + F_{0,1}$, where, for instance, $F_{0,1}$ corresponds to a slow-varying many-body force,
that cannot be treated efficiently with velocity jumps, but that doesn't need to be computed as often as the fast-varying forces gathered in $F_{0,0}$. In that case, let $Q$ be the transition kernel corresponding to a splitting scheme (like BJAOAJB, CAOAC or C'J'AOAJ'C') associated to $\nabla U - F_{0,1}$ with a certain time-step $\delta$. Let $\Delta = n\delta$ (with $n \in 2\mathbb{N}$) be a larger time step. The principle of the JUMP-RESPA scheme is to use the following splitting:

1. $v \leftarrow v - \frac{\Delta}{2} M^{-1} F_{0,1}(x)$.

2. Perform $n$ times the step $Q$ with time step $\delta$.

3. $v \leftarrow v - \frac{\Delta}{2} M^{-1} F_{0,1}(x)$.

This procedure can itself be iterated to include more different time steps, such as BAOAB-RESPA1 [34] that has three layers of time steps. By analogy, We will refer to this situation as the JUMP-RESPA1 integrator. More details are given in SI, section 2.

Let us now see how to simulate efficiently the processes $(W_s)_{s \ge 0}$ and $(\tilde{W}_s)_{s \ge 0}$, and how this gives rise to a computational advantage.

## 2.4 Efficient simulation of jumps

Let us start with the process $(\tilde{W}_s)_{s \ge 0}$ (corresponding to the (C) step in the CAOAC scheme). We suppose that each jump rate $\lambda_i$ is bounded by a constant $\bar{\lambda}_i$, and denote $\bar{\lambda} = \sum_i \bar{\lambda}_i$. Then an easy computation on the generator (see SI, section 3) shows that the jumps can be simulated exactly this way: starting from $(x, v)$,

1. Draw $\mathcal{E}$ a standard exponential random variable, and let $T = \mathcal{E}/\bar{\lambda}$ be the next jump time proposal.

2. Draw $I$ such that $\mathbb{P}(I = i) = \bar{\lambda}_i/\bar{\lambda}$. Propose a jump of type $I$ at time $T$.

3. The jump is accepted with probability $\lambda_I(x, v - M^{-1}F_0(x)T)/\bar{\lambda}_I$, in which case the velocity is resampled at time $T$ according to $q_I(x, v - M^{-1}F_0(x)T, dv')$, otherwise the velocity at time $T$ is simply $v - M^{-1}F_0(x)T$ (there is no jump).

The construction is then repeated by induction over the jump time proposals. The case of $(W_s)_{s \ge 0}$ is very similar: the only difference is that since the process is constant between jumps, in the step 3, a jump of type $I$ is accepted with probability $\lambda_I(x, v)/\bar{\lambda}_I$, in which case the velocity is resampled according to $q_I(x, v, dv')$. This method, called Poisson thinning [60, 61] (which is exact, in the sense that it does not induce any further time discretization error) is the key to the computational advantage of the algorithm: contrary to a classical numerical scheme of the Langevin diffusion (such as BAOAB) where the gradient $\nabla U$ is computed at each time step, here we only need to evaluate $F_i$ when a jump of type $i$ is proposed (at the step 3 above), which does not happen at each time step for each $i$. However, notice that similarly to BAOAB and contrary to multi-timestep splitting methods, there is still a unique time step in the discretization and no additional parameters to tune: the frequency at which $F_i$ is evaluated is random and is adapted to each force.

In many cases, we will use the same bound for all the forces $F_i$, so the $\bar{\lambda}_i$ terms will be uniform in $i: \bar{\lambda}_i = \lambda^*$ (so that $\bar{\lambda} = K\lambda^*$). In those cases, the step 2 above is simply a uniform draw among all the $i$. The Algorithm 1 describes the (C) step of the CAOAC scheme, and the Algorithm 2 describes the (J) step of the BJAOAJB scheme.

In fact, as shown in SI, section 6, in the general case, the bounds $\bar{\lambda}_i$ are local (they depend the velocity $v$) and need to be re-evaluated after each jump. This is not a problem, since the only requirement in order to have a computational advantage is that those bounds do not depend on $F_i$. However, in the case of the Bouncy Particle Sampler (i.e. the jump mechanism defined in (4) and (5)), the norm of the velocity doesn't change after a jump, so $\bar{\lambda}$ stays constant throughout the time step. In those cases, when we simulate the process up to a time $\delta$, we don't need to know exactly the jump times, but only the number of jumps and the order in which they occur. Classical properties of the exponential distribution ensure that during the time interval $[0, \delta]$, the number of jump proposals (step 1 above) follows a Poisson distribution of parameter $\bar{\lambda}\delta$. Then the simulation procedure for the (J) part of BJAOAJB is the following: starting from $(x, v)$,

1. Draw $M \sim \text{Poisson}(\bar{\lambda}\delta)$ (the total number of jump proposals).

2. For each $j \in [1, M]$

(a) Draw $I$ such that $\mathbb{P}(I = i) = \bar{\lambda}_i/\bar{\lambda}$. The $j$-th proposal is a jump of type $I$.
**Algorithm 1** Simulation of (C) part.

Simulates the jump process $(\bar{W}_s)_{s\ge0}$ for a time $\delta/2$, starting from $(x, v)$.

1: **procedure** VELOCITY-UPDATE($\delta, x, v$)
2:     $T \leftarrow 0, t \leftarrow 0$ &rarr; $t$ is the current time, $T$ is the time of jump proposals
3:     $a \leftarrow -M^{-1}F_0(x)$ &rarr; $a$ is the acceleration due to the short range forces
4:     **while** $t \le \delta/2$ **do**
5:         $T \leftarrow \mathcal{E}(K\lambda^*)$ &rarr; Draw the next jump proposal time: an exponential law of parameter $\bar{\lambda}$
6:         **if** $T + t > \delta/2$ **then**
7:             $v \leftarrow v + (\delta/2 - t)a$
8:             **else**
9:                 $t \leftarrow t + T$ &rarr; Update the current time
10:             $v \leftarrow v + Ta$ &rarr; Update the velocity until the jump proposal
11:             $I \leftarrow \text{RANDOM}([1, K])$ &rarr; Choose a jump type
12:             $U \leftarrow \text{RANDOM}([0, 1])$
13:             $\lambda \leftarrow \lambda_I(x, v)$ &rarr; compute the I-th jump rate at the current state
14:             **if** $U \le \lambda/\lambda^*$ **then**
15:                 $v \leftarrow \text{SAMPLING}(q_I(x, v, dv'))$ &rarr; the velocity is resampled
16:             **end if**
17:         **end if**
18:     **end while**
19: **end procedure**

**Algorithm 2** Simulation of (J) part (general case).

Simulates the piecewise-constant jump process $(W_s)_{s\ge0}$ for a time $\delta/2$, starting from $(x, v)$.

1: **procedure** JUMP($\delta, x, v$)
2:     $T \leftarrow 0, t \leftarrow 0$ &rarr; $t$ is the current time, $T$ is the time of jump proposals
3:     **while** $t \le \delta/2$ **do**
4:         $T \leftarrow \mathcal{E}(K\lambda^*)$ &rarr; Draw the next jump proposal time: an exponential law of parameter $\bar{\lambda}$
5:         **if** $T + t > \delta/2$ **then**
6:             $t \leftarrow \delta/2$
7:             **else**
8:                 $t \leftarrow t + T$ &rarr; Update the current time
9:                 $I \leftarrow \text{RANDOM}([1, K])$ &rarr; Choose a jump type
10:             $U \leftarrow \text{RANDOM}([0, 1])$
11:             $\lambda \leftarrow \lambda_I(x, v)$ &rarr; compute the I-th jump rate at the current state
12:             **if** $U \le \lambda/\lambda^*$ **then**
13:                 $v \leftarrow \text{SAMPLING}(q_I(x, v, dv'))$ &rarr; the velocity is resampled
14:             **end if**
15:         **end if**
16:     **end while**
17: **end procedure**
**Algorithm 3** Simulation of (J) part (constant bound case)

Special case: a jump does not affect the bound on the jump rate.

Simulates the piecewise-constant jump process $(W_s)_{s \ge 0}$ for a time $\delta/2$, starting from $(x, v)$.

1: **procedure** JUMP($\delta, x, v$)

2: $M \leftarrow POISSON(K\lambda^*\delta/2)$

3: **for** $j \in [1, M]$ **do**

4: $I \leftarrow RANDOM([1, K])$

5: $U \leftarrow RANDOM([0, 1])$

6: $\lambda \leftarrow \lambda_I(x, v)$

7: **if** $U \le \lambda/\lambda^*$ **then**

8: $v \leftarrow SAMPLING(q_I(x, v, dv'))$

9: **end if**

10: **end for**

11: **return**

12: **end procedure**

▷ For each jump proposal $j$, choose randomly a jump type $I$

▷ compute the I-th jump rate at state $(x, v)$

▷ the velocity is resampled

(b) Accept the proposal (i.e. resample the velocity according to $q_I(x, v, dv')$) with probability $\lambda_I(x, v)/\lambda_I$.

This situation is described in Algorithm 3.

Regarding the simulation of the jump itself, it is shown in [50] that the random variable $\tilde{G}$ that appears in (7) can be easily simulated with rejection sampling. Assuming that we can sample a distribution with density $g$ that satisfies $f(x) \le Cg(x)$ and such that the ratio $f(x)/Cg(x)$ is easy to compute for any $x$, the procedure is the following:

1. Draw a proposal $Y$ with density $g$.

2. Draw $U$ uniform in $[0, 1]$. If $U \le \frac{f(Y)}{Cg(Y)}$ then the proposal is accepted.

3. Otherwise, repeat from step 1.

The mean number of proposals before acceptation is $C$, so the smaller the constant the faster the procedure is. Hence, following [50] (to which we refer for details), in our case, $g$ is a:

$$ \left\{ \begin{array}{ll} \text{Gamma law} & \text{for } \eta < -2.5 \\ \text{Exponential law} & \text{for } -2.5 \le \eta < -1 \\ \text{Rayleigh law} & \text{for } -1 \le \eta \le 0 \\ \text{Mixed Rayleigh-Gaussian law} & \text{for } \eta > 0. \end{array} \right. $$

## 3 Application to electrostatic and van der Waals interactions

In a classical (non-polarizable) force field such as OPLS [14], AMBER [62] or CHARMM [12], the non-bonded parts of the interaction potential are the Lennard-Jones and electrostatic contributions. Since the Lennard-Jones term decreases very fast, in periodic boundary conditions one can simply compute the forces in a radius up to a certain cutoff. The electrostatic part, decreasing at a much slower rate, is more problematic. A classical way of treating the long range interactions is by using Ewald summation techniques, for instance through the SPME method [63] (although other approaches exist, see for instance [64]). In this method, the Coulomb potential is decomposed in a direct and a reciprocal part, the latter being computed in the Fourier space with a fast Fourier transform. In practice, the algorithm BOUNCE described in [48], i.e. a BJAOAJB splitting on the process with the BPS jump mechanism, has been implemented in Tinker-HP [29] and tested in periodic boundary conditions with jumps on the long-range parts of the van der Waals interactions and the direct part of the electrostatic interactions. The reciprocal part was treated with a classical multi-timestep method, the reference system propagator algorithm (RESPA)[35]. Since electrostatics are much stronger than van der Waals at a long range, the BPS electrostatic jumps can only be done on a very small range. Jumping on a larger range yields too many velocity bounces, in the sense that the numerical scheme with the classical timestep of $1f$s is very unstable. One solution would be
to decrease the timestep, which is of course counterproductive, since the goal is to reduce the final cost of the calculations.

On the contrary, the jump mechanism defined in (6) and (7) yields milder jumps than the orthogonal reflexions of the BPS. Moreover, as discussed earlier, the process can be seen as an interpolation between bounces and diffusions: in the situations where BOUNCE is numerically unstable, choosing a small parameter $\epsilon$ makes the process closer to classical Hamiltonian dynamics. In addition to that, combining the (B) and (J) parts (i.e. doing the CAOAC scheme instead of BJAOAJB) stabilizes further the discretized process.

These two elements combined, namely a more general jump process and the fusion of the “jump” and “force” steps, allow to replace a larger range of the potential by jumps than in the BOUNCE algorithm, yielding more computational gain, while preserving dynamical properties of the system such as the diffusion coefficients. Finally, the BPS jumps of [48] on the van der Waals forces can also be integrated in the general procedure, and two multi-timestep versions of the algorithm, namely JUMP-RESPA and JUMP-RESPA1 combine all type of jumps and yields a computational gain over the BAOAB-RESPA and BAOAB-RESPA1 [34] algorithm.

## 3.1 Direct electrostatic jumps

As reminded in the SI, section 5, the Ewald summation splits the electrostatic interactions into a direct and a reciprocal part. In Periodic Boundary Conditions (PBC), this decomposition can be chosen in a way that enforces the “minimum image convention”: in the direct space, each atom interacts only with the closest image of an other atom. We assume in the following that this condition is satisfied. By denoting $r_{ij} = |x_i - x_j|$ the distance between the atoms i and j, the direct electrostatic energy associated to each pair $(i, j)$ is:

$$U_{ij}(x) = \frac{q_i q_j}{r_{ij}} \text{erfc}(\alpha r_{ij}).$$

In order to decompose this term into a short and a long range contribution, let $r_c$ be a cutoff (smaller than the “Ewald cutoff” that separates the direct and reciprocal part, whose value is enforced by the choice of the parameter $\alpha$), $\chi(r)$ a switching function that goes smoothly from 1 to 0 around the cutoff (its precise formula is given in SI, section 6) and let

$$U_{ij} = \chi(r_{ij})U_{ij} + (1 - \chi(r_{ij}))U_{ij}.$$

We now denote $F_{ij} = \nabla_{x_i}(1 - \chi(r_{ij}))U_{ij} \in \mathbb{R}^3$, the component of the long-range force of the atom j acting on the atom i, Therefore, the direct part of electrostatic interactions is decomposed into

$$\nabla U_{direct} = F_0 + \sum_{i \neq j} A_i F_{ij},$$

where $A_i \in \mathcal{M}_{3N \times 3}$ is the matrix with all coefficients are equal to zero except $A(3i - 2, 1) = A(3i - 1, 2) = A(3i, 3) = 1$. Therefore, a jump mechanism is associated to each pair interaction $A_i F_{ij}$. In this context, (6) and (7) have simpler expressions than in the general case described in Section 2. Here, $|M^{-1/2}F_{ij}(x)| = \frac{|F_{ij}(x)|}{\sqrt{m_i}}$ (with $m_i$ the mass of the atom i), which yields the following jump rate:

$$\lambda_{ij}(x, v) = \frac{\sqrt{\beta}}{\varepsilon(x)\sqrt{m_i}} |F_{ij}(x)| \Theta \left( \varepsilon(x) \sqrt{\beta m_i} \frac{v_i \cdot F_{ij}(x)}{|F_{ij}(x)|} \right),$$

and a jump only modifies the atom i:

$$v_i \leftarrow v_i - \frac{2\varepsilon(x)}{1 + \varepsilon^2(x)} \left( \varepsilon(x) \frac{v_i \cdot F_{ij}(x)}{|F_{ij}(x)|} + \frac{\tilde{G}}{\sqrt{\beta m_i}} \right) \frac{F_{ij}(x)}{|F_{ij}(x)|},$$

where $\tilde{G}$ has been generated from the distribution $f_\eta$ with $\eta = \varepsilon(x)\sqrt{\beta m_i \frac{v_i \cdot F_{ij}(x)}{|F_{ij}(x)|}}$. As detailed in the SI, section 6, in this case it is easy to obtain explicit analytic bounds on $|F_{ij}|$ that are required in the thinning method. Moreover, for each $i$, those bounds are uniform in $j$ which, as explained in Section 2.4, simplifies the thinning procedure.

Since the jump mechanisms associated with $F_{ij}$ only involve the i-th atom, the jump processes associated to the different atoms can be seen as independent Markov chains that can be simulated in parallel, which in particular allows massively parallel implementations on GPU, see Section 4.
Finally, as mentioned in the introduction of this section, in this case it is more convenient to perform a CAOAC scheme instead of BJAOAJB. If we treat a large portion of the interactions with jumps, the rate $\lambda_{ij}$ will be higher, yielding many velocity jumps, so merging the (J) and (B) steps of the splitting scheme has the effect of stabilizing the dynamics. The Algorithm 4 below describes precisely the procedure when $\varepsilon$ is constant. The case of an adaptive $\varepsilon(x) = \varepsilon_0|F_i(x)|$ is exactly the same but with a slightly different expression for the jump rate, the bound and the kernel.

**Algorithm 4 Step (C): a general velocity update (with constant $\varepsilon$)**

Simulates during a time $\delta/2$ the PDMP treating the acceleration due to the short range forces and the jumps due to the long range forces. The position $x$ is fixed.

**Input:** N the number of atoms, $(x, v)$ their positions and velocities, $\delta$ the time step

1: **procedure** VELOCITYUPDATE($N, \delta, x, v$)

2: **for** $i \in [1, N]$ **do**

3: $T \leftarrow 0, t \leftarrow 0$

▷ $t$ is the current time, $T$ is the time of jump proposals

4: $v_i \leftarrow v_i^0$

▷ $v_i^0$ is the initial velocity of atom $i$

5: $a_i \leftarrow -\nabla U_{\text{short}}(x)/m_i$

▷ $a$ is the acceleration due to the short range forces

6: $L \leftarrow L_i$

▷ The uniform bound on the $|F_{ij}|, 1 \le j \le N$

7: **while** $t \le \delta/2$ **do**

8: $B_v \leftarrow (|v_i|^2 + (\delta/2 - t)^2|a_i|^2 + 2(\delta/2 - t)(v_i \cdot a_i)_+)^{1/2}$

▷ bound on $|v_i|$ during the time $[t, \delta/2]$

9: $\bar{\lambda} \leftarrow L \times N \left( \beta B_v + \frac{1}{\varepsilon} \sqrt{\frac{\beta}{2\pi m_i}} \right)$

▷ The bound on the jump rate

10: $T \leftarrow \mathcal{E}(\bar{\lambda})$

▷ Draw an exponential law of parameter $\bar{\lambda}$

11: **if** $T + t > \delta/2$ **then**

12: $v_i \leftarrow v_i + (\delta/2 - t)a_i$

13: **else**

14: $t \leftarrow t + T$

15: $v_i \leftarrow v_i + T a_i$

16: $\eta \leftarrow \varepsilon \sqrt{\beta m_i \frac{v_i \cdot F_{ij}(x)}{|F_{ij}(x)|}}$

17: $J \leftarrow \text{RANDOM}([1, N])$

18: $U \leftarrow \text{RANDOM}([0, 1])$

19: $p \leftarrow \lambda_{ij}(x, v)/\bar{\lambda} = \frac{|F_{ij}(x)|\Theta(\eta)}{L(B_v \varepsilon \sqrt{m_i} + 1/\sqrt{2\pi})}$

▷ the probability of acceptation

20: **if** $U \le p$ **then**

21: $\tilde{G} \leftarrow \text{REJECTION SAMPLING}(\eta)$

▷ Generate $\tilde{G}$

22: $v_i \leftarrow v_i - \frac{2\varepsilon}{\sqrt{\beta m_i(1+\varepsilon^2)}} \left( \eta + \tilde{G} \right) \frac{F_{i,j}(x)}{|F_{i,j}(x)|}$

▷ Performs a jump

23: **end if**

24: **end if**

25: **end while**

26: **end for**

27: **return**

28: **end procedure**

### 3.1.1 Optimization of the bound: the ring technique

Let $r_E$ be the Ewald cutoff, i.e. the maximum radius that is taken into account in the direct part and $r_c$ the cutoff between the short-range interactions treated in the acceleration part and the long-range interactions treated with the jumps. In order to improve the bound on the jump rate (and therefore reduce the number of jump proposals), we can introduce an intermediate cutoff $r_M$ ($r_c < r_M < r_E$) and compute two different bounds: a bound on the $|F_{ij}|$ such that $r_c - h \le r_{ij} \le r_M$ ($h$ is the switching parameter) denoted $L_M$ and a bound on the $|F_{ij}|$ such that $r_M \le r_{ij} \le r_E$ denoted $L_L$. For each $i$, we denote $N_M(i)$ and $N_L(i)$ the number of atoms $j$ in these two cases respectively. Now, the term $N \times L_i$ in the bound on the jump rate $\bar{\lambda}_i$ can be replaced by $N_M(i)L_M + N_L(i)L_L$. The algorithm (again in the case where $\varepsilon$ is constant) can then be modified the following way:
1. Draw a proposition at rate $\bar{\lambda}_i$.

2. Draw $Y \in \{M, L\}$ such that $P(Y = M) = \frac{L_M N_M(i)}{N_M(i)L_M+N_L(i)L_L}$ and $P(Y = L) = \frac{L_L N_L(i)}{N_M(i)L_M+N_L(i)L_L}$.

3. Draw uniformly $j$ in the right neighbor-lists (namely the “middle-range” or the “long-range” list).

4. Accept a jump of type $(i, j)$ with probability $\frac{|F_{ij}(x)|\Theta(\eta)}{L_Y(B_v\epsilon\sqrt{m_i}+1/\sqrt{2\pi})}$.

**Remark.** This modification requires to build new “middle” and “long” neighbor-lists, in order to know $N_L(i)$ and $N_M(i)$ for each atom and to draw uniformly among them.

## 3.2 Lennard-Jones jumps

The same procedure can be used for van der Waals forces. In this case, the potential is given by:

$$U_{\text{vdW}}(x_1, \dots, x_N) = \sum_{i \neq j} \varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right].$$

for some parameters $\sigma_{ij}, \varepsilon_{ij} > 0$. Since each term of this sum decreases very fast in $r_{ij}$, there is no splitting in direct and reciprocal parts and all the terms beyond a certain cutoff $r_{\text{vdW}}$ (i.e. for $r_{ij} > r_{\text{vdW}}$) are neglected. Again, we can define the pairwise term

$$U_{ij}^{\text{vdW}}(x) = \varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right],$$

choose a cutoff $r_c < r_{\text{vdW}}$ and use the switching function $\chi$ to separate each term into a short-range and a long-range part:

$$U_{ij}^{\text{vdW}} = \chi(r_{ij})U_{ij}^{\text{vdW}} + (1 - \chi(r_{ij}))U_{ij}^{\text{vdW}}.$$

From now, we define $F_{ij} = \nabla_{x_i}(1 - \chi(r_{ij}))U_{ij}$ and build a jump mechanism associated to each couple $(i, j)$. Here, since the van der Waals forces are weaker than the Coulomb ones at long range, there are less stability issues and we can choose the parameter giving the least jumps: $\varepsilon = \infty$, which is the case that corresponds to the BPS. The jump rate associated to the force $F_{ij}$ is:

$$\lambda_{ij}(x, v) = \beta (v_i \cdot F_{ij}(x))_+ ,$$

and the jumps are deterministic: $q_{ij}(x, v, dv') = \delta_{R_{ij}(x, v)}(dv')$ where

$$R_{ij}(x, v) = v - 2 \frac{v \cdot F_{ij}(x)}{|F_{ij}(x)|^2} F_{ij}(x).$$

If we denote $L_i$ a bound on all the $|F_{ij}|$ for $1 \le j \le N$, the jump rate can be bounded by:

$$\lambda_{ij}(x, v) \le \beta |v| L_i := \lambda_{ij}^*.$$

A jump being an orthogonal reflexion of $v$ against $F_{ij}(x)^{\perp}$, it does not modify $|v|$ and the bound $\lambda_{ij}^*$ can stay the same during all the time step. It is therefore advantageous to use a BJAOAJB splitting with the Algorithm 3 described in 2.4. Moreover, as in the case of jumps on the direct electrostatic part, $\lambda_{ij}^*$ is the same for all $j$, we can define $\bar{\lambda}_i = N\lambda_{ij}^*$ and the process associated to the velocity of each atom can be simulated independently. This yields the Algorithm 5, that was also described in [48].

## 3.3 Conclusion: three new integrators

To conclude, using the jump mechanisms described in this section gives rise to three new integrators. First, the JUMP integrator consists in running either a BJAOAJB scheme, where the (J) part treats the Lennard-Jones jumps, a CAOAC scheme, where the (C) part treats the direct electrostatic jumps (along with the acceleration due to the rest of the forces) or, combining both, a C'J'AOAJ'C' scheme, where (J') treats the Lennard-Jones jumps, and (C') the electrostatic direct jumps. Regarding the multi-timesteps versions, JUMP-RESPA simply consists in treating in a small time-step $\delta$ all the intramolecular forces of the potential (bond stretching, angle
**Algorithm 5** Step (J): BPS for van de Waals interactions

Simulates the piecewise-constant jump process for a time $\delta/2$, starting from $(x, v)$.

**Input**: N the number of atoms, $(x, v)$ their positions and velocities, $\delta$ the time step

1: **procedure** JUMP($N, \delta, x, v$)

2: **for** $i \in [1, N]$ **do**

3: $M_i \leftarrow \text{POISSON}(N\beta|v_i|L_i\delta/2)$

▷ Draw the number of jump proposals

4: **for** $k \in [1, M_i]$ **do**

5: $J \leftarrow \text{RANDOM}([1, N])$

▷ For each jump proposal $k$, choose randomly a jump of type $J$

6: $U \leftarrow \text{RANDOM}([0, 1])$

7: $p \leftarrow \frac{(v_i \cdot F_{iJ}(x))_{+}}{|v_i|L_i}$

▷ compute the probability of accepting the jump

8: **if** $U \le p$ **then**

9: $v_i \leftarrow R_{ij}(x, v)$

▷ the velocity bounces

10: **end if**

11: **end for**

12: **end for**

13: **return**

14: **end procedure**

bending, bond-angle cross terms, out-of-plane bending and torsional rotations), and leaving all the non-bonded forces (namely the van der Waals and Coulomb interactions) in a larger time step $\Delta = n\delta$ (with $n \in 2\mathbb{N}$), including the jumps. However, when using jumps, this is not very natural, as it increases the issues related to time discretization (unstability or numerical bias) while having no effect on the complexity of the jump parts, since with the thinning method, a step-size $n\delta$ requires on average $n$ times more jumps (hence computations) than a step-size $\delta$.

This remark leads to the following scheme: in JUMP-RESPA1, a small time step $\delta$ treats all the intramolecular forces, an intermediate time step $\kappa = m\delta$ treats the short-range van der Waals, the short-range direct electrostatics and the jumps (that take into account the long range van der Waals forces and the long range direct electrostatics), and a large time step $\Delta = n\kappa$ treats the reciprocal electrostatics, the many-body force that is, in practice (especially in large systems) the most numerically expensive force to compute. If we applied the algorithm to a polarizable force field such as AMOEBA [16], the expensive polarization part would also be treated in this largest time step.

In comparison, in BAOAB-RESPA1, the long range van der Waals and long range direct space electrostatics are treated in the larger time step. Here, in JUMP-RESPA1, those are treated by the jumps, in the intermediate time step. As discussed above, indeed, the jump process has the property of not requiring more computations when put in a smaller time step: the jump rate stays exactly the same.

# 4 Numerical experiments

All variations of the algorithms described in Section 3 have been implemented in Tinker-HP software for molecular dynamics, both on the CPU and the parallel GPU versions. The software will be made freely available to academic users within the next release of the Tinker-HP code [65, 66, 21, 29]. Computations were done on the computing cluster of the Theoretical Chemistry Laboratory (LCT) of Sorbonne University, Paris, on a 22 cores Intel(R) Xeon(R) Silver 4116 CPU @ 2.10GHz for the CPU version, and a Nvidia Quadro GV100 for the GPU version. All numerical tests were done on water systems in the NVT ensemble with temperature $T = 300K$, of sizes ranging from 216 to 96000 molecules, with SPC/Fw model in the OPLS-AA force-field, in periodic boundary conditions. Simulation parameters are given at the end of this section, and SI, section 7, gives more details on the implementation.

**Tuning of parameters.** As explained in Section 2.4, treating certain forces with jumps allows to reduce the total number of computations, and therefore accelerates the simulation. However, replacing a too large range of interactions by jumps yields either numerical instabilities (as it is the case for the BOUNCE algorithm described
in [48]) or a significant loss in the auto-diffusion constant, which indicates a reduction of the sampling rate. As mentioned in [34] the ideal situation is when the acceleration rate is higher than the loss rate in the diffusion constant.

In order to do that, the choice of an adaptive $\epsilon$ helps the conservation of the diffusion, as the jump process will yield milder jumps in the shorter-range regions (where the forces are stronger). Moreover, the section of pairwise interaction that are replaced by jumps should not be too large, or in too short-range regions. In other words, there should not be too many jumps, either because of too many jump processes running in parallel (if too many interactions are treated by jumps) or because the jump rates are too high (if those interactions are too strong). Those settings, although they can push further the computational speedup, always yields a loss in auto-diffusion.

In other words, the best choice is to always choose an adaptive $\epsilon$ parameter, with a small enough $\epsilon_0$ (so that the trajectories stay close to the Langevin process), but not too small so there is still a computational advantage, and to replace either a small portion of relatively short-range electrostatic interactions (for example between 5 and 7 Angstroms), or a larger portion of longer-range interactions (for example between 7 and 11 Angstroms).

The first setting seems to be the best for small systems, and the second one is more appropriate for larger systems, especially on the GPU version. Indeed, treating very long range interactions with jumps allows to reduce the reciprocal Coulomb force, which is very expensive in large systems (it can take up to 50% of the total computation time). Moreover, since the reciprocal Coulomb force is partly responsible for resonance issues in multi-timestep integrators, reducing this part allows to extend the stability limit of the external time step for the JUMP-RESPA1 algorithm, which allows for a further increase in simulation speed.

**Precision in sampling and dynamical properties.** Figure 1 shows the oxygen-oxygen radial distributions, the mean temperature and mean potential energy of 500 water molecules, comparing the new JUMP, JUMP-RESPA and JUMP-RESPA1 integrators to its classical BAOAB, BAOAB-RESPA and BAOAB-RESPA1 counterparts, as well as the VERLET integrator for reference, with $2ns$ of simulation. To illustrate the non negligible influence of the forces that are replaced by jumps, this figure also contains a radial distribution where long-range direct electrostatic interactions between 5 and 9 Angström are removed from the potential. Table 1 shows a comparison in hydration free energy of a water molecule, of a sodium ion Na⁺ and of a potassium ion K⁺ in a system of 216 water molecules simulated for $2ns$, using the Bennet Acceptance Ratio [67] (BAR) method. Note that the JUMP integrators are directly compatible with more advanced free energy methods [68]. The canonical measure is therefore similarly sampled by the JUMP integrators and the classical ones. Figure 2 shows that the auto-diffusion constants (computed using the Einstein formula) are well preserved, thanks to the proximity of the process with the Langevin dynamics, as explained in Section 2.1. Again, the tests are done on a system of 500 water molecules, simulated for $2ns$.

**Resonance effects.** One of the well known issues of multi-timestep methods are resonance effects between internal periodicities of the molecular system and the non-physical periodicities created by the various time steps, that can induce numerical instabilities (such as abnormal fluctuations, or even explosion of the kinetic energy) and limit the largest time step that can be taken in the simulation. The reciprocal Coulomb part of the electrostatic interactions is partly responsible for those resonances (recall that this force is not purely long-range, since the erf function is not equal to zero around zero, and that the Ewald sum indeed implies all the atoms of the system). For this reason, increasing the real-space cutoff and replacing a part of the reciprocal electrostatic interactions by jumps, create much less oscillations in the system, thanks to the intrinsically random nature of velocity jumps, which allows to take a larger external time step in JUMP-RESPA1 than in BAOAB-RESPA1 and accelerate further the simulation.

These phenomena can be well observed on velocity autocorrelation spectra (computed from the MD trajectory at the level of the smallest time step using the implementation from ref. [69]), illustrated in Figures 3 and 4. On each graph, from left to right, the first band (below $1000 \text{ cm}^{-1}$) corresponds to molecule libration and slow fluctuations in the hydrogen bond network, the first peak (around $2000 \text{ cm}^{-1}$) corresponds to intra-molecular bending oscillations, the second and third peaks (around $3500 \text{ cm}^{-1}$) correspond to bond stretching motions. All the peaks at the right of those are non-physical artifacts due to the multi time-step. On the zoomed graphs, in BAOAB-RESPA1, the main non-physical peak moves to the left when the external time step increases, coming closer to the physical stretching peak, and even couples to it when the time step is too large. In JUMP-RESPA1, this phenomenon is largely suppressed, resulting in better numerical stability.

**Performances of the algorithm.** In terms of simulation speed, the CPU version gives its best results on small systems (35.5% acceleration of JUMP with respect to BAOAB, and 21.6% of acceleration of JUMP-RESPA1 with
respect to BAOAB-RESPA1 on the system of 1500 atoms), but this advantage disappears when looking at larger systems (with almost no difference between the algorithms when simulating 288000 atoms). On the highly parallel GPU implementation however, it is the opposite: while there are almost no differences between JUMP and BAOAB on small systems (mostly because the GPU is not fully exploited on those, since the computational resources are not saturated, so tests on large molecular systems are more appropriate to this setting), the advantage appears on larger systems, and even is the best on the largest tested system (21.3% acceleration of JUMP with respect to BAOAB and 19.3% acceleration of JUMP-RESPA1 with respect to BAOAB-RESPA1 on the system of 288000 atoms). Table 2 shows the speed of the algorithms (in nanoseconds of simulation per day) on CPU and GPU implementations.

Those differences between GPU and CPU have different reasons. First, the neighbor lists are treated differently, see the end of this section for more details. Then, the proportion of time spent in the different routines of the integrator is different on CPU and GPU. For instance, the parallel architecture of GPUs allows for a more efficient treatment of the van der Waals and direct electrostatic interactions, while the many-body Ewald sum takes almost 50% of the total computation time, especially on larger systems. As a consequence, on a large system, reducing the reciprocal part of electrostatics and treating a larger portion of the potential with jumps yields a significant computational advantage. For reference, Table 3 gives the proportion of time spent in each of the main routines of the CPU and GPU implementations, comparing two system sizes.

As explained in [27], when increasing the real-space cutoff $r_c$, the Fourier space spacing can be multiplied by the same factor, which reduces the cost of its computation. It is particularly advantageous to do that on the large systems, where the interpolation of the structure factors of the Ewald sum is done on large grids.

<table><caption>Table 1: Hydration free energies of water, sodium and potassium in 216 water molecules.</caption><thead><tr><th>Integrator</th><th>BAOAB</th><th>JUMPS</th></tr></thead><tbody><tr><td>ΔG<sub>hydrat</sub> H<sub>2</sub>O (Kcal/Mol)</td><td>6.73 ± 0.09</td><td>6.77 ± 0.08</td></tr><tr><td>ΔG<sub>hydrat</sub> Na+ (Kcal/Mol)</td><td>75.48 ± 0.09</td><td>75.61 ± 0.03</td></tr><tr><td>ΔG<sub>hydrat</sub> K+ (Kcal/Mol)</td><td>58.81 ± 0.04</td><td>58.82 ± 0.02</td></tr></tbody></table>

<table><caption>Table 2: Speed performances of the algorithms, expressed in nanoseconds of simulation per day</caption><thead><tr><th></th><th colspan="4">CPU</th><th colspan="4">GPU</th></tr><tr><th>Number of atoms</th><th>1500</th><th>12000</th><th>96000</th><th>288000</th><th>1500</th><th>12000</th><th>96000</th><th>288000</th></tr></thead><tbody><tr><th>BAOAB (1 fs)</th><td>3.50</td><td>3.95e-1</td><td>3.89e-2</td><td>8.15e-3</td><td>178.0</td><td>98.19</td><td>16.91</td><td>5.13</td></tr><tr><th>JUMP (1 fs)</th><td>4.74</td><td>4.94e-1</td><td>4.52e-2</td><td>8.74e-3</td><td>175.3</td><td>1.769</td><td>20.14</td><td>6.22</td></tr><tr><th>Acceleration</th><td>35.5%</td><td>25.1%</td><td>16.2%</td><td>7.24%</td><td>-1.52%</td><td>9.68%</td><td>19.1%</td><td>21.3%</td></tr><tr><th>BAOAB-RESPA (1:2 fs)</th><td>5.68</td><td>6.29e-1</td><td>6.25e-2</td><td>1.46e-2</td><td>158.4</td><td>110.78</td><td>27.61</td><td>8.59</td></tr><tr><th>JUMP-RESPA (1:2 fs)</th><td>6.62</td><td>6.98e-1</td><td>6.24e-2</td><td>1.38e-2</td><td>154.7</td><td>115.93</td><td>30.58</td><td>9.82</td></tr><tr><th>Acceleration</th><td>16.5%</td><td>11.0%</td><td>-0.16%</td><td>-5.48%</td><td>-2.34%</td><td>4.65%</td><td>10.8%</td><td>14.3%</td></tr><tr><th>BAOAB-RESPA1 (1:2:5 fs)</th><td>6.75</td><td>7.13e-1</td><td>6.73e-2</td><td>1.82e-2</td><td>276.6</td><td>141.3</td><td>34.69</td><td>11.93</td></tr><tr><th>JUMP-RESPA1 (1:2:6 fs)</th><td>8.21</td><td>8.30e-1</td><td>6.58e-2</td><td>1.81e-2</td><td>308.9</td><td>158.6</td><td>41.48</td><td>14.23</td></tr><tr><th>Acceleration</th><td>21.6%</td><td>16.4%</td><td>-2.22%</td><td>-0.55%</td><td>11.7%</td><td>12.2%</td><td>19.6%</td><td>19.3%</td></tr></tbody></table>

<table><caption>Table 3: Time spent in each routine</caption><thead><tr><th rowspan="2">Nb of atoms</th><th colspan="3">BAOAB on CPU</th><th colspan="4">JUMP on CPU</th></tr><tr><th>vdW</th><th>Direct Coulomb</th><th>SPME</th><th>vdW</th><th>Direct Coulomb</th><th>SPME</th><th>Jumps</th></tr></thead><tbody><tr><td>12000</td><td>8.34%</td><td>41.23%</td><td>37.44%</td><td>5.24%</td><td>21.79%</td><td>43.63%</td><td>2.41%</td></tr><tr><td>288000</td><td>13.64%</td><td>33.25%</td><td>40.45%</td><td>7.71%</td><td>38.48%</td><td>23.20%</td><td>3.37%</td></tr><tr><th colspan="8">BAOAB on GPU</th></tr><tr><td>12000</td><td>25.5%</td><td>34.0%</td><td>17.7%</td><td>26.2%</td><td>21.9%</td><td>18.6%</td><td>8.0%</td></tr><tr><td>280000</td><td>8.5%</td><td>32.3%</td><td>52.0%</td><td>6.6%</td><td>46.7%</td><td>28.0%</td><td>7.4%</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>(a) Oxygen-Oxygen radial distributions</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Oxygen-Oxygen radial distributions without Coulomb forces between 5 and 9 Angström</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Mean temperature</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(d) Mean potential energy</figcaption></figure>

<caption>Figure 1: Sampling properties of the JUMP algorithms for a system of 500 water molecules.</caption>
<figure><img src="image_6.png" /><figcaption>Figure 2: Diffusion coefficients for a system of 500 water molecules.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: Velocity autocorrelation spectra of BAOAB-RESPA1 with various external time-steps, for a system of 500 water molecules. From left to right, the first bump corresponds to some many-body frequencies of the system, the first peak corresponds to intra-molecular angle oscillations, the second and third peaks correspond to bond oscillations. All the peaks at the right of those are non-physical artifacts due to the multi time-step. On the zoomed graph, the main non-physical peak moves to the left as the external time step increases, coupling with the bond-oscillation peak at the stability limit.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4: Velocity autocorrelation spectra of BAOAB-RESPA1 and JUMP-RESPA1 with various external time-steps. The non-physical artifacts are largely suppressed in JUMP-RESPA1.</figcaption></figure>
Simulation details

* Tinker-HP version: 1.2

* Force-field: OPLS-AA

* Water model: SPC/Fw

* switching parameter: 0.5Å

* Friction: $1.0ps^{-1}$

* Temperature: 300K.

Parameters for BAOAB, BAOAB-RESPA and BAOAB-RESPA1 simulations:

* Time step for BAOAB: $\delta = 1fs$.

* Time steps for BAOAB-RESPA: $\delta = 1fs$, $\Delta = 2fs$.

* Time steps for BAOAB-RESPA1: $\delta = 1fs$, $\kappa = 2fs$, $\Delta = 5fs$.

* Charge real-space cutoff: 7Å.

* Ewald parameter: $\alpha = 0.5446$.

* Van de Waals long range cutoff: 10Å.

Parameters for JUMP, JUMP-RESPA and JUMP-RESPA1 simulations:

* Time step for JUMP: $\delta = 1fs$.

* Time steps for JUMP-RESPA: $\delta = 1fs$, $\Delta = 2fs$.

* Time steps for JUMP-RESPA1: $\delta = 1fs$, $\kappa = 2fs$, $\Delta = 6.5fs$.

* Van der Waals jumps only on CPU.

* Van de Waals long range cutoff: 10Å.

* Van de Waals short range cutoff (if activated): 7Å.

* Type of jump parameter: adaptive $\epsilon$.

<table><caption>Table 4: JUMP and PME-grid parameters.</caption><tr><td>Number of atoms</td><td>648</td><td>1500</td><td>12000</td><td>96000</td><td>288000</td></tr><tr><td>Box size (in Å)</td><td>18.643</td><td>24.662</td><td>49.232</td><td>98.5</td><td>142.27</td></tr><tr><td>default PME-grid</td><td>24x24x24</td><td>30x30x30</td><td>60x60x60</td><td>120x120x120</td><td>180x180x180</td></tr><tr><td>PME-grid in JUMP</td><td>24x24x24</td><td>30x30x30</td><td>60x60x60</td><td>96x96x96</td><td>128x128x128</td></tr><tr><td>α parameter in JUMP</td><td>0.554</td><td>0.554</td><td>0.554</td><td>0.42</td><td>0.3767</td></tr><tr><td>Real-space cutoff JUMP (in Å)</td><td>7</td><td>7</td><td>7</td><td>9</td><td>10</td></tr><tr><td>Short-range charge cutoff (in Å)</td><td>5</td><td>5</td><td>5</td><td>6</td><td>7</td></tr><tr><td>Jump parameter ε<sub>0</sub></td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.01</td><td>0.01</td></tr></table>
# 5 Conclusion and perspectives

To conclude, our work presents a new framework to construct velocity jumps (or JUMP) integrators for molecular dynamics simulations that can be parallelized on GPU and integrated into multi-timestep methods. It allows to accelerate simulations while preserving sampling precision and dynamical properties, and avoiding certain classical resonance issues. The JUMP approaches provide significant acceleration over their BAOAB, BAOAB-RESPA and BAOAB-RESPA1 counterparts. If the present implementation has been performed in the Tinker-HP package, the approach could be easily ported to other molecular dynamics software. Future work will deal with the extension of the approach to polarizable force fields [66, 70, 19, 20] and with its combination with advanced periodic boundary conditions treatment of electrostatics such as ANKH [72].

**Supporting information** The supporting information contains rigorous definitions on the continuous-time process and splitting schemes, proofs of the mathematical statements as well as technical details on the TinkerHP implementation of the algorithms.

**Conflicts of interest** Louis Lagardère and Jean-Philip Piquemal are co-founders and shareholders of Qubit Pharmaceuticals. The remaining authors declare no competing interests.

**Acknowledgments** This work was made possible thanks to funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810367), project EMC2 (JPP). Computations have been performed at IDRIS, GENCI (Jean Zay machine, France) on grant no A0130712052 (JPP).
References

[1] M. P. Allen and D. J. Tildesley. *Computer Simulation of Liquids*. Oxford University Press, 2017.

[2] B. Leimkuhler and C. Matthews. *Molecular Dynamics*. Springer, 2015.

[3] B. C. Goh, J. A. Hadden, R. C. Bernardi, A. Singharoy, R. McGreevy, T. Rudack, C. K. Cassidy, and K. Schulten. Computational methodologies for real-space structural refinement of large macromolecular complexes. *Annual Review of Biophysics*, 45(1):253–278, 2016.

[4] J. Mortier, C. Rakers, M. Bermudez, M. S Murgueitio, S. Riniker, and G. Wolber. The impact of molecular dynamics on drug design: applications for the characterization of ligand–macromolecule complexes. *Drug discovery today*, 20(6):686–702, 2015.

[5] T. R. Walsh and M. R. Knecht. Biointerface structural effects on the properties and applications of bioinspired peptide-based nanomaterials. *Chemical reviews*, 117(20):12641–12704, 2017.

[6] M. Blazhynska, L. Lagardère, C. Liu, O. Adjoua, P. Ren, and J.-P. Piquemal. Water–glycan interactions drive the SARS-CoV-2 spike dynamics: insights into glycan-gate control and camouflage mechanisms. *Chem. Sci.*, 15:14177–14187, 2024.

[7] P. Freddolino, C. Harrison, Y. Liu, and K. Schulten. Challenges in protein folding simulations: Timescale, representation, and analysis. *Nature physics*, 6:751–758, 2010.

[8] J. Jung, W. Nishima, Ma. Daniels, G. Bascom, C. Kobayashi, A. Adedoyin, M. Wall, A. Lappala, D. Phillips, W. Fischer, C.-S. Tung, T. Schlick, Y. Sugita, and K. Y. Sanbonmatsu. Scaling molecular dynamics beyond 100,000 processor cores for large-scale biophysical simulations. *Journal of Computational Chemistry*, 40(21):1919–1930, 2019.

[9] D. E. Shaw, R. O. Dror, J. K. Salmon, J. P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. Deneroff, B. Batson, K. J. Bowers, E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson, K. Lindorff-Larsen, P. Maragakis, M. A. Moraes, S. Piana, Y. Shan, and B. Towles. Millisecond-scale Molecular Dynamics Simulations on Anton. In *Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis*, SC '09, pages 39:1–39:11, New York, NY, USA, 2009. ACM.

[10] Y. Shibuta, S. Sakane, E. Miyoshi, T. Takaki, and M. Ohno. Micrometer-scale molecular dynamics simulation of microstructure formation linked with multi-phase-field simulation in same space scale. *Modelling and Simulation in Materials Science and Engineering*, 27(5):054002, may 2019.

[11] I. Yu, T. Mori, T. Ando, R. Harada, J. Jung, Y. Sugita, and M. Feig. Biomolecular interactions modulate macromolecular structure and dynamics in atomistic model of a bacterial cytoplasm. *eLife*, 5:e19274, nov 2016.

[12] J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. de Groot, H. Grubmueller, and A. D. MacKerell, Jr. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. *NATURE METHODS*, 14(1):71–73, JAN 2017.

[13] J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser, and C. Simmerling. ff14sb: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99sb. *Journal of Chemical Theory and Computation*, 11(8):3696–3713, Aug 2015.

[14] M. J. Robertson, J. Tirado-Rives, and W. L. Jorgensen. Improved Peptide and Protein Torsional Energetics with the OPLS-AA Force Field. *Journal of Chemical Theory and Computation*, 11(7):3499–3509, Jul 2015.

[15] M. M. Reif, P. H. Hünenberger, and C. Oostenbrink. New Interaction Parameters for Charged Amino Acid Side Chains in the GROMOS Force Field. *Journal of Chemical Theory and Computation*, 8(10):3705–3723, Oct 2012.

[16] P. Ren and J. W. Ponder. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. *The Journal of Physical Chemistry B*, 107(24):5933–5947, 2003.
[17] Sehr Naseem-Khan, Louis Lagardère, Christophe Narth, G. Andrés Cisneros, Pengyu Ren, Nohad Gresh, and Jean-Philip Piquemal. Development of the quantum-inspired sibfa many-body polarizable force field: Enabling condensed-phase molecular dynamics simulations. *Journal of Chemical Theory and Computation*, 18(6):3607-3621, 2022. PMID: 35575306.

[18] J. W. Ponder and D. A. Case. Force Fields for Protein Simulations. In *Protein Simulations*, volume 66 of *Advances in Protein Chemistry*, pages 27–85. Academic Press, 2003.

[19] Josef Melcr and Jean-Philip Piquemal. Accurate biomolecular simulations account for electronic polarization. *Frontiers in Molecular Biosciences*, 6:143, 2019.

[20] Zhifeng Jing, Chengwen Liu, Sara Y. Cheng, Rui Qi, Brandon D. Walker, Jean-Philip Piquemal, and Pengyu Ren. Polarizable force fields for biomolecular simulations: Recent advances and applications. *Annual Review of Biophysics*, 48(1):371–394, 2019. PMID: 30916997.

[21] O. Adjoua, L. Lagardère, L.-H. Jolly, A. Durocher, T. Very, I. Dupays, Z. Wang, T. Jaffrelot Inizan, F. Célerse, P. Ren, J. Ponder, and J.-P. Piquemal. Tinker-hp: Accelerating Molecular Dynamics Simulations of Large Complex Systems with Advanced Point Dipole Polarizable Force Fields Using GPUs and Multi-GPU Systems. *Journal of Chemical Theory and Computation*, 17:2034–2053, 03 2021.

[22] P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, et al. Openmm 7: Rapid development of high performance algorithms for molecular dynamics. *PLoS computational biology*, 13(7):e1005659, 2017.

[23] C. Kutzner, S. Páll, M. Fechner, A. Esztermann, B. L. de Groot, and H. Grubmüller. Best bang for your buck: GPU nodes for GROMACS biomolecular simulations, 2015.

[24] R. Salomon-Ferrer, A. W. Gotz, D. Poole, S. Le Grand, and R. C. Walker. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. *Journal of chemical theory and computation*, 9(9):3878–3888, 2013.

[25] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. Dac Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton. Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Computer Physics Communications*, 271:108171, 2022.

[26] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, and K. Schulten. Scalable molecular dynamics with NAMD. *Journal of Computational Chemistry*, 26(16):1781–1802, 2005.

[27] M. J. Abraham, T. Murtola, R. Schulz, S. Pall, J. C. Smith, B. Hess, and E. Lindahl. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX*, 1-2:19 – 25, 2015.

[28] R. Salomon-Ferrer, D. A. Case, and R. C. Walker. An overview of the Amber biomolecular simulation package. *WIREs Computational Molecular Science*, 3(2):198–210, 2013.

[29] L. Lagardère, L.-H. Jolly, F. Lipparini, F. Aviat, B. Stamm, Z. F. Jing, M. Harger, H. Torabifard, G. A. Cisneros, M. J. Schnieders, N. Gresh, Y. Maday, P. Y. Ren, J. W. Ponder, and J.-P. Piquemal. Tinker-hp: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields. *Chem. Sci.*, 9:956–972, 2018.

[30] A.-P. Hynninen and M. F. Crowley. New faster CHARMM molecular dynamics engine. *Journal of Computational Chemistry*, 35(5):406–413, 2014.

[31] C. Kobayashi, J. Jung, Y. Matsunaga, T. Mori, T. Ando, K. Tamura, M. Kamiya, and Y. Sugita. Gen- esis 1.1: A hybrid-parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms. *Journal of Computational Chemistry*, 38(25):2193–2206, 2017.

[32] S. Plimpton. Fast Parallel Algorithms for Short-Range Molecular Dynamics. *Journal of Computational Physics*, 117(1):1–19, 1995.
[33] D.A. Gibson and E.A. Carter. Time-reversible multiple time scale ab initio molecular dynamics. *J. Phys. Chem.*, pages 13429–13434, 1993.

[34] L. Lagardère, F. Aviat, and J.-P. Piquemal. Pushing the limits of Multiple-Timestep Strategies for Polarizable Point Dipole Molecular Dynamics. *Journal of Physical Chemistry Letters*, 10(10):2593–2599, 2019.

[35] M.E. Tuckerman, B.J Berne, and Rossi A. Molecular dynamics algorithm for multiple time scales: Systems with disparate masses. *J. Chem. Phys.*, 1991.

[36] M. Tuckerman, B. J. Berne, and G. J. Martyna. Reversible multiple time scale molecular dynamics. *J. Chem. Phys.*, 97(3):1990–2001, 1992.

[37] N. Bou-Rabee. Time Integrators for Molecular Dynamics. *Entropy*, 16:138–162, 2014.

[38] B. Leimkuhler and C. Matthews. Robust and efficient configurational molecular sampling via Langevin dynamics. *The Journal of Chemical Physics*, 138(17):174102, 2013.

[39] B. Leimkuhler, C. Matthews, and G. Stoltz. The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics. *IMA J. Numer. Anal.*, 36(1):13–79, 2016.

[40] C. Matthews and B. Leimkuhler. Rational Construction of Stochastic Numerical Methods for Molecular Sampling. *Applied Mathematics Research eXpress*, 2013(1):34–56, 06 2012.

[41] J. J. Biesiadecki and R. D. Skeel. Dangers of Multiple Time Step Methods. *Journal of Computational Physics*, 109(2):318–328, 1993.

[42] Q. Ma, J. Izaguirre, and R. Skeel. Verlet-i/r-respa/impulse is Limited by Nonlinear Instabilities. *SIAM Journal on Scientific Computing*, 24(6):1951–1973, 2003.

[43] Joseph A. Morrone, Thomas E. Markland, Michele Ceriotti, and B. J. Berne. Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances. *The Journal of Chemical Physics*, 134(1):014103, 2011.

[44] A. Albaugh, M. E. Tuckerman, and T. Head-Gordon. Combining Iteration-Free Polarization with Large Time Step Stochastic-Isokinetic Integration. *Journal of Chemical Theory and Computation*, 15(4):2195–2205, 2019.

[45] B. Leimkuhler, D. T. Margul, and M. E. Tuckerman. Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps. *Molecular Physics*, 111(22-23):3579–3594, 2013.

[46] D. T. Margul and M. E. Tuckerman. A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps. *Journal of Chemical Theory and Computation*, 12(5):2170–2180, 2016.

[47] K. A. Feenstra, B. Hess, and H. J. C. Berendsen. Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich systems. *Journal of Computational Chemistry*, 20(8):786–798, 1999.

[48] P. Monmarché, J. Weisman, L. Lagardère, and J.-P. Piquemal. Velocity jump processes: An alternative to multi-timestep methods for faster and accurate molecular dynamics simulations. *The Journal of Chemical Physics*, 153(2), 07 2020. 024101.

[49] N. Gouraud, L. Journel, and P. Monmarché. The velocity jump Langevin process and its splitting scheme: long time convergence and numerical accuracy. *arXiv e-prints*, 2024.

[50] P. Monmarché, M. Rousset, and P.-A. Zitt. Exact targeting of Gibbs distributions using velocity-jump processes. *Stochastics and Partial Differential Equations: Analysis and Computations*, 2022.

[51] P. Monmarché. Piecewise deterministic simulated annealing. *ALEA Lat. Am. J. Probab. Math. Stat.*, 13(1):357–398, 2016.

[52] E. A. J. F. Peters and G. de With. Rejection-free Monte Carlo sampling for general potentials. *Phys. Rev. E* 85, 026703, 2012.

[53] P. Vanetti, A. Bouchard-Côté, G. Deligiannidis, and A. Doucet. Piecewise Deterministic Markov Chain Monte Carlo. *ArXiv e-prints*, July 2017.
[54] A. Bertazzi, J. Bierkens, and P. Dobson. Approximations of Piecewise Deterministic Markov Processes and their convergence properties. *Stochastic Processes and their Applications*, 154:91–153, 2022.

[55] A. Bertazzi, P. Dobson, and P. Monmarché. Splitting schemes for second order approximations of piecewise-deterministic Markov processes. *arXiv e-prints*, 2023.

[56] J. Bierkens, P. Fearnhead, and G. Roberts. The Zig-Zag process and super-efficient sampling for Bayesian analysis of big data. *The Annals of Statistics*, 47(3):1288–1320, 2019.

[57] A. Corbella, S. E. F. Spencer, and G. O. Roberts. Automatic Zig-Zag sampling in practice. *Statistics and Computing*, 32(6):107, Nov 2022.

[58] F. Pagani, A. Chevallier, S. Power, T. House, and S. Cotter. Nuzz: Numerical Zig-Zag for general models. *Statistics and Computing*, 34(1):61, Jan 2024.

[59] A. Bouchard-Côté, S. J. Vollmer, and A. Doucet. The bouncy particle sampler: a nonreversible rejection-free Markov chain Monte Carlo method. *J. Am. Stat. Assoc.*, 113(522):855–867, 2018.

[60] V. Lemaire, M. Thieullen, and N. Thomas. Exact Simulation of the Jump Times of a Class of Piecewise Deterministic Markov Processes. *Journal of Scientific Computing*, 75(3):1776–1807, 2017.

[61] P. A. W. Lewis and G. S. Shedler. Simulation of nonhomogeneous poisson processes by thinning. *Naval Research Logistics Quarterly*, 26(3):403–413, September 1979.

[62] J. Wang, P. Cieplak, and P. A. Kollman. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? *Journal of computational chemistry*, 21(12):1049–1074, 2000.

[63] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen. A smooth particle mesh Ewald method. *The Journal of Chemical Physics*, 103(19):8577–8593, 11 1995.

[64] J. A. Morrone, R. Zhou, and B. J. Berne. Molecular Dynamics with Multiple Time Scales: How to Avoid Pitfalls. *Journal of Chemical Theory and Computation*, 6(6):1798–1804, 2010.

[65] Tinker-HP's github. https://github.com/TinkerTools/Tinker-HP. Accessed:2024-24-09.

[66] Tinker-HP's website. https://tinker-hp.org/. Accessed: 2024-24-09.

[67] C. H. Bennett. Efficient estimation of free energy differences from Monte Carlo data. *Journal of Computational Physics*, 22(2):245–268, 1976.

[68] Louis Lagardère, Lise Maurin, Olivier Adjoua, Krystel El Hage, Pierre Monmarché, Jean-Philip Piquemal, and Jérôme Hénin. Lambda-abf: Simplified, portable, accurate, and cost-effective alchemical free-energy computation. *Journal of Chemical Theory and Computation*, 2024.

[69] Thomas Plé, Nastasia Mauger, Olivier Adjoua, Théo Jaffrelot Inizan, Louis Lagardère, Simon Huppert, and Jean-Philip Piquemal. Routine molecular dynamics simulations including nuclear quantum effects: From force fields to machine learning potentials. *Journal of Chemical Theory and Computation*, 19(5):1432–1445, 2023.

[70] Nohad Gresh, G. Andrés Cisneros, Thomas A. Darden, and Jean-Philip Piquemal. Anisotropic, polarizable molecular mechanics studies of inter- and intramolecular interactions and ligand- macromolecule complexes. a bottom-up strategy. *Journal of Chemical Theory and Computation*, 3(6):1960–1986, 2007. PMID: 18978934.

[71] Yue Shi, Pengyu Ren, Michael Schnieders, and Jean-Philip Piquemal. *Polarizable Force Fields for Biomolecular Modeling*, chapter 2, pages 51–86. John Wiley & Sons, Ltd, 2015.

[72] Igor Chollet, Louis Lagardère, and Jean-Philip Piquemal. Ankh: A generalized o(n) interpolated ewald strategy for molecular dynamics simulations. *Journal of Chemical Theory and Computation*, 19(10):2887–2905, 2023. PMID: 37134146.