<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">% PRL look and style (easy on the eyes)
\documentclass[aps,prl,twocolumn,superscriptaddress,nofootinbib]{revtex4-1}
% Two-column style (for submission/review/editing)
%\documentclass[aps,preprint,prl,superscriptaddress,showpacs]{revtex4}

\usepackage{palatino}
\usepackage{amsmath}
\usepackage{amssymb}

\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{boxedminipage}


\usepackage{verbatim}
\usepackage[colorlinks=true,citecolor=blue,linkcolor=blue]{hyperref}

% The figures are in a figures/ subdirectory.
\graphicspath{{figures/}}

% italicized boldface for math (e.g. vectors)
\newcommand{\bfv}[1]{{\mbox{\boldmath{$#1$}}}}
% non-italicized boldface for math (e.g. matrices)
\newcommand{\bfm}[1]{{\bf #1}}

%\newcommand{\bfm}[1]{{\mbox{\boldmath{$#1$}}}}
%\newcommand{\bfm}[1]{{\bf #1}}
\newcommand{\expect}[1]{\left \langle #1 \right \rangle}                % &lt;.&gt; for denoting expectations over realizations of an experiment or thermal averages

% vectors
\newcommand{\x}{\bfv{x}}
\newcommand{\y}{\bfv{y}}
\newcommand{\f}{\bfv{f}}

\newcommand{\bfc}{\bfm{c}}
\newcommand{\hatf}{\hat{f}}

\newcommand{\bTheta}{\bfm{\Lambda}}
\newcommand{\btheta}{\bfm{\theta}}
\newcommand{\bhatf}{\bfm{\hat{f}}}
\newcommand{\Cov}[1] {\mathrm{cov}\left( #1 \right)}
\newcommand{\Ept}[1] {{\mathrm E}\left[ #1 \right]}
\newcommand{\Eptk}[2] {{\mathrm E}_{#1}\left[ #2\right]}
\newcommand{\T}{\mathrm{T}}                                % T used in matrix transpose

% Commands to force LaTeX to put two figures on a page:
\renewcommand{\textfraction}{0.05}
\renewcommand{\topfraction}{0.95}
\renewcommand{\bottomfraction}{0.95}
\renewcommand{\floatpagefraction}{0.35}
\setcounter{totalnumber}{5}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DOCUMENT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TITLE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Nonequilibrium candidate Monte Carlo:\\
A new tool for efficient equilibrium simulation}

\author{Jerome P. Nilmeier}
\email[Electronic Address: ]{nilmeier1@llnl.gov}
%\affiliation{Postdoctoral Fellow at Lawrence Livermore National Laboratory, Physical and Life Sciences Directorate, Livermore, California 94550, USA}
\affiliation{Biosciences and Biotechnology Division, Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California 94550, USA}

\author{Gavin E. Crooks}
\email[Electronic Address: ]{gecrooks@lbl.gov}
\affiliation{Physical Biosciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA}

\author{David D. L. Minh}
\email[Electronic Address: ]{daveminh@anl.gov}
\affiliation{Biosciences Division, Argonne National Laboratory, Argonne, Illinois 60439, USA}

\author{John D. Chodera}
\thanks{Corresponding author}
\email[Electronic Address: ]{jchodera@berkeley.edu}
\affiliation{California Institute of Quantitative Biosciences (QB3), University of California, Berkeley, 260J Stanley Hall, Berkeley, California 94720, USA}

\date{\today}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{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 \emph{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.

\emph{Keywords: Metropolis-Hastings; Markov chain Monte Carlo; molecular dynamics; expanded ensembles}
\end{abstract}

\maketitle

In this paper, we describe a new technique for constructing efficient Markov chain Monte Carlo (MCMC)~\cite{jun-s-liu:MCMC} 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~\cite{metropolis:jcp:1953:metropolis-monte-carlo, hastings:biometrika:1970:metropolis-hastings} sampling procedure is generalized by using \emph{nonequilibrium} processes to generate candidates for \emph{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)~\cite{metropolis:jcp:1953:metropolis-monte-carlo, hastings:biometrika:1970:metropolis-hastings} and molecular dynamics (MD)~\cite{rahman:phys-rev:1964:molecular-dynamics} 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~\cite{swendsen-wang:prl:1987:swendsen-wang,wolff:prl:1989:cluster-moves}.

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.

One approach that can allow unperturbed degrees of freedom to relax, and hence maintain a reasonable acceptance rate, is to use a \emph{nonequilibrium} process to generate candidate configurations.
Using the appropriate acceptance criterion for the final configuration will preserve the equilibrium distribution.
This idea was suggested in a statistical context by Neal~\cite{neal:2004:tech-report:dragging-fast-variables} and first applied to a chemical system by Stern~\cite{stern:jcp:2007:constant-pH}, who used it to sample an equilibrium mixture of protonation states at constant pH in explicit solvent (an inexact variant was proposed earlier~\cite{buergi:proteins:2002:constant-pH}).
Nilmeier et al.~\cite{nilmeier-jacobson:2009:jctc:posh} proposed the driving of a subset of degrees of freedom to enhance acceptance rates, but derived only an approximate acceptance criteria.
Nonequilibrium processes have also been used to generate configurations for parallel tempering simulations ~\cite{opps-schofield:pre:2001:extended-state-space,brown-head-gordon:jcc:2003:cool-walk,jarzynski:pnas:2009:nonequilibrium-swaps}.
%, and Wagoner and Pande~\cite{wagoner-pande:2010:hybrid-solvent-model} used nonequilibrium insertion/deletion moves to improve the efficiency of semigrand canonical sampling.
% Siepmann and Frenkel~\cite{Siepmann-Frenkel:molphys:1992:configurational-bias}, Brown and Head-Gordon~\cite{brown-head-gordon:jcc:2003:cool-walk}

Here, we unify these ideas and significantly extend the application of nonequilibrium moves in physical simulations.
We present a theoretical framework, \emph{nonequilibrium candidate Monte Carlo} (NCMC), that is applicable to both single thermodynamic states (e.g.~NVT, NpT, $\mu$VT ensembles) as well as \emph{mixtures} of thermodynamic states (e.g.~expanded ensemble~\cite{lyubartsev:jcp:1992:expanded-ensembles,park:pre:2008:simulated-tempering} simulations).
Nonequilibrium proposals may drive a subset of degrees of freedom, the thermodynamic parameters characterizing the equilibrium distribution, or \emph{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 ENSEMBLES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Equilibrium and Expanded Thermodynamic Ensembles}

For physical systems in equilibrium, the probability of observing a microstate is given by the Boltzmann distribution,
\begin{eqnarray}
\pi_\lambda(x) &amp;=&amp; Z_{\lambda}^{-1} \, e^{-u_\lambda(x)} \:\:;\:\: Z_{\lambda}=\int_\Gamma dx \, e^{-u_\lambda(x)}, \label{equation:fixed lambda_probability}
\end{eqnarray}
where $x \in \Gamma$ denotes a microstate of the system (which may include coordinates, momenta, and other dynamical variables, such as simulation box dimensions), $\lambda$ denotes a set of thermodynamic parameters whose values define a \emph{thermodynamic state}, and $Z_{\lambda}$ is a normalizing constant known as the \emph{partition function}.

The reduced potential $u_\lambda(x)$ depends on the thermodynamic ensemble under consideration~\cite{shirts-chodera:jcp:2008:mbar}.
For instance, in an isothermal-isobaric ($NpT$) ensemble, the reduced potential will assume the form,
\begin{eqnarray}
u_\lambda(x) &amp;=&amp; \beta [ H(x) + p V(x) ], \label{equation:reduced-potential}
\end{eqnarray}
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 $\lambda \equiv \{\beta, H, p\}$ includes the inverse temperature $\beta$, 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 \emph{multiple} thermodynamic states within a single simulation, we also define an \emph{expanded ensemble}~\cite{lyubartsev:jcp:1992:expanded-ensembles,park:pre:2008:simulated-tempering}, which specifies a joint distribution for $(x,\lambda)$ in a weighted mixture of thermodynamic states,
\begin{eqnarray}
\pi(x,\lambda) &amp;=&amp; \frac{Z_\lambda \pi_\lambda(x) \, \omega_\lambda}{\sum\limits_{\nu \in \mathcal{G}} \int_\Gamma dy \, Z_\nu \pi_\nu(y) \, \omega_\nu}, \label{equation:expanded-ensemble}
\end{eqnarray}
where $\omega_\lambda &gt; 0$ specifies an externally-imposed weight for state $\lambda$.
Here, $\lambda \in \mathcal{G}$ may assume values in a discrete or continuous space $\mathcal{G}$.
If the set $\mathcal{G}$ consists of a single value of $\lambda$, a single thermodynamic state is sampled, and $\pi(x,\lambda) = \pi_\lambda(x)$. 
These thermodynamic states may correspond to a variety of different states of interest, such as temperatures in a simulated tempering simulation~\cite{marinari-parisi:europhys-lett:1992:simulated-tempering}, alchemical states in a simulated scaling simulation~\cite{yang:jcp:2007:simulated-scaling}, or protonation states in a constant-pH simulation~\cite{mongan:j-comput-chem:2004:constant-ph}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACCEPTANCE CRITERIA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{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)},\lambda^{(n)})$, which is assumed to be drawn from $\pi(x,\lambda)$, is used to initialize a trajectory, $(x_0, \lambda_0) = (x^{(n)},\lambda^{(n)})$.
A \emph{candidate} configuration $(x_T,\lambda_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~\cite{jarsynzki:j-stat-phys:2000:fluctuation-theorem} selected with a probability dependent only on $(x_0,\lambda_0)$.
Even if we only wish to sample from a single thermodynamic state $\lambda$, we may use a protocol that transiently drives the thermodynamic parameters away from $\lambda$ and back again.
Finally, an acceptance probability is computed and used to decide whether the next sample in the Markov chain, $(x^{(n+1)},\lambda^{(n+1)})$, is the candidate, $(x_T,\lambda_T)$, or the momentum reversal of the initial sample, $(\tilde{x}^{(n)},\lambda^{(n)})$.

Our nonequilibrium move begins by choosing a protocol $\Lambda$ from a set of possible protocols with probability $P(\Lambda | x_0, \lambda_0)$, such that there exists a reverse protocol labeled as $\tilde{\Lambda}$ (to be defined momentarily) with $P(\tilde{\Lambda} | \tilde{x}_T, \lambda_T) &gt; 0$.
A protocol $\Lambda$ specifies both a series of $T$ \emph{perturbation kernels} $\alpha_t(x,y)$ and \emph{propagation kernels} $K_t(x,y)$, arranged in an alternating pattern such that $\Lambda \equiv \{ \alpha_1, K_1, \alpha_2, K_2, \ldots, \alpha_T, K_T \}$.
Both $\alpha_t(x,y)$ and $K_t(x,y)$ are conditional probabilities of $y \in \Gamma$ given any $x \in \Gamma$, and must satisfy the requirement that if $p(x,y) &gt; 0$, then $p(y,x) &gt; 0$, for $p$ substituted by $\alpha_t$ and $K_t$.

Each perturbation kernel $\alpha_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 $K_t$ propagates some or all of the coordinates of the system at fixed $\lambda_t$ according to some form of MCMC or MD (e.g.~Metropolis Monte Carlo~\cite{metropolis:jcp:1953:metropolis-monte-carlo,hastings:biometrika:1970:metropolis-hastings}, velocity Verlet~\cite{swope:jcp:1982:velocity-verlet} deterministic dynamics, or overdamped Langevin stochastic dynamics~\cite{ermak-yeh:cpl:1974:brownian-dynamics,ermak:jcp:1975:brownian-dynamics}) 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 \equiv (x_0, x_1, \ldots, x_T)$ is generated from an initial microstate $x_0$ according to a protocol $\Lambda$ can be illustrated by the scheme,
\begin{eqnarray}
x_0 \stackrel{\alpha_1}{\longrightarrow} x_1^* \stackrel{K_1}{\longrightarrow} x_1 \rightarrow \cdots \longrightarrow x_{T-1} \stackrel{\alpha_T}{\longrightarrow} x_T^* \stackrel{K_T}{\longrightarrow} x_T
\end{eqnarray}
Application of the perturbation $\alpha_t$ to $x_{t-1}$ generates a perturbed configuration $x_t^*$, which is then propagated by the kernel $K_t$ to obtain $x_t$.

The reverse protocol $\tilde{\Lambda} \equiv \{K_T, \alpha_T, \ldots, K_0, \alpha_0\}$ reverses the order in which the perturbation and propagation steps are applied, generating the time-reversed trajectory $\tilde{X} \equiv \{\tilde{x}_T, \ldots, \tilde{x}_0\}$, where $\tilde{x}$ denotes $x$ with inverted momenta,
\begin{eqnarray}
\tilde{x}_T \stackrel{K_T}{\longrightarrow} \tilde{x}_T^* \stackrel{\alpha_T}{\longrightarrow} \tilde{x}_{T-1} \rightarrow \cdots \rightarrow \tilde{x}_1 \stackrel{K_1}{\longrightarrow} \tilde{x}_1^*  \stackrel{\alpha_1}{\longrightarrow} \tilde{x}_0 
\end{eqnarray}

The next step in NCMC is to accept or reject $(x_T,\lambda_T)$ as the next sample in the chain.  
To ensure that the stationary distribution $\pi(x,\lambda)$ is preserved, we enforce a strict \emph{pathwise} form of detailed balance (see Appendix for proof),
\begin{eqnarray}
\lefteqn{A(X | \Lambda) \, \Pi(X | x_0, \Lambda) \, P(\Lambda | x_0, \lambda_0) \, \pi(x_0, \lambda_0)} \nonumber \\
&amp;=&amp; A(\tilde{X} | \tilde{\Lambda}) \, \Pi(\tilde{X} | \tilde{x}_T, \tilde{\Lambda}) \, P(\tilde{\Lambda} | \tilde{x}_T, \lambda_T) \, \pi(\tilde{x}_T, \lambda_T) . \label{equation:pathwise-detailed-balance}
\end{eqnarray}
The quantity $A(X | \Lambda)$ is the NCMC acceptance probability, while $\Pi(X | x_0, \Lambda)$ and $\Pi(\tilde{X} | \tilde{x}_T, \tilde{\Lambda})$ denote the probability of generating trajectory $X$ from initial configuration $x_0$ using protocol $\Lambda$, or $\tilde{X}$ from final configuration $\tilde{x}_T$ with protocol $\tilde{\Lambda}$, respectively,
\begin{eqnarray}
\Pi(X | x_0) &amp;=&amp; \prod_{1 \le t \le T} \alpha_t(x_{t-1}, x_t^*) \, K_t(x_t^*, x_t) \\
\Pi(\tilde{X} | \tilde{x}_T) &amp;=&amp; \prod_{T \ge t \ge 1} \alpha_t(\tilde{x}_t^*, \tilde{x}_{t-1}) \, K_t(\tilde{x}_t, \tilde{x}_t^*) .
\label{equation:forward-transition-probability}
\end{eqnarray}

We define the ratio of proposal kernels as,
\begin{eqnarray}
\frac{\alpha(\tilde{X}|\tilde{\Lambda})}{\alpha(X|\Lambda)} &amp;\equiv&amp; \prod_{t=1}^{T} \frac{\alpha_t(\tilde{x}_t^*, \tilde{x}_{t-1})}{\alpha_t(x_{t-1}, x_t^*)} \label{equation:ratio-of-perturbation-kernels},
\end{eqnarray}
and the ratio of propagation kernels as the exponentiated difference in forward and backward conditional path actions as,
\begin{eqnarray}
e^{-\Delta \mathcal{S}(X|\Lambda)} &amp;\equiv&amp; \prod_{t=1}^T \frac{K_t(\tilde{x}_t, \tilde{x}_t^*)}{K_t(x_t^*, x_t)} . \label{equation:ratio-of-propagation-kernels}
\end{eqnarray}
Using the above expressions and the momentum invariance property $\pi(x,\lambda) = \pi(\tilde{x},\lambda)$, we may write the ratio of acceptance probabilities as,
\begin{eqnarray}
\lefteqn{\frac{A(X|\Lambda)}{A(\tilde{X}|\tilde{\Lambda})} =\frac{\pi(\tilde{x}_T, \lambda_T)}{\pi(x_0, \lambda_0)} \frac{P(\tilde{\Lambda} | \tilde{x}_T, \lambda_T)}{P(\Lambda|x_0,\lambda_0)} \frac{\Pi(\Tilde{X}|\tilde{x}_T)}{\Pi(X | x_0)}} \nonumber \\
&amp;=&amp;\frac{\pi(x_T,\lambda_T)}{\pi(x_0,\lambda_0)} \frac{P(\tilde{\Lambda} | \tilde{x}_T, \lambda_T)}{P(\Lambda|x_0,\lambda_0)} \prod_{t=1}^{T} \frac{\alpha_t(\tilde{x}_t^*, \tilde{x}_{t-1})}{\alpha_t(x_{t-1}, x_t^*)} \frac{K_t(\tilde{x}_t, \tilde{x}_t^*)}{K_t(x_t^*, x_t)} \nonumber \\
&amp;\equiv&amp; \frac{\omega_T}{\omega_0} \frac{P(\tilde{\Lambda}|\tilde{x}_T,\lambda_T)}{P(\Lambda | x_0, \lambda_0)}  \frac{\alpha(\tilde{X}|\tilde{\Lambda})}{\alpha(X|\Lambda)} e^{-\Delta \mathcal{S}(X|\Lambda) - \Delta u(X|\Lambda)} 
\label{equation:acceptance-criteria-relationship}
\end{eqnarray}
where $\Delta u(X|\Lambda) \equiv u_T(x_T) - u_0(x_0)$ is the energy difference.
Eq.~\ref{equation:acceptance-criteria-relationship} is the main result of this paper, and is highly general with regard to both the choice of protocol for perturbation and propagation.
%Eq.~\ref{equation:acceptance-criteria-relationship} 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|\Lambda)$ that satisfy Eq.~\ref{equation:acceptance-criteria-relationship} are possible, including the well-known Metropolis-Hastings criterion~\cite{metropolis:jcp:1953:metropolis-monte-carlo,hastings:biometrika:1970:metropolis-hastings},
\begin{eqnarray}
\lefteqn{A(X|\Lambda) =} \label{equation:modified-metropolis-criteria} \\
&amp;&amp;\min \left\{ 1, \frac{\omega_T}{\omega_0} \frac{P(\tilde{\Lambda}|\tilde{x}_T,\lambda_T)}{P(\Lambda | x_0, \lambda_0)}  \frac{\alpha(\tilde{X}|\tilde{\Lambda})}{\alpha(X|\Lambda)} e^{-\Delta \mathcal{S}(X|\Lambda) - \Delta u(X|\Lambda)}  \right\} .  \nonumber
\end{eqnarray}
After generating $(x_T,\lambda_T)$ and evaluating $A(X|\Lambda)$, we generate a uniform random variate $U$.  If $A(X|\Lambda) &lt; U$, then the candidate becomes the next value in the chain, $(x^{(n+1)},\lambda^{(n+1)}) = (x_T,\lambda_T)$.  Otherwise, it  is rejected, we perform a momentum flip, and the next value becomes $(x^{(n+1)},\lambda^{(n+1)})=(\tilde{x}_0,\lambda_0)$.
Alternately, we may perform a momentum flip upon acceptance, $(x^{(n+1)},\lambda^{(n+1)}) = (\tilde{x}_T,\lambda_T)$ and preserve the momentum upon rejection, $(x^{(n+1)},\lambda^{(n+1)})=(x_0,\lambda_0)$.
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 $\pi(x,\lambda)$, but can be mixed with other MCMC moves or with MD~\cite{jun-s-liu:MCMC}.
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.

%Various useful choices for the perturbation and propagation kernels are discussed in the Appendix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PERTURBATION KERNELS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Perturbation Kernels}

A large variety of choices are available for the perturbation kernels $\alpha_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.

\subsection{Stochastically Driven Degrees of Freedom}
Suppose we wish to drive a  torsion angle $\phi$ (an angle subtended by four bonded atoms) \emph{stochastically} by rotating it to a new torsion angle $\phi'$ according to some probability, such as the von Mises circular distribution centered on $\phi$,
\begin{equation}
\eta(\phi' | \phi) = [2 \pi I_0(\kappa)]^{-1} e^{\kappa \cos(\phi' - \phi)} ,
\end{equation}
with $I_0(\kappa)$ denoting the modified Bessel function of order zero and $\kappa &gt; 0$ taking the role of a dimensionless force constant.
Because the stochastic perturbation is made in a non-Cartesian coordinate, a Jacobian $J(\phi)$ must be included to compute $\alpha(x,y)$ in Cartesian coordinates, resulting in the ratio,
\begin{equation}
\frac{\alpha_t(\tilde{y},\tilde{x})}{\alpha_t(x,y)} = \frac{\eta(\phi|\phi') J_{\phi'}(\phi)}{\eta(\phi'|\phi) J_\phi(\phi')} = e^{\kappa \left[ \cos(\phi - \phi') - \cos(\phi' - \phi) \right]}.
\end{equation}
where $J_\phi(\phi') = J_{\phi'}(\phi) =1$ because the transformation (a rotation about a bond vector) preserves the Cartesian phase space volume.

\subsection{Deterministically Driven Degrees of Freedom}

Instead of perturbing the torsion angle stochastically, we can deterministically drive it in small, fixed increments $\Delta \phi$.
In this case, we effectively define an invertible map $\mathcal{M}$ that takes $x \rightarrow y$, such that $y = \mathcal{M} x$ differs from $x$ only in the rotation of the specified torsion $\phi$ by $\Delta \phi$.
To implement this, we may choose a perturbation $\Delta \phi$ from a distribution where $\pm \Delta \phi$ have equal probability, and drive $\phi(x)$ from its current value $\phi_0$ to a final value $\phi_T = \phi_0 + \Delta \phi$ over $T$ steps in equal increments, such that $\phi(x_t)$ is constrained to $\phi_t \equiv (1-t/T) \phi_0 + (t/T) \phi_T$.
In this case, $\alpha_t(x,y) = \delta(y - \mathcal{M} x) J(x)$, where the Jacobian $J(x)$ represents the factor by which Cartesian phase space is compressed on the application of the map $\mathcal{M}$, which is again unity for rotation about a torsion angle by $\Delta \phi$, and, due to the invertibility of the map, the ratio $\alpha_t(\tilde{y}, \tilde{x}) / \alpha_t(x,y) = 1$.

\subsection{Simulation Box Scaling}

%Another possible perturbation is 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)+\Delta V)/V(x)]^{1/3}$ with $\Delta V$ chosen uniformly from $[V-\Delta V_0,V+\Delta V_0]$ applied as a factor of $s^{1/T}$ over the course of $T$ steps.
In this case, the perturbation kernel $\alpha_t(x,y)$ 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 $\alpha(\tilde{X}|\tilde{\Lambda}) / \alpha(X|\Lambda) = s^{3 N}$, where $N$ denotes the number of molecular centers.

\subsection{Thermodynamic Perturbation}

In many driven nonequilibrium processes, there is no direct perturbation to the coordinates, such that $\alpha_t(x,y) = \delta(x-y)$ and the ratio $\alpha(\tilde{X}|\tilde{\Lambda}) / \alpha(X|\Lambda) = 1$.
Instead, only the thermodynamic parameters $\lambda$ are varied in time, carrying the system out of equilibrium through action of the $K_t$ propagation kernels. 
We recover Neal's method~\cite{neal:2004:tech-report:dragging-fast-variables} if the reduced potential $u_t$ is a simple linear interpolation such that $u_t(x) = (1-t/T) u_0(x) + (t/T) u_T(x)$, the probability of choosing protocol $\Lambda$ is symmetric with $\tilde{\Lambda}$, and MC~\cite{metropolis:jcp:1953:metropolis-monte-carlo, hastings:biometrika:1970:metropolis-hastings} is used for the propagation kernel $K_t$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PROPAGATION KERNELS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Propagation Kernels}

The choice of propagation kernels available is also very broad.
If strong driving is performed in selection of $\alpha$, one may elect to choose a time-independent propagation kernel $K_t(x,y) \equiv K(x,y)$ that samples from a stationary distribution $\pi(x)$ of interest.
Alternatively, a strongly time-dependent $K_t$ 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Monte Carlo dynamics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{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 $\alpha_t$) by a method that satisfies detailed balance in $\pi_t$,
\begin{eqnarray}
\pi_t(x) K_t(x,y) &amp;=&amp; \pi_t(\tilde{y}) K_t(\tilde{y},\tilde{x}) ,
\label{eq:detailed-balance}
\end{eqnarray}
where $\pi_t \equiv Z_t^{-1} e^{-u_t(x)}$ for a specified $u_t(x)$.
Many MCMC methods~\cite{jun-s-liu:MCMC}, including Metropolis~\cite{metropolis:jcp:1953:metropolis-monte-carlo, hastings:biometrika:1970:metropolis-hastings} and various hybrid Monte Carlo (HMC) algorithms that combine discrete-timestep integrators with Monte Carlo acceptance/rejection steps~\cite{duane:1987:phys-lett-b:hybrid-monte-carlo,lelievre-stoltz-rousset:2010:GHMC}, obey detailed balance and are commonly used for the simulation of physical systems.

By analogy with Crooks~\cite{crooks:jsp:1998:neq}, we define a \emph{work} $w$ and \emph{heat} $q$ for the nonequilibrium driven process,
\begin{eqnarray}
w(X|\Lambda) &amp; = &amp; \sum_{t=1}^T \left[ u_t(x_t^*) - u_{t-1}(x_{t-1}) \right] \label{eq:work} \\
q(X|\Lambda) &amp; = &amp; \sum_{t=1}^T \left[ u_t(x_t) - u_t(x_t^*) \right] \label{eq:heat}
\end{eqnarray}
such that $w(X|\Lambda)+q(X|\Lambda) = \Delta u(X|\Lambda)$, 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|\Lambda)$,
\begin{eqnarray}
\Delta \mathcal{S}(X | \Lambda)
= -\ln \prod_{t=1}^{T} \frac{\pi_t(x_t^*)}{\pi_t(\tilde{x}_t)}
= -q(X|\Lambda),
\end{eqnarray}
leading to an acceptance probability similar to standard MC, except that the work, $w(X|\Lambda)$, replaces the instantaneous potential energy difference,
\begin{eqnarray}
\frac{A(X|\Lambda)}{A(\tilde{X}|\tilde{\Lambda})} &amp;=&amp; \frac{\omega_T}{\omega_0} \frac{P(\tilde{\Lambda}|\tilde{x}_T,\lambda_T)}{P(\Lambda | x_0, \lambda_0)}  \frac{\alpha(\tilde{X}|\tilde{\Lambda})}{\alpha(X|\Lambda)} e^{-w(X|\Lambda)}.
\label{equation:detailed-balance-acceptance-criteria}
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Deterministic integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\subsection{Deterministic MD}

\subsection{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 $K_t(x,y) = K_t(\tilde{y},\tilde{x})$.  
Hence, $\Delta \mathcal{S}(X|\Lambda) = 0$ and the acceptance ratio is,
\begin{eqnarray}
\frac{A(X|\Lambda)}{A(\tilde{X}|\tilde{\Lambda})} &amp;=&amp; \frac{\omega_T}{\omega_0} \frac{P(\tilde{\Lambda}|\tilde{x}_T,\lambda_T)}{P(\Lambda | x_0, \lambda_0)}  \frac{\alpha(\tilde{X}|\tilde{\Lambda})}{\alpha(X|\Lambda)} e^{-\Delta u(X|\Lambda)}
\label{equation:deterministic-acceptance-criteria},
\end{eqnarray}
where $\Delta u(X | \Lambda) \equiv u_T(x_T) - u_0(x_0)$ 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~\cite{jarzynski:prl:1997:neq-fe, lechner:2006:jcp:fast-switching}.
%by Jarzynski~\cite{jarzynski:prl:1997:neq-fe} and Lechner et al.~\cite{lechner:2006:jcp:fast-switching}.

Symplectic integrators include velocity Verlet~\cite{swope:jcp:1982:velocity-verlet}.
These integrators are also symplectic when utilizing constraints through the application of algorithms such as RATTLE~\cite{andersen:1983:j-comput-phys:RATTLE}, provided that the constraints are iterated to convergence each timestep~\cite{leimkuhler-skeel:1994:j-comput-phys:symplectic-integrators-in-constrained-systems}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\subsection{Stochastic integrators}

\subsection{Stochastic Dynamics}
Stochastic integrators such as velocity Verlet discretizations of Langevin dynamics~\cite{paterlini-ferguson:chem-phys:1998:velocity-verlet-langevin} satisfy detailed balance only with respect to a modified distribution $\hat{\pi}_t$ that differs from $\pi_t$ in a timestep-dependent manner,
\begin{eqnarray}
\hat{\pi}_t(x) K_t(x,y) &amp;=&amp; \hat{\pi}_t(\tilde{x}) K_t(\tilde{y},\tilde{x}) .
\end{eqnarray}
In this case the relative action $\Delta \mathcal{S}(X|\Lambda)$ may not have a simple form, but is still relatively straightforward to compute (see, e.g.~Ref.~\cite{chodera:2011:jcp:dynamical-reweighting}).
As examples, we compute $\Delta \mathcal{S}(X | \Lambda)$ for the overdamped Langevin (Brownian) dynamics integrator of Ermak and Yeh~\cite{ermak-yeh:cpl:1974:brownian-dynamics,ermak:jcp:1975:brownian-dynamics} and the Br\"{u}nger-Brooks-Karplus (BBK) Langevin integrator~\cite{brunger-brooks-karplus:cpl:1984:bbk-integrator,pastor-brooks-szabo:mol-phys:1988:bbk-integrator,schlick} in the Appendix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stochastic integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\color{blue}
%
%\subsection{Nonequilibrium Hybrid Monte Carlo Trials}
%
%The standard hybrid Monte Carlo (HMC) scheme~\cite{duane:1987:phys-lett-b:hybrid-monte-carlo} can be considered in the context of a NCMC trial move.  The standard procedure for generating an HMC trial move is to assign velocities from a Maxwell Boltzmann distribution, propagate for a number of steps, and accept the final move with the change in the Hamiltonian.  Here we derive a more general form of this procedure to include nonequilibrium trajectories.  It can be considered a generalization of \ref{equation:acceptance-criteria-relationship-symplectic}. The derivation of the HMC criterion, however uses a slightly different convention, and is presented here {\color{blue}[JN: or maybe in an appendix, along with the BD derivation?]}
%for completeness.  We consider the HMC trial move as a single perturbation of momentum, followed by a deterministic trajectory (as in \ref{equation:acceptance-criteria-relationship-symplectic}).  To to this, we revisit Eq. \ref{equation:acceptance-criteria-relationship},
%
%\begin{eqnarray}
%\frac{acc(X)_{HMC}}{acc(\Tilde{X})_{HMC}}&amp;=&amp;\frac{\pi_H(x_T)}{\pi_H(x_0)}\prod_{t=1}^{T} \frac{\alpha_t(x_t^*, x_{t-1}) K_t(x_t, x_t^*)}{\alpha_t(x_{t-1}, x_t^*) \, K_t(x_t^*, x_t)}
%\nonumber\\
%&amp;=&amp;\frac{\pi_H(x_T)}{\pi_H(x_0)}\frac{\alpha_1(x_1^*, x_{0})}{\alpha_t(x_{t-1}, x_t^*)}
%\nonumber
%%\\
%%&amp;=&amp;\frac{\pi_H(x_T)}{\pi_H(x_0)}\frac{\alpha^{(p)}_1(-p_1^*, -p_{0})}{\alpha^{(p)}_1(\\p_{0},\\ p_1^*)}
% \label{equation:HMC-criteria-part-A}
%\end{eqnarray}
%
%Since the initial position and thermodynamic parameter remain unchanged, we consider only the probability of selecting the trial momentum $\alpha^{(p)}(p,p')=g(p')$.  We also notice that the reverse selection probability invokes time reversibilty to assign reverse velocities to the selection probability.  Finally, noting that $g(p)=g(-p)$, we obtain
%
%\begin{eqnarray}
%\frac{acc(X)_{HMC}}{acc(\Tilde{X})_{HMC}}&amp;=&amp;\frac{\pi(q_T|\lambda_T)g(p_T)g(p_0)}{\pi(q_0|\lambda_0)g(p_0)g(p_1^*)}\nonumber\\
%acc(X)_{HMC}&amp;=&amp; e^{-h(x_T)+h(x_1^*)}
%\label{equation:HMC-criteria-part-B}
%\end{eqnarray}
%
%This result is the standard HMC criterion~\cite{duane:1987:phys-lett-b:hybrid-monte-carlo}, allowing for the trial trajectory to be a driven nonequilibrium trajectory.  If the $\lambda_0=\lambda_T$, we recover the exact criterion.  Note that the Hamiltonian is not expected to be conserved for the cases where $\lambda$ varies over the course of the trajectory, whereas it is conserved for the fixed $\lambda$ case.  The nonequilibrium version of the HMC criterion is interesting because it is derived using the microscopic reversibility condition.
%
%An analogous, but more complex, acceptance criteria can be derived for stochastic integrators (see Appendix~\ref{appendix:stochastic-integrators}).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OPTIMAL ACCEPTANCE CRITERIA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\section{Optimal acceptance criteria}
%
%{\color{blue} [JDC: Here, we discuss the optimality of modified Metropolis-Hastings, in terms of achieving the maximum possible acceptance rate.  If we can work this out, maybe we should ask Phillip Geissler, who originally gave us this idea, to be a coauthor.]}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ILLUSTRATIVE APPLICATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\section{Illustrative scenarios}
%
%To illustrate the use of nonequilibrium candidate Monte Carlo, we briefly outline a number of scenarios in which it may be applied and the corresponding acceptance criteria that would generate an equilibrium distribution.
%
%\subsection{Charge inversion}
%
%Consider a solvated monoatomic, monovalent ion in water.
%We wish to estimate the free energy difference between the anionic and cationic forms of this ion in an expanded ensemble simulation.
%The allowed values for the thermodynamic parameter $\lambda$ in our state space therefore consist of two states, the cation form and anion forms, each with distinct $u_\lambda(x)$.
%In a standard expanded ensemble Monte Carlo simulation, we would propose Monte Carlo moves that simply invert the charge on the ion, instantaneously converting a cation into an anion, or vice-versa.
%Because the water dipoles arrange themselves in a favorable manner around either species, this dipole configuration would strongly disfavor charge inversion, leading to an extremely low acceptance rate for charge inversion.
%
%To increase the acceptance rate through the use of nonequilibrium candidates, we can propose a protocol $\Lambda$ that inverts the charge through a series of $N$ steps, in which the charge is modified by an increment of magnitude $2/N$ in the same direction each step until the inverted charge is reached.
%The propagation kernels ($K_t$) allow all of the degrees of freedom to evolve by Hamiltonian dynamics under the modified Hamiltonian corresponding to the current charge state, after drawing the initial velocity vector from the Maxwell-Boltzmann distribution.
%In the cation form, we only propose protocols $\Lambda$ that carry us to the anion form, and the perturbation kernel is the identity ($\alpha_t(x,y) = \delta(x - y)$).
%
%The modified Metropolis-Hastings acceptance probability for this type of nonequilibrium candidate move then simplifies to
%\begin{eqnarray}
%A_\mathrm{MH}(X|\Lambda) &amp;=&amp;  \min\left\{1, e^{-w(X | \Lambda)} \right\}
%\end{eqnarray}
%where the work $w(X | \Lambda)$ is simply the work accumulated during the application of charge increments, evaluated at fixed configuration
%\begin{eqnarray}
%w(X | \lambda) &amp;=&amp; \beta \sum_{t=1}^T (q_t - q_{t-1}) \phi_\mathrm{ion}(x_{t-1})
%\end{eqnarray}
%where $q_t$ is the time-dependent charge on the ion (equal to $1 - 2 t/T$ for transformation from a cation to an anion) and $\phi_\mathrm{ion}(x_t)$ is the electric field at the ion center due to the waters, as a function of configuration $x$.
%
%\subsection{Particle insertion}
%
%\subsection{Volume perturbation}
%
%\subsection{Protein sidechain sampling}
%
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NUMERICAL ILLUSTRATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Illustrative Application: Bistable Dimer in a WCA Fluid}
\label{section:wca-application}

\begin{figure}[h]
\noindent
\resizebox{\columnwidth}{!}{\includegraphics{wca-fluid-diagram}}
\caption{{\bf Bistable dimer potential and instantaneous MC moves in WCA fluid.}
An extension move increases the dimer extension by $\Delta r = +r_0$, while a compaction move decreases the dimer extension by $\Delta r = -r_0$. 
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 illustration---the actual simulation described here is of a 3D system.
\label{figure:wca-fluid-diagram}}
\end{figure}

\begin{figure}[h]
\noindent
\resizebox{\columnwidth}{!}{\includegraphics{wca-dimer-trajectory}}
\caption{{\bf Trajectories of WCA dimer system in vacuum and solvent.}
\emph{Left:} The dimer extension $r$ as a function of simulation iteration.
The dotted horizontal line denotes division between compact and extended configurations.
The quantity $\tau$ printed above each plot indicates the estimated integrated autocorrelation time for the dimer extension $r$.
\emph{Right:} Histogram accumulated over trajectory (black), with true equilibrium distribution (red).
Plot titles denote whether simulation was run in vacuum (\emph{vacuum}) or dense WCA fluid (\emph{solvent}), and whether the simulation utilized only 500 GHMC steps per iteration (\emph{MD}) or included instantaneous MC (\emph{MC}) or NCMC moves (\emph{NCMC}) following each iteration.
\label{figure:wca-trajectory}}
\end{figure}

\begin{figure}[h]
\noindent
\resizebox{\columnwidth}{!}{\includegraphics{wca-dimer-efficiency}}
\caption{{\bf Acceptance probabilities and efficiency gain of NCMC proposals over instantaneous MC proposals.}
\emph{Top:} Acceptance probability of NCMC proposals as a function of length of nonequilibrium proposal trajectory (black dots), compared with instantaneous MC proposal (red line).
\emph{Bottom:} Relative acceptance efficiency of NCMC proposals compared to instantaneous MC proposal (black dots), with no efficiency improvement (gain factor of unity) denoted as red line.
Estimated 95\% confidence intervals are shown as vertical lines.
%{\color{red}[JDC: Note that error bars are probably large on this log-plot.  Change this to linear scale instead of log scale?  Or should we perform umbrella sampling calculations to compute these log-speedup factors very precisely?]}
%DM: Can you show them?
\label{figure:wca-efficiency}}
\end{figure}

To illustrate the use of NCMC for enhancing sampling efficiency in a condensed-phase system, we ran simulations with a bistable dimer (adapted from Section 1.3.2.4 of Ref.~\cite{lelievre-stoltz-rousset:2010:free-energy-computations}) 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 = r_0$ (the \emph{compact} configuration, where the particles are in contact) and $r = 2 r_0$ (the \emph{extended} configuration), and an interatomic potential barrier height of $5 \, k_\mathrm{B} T$.
In the solvated simulations, the dimer was immersed in a dense bath of particles that interact with the bonded particles and each other via the Weeks-Chandler-Andersen (WCA) soft repulsive potential~\cite{weeks-chandler-andersen:1971:jcp:wca-fluid} under dense fluid conditions (reduced density $\rho \sigma^3 = 0.96$).
Each simulation ``iteration'' consisted of velocity reassignment from the Maxwell-Boltzmann distribution, 500 steps of generalized hybrid Monte Carlo (GHMC) dynamics~\cite{gustafson:statistics-and-computing:1998:ghmc,jun-s-liu:MCMC,lelievre-stoltz-rousset:2010:free-energy-computations,lelievre-stoltz-rousset:2010:GHMC} (essentially, a Metropolis-corrected form of Langevin dynamics, henceforth referred to here as MD), optionally followed by either an instantaneous MC move or NCMC move prior to recording the dimer extension distance $r$.

In vacuum, the MD simulation of the dimer shows slow hopping between compact and extended configurations (Figure~\ref{figure:wca-trajectory}, top trace), resulting in slow convergence of estimates of quantities that are sensitive to transitions across the barrier, such as the ratio of populations of compact and extended states.
The average number of simulation iterations required to sample a new uncorrelated equilibrium configuration can be quantified by computing the integrated autocorrelation time $\tau$ for the dimer extension $r$  (Figure~\ref{figure:wca-trajectory}, subplot titles), which here is $\approx$ 59.2 iterations.
The convergence rate in \emph{vacuum} is greatly increased by the introduction of periodic instantaneous MC attempts that modify the interparticle distance $r$ by $\pm r_0$ (Figure~\ref{figure:wca-trajectory}, second trace from top).
When the dimer is immersed in a dense fluid of WCA particles, however, MD alone results in extremely slow barrier crossings (Figure~\ref{figure:wca-trajectory}, middle trace), and hence extremely poor convergence rates.
Unlike the vacuum case, introduction of instantaneous MC moves does not significantly reduce the correlation time $\tau$ for barrier crossing (Figure~\ref{figure:wca-trajectory}, second trace from bottom).
However, performing these same dimer expansion and contraction moves over 2048-step NCMC moves (Figure~\ref{figure:wca-trajectory}, bottom trace) once again allows the system to rapidly mix between both compact and extended states, requiring only $\approx 5 \times$ more evaluations of the force per iteration than the standard MD simulation for an effective $\approx 75 \times$ reduction in the correlation time for dimer barrier crossing, resulting a real efficiency gain of more than an order of magnitude.

Even more remarkable is the degree to which the NCMC procedure enhances the \emph{acceptance probability} of the proposed Monte Carlo move in dense solvent.
Here, the probability that an instantaneous MC proposal of $\pm 2 w$ is accepted is estimated to be $\approx 10^{-18}$.
Through the use of NCMC proposals, the acceptance probability can be increased to $\approx 10^{-1}$, achieving an acceptance rate of 12.6\% using nonequilibrium proposal trajectories of 2048 steps (Figure~\ref{figure:wca-efficiency}, top).
We can quantify the \emph{acceptance efficiency gain} relative to instantaneous moves as the ratio of acceptance probabilities of NCMC to instantaneous MC over the relative cost (in number of whole-system force or energy calculations), $(1/T) (\expect{A}_T / \expect{A}_0)$.  Here, $\expect{A}_n$ denotes the equilibrium expectation of the acceptance probability for $n$ steps of NCMC switching, with $n=0$ denoting instantaneous MC. 
%{\color{blue}(This accounting does not consider optimizations that can be made with pairwise-additive force fields, though the impact of these optimizations is at most a factor of $N$.)}
%A gain of unity would mean that the acceptance probability was increased only by a factor of the number of the relative cost, resulting in no net improvement.
NCMC moves provide a vast \emph{gain} in acceptance efficiency over instantaneous MC, reaching an efficiency boost of $\approx 10^{14}$ in this particular example (Figure~\ref{figure:wca-efficiency}, bottom).
Though the computed statistical error in the estimated instantaneous acceptance probability is likely underestimated due to the extremity of this quantity (see Methods), it is clear that the efficiency gain is no doubt remarkable.

%{\color{red}[JDC: We should temper this amazing efficiency gain by noting that the efficiency gain over MC is determined almost exclusively by the properties of the solvent, while the real efficiency gain is determined almost entirely by the barrier height.  By manipulating the solvent density or the barrier height, we can mostly independently modulate either the MC speedup or real speedup.  In an extreme case of the vacuum system, NCMC offers no MC speedup, since there is no solvent.}

As seen in Fig.~\ref{figure:wca-efficiency}, there is a tradeoff between the switching time and the gain in efficiency of NCMC over the equivalent instantaneous Monte Carlo step, which can be considered a special case of NCMC where $T=1$ steps are taken and $K_1(x,y) = \delta(x-y)$.
For very short switching times, as in the illustrative example, there may actually be a \emph{loss} in efficiency, followed by superlinear \emph{gain} in efficiency, until a plateau is reached where longer switching times result in no further improvement in acceptance probability; switching is nearly reversible.
(A similar plateau is observed in the modified parallel tempering protocol of Ref.~\cite{jarzynski:pnas:2009:nonequilibrium-swaps}.)
In the intermediate switching time regime, it should be possible to achieve useful acceptance probabilities without paying the full price of quasi-reversible switching, but precise details of these regimes is expected to be problem-dependent.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISCUSSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Epilogue}

%The described procedure provides a simple but effective way to borrow Monte Carlo moves used in gas-phase or implicit solvent simulations and apply them to simulations in explicit condensed phases.



We note that switching trajectories contain potentially useful information.
Indeed, several methods~\cite{minh-adib:2008:prl:minh-adib,minh-chodera:2009:jcp:minh-chodera,minh-chodera:2011:jcp:mis-estimator} have recently been developed to estimate \emph{equilibrium} properties from nonequilibrium samples through the application of statistical estimator theory to nonequilibrium fluctuation theorems~\cite{jarzynski:prl:1997:neq-fe,crooks:pre:1999:cft,crooks:pre:2000:pathensemble}; these are particularly relevant to switching between multiple thermodynamic states.
Though the development of efficient estimators that utilize \emph{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~\cite{frenkel:2004:pnas:rejected-states}.
Recent work has also suggested that estimation of the thermodynamic metric tensor along the nonequilibrium parameter switching path~\cite{salamon-berry:prl:1983:thermodynamic-length,crooks:2007:prl:thermodynamic-length,feng-crooks:2009:pre:thermodynamic-length,minh-chodera:2011:jcp:mis-estimator}, could prove useful in adaptively optimizing the switching protocol~\cite{shenfeld-xu-eastwood-dror-shaw:2009:pre:optimizing-thermodynamic-length}.

%{\color{red}Additional points we could cover in discussion:
%\begin{itemize}
%  \item Discuss relation with Jerome's POSH algorithm, and why it works so well despite not being exact?
%  \item Discuss tradeoff between switching time and efficiency.
%  \item Allude to Jerome's idea of \emph{also} using the nonequilibrium portions of the trajectories through some form of reweighting.
%  \item Relationship to Cool-Walking and Opps-Schofield?
%  \item Discuss how we expect these to be useful for biomolecules?  Give ideas of switching protocols we'd like to see?
%  \item Optimization with thermodynamic length.
%  \item How close do we have to be to the ``entropically reversible limit'' for switching?  Could this be prohibitively slow for some systems?
%\end{itemize}
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATERIALS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{materials}

\begin{figure}[tbp]
\noindent
\resizebox{\columnwidth}{!}{\includegraphics{wca-dimer-umbrella}}
\caption{{\bf Umbrella sampling simulation of the dimer in WCA solvent.}
\emph{Left:} The dimer extension $r$ as a function of simulation iteration, with the dotted line denoting the division between compact and extended configurations.
\emph{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.
\label{figure:wca-dimer-umbrella}}
\end{figure}

\section{WCA Dimer Simulations}
The dimer system considered here consists of two particles that interact via a double-well ``bonded'' potential in the interatomic distance~$r$,
\begin{eqnarray}
U_\mathrm{bond}(r) &amp;=&amp; h \left[ 1 - \frac{(r - r_0 - s)^2}{s^2}\right]^2
\end{eqnarray}
with $h = 5 \, k_\mathrm{B} T$, $r_0 = r_\mathrm{WCA}$, and $s = r_\mathrm{WCA}/2$, where $r_\mathrm{WCA} \equiv 2^{1/6} \sigma$.
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~\cite{weeks-chandler-andersen:1971:jcp:wca-fluid},
\begin{eqnarray}
U_\mathrm{WCA}(r) &amp;=&amp; \begin{cases}
4 \epsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right] + \epsilon , &amp; r &lt; r_\mathrm{WCA} \\
0 &amp; r \ge r_\mathrm{WCA}
\end{cases},
\end{eqnarray}
with mass $m = 39.9$ amu, $\sigma = 3.4$ \AA, and $\epsilon = 120 \, k_\mathrm{B} T$.
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 $\rho \sigma^3 = 0.96$.
For all simulations, the reduced temperature is $k_\mathrm{B} T / \epsilon = 0.824$.
A custom Python code making use of the GPU-accelerated OpenMM package~\cite{friedrichs:2009:j-comput-chem:openmm,eastman:2010:comp-sci-eng:openmm,eastman:2010:j-comput-chem:openmm} and the PyOpenMM Python wrapper~\cite{pyopenmm} was used to conduct the simulations.
All scripts are available for download from \url{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)~(\cite{gustafson:statistics-and-computing:1998:ghmc,jun-s-liu:MCMC,lelievre-stoltz-rousset:2010:free-energy-computations,lelievre-stoltz-rousset:2010:GHMC}) for all our simulations.
GHMC is based on a velocity Verlet discretization~\cite{swope:jcp:1982:velocity-verlet} 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 \tau$, where $\tau = \sqrt{\sigma^2 m / \epsilon}$, and the collision rate was set to $\tau^{-1}$.
With this timestep, the acceptance probability is $99.929\pm0.001$\%; the resulting dynamics closely approximates Langevin dynamics.

In simulations employing instantaneous Monte Carlo moves, a perturbation $\Delta r$ to the interatomic distance $r$ of the dimer was chosen according to,
\begin{eqnarray}
\Delta r &amp;=&amp; \begin{cases}
+2 w &amp; \mathrm{if} \:\: r &lt; 1.5 r_0 \\
-2 w &amp; \mathrm{if} \:\: 1.5 r_0 \le r \le 3 r_0 \\
0 &amp; \mathrm{otherwise}
\end{cases} .
\end{eqnarray}
The dimer was contracted or expanded about the bond midpoint to generate a new configuration $x_\mathrm{new}$ with dimer extension $r_\mathrm{new}$ from the old configuration $x_\mathrm{old}$ with dimer extension $r_\mathrm{old}$, and the move accepted or rejected with the Metropolis-Hastings criterion~\cite{hastings:biometrika:1970:metropolis-hastings},
\begin{eqnarray}
A(x_\mathrm{new} | x_\mathrm{old}) &amp;=&amp; \mathrm{min}\left\{1, e^{-\beta[U(x_\mathrm{new})+U(x_\mathrm{old})]} J_r(x_\mathrm{old},x_\mathrm{new}) \right\} \nonumber \\
\end{eqnarray}
where the Jacobian ratio term $J_r(x_\mathrm{old}, x_\mathrm{new}) = (r_\mathrm{new} / r_\mathrm{old})^2$ accounts for the expansion and contraction of phase space due to the Monte Carlo proposals.

For simulations employing NCMC moves, proposals were made by selecting a new velocity vector from the Maxwell-Boltzmann distribution, integrating $n$ steps of velocity Verlet dynamics~\cite{swope:jcp:1982:velocity-verlet} for all bath atoms as the dimer extension was driven from $r_o$ to $r_n$ in equal steps of size $\Delta r / n$, and accepting or rejecting based on the modified Metropolis criteria for symplectic integrators (Eqs.~\ref{equation:modified-metropolis-criteria} and \ref{equation:deterministic-acceptance-criteria}),
\begin{eqnarray}
A(X) &amp;=&amp; \mathrm{min}\left\{1, e^{-\beta[H(x_t)+H(x_0)]} J_r(x_o,x_n) \right\} .
\end{eqnarray}
The Jacobian ratio is also $(r_n/r_o)^2$.
The final configuration after the MC or NCMC rejection procedure was stored and plotted to generate Figure~\ref{figure:wca-trajectory}.

The mean acceptance probabilities for each switching time $\tau$ can be estimated via the sample mean
\begin{eqnarray}
\expect{A}_\tau &amp;\approx&amp; \frac{1}{N} \sum_{n=1}^N A(X_n)  .
\end{eqnarray}
For numerical stability, logarithms of $A(X_n)$ were stored, as $a_n \equiv \ln A(X_n)$. 
We then estimated $\ln \expect{A}_\tau$ (shown in Figure~\ref{figure:wca-efficiency}) as
\begin{eqnarray}
\ln \expect{A}_\tau &amp;\approx&amp; - \ln N + \ln b + \sum_{n=1}^N e^{a_n - b} 
\end{eqnarray}
where $b \equiv \max_n a_n$.

Integrated autocorrelation times were estimated using the rapid scheme described in Section 5.2 of Ref.~\cite{chodera:jctc:2007:ptwham}.

%The statistical error in the estimate of $\ln \expect{A}_\tau$ can be estimated in the standard way,
%\begin{eqnarray}
%\delta^2 \ln \expect{A}_\tau &amp;\approx&amp; \frac{\delta^2 \expect{A}_\tau}{\expect{A}_\tau^2} \\
%\delta^2 \expect{A}_\tau &amp;\approx&amp; \frac{\mathrm{var}A(X_n)}{N} \nonumber 
%\end{eqnarray}
%We note that this is likely to be an \emph{underestimate} of the true statistical uncertainty in the estimated $\ln \expect{A}_\tau$ because of the deficiencies of first-order Taylor expansion propagation of error and the heavy asymmetric tails in the distribution of $A(X_n)$.

The acceptance probabilities plotted in Figure~\ref{figure:wca-efficiency} 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 $\ln \expect{A}_\tau$ 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.

The reference distribution for the interparticle distribution $\mathcal{P}(r)$ plotted in red on the right side of Figure~\ref{figure:wca-trajectory} was computed analytically for the vacuum system,
\begin{eqnarray}
\mathcal{P}_\mathrm{vac}(r) &amp;\propto&amp; 4 \pi r^2 e^{-\beta U_\mathrm{bond}(r)} .
\end{eqnarray}
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, 
\begin{eqnarray}
U_\mathrm{umbrella}(r) &amp;=&amp; k_\mathrm{B} T \ln r^2 + \theta(r_\mathrm{min} - r) (K/2) [r-r_\mathrm{min}]^2 \nonumber \\
&amp;&amp; \mbox{} + \theta(r - r_\mathrm{max}) (K/2) [r-r_\mathrm{max}]^2 ,
\end{eqnarray}
where $r_\mathrm{min} = r_0$, $r_\mathrm{max} = 2.05 r_0$, and $K = k_\mathrm{B} T / \eta^2$, with $\eta = 0.3$ \AA, and $\theta(r)$ is the Heaviside function that assumes a value of unity for $r \ge 0$ and zero otherwise.
The true solvated interparticle distribution $p(r)$ was estimated by reweighting the data produced from this simulation, using the relationship
\begin{equation}
\mathcal{P}_\mathrm{sol}(r) \propto \frac{\sum_{n=1}^N \delta(r - r_n) \, e^{-\beta [U_\mathrm{bond}(r_n) - U_\mathrm{umbrella}(r_n)]}}{\sum_{n=1}^N e^{-\beta [U_\mathrm{bond}(r_n) - U_\mathrm{umbrella}(r_n)]}}
\end{equation}
where $r_n$ denotes the bond separation for sample $n$, and a finite-width histogram bin was used instead of the delta function $\delta(r)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACKNOWLEDGMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{acknowledgments}
The authors thank Gabriel Stoltz (CERMICS, \'{E}cole des Ponts ParisTech); Jed W.~Pitera (IBM Almaden Research Center); Yael Elmatad, Anna Schneider, Paul Nerenberg, Todd Gingrich, David Chandler, David Sivak, Phillip Geissler, Michael Gr\"{u}nwald, and Ulf R{\"o}rb{\ae}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.
JDC was supported through a distinguished postdoctoral fellowship from the California Institute for Quantitative Biosciences (QB3) at the University of California, Berkeley.
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.
Additionally, the authors are grateful to {\sc OpenMM} developers Peter Eastman, Mark Friedrichs, Randy Radmer, and Christopher Bruns for their generous help with the {\sc OpenMM} GPU-accelerated computing platform and associated {\sc PyOpenMM} Python wrappers.
The submitted manuscript has been created in part 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIBLIOGRAPHY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliography{neq-mc}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\eject
\appendix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pathwise detailed balance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{widetext}

\section{Proof that NCMC preserves the equilibrium distribution}
\label{appendix:NCMC-proof}

Following the proof for GHMC in Ref.~\cite{lelievre-stoltz-rousset:2010:free-energy-computations}, here we show that NCMC preserves the equilibrium distribution.  The expected acceptance rate for NCMC moves initiated from $(x,\lambda)$ is,
\begin{eqnarray}
\alpha(x,\lambda) &amp;\equiv&amp; \int d\Lambda \int dX \, \rho(X, \Lambda | x, \lambda) \, A(X|\Lambda) .
\end{eqnarray}

Suppose that we have a variate $(x^{(n)},\lambda^{(n)})$ drawn from the equilibrium distribution $\pi(x,\lambda)$.
The probability density of the next value in the chain, $p(x^{(n+1)},\lambda^{(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 $(x^{(n+1)},\lambda^{(n+1)}) = (\tilde{x}^{(n)},\lambda^{(n)})$ is,
\begin{eqnarray}
\int dx \sum_\lambda \pi(x,\lambda) [1 - \alpha(x,\lambda)] \delta(\tilde{x}-x^{(n+1)}) \delta(\lambda - \lambda^{(n+1)}) &amp; = &amp; \pi(\tilde{x}^{(n+1)}, \lambda^{(n+1)}) [1 - \alpha(\tilde{x}^{(n+1)},\lambda^{(n+1)})] .
\label{eq:rejection-term}
\end{eqnarray}
The latter contribution from accepting the candidate such that $(x^{(n+1)},\lambda^{(n+1)}) = (x_T,\lambda_T)$ is,
\begin{eqnarray}
&amp; &amp; \int dx \sum_{\lambda}  \pi(x,\lambda) \int dX \int d\Lambda \, \rho(X,\Lambda | x, \lambda) \, A(X|\Lambda) \delta(x_T-x^{(n+1)}) \delta(\lambda_T - \lambda^{(n+1)}) \nonumber \\
&amp; = &amp; \int dx_0 \sum_{\lambda_0} \int dX \int d\Lambda \, \left[ \pi(x_0,\lambda_0) \, \rho(X,\Lambda | x_0, \lambda_0) \, A(X|\Lambda)  \right] \delta(x_T-x^{(n+1)}) \delta(\lambda_T - \lambda^{(n+1)}) \nonumber \\
&amp; = &amp; \int dx_T \sum_{\lambda_T} \int d\tilde{X} \int d\tilde{\Lambda} \, \left[ \pi(\tilde{x}_T,\lambda_T) \, \rho(\tilde{X},\tilde{\Lambda} | \tilde{x}_T, \lambda_T) \, A(\tilde{X}|\tilde{\Lambda})  \right] \delta(x_T-x^{(n+1)}) \delta(\lambda_T - \lambda^{(n+1)}) \nonumber \\
&amp; = &amp; \pi(\tilde{x}^{(n+1)},\lambda^{(n+1)}) \alpha(\tilde{x}^{(n+1)}, \lambda^{(n+1)}),
\label{eq:acceptance-term}
\end{eqnarray}
where the pathwise detailed balance condition (Eq. \ref{equation:pathwise-detailed-balance}) is used.

Taking the sum of Eqs. \ref{eq:rejection-term} and \ref{eq:acceptance-term}, we find that the equilibrium distribution is preserved,
\begin{eqnarray}
p(x^{(n+1)},\lambda^{(n+1)}) = \pi(x^{(n+1)},\lambda^{(n+1)})
\end{eqnarray}

By analogous reasoning, maintaining the momentum upon rejection, $(x^{(n+1)},\lambda^{(n+1)} = (x^{(n)},\lambda^{(n)})$, and flipping it upon acceptance, $(x^{(n+1)},\lambda^{(n+1)} = (\tilde{x}_T,\lambda_T)$ will also preserve the equilibrium distribution.

%\end{widetext}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Acceptance criteria for Brownian dynamics integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{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~\cite{ermak-yeh:cpl:1974:brownian-dynamics,ermak:jcp:1975:brownian-dynamics}
\begin{eqnarray}
%x_t &amp;=&amp; x_t^* - \eta \, \Delta t \, \nabla_x H(x_t) + \sqrt{2 \, \eta \, \Delta t} \,\xi_t
x_t &amp;=&amp; x_t^* - \frac{\Delta t}{\gamma m} F_t(x_t) + \sqrt{2} \left( \frac{\Delta t}{\gamma m}\right)^{1/2} \,\xi_t
\end{eqnarray}
where $m$ is the particle mass, $F_t(x) \equiv -(\partial/\partial x) H_t(x)$ is the (potentially time-dependent) systematic force, and $\gamma$ is an effective collision frequency or friction coefficient with units of inverse time.
The noise history $\xi_t$ for each degree of freedom is a normal random variate with zero mean and variance $\beta^{-1}$, drawn from the distribution
\begin{eqnarray}
\phi(\xi_t) &amp;=&amp; \frac{1}{\sqrt{2 \pi \beta^{-1}}} \exp\left[-\frac{\beta}{2} \xi_t^2\right]
\end{eqnarray}

In NCMC, every application of the propagation kernel $K_t$ produces a transition $x_t^* \rightarrow x_t$ determined by the noise history variable $\xi_t$, there is a corresponding $\tilde{\xi}_t$ that generates the opposite step, $x_t \rightarrow x_t^*$.
By noting
\begin{eqnarray}
x_t  &amp;=&amp; x_t^* - \frac{\Delta t}{\gamma m} F_t(x_t) + \sqrt{2} \left( \frac{\Delta t}{\gamma m}\right)^{1/2} \,\xi_t \nonumber \\
x_t^* &amp;=&amp; x_t - \frac{\Delta t}{\gamma m} F_t(x_t^*) + \sqrt{2} \left( \frac{\Delta t}{\gamma m}\right)^{1/2} \,\tilde{\xi}_t
\end{eqnarray}
we can compute the relationship,
\begin{eqnarray}
\tilde{\xi}_t &amp;=&amp;  \frac{1}{\sqrt{2}} \left(\frac{\Delta t}{\gamma m}\right)^{1/2} \left[F_t(x_t) + F_t(x_t^*)\right] - \xi_t,
\end{eqnarray}

Then, the ratio of transition kernels appearing in Eq.~\ref{equation:ratio-of-propagation-kernels} can be written in terms of noise history $\xi_t$ and the computed reverse noise history $\xi_t^*$,
%\begin{eqnarray}
%\Delta \mathcal{S}(X) &amp;=&amp; - \ln \prod_{t=1}^T \frac{K_t(x_t, x_t^*)}{K_t(x_t^*, x_t)} \nonumber \\
%&amp;=&amp; - \ln \prod_{t=1}^T \frac{\phi(\xi_t^*) \left|\frac{\partial x_t^*}{\partial \xi_t^*}\right|}{\phi(\xi_t) \left|\frac{\partial x_t}{\partial \xi_t}\right|} \nonumber \\
%&amp;=&amp; - \ln \prod_{t=1}^T \frac{\phi(\xi_t^*)}{\phi(\xi_t)} \nonumber \\
%&amp;=&amp; - \ln \prod_{t=1}^T \exp\left[ -\frac{\beta}{2} ({\xi_t^*}^2 - {\xi_t}^2)\right] \nonumber \\
%&amp;=&amp; \frac{\beta}{2} \sum_{t=1}^T ({\xi_t^*}^2 - {\xi_t}^2)
%\end{eqnarray}
\begin{eqnarray}
\Delta \mathcal{S}(X) = - \ln \prod_{t=1}^T \frac{K_t(x_t, x_t^*)}{K_t(x_t^*, x_t)}  = - \ln \prod_{t=1}^T \frac{\phi(\tilde{\xi}_t) \left|\frac{\partial x_t^*}{\partial \tilde{\xi}_t}\right|}{\phi(\xi_t) \left|\frac{\partial x_t}{\partial \xi_t}\right|} = - \ln \prod_{t=1}^T \frac{\phi(\tilde{\xi}_t)}{\phi(\xi_t)} = - \ln \prod_{t=1}^T \exp\left[ -\frac{\beta}{2} ({\tilde{\xi}_t}^2 - {\xi_t}^2)\right] = \frac{\beta}{2} \sum_{t=1}^T ({\tilde{\xi}_t}^2 - {\xi_t}^2)
\label{equation:Ermak-action}
\end{eqnarray}
where the tildes are dropped because the microstate $x$ contains no momenta.
The quantity $\left|\partial x_t/\partial \xi_t\right|$ represents the Jacobian for the change of variables from the $\xi_t$ to $x_t$, and the Jacobians in the numerator and denominator cancel.
The quantity in Eq. \ref{equation:Ermak-action} can easily be computed during integration.
%For integrators of Langevin dynamics, such as the Leapfrog Langevin integrator~\cite{paterlini-ferguson:chem-phys:1998:velocity-verlet-langevin}, the computation is similar.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Acceptance criteria for Langevin integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Acceptance criteria for Langevin integrator of Brooks, Br\"{u}nger, and Karplus (BBK)}

The Br\"{u}nger-Brooks-Karplus (BBK) stochastic integrator~\cite{brunger-brooks-karplus:cpl:1984:bbk-integrator,pastor-brooks-szabo:mol-phys:1988:bbk-integrator} is a popular integrator for simulating Langevin dynamics.
The application of the corresponding propagator $K_t(x_t^*,x_t)$ can be implemented via a velocity Verlet scheme~\cite{schlick,lelievre-stoltz-rousset:2010:free-energy-computations},
\begin{eqnarray}
v'_{t} &amp;=&amp; v^*_t + \frac{\Delta t}{2 m} \left(F_t(r^*_t) - \gamma m v^*_t + \sqrt{\frac{2 \gamma m}{\Delta t}} \xi_t\right) \nonumber \\
r_{t} &amp;=&amp; r^*_t + \Delta t \, v'_t \nonumber \\
v_{t} &amp;=&amp; \frac{1}{1 + \frac{\gamma \Delta t}{2}} \left[ v'_{t} + \frac{\Delta t}{2 m} \left( F_t(r_t) + \sqrt{\frac{2 \gamma m}{\Delta t}} \xi'_{t} \right) \right]
\end{eqnarray}
where $r_t$ and $v_t$ denote the respective Cartesian position and velocity components of the microstate $x_t$,  $\gamma$ the effective collision frequency with units of inverse time, and $m$ the particle mass.  
$v'_t$ is an auxiliary variable used only in simplifying the mathematical representation of the integration scheme.
Note that we require \emph{two} random variates, $\xi_t$ and $\xi'_t$, 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 $\tilde{X}$ (see, e.g., Section 2.2.3.2 of \cite{lelievre-stoltz-rousset:2010:free-energy-computations}).

The noise history terms $\xi_t$ and $\xi'_t$ are normal random variates with zero mean and variance $\beta^{-1}$.
Their joint distribution can therefore be written,
\begin{eqnarray}
\psi(\xi_t,\xi'_t) &amp;=&amp; \frac{1}{2 \pi \beta^{-1}} \exp\left[-\frac{\beta}{2} \left(\xi_t^2 + {\xi'_t}^2\right) \right] .
\end{eqnarray}
%\begin{eqnarray}
%\psi(\xi_t,\xi_{t+1/2}) &amp;=&amp; \frac{1}{2 \pi \beta^{-1}} \exp\left[-\frac{\beta}{2} \left(\xi_t^2 + \xi_{t+1/2}\right)^2 \right] .
%\end{eqnarray}

For every step $x_t^* \rightarrow x_t$, the positions and velocities undergo a transition $(r_t^*,v_t^*) \rightarrow (r_t,v_t)$ determined by the noise variables $(\xi_t, \xi'_t)$.
A corresponding choice of noise variables {\color{blue}$(\tilde{\xi}_t,\tilde{\xi}'_t)$} will generate the reverse step, $\tilde{x}_t \rightarrow \tilde{x}_t^*$, carrying $(r_t,-v_t) \rightarrow (r_t^*,-v_t^*)$.
With a little algebra, it is seen that,
\begin{eqnarray}
\tilde{\xi}_t &amp;=&amp; \xi'_t - \sqrt{2 \gamma \, m \, \Delta t} \, v_t \nonumber \\
\tilde{\xi}'_t &amp;=&amp; \xi_t - \sqrt{2 \gamma \, m \, \Delta t} \, v_t^* .
\end{eqnarray}
In order to write the ratio of transition kernels appearing in Eq.~\ref{equation:ratio-of-propagation-kernels} in terms of noise variables $(\xi_t,\xi'_t)$ and the computed reverse noise variables $(\tilde{\xi}_t, \tilde{\xi}'_t)$, we must first compute the Jacobian $J(\xi_t, \xi'_t)$ because the random variates are not in Cartesian space,
\begin{eqnarray}
J(\xi_t, \xi'_t) &amp;\equiv&amp; \det \left[ 
\begin{array}{cc}
\frac{\partial r_t}{\partial \xi_t} &amp; \frac{\partial v_t}{\partial \xi_t} \\
\frac{\partial r_t}{\partial \xi'_t} &amp; \frac{\partial v_t}{\partial \xi'_t}
\end{array}
\right] ,
\end{eqnarray}
which can be shown to be independent of $\xi_t$ and $\xi'_t$.
The conditional path action difference can now be computed,
\begin{eqnarray}
\Delta \mathcal{S}(X) = - \ln \prod_{t=1}^T \frac{K_t(\tilde{x}_t, \tilde{x}_t^*)}{K_t(x_t^*, x_t)}  = - \ln \prod_{t=1}^T \frac{\psi(\tilde{\xi}_t,\tilde{\xi}'_t) J(\tilde{\xi}_t,\tilde{\xi}'_t)}{\psi(\xi_t,\xi'_t) J(\xi_t,\xi'_t)}  = \frac{\beta}{2} \sum_{t=1}^T \left[ \left(\left.{\tilde{\xi}_t}\right.^2 + \left.{\tilde{\xi}'_t}\right.^2\right) - \left({\xi_t}^2 + {\xi'_t}^2\right)\right]
\end{eqnarray}
where the ratio of Jacobians $J(\tilde{\xi}_t, \tilde{\xi}'_t) / J(\xi_t, \xi'_t)$ cancels because the Jacobians are independent of the noise variates.

\end{widetext}

\end{document}</pre></body></html>