Time Integrators for Molecular Dynamics

Abstract

This paper invites the reader to learn more about time integrators for Molecular Dynamics simulation through a simple MATLAB implementation. An overview of methods is provided from an algorithmic viewpoint that emphasizes long-time stability and finite-time dynamic accuracy. The given software simulates Langevin dynamics using an explicit, second-order (weakly) accurate integrator that exactly reproduces the Boltzmann-Gibbs density. This latter feature comes from adding a Metropolis acceptance-rejection step to the integrator. The paper discusses in detail the properties of the integrator. Since these properties do not rely on a specific form of a heat or pressure bath model, the given algorithm can be used to simulate other bath models including, e.g., the widely used v-rescale thermostat.

1. Introduction

Molecular Dynamics (MD) simulation refers to the time integration of Hamilton’s equations often coupled to a heat or pressure bath [1–5]. From its early use in computing equilibrium dynamics of homogeneous molecular systems [6–13] and pico- to nano-scale protein dynamics [14–23], the method has evolved into a general purpose tool for simulating statistical properties of heterogeneous molecular systems [24]. Accessible time horizons have increased remarkably: the time line in Figure 1 attempts to capture this nearly billion-fold improvement in capability over the last forty or so years. To put this speedup in perspective, though, computing power has increased by about eight powers of ten over this time period as predicted by Moore’s law.

To be clear, the selection of applications and methods shown in Figure 1 is not comprehensive and heavily biased towards the specific ideas and methods that inform this paper. The applications highlighted are simulations of liquid argon [6], water [11], protein dynamics without solvent [14,15] and biopolymer dynamics with solvent [25–31]. The methods include the following “upgrades” to MD simulation: Verlet integrator and neighbor lists [7], cell linked list [32], the SHAKE integrator for constraints [33], stochastic heat baths via Langevin dynamics [34,35], a library of empirical potentials [36], a deterministic heat bath via Nosé-Hoover dynamics [37,38], the fast multipole method [39], multiple time steps [40], splitting methods for Langevin dynamics [41–43], quasi-symplectic integrators [44,45], (fast) combined neighbor and cell lists [46], the v-rescale thermostat [47] and the stochastic Nosé-Hoover Langevin thermostat [48–50].

Near future applications of MD simulation include micro- to milli-scale simulations of biomolecular processes, like protein folding, ligand binding, membrane transport and biopolymer conformational changes [51–53]. In addition, atomistic MD simulations are used more sparingly in multiscale models [54–58] and rare event simulation, such as the finite temperature string method and milestoning [59–62]. Given this continuous development and generalization of MD, it is not a stretch to suppose that MD will play a transformative role in medicine, technology and education in the twenty-first century.

In its standard form, the method inputs a random initial condition, physical and numerical parameters and outputs a long discrete path of the molecular system. Statistical quantities, like velocity correlation or mean radius of gyration, are usually computed online, i.e., as points along this trajectory are produced. MD simulation is built atop a cheap forward Euler-like integrator that requires only a single interactomic force field evaluation per step. Even though MD seems straightforward, software implementations of MD are typically optimized for performance [36,63,64], and as a side effect, make it cumbersome for non-experts to learn and modify.

Also, besides this issue, due to the interplay between stochastic Brownian and molecular forces, infinitely long trajectories of existing MD integrators do not have the right distribution. What happens is that the Brownian force can cause the integrator to enter regions where its approximation to the molecular force is inaccurate and possibly destabilizing. In the latter case, the approximation spends a disproportionate amount of time at higher energies, and thus, the invariant measure of the approximation, if it even exists, is not correct. This phenomenon is a well-known shortcoming of explicit integrators for nonlinear diffusions [65–69].

Recently, a probabilistic approach was proposed to solve this problem, which questions the notion that Monte Carlo methods and MD have different aims: the former strictly samples probability distributions, and the latter estimates dynamics. The basic idea is to combine a standard MD integrator with a Metropolis-Hastings algorithm targeted to the Boltzmann-Gibbs distribution [70–72]. Because the scheme is a Monte Carlo method, it exactly preserves the desired distribution [71,72]. This property implies numerical stability over long-time simulations. However, the price to be paid for this stability is a loss of accuracy whenever a move is rejected and some overhead in evaluating the Metropolis acceptance-rejection step. Still, a Metropolized integrator is dynamically accurate on finite-time intervals [72,73], and so, even though a Metropolized integrator involves a Monte Carlo step, its aim and philosophy are very different from Monte Carlo methods, whose only goal is to sample a target distribution with no concern for the dynamics [71,74–82]. In principle, this approach offers a simple alternative to costly implicit integrators, but are Metropolized integrators ready for daily use in MD simulation? The answer to this question is unclear, since this approach is new and has not been tested on enough examples.

Motivated by these issues, this paper builds a software system for MD simulation with a Metropolis step built in and applies it to a homogeneous molecular system. The algorithm and its properties are introduced in a step-by-step fashion. In particular, we show that the integrator is second-order weakly accurate on finite-time intervals and converges to the Boltzmann-Gibbs distribution in the long-time limit. The software version of the algorithm is written in the latest version of MATLAB with plenty of comments, variables that are descriptively named and operations that can be easily translated into mathematical expressions [83]. Since MATLAB is widely available, this design ensures that the software will be easy-to-use and cross-platform. The following MATLAB-specific file formats will be used.

• (F1) MATLAB script and function files are written in the MATLAB language and can be run from the MATLAB command line without ever compiling them.

• (F2) MATLAB executable (MEX) files are written in the “C” language and compiled using the MATLAB mex function. The resulting executable is comparable in efficiency to a “C” code and can be called directly from the MATLAB command line. We will use MEX-files for performance-critical routines [84].

• (F3) MATLAB binary (MAT) files will be used to store simulation data.

The paper is organized as follows. We begin with an overview of integrators that have been proposed in MD simulation in Section 2. We explain how to Metropolize each of these schemes to make them long-time stable in Section 3, and as an application, we use a Metropolized scheme to generate a long trajectory of a Lennard-Jones fluid in Section 4. Generalizations of corrected MD integrators to other molecular models are discussed in Section 5. The paper closes by discussing some potential pitfalls in high dimension and tricks to get the integrator to scale well in Section 6.

2. Algorithmic Introduction to Time Integrators for MD Simulation

For pedagogical reasons, we will start with Langevin dynamics of a system of N molecules. Then, we show in Section 5 how to simulate more general models of molecular systems. Denote by m_j > 0 and q_j the mass and position of the j-th molecule, respectively. The governing Langevin equation is given by:

$$\{\begin{array}{l}\frac{d{\mathit{q}}_j}{dt}(t)=m_j^{-1}{\mathit{p}}_j(t),\\ d{\mathit{p}}_j(t)=-\frac{\partial U}{\partial {\mathit{q}}_j}(\mathit{q}(t))dt-\gamma {\mathit{p}}_j(t)dt+\sqrt{2kT\gamma m_j}d{\mathit{w}}_j,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}j=1,\cdots ,N$$

(1)

where q = (q₁, · · ·, q_N) and p = (p₁, · · ·, p_N) denote the positions and momenta of the particles, kT is the temperature factor, and ${\left\{{\mathit{w}}_j\right\}}_{j=1}^N$ are N-independent Brownian motions. The last two terms in the second equation in (1) represent the effect of a heat bath with parameter γ. In Langevin dynamics, positions are differentiable, and due to the irregularity of the Brownian force, momenta are just continuous, but not differentiable. This difference in regularity explains why the first equation in (1) is written as an ordinary differential equation (ODE) and the second equation is written as a stochastic differential equation (SDE).

The bath-free dynamics is a Hamiltonian system with the following Hamiltonian energy function:

$$H(\mathit{q},\mathit{p})=\sum _{j=1}^N\frac{1}{2m_j}{\left|{\mathit{p}}_j\right|}^2+U(\mathit{q})$$

(2)

Since the masses are constant, this Hamiltonian nicely separates into a kinetic and potential energy that are purely functions of p and q, respectively. The stationary probability density of the solution to Equation (1) is the Boltzmann-Gibbs density given by:

$$\nu (\mathit{q},\mathit{p})=Z^{-1}\text{exp}\left(-\frac{1}{kT}H(\mathit{q},\mathit{p})\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Z=\int \text{exp}\left(-\frac{1}{kT}H(\mathit{q},\mathit{p})\right)\hspace{0.17em}\hspace{0.17em}d\mathit{q}d\mathit{p}$$

(3)

Let h be a given time step size and m = diag(m₁, · · ·, m_N). Let (Q₀, P₀) denote the position and momentum of the molecular system at time t > 0. The simplest approximation to Equation (1) is a forward Euler discretization or Euler-Maruyama scheme [85] that computes an updated position and momentum (Q₁, P₁) at t + h using:

$$\begin{array}{l}{\mathit{Q}}_1={\mathit{Q}}_0+h{\mathit{m}}^{-1}{\mathit{P}}_0\\ {\mathit{P}}_1={\mathit{P}}_0-h\nabla U{\mathit{Q}}_0-h\gamma {\mathit{P}}_0+\sqrt{h}\sqrt{2kT\gamma }{\mathit{m}}^{1/2}\mathit{\xi }\end{array}$$

(forward Euler)

Here, ξ ∈ ℝⁿ denotes a Gaussian random vector with mean zero and covariance 𝔼(ξ_iξ_j) = δ_ij. The problem with this approximation is that the forward Euler method is known to diverge in finite-time when the derivatives of the potential are unbounded, which is the norm in MD simulation. The precise statement and proof of divergence in a general setting can be found in [86]. By far the most computationally intensive part of the time-stepping algorithm is the evaluation of the potential force. Thus, we will restrict our discussion to schemes that, like Euler, only require a single force field evaluation per step.

An improvement to the forward Euler method is the following two-step scheme:

$${\mathit{Q}}_2=(1+e^{-\gamma h}){\mathit{Q}}_1-e^{-\gamma h}{\mathit{Q}}_0+\frac{1-e^{-\gamma h}}{\gamma }{\mathit{m}}^{-1}\left(-h\nabla U({\mathit{Q}}_1)+\sqrt{h}\sqrt{2kT\gamma }{\mathit{m}}^{1/2}\mathit{\xi }\right)$$

(BBK)

In the limit, γ → 0, this scheme reduces to the well-known Verlet integrator for MD simulation [7]. Just like Verlet, this integrator defines a map on pairs of molecular system configurations. Substituting the approximation, e^−γh ≈ (1−γh/2)/(1+γh/2), into the above yields the Brünger-Brooks-Karplus (BBK) scheme, as appearing in [35]. Like the forward Euler method, this method is explicit and only requires one new force evaluation per step.

Second-order accurate schemes that generalize the Velocity Verlet integrator to Langevin dynamics were proposed in a sequence of papers [42–44,87,88]. Here, we mention two of these schemes that are both Strang splittings of Equation (1). The first was proposed by Ricci and Ciccotti [42] and consists of the following sub-steps:

$$\underset{\text{exactly}\hspace{0.17em}\text{evolve}\hspace{0.17em}\text{by}\hspace{0.17em}1/2\hspace{0.17em}\text{a}\hspace{0.17em}\text{step}}{\underbrace{\left(\begin{array}{c}\dot{\mathit{q}}(t)={\mathit{m}}^{-1}\mathit{p}(t)\\ d\mathit{p}(t)=0\end{array}\right)}}°\underset{\text{exactly}\hspace{0.17em}\text{evolve}\hspace{0.17em}\text{by}\hspace{0.17em}\hspace{0.17em}\text{a}\hspace{0.17em}\text{step}}{\underbrace{\left(\begin{array}{c}\dot{\mathit{q}}(t)=0\\ d\mathit{p}(t)=-\nabla U(\mathit{q}(t))dt-\gamma \mathit{p}(t)dt+\sqrt{2kT\gamma }{\mathit{m}}^{1/2}d\mathit{W}\end{array}\right)}}°\underset{\text{exactly}\hspace{0.17em}\text{evolve}\hspace{0.17em}\text{by}\hspace{0.17em}1/2\hspace{0.17em}\text{a}\hspace{0.17em}\text{step}}{\underbrace{\left(\begin{array}{c}\dot{\mathit{q}}(t)={\mathit{m}}^{-1}\mathit{p}(t)\\ d\mathit{p}(t)=0\end{array}\right)}}$$

Each step in this decomposition can be exactly solved. Clearly, the half-steps are easy to solve, since momentum is constant over each of these half-steps. The SDE appearing in the inner step can also be exactly solved, since it is linear in momentum (see Chapter 5 in [89]). This splitting is quite natural, since it treats the heat bath forces in the same way as the potential forces.

A related, but different, splitting method was proposed by Bussi and Parinello in [43] and is given by:

$$\underset{\text{exactly}\hspace{0.17em}\text{evolve}\hspace{0.17em}\text{by}\hspace{0.17em}1/2\hspace{0.17em}\text{a}\hspace{0.17em}\text{step}}{\underbrace{\left(\begin{array}{c}\dot{\mathit{q}}(t)=0\\ d\mathit{p}(t)=-\gamma \mathit{p}(t)dt+\sqrt{2kT\gamma }{\mathit{m}}^{1/2}d\mathit{W}\end{array}\right)}}°\underset{\begin{array}{l}\text{approximately}\hspace{0.17em}\text{evolve}\\ \text{using}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{a}\hspace{0.17em}\text{step}\hspace{0.17em}\text{of}\hspace{0.17em}\text{Verlet}\end{array}}{\underbrace{\left(\begin{array}{c}\dot{\mathit{q}}(t)={\mathit{m}}^{-1}\mathit{p}(t)\\ \dot{\mathit{p}}(t)=-\nabla U(\mathit{q}(t))\end{array}\right)}}°\underset{\text{exactly}\hspace{0.17em}\text{evolve}\hspace{0.17em}\text{by}\hspace{0.17em}1/2\hspace{0.17em}\text{a}\hspace{0.17em}\text{step}}{\underbrace{\left(\begin{array}{c}\dot{\mathit{q}}(t)=0\\ d\mathit{p}(t)=-\gamma \mathit{p}(t)dt+\sqrt{2kT\gamma }{\mathit{m}}^{1/2}d\mathit{W}\end{array}\right)}}$$

Notice that this decomposition splits the Langevin dynamics into its Hamiltonian and heat bath parts, which makes it easy to analyze the structural properties of the scheme. A Velocity Verlet integrator is used to approximate the Hamiltonian dynamics. This approximation exactly preserves phase space volume and preserves energy to third-order accuracy per step. Moreover, the solution to the SDE appearing in the half-steps exactly preserves the Boltzmann-Gibbs density.

Since the Velocity Verlet integrator does not exactly preserve energy, the composition above does not exactly preserve the stationary distribution with density in Equation (3). In [90], it was shown that if the derivatives of the potential are all bounded, the Bussi and Parinello integrator possesses an invariant measure that is $𝒪$(h²) close to the Boltzmann-Gibbs distribution. In this same context, the leading order error term in the integrator’s approximation to the invariant measure was explicitly determined [91]. Technically speaking, however, these results do not directly apply to MD simulation, since real MD simulation involves potentials whose derivatives are unbounded, e.g., Lennard-Jones forces. As a consequence of this irregularity in the force fields and discretization error, explicit schemes, like this one, may either not detect features of the potential energy properly, which leads to unnoticed, but large errors in dynamic quantities such as the mean first passage time, or may mishandle soft- or hard-core potentials, which leads to numerical instabilities; see the numerical examples in [92]. These numerical artifacts motivate adding a Metropolis accept/refusal sub-step to the integrator. In the next section, we show how to Metropolize all of the MD integrators presented in this section. In Section 5, we explain how to generalize the Metropolis-corrected Bussi and Parinello algorithm to a larger class of diffusion processes.

3. Metropolis-Corrected MD Integrators

Here, we show how to add a Metropolis acceptance-rejection step to a BBK-type scheme and the Bussi and Parinello splitting scheme and then precisely state the properties of these integrators. We start with a detailed description of each algorithm. Both algorithms require evaluating the acceptance probability given by the usual Metropolis ratio:

$$\alpha (\mathit{q},\mathit{p},\mathit{Q},\mathit{P})=\text{min}\left(1,\text{exp}\left(-\frac{1}{kT}(H(\mathit{Q},\mathit{P})-H(\mathit{q},\mathit{p},))\right)\right)$$

(4)

The procedure to Metropolize the Ricci and Ciccotti scheme can be found in Section 2 of [70].

Algorithm 3.1 (First-order BBK-type integrator). Given the current state (Q₀, P₀) at time t, the algorithm proposes a new state ( ${\mathit{Q}}_1^{★}$, ${\mathit{P}}_1^{★}$) at time t + h for some time step h > 0 via:

$$\left(\begin{array}{c}{\mathit{Q}}_1^{★}\\ {\mathit{P}}_1^{★}\end{array}\right)=\left(\begin{array}{c}{\mathit{Q}}_0+{\mathit{m}}^{-1}\left(h{\mathit{P}}_0-\frac{h^2}{2}\nabla U({\mathit{Q}}_0)\right)\\ {\mathit{P}}_0-\frac{h}{2}\left(\nabla U({\mathit{Q}}_0)+\nabla U({\mathit{Q}}_1^{★})\right)\end{array}\right)$$

(Step 1)

This “proposal move” ( ${\mathit{Q}}_1^{★}$, ${\mathit{P}}_1^{★}$) is then accepted or rejected:

$$\left(\begin{array}{c}{\tilde{\mathit{Q}}}_1\\ {\tilde{\mathit{P}}}_1\end{array}\right)=x\left(\begin{array}{c}{\mathit{Q}}_1^{★}\\ {\mathit{P}}_1^{★}\end{array}\right)+(1-x)\hspace{0.17em}\hspace{0.17em}\left(\begin{array}{c}{\mathit{Q}}_0\\ -{\mathit{P}}_0\end{array}\right)$$

(Step 2)

where x is a Bernoulli random variable with parameter α(Q₀, P₀, ${\mathit{Q}}_1^{★}$, ${\mathit{P}}_1^{★}$) given by Equation (4). The actual update of the system is taken to be:

$$\left(\begin{array}{c}{\mathit{Q}}_1\\ {\mathit{P}}_1\end{array}\right)=\left(\begin{array}{c}{\tilde{\mathit{Q}}}_1\\ \text{exp}(-\gamma h){\tilde{\mathit{P}}}_1+\sqrt{kT}\sqrt{1-\text{exp}(-2\gamma h)}{\mathit{m}}^{1/2}\mathit{\xi }\end{array}\right)$$

(Step 3)

Here, ξ ∈ ℝⁿ denotes a Gaussian random vector with mean zero and covariance 𝔼(ξ_iξ_j) = kTδ_ij.

The momenta of the molecules gets reversed if a move is rejected in Step 2 of Algorithm 3.1. This momentum flip is necessary for the algorithm to preserve the correct stationary distribution [70,71], but results in an O(1) error in dynamics. High acceptance rates are therefore needed to ensure that the time lag between successive rejections is frequently long enough for the approximation to capture the desired dynamics. Since the acceptance rate in Equation (4) is related to how well the Verlet integrator in (Step 1) preserves energy after a single step, this rejection rate is O(h³). Thus, in practice, we find that the time step required to obtain a sufficiently high acceptance rate is often automatically fulfilled by a time step that sufficiently resolves the desired dynamics. Each step of this algorithm requires: evaluating the atomic force field once in the third equation of (Step 1), generating a Bernoulli random variable with parameter α in (Step 2) and generating an n-dimensional Gaussian vector in (Step 3). We stress that (Step 2) in Algorithm 3.1 is all that is needed to get MD integrators to exactly preserve the Boltzmann-Gibbs density in Equation (3).

Next, we show how to Metropolize the Bussi and Parinello splitting integrator.

Algorithm 3.2 (Second-order Bussi and Parinello integrator). Let ξ, η ∈ ℝⁿ be two independent Gaussian random vectors with mean zero and covariance 𝔼(ξ_iξ_j) = 𝔼 (η_iη_j) = δ_ij. Given a time step size h and the current state (Q₀, P₀) at time t, the algorithm takes a half-step of the heat bath dynamics:

$$\left(\begin{array}{c}{\tilde{\mathit{Q}}}_0\\ {\tilde{\mathit{P}}}_0\end{array}\right)=\left(\begin{array}{c}{\mathit{Q}}_0\\ \text{exp}(-\gamma h/2){\tilde{\mathit{P}}}_0+\sqrt{kT}\sqrt{1-\text{exp}(-\gamma h)}{\mathit{m}}^{1/2}\mathit{\xi }\end{array}\right)$$

(Step 1)

Followed by a full step of Verlet to compute a proposal move ( ${\tilde{\mathit{Q}}}_1^{★}$, ${\tilde{\mathit{P}}}_1^{★}$):

$$\left(\begin{array}{c}{\tilde{\mathit{Q}}}_1^{★}\\ {\tilde{\mathit{P}}}_1^{★}\end{array}\right)=\left(\begin{array}{c}{\tilde{\mathit{Q}}}_0+{\mathit{m}}^{-1}\left(h{\tilde{\mathit{P}}}_0-\frac{h^2}{2}\nabla U({\tilde{\mathit{Q}}}_0)\right)\\ {\mathit{P}}_0-\frac{h}{2}\left(\nabla U({\tilde{\mathit{Q}}}_0)+\nabla U({\tilde{\mathit{Q}}}_1^{★})\right)\end{array}\right)$$

(Step 2)

This proposal move ( ${\tilde{\mathit{Q}}}_1^{★}$, ${\tilde{\mathit{P}}}_1^{★}$) is then accepted or rejected:

$$\left(\begin{array}{c}{\tilde{\mathit{Q}}}_1\\ {\tilde{\mathit{P}}}_1\end{array}\right)=x\left(\begin{array}{c}{\tilde{\mathit{Q}}}_1^{★}\\ {\tilde{\mathit{P}}}_1^{★}\end{array}\right)+(1-x)\hspace{0.17em}\hspace{0.17em}\left(\begin{array}{c}{\tilde{\mathit{Q}}}_0\\ -{\tilde{\mathit{P}}}_0\end{array}\right)$$

(Step 3)

where x is a Bernoulli random variable with parameter α(Q̃₀, P̃₀, ${\tilde{\mathit{Q}}}_1^{★}$, ${\tilde{\mathit{P}}}_1^{★}$) given by Equation (4). The actual update of the system t + h is taken to be:

$$\left(\begin{array}{c}{\mathit{Q}}_1\\ {\mathit{P}}_1\end{array}\right)=\left(\begin{array}{c}{\tilde{\mathit{Q}}}_1\\ \text{exp}(-\gamma h/2){\tilde{\mathit{P}}}_1+\sqrt{kT}\sqrt{1-\text{exp}(-\gamma h)}{\mathit{m}}^{1/2}\mathit{\eta }\end{array}\right)$$

(Step 4)

This algorithm requires generating two independent n-dimensional Gaussian vectors per step. Thus, it is more costly than Algorithm 3.1. However, the advantage of doing this is that the resulting Metropolis corrected algorithm is second-order weakly accurate, as the following Proposition states.

Proposition 3.3. Let (Q_n, P_n) represent the numerical approximation produced by Algorithm 3.2 at time nh with the same initial condition as the true solution: (Q₀, P₀) = (q(0), p(0)). For every time interval T > 0 and for suitable observables f(q, p), there exists a C(T) > 0, such that:

$$\left|𝔼f(\mathit{q}(\lfloor t/h\rfloor h),\hspace{0.17em}\hspace{0.17em}\mathit{p}(\lfloor t/h\rfloor h))-𝔼f\left({\mathit{Q}}_{\lfloor t/h\rfloor },{\mathit{P}}_{\lfloor t/h\rfloor }\right)\right|\le C(T)h^2$$

(5)

for all t < T.

This accuracy concept is sufficient for computing means and correlation functions at finite-time and equilibrium correlations. Figure 2 verifies this Proposition by checking the weak accuracy of Algorithms 3.1 and 3.2 on a harmonic oscillator test problem.

To be specific, Figure 2 plots the weak accuracy of the Metropolis-corrected MD integrators with respect to the true solution of the Langevin dynamics of a harmonic oscillator: q̇(t) = p(t), $dp(t)=-q(t)-p(t)+\sqrt{2}dw(t)$, with initial condition q(0) = 1.0, p(0) = 0. The time steps tested are h = 2⁻ⁿ, where n is given on the x-axis. The quantity monitored for the error is the estimate of 𝔼(q(1)² + p(1)²) = 1.699445410 computed analytically. The dashed and solid curves are the graphs of 2⁻ⁿ(= h) and 2⁻²ⁿ(= h²) versus n, respectively.

Proof. The desired single-step error estimate can be obtained from an application of the triangle inequality:

$$\left|𝔼f(\mathit{q}(h),\hspace{0.17em}\hspace{0.17em}\mathit{p}(h))-𝔼f\left({\mathit{Q}}_1,{\mathit{P}}_1\right)\right|\le \left|𝔼f(\mathit{q}(h),\hspace{0.17em}\hspace{0.17em}\mathit{p}(h))-𝔼f\left({\hat{\mathit{Q}}}_1,{\hat{\mathit{P}}}_1\right)\right|+\left|𝔼f\left({\hat{\mathit{Q}}}_1,{\hat{\mathit{P}}}_1\right)-𝔼f\left({\mathit{Q}}_1,{\mathit{P}}_1\right)\right|$$

(6)

where (Q̂₁, P̂₁) denotes one step of the uncorrected Bussi and Parinello scheme with (Q̂₀, P̂₀) = (q(0), p(0)). The first term in the upper bound in Equation (6) is O(h³), since the unadjusted scheme is a Strang splitting of Equation (1). To bound the second term in Equation (6), note that:

$$𝔼f\left({\mathit{Q}}_1,{\mathit{P}}_1\right)-𝔼f\left({\hat{\mathit{Q}}}_1,{\hat{\mathit{P}}}_1\right)=𝔼\left\{\left(\overline{f}({\tilde{\mathit{Q}}}_1^{★},{\tilde{\mathit{P}}}_1^{★})-\overline{f}\left({\tilde{\mathit{Q}}}_0,-{\tilde{\mathit{P}}}_0\right)\right)\hspace{0.17em}\hspace{0.17em}\left(\alpha \left({\tilde{\mathit{Q}}}_0,{\tilde{\mathit{P}}}_0,{\tilde{\mathit{Q}}}_1^{★},{\tilde{\mathit{P}}}_1^{★}\right)-1\right)\right\}$$

where we have introduced the auxilary function:

$$\overline{f}(\mathit{q},\mathit{p})=𝔼f(\mathit{q},\text{exp}(-\gamma h/2)\mathit{p}+\sqrt{kT}\sqrt{1-\text{exp}(-\gamma h)}{\mathit{m}}^{1/2}\mathit{\eta })$$

Since the rejection rate is $𝒪$(h³), it follows from the above expression that the second term in the upper bound of Equation (6) is also O(h³). Standard results in numerical analysis for SDEs then imply that the algorithm converges weakly on finite-time intervals with global order two; see, for instance, [93] (Chapter 2.2).

For completeness sake, we also provide a statement that both algorithms are ergodic.

Proposition 3.4. Let (Q_n, P_n) be the numerical approximation produced by Algorithms 3.1 or 3.2 at time nh. Then, for suitable observables f(q, p):

$$\underset{T\to \infty }{\text{lim}}\frac{1}{T}{\int }_0^{T}f\left({\mathit{Q}}_{\lfloor t/h\rfloor },{\mathit{P}}_{\lfloor t/h\rfloor }\right)dt\to {\int }_{{ℝ}^{2n}}f(\mathit{q},\mathit{p})\nu (\mathit{q},\mathit{p})d\mathit{q}d\mathit{p}$$

(7)

Here, ν(q, p) denotes the Boltzmann-Gibbs density defined in Equation (3).

A proof of this Proposition can be found in [72].

4. Application to Lennard-Jones Fluid

Listing 1 translates Algorithm 3.2 into the MATLAB language. Intrinsically defined MATLAB functions appear in boldface. The algorithm uses MATLAB’s built in random number generators to carry out Step 1, Step 3 and Step 4. In particular, the Bernoulli random variable, x, in Step 3 is generated in Line 20, and the Gaussian vectors in Step 1 and Step 4 are generated on Line 9 and Line 29, respectively. In addition to updating the positions and momenta of the system, the program also stores the previous value of the potential energy and force, so that the force and potential energy is evaluated in Line 15 just once per simulation step. This evaluation calls a MEX function, which inputs the current position of the molecular system and outputs the force field and potential energy at that position. We use a MEX function, because the atomistic force field evaluation cannot be easily vectorized and is, by far, the most computationally demanding step in MD. The PreProcessing script file called in Line 2 defines the physical and numerical parameters, sets the initial condition and allocates space for storing simulation data. Sample averages are updated as new points on the trajectory are produced in the UpdateSampleAverages script file invoked in Line 35. Finally, the outputs produced by the algorithm are handled by the PostProcessing script file in Line 39.

Let us consider a concrete example: a Lennard-Jones fluid that consists of N identical atoms [1–3]. The configuration space of this system is a fixed cubic box with periodic boundary conditions. The distance between the i-th and j-th particle is defined according to the minimum image convention, which states that the distance between q_i and q_j in a cubic box of length ℓ is:

$$d_{MD}({\mathit{q}}_i,{\mathit{q}}_j)\stackrel{\text{def}}{=}\left|({\mathit{q}}_i,{\mathit{q}}_j)-\ell \lfloor ({\mathit{q}}_i,{\mathit{q}}_j)/\ell \rceil \right|$$

(8)

where ⌊·⌉ is the nearest integer function. In terms of this distance, the total potential energy is a sum over all pairs:

$$U(\mathit{q})=\sum _{i=1}^{n-1}\sum _{j=i+1}^n{U}_{LJ}\left(d_{MD}({\mathit{q}}_i,{\mathit{q}}_j)\right)$$

(9)

where U_LJ(r) is the following truncated Lennard-Jones potential function:

$$U_{LJ}(r)=\{\begin{array}{ll}f(r)-f(r_c),\hfill & r<r_c\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$$

(10)

Listing 1. Metropolized MD Integrator: MDintegrator.m

Listing 1. Metropolized MD Integrator: MDintegrator.m

Listing 2. Metropolized MD Integrator: PreProcessing.m

Listing 2. Metropolized MD Integrator: PreProcessing.m

Here, f(r) = 4(1/r¹² − 1/r⁶) and r_c is the cutoff radius, which is bounded above by the size of the simulation box; and we have used dimensionless units to describe this system, where energy is rescaled by the depth of the Lennard-Jones potential energy and length by the point where the potential energy is zero. The error introduced by the truncation in Equation (10) is proportional to the density of the molecular system and can be made arbitrarily small by selecting the cutoff distance to be sufficiently large. A direct evaluation of the potential force, ∇U(q), scales like O(N²), and typically dominates the total computational cost. In practice, neighbor/cell lists, also called Verlet lists, are used in order to obtain a force evaluation that scales linearly with system size. Since the system we consider will have just a few hundred atoms, there is, however, little advantage to using these data structures, or using a fast force field evaluation, and thus, ForceFieldmex evaluates the force and energy using a sum over all particle pairs.

Listing 2 shows the PreProcessing script, which sets the parameters provided in

Table 1. Simulation parameters.

	Parameter	Description	Value
Physical Parameters	ρ	density	{0.6, 0.7, 0.8, 0.9, 1.0, 1.1}
	kT	temperature factor	0.5
	γ	heat bath parameter	0.01
	Nm	# of molecules	512
	T	time-span for autocorrelation	2
Numerical Parameters	h	time step	0.005
	Ns	# of simulation steps	105
	rc	Lennard-Jones force cutoff radius	21/6

and constructs the initial condition, where the N atoms are assumed to be at rest and on the sites of a face-centered cubic lattice. The command, rng(123), on Line 3 sets the seed of the random number generator functions, RAND and RANDN. The acceptance rates at every step and the velocity autocorrelation are updated in the UpdateSampleAverages script shown in Listing 3. The mean acceptance rate, which is outputted in the PostProcessing script shown in Listing 4, must be high enough to ensure that the dynamics is accurately represented. To compute the autocorrelation of an observable over a time interval of length T, the value of that observable along the entire trajectory is not needed. In fact, it suffices to use the values of this observable along a piece of trajectory over a moving time-window [t_i, t_i + T ], where t_i = i × h. This storage space is allocated in PreProcessing and is updated in UpdateSampleAverages. More precisely, the molecular velocities are stored in the pivot array from i − N_a to i, where i is the index of the current position and N_a = ⌈T/h⌉ + 1. Notice that velocity autocorrelations are not computed until after the index, i, exceeds 10⁴. This equilibration time removes some of the statistical bias that may arise from using a non-random initial condition. Short-time trajectories of this molecular system are plotted in Figure 3 from an initial condition where atoms are placed on the sites of a face-centered cubic lattice and at rest. The trajectory is computed using the numerical and physical parameters indicated in Table 1, with the exception of the number of steps, which is set equal to N_s = 1000. Notice that at lower densities particle trajectories are more diffusive and less localized. Using the parameters provided in Table 1, we compute velocity autocorrelations for a range of density values in Figure 4. Since the heat bath parameter is set to a small value, these figures are in qualitative agreement with those obtained by simulating the molecular system with no heat bath as shown in Figure 5.2 of [3].

Listing 3. Metropolized MD Integrator: UpdateSampleAverages.m

Listing 3. Metropolized MD Integrator: UpdateSampleAverages.m

Listing 4. Metropolized MD Integrator: PostProcessing.m

Listing 4. Metropolized MD Integrator: PostProcessing.m

5. General Case

Here, we show how the preceding ideas extend to other molecular systems that obey stochastic differential equations. In the process, we generalize the Metropolized Bussi and Parinello integrator (Algorithm 3.2) to a big class of diffusion processes, including the v-rescale thermostat. We begin with the underlying Hamiltonian dynamics of a molecular system.

5.1. Bath-Free Dynamics

MD is based on Hamilton’s equations for a Hamiltonian H: ℝ^2d → ℝ:

$$\dot{\mathit{z}}(t)=\mathit{J}\nabla H(\mathit{z}(t)),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{z}(0)\in {ℝ}^{2d}$$

(11)

where z(t) = (q(t), p(t)) is a vector of molecular positions q(t) ∈ ℝ^d and momenta p(t) ∈ ℝ^d and J is the 2d × 2d skew-symmetric matrix defined as:

$$\mathit{J}=\left(\begin{array}{cc}{\mathbf{0}}_{d\times d}& {\mathit{I}}_{d\times d}\\ -{\mathit{I}}_{d\times d}& {\mathbf{0}}_{d\times d}\end{array}\right)$$

(12)

The Hamiltonian, H(z), represents the total energy of the molecular system and is typically “separable”, meaning that it can be written as:

$$H(\mathit{z})=K(\mathit{p})+U(\mathit{q}),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{z}=(\mathit{q},\mathit{p})$$

(13)

where K(p) and U(q) are the kinetic and potential energy functions, respectively [94]. In MD, the kinetic energy function is a positive definite quadratic form, and the potential energy function involves “fudge factors” determined from experimental or quantum mechanical studies of pieces of the molecular system of interest [36]. The accuracy of the resulting energy function must be systematically verified by comparing MD simulation data to experimental data [95]. The flow that Equation (11) determines has the following structure:

• (S1) volume-preserving (since the vector-field in Equation (11) is divergenceless); and

• (S2) energy-preserving (since J is skew-symmetric and constant).

Explicit symplectic integrators, like the Verlet scheme, exploit these properties to obtain long-time stable schemes for Hamilton’s equations [96,97].

5.2. Governing Stochastic Dynamics

In order to mimic experimental conditions, Equation (11) is often coupled to a bath that puts the system at constant temperature and/or pressure. The standard way to do this is to assume that the system with a bath is governed by a stochastic ordinary differential equation (SDE) of the type:

$$d\mathit{Y}(t)=\underset{\text{deterministic}\hspace{0.17em}\text{drift}}{\underbrace{\mathit{A}(\mathit{Y}(t))dt}}+\underset{\text{heat}\hspace{0.17em}\text{bath}}{\underbrace{(\text{div}\hspace{0.17em}\mathit{D})(\mathit{Y}(t))dt+\sqrt{2kT}\mathit{B}(\mathit{Y}(t))d\mathit{W}(t)}}$$

(14)

Here, we have introduced the following notation.

Y(t) ∈ ℝn	state of the (extended) system
A(x) ∈ ℝn	deterministic drift vector field
B(x) ∈ ℝn×n	noise-coefficient matrix
D(x) ∈ ℝn×n	diffusion matrix
W (t) ∈ ℝn	n-dimensional Brownian motion
kT	temperature factor

The n × n diffusion matrix, D(x), is defined in terms of the noise coefficient matrix, B(x), as:

$$\mathit{D}(\mathit{x})\stackrel{\text{def}}{=}kT\mathit{B}(\mathit{x})\mathit{B}{(\mathit{x})}^T,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{all}\hspace{0.17em}\hspace{0.17em}\mathit{x}\in {ℝ}^n$$

(15)

where B(x)^T denotes the transpose of the real matrix, B(x). The diffusion matrix is symmetric and nonnegative definite. Depending on the particular bath that is used, the dimension, n, of Y (t) in Equation (14) is related to the dimension, 2d, of z(t) in Equation (11) by the inequality: n ≥ 2d. For example, in Nosé-Hoover Langevin dynamics, a single bath degree of freedom is added to Equation (11), so that n = 2d + 1, while in Langevin dynamics, the effect of the bath is modeled by added friction and Brownian forces that keep n = 2d. The Langevin Equation (1) can be put in the form of Equation (14) by letting x = (q, p),

$$\mathit{A}(\mathit{x})=\left(\begin{array}{c}{\mathit{m}}^{-1}\mathit{p}\\ -\nabla U(\mathit{q})-\gamma \mathit{p}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{B}=\sqrt{\gamma }\left(\begin{array}{cc}0& 0\\ 0& m^{1/2}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{W}=({\mathit{w}}_1,\cdots ,{\mathit{w}}_N)$$

(16)

where m = diag(m₁, · · ·, m_N).

Equation (14) generates a stochastic process, Y (t), that is a Markov diffusion process. We assume that this diffusion process admits a stationary distribution μ(dx), i.e., a probability distribution preserved by the dynamics [98,99]. We denote by ν(x) the density of this distribution. Even though the diffusion matrix in Equation (15) is not necessarily positive definite, one can use the Hörmander’s condition to prove that the process, Y (t), is an ergodic process with a unique stationary distribution [100,101]. By the ergodic theorem, it then follows that:

$$\frac{1}{T}{\int }_0^{T}f(\mathit{Y}(t))dt\to {\int }_{{ℝ}^n}f(\mathit{x})\nu (\mathit{x})d\mathit{x},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{as}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}T\to \infty ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{a}.\text{s}.$$

(17)

where f(x) is a suitable test function.

The evolution of the probability density of the law of Y (t) at time t, ρ(t, x), satisfies the Fokker-Planck equation:

$$-\frac{\partial \rho }{\partial t}+L\rho =0$$

(18)

where ρ(0, ·) is the density of the initial distribution, Y (0) ∼ ρ(0, ·), and L is defined as the following second-order partial differential operator:

$$(Lf)(\mathit{x})\stackrel{\text{def}}{=}\text{div}\left(\text{div}\left(\mathit{D}(\mathit{x})f(\mathit{x})\right)-\mathit{A}(\mathit{x})f(\mathit{x})\right)$$

(19)

Since μ(dx) = ν(x)dx is a stationary distribution of Y (t), the probability density, ν(x), is a steady-state solution of Equation (18), i.e., it satisfies:

$$(L\nu )(\mathit{x})=0$$

(20)

Define the probability current as the vector field:

$$\mathit{j}(\mathit{x})\stackrel{\text{def}}{=}\text{div}\left(\mathit{D}(\mathit{x})\nu (\mathit{x})\right)-\mathit{A}(\mathit{x})\nu (\mathit{x})$$

(21)

The stationarity condition in Equation (20) implies that j(x) is divergenceless. In the zero-current case, the diffusion process, Y (t), is reversible, and the stationary density ν(x) is called the equilibrium probability density of the diffusion [102].

In this case, the operator, L, is self-adjoint, in the sense that:

$${\langle Lf,g\rangle }_{\nu }={\langle f,Lg\rangle }_{\nu }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\text{all}\hspace{0.17em}\hspace{0.17em}\text{suitable}\hspace{0.17em}\text{test}\hspace{0.17em}\text{functions}\hspace{0.17em}f,g$$

(22)

where 〈·, ·〉_ν denotes an L² inner product weighted by the density, ν(x). This property implies that the diffusion is ν-symmetric [103]:

$$\nu (\mathit{x})p_t(\mathit{x},\mathit{y})=\nu (\mathit{y})p_t(\mathit{y},\mathit{x})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\hspace{0.17em}\text{all}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t>0$$

(23)

where p_t(x, y) denotes the transition probability density of Y (t). Indeed, Equation (22) is simply an infinitesimal version of Equation (23), which is referred to as the detailed balance condition. In the self-adjoint case, the drift is uniquely determined by the diffusion matrix and the stationary density ν(x):

$$\mathit{j}(\mathit{x})=\mathbf{0}\Rightarrow \mathit{A}(\mathit{x})=\frac{1}{\nu (\mathit{x})}\text{div}\left(\mathit{D}(\mathit{x})\nu (\mathit{x})\right)$$

(24)

Long-time stable explicit schemes adapted to this structure have been recently developed [92].

5.3. Splitting Approach to MD Simulation

We are now in a position to explain our general approach for deriving a long-time stable scheme for Equation (14). Crucial to our approach is that in MD simulation, we usually have a formula for a function proportional to the stationary density ν(x). Following [90], we can split Equation (14) into:

$$d\mathit{Y}=-\mathit{D}(\mathit{Y})\nabla H_{\nu }(\mathit{Y})dt+\text{div}\mathit{D}(\mathit{Y})dt+\sqrt{2kT}B(\mathit{Y})dW$$

(25)

$$\dot{\mathit{Y}}=\mathit{A}(\mathit{Y})+\mathit{D}(\mathit{Y})\nabla H_{\nu }(\mathit{Y})$$

(26)

where we have introduced H_ν(x) = −(log ν)(x). An exact splitting method preserves μ(dx). It is formed by taking the exact solution (in law) of Equation (25) in composition with the exact flow of Equation (26). The process produced by Equation (25) is self-adjoint with respect to ν(x). Moreover, the stationarity of ν(x) implies that the flow of the ODE (26) preserves it. Since each step is preservative, their composition is, too.

In place of the exact splitting, a Metropolized explicit integrator can be used for Equation (25) [92], and a measure-preserving scheme can be designed to solve the ODE [72,104]. In [92], explicit schemes are introduced for Equation (25) that: (i) sample the exact equilibrium probability density of the SDE when this density exists (i.e., whenever ν(x) is normalizable); (ii) generates a weakly accurate approximation to the solution of Equation (14) at constant kT ; (iii) acquire higher order accuracy in the small noise limit, kT → 0; and (iv) avoid computing the divergence of the diffusion matrix D(x). Compared to the methods in [72], the main novelty of these schemes stems from (iii) and (iv). The resulting explicit splitting method is accurate, since it is an additive splitting of Equation (14); and typically ergodic when the continuous process is ergodic [72].

This type of splitting of Equation (14) is quite natural and has been used before in MD [43,87], dissipative particle dynamics [105,106] and the simulation of inertial particles [107]. Other closely related schemes for Equation (14) include Brünger-Brooks-Karplus (BBK) [35], van Gunsteren and Berendsen (vGB) [108] and the Langevin-Impulse (LI) methods [41] and quasi-symplectic integrators [44]. However, for general MD force fields, none of these explicit integrators are long-time stable. Our framework to stabilize explicit MD integrators is the Metropolis-Hastings algorithm.

5.4. Metropolis-Hastings Algorithm

A Metropolis-Hastings method is a Monte Carlo method for producing samples from a probability distribution, given a formula for a function proportional to its density [74,75]. The algorithm consists of two sub-steps: firstly, a proposal move is generated according to a transition density, g(x, y); and secondly, this proposal move is accepted or rejected with a probability:

$$\alpha (\mathit{x},\mathit{y})=1\wedge \frac{g(\mathit{x},\mathit{y})\nu (\mathit{y})}{g(\mathit{x},\mathit{y})\nu (\mathit{x})}$$

(27)

Standard results on Metropolis-Hastings methods can be used to classify this algorithm as ergodic [100,109,110].

6. Conclusions

This paper provided an algorithmic introduction to time integrators for MD simulation. A quick overview of existing algorithms was given. When the derivatives of the potential are bounded, it is well known that these integrators work just fine: they are convergent on finite-time intervals and possess an invariant measure that is nearby the Boltzmann-Gibbs density. However, in realistic MD simulation, the derivatives of the potential are unbounded. This lack of regularity can cause numerical instabilities or artifacts in explicit integrators. The paper demonstrated how a Metropolis acceptance-rejection step can be added to explicit MD integrators to mitigate some of these problems and, in principle, obtain long-time stable and finite-time accurate schemes. A MATLAB implementation of Metropolis-corrected MD integrators was provided and used to compute the velocity autocorrelation of a sea of Lennard-Jones particles at various densities between the solid and liquid phases. The paper did not provide an in-depth review of the theory of Metropolis integrators, which can be found elsewhere [72,73].

Calculating the force field at every step dominates the overall computational cost of MD simulation. These force fields involve: bonded interactions and non-bonded Lennard-Jones and electrostatic interactions. The calculation of bonded interactions is straightforward to vectorize and scales like O(N). In addition, Lennard-Jones forces rapidly decay with interatomic distance. To a good approximation, every atom interacts only with neighbors within a sufficiently large ball. By using data structures, like neighbor lists and cell linked lists, these interactions can be calculated in O(N) steps, and therefore, the Lennard-Jones interactions can be calculated in O(N) steps [46]. On the other hand, the electrostatic energy between particles decays, like 1/r, where r denotes an interatomic distance, which leads to long-range interactions between atoms. Unlike Lennard-Jones interaction, this interaction cannot be cutoff without introducing large errors. In this case, one can use sophisticated techniques, like the fast multipole method, to rigorously handle such interactions in $𝒪$(N) steps [39,58].

However, the effect of these ‘mathematical tricks’ for fast calculation of the force field can become muted if the time step requirement for stability or accuracy becomes more severe in high dimension. This can happen in the Metropolis integrator, if the acceptance probability in Step 2 of Algorithm 3.1 or Step 3 of Algorithm 3.2 deteriorates in high dimension. The scaling of Metropolis algorithms has been quantified for the random walk Metropolis, hybrid Monte Carlo and Metropolis-adjusted Langevin algorithm (MALA) [111–115]. Since the acceptance probability is a function of an extensive quantity, the acceptance rate can artificially deteriorate with increasing system size, unless the time step is reduced. Because high acceptance rates are required to maintain dynamic accuracy, the dependence of the time step on system size limits the application of Metropolized schemes to large-scale systems. Fortunately, this scalability issue can often be resolved by using local, rather than global proposal moves, because the change in energy induced by a local move is typically an intensive quantity. For molecular dynamics calculations, this approach was pursued in [73]. Using dynamically consistent local moves (a so-called J-splitting [116]), it was shown that in certain situations, a scalable Metropolis integrator can be designed; however, the extent to which this strategy remedies the issue of high rejection rate in high dimension is not clear at this point and should be tested in applications.