Path-accelerated molecular dynamics: Parallel-in-time integration using path integrals

TL;DR

The paper introduces path-accelerated molecular dynamics (PAMD), a parallel-in-time method leveraging path integrals to significantly speed up long-timescale MD simulations on parallel architectures.

Contribution

It presents a novel parallelization approach for MD using path integrals, enabling substantial reductions in simulation time compared to traditional methods.

Findings

Achieved over 1000x speedup in Brownian dynamics simulations.

Demonstrated effectiveness on harmonic oscillator and Lennard-Jones systems.

Method is generalizable to other stochastic equations of motion.

Abstract

Massively parallel computer architectures create new opportunities for the performance of long-timescale molecular dynamics (MD) simulations. Here, we introduce the path-accelerated molecular dynamics (PAMD) method that takes advantage of distributed computing to reduce the wall-clock time of MD simulation via parallelization with respect to MD timesteps. The marginal distribution for the time evolution of a system is expressed in terms of a path integral, enabling the use of path sampling techniques to numerically integrate MD trajectories. By parallelizing the evaluation of the path action with respect to time and by initializing the path configurations from a non-equilibrium distribution, the algorithm enables significant speedups in terms of the length of MD trajectories that can be integrated in a given amount of wall-clock time. The method is demonstrated for Brownian dynamics, although it is generalizable to other stochastic equations of motion including open systems. We apply the method to two simple systems, a harmonic oscillator and a Lennard-Jones liquid, and we show that in comparison to the conventional Euler integration scheme for Brownian dynamics, the new method can reduce the wall-clock time for integrating trajectories of a given length by more than three orders of magnitude in the former system and more than two in the latter. This new method for parallelizing MD in the dimension of time can be trivially combined with algorithms for parallelizing the MD force evaluation to achieve further speedup.

Tables

Table 1. Table 1 .: Summary of notation employed to describe the path-based MD integration scheme introduced in the current work. The gray region corresponds to the definitions of parameters specific to the path sampling algorithm used.

$dt$	timestep for discretization of the sampled path
$N_{\mathrm{path}}$	number of timesteps in the sampled path
$N_{\mathrm{procs}}$	number of processors for parallel-in-time force evaluations
$N_{\mathrm{force}}$	number of force evaluations per shifting event
$N_{\mathrm{sample}}$	number of MC steps per shifting event
$N_{\mathrm{shift}}$	number of timesteps shifted
$L$	total number of levels in the sampled path
$l_{\min }$	finest sampled level of the path
$l_{\max }$	coarsest sampled level of the path

Table 2. Table 2 .: Summary of PAMD simulation parameters used for the harmonic oscillator application. The gray region corresponds to the parameters specific to the path sampling algorithm used.

	Euler
$d{t}_E$	$0.025$
$E_{\mathrm{eq}}$	$0.3\%$
$E_{\mathrm{dyn}}$	$1.0\%$
	PAMD
	Simulation $1$	Simulation $2$	Simulation $3$
$dt$	$0.4$	$0.4$	$0.4$
$N_{\mathrm{path}}$	$16$	$256$	$4096$
$N_{\mathrm{sample}}$	$2$	$1$	$1$
$N_{\mathrm{shift}}$	$2$	$8$	$64$
$L$	$4$	$8$	$12$
$l_{\min }$	$3$	$4$	$4$
$l_{\max }$	$3$	$5$	$8$
$\chi$	$16$	$128$	$1024$
$E_{\mathrm{eq}}$	$0.3\%$	$0.3\%$	$0.5\%$
$E_{\mathrm{dyn}}$	$0.3\%$	$1.1\%$	$2.9\%$

Table 3. Table 3 .: Summary of PAMD simulation parameters used for the application to the Lennard-Jones liquid. The gray region corresponds to the parameters specific to the path sampling algorithm used.

	Euler
$d{t}_E$	$2.5\times {10}^{-4}$
$E_{\mathrm{eq}}$	$1.4\%$
$E_{\mathrm{dyn}}$	$0.1\%$
	PAMD
	Simulation $1$	Simulation $2$
$dt$	$5\times {10}^{-4}$	$5\times {10}^{-4}$
$N_{\mathrm{path}}$	$128$	$4096$
$N_{\mathrm{sample}}$	$1$	$1$
$N_{\mathrm{shift}}$	$8$	$64$
$L$	$7$	$12$
$l_{\min }$	$3$	$4$
$l_{\max }$	$5$	$7$
$\chi$	$16$	$128$
$E_{\mathrm{eq}}$	$3.5\%$	$3.7\%$
$E_{\mathrm{dyn}}$	$1.0\%$	$1.9\%$

Equations

$$\dot{x}(t)=-{\gamma }^{-1}{V}^{\prime }(x(t))+\sqrt{2D}\dot{w}(t),$$

$$x(t+dt)-x(t)=-{\gamma }^{-1}{V}^{\prime }(x(t))dt+\sqrt{2Ddt}\xi ,$$

$$K($$

$$\exp \left\{-\frac{dt}{4D}{\left(\frac{x(t+dt)-x(t)}{dt}+\frac{V^{\prime }(x(t))}{\gamma }\right)}^2\right\},$$

$$\prod _{n=0}^{N-1}K(x(t_{n+1})|x(t_n);dt)\equiv e^{-S[\mathbf{X}]},$$

$$K(x(T+t)|x(t);T)\propto {\int }_{\mathbb{R}}\mathrm{d}{x}_1\cdots {\int }_{\mathbb{R}}\mathrm{d}{x}_{N-1}{e}^{-S[\mathbf{X}]},$$

$$D_{\mathrm{KL}}(t)={⟨{\int }_{\mathbb{R}}\mathrm{d}{x}_{t}P(x_t|x_0;t)\log \frac{P(x_t|x_0;t)}{Q(x_t|x_0;t)}⟩}_{x_0},$$

$$P(x_t|x_0;t)\propto \exp \left\{-\frac{\beta k}{2}\frac{(x_t-e^{-{\gamma }^{-1}kt}{x}_0{)}^2}{(1-e^{-2{\gamma }^{-1}kt})}\right\},$$

$$\chi =N_{\mathrm{shift}}\cdot \frac{N_{\mathrm{procs}}}{N_{\mathrm{force}}}\cdot \frac{dt}{d{t}_{\mathrm{E}}}.$$

$$N_{\mathrm{force}}=N_{\mathrm{sample}}\cdot N_{\mathrm{path}}.$$

$$N_{\mathrm{procs}}=N_{\mathrm{path}}.$$

$$\chi =\frac{N_{\mathrm{shift}}}{N_{\mathrm{sample}}}\cdot \frac{dt}{d{t}_{\mathrm{E}}}.$$

$$E_{\mathrm{eq}}=\frac{\sqrt{{\int }_{\mathbb{R}}\mathrm{d}x|P(x)-P_{\mathrm{s}}(x){|}^2}}{Z}$$

$$C(t)=⟨x(t^{\prime })x(t^{\prime }+t)⟩=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\mathrm{d}{t}^{\prime }x(t^{\prime })x(t^{\prime }+t),$$

$$C(t)=⟨x^2⟩\exp (-\kappa t),$$

$$E_{\mathrm{dyn}}=\frac{\sqrt{(\kappa -{\kappa }_{\mathrm{s}}{)}^2}}{k},$$

$$U(r)=\left\{\begin{array}[]{rr}u(r)-u(r_{\mathrm{c}})-(r-r_{\mathrm{c}})\,u^{\prime}(r_{\mathrm{c}}),&r\leq r_{\mathrm{c}}\\ 0,&r>r_{\mathrm{c}}\end{array}\right.,$$

$$V(r;{\sigma }_{\mathrm{HS}}(\tau ))=U(r)+U_{\mathrm{HS}}(r;{\sigma }_{\mathrm{HS}}(\tau )),$$

$$U_{\mathrm{HS}}(r;\sigma_{\mathrm{HS}})=\left\{\begin{array}[]{rr}+\infty,&r\leq\sigma_{\mathrm{HS}}\\ 0,&r>\sigma_{\mathrm{HS}}\end{array}\right..$$

$$E_{\mathrm{eq}}=\frac{\sqrt{{\int }_0^{r_{\mathrm{c}}}\mathrm{d}r|g(r)-g_{\mathrm{s}}(r){|}^2}}{{\int }_0^{r_{\mathrm{c}}}\mathrm{d}rg(r)}$$

$$E_{\mathrm{dyn}}=\frac{\sqrt{(D-D_{\mathrm{s}}{)}^2}}{D},$$

$$⟨R^2(t)⟩=\underset{T\to \infty }{\lim }{\int }_0^T\mathrm{d}{t}^{\prime }\frac{1}{N}\sum _{i=1}^N|{\mathbf{r}}_i(t^{\prime }+t)-{\mathbf{r}}_i(t^{\prime }){|}^2$$

Full text

I Introduction

Molecular dynamics (MD) Frenkel and Smit (2002); Allen and Tildesley (2017) is the central tool for simulating chemical, biological, and materials systems, with new algorithms and hardware expanding the range of accessible timescales and lengthscales Durrant and McCammon (2011); Dror et al. (2012); Voter, Montalenti, and Germann (2002). Faster processors have played an important role in this expansion, although the most dramatic improvements in recent years have come from the number of available processors, rather than the clock-speed of the individual cores Keyes (2007); Dongarra et al. (2011). In particular, highly multi-threaded computer architectures have been used to parallelize the MD force evaluation, greatly reducing the wall-clock time needed to perform an individual MD step Plimpton (1995); Phillips et al. (2005); Shaw et al. (2009); Salomon-Ferrer et al. (2013); Páll et al. (2015); Grossman et al. (2015). However, despite this progress in the parallelization of MD simulations with respect to the force evaluations (i.e., in space), less attention has been dedicated to the notion of parallelization with respect to the MD timesteps (i.e., in time).

The sequential nature of MD (i.e., the need to have access to a given timestep before the next timestep can be computed) would seem to discount the possibility of exploiting parallelization in time; nonetheless, methods for parallel-in-time integration are being developed and applied to MD simulation. Most approaches Lions, Maday, and Turinici (2001); Farhat and Chandesris (2003); Garrido et al. (2006); Emmett and Minion (2012) are based on a prediction-correction paradigm that combines fine (i.e., accurate and expensive) and coarse (i.e., inaccurate and inexpensive) solvers to iteratively refine approximations of a trajectory in a convergent and parallel-in-time fashion. A range of coarse solvers and iteration schemes have been employed to evaluate MD trajectories of molecular systems with parallelization in the time domain Baffico et al. (2002); Yanan, Srinivasan, and Chandra (2006); Speck et al. (2012); Bylaska, Weare, and Weare (2013); Blumers, Li, and Karniadakis (2019), leading to order-of-magnitude reductions in the wall-clock time-to-solution with respect to sequential integration at the fine level of accuracy. Schemes for approximate long-timescale integration via trajectory splicing are an alternative route to parallelization in time, yielding accurate time evolution for systems that exhibit strong timescale separation on well-characterized regions of the potential energy landscape Perez et al. (2016).

The current work takes a different approach to parallelizing MD in time. We demonstrate that by working with ensembles of trajectories in a path-integral framework, multiple processors can be employed to reduce the wall-clock time needed to evolve an MD trajectory of arbitrary length, without resorting to parallelization of the MD force evaluation. This method of parallelization for MD trajectories is independent of, and thus entirely complementary to, parallelization of the MD force evaluations, and it creates new opportunities to harness large numbers of available computer processors for the generation of long-timescale MD trajectories.

II Method

II.1 MD integration based on path distributions

In this work, we focus on the MD equation of motion governing Brownian (i.e., overdamped Langevin) dynamics under potential $V$ at temperature $\beta^{-1}$,

[TABLE]

where the diffusion coefficient $D$ and the friction coefficient $\gamma$ are related by the Einstein relation $D=(\beta\gamma)^{-1}$, and $w(t)$ is the standard Wiener process. MD trajectories can be generated by discretizing Eq. 1 with various numerical integration schemes Brünger, Brooks, and Karplus (1984); Brańka and Heyes (1998); Ricci and Ciccotti (2003); Bussi and Parrinello (2007); Bou-Rabee (2014), such as the forward Euler algorithm Allen and Tildesley (2017)

[TABLE]

where $dt$ is the discretization timestep, and $\xi$ a standard Gaussian random variate. The marginal distribution associated with time evolution of the system by $dt$ according to Eq. 2 is Risken and Frank (1996)

[TABLE]

such that the likelihood of a MD trajectory of length $T=N\,dt$ that evolves the system along positions $\bm{X}=\{x(t_{0}),x(t_{1}),\ldots,x(t_{N})\}$ at times $t_{n}=t+n\,dt$ is

[TABLE]

where $S[\bm{X}]$ is the action associated with the MD trajectory. From Eq. 4, the position of the time-evolved system at time $T$ has a marginal distribution given by the path integral

[TABLE]

where $x_{n}=x(t_{n})$. It is clear that this path-integral formulation of the ensemble of MD trajectories provides an equivalent description of the time evolution of the system as Eq. 2. Numerous studies have explored this path-integral formulation with variations of the underlying equation of motion and of the discretization of the action Pratt (1986); Olender and Elber (1996); Bolhuis et al. (2002); Miller III and Predescu (2007); Sivak, Chodera, and Crooks (2014).

Setting aside issues of efficiency until section II.2, we note that the path-integral formulation of the marginal distribution for the time-evolved system offers a simple MD integration scheme, illustrated in Fig. II.1. First, sampling from the distribution of paths of length $T$, with likelihood given by Eq. 4, is performed using Monte Carlo (MC) or related methods (Fig. II.1A) Ceperley (1995); Stuart, Voss, and Wilberg (2004); Stoltz (2007); by drawing a realization from this distribution, we obtain a segment of MD trajectory from time [math] to time $T$ (illustrated by the heavy orange path in Fig. II.1A). Then, by shifting from $x(0)$ to $x(T)$ along the sampled path, we resolve a trajectory from $x(0)$ to $x(T)$ (represented by the heavy green path in Fig. II.1B) that is statistically equivalent to a realization from the Euler algorithm defined in Eq. 2. After shifting the tail of the path to $x(T)$, we restart the path sampling to extend the trajectory from time $T$ to time $2T$. Iteration of this scheme will lead to the numerical integration of a MD trajectory of arbitrary length in time.