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. Whereas 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.
AUTHOR SUMMARY
We expect that the theoretical framework developed in this study will extend the repertoire of useful Monte Carlo moves that can be applied to condensed-phase systems. We also anticipate NCMC to enable the sampling of widely disparate chemical states in expanded ensemble frameworks, which will be vital for exploring large chemical spaces, such as the space of small organic chemical compounds—a key aspect of molecular design problems.
The situation differs, however, when the dimer is placed in a dense solvent. As seen in the third panel of Fig. P1, transitions between compact and extended states are substantially hindered, and the autocorrelation time for MD simulation is increased to 299.8 iterations. Because the number of statistically uncorrelated configurations sampled in a simulation of 10,000 iterations is greatly reduced, the histogram from this simulation significantly differs from the true distribution (shown in red). Moreover, sampling is not improved by attempting the same instantaneous MC moves as in vacuum (fourth panel, Fig. P1). However, when these MC moves are divided into 2,048 steps in our new NCMC framework, the autocorrelation time is reduced to 4.0 iterations, and it is possible to accurately reconstruct the equilibrium distribution of dimer extensions (the bottom panel of Fig. P1).
As an illustrative example, consider a bistable dimer of “bonded” particles interacting via a double-well potential, with minima at r = r0 (the compact configuration, where the particles are in contact) and r = 2r0 (the extended configuration), and an interatomic potential barrier height of 5 kBT. If the dimer is simulated in vacuum, thermostatted MD is capable of sampling dimer configurations with an autocorrelation time τ—a measure of the simulation time required to generate a statistically uncorrelated configuration—of 59.2 iterations of 500 MD steps per iteration (see the top panel of Fig. P1). For this system, a logical MC move is to extend the dimer from a compact configuration, and vice versa. Indeed, attempting this move after each iteration of MD vastly decreases the autocorrelation time, as seen in the second panel of Fig. P1.
Fig. P1.
Trajectories of dimer system in vacuum and solvent. (Left) The dimer extension r 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 r. (Right) Histogram accumulated over trajectory (black), with true equilibrium distribution (red). Plot titles denote whether simulation was run in vacuum (vacuum) or dense Weeks–Chandler–Andersen fluid (solvent), and whether the simulation utilized only 500 MD steps per iteration (actually Generalized Hybrid Monte Carlo; see paper) or included instantaneous MC or NCMC moves following each iteration.
Our approach to enhanced sampling, called nonequilibrium candidate Monte Carlo (NCMC), is inspired by recently discovered fluctuation theorems connecting nonequilibrium driven processes to equilibrium properties (2, 3). The method has some precedents in the statistics and statistical physics literature. Notably, Athènes proposed a work-biased MC scheme for the computation of chemical potentials (4), Ballard and Jarzynski developed a framework for using nonequilibrium driven processes in replica exchange simulations (5), and Stern used nonequilibrium driven processes as MC moves to sample an equilibrium mixture of protonation states for constant-pH simulations in explicit solvent. (6). In our current study, we develop a simple theoretical framework that significantly generalizes these results in a rigorous manner.
Molecules are constantly in motion, dancing between microscopic configurations that can be very distinct. Computer simulations that accurately predict a system’s properties must account for this structural diversity. Unfortunately, the most popular simulation methods, molecular dynamics (MD) and Metropolis Monte Carlo (MC) (1), are best suited for sampling very similar conformations. Whereas small-time steps are necessary for the numerical stability of MD simulations, MC allows for large moves that can more quickly lead to distinct configurations. In condensed-phase systems, however, large MC moves are often useless because good moves for one component of the system (such as a biomolecule) generally lead to unfavorable overlaps with other components in a system (such as solvent). In this article, we introduce a new class of MC moves based on nonequilibrium dynamics: (i) a candidate configuration is generated (sampled) by actively driving a system out of equilibrium, and (ii) the candidate is accepted or rejected with criteria that keep the system in equilibrium. When applied to a simple but informative model—a bistable dimer in solvent—we find that our approach substantially accelerates sampling and thus expands the scope of MC simulations.
This article is a PNAS Direct Submission.
Data deposition: Simulation and analysis scripts, as well as simulation data generated for this paper, are available online at http://simtk.org/home/ncmc/; simtk.org is an online repository for exchanging scientific data produced by SimBios, a National Institutes of Health Center for Biomedical Computing.
N Metropolis, AW Rosenbluth, MN Rosenbluth, AH Teller, E Teller, Equation of state calculations by fast computing machines. J Chem Phys21, 1087–1092 (1953).
In 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 (MC) (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, 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. Whereas 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 nonlocal 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. 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.
Fig. 1.
Bistable dimer potential and instantaneous MC moves in WCA fluid. An extension move increases the dimer extension by Δr = +r0, whereas a compaction move decreases the dimer extension by Δr = -r0. 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.
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. In the case of the bistable dimer immersed in dense solvent, the extension (or 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–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) [although an inexact variant was proposed earlier (10)], and Nilmeier et al. (12) proposed the driving of a subset of degrees of freedom to enhance acceptance rates (using an approximate acceptance criterion). Nonequilibrium processes have also been used to generate configurations for parallel tempering simulations (13–15).
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.
Equilibrium and Expanded Thermodynamic Ensembles
For physical systems in equilibrium, the probability of observing a microstate is given by the Boltzmann distribution,
[1]
where x∈Γ 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 Zλ is a normalizing constant known as the partition function.
The reduced potential uλ(x) depends on the thermodynamic ensemble under consideration (18). For instance, in an isothermal-isobaric (NpT) ensemble, the reduced potential will assume the form,
[2]
which depends on the Hamiltonian H(x) (which may include an external biasing potential, and is presumed to be invariant under momentum inversion) and the system volume V(x). In this ensemble, the vector of controllable thermodynamic parameters λ ≡ {β,H,p} includes the inverse temperature β, the Hamiltonian H(x), and external pressure p. 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 (x,λ) in a weighted mixture of thermodynamic states,
[3]
where ωλ > 0 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 π(x,λ) = πλ(x). 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).
Nonequilibrium Candidate Monte Carlo
We first describe the general form of NCMC. At the start of an iteration, the current sample in the Markov chain, (x(n),λ(n)), which is assumed to be drawn from π(x,λ), is used to initialize a trajectory, (x0,λ0) = (x(n),λ(n)). A candidate configuration (xT,λT) 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 (x0,λ0). 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, (x(n+1),λ(n+1)), is the candidate, (xT,λT), or the momentum reversal of the initial sample, .
An NCMC move begins by selecting a protocol Λ from a set of possible protocols with probability P(Λ|x0,λ0), such that there exists a reverse protocol labeled as (to be defined momentarily) with . A protocol Λ specifies both a series of T perturbation kernels αt(x,y) and propagation kernels Kt(x,y), arranged in an alternating pattern such that Λ ≡ {α1,K1,α2,K2,…,αT,KT}. Both αt(x,y) and Kt(x,y) are conditional probabilities of y∈Γ given any x∈Γ, and must satisfy the requirement that if p(x,y) > 0, then p(y,x) > 0, for p substituted by αt and Kt.
Each perturbation kernel αt drives some or all of the degrees of freedom x 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 Kt propagates some or all of the coordinates of the system at fixed λt 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 t. 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 X ≡ (x0,x1,…,xT) is generated from an initial microstate x0 according to a protocol Λ can be illustrated by the scheme,
[4]
Application of the perturbation αt to xt-1 generates a perturbed configuration , which is then propagated by the kernel Kt to obtain xt.
The reverse protocol reverses the order in which the perturbation and propagation steps are applied, generating the time-reversed trajectory , where denotes x with inverted momenta,
[5]
The next step in NCMC is to accept or reject (xT,λT) as the next sample in the chain. To ensure that the stationary distribution π(x,λ) is preserved, we enforce a strict pathwise form of detailed balance*,
[6]
The quantity A(X|Λ) is the NCMC acceptance probability, whereas Π(X|x0,Λ) and denote the probability of generating trajectory X from initial configuration x0 using protocol Λ, or from final configuration with protocol , respectively,
[7]
[8]
Summation of Eq. 6 over all trajectories starting with x0 and ending with xT recovers the standard detailed balance condition (see Appendix for proof).
We define the ratio of proposal kernels as
[9]
and the ratio of propagation kernels as the exponentiated difference in forward and backward conditional path actions as
[10]
Using the above expressions and the momentum invariance property , we may write the ratio of acceptance probabilities as
[11]
where Δu(X|Λ) ≡ uT(xT) - u0(x0) is the energy difference. Eq. 11 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 A(X|Λ) that satisfy Eq. 11 are possible, including the well-known Metropolis–Hastings criterion (2, 3)
[12]
After generating (xT,λT) and evaluating A(X|Λ), we generate a uniform random variate U. If A(X|Λ) > U, then the candidate becomes the next value in the chain, (x(n+1),λ(n+1)) = (xT,λT). 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, (x(n+1),λ(n+1)) = (x0,λ0). We cannot, however, ignore the momentum flip completely; as explained in Appendix, this flip is necessary to preserve the equilibrium distribution.
We note that NCMC need not be used exclusively to sample from π(x,λ), 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.
Perturbation Kernels
A large variety of choices are available for the perturbation kernels αt(x,y). 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.
Stochastically Driven Degrees of Freedom.
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 angles and bonds fixed) according to some probability, such as the von Mises circular distribution centered on ϕ,
[13]
with I0(κ) denoting the modified Bessel function of order zero and κ > 0 taking the role of a dimensionless force constant. Because the stochastic perturbation is made in a non-Cartesian coordinate, a Jacobian J(ϕ) must be included to compute α(x,y) in Cartesian coordinates, resulting in the ratio,
[14]
where J(ϕ′) = J(ϕ) = 1 because the transformation (a rotation about a bond vector) preserves the Cartesian phase space volume.
Deterministically Driven Degrees of Freedom.
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 x → y, such that differs from x 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 ϕ(x) from its current value ϕ0 to a final value ϕT = ϕ0 + Δϕ over T steps in equal increments, such that ϕ(xt) is constrained to ϕt ≡ (1 - t/T)ϕ0 + (t/T)ϕT. In this case, , where the Jacobian J(x) 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 .
Simulation Box Scaling.
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 s = [(V(x) + ΔV)/V(x)]1/3 with ΔV chosen uniformly from [V - ΔV0,V + ΔV0] applied as a factor of s1/T over the course of T steps. In this case, the perturbation kernel αt(x,y) is a delta function that compresses or expands the molecular centers and box geometry. Because the Jacobian is the ratio of infinitesimal volumes upon scaling, the ratio of perturbation kernels is , where N denotes the number of molecular centers.
Thermodynamic Perturbation.
In many driven nonequilibrium processes, there is no direct perturbation to the coordinates, such that αt(x,y) = δ(x - y) and the ratio . Instead, only the thermodynamic parameters λ are varied in time, carrying the system out of equilibrium through action of the Kt propagation kernels. We recover Neal’s method (7) if the reduced potential ut is a simple linear interpolation such that ut(x) = (1 - t/T)u0(x) + (t/T)uT(x), the probability of choosing protocol Λ is symmetric with , and MC (2, 3) is used for the propagation kernel Kt.
Propagation Kernels
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 Kt(x,y) ≡ K(x,y) that samples from a stationary distribution π(x) of interest. Alternatively, a strongly time-dependent Kt 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.
Reversible Markov Chain Monte Carlo.
One may propagate some or all of a system’s degrees of freedom (e.g., those not affected by the perturbation kernel αt) by a method that satisfies detailed balance in πt,
[15]
where for a specified ut(x). 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 w and heat q for the nonequilibrium driven process,
[16]
[17]
such that w(X|Λ) + q(X|Λ) = Δu(X|Λ), 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, q(X|Λ),
[18]
leading to an acceptance probability similar to standard MC, except that the work, w(X|Λ), replaces the instantaneous potential energy difference,
[19]
Deterministic Dynamics.
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,
[20]
where Δu(X|Λ) ≡ uT(xT) - u0(x0) 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).
Stochastic Dynamics.
Stochastic integrators such as velocity Verlet discretizations of Langevin dynamics (34, 35) sample a modified distribution that differs from the desired equilibrium distribution πt in a timestep-dependent manner (36). Whereas this modified distribution may be difficult or impossible to compute 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 ensemble†. 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–40) in Appendix.
Illustrative Application: Bistable Dimer in a WCA Fluid
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 r = r0 (compact) and r = 2r0 (extended), and a 5kBT barrier (see Fig. 1). In the solvated simulations, the dimer was immersed in a dense bath (reduced density ρσ3 = 0.96) 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 (1, 28, 36, 37) (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.
The rate at which effectively uncorrelated samples are generated can be quantified in terms of the correlation time τ for the dimer extension r(t) (shown in Fig. 2). This time represents the asymptotic decay time constant for the correlation function , which will behave like
[21]
for large t, where and . The correlation time is related to the statistical inefficiency, g = 1 + 2τ, a factor that describes the number of iterations necessary to generate an effectively uncorrelated sample (42).
Fig. 2.
Trajectories of WCA dimer system in vacuum and solvent. (Left) The dimer extension r 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 r. (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 2,048-step NCMC moves (NCMC) following each iteration.
For the MD simulation in vacuum (Fig. 2, top trace), we observe slow hopping between compact and extended configurations with a correlation time τ = 59.2 iterations, resulting in slow convergence of the histogram. Introducing instantaneous MC dimer extension/contraction moves that modify the dimer extension by Δr = ± r0 reduces the correlation time to τ ≈ 0.0, 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).
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 2,048-step NCMC moves (Fig. 2, bottom trace) allows the system to rapidly mix between both compact and extended states with a correlation time of τ = 4.0 iterations. Whereas each iteration requires a fivefold increase in computational effort (500 MD steps + 2,048 NCMC switching steps = 2,548 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. Toward this end, we estimated the acceptance probability of the extension/contraction moves in dense solvent as a function of switching length (Fig. 3). Whereas instantaneous MC proposals of ± r0 are only accepted with probability ≈10-27 (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–8), there is little to no increase in acceptance rate, but for an intermediate number of steps (16–1,024), 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 2,000 steps, achieving an acceptance rate of 12% using nonequilibrium proposal trajectories of 2,048 steps, or 38% for 8,192 steps.
Fig. 3.
Acceptance probabilities of NCMC proposals. (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.
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 τeff, an indirect estimate of the correlation time,
[22]
where τMD and τNCMC 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 τNCMC ≈ -1/ ln(1 - 2γ) (see Appendix for derivation).
As shown in Fig. 4, the effective correlation time τeff is only diminished when the NCMC acceptance probability is large enough such that τNCMC ≈ τMD, which occurs when γ≥0.13% (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 τeff.
Fig. 4.
Statistical efficiency gain of NCMC proposals relative to instantaneous MC proposals. (Upper) Effective correlation time τeff, in iterations, for MD+NCMC (black dots) compared to MD alone (red line). (Lower) 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).
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,
[23]
Here, TMD = 500 steps per iteration, and TNCMC is again varied from 1–8,192 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 86.9% the efficiency of MD alone at 128 steps—followed by a rapid gain in efficiency, plateauing at an efficiency gain of approximately 13× the efficiency of MD alone for 2,048–4,096-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 8,192 steps, the reduction in correlation time is not sufficient to offset the additional cost of these moves.
Epilogue
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–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–52). Recent work has also suggested that estimating the thermodynamic metric tensor along the nonequilibrium parameter switching path (53–56), could prove useful in adaptively optimizing the switching protocol (57).
We note that switching trajectories contain potentially useful information. Indeed, several methods (56, 58, 59) 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
The dimer system considered here consists of two particles that interact via a double-well bonded potential in the interatomic distance r,
[24]
with h = 5kBT, r0 = rWCA, and s = rWCA/2, where rWCA ≡ 21/6σ. Simulations denoted as “vacuum” contain only these two particles, whereas simulations denoted as “solvated” also interact with a dense bath of particles via the WCA nonbonded potential (41),
[25]
with mass m = 39.9 amu, σ = 3.4 Å, and ϵ = 120 kBT. 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 ρσ3 = 0.96. For all simulations, the reduced temperature is kBT/ϵ = 0.824. A custom Python code making use of the graphics processing unit (GPU)-accelerated OpenMM package (63–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 GHMC (1, 28, 36, 37) 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 0.002τ, where , and the collision rate was set to τ-1. With this timestep, the acceptance probability is 99.929 ± 0.001%; the resulting dynamics closely approximates Langevin dynamics.
In simulations employing instantaneous Monte Carlo moves, a perturbation Δr to the interatomic distance r of the dimer was chosen according to
[26]
The dimer was contracted or expanded about the bond midpoint to generate a new configuration xnew with dimer extension rnew from the old configuration xold with dimer extension rold, and the move accepted or rejected with the Metropolis–Hastings criterion (3),
[27]
where the Jacobian ratio term Jr(xold,xnew) = (rnew/rold)2 accounts for the expansion and contraction of phase space due to the Monte Carlo proposals.
For simulations employing T-step NCMC moves, proposals were made by selecting a new velocity vector from the Maxwell–Boltzmann distribution, integrating T steps of velocity Verlet dynamics (23) for all bath atoms as the dimer extension was driven from rold to rnew in equal steps of size Δr/T, and accepting or rejecting based on the modified Metropolis criteria for symplectic integrators (Eqs. 12–20),
[28]
The Jacobian ratio is also (rnew/rold)2. This scheme corresponds to a choice for the perturbation kernel of
[29]
where r(x) denotes the dimer separation of configuration x. The propagation kernel Kt(x,y) corresponds to velocity Verlet dynamics where the dimer atoms are held fixed in space. 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
[30]
For numerical stability, logarithms of A(Xn) were stored, as an ≡ ln A(Xn). We then estimated (shown in Fig. 4) as
[31]
where b ≡ max nan.
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 2,048-step NCMC, with 500 steps of GHMC dynamics in between each NCMC trial, to ensure equilibrium sampling. Prior to each 2,048-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 1,000 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.
The reference distribution for the interparticle distribution plotted in red on the right side of Fig. 2 was computed analytically for the vacuum system,
[32]
For the solvated system, this distribution was estimated from an umbrella sampling simulation (Fig. 5) employing a modified bonded potential intended to remove the barrier in between compact and extended states,
[33]
where rmin = r0, rmax = 2.05r0, and K = kBT/η2, with η = 0.3 Å, and θ(r) is the Heaviside function that assumes a value of unity for r≥0 and zero otherwise. The true solvated interparticle distribution p(r) was estimated by reweighting the data produced from this simulation, using the relationship
[34]
where rn denotes the bond separation for sample n, and a finite-width histogram bin was used instead of the delta function δ(r).
Fig. 5.
Umbrella sampling simulation of the dimer in WCA solvent. (Left) The dimer extension r 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.
Appendix
Proof that NCMC Preserves the Equilibrium Distribution.
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 (x,λ) is
[35]
Suppose that we have a variate (x(n),λ(n)) drawn from the equilibrium distribution π(x,λ). The probability density of the next value in the chain, p(x(n+1),λ(n+1)), 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
[36]
The latter contribution from accepting the candidate such that (x(n+1),λ(n+1)) = (xT,λT) is,
[37]
where ρ(X,Λ|x0,λ0) ≡ Π(X|x0,Λ)P(Λ|x0,λ0) is the probability of generating the trajectory-protocol pair (X,Λ) from (x0,λ0), and the pathwise detailed balance condition (Eq. 6) is used to produce the quantity in brackets.
Taking the sum of Eqs. 36 and 37, we find that the equilibrium distribution is preserved
[38]
By analogous reasoning, maintaining the momentum upon rejection, (x(n+1),λ(n+1)) = (x(n),λ(n)), and flipping it upon acceptance, will also preserve the equilibrium distribution.
Acceptance Criteria for Overdamped Langevin (Brownian) Integrator of Ermak and Yeh.
A common integrator for Brownian dynamics (the overdamped regime of Langevin dynamics), in which only coordinates x are explicitly integrated, is given by Ermak and Yeh (24, 25). In our notation, where the perturbed coordinates are propagated by one step of the stochastic integrator to obtain xt, application of the propagation kernel can be written,
[39]
where m is the particle mass, Ft(x) ≡ -(∂/∂x)Ht(x) is the (potentially time-dependent) systematic force, and γ is an effective collision frequency or friction coefficient with units of inverse time. The noise history ξt for each degree of freedom is a normal random variate with zero mean and variance β-1, drawn from the distribution
[40]
In NCMC, every application of the propagation kernel Kt produces a transition determined by the noise history variable ξt, there is a corresponding that generates the opposite step, . By noting
[41]
we can compute the relationship
[42]
Then, the ratio of transition kernels appearing in Eq. 10 can be written in terms of noise history ξt and the computed reverse noise history
[43]
where the tildes are dropped because the microstate x contains no momenta. The quantity |∂xt/∂ξt| represents the Jacobian for the change of variables from the ξt to xt, and the Jacobians in the numerator and denominator cancel. The quantity in Eq. 43 can easily be computed during integration.
Acceptance Criteria for Langevin Integrator of Brooks, Brünger, and Karplus (BBK).
The Brünger–Brooks–Karplus (BBK) stochastic integrator (38, 39) is a popular integrator for simulating Langevin dynamics. In our notation, where the perturbed coordinates are propagated by one step of the stochastic integrator to obtain xt, application of the propagation kernel can be written
[44]
where we have used a velocity Verlet discretization of the BBK integrator (36, 40). Here rt and vt denote the respective Cartesian position and velocity components of the microstate xt, γ the effective collision frequency with units of inverse time, and m the particle mass. is an auxiliary variable used only in simplifying the mathematical representation of the integration scheme. Note that we require two random variates, ξt and , per degree of freedom per timestep in order for this scheme to be able to generate both the forward trajectory X and its time-reverse (e.g., section 2.2.3.2 of ref. 36).
The noise history terms ξt and are normal random variates with zero mean and variance β-1. Their joint distribution can therefore be written
[45]
For every step , the positions and velocities undergo a transition determined by the noise variables . A corresponding choice of noise variables will generate the reverse step, carrying With a little algebra, it is seen that
[46]
To write the ratio of transition kernels appearing in Eq. 10 in terms of noise variables and the computed reverse noise variables , we must first compute the Jacobian because the random variates are not in Cartesian space
[47]
which can be shown to be independent of ξt and . The conditional path action difference can now be computed
[48]
where the ratio of Jacobians cancels because the Jacobians are independent of the noise variates.
Derivation of Effective Correlation Time for Mixed MD/NCMC Sampling.
For simplicity, we make the assumption that the system of interest has two long-lived conformational states of equal population with dimer extensions rc and re. This assumption holds to good approximation for the WCA dimer example considered here, and may apply to other systems of interest as well. We assume that the correlation time for a fixed number of MD simulation steps per iteration is given by τMD, and describe the probability of finding the system ends up in a given conformational state after one iteration by a 2 × 2 column-stochastic transition matrix TMD
[49]
For a 2 × 2 system whose time evolution is governed by the column stochastic transition matrix T, we can write the autocorrelation function for the dimer extension r as
[50]
where the transition matrix T has unitary eigenvalue decomposition UMU-1, and μ is the nonunit eigenvalue of T. The constants are and C∞ = (1/4)(rc + re)2.
Relating this to the autocorrelation time τ estimated from a timeseries, intended to reflect the fit to
[51]
we can see that τ = -1/ ln μ. We then determine that the correlation time τMD = -1/ ln μMD, with μMD = 1 - 2α.
Similarly, we can write the probability that the NCMC step will carry the system from one conformational state to another in terms of the acceptance probability γ, which we assume to be symmetric,
[52]
where we have correlation time τNCMC = -1/ ln μNCMC and μNCMC = 1 - 2γ.
The effective transition matrix Teff for iterations consisting of MD simulation steps followed by an NCMC trial is given by
[53]
where the effective correlation time τeff = -1/ ln μeff, with μeff = 1 - 2[(α + γ) - 2αγ]. Substituting in α = (1 - e-1/τMD)/2 and γ = (1 - e-1/τNCMC)/2, we obtain
[54]
As a check of the accuracy of Eq. 54, we note that for MD with 2,048-step NCMC switching, we compute τeff ≈ 4.0 iterations, using only τMD = 299.8 iterations and the NCMC acceptance probability γ = 12.1%. The actual correlation time measured from a 10,000 iteration simulation is computed to be τeff = 4.0.
Note
The authors declare no conflict of interest.
Data Availability
Data deposition: Simulation and analysis scripts, as well as simulation data generated for this paper, are available online at http://simtk.org/home/ncmc/; simtk.org is an online repository for exchanging scientific data produced by SimBios, a National Institutes of Health Center for Biomedical Computing.
Acknowledgments.
We thank Gabriel Stoltz (CERMICS, École des Ponts ParisTech); Jed W. Pitera (IBM Almaden Research Center); Manuel Athenès (CEA Saclay); Firas 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, as well as the two anonymous referees for their helpful suggestions for improving clarity. J.P.N. was supported by the US Department of Energy (DOE) by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. G.E.C. 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 DOE under Contract DE-AC02-05CH11231. D.D.L.M. was funded by a Director’s Postdoctoral Fellowship from Argonne National Laboratory. J.D.C. 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 DOE by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and by University of Chicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a USDE Office of Science laboratory, is operated under Contract DE-AC02-06CH11357.
References
1
JS Liu Monte Carlo Strategies in Scientific Computing (Springer, 2nd Ed, New York, 2002).
N Metropolis, AW Rosenbluth, MN Rosenbluth, AH Teller, E Teller, Equation of state calculations by fast computing machines. J Chem Phys21, 1087–1092 (1953).
RM Neal, Taking bigger Metropolis steps by dragging fast variables. (Department of Statistics, University of Toronto, Toronto, Technical Report 0411. (2004).
R Bürgi, PA Kollman, WF van Gunsteren, Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation. Proteins47, 469–480 (2002).
AP Lyubartsev, AA Martsinovski, SV Shevkunov, PN Vorontsov-Velyaminov, New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles. J Chem Phys96, 1776–1783 (1992).
H Li, M Fajer, W Yang, 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 Phys126, 024106 (2007).
WC Swope, HC Andersen, PH Berens, KR Wilson, 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 Phys76, 637–649 (1982).
DL Ermak, Y Yeh, Equilibrium electrostatic effects on the behavior of polyions in solution: Polyion-mobile ion interaction. Chem Phys Lett24, 243–248 (1974).
T Lelièvre, G Stoltz, M Rousset, Langevin dynamics with constraints and computation of free energy differences. (Imperial College Press, London, 2010).
W Lechner, H Oberhofer, C Dellago, PL Geissler, Equilibrium free energies from fast-switching trajectories with large time steps. J Chem Phys124, 044113 (2006).
MG Paterlini, DM Ferguson, Constant temperature simulations using the Langevin equation with velocity Verlet integration. Chem Phys236, 243–252 (1998).
JD Weeks, D Chandler, HC Andersen, Role of repulsive forces in determining the equilibrium structure of simple liquids. J Chem Phys54, 5237–5247 (1971).
JD Chodera, WC Swope, JW Pitera, C Seok, KA Dill, Use of the weighted histogram analysis method for the analysis of simulated and parallel tempering simulations. J Chem Theor Comput3, 26–41 (2007).
PG Bolhuis, D Chandler, C Dellago, PL Geissler, Transition path sampling: Throwing ropes over rough mountain passes, in the dark. Annu Rev Phys Chem53, 291–318 (2002).
A Berezhkovskii, A Szabo, One-dimensional reaction coordinates for diffusive activated rate processes in many dimensions. J Chem Phys122, 014503 (2005).
B Ensing, M de Vivo, Z Liu, P Moore, ML Klein, Metadynamics as a tool for exploring free energy landscapes of chemical reactions. Acc Chem Res39, 73–81 (2006).
A Gomez-Marin, T Schmiedl, U Seifert, Optimal protocols for minimal work processes in underdamped statistical thermodynamics. J Chem Phys129, 024114 (2008).
DK Shenfeld, H Xu, MP Eastwood, RO Dror, DE Shaw, Minimizing thermodynamic length to select intermediate states for free-energy calculations and replica-exchange simulations. Phys Rev E80, 046705 (2009).
Data deposition: Simulation and analysis scripts, as well as simulation data generated for this paper, are available online at http://simtk.org/home/ncmc/; simtk.org is an online repository for exchanging scientific data produced by SimBios, a National Institutes of Health Center for Biomedical Computing.
We thank Gabriel Stoltz (CERMICS, École des Ponts ParisTech); Jed W. Pitera (IBM Almaden Research Center); Manuel Athenès (CEA Saclay); Firas 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, as well as the two anonymous referees for their helpful suggestions for improving clarity. J.P.N. was supported by the US Department of Energy (DOE) by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. G.E.C. 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 DOE under Contract DE-AC02-05CH11231. D.D.L.M. was funded by a Director’s Postdoctoral Fellowship from Argonne National Laboratory. J.D.C. 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 DOE by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and by University of Chicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a USDE Office of Science laboratory, is operated under Contract DE-AC02-06CH11357.
Notes
This article is a PNAS Direct Submission.
Data deposition: Simulation and analysis scripts, as well as simulation data generated for this paper, are available online at http://simtk.org/home/ncmc/; simtk.org is an online repository for exchanging scientific data produced by SimBios, a National Institutes of Health Center for Biomedical Computing.
See Author Summary on page 18199 (volume 108, number 45).
*The described pathwise detailed balance condition is closely related to “super-detailed balance” (e.g., ref. 26), but also accounts for momentum reversal to extend its definition to include molecular dynamics integrators.
†
Alternatives 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 (1, 28, 36, 37).
Authors
Affiliations
Jerome P.Nilmeier
Biosciences and Biotechnology Division, Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, CA 94550;
To whom correspondence should be addressed. E-mail: [email protected].
Author contributions: J.P.N., G.E.C., D.D.L.M., and J.D.C. designed research; D.D.L.M. and J.D.C. performed research; J.P.N., G.E.C., D.D.L.M., and J.D.C. contributed new analytic tools; J.D.C. analyzed data; and J.P.N., G.E.C., D.D.L.M., and J.D.C. wrote the paper.
Competing Interests
The authors declare no conflict of interest.
Metrics & Citations
Metrics
Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.
Citation statements
Altmetrics
Citations
Cite this article
J.P. Nilmeier, G.E. Crooks, D.D.L. Minh, J.D. Chodera, Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation, Proc. Natl. Acad. Sci. U.S.A.
Collective cognition is often mentioned as one of the advantages of group living.
But which factors actually facilitate group smarts? To answer this, we compared how
individuals and groups of either ants or people tackle an identical ...Biological ensembles use collective intelligence to tackle challenges together, but
suboptimal coordination can undermine the effectiveness of group cognition. Testing
whether collective cognition exceeds that of the individual is often impractical since
...
US municipal wastewater facilities are major per- and polyfluoroalkyl substances (PFAS)
sources known to affect drinking water quality. Among eight large wastewater treatment
facilities with comparable treatment technologies and sizes to those ...Wastewater receives per- and polyfluoroalkyl substances (PFAS) from diverse consumer
and industrial sources, and discharges are known to be a concern for drinking water
quality. The PFAS family includes thousands of potential chemical structures ...
One of the largest volcanic eruptions of the nineteenth century took place in 1831
CE. Although this event led to significant Northern Hemisphere climate cooling, the
source of this eruption remains a mystery. Using evidence from well-dated ...Polar ice cores and historical records evidence a large-magnitude volcanic eruption
in 1831 CE. This event was estimated to have injected ~13 Tg of sulfur (S) into the
stratosphere which produced various atmospheric optical phenomena and led to Northern
...
where x∈Γ 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 Zλ is a normalizing constant known as the partition function.
which depends on the Hamiltonian H(x) (which may include an external biasing potential, and is presumed to be invariant under momentum inversion) and the system volume V(x). In this ensemble, the vector of controllable thermodynamic parameters λ ≡ {β,H,p} includes the inverse temperature β, the Hamiltonian H(x), and external pressure p. Other thermodynamic parameters and their conjugate variables can be included or excluded to generate alternative physical (or unphysical) ensembles.
where ωλ > 0 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 π(x,λ) = πλ(x). 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).
.
(to be defined momentarily) with 
. A protocol Λ specifies both a series of T perturbation kernels αt(x,y) and propagation kernels Kt(x,y), arranged in an alternating pattern such that Λ ≡ {α1,K1,α2,K2,…,αT,KT}. Both αt(x,y) and Kt(x,y) are conditional probabilities of y∈Γ given any x∈Γ, and must satisfy the requirement that if p(x,y) > 0, then p(y,x) > 0, for p substituted by αt and Kt.
Application of the perturbation αt to xt-1 generates a perturbed configuration 
, which is then propagated by the kernel Kt to obtain xt.
reverses the order in which the perturbation and propagation steps are applied, generating the time-reversed trajectory 
, where 
denotes x with inverted momenta,
The quantity A(X|Λ) is the NCMC acceptance probability, whereas Π(X|x0,Λ) and 
denote the probability of generating trajectory X from initial configuration x0 using protocol Λ, or 
from final configuration 
with protocol 
, respectively,
Summation of Eq. 6 over all trajectories starting with x0 and ending with xT recovers the standard detailed balance condition (see Appendix for proof).
and the ratio of propagation kernels as the exponentiated difference in forward and backward conditional path actions as
Using the above expressions and the momentum invariance property 
, we may write the ratio of acceptance probabilities as
where Δu(X|Λ) ≡ uT(xT) - u0(x0) is the energy difference. Eq. 11 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.
After generating (xT,λT) and evaluating A(X|Λ), we generate a uniform random variate U. If A(X|Λ) > U, then the candidate becomes the next value in the chain, (x(n+1),λ(n+1)) = (xT,λT). 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, (x(n+1),λ(n+1)) = (x0,λ0). We cannot, however, ignore the momentum flip completely; as explained in Appendix, this flip is necessary to preserve the equilibrium distribution.
with I0(κ) denoting the modified Bessel function of order zero and κ > 0 taking the role of a dimensionless force constant. Because the stochastic perturbation is made in a non-Cartesian coordinate, a Jacobian J(ϕ) must be included to compute α(x,y) in Cartesian coordinates, resulting in the ratio,
where J(ϕ′) = J(ϕ) = 1 because the transformation (a rotation about a bond vector) preserves the Cartesian phase space volume.
that takes x → y, such that 
differs from x 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 ϕ(x) from its current value ϕ0 to a final value ϕT = ϕ0 + Δϕ over T steps in equal increments, such that ϕ(xt) is constrained to ϕt ≡ (1 - t/T)ϕ0 + (t/T)ϕT. In this case, 
, where the Jacobian J(x) 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 
.
, where N denotes the number of molecular centers.
. Instead, only the thermodynamic parameters λ are varied in time, carrying the system out of equilibrium through action of the Kt propagation kernels. We recover Neal’s method (7) if the reduced potential ut is a simple linear interpolation such that ut(x) = (1 - t/T)u0(x) + (t/T)uT(x), the probability of choosing protocol Λ is symmetric with 
, and MC (2, 3) is used for the propagation kernel Kt.
where 
for a specified ut(x). 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.
such that w(X|Λ) + q(X|Λ) = Δu(X|Λ), a restatement of the first law of thermodynamics.
leading to an acceptance probability similar to standard MC, except that the work, w(X|Λ), replaces the instantaneous potential energy difference,
. Hence, 
and the acceptance ratio is,
where Δu(X|Λ) ≡ uT(xT) - u0(x0) 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).
is relatively straightforward, and the NCMC acceptance criteria ensures that the NCMC-sampled configurations are distributed according to the desired equilibrium ensemble†. 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–40) in Appendix.
, which will behave like
for large t, where 
and 
. The correlation time is related to the statistical inefficiency, g = 1 + 2τ, a factor that describes the number of iterations necessary to generate an effectively uncorrelated sample (42).
where τMD and τNCMC 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 τNCMC ≈ -1/ ln(1 - 2γ) (see Appendix for derivation).
Here, TMD = 500 steps per iteration, and TNCMC is again varied from 1–8,192 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 86.9% the efficiency of MD alone at 128 steps—followed by a rapid gain in efficiency, plateauing at an efficiency gain of approximately 13× the efficiency of MD alone for 2,048–4,096-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 8,192 steps, the reduction in correlation time is not sufficient to offset the additional cost of these moves.
with h = 5kBT, r0 = rWCA, and s = rWCA/2, where rWCA ≡ 21/6σ. Simulations denoted as “vacuum” contain only these two particles, whereas simulations denoted as “solvated” also interact with a dense bath of particles via the WCA nonbonded potential (41),
with mass m = 39.9 amu, σ = 3.4 Å, and ϵ = 120 kBT. 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 ρσ3 = 0.96. For all simulations, the reduced temperature is kBT/ϵ = 0.824. A custom Python code making use of the graphics processing unit (GPU)-accelerated OpenMM package (63–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.
, and the collision rate was set to τ-1. With this timestep, the acceptance probability is 99.929 ± 0.001%; the resulting dynamics closely approximates Langevin dynamics.
The dimer was contracted or expanded about the bond midpoint to generate a new configuration xnew with dimer extension rnew from the old configuration xold with dimer extension rold, and the move accepted or rejected with the Metropolis–Hastings criterion (3),
where the Jacobian ratio term Jr(xold,xnew) = (rnew/rold)2 accounts for the expansion and contraction of phase space due to the Monte Carlo proposals.
The Jacobian ratio is also (rnew/rold)2. This scheme corresponds to a choice for the perturbation kernel of
where r(x) denotes the dimer separation of configuration x. The propagation kernel Kt(x,y) corresponds to velocity Verlet dynamics where the dimer atoms are held fixed in space. The final configuration after the MC or NCMC rejection procedure was stored and plotted to generate Fig. 2.
For numerical stability, logarithms of A(Xn) were stored, as an ≡ ln A(Xn). We then estimated 
(shown in Fig. 4) as
where b ≡ max nan.
and the computed relative efficiency over instantaneous Monte Carlo was estimated by 1,000 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.
plotted in red on the right side of Fig. 2 was computed analytically for the vacuum system,
For the solvated system, this distribution was estimated from an umbrella sampling simulation (Fig. 5) employing a modified bonded potential intended to remove the barrier in between compact and extended states,
where rmin = r0, rmax = 2.05r0, and K = kBT/η2, with η = 0.3 Å, and θ(r) is the Heaviside function that assumes a value of unity for r≥0 and zero otherwise. The true solvated interparticle distribution p(r) was estimated by reweighting the data produced from this simulation, using the relationship
where rn denotes the bond separation for sample n, and a finite-width histogram bin was used instead of the delta function δ(r).
is
The latter contribution from accepting the candidate such that (x(n+1),λ(n+1)) = (xT,λT) is,
where ρ(X,Λ|x0,λ0) ≡ Π(X|x0,Λ)P(Λ|x0,λ0) is the probability of generating the trajectory-protocol pair (X,Λ) from (x0,λ0), and the pathwise detailed balance condition (Eq. 6) is used to produce the quantity in brackets.
will also preserve the equilibrium distribution.
are propagated by one step of the stochastic integrator to obtain xt, application of the propagation kernel 
can be written,
where m is the particle mass, Ft(x) ≡ -(∂/∂x)Ht(x) is the (potentially time-dependent) systematic force, and γ is an effective collision frequency or friction coefficient with units of inverse time. The noise history ξt for each degree of freedom is a normal random variate with zero mean and variance β-1, drawn from the distribution
determined by the noise history variable ξt, there is a corresponding 
that generates the opposite step, 
. By noting
we can compute the relationship
where the tildes are dropped because the microstate x contains no momenta. The quantity |∂xt/∂ξt| represents the Jacobian for the change of variables from the ξt to xt, and the Jacobians in the numerator and denominator cancel. The quantity in Eq. 43 can easily be computed during integration.
are propagated by one step of the stochastic integrator to obtain xt, application of the propagation kernel 
can be written
where we have used a velocity Verlet discretization of the BBK integrator (36, 40). Here rt and vt denote the respective Cartesian position and velocity components of the microstate xt, γ the effective collision frequency with units of inverse time, and m the particle mass. 
is an auxiliary variable used only in simplifying the mathematical representation of the integration scheme. Note that we require two random variates, ξt and 
, per degree of freedom per timestep in order for this scheme to be able to generate both the forward trajectory X and its time-reverse 
(e.g., section 2.2.3.2 of ref. 36).
are normal random variates with zero mean and variance β-1. Their joint distribution can therefore be written
For every step 
, the positions and velocities undergo a transition 
determined by the noise variables 
. A corresponding choice of noise variables 
will generate the reverse step, 
carrying 
With a little algebra, it is seen that
To write the ratio of transition kernels appearing in Eq. 10 in terms of noise variables 
and the computed reverse noise variables 
, we must first compute the Jacobian 
because the random variates are not in Cartesian space
which can be shown to be independent of ξt and 
. The conditional path action difference can now be computed
where the ratio of Jacobians 
cancels because the Jacobians are independent of the noise variates.
where the transition matrix T has unitary eigenvalue decomposition UMU-1, and μ is the nonunit eigenvalue of T. The constants are 
and C∞ = (1/4)(rc + re)2.
we can see that τ = -1/ ln μ. We then determine that the correlation time τMD = -1/ ln μMD, with μMD = 1 - 2α.
where we have correlation time τNCMC = -1/ ln μNCMC and μNCMC = 1 - 2γ.
where the effective correlation time τeff = -1/ ln μeff, with μeff = 1 - 2[(α + γ) - 2αγ]. Substituting in α = (1 - e-1/τMD)/2 and γ = (1 - e-1/τNCMC)/2, we obtain
As a check of the accuracy of Eq. 54, we note that for MD with 2,048-step NCMC switching, we compute τeff ≈ 4.0 iterations, using only τMD = 299.8 iterations and the NCMC acceptance probability γ = 12.1%. The actual correlation time measured from a 10,000 iteration simulation is computed to be τeff = 4.0.