Nonequilibrium candidate Monte Carlo:A new tool for efficient equilibrium simulation

Nonequilibrium candidate Monte Carlo:
A new tool for efficient equilibrium simulation

Jerome P. Nilmeier , Gavin E. Crooks , David D. L. Minh ,    John D. Chodera  To whom correspondence should be addressed. E-mail: jchodera@berkeley.edu
Abstract

Metropolis Monte Carlo simulation is a powerful tool for studying the equilibrium properties of matter. In complex condensed-phase systems, however, it is difficult to design Monte Carlo moves with high acceptance probabilities that also rapidly sample uncorrelated configurations. Here, we introduce a new class of moves based on nonequilibrium dynamics: candidate configurations are generated through a finite-time process in which a system is actively driven out of equilibrium, and accepted with criteria that preserve the equilibrium distribution. The acceptance rule is similar to the Metropolis acceptance probability, but related to the nonequilibrium work rather than the instantaneous energy difference. Our method is applicable to sampling from both a single thermodynamic state or a mixture of thermodynamic states, and allows both coordinates and thermodynamic parameters to be driven in nonequilibrium proposals. While generating finite-time switching trajectories incurs an additional cost, driving some degrees of freedom while allowing others to evolve naturally can lead to large enhancements in acceptance probabilities, greatly reducing structural correlation times. Using nonequilibrium driven processes vastly expands the repertoire of useful Monte Carlo proposals in simulations of dense solvated systems.

Metropolis-Hastings      Markov chain Monte Carlo      molecular dynamics      expanded ensembles

Abbreviations: MC, Metropolis Monte Carlo; MD, molecular dynamics; MCMC, Markov chain Monte Carlo; NCMC, Nonequilibrium candidate Monte Carlo

I n this paper, we describe a new technique for constructing efficient Markov chain Monte Carlo (MCMC) [1] moves that both have high acceptance rates and also allow rapid transit through configuration space, greatly enhancing convergence rates in simulations of dense solvated systems. The Metropolis Monte Carlo [2, 3] sampling procedure is generalized by using nonequilibrium processes to generate candidates for equilibrium simulations. Within this framework, moves that are efficient for an isolated part of a system but lead to near-universal rejection in standard Monte Carlo simulations of dense mixtures can be converted to nonequilibrium processes that generate candidates with higher acceptance probabilities. In this new procedure, the acceptance criteria is related to the nonequilibrium work, rather than the potential energy difference used in traditional Monte Carlo moves.

Since their introduction in the mid-twentieth century, Metropolis Monte Carlo (MC) [2, 3] and molecular dynamics (MD) [4] simulations have become favored tools for sampling from complex multidimensional distributions, such as configurations of microscopic physical systems in thermodynamic ensembles. However, these methods can produce highly correlated samples, leading to slow convergence of estimated expectations. While MD requires the use of small timesteps for numerical stability and to approximate sampling from the desired distribution, MC simulations can, in principle, make use of non-local moves that accelerate mixing of the Markov chain. Indeed, vast improvements in efficiency have been obtained by applying cleverly constructed move sets that exploit physical intuition about the system under study, such as cluster moves in Potts and Ising model simulations [5, 6].

Designing efficient moves requires striking a balance between rapid traversal of phase space and ensuring reasonable acceptance probabilities. For complex heterogeneous systems such as solvated biomolecules, achieving this balance remains challenging. Typically, efficient moves exploit physical insight into kinetically slow processes and energetically favorable configurations. Often, the experimenter may possess physical insight about one component in the system (e.g. a biomolecule) that permits the design of moves that would be efficient in the absence of other components (e.g. solvent), but encounter energetically unfavorable interactions in their presence, reducing acceptance rates to levels where standard MC provides no benefit. \colorblack As an illustrative example, consider a bistable dimer—a pair of particles interacting with a potential with minima in compact or extended configurations, separated by a high barrier (see Fig. 1). For simulations of this system in a vacuum, a simple and effective standard MC move is to instantaneously increase the interparticle distance from a compact to extended configuration (or conversely, to decrease the distance from an extended to compact configuration). When the dimer is immersed in a dense solvent, however, this move is met with near-universal rejection because solvent molecules overlap with proposed configurations.

\color

black One approach that can allow unperturbed degrees of freedom to relax, and hence maintain a reasonable acceptance rate, is to use a nonequilibrium process to generate candidate configurations. Using the appropriate acceptance criterion for the final configuration will preserve the equilibrium distribution. \colorblack In the case of the bistable dimer immersed in dense solvent, the extension (contraction) may be carried out over a finite number of increments interleaved with standard Metropolis Monte Carlo or molecular dynamics steps that allow the solvent to reorganize to avoid overlap with the dimer particles.

The basic idea of using nonequilibrium driven processes as Monte Carlo moves has precedents in both the statistical [7, 8] and chemical [9, 10, 11, 12] literature. Among the latter, Athènes developed “work-bias Monte Carlo” to enhance acceptance rates in grand canonical Monte Carlo simulations [9], Stern presented a scheme to sample an equilibrium mixture of protonation states at constant pH in explicit solvent [11] (though an inexact variant was proposed earlier [10]), and \colorblack Nilmeier et al. [12] proposed the driving of a subset of degrees of freedom to enhance acceptance rates (using an approximate acceptance criteria). Nonequilibrium processes have also been used to generate configurations for parallel tempering simulations [13, 14, 15].

\color

black Here, we unify these ideas and significantly extend the application of nonequilibrium moves in physical simulations. We present a theoretical framework, nonequilibrium candidate Monte Carlo (NCMC), that is applicable to both single thermodynamic states (e.g. NVT, NpT, VT ensembles) as well as mixtures of thermodynamic states (e.g. expanded ensemble [16, 17] simulations). Nonequilibrium proposals may drive a subset of degrees of freedom, the thermodynamic parameters characterizing the equilibrium distribution, or both, significantly expanding the repertoire of Monte Carlo moves that lead to high acceptance and efficient mixing in dense condensed-phase systems.

\@xsect

For physical systems in equilibrium, the probability of observing a microstate is given by the Boltzmann distribution,

(0)

where denotes a microstate of the system (which may include coordinates, momenta, and other dynamical variables, such as simulation box dimensions), denotes a set of thermodynamic parameters whose values define a thermodynamic state, and is a normalizing constant known as the partition function.

The reduced potential depends on the thermodynamic ensemble under consideration [18]. For instance, in an isothermal-isobaric () ensemble, the reduced potential will assume the form,

(0)

which depends on the Hamiltonian (which may include an external biasing potential, and is presumed to be invariant under momentum inversion) and the system volume . In this ensemble, the vector of controllable thermodynamic parameters includes the inverse temperature , the Hamiltonian , and external pressure . Other thermodynamic parameters and their conjugate variables can be included or excluded to generate alternative physical (or unphysical) ensembles.

To allow sampling from multiple thermodynamic states within a single simulation, we also define an expanded ensemble [16, 17], which specifies a joint distribution for in a weighted mixture of thermodynamic states,

(0)

where specifies an externally-imposed weight for state . Here, may assume values in a discrete or continuous space . If the set consists of a single value of , a single thermodynamic state is sampled, and . These thermodynamic states may correspond to a variety of different states of interest, such as temperatures in a simulated tempering simulation [19], alchemical states in a simulated scaling simulation [20], or protonation states in a constant-pH simulation [21].

Figure 1: Bistable dimer potential and instantaneous MC moves in WCA fluid. An extension move increases the dimer extension by , while a compaction move decreases the dimer extension by . Both move types meet with near-universal rejection when implemented as instantaneous MC moves in a dense WCA fluid. Note that the lower panel is only a cartoon — the actual described simulation is of a 3D system.
\@xsect

We first describe the general form of NCMC. At the start of an iteration, the current sample in the Markov chain, , which is assumed to be drawn from , is used to initialize a trajectory, . A candidate configuration is then proposed through a nonequilibrium process in which a set of degrees of freedom and/or thermodynamic parameters may be driven according to some protocol [22] selected with a probability dependent only on . Even if we only wish to sample from a single thermodynamic state , we may use a protocol that transiently drives the thermodynamic parameters away from and back again (as in Ref. [14]). Finally, an acceptance probability is computed and used to decide whether the next sample in the Markov chain, , is the candidate, , or the momentum reversal of the initial sample, .

An NCMC move begins by selecting a protocol from a set of possible protocols with probability , such that there exists a reverse protocol labeled as (to be defined momentarily) with . A protocol specifies both a series of perturbation kernels and propagation kernels , arranged in an alternating pattern such that . Both and are conditional probabilities of given any , and must satisfy the requirement that if , then , for substituted by and .

Each perturbation kernel drives some or all of the degrees of freedom in a stochastic or deterministic way (e.g. by driving a torsion angle, a distance between two atoms, or the volume of the simulation cell). Similarly, each propagation kernel propagates some or all of the coordinates of the system at fixed according to some form of MCMC or MD (e.g. Metropolis Monte Carlo [2, 3], velocity Verlet [23] deterministic dynamics, or overdamped Langevin stochastic dynamics [24, 25]) that may also depend on the time index . Interleaving perturbation and propagation allows for energetically unfavorable interactions introduced by perturbation to be relaxed during propagation, potentially increasing acceptance rates relative to the instantaneous perturbations of standard Metropolis Monte Carlo.

The procedure by which a trajectory is generated from an initial microstate according to a protocol can be illustrated by the scheme,

(0)

Application of the perturbation to generates a perturbed configuration , which is then propagated by the kernel to obtain .

The reverse protocol reverses the order in which the perturbation and propagation steps are applied, generating the time-reversed trajectory , where denotes with inverted momenta,

(0)

The next step in NCMC is to accept or reject as the next sample in the chain. To ensure that the stationary distribution is preserved, we enforce a strict pathwise form of detailed balance,555\colorblackThe described pathwise detailed balance condition is closely related to “super-detailed balance” (see, e.g. [26]), but also accounts for momentum reversal to extend its definition to include molecular dynamics integrators.

(0)

The quantity is the NCMC acceptance probability, while and denote the probability of generating trajectory from initial configuration using protocol , or from final configuration with protocol , respectively,

(0)
(0)
\color

blackSummation of Eq. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation over all trajectories starting with and ending with recovers the standard detailed balance condition (see Appendix for proof).

We define the ratio of proposal kernels as,

(0)

and the ratio of propagation kernels as the exponentiated difference in forward and backward conditional path actions as,

(0)

Using the above expressions and the momentum invariance property , we may write the ratio of acceptance probabilities as,

(0)

where is the energy difference. Eq. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation is the main result of this paper, and is highly general with regard to both the choice of protocol for perturbation and propagation. In subsequent sections, we discuss specific choices for these protocols that lead to particularly simple acceptance criteria.

Many choices of acceptance probabilities that satisfy Eq. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation are possible, including the well-known Metropolis-Hastings criterion [2, 3],

After generating and evaluating , we generate a uniform random variate . If , then the candidate becomes the next value in the chain, . Otherwise, it is rejected, we perform a momentum flip, and the next value becomes . Alternately, we may perform a momentum flip upon acceptance, and preserve the momentum upon rejection, . We cannot, however, ignore the momentum flip completely; as explained in the Appendix, it is necessary to preserve the equilibrium distribution.

We note that NCMC need not be used exclusively to sample from , but can be mixed with other MCMC moves or with MD [1]. For example, one may reinitialize velocities from the Maxwell-Boltzmann distribution after each NCMC step; this is a Gibbs sampling MCMC move using the marginal distribution for velocities.

\@xsect

A large variety of choices are available for the perturbation kernels . Through judicious selection of these kernels, a practitioner can design nonequilibrium proposals that carry some component of the system from one high-probability region to another with high acceptance rates. We briefly describe a few possibilities.

\@xsect

Suppose we wish to drive a torsion angle (an angle subtended by four bonded atoms) stochastically by rotating it to a new torsion angle (holding bond lengths and angles fixed)according to some probability, such as the von Mises circular distribution centered on ,

(0)

with denoting the modified Bessel function of order zero and taking the role of a dimensionless force constant. Because the stochastic perturbation is made in a non-Cartesian coordinate, a Jacobian must be included to compute in Cartesian coordinates, resulting in the ratio,

(0)

where because the transformation (a rotation about a bond vector) preserves the Cartesian phase space volume.

\@xsect

Instead of perturbing the torsion angle stochastically, we can deterministically drive it in small, fixed increments . In this case, we effectively define an invertible map that takes , such that differs from only in the rotation of the specified torsion by . To implement this, we may choose a perturbation from a distribution where have equal probability, and drive from its current value to a final value over steps in equal increments, such that is constrained to . In this case, , where the Jacobian represents the factor by which Cartesian phase space is compressed on the application of the map , which is again unity for rotation about a torsion angle by , and, due to the invertibility of the map, the ratio .

\@xsect

Another possible deterministic perturbation kernel is simulation box scaling. A barostat can be implemented by proposing propagation kernels that scale the molecular centers and box geometry by a factor with chosen uniformly from applied as a factor of over the course of steps. In this case, the perturbation kernel is a delta function that compresses or expands the molecular centers and box geometry. Since the Jacobian is the ratio of infinitesimal volumes upon scaling, the ratio of perturbation kernels is , where denotes the number of molecular centers.

\@xsect

In many driven nonequilibrium processes, there is no direct perturbation to the coordinates, such that and the ratio . Instead, only the thermodynamic parameters are varied in time, carrying the system out of equilibrium through action of the propagation kernels. We recover Neal’s method [7] if the reduced potential is a simple linear interpolation such that , the probability of choosing protocol is symmetric with , and MC [2, 3] is used for the propagation kernel .

\@xsect

The choice of propagation kernels available is also very broad. If strong driving is performed in selection of , one may elect to choose a time-independent propagation kernel that samples from a stationary distribution of interest. Alternatively, a strongly time-dependent could be selected to transiently drive the system out of equilibrium, or from the equilibrium distribution at one thermodynamic state to another. Some possible choices are described below.

\@xsect

One may propagate some or all of a system’s degrees of freedom (e.g. those not affected by the perturbation kernel ) by a method that satisfies detailed balance in ,

(0)

where for a specified . Many MCMC methods [1], including Metropolis [2, 3] and various hybrid Monte Carlo (HMC) algorithms that combine discrete-timestep integrators with Monte Carlo acceptance/rejection steps [27, 28], obey detailed balance and are commonly used for the simulation of physical systems.

By analogy with Crooks [29], we define a work and heat for the nonequilibrium driven process,

(0)
(0)

such that , a restatement of the first law of thermodynamics.

The conditional path action difference can then be written in terms of the heat of the process, ,

(0)

leading to an acceptance probability similar to standard MC, except that the work, , replaces the instantaneous potential energy difference,

(0)
\@xsect

When an isolated system is propagated by a symplectic integrator—a reversible, deterministic integrator that preserves phase space volume—the propagation kernels follow . Hence, and the acceptance ratio is,

(0)

where is the energy difference. The equivalence of the work and energy difference for volume-preserving integrators was previously recognized in the context of fluctuation theorem calculations [30, 31].

Symplectic integrators include velocity Verlet [23]. These integrators are also symplectic when utilizing constraints through the application of algorithms such as RATTLE [32], provided that the constraints are iterated to convergence each timestep [33].

\@xsect\color

black Stochastic integrators such as velocity Verlet discretizations of Langevin dynamics [34, 35] sample a modified distribution that differs from the desired equilibrium distribution in a timestep-dependent manner [36]. While this modified distribution may be difficult or impossible to compute in order to recover equilibrium properties by reweighting, computation of the relative action is relatively straightforward, and the NCMC acceptance criteria ensures that the NCMC-sampled configurations are distributed according to the desired equilibrium ensemble666\colorblackAlternatives to using NCMC to correct stochastic integration include introducing a Metropolization correction after each step, as in the generalized hybrid Monte Carlo (GHMC) integrator we use in the example [37, 1, 36, 28].. \colorblack As examples, we compute for the overdamped Langevin (Brownian) dynamics integrator of Ermak and Yeh [24, 25] and the Brünger-Brooks-Karplus (BBK) Langevin integrator [38, 39, 40] in the Appendix.

\@xsect
Figure 2: Trajectories of WCA dimer system in vacuum and solvent. Left: The dimer extension as a function of simulation iteration. The dotted horizontal line denotes division between compact and extended configurations. The quantity printed above each plot indicates the estimated integrated autocorrelation time for the dimer extension . Right: Histogram accumulated over trajectory (black), with true equilibrium distribution (red). Plot titles denote whether simulation was run in vacuum (vacuum) or dense WCA fluid (solvent), and whether the simulation utilized only 500 GHMC steps per iteration (MD) or included instantaneous MC (MC) or 2048-step NCMC moves (NCMC) following each iteration.

To demonstrate NCMC, we ran simulations of a bistable dimer (adapted from Section 1.3.2.4 of Ref. [36]) in vacuum as well as a dense fluid. The dimer consists of a pair of “bonded” particles interacting via a double-well potential, with minima at (compact) and (extended), and a barrier \colorblack(see Fig. 1). In the solvated simulations, the dimer was immersed in a dense bath (reduced density ) of particles that interact with the bonded particles and each other via the Weeks-Chandler-Andersen (WCA) soft repulsive potential [41]. Each simulation “iteration” consisted of velocity reassignment from the Maxwell-Boltzmann distribution, 500 steps of generalized hybrid Monte Carlo (GHMC) dynamics [37, 1, 36, 28] (essentially, a Metropolis-corrected form of Langevin dynamics, henceforth referred to here as MD), optionally followed by either an instantaneous MC move or an NCMC move.

\color

black The rate at which effectively uncorrelated samples are generated can be quantified in terms of the correlation time for the dimer extension (shown in Fig. 2). This time represents the asymptotic decay time constant for the correlation function , which will behave like

(0)

for large , where and . The correlation time is related to the statistical inefficiency, , a factor that describes the number of iterations necessary to generate an effectively uncorrelated sample [42].

For the MD simulation in vacuum (Fig. 2, top trace), we observe slow hopping between compact and extended configurations with a correlation time iterations, resulting in slow convergence of the histogram. Introducing instantaneous MC dimer extension/contraction moves that modify the dimer extension by reduces the correlation time to , such that an uncorrelated configuration is generated after each iteration of 500 MD steps and one instantaneous MC step (Fig. 2, second trace from top).

Figure 3: Acceptance probabilities of NCMC proposals. \colorblack Top: Acceptance probability of NCMC proposals as a function of length of nonequilibrium proposal trajectory (black dots), compared with instantaneous MC proposal (red line). Inset: Enlarged region with acceptance probability shown on linear scale. Estimated 95% confidence intervals are shown as vertical lines. \colorblack

When the dimer is immersed in a dense fluid of WCA particles, however, iterations consisting of 500 MD steps alone result in extremely slow barrier crossings, requiring g 600 iterations to produce an uncorrelated sample (Fig. 2, middle trace). Unlike in vacuum, the introduction of instantaneous MC moves does not significantly reduce the correlation time (Fig. 2, second trace from bottom). However, performing these same dimer expansion and contraction moves over 2048-step NCMC moves (Fig. 2, bottom trace) allows the system to rapidly mix between both compact and extended states with a correlation time of iterations. While each iteration requires a 5-fold increase in computational effort (500 MD steps + 2048 NCMC switching steps = 2548 force evaluations, versus 500 force evaluations for MD alone), a 67-fold reduction in correlation time is achieved, resulting in a remarkable order-of-magnitude gain in overall efficiency.

The length of the NCMC switching process is a free parameter that may be tuned to further improve efficiency. Towards this end, we estimated the acceptance probability of the extension/contraction moves in dense solvent as a function of switching length (Fig. 3). While instantaneous MC proposals of are only accepted with probability (the error is this quantity is likely underestimated due to its extremity), dividing the move into smaller steps boosts the acceptance rate to a level useful in condensed-phase simulation. If the move is divided into a small number of steps (1 to 8), there is little to no increase in acceptance rate, but for an intermediate number of steps (16 to 1024), there is a superlinear boost in the acceptance probability relative to the length of the switching process. The acceptance probability finally reaches useful levels around 2000 steps, achieving an acceptance rate of 12% using nonequilibrium proposal trajectories of 2048 steps, or 38% for 8192 steps.

In general, there is no direct relationship between acceptance probability and efficiency. Under certain assumptions relevant to the bistable dimer, however, it is possible to link the NCMC acceptance probability to , an indirect estimate of the correlation time,

(0)

where and are correlation times for iterations consisting solely of MD or NCMC moves, respectively. The latter correlation time may be estimated from the average acceptance probability by (see Appendix for derivation).

As shown in Fig. 4, the effective correlation time is only diminished when the NCMC acceptance probability is large enough such that , which occurs when (about 256 switching steps or more). For shorter switching times, even though the acceptance probability is high relative to instantaneous MC, it is still too small to significantly reduce .

Figure 4: Statistical efficiency gain of NCMC proposals relative to instantaneous MC proposals. \colorblack Top: Effective correlation time , in iterations, for MD+NCMC (black dots) compared to MD alone (red line). Bottom: Relative statistical efficiency of MD+NCMC, in terms of number of uncorrelated configurations generated for a fixed amount of computational effort, for MD+NCMC (black dots) relative to MD alone (red line). \colorblack

When comparing efficiency, we are most interested in the rate of generating uncorrelated configurations for a given amount of computational effort. Relative to MD alone, this rate is described by the efficiency gain,

(0)

Here, steps per iteration, and is again varied from 1 to 8192 steps. The results are shown in the bottom panel of Fig. 4. Surprisingly, there is actually a slight loss in efficiency at short switching times—dropping to a minimum of the efficiency of MD alone at 128 steps—followed by a rapid gain in efficiency, plateauing at an efficiency gain of the efficiency of MD alone for 2048–4096-step NCMC proposals. (A similar plateau behavior is observed in the modified parallel tempering protocol of Ref. [15].) After this point, longer switching times do not achieve as high of an efficiency gain; even though the acceptance rate continues to increase as the number of NCMC switching steps is doubled again to 8192 steps, the reduction in correlation time is not sufficient to offset the additional cost of these moves. \colorblack

\@xsect\color

black We have described a procedure—nonequilibrium candidate Monte Carlo (NCMC)—by which nonequilibrium proposals can be used within MCMC simulations to enhance acceptance rates. In our illustration, we have demonstrated how its use can lead to large improvements in statistical efficiency—the rate at which uncorrelated samples are generated for a fixed amount of computational effort. In other applications, whether similarly large efficiency gains are achieved will depend on the precise nature of the system under study and the nonequilibrium proposals introduced. The most straightforward approach—borrowing Metropolis Monte Carlo proposals that are reasonable for one component of the system in isolation, and converting them to nonequilibrium proposals—is likely to be a fruitful avenue for generating efficient schemes, as was demonstrated here for a simple system.

More generally, the problem of selecting efficient nonequilibrium proposals is similar to the problem of choosing good reaction coordinates, in that it is desirable to drive the system along (possibly complex) slow collective coordinates where orthogonal degrees of freedom relax quickly. The search for such collective coordinates is a topic of active research [43, 44, 45, 46, 47, 48, 49]. Given an initial nonequilibrium protocol, the issue of optimizing such a protocol to minimize dissipation (and maximize acceptance) is also a topic of active study, led by forays into the world of single-molecule measurement [50, 51, 52]. Recent work has also suggested that estimating the thermodynamic metric tensor along the nonequilibrium parameter switching path [53, 54, 55, 56], could prove useful in adaptively optimizing the switching protocol [57].

\color

black We note that switching trajectories contain potentially useful information. Indeed, several methods [58, 59, 56] have recently been developed to estimate equilibrium properties from nonequilibrium samples through the application of statistical estimator theory to nonequilibrium fluctuation theorems [30, 60, 61]; these are particularly relevant to switching between multiple thermodynamic states. Though the development of efficient estimators that utilize both nonequilibrium switching trials and sampled equilibrium data generated in NCMC simulations remains an open challenge, it is at least straightforward to incorporate information from rejected NCMC proposals in the estimation of equilibrium averages [26, 62].

Materials and Methods

\@xsect

The dimer system considered here consists of two particles that interact via a double-well “bonded” potential in the interatomic distance ,

(0)

with , , and , where . Simulations denoted as “vacuum” contain only these two particles, while simulations denoted as “solvated” also interact with a dense bath of particles via the WCA nonbonded potential [41],

(0)

with mass amu, Å, and . The nonbonded WCA interaction is excluded between the two bonded particles. The solvated system contains a total of 216 WCA particles at a reduced density of . For all simulations, the reduced temperature is . A custom Python code making use of the GPU-accelerated OpenMM package [63, 64, 65] and the PyOpenMM Python wrapper [66] was used to conduct the simulations. All scripts are available for download from http://simtk.org/home/ncmc.

To ensure that observed differences were not due to changes in the stationary distribution of the integrator, we used generalized hybrid Monte Carlo (GHMC) ([37, 1, 36, 28]) for all our simulations. GHMC is based on a velocity Verlet discretization [23] of Langevin dynamics—the two are equivalent in the limit of small timesteps - but includes an acceptance/rejection step to correct for errors introduced by finite timesteps so that the stationary distribution is the exact equilibrium distribution. We used a timestep of , where , and the collision rate was set to . With this timestep, the acceptance probability is %; the resulting dynamics closely approximates Langevin dynamics.

In simulations employing instantaneous Monte Carlo moves, a perturbation to the interatomic distance of the dimer was chosen according to,

(0)

The dimer was contracted or expanded about the bond midpoint to generate a new configuration with dimer extension from the old configuration with dimer extension , and the move accepted or rejected with the Metropolis-Hastings criterion [3],

where the Jacobian ratio term accounts for the expansion and contraction of phase space due to the Monte Carlo proposals.

\color

black For simulations employing -step NCMC moves, proposals were made by selecting a new velocity vector from the Maxwell-Boltzmann distribution, integrating steps of velocity Verlet dynamics [23] for all bath atoms as the dimer extension was driven from to in equal steps of size , and accepting or rejecting based on the modified Metropolis criteria for symplectic integrators (Eqs. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation and Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation),

(0)

The Jacobian ratio is also . \colorblack This scheme corresponds to a choice for the perturbation kernel of,

(0)

where denotes the dimer separation of configuration . The propagation kernel corresponds to velocity Verlet dynamics where the dimer atoms are held fixed in space. \colorblack The final configuration after the MC or NCMC rejection procedure was stored and plotted to generate Fig. 2.

The mean acceptance probabilities for each switching time can be estimated via the sample mean

(0)

For numerical stability, logarithms of were stored, as . We then estimated (shown in Fig. 4) as

(0)

where .

Integrated autocorrelation times were estimated using the rapid scheme described in Section 5.2 of Ref. [42].

The acceptance probabilities plotted in Fig. 4 were estimated from a trajectory consisting of 10 000 iterations of 2048-step NCMC, with 500 steps of GHMC dynamics in between each NCMC trial, to ensure equilibrium sampling. Prior to each 2048-step NCMC iteration, a velocity from the Maxwell-Boltzmann distribution was selected, and NCMC trial moves with varying switching times were conducted solely to accumulate statistics. The statistical error in the estimate of and the computed relative efficiency over instantaneous Monte Carlo was estimated by 1000 bootstrap trials, in which the dataset of 10 000 work samples was resampled with replacement in each bootstrap trial and 95% confidence intervals computed from the distribution over bootstrap replicates.

Figure 5: Umbrella sampling simulation of the dimer in WCA solvent. Left: The dimer extension as a function of simulation iteration. Right: The histogram accumulated over the trajectory, with the observed histogram in black and the reweighted histogram (corrected for the applied umbrella bias potential) in red.

The reference distribution for the interparticle distribution plotted in red on the right side of Fig. 2 was computed analytically for the vacuum system,

(0)

For the solvated system, this distribution was estimated from an umbrella sampling simulation employing a modified bonded potential intended to remove the barrier in between compact and extended states,

(0)

where , , and , with Å, and is the Heaviside function that assumes a value of unity for and zero otherwise. The true solvated interparticle distribution was estimated by reweighting the data produced from this simulation, using the relationship

(0)

where denotes the bond separation for sample , and a finite-width histogram bin was used instead of the delta function .


ACKNOWLEDGMENTS. The authors thank Gabriel Stoltz (CERMICS, École des Ponts ParisTech); Jed W. Pitera (IBM Almaden Research Center); \colorblackManuel Athenès (CEA Saclay); \colorblackFiras Hamze (D-Wave Systems); Yael Elmatad, Anna Schneider, Paul Nerenberg, Todd Gingrich, David Chandler, David Sivak, Phillip Geissler, Michael Grünwald, and Ulf Rörbæck Pederson (University of California, Berkeley); Vijay S. Pande (Stanford University); and Suriyanarayanan Vaikuntanathan, Andrew J. Ballard, and Christopher Jarzynski (University of Maryland), and Huafeng Xu (D. E. Shaw Research) for enlightening discussions on this topic and constructive feedback on this manuscript, \colorblackas well as the two anonymous referees for their helpful suggestions for improving clarity. JPN was supported by the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. GEC was funded by the Helios Solar Energy Research Center, which is supported by the Director, Office of Science, Office of Basic Energy Sciences of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. DDLM was funded by a Director’s Postdoctoral Fellowship from Argonne National Laboratory. JDC was supported through a distinguished postdoctoral fellowship from the California Institute for Quantitative Biosciences (QB3) at the University of California, Berkeley. Additionally, the authors are grateful to OpenMM developers Peter Eastman, Mark Friedrichs, Randy Radmer, and Christopher Bruns for their generous help with the OpenMM GPU-accelerated computing platform and associated PyOpenMM Python wrappers. This research was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and by UChicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the government.

References

  • 1. Liu, J. S. (2002) Monte Carlo strategies in scientific computing. (Springer-Verlag, New York), 2nd ed. edition.
  • 2. Metropolis, N, Rosenbluth, A. W, Rosenbluth, M. N, Teller, A. H, & Teller, E. (1953) Equation of state calculations by fast computing machines. J. Chem. Phys. 21:1087–1092.
  • 3. Hastings, W. K. (1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57:97–109.
  • 4. Rahman, A. (1964) Correlations in the motion of atoms in liquid argon. Phys. Rev. 136:A405–A411.
  • 5. Swendsen, R. H & Wang, J.-S. (1987) Nonuniversal critical dynamics in Monte Carlo simulations. Phys. Rev. Lett. 58:86–88.
  • 6. Wolff, U. (1989) Collective Monte Carlo updating for spin systems. Phys. Rev. Lett. 62:361–364.
  • 7. Neal, R. M. (2004) Taking bigger Metropolis steps by dragging fast variables, (Department of Statistics, University of Toronto), Technical Report 0411.
  • 8. Andrieu, C, Doucet, A, & Holenstein, R. (2010) Particle markov chain monte carlo methods. J. Royal Stat. Soc. B 72:269–342.
  • 9. Athènes, M. (2002) Computation of a chemical potential using a residence weight algorithm. Phys. Rev. E 66:046705.
  • 10. Bürgi, R, Kollman, P. A, & van Gunsteren, W. F. (2002) Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation. Proteins 47:469–480.
  • 11. Stern, H. A. (2007) Molecular simulation with variable protonation states at constant pH. J. Chem. Phys. 126:164112.
  • 12. Nilmeier, J & Jacobson, M. P. (2009) Monte carlo sampling with hierarchical move sets: POSH Monte Carlo. J. Chem. Theor. Comput. pp. 1968–1984.
  • 13. Opps, S. B & Schofield, J. (2001) Extended state-space monte carlo methods. Phys. Rev. E 63:56701.
  • 14. Brown, S & Head-Gordon, T. (2003) Cool walking: A new Markov chain Monte Carlo sampling method. J. Comput. Chem. 24:68–76.
  • 15. Ballard, A. J & Jarzynski, C. (2009) Replica exchange with nonequilibrium switches. Proc. Natl. Acad. Sci. USA 106:12224–12229.
  • 16. Lyubartsev, A. P, Martsinovski, A. A, Shevkunov, S. V, & Vorontsov-Velyaminov, P. N. (1992) New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles. J. Chem. Phys. 96:1776–1783.
  • 17. Park, S. (2008) Comparison of the serial and parallel algorithms of generalized ensemble simulations: An analytical approach. Phys. Rev. E 77:016709.
  • 18. Shirts, M. R & Chodera, J. D. (2008) Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 129:124105.
  • 19. Marinari, E & Parisi, G. (1992) Simulated tempering: a new Monte Carlo scheme. Europhys. Lett. 19:451–458.
  • 20. Li, H, Fajer, M, & Yang, W. (2007) Simulated scaling method for localized enhanced sampling and simultaneous “alchemical” free energy simulations: A general method for molecular mechanical, quantum mechanical, and quantum mechanical/molecular mechanical simulations. J. Chem. Phys. 126:024106.
  • 21. Mongan, J, Case, D. A, & McCammon, J. A. (2004) Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem. 25:2038–2048.
  • 22. Jarzynski, C. (2000) Hamiltonian derivation of a detailed fluctuation theorem. J. Stat. Phys. 98:77–102.
  • 23. Swope, W. C, Andersen, H. C, Berens, P. H, & Wilson, K. R. (1982) A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J. Chem. Phys. 76:637–649.
  • 24. Ermak, D. L & Yeh, Y. (1974) Equilibrium electrostatic effects on the behavior of polyions in solution: Polyion-mobile ion interaction. Chem. Phys. Lett. 24:243–248.
  • 25. Ermak, D. L. (1975) A computer simulation of charged particles in solution. I. Technique and equilibrium properties. J. Chem. Phys. 62:4189–4196.
  • 26. Frenkel, D. (2004) Speed-up of Monte Carlo simulations by sampling of rejected states. Proc. Natl. Acad. Sci. USA 101:17571–17575.
  • 27. Duane, S, Kennedy, A. D, Pendleton, B. J, & Roweth, D. (1987) Hybrid Monte Carlo. Phys. Lett. B 195:216.
  • 28. Lelièvre, T, Stoltz, G, & Rousset, M. (2010) Langevin dynamics with constraints and computation of free energy differences. ariv:1006.4914.
  • 29. Crooks, G. E. (1998) Nonequilibrium measurements of free energy differences for microscopically reversible markovian systems. J. Stat. Phys. 90:1481–1487.
  • 30. Jarzynski, C. (1997) Nonequilibrium equality for free energy differences. Phys. Rev. Lett. 78:2690–2693.
  • 31. Lechner, W, Oberhofer, H, Dellago, C, & Geissler, P. L. (2006) Equilibrium free energies from fast-switching trajectories with large time steps. J. Chem. Phys. 124:044113.
  • 32. Andersen, H. C. (1983) Rattle: A “velocity” version of the Shake algorithm for molecular dynamics calculations. J. Comput. Phys. 52:24–34.
  • 33. Leimkuhler, B & Skeel, R. D. (1994) Symplectic numerical integrators in constrained Hamiltonian systems. J. Comput. Phys. 112:117–125.
  • 34. Paterlini, M. G & Ferguson, D. M. (1998) Constant temperature simulations using the Langevin equation with velocity Verlet integration. Chem. Phys. 236:243–252.
  • 35. Izaguirre, J. A, Sweet, C. R, & Pande, V. S. (2010) Multiscale dynamics of macromolecules using normal mode Langevin. Pac. Symp. Biocomp. 15:240–251.
  • 36. Lelièvre, T, Stoltz, G, & Rousset, M. (2010) Free energy computations: A mathematical perspective. (Imperial College Press), 1st edition.
  • 37. Gustafson, P. (1998) A guided walk Metropolis algorithm. Stat. Comput. 8:357–364.
  • 38. Brünger, A, Brooks III, C. L, & Karplus, M. (1984) Stochastic boundary conditions for molecular dynamics simulations of ST2 water. Chem. Phys. Lett. 105:495–500.
  • 39. Pastor, R. W, Brooks, B. R, & Szabo, A. (1988) An analysis of the accuracy of Langevin and molecular dynamics algorithms. Mol. Phys. 65:1409–1419.
  • 40. Schlick, T. (2002) Molecular Modeling and Simulation: An Interdisciplinary Guide. (Springer-Verlag), 1st edition.
  • 41. Weeks, J. D, Chandler, D, & Andersen, H. C. (1971) Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 54:5237–5247.
  • 42. Chodera, J. D, Swope, W. C, Pitera, J. W, Seok, C, & Dill, K. A. (2007) Use of the weighted histogram analysis method for the analysis of simulated and parallel tempering simulations. J. Chem. Theor. Comput. 3:26–41.
  • 43. Bolhuis, P. G, Chandler, D, Dellago, C, & Geissler, P. L. (2002) Transition path sampling: Throwing ropes over rough mountain passes, in the dark. Annu. Rev. Phys. Chem. 53:291–318.
  • 44. Best, R. B & Hummer, G. (2005) Reaction coordinates and rates from transition paths. Proc. Natl. Acad. Sci. USA 102:6732–6737.
  • 45. Ma, A & Dinner, A. R. (2005) Automatic method for identifying reaction coordinates in complex systems. J. Phys. Chem. B 109:6769–6779.
  • 46. E, W, Ren, W, & Vanden-Eijnden, E. (2005) Finite temperature string method for the study of rare events. J. Phys. Chem. B 109:6688–6693.
  • 47. Berezhkovskii, A & Szabo, A. (2005) One-dimensional reaction coordinates for diffusive activated rate processes in many dimensions. J. Chem. Phys. 122:014503.
  • 48. Peters, B & Trout, B. L. (2006) Obtaining reaction coordinates by likelihood maximization. J. Chem. Phys. 125:054108.
  • 49. Ensing, B, de Vivo, M, Liu, Z, Moore, P, & Klein, M. L. (2006) Metadynamics as a tool for exploring free energy landscapes of chemical reactions. Acc. Chem. Res. 39:73–81.
  • 50. Schmiedl, T & Seifert, U. (2007) Optimal finite-time processes in stochastic thermodynamics. Phys. Rev. Lett. 98:108301.
  • 51. Then, H & Engel, A. (2008) Computing the optimal protocol for finite-time processes in stochastic thermodynamic. Phys. Rev. E 77:041105.
  • 52. Gomez-Marin, A, Schmiedl, T, & Seifert, U. (2008) Optimal protocols for minimal work processes in underdamped statistical thermodynamics. J. Chem. Phys. 129:024114.
  • 53. Salamon, P & Berry, R. S. (1983) Thermodynamic length and dissipated availability. Phys. Rev. Lett. 51:1127.
  • 54. Crooks, G. E. (2007) Measuring thermodynamic length. Phys. Rev. Lett. 99:100602.
  • 55. Feng, E. H & Crooks, G. E. (2009) Far-from-equilibrium measurements of thermodynamic length. Phys. Rev. E 79:012104.
  • 56. Minh, D. D. L & Chodera, J. D. (2011) Multiple-timeslice nonequilibrium estimators: Estimating equilibrium ensemble averages from nonequilibrium experiments using multiple timeslices, with application to thermodynamic length. J. Chem. Phys. 134:024111.
  • 57. Shenfeld, D. K, Xu, H, Eastwood, M. P, Dror, R. O, & Shaw, D. E. (2009) Minimizing thermodynamic length to select intermediate states for free-energy calculations and replica-exchange simulations. Phys. Rev. E 80:046705.
  • 58. Minh, D. D. L & Adib, A. B. (2008) Optimized free energies from bidirectional single-molecule force spectroscopy. Phys. Rev. Lett. 100:180602.
  • 59. Minh, D. D. L & Chodera, J. D. (2009) Optimal estimators and asymptotic variances for nonequilibrium path-ensemble averages. J. Chem. Phys. 131:134110.
  • 60. Crooks, G. E. (1999) Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences. Phys. Rev. E 60:2721–2726.
  • 61. Crooks, G. E. (2000) Path-ensemble averages in systems driven far from equilibrium. Phys. Rev. E 61:2361–2366.
  • 62. Athènes, M & Marinica, M.-C. (2010) Free energy reconstruction from steered dynamics without post-processing. J. Chem. Phys. 229:7129–7146.
  • 63. Friedrichs, M. S, Eastman, P, Vaidyanathan, V, Houston, M, LeGrand, S, Beberg, A. L, Ensign, D. L, Bruns, C. M, & Pande, V. S. (2009) Accelerating molecular dynamic simulation on graphics processing units. J. Comput. Chem. 30:864–872.
  • 64. Eastman, P & Pande, V. S. (2010) OpenMM: A hardware-independent framework for molecular simulations. Computing in Science and Engineering 12:34–39.
  • 65. Eastman, P & Pande, V. S. (2010) Efficient nonbonded interactions for molecular dynamics of a graphics processing unit. J. Comput. Chem. 31:1268.
  • 66. Bruns, C. M, Radmer, R. A, Chodera, J. D, & Pande, V. S. (2010) PyOpenMM. http://simtk.org/home/pyopenmm.

Appendix

\@xsect

Following the proof for GHMC in Ref. [36], here we show that NCMC preserves the equilibrium distribution. The expected acceptance rate for NCMC moves initiated from is,

(0)

Suppose that we have a variate drawn from the equilibrium distribution . The probability density of the next value in the chain, , has contributions from two scenarios: when the candidate variate is rejected and when it is accepted. The contribution from rejecting the candidate and flipping the momentum such that is,

(0)

The latter contribution from accepting the candidate such that is,

(0)
\color

black where is the probability of generating the trajectory-protocol pair from , and the pathwise detailed balance condition (Eq. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation) is used to produce the quantity in brackets. \colorblack

Taking the sum of Eqs. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation and Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation, we find that the equilibrium distribution is preserved,

(0)

By analogous reasoning, maintaining the momentum upon rejection, , and flipping it upon acceptance, will also preserve the equilibrium distribution.

\@xsect

A common integrator for Brownian dynamics (the overdamped regime of Langevin dynamics), in which only coordinates are explicitly integrated, is given by Ermak and Yeh [24, 25]. \colorblackIn our notation, where the perturbed coordinates are propagated by one step of the stochastic integrator to obtain , application of the propagation kernel can be written,

(0)
\color

black where is the particle mass, is the (potentially time-dependent) systematic force, and is an effective collision frequency or friction coefficient with units of inverse time. The noise history for each degree of freedom is a normal random variate with zero mean and variance , drawn from the distribution

(0)

In NCMC, every application of the propagation kernel produces a transition determined by the noise history variable , there is a corresponding that generates the opposite step, . By noting \colorblack

(0)
\color

black we can compute the relationship,

(0)

Then, the ratio of transition kernels appearing in Eq. Nonequilibrium candidate Monte Carlo: A new tool for efficient equilibrium simulation can be written in terms of noise history and the computed reverse noise history ,