% PRL look and style (easy on the eyes)
\documentclass[aps,pre,twocolumn,showkeys,nofootinbib,superscriptaddress]{revtex4-1}
% Two-column style (for submission/review/editing)
%\documentclass[aps,prl,preprint,nofootinbib,superscriptaddress]{revtex4-1}

\usepackage{palatino}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
%\usepackage[mathbf,mathcal]{euler}
%\usepackage{citesort}
\usepackage{dcolumn}
\usepackage{boxedminipage}
\usepackage{verbatim}
\usepackage[colorlinks=true,citecolor=blue,linkcolor=blue]{hyperref}

%\graphicspath{{figures/}}

%\bibliographystyle{apsrevlong}

% 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}                % <.> 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{\z}{\bfv{z}}
%\newcommand{\p}{\bfv{p}}
%\newcommand{\q}{\bfv{q}}
%\newcommand{\T}{\mathrm{T}}

\newcommand{\op}[1]{\mathcal{#1}}
%\newcommand{\vektor}[1]{{\bf #1}}
%\renewcommand{\vec}[1]{\bfm{#1}} % italicized boldface
%\renewcommand{\vec}[1]{{\bf #1}}
\renewcommand{\vec}[1]{{\mbox{\boldmath{$#1$}}}}
\newcommand{\q}{\vec{q}} % coordinates
\newcommand{\p}{\vec{p}} % momenta
\newcommand{\z}{\vec{z}} % phase space point
\newcommand{\x}{\vec{x}} % general vector x
\newcommand{\M}{\bfm{M}} % diagonal mass matrix
\newcommand{\grad}{\nabla}
\newcommand{\timeavg}[1]{\overline{#1}}                    % time average over a trajectory

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

\newcommand{\bTheta}{\bfm{\Theta}}
\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}\left[ #2 \,|\, #1\right]}
\newcommand{\T}{\mathrm{T}}                                % T used in matrix transpose
\newcommand{\qed}{\hfill \mbox{\raggedright \rule{.07in}{.1in}}}  % QED mark
\newcommand{\var}[1]{\mathrm{var} \, #1}
\newcommand{\cov}[2]{\mathrm{cov}\left( {#1}, {#2} \right)}

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

%% TITLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Dynamical reweighting:\\
Improved estimates of dynamical properties from simulations at multiple temperatures}

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

\author{William C. Swope}
 \email{swope@us.ibm.com}
 \affiliation{IBM Almaden Research Center, San Jose, CA 95120}

\author{Frank No\'{e}}
 \email{frank.noe@fu-berlin.de}
 \affiliation{DFG Research Center Matheon, FU Berlin, Arnimallee 6, 14195 Berlin, Germany}

\author{Jan-Hendrik Prinz}
 \email{jan.prinz@fu-berlin.de}
 \affiliation{IWR, University of Heidelberg, INF 368, 69120 Heidelberg, Germany}
 \affiliation{DFG Research Center Matheon, FU Berlin, Arnimallee 6, 14195 Berlin, Germany}

\author{Michael R. Shirts}
 \email{michael.shirts@virginia.edu}
 \affiliation{Department of Chemical Engineering, University of Virginia, VA 22904}

\author{Vijay S. Pande}
 \email{pande@stanford.edu}
 \affiliation{Department of Chemistry, Stanford University, Stanford, CA 94305}

\date{\today}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT/pacs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}

Dynamical averages based on \emph{functionals} of dynamical trajectories, such as time-correlation functions, play an important role in determining kinetic or transport properties of matter.
At temperatures of interest, the expectations of these quantities are often dominated by contributions from rare events, making the precise calculation of these quantities by molecular dynamics simulation difficult.
Here, we present a reweighting method for combining simulations from multiple temperatures (or from simulated or parallel tempering simulations) to compute an optimal estimate of the dynamical properties at the temperature of interest without the need to invoke an approximate kinetic model (such as the Arrhenius law).
Continuous and differentiable estimates of these expectations at any temperature in the sampled range can also be computed, along with an assessment of the associated statistical uncertainty.
For rare events, aggregating data from multiple temperatures can produce an estimate of the desired precision at greatly reduced computational cost compared with simulations conducted at a single temperature.
Here, we describe use of the method for the canonical (NVT) ensemble using four common models of dynamics (canonical distribution of Hamiltonian trajectories, Andersen thermostatting, Langevin, and overdamped Langevin or Brownian dynamics), but it can be applied to \emph{any} thermodynamic ensemble provided the ratio of path probabilities at different temperatures can be computed.
To illustrate the method, we compute a time-correlation function for solvated terminally-blocked alanine peptide across a range of temperatures using trajectories harvested using a modified parallel tempering protocol.

\emph{Keywords: Weighted histogram analysis method (WHAM); multiple histogram reweighting; multistate Bennett acceptance ratio (MBAR); extended bridge sampling (EBS); temperature reweighting; temperature-dependent dynamics; parallel tempering; molecular dynamics.}

\end{abstract}

%\pacs{02.30.Cj,05.70.Ce,82.20.Wt}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTRODUCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{section:introduction}

Dynamical properties, such as diffusion constants, position and velocity autocorrelation functions, rotational correlation times, frequency-dependent dielectric constants, and reaction or isomerization rates play a critical role in our understanding of various phenomena in chemistry and biology.
Besides providing physical insight, calculation of these properties from simulations is often necessary for making comparison with spectroscopic experiments (such as FTIR, 2DIR, NMR, dynamic light scattering, or neutron correlation spectroscopy) or to make predictions about material or structural properties under conditions difficult to access experimentally or of substances yet to be created in physical form.

When these phenomena involve enthalpic barriers or entropic bottlenecks in systems obeying classical statistical mechanics, averages of these properties can be extremely slow to converge in standard molecular dynamics simulations, requiring tediously long simulations to produce estimates with the desired statistical precision.
Additionally, the temperature-dependence of these properties is often of interest, requiring either the use of simulations at numerous temperatures, where each simulation must be long enough to ensure the difference in estimates between different temperatures is statistically significant, or the assumption of a kinetic model, such as Arrhenius behavior for rate processes.

In recent years, a number of algorithmic advances have helped ameliorate difficulties for the computation of equilibrium expectations or thermodynamic properties caused by the presence of significant enthalpic and entropic barriers.
Chief among these have been \emph{simulated}~\cite{marinari:europhys-lett:1992:simulated-tempering,lyubartsev:jcp:1992:expanded-ensembles,mitsutake:2000a} and \emph{parallel}~\cite{geyer-thompson:j-am-stat-assoc:1995:parallel-tempering,hukushima-nemoto:j-phys-soc-jpn:1996:parallel-tempering,hansmann:chem-phys-lett:1997:parallel-tempering-monte-carlo,sugita:chem-phys-lett:1999:REMD} \emph{tempering}, in which the time required to cross enthalpic (entropic) barriers is reduced by allowing the system to access higher (lower) temperatures during the simulation.
The use of random temperature-switching proposals and a Metropolis-like criterion for their acceptance embeds the simulations in a Markov chain that ensures, in the long run, the stationary distribution at thermal equilibrium is sampled at each temperature~\cite{jun-s-liu:mcmc}.
The advantages of this procedure are twofold: convergence times for averages at a \emph{single} temperature can often be reduced by appropriate choice of temperatures to enhance mixing~\cite{jun-s-liu:mcmc,hansmann:jcp:2006:optimized-parallel-tempering,machta:pre:2009:strengths-and-weaknesses-of-parallel-tempering,hummer:jcp:2009:error-and-efficiency-parallel-tempering}, and data from all temperatures can be combined to produce superior estimates of equilibrium expectations over a range of temperatures using histogram reweighting~\cite{ferrenberg:1989a,kumar:1992a,bartels:1997a,gallicchio:2005a} or histogram-free~\cite{bartels:cpl:2000:wham,shirts-chodera:jcp:2008:mbar} statistical analysis methods.

While helpful in the calculation of equilibrium expectations of mechanical observables that are functions only of the atomic coordinates and momenta, the data produced by simulated and parallel tempering simulations are generally not directly useful for the computation of dynamical quantities.
This is due to the fact that replica trajectories include unphysical changes in temperature, and permuting the data to separate trajectories by temperature results in discontinuous trajectories for each temperature, with a complicated correlation structure entangling the temperatures~\cite{gallicchio:2005a,chodera:jctc:2007:parallel-tempering-wham}.
The short trajectory segments generated between exchanges, however, are valid dynamical trajectories that can be used to estimate dynamical properties, as in~\cite{hummer:pre:2008:kinetics-from-replica-exchange}.

Here, we show how trajectories from multiple temperatures (including those harvested from tempering simulations) can be reweighted to produce optimal estimates of dynamical quantities at the temperature(s) of interest, and how these estimates might be superior to those from single-temperature simulation if the crossing of enthalpic barriers or entropic bottlenecks is accelerated at some of the replica temperatures.
Application of this reweighting scheme to simulated and parallel tempering simulations requires little, if any, modification of the simulation protocol.
We require only that the time between exchanges is long enough to compute dynamical expectations of interest, the model of dynamics used to propagate the replica samplers in between exchanges is of a form amenable to reweighting, and the temperature spacing is close enough to permit estimation of appropriate normalization constants.
We illustrate this approach by estimating a slowly-convergent property---the normalized fluctuation autocorrelation function for a conformational state of terminally-blocked solvated alanine peptide that is only sparsely populated at 300 K---and show that dynamical reweighting provides significant advantages over standard estimates from single temperatures alone.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THEORY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{section:theory}

We now lay out the main theoretical tools necessary for estimating dynamical expectations by making use of simulations at multiple temperatures.
In Section~\ref{section:dynamical-expectations}, we provide a precise definition of the dynamical expectations we can estimate through reweighting schemes.
Next, in Section~\ref{section:dynamical-reweighting}, we review the generalized path estimator for making optimal use of trajectories sampled from multiple ensembles.
As this estimator requires we compute path action differences at different temperatures for every trajectory sampled, Section~\ref{section:models-of-dynamics} presents simplified, convenient forms of these quantities for several common models of dynamics within the canonical ensemble.
Finally, Section~\ref{section:modified-parallel-tempering-protocol} describes a modified parallel tempering protocol that can be used to easily sample trajectories from multiple temperatures in a way that they can easily be used with this estimation-by-reweighting procedure.

\subsection{Dynamical expectations}
\label{section:dynamical-expectations}

The equilibrium thermodynamic expectation of some static (non-kinetic) property $A$ for a system obeying classical statistical mechanics can be written as
\begin{eqnarray}
\expect{A} \equiv \int dx \, p(x) \, A(x) \label{equation:thermodynamic-expectation}
\end{eqnarray}
where $x$ denotes the instantaneous configuration of the system (such as the coordinates and momenta of all particles), $A(x)$ is often referred to as a \emph{phase function} or \emph{mechanical observable}, and $p(x)$ is the equilibrium probability density, given by
\begin{eqnarray}
p(x) = Z^{-1} \, q(x) \:\:;\:\: Z \equiv \int dx \, q(x)
\end{eqnarray}
where $Z$ is a normalizing constant or \emph{partition function}, and $q(x) > 0$ is an unnormalized probability density.
(Here, we will choose $q(x)$ to contain terms that depend explicitly on $x$; $x$-independent multiplicative terms will generally be absorbed into $Z$ throughout.)
In the canonical ensemble, for example, the system is in contact with a heat bath and held at fixed volume, and we have $q(x) = e^{-\beta H(x)}$, where $H(x)$ is the Hamiltonian and $\beta = (k_B T)^{-1}$ is the inverse temperature.

We can write an analogous expression for equilibrium \emph{dynamical expectations} of a kinetic property $\mathcal{A}$ as
\begin{eqnarray}
\expect{\mathcal{A}} &=& \int dX \, p[X] \, \mathcal{A}[X] \label{equation:dynamical-expectation} 
\end{eqnarray}
where $\mathcal{A}[X]$ now denotes a \emph{functional} of a trajectory $X \equiv x(t)$, $p[X]$ is the probability density of trajectories, and the integral is taken over the space of all such trajectories with respect to an appropriate measure $dX$.
Analogous to the case of phase space probability densities, $p[X]$ can also be written in terms of an unnormalized density $q[X] > 0$,
\begin{eqnarray}
p[X] = Z^{-1} \, q[X] \:\:;\:\: Z \equiv \int dX \, q[X] \label{equation:trajectory-probability-density}
\end{eqnarray}
As we will see in Section~\ref{section:models-of-dynamics}, the precise definition of $p[X]$ will depend on the dynamical model under consideration.

While the expectation of \emph{any} trajectory functional can, at least formally, be computed this way, the most common dynamical quantities of interest are \emph{time-correlation functions} of the form
\begin{eqnarray}
C_{AB}(t) &\equiv& \expect{A(0) \, B(t)}
\end{eqnarray} 
for some pair of phase functions $A$ and $B$.
This corresponds to the choice
\begin{eqnarray}
\mathcal{A}[X] &\equiv& A(x(0)) \, B(x(t)) \label{equation:time-correlation-function}
\end{eqnarray}
For practical purposes, we will henceforth presume that the functional $\mathcal{A}[X]$ is \emph{temporally local}, in that it can be expressed in a way that operates on a finite time interval that can be bounded by some fixed duration $\tau$.
We can then restrict ourselves to considering trajectory objects $X \equiv x(t)$ of fixed length $\tau$, so that $x(t)$ is defined only on the interval $t \in [0,\tau]$.
In the case of time-correlation functions $C_{AB}(t)$ above, for example, $t \le \tau$.

Expectations with respect to some altered trajectory probability density $p_*[X]$ can also be computed within this equilibrium expectation framework by simply incorporating the weighting factor $p_*[X] / p[X]$ into the trajectory functional $\mathcal{A}$.
For example, to compute time-correlation functions such as Eq.~\ref{equation:time-correlation-function} with respect to \emph{nonequilibrium} initial conditions, $\mathcal{A}[X]$ could be redefined to include a factor of $\rho[X]/ p[X]$ to account for non-equilibrium (or even non-canonical) initial phase space density $\rho[X]$~\cite{skeel:siam-j-sci-comput:2008:what-makes-md-work}.

\subsection{Dynamical reweighting}
\label{section:dynamical-reweighting}

The groundwork for reweighting \emph{trajectories} sampled from multiple thermodynamic states in order to produce an optimal estimate of some dynamical expectation was laid in work by Minh and Chodera~\cite{minh-chodera:asymptotic-variance}.
Though presented in the context of nonequilibrium experiments (in which the consequences of the Crooks fluctuation theorem~\cite{crooks:pre:1999:fluctuation-theorem,crooks:2000a} were explored), the estimator is sufficiently general that it applies to equilibrium trajectories sampled from different equilibrium thermodynamic states within the same thermodynamic ensemble, producing an asymptotically optimal estimate of properties within the thermodynamic state of interest.
The generalized path ensemble estimator~\cite{minh-chodera:asymptotic-variance} is in turn based on the statistical inference framework of extended bridge sampling~\cite{meng-wong:1996:statistica-sinica:estimating-normalizing-constants,kong:j-royal-stat-soc-b:2003:bridge-sampling,tan:j-am-stat-assoc:2004:extended-bridge-sampling}, which provides a solid statistical foundation for earlier estimation and reweighting schemes found in statistical physics and chemistry~\cite{bennett:1977a,ferrenberg:1989a,kumar:1992a,bartels:1997a,bartels:cpl:2000:wham}.

Here, we briefly review the general estimator formalism and examine its application to common schemes used to model dynamics at constant temperature.
While we restrict our consideration to the canonical (NVT) ensemble, extension to other thermodynamic ensembles (such as the isobaric and semigrand-canonical ensembles) and other schemes for generating trajectories within the canonical and other thermodynamic ensembles is operationally straightforward.

Suppose we have $K$ path ensembles at different thermodynamic conditions, indexed by $i \in \{1,\ldots,K\}$, characterized by trajectory probability densities
\begin{eqnarray}
p_i[X] = Z_i^{-1} q_i[X] \:\:;\:\:Z_i \equiv \int dX \, q_i[X]
\end{eqnarray}
where $q_i[X] > 0$ is an unnormalized density and $Z_i$ an unknown normalization constant, from which we have collected $N_i$ trajectories of duration $\tau$.
We denote these trajectories $X_n$, where the trajectory index $n$ runs from $1$ to $N \equiv \sum_{k=1}^K N_k$, with the trajectories from different path ensembles indexed in \emph{arbitrary} order.
The association of a trajectory with the path ensemble from which it was sampled will not be relevant in the estimating equations.

The optimal estimator for a dynamical expectation
\begin{eqnarray}
\expect{\mathcal{A}}_i &=& \int dX \, p_i[X] \, \mathcal{A}[X]
\end{eqnarray}
was shown by Minh and Chodera~\cite{minh-chodera:asymptotic-variance} to be given by 
\begin{eqnarray}
\hat{\mathcal{A}}_i &=& \sum_{n=1}^N w_{ni} \, \mathcal{A}[X_{n}] \label{equation:dynamical-expectation-estimator}
\end{eqnarray}
where the $N \times K$ weight matrix $W \equiv (w_{ni})$ containing the appropriate trajectory weights for all trajectories $n$ in all path ensembles $i$ is given by
\begin{eqnarray}
w_{ni} &=& \hat{Z}_i^{-1} \left[ \sum_{k=1}^K N_k \hat{Z}_k^{-1} q_k[X_n] / q_i[X_n] \right]^{-1}. \label{equation:trajectory-weights}
\end{eqnarray}
Here, $w_{ni}$ denotes the weight contribution from trajectory $n$ in the aggregated pool of $N$ trajectories for estimating expectations for state $i$.
The presence of the leading normalizing factor $\hat{Z}_i^{-1}$ ensures that the weights are normalized such that $\sum_{n=1}^N w_{ni} = 1$ for any state $i$.
The unnormalized path densities $q_k[X_n] / q_i[X_n]$ represent the ratio of how likely a particular trajectory $X_n$ is to appear in path ensemble $k$ over path ensemble $i$, up to some ratio of normalizing constants $Z_i/Z_k$ that is trajectory-independent.
Importantly, the estimator expression in Eqs.~\ref{equation:dynamical-expectation-estimator}--\ref{equation:trajectory-weights} holds even if no trajectories are sampled from some of the path ensembles, such that $N_i = 0$ for these unsampled ensembles, but there are still samples from other path ensembles such  that $N > 0$).

The normalizing constants $\{\hat{Z}_i\}$, determined only up to an arbitrary multiplicative constant, are determined by solving the coupled set of $K$ nonlinear equations for $i = 1,\ldots,K$, under the constraint that $\sum_{n=1}^N w_{ni} = 1$:
\begin{eqnarray}
\hat{Z}_i &=& \sum_{n=1}^N \left[ \sum_{k=1}^K N_k \hat{Z}_k^{-1} q_k[X_n] / q_i[X_n] \right]^{-1} \label{equation:normalization-constants}
\end{eqnarray}
This can be done efficiently through a number of methods~\cite{minh-chodera:asymptotic-variance,shirts-chodera:jcp:2008:mbar}.
For example, the simplest such approach is to iterate a form of Eq.~\ref{equation:normalization-constants} to self-consistency.
Suppose we choose an initial guess $\hat{Z}^{(0)}_i = 0$, $i = 1,\ldots,K$.
We can employ an iterative update procedure,
\begin{eqnarray}
\hat{Z}^{(n+1)}_i &=& \sum_{n=1}^N \left[ \sum_{k=1}^K N_k [\hat{Z}^{(n)}_k]^{-1} q_k[X_n] / q_i[X_n] \right]^{-1} \label{equation:normalization-constant-iteration}
\end{eqnarray}
to generate a new set of estimates $\{\hat{Z}^{(n+1)}_i\}$ from a previous set of estimates $\{\hat{Z}^{(n)}_i\}$.
Each iteration of the update procedure requires only that we be able to compute the ratio $q_k[X_n] / q_i[X_n]$ for all trajectories in the pool of $N$ trajectories; these ratios can be precomputed and stored.
This iterative procedure is continued until these estimates converge to within some specified tolerance~\cite{minh-chodera:asymptotic-variance,shirts-chodera:jcp:2008:mbar}.
For numerical stability, it is convenient to work with $\ln Z_i^{(n)}$ instead of $Z_i^{(n)}$ directly~\cite{shirts-chodera:jcp:2008:mbar}.

In the canonical ensemble, the probability distributions are parameterized by a temperature $T$, or equivalently, the inverse temperature $\beta \equiv (k_B T)^{-1}$.
The expectation of observable $\mathcal{A}$ at some \emph{arbitrary} inverse temperature $\beta$, $\expect{\mathcal{A}}_\beta$, estimated from simulations at fixed inverse temperatures $\beta_1, \ldots, \beta_K$, can then be estimated as
\begin{eqnarray}
\hat{\mathcal{A}}(\beta) &=& \sum_{n=1}^N w_n(\beta) \, \mathcal{A}[X_n] \label{equation:reweighted-expectation}
\end{eqnarray}
where the temperature-dependent trajectory weights $w_n(\beta)$ are
\begin{eqnarray}
w_n(\beta) &=& \hat{Z}(\beta)^{-1} \left[ \sum_{k=1}^K N_k \hat{Z}_k^{-1} q[X_n | \beta_k] / q[X_n | \beta] \right]^{-1} \label{equation:temperature-dependent-trajectory-weights}
\end{eqnarray}
with normalization constants $Z_k \equiv Z(\beta_k)$, and 
\begin{eqnarray}
\hat{Z}(\beta) &=& \sum_{n=1}^N \left[ \sum_{k=1}^K N_k \hat{Z}_k^{-1} q[X_n | \beta_k] / q[X_n | \beta] \right]^{-1}
\end{eqnarray}
Once the normalizing constants $\{\hat{Z}_k\}$ have been determined by solving the coupled nonlinear equations in Eq.~\ref{equation:normalization-constants}, no further nonlinear equation solution iterations are necessary to estimate $\hat{\mathcal{A}}(\beta)$ at other temperatures of interest.

Perhaps surprisingly, Eqs.~\ref{equation:trajectory-weights} and \ref{equation:temperature-dependent-trajectory-weights} do not contain any information linking a trajectory $X_n$ with the path ensemble $p_j[X]$ from which it was sampled; this fact is a direct, if unobvious, consequence of the extended bridge sampling formalism~\cite{kong:j-royal-stat-soc-b:2003:bridge-sampling,tan:j-am-stat-assoc:2004:extended-bridge-sampling}.

The unnormalized trajectory probability density $q[X | \beta]$ will depend on the model of dynamics within the canonical ensemble used in the simulations; we give expressions for several popular models in Section~\ref{section:models-of-dynamics}.
In all the cases treated here, however, $q[X | \beta]$ is a continuous and differentiable function of $\beta$, meaning that $\mathcal{A}(\beta)$ will also be continuous and differentiable.
The utility of this estimate will depend on a number of contributing circumstances; the estimate is only expected to be reliable within the range of temperatures sampled, but estimates can, in principle, be obtained for any temperature.

The statistical uncertainty in $\hat{\mathcal{A}}(\beta)$ can be computed in a straightforward manner from an estimate of the asymptotic variance of $\hat{\mathcal{A}}(\beta)$~\cite{minh-chodera:asymptotic-variance,shirts-chodera:jcp:2008:mbar}.
Briefly, the estimating equations are linearized about the optimal estimator and the variance in the estimating equations is propagated into the corresponding variance in the estimator (see, for example, the Appendix in~\cite{tan:j-am-stat-assoc:2004:extended-bridge-sampling}.)
Operationally, the $N \times K$ weight matrix $W$ is augmented by two additional columns, indexed (for convenience) by $a$ and $A$, consisting of
\begin{eqnarray}
w_{na} = w_n(\beta) \:\: ; \:\: w_{nA} = \frac{\mathcal{A}[X_n]}{\hat{\mathcal{A}}(\beta)} w_{na}
\end{eqnarray}
where $\beta$ is the temperature of interest.

The variance in $\hat{A}(\beta)$ can be computed from the asymptotic covariance matrix $\bfm{\Theta}$, which is estimated by
\begin{eqnarray}
\hat{\bfm{\Theta}} &\equiv& \bfm{W}^T \cdot [\bfm{I}_N - \bfm{W} \bfm{N} \bfm{W}^\T]^+ \cdot \bfm{W} \label{equation:asymptotic-covariance}
\end{eqnarray}
where $\bfm{I}_N$ is the $N \times N$ identity matrix, $\bfm{N} = \mathrm{diag}(N_1, N_2, \ldots, N_K, 0, 0)$, and $\bfm{A}^+$ denotes the pseudoinverse of a matrix $\bfm{A}$.
For cases where the columns of $\bfm{W}$ are linearly independent (which will be the case when $\beta \neq \beta_i$ for all $i \in \{1,\ldots,K\}$ and $\beta_i \ne \beta_j$ for $i \ne j$), a simpler expression that requires pseudoinversion of only $\mathcal{O}(K) \times \mathcal{O}(K)$ can be used:
\begin{eqnarray}
\hat{\bfm{\Theta}} &=& [(\bfm{W}^\T \bfm{W})^+ - \bfm{N}]^+
\end{eqnarray}
Efficient computation of $\hat{\bfm{\Theta}}$ by other means is discussed in the Appendix of~\cite{shirts-chodera:jcp:2008:mbar}. 

The uncertainty in $\hat{\mathcal{A}}(\beta)$ is then estimated as
\begin{eqnarray}
\var \hat{\mathcal{A}}(\beta) &\equiv& \expect{\left(\hat{\mathcal{A}}(\beta) - \expect{\hat{\mathcal{A}}(\beta)}\right)^2} \nonumber \\
&\approx& [\hat{\mathcal{A}}(\beta)]^2 \left[ \hat{\Theta}_{aa}  - 2 \hat{\Theta}_{aA} + \hat{\Theta}_{AA} \right] \label{equation:uncertainty-estimate}
\end{eqnarray}
Near $\hat{A}(\beta) \equiv 0$, the uncertainty estimation in Eq.~\ref{equation:reweighted-expectation} may run into numerical issues; in this case, it is recommended that the relative uncertainty $\var{\hat{A}(\beta)} / [\hat{A}(\beta)]^2$ be computed instead.

Similarly, the covariance between two estimates $\hat{\mathcal{A}}(\beta)$ and $\hat{\mathcal{B}}(\beta')$ at respective temperatures $\beta$ and $\beta'$ (which may be the same or different) can be computed by augmenting $\mathrm{W}$ with additional columns~\cite{minh-chodera:asymptotic-variance,shirts-chodera:jcp:2008:mbar},
\begin{eqnarray}
w_{nb} = w_n(\beta') \:\: ; \:\: w_{nB} = \frac{\mathcal{B}[X_n]}{\hat{\mathcal{B}}(\beta')} w_{nb}
\end{eqnarray}
yielding
\begin{eqnarray}
\lefteqn{\mathrm{cov} \left( \hat{\mathcal{A}}(\beta),  \hat{\mathcal{B}}(\beta') \right) \:\: } \nonumber \\
&\equiv& \expect{\left(\hat{\mathcal{A}}(\beta) - \expect{\hat{\mathcal{A}}(\beta)}\right) \left(\hat{\mathcal{B}}(\beta') - \expect{\hat{\mathcal{B}}(\beta')}\right)} \nonumber \\
&\approx& \hat{\mathcal{A}}(\beta) \, \hat{\mathcal{B}}(\beta') \left[ \hat{\Theta}_{ab}  - \hat{\Theta}_{aB} - \hat{\Theta}_{Ab} + \hat{\Theta}_{AB} \right] \label{equation:correlated-uncertainty}
\end{eqnarray}
This can be extended to more than two expectations if the covariance of many such pairs is desired simultaneously~\cite{minh-chodera:asymptotic-variance,shirts-chodera:jcp:2008:mbar}.

For purposes of numerical stability, it is convenient to work with the logarithms of quantities like $q[X | \beta_i] / q[X | \beta_j]$ and $\hat{Z}_i / \hat{Z}_j$.
We therefore define a trajectory \emph{action} functional $S[X | \beta]$ as
\begin{eqnarray}
S[X | \beta] &\equiv& - \ln q[X | \beta]
\end{eqnarray}
Note that this definition differs from standard definitions (e.g.~\cite{chandler:1998:cqdcps:tps}) in that only terms that depend on \emph{both} trajectory $X$ and temperature $\beta$ are included in $q[X | \beta]$ and $S[X | \beta]$---any remaining multiplicative terms in the path probability density $p[X | \beta]$ are subsumed by the normalization constant $Z(\beta)$ to simplify the subsequent development.
Inclusion of explicit multiplicative terms that depend on $\beta$ but not $X$ could reduce the variance of dynamical reweighting estimates, but this is not done here for simplicity.

For all of the models of dynamics considered in Section~\ref{section:models-of-dynamics}, the trajectory- and temperature-dependent action $S[X | \beta]$ is \emph{linear} in the inverse temperature $\beta$, having the form
\begin{eqnarray}
S[X | \beta] &=& \beta \mathcal{H}[X]
\end{eqnarray}
where the quantity $\mathcal{H}[X]$ plays the role of a generalized path Hamiltonian that operates on trajectory space, analogous to the role of the standard Hamiltonian $H(x)$ that operates on microstates in phase space.
We note that the path Hamiltonian can be trivially computed from the action $S[X|\beta]$ or unnormalized path density $q[X|\beta]$ as
\begin{eqnarray}
\mathcal{H}[X] &=& \beta^{-1} S[X | \beta] = - \beta^{-1} \ln q[X | \beta] \label{equation:path-hamiltonian}
\end{eqnarray}

We also define the dimensionless path free energy $f_i$ of path ensemble $i$ as
\begin{eqnarray}
f_i &\equiv& - \ln Z_i
\end{eqnarray}
The logarithms of the trajectory weights $w_n(\beta)$ can then be written in terms of the path Hamiltonians $\mathcal{H}[X]$ and dimensionless free energies $f_i$ as
\begin{eqnarray}
\ln w_n(\beta) &=& \hat{f}(\beta) - \ln \sum_{k=1}^K N_k e^{ \hat{f}_k - (\beta_k - \beta) \mathcal{H}[X_n] }
\end{eqnarray}
where the temperature-dependent dimensionless free energy $\hat{f}(\beta)$ is given by
\begin{eqnarray}
\hat{f}(\beta) &=& - \ln \sum_{n=1}^N \left[\sum_{k=1}^K N_k e^{ \hat{f}_k - (\beta_k - \beta) \mathcal{H}[X_n] }\right]^{-1}
\end{eqnarray}
and the $\hat{f}_k \equiv \hat{f}(\beta_k)$.
Practical notes about working with the estimation equations in a numerically stable and efficient manner are given in~\cite{minh-chodera:asymptotic-variance,shirts-chodera:jcp:2008:mbar}.
Due to the analogous nature of estimating static equilibrium expectations using the Bennett acceptance ratio (MBAR) method~\cite{shirts-chodera:jcp:2008:mbar}, {\tt pymbar}, the Python package for performing MBAR analysis (available at \url{http://simtk.org/home/pymbar}), can be used for the computation of dimensionless path free energies $\hat{f}_i$, expectations $\hat{A}$, and covariances $\mathrm{cov} (\hat{A},\hat{B})$ by using path actions $S[X]$ in place of the reduced potential energy $u(x)$.
This scheme was used for the calculations presented here.

\subsection{Models of dynamics in the canonical ensemble}
\label{section:models-of-dynamics}

While the path ensemble estimator (described for temperature-dependent dynamics in Section~\ref{section:dynamical-reweighting}, and more generally in~\cite{minh-chodera:asymptotic-variance}) can be applied to dynamics in \emph{any} thermodynamic ensemble in which the unnormalized trajectory probability density ratios $q_i[X] / q_j[X]$ are finite and nonzero over the same domain, we restrict our consideration to the canonical (NVT) ensemble in classical statistical mechanics.
As the concept of a thermodynamic ensemble is a purely \emph{equilibrium} construct, it specifies only the relative probability $p(x | \beta)$ with which a configuration or phase space point $x$ is observed at equilibrium.
To write a probability density over equilibrium trajectories at inverse temperature $\beta$, $p[X | \beta]$, we must further specify a particular scheme for modeling dynamics within this ensemble.
The choice of model used to simulate dynamics within this ensemble is not unique---many schemes can be employed that will generate the same equilibrium density in phase space $p(x | \beta)$, but will have different trajectory probability densities $p[X | \beta]$.
This choice must therefore be guided by both a desire to mimic the relevant physics for the system of interest (such as whether the system must be allowed to exchange energy with an external heat reservoir during the course of dynamical evolution), balanced with computational convenience (e.g.~the use of a stochastic thermostat in place of explicitly simulating a large external reservoir).

We consider four common models for simulating dynamics within the NVT ensemble for which the ratio of unnormalized trajectory probability densities $q[X | \beta_i] / q[X | \beta_j]$ can be computed: \emph{(i)} Hamiltonian dynamics with canonically-distributed initial conditions, \emph{(ii)} an Andersen thermostat~\cite{andersen:jcp:1980:andersen-thermostat}, \emph{(iii)} Langevin dynamics, and \emph{(iv)} overdamped Langevin (sometimes called \emph{Brownian}) dynamics.
Notably, the deterministic nature of the Nos\'{e}-Hoover thermostat~\cite{hoover:pra:1985:nose-hoover,evans:jcp:1985:nose-hoover} and the way that the inverse temperature $\beta$ appears in its equations of motion means a trajectory $X \equiv x(t)$ generated at inverse temperature $\beta_i$ will have identically zero probability $p[X | \beta_j]$ at inverse temperature $\beta_j \ne \beta_i$, preventing its use in the reweighting scheme described here.
The same is true of the Berendsen weak-coupling algorithm~\cite{berendsen:1984:jcp:weak-coupling-thermostat}, though any usage of this form of thermal control is highly discouraged because of its failure to generate a canonical ensemble~\cite{morishita:2000:jcp:berendsen-analysis,schmidler:2008:jcp:boltzmann-in-replica-exchange,hummer:2009:jctc:thermostat-artifacts} or be ergodic \cite{schmidler:2008:jcp:boltzmann-in-replica-exchange}.

\subsubsection{Canonical distribution of Hamiltonian trajectories}
\label{section:canonical-distribution-of-hamiltonian-trajectories}

Consider a collection of trajectories of length $\tau$ in which initial phase space points $x_0 \equiv x(0)$ are sampled from the canonical density $p(x_0 | \beta)$, and whose dynamical evolution is governed by Hamiltonian mechanics.
Physically, this corresponds to the statistics of a situation where the system is initially in contact with a thermal reservoir at inverse temperature $\beta$, but the evolution time $\tau$ is short enough that effectively no heat is exchanged with the reservoir during this time.
Practically, an ensemble of these trajectories can be generated in a number of ways: 
For example, the initial configurations $x_0$ could be generated from the desired canonical ensemble $p(x_0 | \beta) \propto e^{-\beta H(x_0)}$ by a form of Metropolis-Hastings Monte Carlo~\cite{hastings:1970:biometrika:metropolis-hastings} or hybrid Monte Carlo~\cite{duane:1987:phys-lett-b:hmc}.

The trajectory probability density $p[X | \beta]$ is given by
\begin{eqnarray}
p[X | \beta] &=& Z(\beta)^{-1} \, e^{-\beta H(x_0)}
\end{eqnarray}
where $H(x)$ denotes the Hamiltonian for the system and $x_0 \equiv x(0)$ is the initial phase space point.
While the normalization constant can be written
\begin{eqnarray}
Z(\beta) &=&\int dx \, e^{-\beta H(x)} \label{equation:canonical-ensemble}
\end{eqnarray}
we do not need to compute it for use in dynamical reweighting.

The path Hamiltonian $\mathcal{H}[X]$ can then be identified (via Eq.~\ref{equation:path-hamiltonian}) as being identical to the Hamiltonian 
\begin{eqnarray}
\mathcal{H}[X] &=& H(x_0) \label{equation:hamiltonian-action} .
\end{eqnarray}
We note that, while the Hamiltonian is invariant along the trajectory for true Hamiltonian dynamics, numerical integrators will not exactly preserve this property.

\subsubsection{Andersen thermostat}
\label{section:andersen-thermostat}

The Andersen thermostat~\cite{andersen:jcp:1980:andersen-thermostat,e:commun-pure-appl-math:2008:andersen-thermostat} allows  the system to exchange heat with an external thermal reservoir through stochastic collisions with virtual bath particles.
These ``collisions'' simply cause the real particle's velocity to be reassigned from a Maxwell-Boltzmann distribution corresponding to the thermostat temperature, ensuring that the canonical ensemble is sampled~\cite{andersen:jcp:1980:andersen-thermostat}.

Two schemes are commonly used.
In the first scheme, the system undergoes so-called \emph{massive collisions}, in which all momenta are reassigned from the Maxwell-Boltzmann distribution at fixed intervals; evolution between collisions occurs by normal Hamiltonian dynamics.
This collision interval can have a large effect on the kinetic properties of the system, so this scheme is usually used only for studying static thermodynamic properties.

In the second scheme, each \emph{particle} has a fixed probability per unit time $\nu$ of having its momenta reassigned.
This corresponds to a scenario in which each particle undergoes independent stochastic collisions with fictitious bath particles, with the system evolving according to Hamiltonian dynamics between collisions.
Andersen suggested a physically-motivated collision frequency $\nu$ could be selected based on the thermal conductivity $\kappa$ of the sample~\cite{andersen:jcp:1980:andersen-thermostat}
\begin{eqnarray}
\nu &=& \frac{2 a \kappa}{3 k_B \rho^{1/3} {N_p}^{2/3}}
\end{eqnarray}
where $a$ is a dimensionless constant that depends on the shape of the system, $\rho$ is the density, and $N_p$ is the number of particles.
Note that, as the system size increases with $N_p$, the collision frequency $\nu$ decreases as $N_p^{-2/3}$, such that the dynamics will approach that of a Hamiltonian system for large systems.
Further elaborations of this thermostatting scheme attempt to model physically realistic dynamics of large solutes, such as solvated biological macromolecules, by coupling the thermostat only to solvent degrees of freedom (see, e.g., Ref.~\cite{bolhuis:pnas:2006:trp-cage}).
Systems with constraints will require the velocity components along the constraints to be removed prior to integration by an algorithm such as RATTLE~\cite{andersen:1983:j-comput-phys:rattle}.

For either thermostatting scheme, the probability of sampling a trajectory $X \equiv x(t)$ at inverse temperature $\beta$ is therefore determined by the probability of selecting the initial configuration $x_0 \equiv x(0)$ from equilibrium and the probability of generating the resulting sequence of random collisions.
Up to irrelevant temperature-independent multiplicative constants that simply determine the probability that a given sequence of $N_c$ particle collisions occurred, this probability is
\begin{eqnarray}
p[X | \beta] &\propto& e^{-\beta H(x_0)} \prod_{n=1}^{N_c} \pi_v(\bfm{v}_n | m_n, \beta)
\end{eqnarray}
where the index $n$ in the product runs over all particle collisions that occurred during integration (either at regular intervals or stochastically determined), indexed in arbitrary order.
$\bfm{v}_n$ is the velocity of the particle after the $n$th massive collision, $m_n$ its mass, and $\pi_v(\bfm{v} | m, \beta)$ is the Maxwell-Boltzmann distribution for the velocity vector of a particle of mass $m$ at inverse temperature $\beta$:
\begin{eqnarray}
\pi_v(\bfm{v} | \beta, m) &\propto& e^{-\beta m v^2 / 2}
\end{eqnarray}

The path Hamiltonian $\mathcal{H}[X]$ can therefore be identified as the Hamiltonian for the initial phase space point $x_0$ plus the new kinetic energies of all particles that have undergone collisions during the course of integration:
\begin{eqnarray}
\mathcal{H}[X] &=& H(x_0) + \frac{1}{2} \sum_{n=1}^{N_c} m_n v_n^2
\end{eqnarray}
where again, $\mathrm{v}_n$ is the new velocity of the particle that had undergone collision indexed by $n$, and $m_n$ is its mass.

\subsubsection{Langevin dynamics}
\label{section:langevin-dynamics}

In Langevin dynamics, stochastic forces are used to mimic the influence of degrees of freedom that are not explicitly represented by particles in the system, such as solvent molecules whose influence is incorporated into the potential by mean-field interactions.
The memoryless Langevin equations of motion for phase space point $x \equiv (\bfm{r},\bfm{v})$, consisting of Cartesian atomic positions $\bfm{r}$ and velocities $\bfm{v}$, are 
\begin{eqnarray}
\dot{r}_i(t) &=& v_i(t) \nonumber \\
m_i \dot{v}_i(t) &=& F_i(\bfm{v}(t)) - \gamma_i v_i(t) + R_i(t)
\end{eqnarray}
where $m_i$ is the mass associated with degree of freedom $i$, $F_i(\bfm{q}) \equiv -\partial H / \partial q_i$ is the systematic force component, and $\gamma_i$ is an effective collision frequency or friction coefficient for that degree of freedom, with units of inverse time (e.g.~91 ps$^{-1}$ for water-like viscosity).
$\gamma_i$ is related to the temperature-dependent diffusion constant $D_i(\beta)$ by 
\begin{eqnarray}
D_i(\beta) &=& (\beta \gamma_i m_i)^{-1} \label{equation:diffusion-constant}
\end{eqnarray}
The stochastic force $R_i(t)$ has zero mean and satisfies the fluctuation-dissipation theorem,
\begin{eqnarray}
\expect{R_i(t) \, R_j(t')} &=& 2 \beta^{-1} m_i \gamma_i \delta_{ij} \delta(t - t') 
\end{eqnarray}

The Langevin Leapfrog integrator \cite{izaguirre:2010:psb:normal-mode-langevin} is an accurate and stable algorithm for integrating the Langevin equations of motion.
Updating of positions $\bfm{r}_t$ and velocities $\bfm{v}_t$ by a discrete timestep $\Delta t$ is handled by the following scheme:
\begin{eqnarray}
&v_{i,t+\frac{1}{2}} &= a_i v_{i,t} + b_i \Delta t \, m_i^{-1} F_i(X_t) + c_i m_i^{-1/2} \xi_{it} \nonumber \\
&r_{i,t+1} &= r_{i,t} + \Delta t \, v_{i,t+\frac{1}{2}} \label{equation:langevin-leapfrog} \\
&v_{i,t+1} &= a_i v_{i,t+\frac{1}{2}} + b_i \Delta t \, m_i^{-1} F_i(X_{t+1}) + c_i m_i^{-1/2} \xi_{i(t+1)} \nonumber
\end{eqnarray}
where we have defined the dimensionless constants $a_i$, $b_i$, and $c_i$ associated with each degree of freedom $i$ as
\begin{eqnarray}
a_i &\equiv& e^{-\gamma_i \Delta t / 2} \nonumber \\
b_i &\equiv& (1 - e^{-\gamma_i \Delta t / 2}) / (\gamma_i \Delta t) \nonumber \\
c_i &\equiv& (1 - e^{-\gamma_i \Delta t})^{1/2}  \nonumber
\end{eqnarray}
As the collision rate $\gamma_i \rightarrow 0$, we have $a_i \rightarrow 1$, $b_i \rightarrow 1/2$, and $c_i \rightarrow 0$, resulting in the standard leapfrog integrator scheme.

$\xi_{it}$ is a unit-bearing random normal variate with zero mean and variance $\beta^{-1} = k_B T$, simply denoted
\begin{eqnarray}
\xi_{it} &\sim& \mathcal{N}(0, \beta^{-1}) .
\end{eqnarray}
Note that $\xi_{it}$ has units of the square-root of energy.

Propagation for $N_t$ steps of Langevin leapfrog integration from initial configuration $x_0 \equiv (\bfm{r}_0, \bfm{v}_0)$ for a system with $N_d$ degrees of freedom requires generating the sequence of random variates $\bfm{\xi}_{it}$.
Ignoring multiplicative factors not involving both trajectory $X$ and temperature $\beta$, the unnormalized equilibrium probability of sampling a trajectory $X$ which originates at $x_0$ and has noise history $\xi_{it}$ at inverse temperature $\beta$ is given by
\begin{eqnarray}
q[X | \beta] &=& e^{-H(x_0)} \prod_{t=0}^{N_t} \prod_{i=1}^{N_d} \phi(\xi_{it} ; 0, \beta^{-1})
\end{eqnarray}
where $\phi(x; \mu, \sigma^2)$ is the normal probability density
\begin{eqnarray}
\phi(y; \mu, \sigma^2) &\propto& \exp\left[ -\frac{1}{2 \sigma^2}(y - \mu)^2 \right]
\end{eqnarray}

The path Hamiltonian can therefore be computed from the initial phase space point $x_0$ and unit-bearing noise history $\xi_{it}$, as
\begin{eqnarray}
\mathcal{H}[X] &=& H(x_0) + \frac{1}{2} \sum_{t=0}^{N_t} \sum_{i=1}^{N_d} \xi_{ti}^2
\end{eqnarray}
which can easily be accumulated during integration.

Note that, if $\gamma_i = 0$ for a degree of freedom $i$, the corresponding noise term $\xi_{ti}^2$ does not contribute to the integration in Eq.~\ref{equation:langevin-leapfrog} and can be omitted from the path Hamiltonian above.

%{\color{red} [JDC: Include Tony Lelievre's trick where we hold $\gamma/\beta$ constant so that we can approach small timestep limit without the generalized heat capacity $\mathcal{C}$ blowing up.]}

\subsubsection{Overdamped Langevin (Brownian) dynamics}
\label{section:brownian-dynamics}

At high friction, momentum relaxation becomes fast compared to particle motion, such that inertial motion may be neglected.
A common integrator for this overdamped regime, in which only coordinates $\bfm{q}$ are explicitly integrated, is given by Ermak and Yeh~\cite{ermak-yeh:cpl:1974:brownian-dynamics,ermak:jcp:1975:brownian-dynamics}
\begin{eqnarray}
r_{t+1,i} &=& r_{ti} + \frac{\Delta t}{\gamma_i m_i} F_i(x_t) + \sqrt{2} \left( \frac{\Delta t}{\gamma_i m_i} \right)^{1/2} \xi_{ti}
\end{eqnarray}
where the stochastic component $\xi_{ti}$ is again a unit-bearing normal random variate with zero mean and variance $\beta^{-1} = k_B T$,
\begin{eqnarray}
\xi_{ti} \sim \mathcal{N}\left(0, \beta^{-1}\right)
\end{eqnarray}
with units of the square-root of energy, and $\gamma_i$ is a collision rate of friction constant for degree of freedom $i$, with units of inverse time, where the temperature-dependent diffusion constant $D_i$ is again related by Eq.~\ref{equation:diffusion-constant}.

Ignoring irrelevant multiplicative factors that do not depend on both trajectory $X$ and temperature $\beta$, the unnormalized probability functional for a realization of $N_t$ steps of this process is then
\begin{eqnarray}
q[X | \beta] &=& e^{-\beta U(x_0)} \prod_{t=1}^{N_t} \prod_{i=1}^{N_d} \phi(\xi_{ti}; 0, \beta^{-1}) 
\end{eqnarray}
where $U(x) = U(\bfm{r})$ denotes the potential energy and $\phi(y ; \mu, \sigma^2)$ again denotes the normal probability density.
This allows us to identify the path Hamiltonian as
\begin{eqnarray}
\mathcal{H}[X] &=& U(x_0) + \frac{1}{2} \sum_{t=1}^{N_t} \sum_{i=1}^{N_d} \xi_{ti}^2
\end{eqnarray}

\subsection{Modified parallel tempering protocol}
\label{section:modified-parallel-tempering-protocol}

An especially convenient way to harvest trajectories from multiple temperatures is through the use of parallel tempering~\cite{geyer-thompson:j-am-stat-assoc:1995:parallel-tempering,hukushima-nemoto:j-phys-soc-jpn:1996:parallel-tempering,hansmann:chem-phys-lett:1997:parallel-tempering-monte-carlo,sugita:chem-phys-lett:1999:REMD}.
As illustrated by Buchete and Hummer~\cite{hummer:jpcb:2008:coarse-master-equations}, the trajectory segments generated in between exchange attempts are valid dynamical trajectories sampled from the corresponding replica temperatures.
Simulated tempering can likewise be employed, though we will not discuss this here.

Below, we enumerate several considerations that must be taken into account for the use of a parallel tempering protocol to collect trajectories for use in dynamical reweighting.

First, the dynamics in between exchange attempts must use a scheme amenable to reweighting, such as one of the models of dynamics within the canonical ensemble described in Section~\ref{section:models-of-dynamics}.
While this means some popular choices, such as Nos\'{e}-Hoover dynamics~\cite{hoover:pra:1985:nose-hoover,evans:jcp:1985:nose-hoover} and the Berendsen weak-coupling algorithm~\cite{berendsen:1984a} cannot be used, considerable flexibility remains in which dynamical schemes are permitted, including Hamiltonian, Andersen, Langevin, and overdamped Langevin (Brownian) dynamics, as discussed in the previous section.

Second, the time between exchange attempts should be at least as long as $\tau$, the time over which the computed dynamical observable $\mathcal{A}[X]$ is temporally local.
For example, if a correlation function $C_{AB}(t) \equiv \expect{A(0) B(t)}$ is desired for $t \in [0, \tau]$, the time between exchanges should be at least $\tau$.
This may, in some circumstances, impact the efficiency of parallel tempering in producing uncorrelated samples.
For functions $\mathcal{A}$ that require a large $\tau$, much of the benefit of the exchanges in parallel tempering may be lost if few exchange attempts occur during the simulation.
Often, the best choice is to make this time \emph{exactly} $\tau$, since the enhanced mixing properties of parallel tempering simulations diminish as the time between exchanges grows (for fixed total simulation time)~\cite{roitberg:jcp:2008:exchange-frequency}.

Thirdly, it is necessary that there is sufficient overlap in the distribution of path Hamiltonians sampled from neighboring temperatures for dynamical reweighting to be effective.
While the standard replica-exchange acceptance criteria~\cite{sugita:chem-phys-lett:1999:REMD} will lead to the correct distribution of trajectories at each temperature, the fraction of accepted exchange attempts reports on the quality of overlap of potential energies, rather than the overlap of path Hamiltonians; as a result, poor reweighting performance may be obtained despite a high exchange acceptance fraction.
A simple way to ensure both good overlap in path Hamiltonian between neighboring temperatures and correct sampling of the equilibrium distribution is to use a modified exchange criterion based on the trajectory action:
\begin{eqnarray}
%\lefteqn{P_\mathrm{exch}(\beta_i, \beta_j) = \min \left\{ 1, \frac{p[X_i | \beta_j] p[X_j | \beta_i]}{p[X_i | \beta_i] p[X_j | \beta_j]} \right\}} \nonumber \\
%&=&\min\left\{ 1, e^{-(\beta_j - \beta_i) (\mathcal{H}[X_i] - \mathcal{H}[X_j])} \right\} \label{equation:modified-acceptance-criterion}
P_\mathrm{exch}(\beta_i, \beta_j) &=& \min \left\{ 1, \frac{p[X_i | \beta_j] p[X_j | \beta_i]}{p[X_i | \beta_i] p[X_j | \beta_j]} \right\} \nonumber \\
&=&\min\left\{ 1, e^{-(\beta_j - \beta_i) (\mathcal{H}[X_i] - \mathcal{H}[X_j])} \right\} \label{equation:modified-acceptance-criterion}
\end{eqnarray}
where, again, $\mathcal{H}[X] = \beta^{-1} S[X | \beta] = - \beta^{-1} \ln q[X | \beta]$ denotes the path Hamiltonian for trajectory $X$ defined in Eq.~\ref{equation:path-hamiltonian}.
The resulting procedure can be considered a form of parallel tempering transition path sampling \cite{vlugt-smit:2001:phys-chem-comm:parallel-tempering-tps}, with the difference that new trajectories are generated from old ones in the more conventional fashion of integrating equations of motion from the final configuration $x_T$ of the previous trajectory, rather than some form of transition path sampling~\cite{chandler:1998:cqdcps:tps,bolhuis:annu-rev-phys-chem:2002:transition-path-sampling}.
Note that, even for the case of a canonical distribution of Hamiltonian trajectories, the modified acceptance criteria of Eq.~\ref{equation:modified-acceptance-criterion} differs from the standard acceptance criteria~\cite{sugita:chem-phys-lett:1999:REMD}, as the total Hamiltonian in the modified criteria replaces the potential energy in the standard criteria.

Attempting to swap temperatures $\beta_i$ and $\beta_j$ for an arbitrarily chosen pair of replicas $(i,j)$ will ensure that the trajectories are correctly distributed from equilibrium at their new (or old) temperatures after the exchange attempt guarantees each replica visiting temperature $\beta_i$ samples from $p[X | \beta_i]$ in the long run~\cite{geyer:1991:parallel-tempering,hukushima-nemoto:1996:j-phys-soc-japan:parallel-tempering,jun-s-liu:mcmc}; as a result, their individual timeslices $x_t$ will also be distributed from the corresponding equilibrium distribution $p(x | \beta_i)$.
If there is poor overlap between the trajectory action distributions, then the exchange attempts between these temperatures will rarely succeed; on the other hand, if the overlap is good, then the exchange rate will be accepted with significant probability.
Because the cost of a single exchange attempt is insignificant compared to the computational cost of propagating all $K$ replicas by another time $\tau$, it is recommended that many swaps between pairs $(i,j)$ chosen with uniform probability be attempted in order to ensure the replicas are well-mixed.
Monitoring of attempted and accepted exchanges early in the simulation can help diagnose whether there is sufficient coupling between all temperatures for successful reweighting.

While it is not strictly necessary to employ this modified acceptance probability, some care must be taken to ensure that a canonical distribution is actually obtained.  
For example, swapping using an exchange criteria based only on potential energies (as suggested by Sugita et al.~\cite{sugita:chem-phys-lett:1999:REMD}) and failing to rescale velocities may lead to an incorrect distribution.
On the other hand, rescaling the velocities by $(\beta_\mathrm{old} / \beta_\mathrm{new})^{1/2}$ as recommended~\cite{sugita:chem-phys-lett:1999:REMD}, or subsequently redrawing the initial velocities from a Maxwell-Boltzmann distribution at the new temperature~\cite{hansmann:chem-phys-lett:1997:parallel-tempering-monte-carlo,pitera:2006a}---a process equivalent to ``massive collisions'' for the Andersen thermostat, is not problematic.
However, extra care must be taken to ensure the temperatures are spaced such that there is sufficient overlap in trajectory action between neighboring temperatures.

Finally, we note that the temporally sequential trajectories produced by each replica will have some degree of correlation in the observable $\mathcal{A}[X]$.  
If this correlation is significant, the uncertainty estimate produced by Eqs.~\ref{equation:uncertainty-estimate} will be an underestimate of the true statistical error, as the samples used in the estimation are assumed to be \emph{uncorrelated}.
A simple way to remove this correlation is to construct the timeseries $A_t \equiv \mathcal{A}[X_t]$ over the sequentially-sampled trajectories $X_t$ from a single replica (\emph{without} having permuted the replicas to collect trajectory data by temperature) and estimate the statistical inefficiency $g_\mathcal{A}$ from the integrated autocorrelation time of the timeseries $A_t$, $t = 0, 1, \ldots$~\cite{swope:1982a,flyvbjerg:1989a,janke:2002a,chodera:jctc:2007,shirts-chodera:jcp:2008:mbar}.
If an equally-spaced subset of $N / g_\mathcal{A}$ trajectories from the replicas is generated, these samples will be effectively uncorrelated.
Alternatively, the time between exchanges $\tau$ can be increased to ensure that sequential $A_t$ are decorrelated.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPLICATION TO ALANINE DIPEPTIDE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Correlation functions of a solvated terminally-blocked alanine peptide}
\label{section:application-to-alanine-dipeptide}

\begin{figure}[tb]
\resizebox{\columnwidth}{!}{\includegraphics{figure1}}
\caption{{\bf Terminally-blocked alanine in explicit solvent and potential of mean force of $\psi$ torsion at 300 K.}
The dark shaded region denotes a 95\% confidence interval (two standard errors) in the free energy difference between the plotted $\psi$ torsion and the lowest point on the potential of mean force (near $\psi = 150$ degrees).
The labeled light shaded region denotes the $\alpha_R$ conformation.
The PMF and error bars were estimated from the parallel tempering dataset using MBAR as described in the text.
\label{figure:alanine-dipeptide-psi-pmf}} 
\end{figure}

%\begin{figure}[tb]
%%\resizebox{\columnwidth}{!}{\includegraphics{figures/alanine-dipeptide/mixing-statistics}}
%\caption{{\bf Mixing statistics among temperatures.}
%The probability that a replica at temperature index $i$ is found in temperature index $j$ on the subsequent iteration of the parallel tempering simulation is depicted.
%\label{figure:alanine-dipeptide-mixing-statistics}} 
%\end{figure}

%\begin{figure}[tb]
%\resizebox{\columnwidth}{!}{\includegraphics{figures/alanine-dipeptide/2d-pmf-placeholder}}
%\caption{{\bf Terminally-blocked alanine in explicit solvent and potential of mean force in $(\phi,\psi)$ torsions at 300 K.}
%{\color{red}[JDC: The current figure is just a placeholder, and must be regenerated with new dataset.  Only the $C_7^\mathrm{ax}$ state boundaries (state 5) will be marked.]}
%\label{figure:alanine-dipeptide-2d-pmf}} 
%\end{figure}

To illustrate the application of dynamical reweighting to a condensed-phase system, we consider the estimation of the normalized fluctuation autocorrelation function for a long-lived conformational state of terminally-blocked alanine peptide in water that is sparsely populated at 300 K.
%This system was previously studied by Bolhuis et al.~\cite{bolhuis:pnas:2000:TPS-alanine-dipeptide} using transition path sampling~\cite{chandler:jcp:1998:transition-path-sampling,bolhuis:annu-rev-phys-chem:2002:transition-path-sampling}.
%We use the same solvated, terminally-blocked alanine peptide system described in Ref.~\cite{chodera:mms:2006}.
Additional simulation details are given in Appendix~\ref{appendix:simulation-details}.
All code, datasets, analysis scripts, and plotting scripts used for this application are made available for download at \url{http://simtk.org/home/dynamical-reweighting}.

The modified parallel tempering protocol described in Section~\ref{section:modified-parallel-tempering-protocol} was used to sample many 10 ps trajectories from each of 40 exponentially-spaced temperatures spanning 300--600 K, using the ``canonical distribution of Hamiltonian trajectories'' model of dynamics described in Section~\ref{section:canonical-distribution-of-hamiltonian-trajectories}.
Replica temperatures were chosen to be exponentially spaced so that $T_k = (T_\mathrm{max} - T_\mathrm{min}) \exp[(k-1)/(K-1)] + T_\mathrm{min}$ for $k = 1,\ldots,K$; this is equivalent to geometrically-spaced replicas, which would achieve equal acceptance rates among replicas if the system being simulated was an ideal gas.
\color{red}
These temperatures may not be optimal in terms of reducing replica correlation times or equalizing acceptance probabilities for the system studied here.
While numerous schemes for intelligently choosing or adapting temperatures exist~\cite{kofke:2002:jcp:acceptance-probability,katzberger-trebst-huse-troyer:j-stat-mech:2006:feedback-optimized-parallel-tempering,trebst-troyer:2006:j-stat-mech:feedback-optimized-monte-carlo,trebst-troyer-hansmann:jcp:2006:optimized-replica-selection,nadler-hansmann:pre:2007:generalized-ensemble,gront-kolinski:j-phys-condens-matter:2007:optimized-replica-selection,park-pande:pre:2007:choosing-weights-simulated-tempering}, no attempt was made to do so here.
\color{black}
The modified protocol makes use of the exchange criteria of Eq.~\ref{equation:modified-acceptance-criterion}, where path Hamiltonians instead of potential energies are used in computing the probability of accepting a proposed temperature swap between replicas; this amounts to simply using the \emph{total} energies $H(x)$ in the exchange criteria, the sum of potential and kinetic energies, instead of only the potential energies in standard parallel tempering protocols.
A total of 6720 iterations were conducted (67 ns/replica, or 2.7 $\mu$s aggregate), where each iteration consisted of a temperature exchange phase, assignment of new velocities from the Maxwell-Boltzmann distribution at the appropriate replica temperature, and generation of a 10 ps trajectory segment.

\begin{figure}[tb]
\resizebox{\columnwidth}{!}{\includegraphics{figure2}}
\caption{{\bf Normalized fluctuation autocorrelation functions for indicator function on $\alpha_R$ conformation.}
In both plots, single-temperature estimates are shown as points with error bars, and the dynamical reweighting estimate is shown as a solid line with shaded error interval.
Error bars denote 95\% confidence intervals (two standard errors).
\label{figure:alanine-dipeptide-correlation-functions}} 
\end{figure}

To ensure thorough mixing of the replica associations with temperatures, $40^3 =$ 64 000 exchange attempts among randomly-selected pairs of replicas were attempted each parallel tempering iteration.
Because exchanges are attempted among \emph{all} temperatures, and not just neighboring ones, the probability any two temperatures are exchanged during the simulation is diagnostic of whether there is sufficient overlap among temperatures for successful reweighting; in particular, if no swaps occur between specified subsets of temperature, reweighting will produce highly uncertain estimates.
Over the course of the simulation, 31.5\% of proposed neighbor swaps were accepted, 4.5\% of next-neighbor swaps, and 1.1\% of more remote proposed exchanges accepted.
The second-largest eigenvalue of the empirical temperature exchange matrix was less than unity---estimated to be 0.99801---indicating that the set of temperatures did not decompose into two or more subsets with extremely poor exchange properties between them; were there two subsets of temperatures with no transitions in between, this eigenvalue would have been unity.
This indicates that a sufficient number of replicas were chosen to span the range of temperatures with reasonably good exchange among temperatures.

A statistically uncorrelated subset of trajectories generated from the modified parallel tempering simulation were extracted according to the procedure described in Section~\ref{section:modified-parallel-tempering-protocol}.
The statistical inefficiency $g$ of the timeseries $u_n$ was computed, where $u_n$ is defined as the \emph{effective reduced potential} for the replica ensemble, defined in terms of the instantaneous joint configuration of all replica conformations $\mathcal{X} \equiv \{X_1, \ldots, X_K\}$:
\begin{eqnarray}
u_n &\equiv& u(\mathcal{X}_n) = \sum_{k=1}^K \beta_k \mathcal{H}[X_{kn}]
\end{eqnarray}
Here, $X_{kn}$ denotes the trajectory sampled from the replica at temperature $\beta_k$ at iteration $n$ of the parallel tempering simulation, and $\mathcal{H}[X]$ the path Hamiltonian.
As the quantity $u_n$ denotes the log-probability of the overall replica-exchange joint probability distribution function $P(X_1, \ldots, X_K | \beta_1, \ldots, \beta_K)$, and the statistical inefficiency of the timeseries $u_n$, $n = 1, \ldots, N$ roughly corresponds to the number of parallel tempering iterations required to generate an independent set of samples.
After discarding the first 50 iterations to equilibration, the statistical inefficiency was estimated to be 3.1 iterations.

We define the $\alpha_R$ conformation as $\psi \in [-124, 28)$ degrees, based on examination of the one-dimensional potential of mean force (PMF) in the $\psi$ torsion at 300 K (Figure~\ref{figure:alanine-dipeptide-psi-pmf}), which was estimated using the multistate Bennett acceptance ratio (MBAR)~\cite{shirts-chodera:jcp:2008:mbar} and a Gaussian kernel density estimator with a bandwidth of 20 degrees.
This state has a low population at 300 K ($9.6 \pm 0.5$\% as estimated by MBAR), but a relatively long lifetime compared to other conformational transitions in this peptide model.
As an observable, we chose the indicator function for the $\alpha_R$ state:
\begin{eqnarray}
h(x) &=& \begin{cases}
1 & \:\:\psi(x) \in [-124,28) \:\mathrm{degrees} \\
0 & \:\:\mathrm{else}
\end{cases}
\end{eqnarray}
We compute the normalized fluctuation correlation function $C(\tau; \beta)$,
\begin{eqnarray}
C(\tau; \beta) &\equiv& \frac{\expect{\delta h(0) \delta h(\tau)}_\beta}{\expect{\delta h^2}_\beta}
\end{eqnarray}
where we use the shorthand $\delta h(t) \equiv h(x(t)) - \expect{h}_\beta$.
We note that, while the correlation function $C(\tau; \beta)$ describes the relaxation dynamics out of the defined region $\psi \in [-124, 28)$ degrees, this correlation function may contain multiple exponential timescales due to the simplistic definition of the $\alpha_R$ stable state.
Estimation of accurate rate estimates would require calculation of a quantity robust to imperfect definition of dividing surface, such as examination of the reactive flux correlation function~\cite{chandler:jcp:1978:reactive-flux-theory,chandler-berne:1979:jcp:reactive-flux-theory} or the Markov model described in a companion paper~\cite{prinz:2010:temperature-dependent-markov-models}.

Because $C(\tau; \beta)$ is not easily expressed as a single expectation, we write it as a combination of two elementary path expectations.
For a fixed choice of $\tau$, we define the path functionals $\mathcal{A}$ and $\mathcal{B}$ as
\begin{eqnarray}
\mathcal{A}[X] &\equiv& h(x(0)) h(x(\tau)) \nonumber \\
\mathcal{B}[X] &\equiv& h(x(0))
\end{eqnarray}
The correlation function $C(\tau; \beta)$ can then be written as
\begin{eqnarray}
C(\tau; \beta) &=& \frac{\expect{\mathcal{A}}_\beta - \expect{\mathcal{B}}_\beta^2}{\expect{\mathcal{B}}_\beta - \expect{\mathcal{B}}_\beta^2}
\end{eqnarray}
and hence its estimator, in terms of $\hat{\mathcal{A}}(\beta)$ and $\hat{\mathcal{B}}(\beta)$, is
\begin{eqnarray}
\hat{C}(\tau; \beta) &=& \frac{\hat{\mathcal{A}}(\beta) - \hat{\mathcal{B}}(\beta)^2}{\hat{\mathcal{B}}(\beta) - \hat{\mathcal{B}}(\beta)^2}
\end{eqnarray}
where we have used the fact that $(\mathcal{B}[X])^2 = \mathcal{B}[X]$ because $[h(x)]^2 = h(x)$.
The statistical uncertainty in $\hat{C}(\tau; \beta)$ is determined by simple propagation of error.
Suppressing the functional dependence on $\beta$, we have
\begin{eqnarray}
\delta^2 \hat{C} &=& \left[\frac{\partial{\hat{C}}}{\partial{\hat{\mathcal{A}}}}\right]^2 \delta^2 \hat{\mathcal{A}} + \left[\frac{\partial{\hat{C}}}{\partial{\hat{\mathcal{B}}}}\right]^2 \delta^2 \hat{\mathcal{B}} +2 \left[\frac{\partial{\hat{C}}}{\partial{\hat{\mathcal{A}}}}\right] \left[\frac{\partial{\hat{C}}}{\partial{\hat{\mathcal{B}}}}\right] \delta \hat{\mathcal{A}} \delta \hat{\mathcal{B}} \nonumber \\
\end{eqnarray}
where the partial derivatives are given by
\begin{eqnarray}
\frac{\partial{\hat{C}}}{\partial{\hat{\mathcal{A}}}} = \frac{1}{\mathcal{B}(1-\mathcal{B})} \:\: ; \:\: \frac{\partial{\hat{C}}}{\partial{\hat{\mathcal{B}}}} = \frac{\mathcal{A}-1}{(\mathcal{B}-1)^2} - \frac{\mathcal{A}}{\mathcal{B}^2}
\end{eqnarray}
Where the quantities $\delta^2 \hat{\mathcal{A}} \equiv \var \hat{\mathcal{A}}$, $\delta^2 \hat{\mathcal{B}} \equiv \var \hat{\mathcal{B}}$, and $\delta \hat{\mathcal{A}} \delta \hat{\mathcal{B}} \equiv \mathrm{cov} (\hat{\mathcal{A}}, \hat{\mathcal{B}})$ represent the statistical errors (variance or covariance) of the estimators $\hat{\mathcal{A}}$ and $\hat{\mathrm{B}}$.

\begin{figure}[tb]
\resizebox{\columnwidth}{!}{\includegraphics{figure3}}
\caption{{\bf Temperature dependence of the normalized fluctuation autocorrelation function $C(\tau; \beta)$ for $\tau = 10$ ps.}
Single-temperature estimates are shown as points with error bars denoting a 95\% confidence interval; the dynamical reweighting estimate is shown as a solid line with shaded confidence interval.
\label{figure:alanine-dipeptide-correlation-temperature-dependence}} 
\end{figure}

Estimates for $C(\tau; \beta)$ were computed using either the pool of uncorrelated trajectories at each replica temperature, or for a superset of temperatures (including the simulation temperatures) by dynamical reweighting as described in Section~\ref{section:dynamical-reweighting}.
For the former case, statistical error estimates were computed using sample (co)variances; for the latter case, statistical errors were computed using covariances obtained using Eqs.~\ref{equation:uncertainty-estimate} and \ref{equation:correlated-uncertainty}.
Figure~\ref{figure:alanine-dipeptide-correlation-functions} compares correlation functions estimated for a few temperatures using single-temperature data and multiple-temperature data using dynamical reweighting.

To demonstrate how dynamical reweighting produces a smooth estimate throughout the temperature range, even at temperatures not sampled, Figure~\ref{figure:alanine-dipeptide-correlation-temperature-dependence} shows the estimate of $C(\tau;\beta)$ for $\tau = 10$ ps.
Clearly, the reweighted estimate produces a smoother, continuous estimate across all temperatures, with a uniformly reduced statistical uncertainty compared to the estimates computed from individual temperatures.
While some features in the reweighted estimate are likely spurious, the 95\% confidence interval envelope (shaded region) still permits a smooth, well-behaved curve within its boundaries.

%\begin{figure}[tb]
%\resizebox{!}{0.35\columnwidth}{\includegraphics{figure4a}}
%\resizebox{!}{0.35\columnwidth}{\includegraphics{figure4b}}
%\caption{{\bf Bias and variance of single-temperature estimate versus dynamical reweighting estimates.}
%The estimated bias (left) and variance (right) of the estimated normalized fluctuation correlation function estimated by dividing the simulation into 20 blocks of equal length.
%The left half of each color block shows the statistic for single-temperature estimates, while the right half shows the statistic for dynamical reweighting estimates employing data from all temperatures.
%Each horizontal row represents the estimated bias or variance for the correlation function $C(\tau; \beta)$ over $\tau \in [0, 10]$ ps for one temperature; temperatures are ordered from lowest (top) to highest (bottom).
%\label{figure:alanine-dipeptide-bias-variance}} 
%\end{figure}

\begin{figure*}[tb]
\resizebox{!}{0.85\columnwidth}{\includegraphics{figure4a}}
\resizebox{!}{0.85\columnwidth}{\includegraphics{figure4b}}
\caption{{\bf Bias and variance of single-temperature estimate versus dynamical reweighting estimates.}
The estimated bias (left) and variance (right) of the estimated normalized fluctuation correlation function estimated by dividing the simulation into 20 blocks of equal length.
The left half of each color block shows the statistic for single-temperature estimates, while the right half shows the statistic for dynamical reweighting estimates employing data from all temperatures.
Each horizontal row represents the estimated bias or variance for the correlation function $C(\tau; \beta)$ over $\tau \in [0, 10]$ ps for one temperature; temperatures are ordered from lowest (top) to highest (bottom).
\label{figure:alanine-dipeptide-bias-variance}} 
\end{figure*}

Finally, we compare the bias and variance of single-temperature estimates with the dynamical reweighting estimate in a manner independent of the statistical error estimator.
The replica-exchange trajectory was divided into 20 contiguous blocks of equal length, and $C(\tau; \beta)$ estimated independently on each block.
Using the dynamical reweighting estimate over the \emph{entire} dataset as a reference, the bias and variance of the single-temperature and dynamical reweighting estimators were assessed by computing the estimators for each of 20 consecutive subsets of the simulation, and computing the mean deviation or standard deviation of the estimates over the blocks.
The bias, which is defined as the expectation over many independent realizations of the experiment of the deviation of the estimate from the true value, was estimated as the sample mean of the block estimates minus the estimate computed using all data; the variance, which represents the expected squared-deviation of the estimator from its expectation over many independent realizations of the experiment, was estimated by the sample variance over the blocks.
Figure~\ref{figure:alanine-dipeptide-bias-variance} depicts the estimated bias and variance in $\hat{C}(\tau; \beta)$ for all times $\tau \in [0,10]$ and all simulation temperatures $\beta_k$.
The advantage of dynamical reweighting is clear.
Bias is minimized by virtue of the fact that contributions to the dynamical property can occur at any temperature and are incorporated with their appropriate weight, rather than relying upon small-number statistics where the events of interest may occur few or no times within the simulation at the temperature of interest. 
Similarly, inclusion of data from multiple temperatures reduces the overall variance component of the statistical error, resulting in improved estimates at all temperatures.

\begin{figure}[tb]
\resizebox{\columnwidth}{!}{\includegraphics{figure5}}
\caption{{\bf Logarithmic total weight contributed from all trajectories at each sampled temperature to various temperatures.}
The base-10 logarithm of the sum of weights $w_n(\beta)$ from all trajectories sampled from each contributing temperature is shown for each target (reweighted) temperature corresponding to inverse temperature $\beta$.
The contribution from the target temperature is normalized in this plot to have a total weight of unity (zero logarithm).
\label{figure:alanine-dipeptide-log-contributions}} 
\end{figure}

Figure~\ref{figure:alanine-dipeptide-log-contributions} shows the total trajectory weight contributing from each sampled temperature as a function of target temperature $\beta$, where the per-trajectory weight $w_n(\beta)$ from Eq.~\ref{equation:temperature-dependent-trajectory-weights} is summed over all trajectories from the contributing sampled temperature and plotted on a base-10 logarithmic scale, normalized so that the target temperature has a total contributing weight of unity (zero on log scale).
After the target temperature, neighboring temperatures contribute the most trajectory weight to averages at the target temperature, with the contribution rapidly falling off for more distant temperatures.
\color{red}
Note that this does not mean that distant temperatures are unhelpful, as they may still serve to reduce the effective replica correlation time if correlation times at high temperatures are especially short~\cite{jun-s-liu:mcmc}.
\color{black}

Normalizing the total trajectory weight contributed from each temperature to the weight contributed from the target temperature allows us to roughly estimate how much useful data in total we are extracting by making use of reweighting in this system by summing this total relative weight over all temperatures.
This quantity, averaged over all target temperatures, is $\sim 2.4$, suggesting that approximately 2.4 times as much data is provided in this system by making use of reweighting instead of just the data collected from the replicas that visit the target temperature.
Increasing the quantity of data collected by a factor of $2.4$ should decrease the standard error by $\sim \sqrt{2.4} \approx 1.6$, an amount consistent with the difference in confidence interval widths depicted in Figure~\ref{figure:alanine-dipeptide-correlation-functions}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISCUSSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion}

The dynamical reweighting scheme outlined here provides a convenient way to estimate equilibrium dynamical properties from simulations at multiple temperatures, complementing reweighting schemes for estimating static expectations (such as the related multistate Bennett acceptance ratio, or MBAR, method~\cite{shirts-chodera:jcp:2008:mbar}, or the histogram-based WHAM~\cite{kumar:1992a,gallicchio:2005a,chodera:jctc:2007:parallel-tempering-wham}).
Dynamical reweighting provides a way to make use of \emph{all} the data from a parallel or simulated tempering simulation, provided the dynamical model employed is amenable to reweighting  by the unnormalized density ratios $q[X | \beta] / q[X | \beta']$ being both finite and nonzero.
This condition is fulfilled for the common models of dynamics in the canonical ensemble---canonical distribution over Hamiltonian trajectories; Andersen, Langevin, and overdamped Langevin / Brownian stochastic dynamics, reviewed in Section~\ref{section:models-of-dynamics}---but not Nos\'{e}-Hoover or Berendsen.

While not all properties will benefit greatly from the use of reweighting, some dynamical expectations---particularly those that involve large contributions from trajectories that are rare at the temperature of interest, but more plentiful at elevated (or reduced) temperatures---will especially benefit in terms of reduced variance.
For example, though not shown here, transitions along the $\phi$ torsion for terminally-blocked alanine involve conformational states with high free energy which may not even be \emph{sampled} at some temperatures in simulated or parallel tempering simulations of typical length, necessitating the use of reweighting to provide an estimate of rates involving these states (illustrated in detail in a companion paper~\cite{prinz:2010:temperature-dependent-markov-models}).
Even in the case where little or no temperature-dependent enhancement of the phenomena is expected, neighboring temperatures from a simulated or parallel tempering simulation often carry a reasonably large amount of information about the temperature of interest, and their inclusion will further reduce the expected statistical error in the estimated expectation if the data is available anyway.

The variety of models of dynamics within the canonical ensemble presented in Section~\ref{section:models-of-dynamics}, all of which lead to the same static equilibrium distribution of configurations, naturally leads one to ask which model is most appropriate.
In some regimes---such as the low-friction regime of Langevin dynamics and the canonical distribution of Hamiltonian trajectories containing large baths of explicit solvent---these models may give nearly identical expectations, while in other regimes---those where the Andersen collision rate $\nu$ or Langevin collision rate $\gamma$ is large---the expectations for dynamical properties may differ considerably.
Which situation is more physically reasonable will undoubtedly depend on both the system under study and the properties of interest, but a detailed investigation of this is beyond the scope of this work.

The efficiency of reweighting---that is, the degree to which temperatures $\beta' \ne \beta$ contribute to expectations at temperature $\beta$---will depend both on the size of the system and, for stochastic dynamics, the trajectory length.
To understand this behavior, consider a standard parallel tempering simulation using the conventional exchange criteria based on potential energies~\cite{sugita:chem-phys-lett:1999:REMD}.
In order for the ensembles at two temperatures to have good potential energy overlap, and hence good exchange acceptance rates, their temperatures should be spaced roughly the half-width of the potential energy distribution, $\delta U \sim \expect{(U - \expect{U})^2}^{1/2} = T \sqrt{k_B C_v}$,.
The heat capacity for the potential energy, $C_v \equiv (\partial/\partial T) \expect{U}_T$, tells us that the temperature shift $\delta T$ corresponding to this energy shift $\delta U$ is roughly $T \sqrt{k_B / C_v}$.

Analogously, good path Hamiltonian overlap (and good exchange acceptance rates using the modified acceptance criterion in Eq.~\ref{equation:modified-acceptance-criterion}) can be assured with a temperature spacing of roughly $T \sqrt{k_B / \mathcal{C}_v}$, where $\mathcal{C}_v$ is a generalized form of the heat capacity for path Hamiltonians $\mathcal{H}[X]$:
\begin{eqnarray}
\mathcal{C}_v &\equiv& \frac{\partial \expect{\mathcal{H}}_T}{\partial T} = k_B \beta^2 \expect{(\mathcal{H} - \expect{\mathcal{H}})^2}_\beta
\end{eqnarray}
For the canonical distribution of Hamiltonian trajectories, this generalized path heat capacity $\mathcal{C}_v$ is identical to the heat capacity for the total energy, 
%$(\partial / \partial T) \expect{H}_T = (\partial / \partial U) \expect{U}_T + N_d k_B / 2$, which grows linearly with the number of degrees of freedom $N_d$.
$(\partial / \partial T) \expect{H}_T = C_v + N_d k_B / 2$, which grows linearly with the number of degrees of freedom $N_d$.

For Andersen, Langevin, and Brownian dynamics, the generalized path heat capacity $\mathcal{C}_v$ differs from the thermodynamic heat capacity $C_v$.
Because the variance of a sum of uncorrelated random variables is simply the sum of the variances, we can compute this difference analytically.
For Andersen dynamics, $\mathcal{C}_v = C_v + N_c k_B / 2$, where $N_c$ is the number of collisions that occur during realization of the trajectory, where roughly $N_c \sim N_d N_t \nu$, with $\nu$ is the per-step collision probability and $N_t$ is the number of integration timesteps in the trajectory. 
For Langevin dynamics, $\mathcal{C}_v = C_v + N_d (N_t+1) k_B / 2$, and for Brownian dynamics, $\mathcal{C}_v = C_v + N_d N_t k_B / 2$ (where $C_v$ here denotes the heat capacity for the \emph{potential energy} contribution only for Brownian dynamics). 
For constant, nonzero collision rates, the generalized heat capacity is therefore extensive in $N_d N_t$; as the systems grow larger or the trajectories grow longer, the efficiency of reweighting is therefore expected to diminish.

The form of the path Hamiltonians suggests a further analogy with the equilibrium statistical mechanics governing single configurations.
For Hamiltonian dynamics, a path is uniquely identified by its initial phase space point.
For the stochastic forms of dynamics (Andersen, Langevin, Brownian), a trajectory given an initial phase space point can be thought of as an ideal polymer consisting of $N_t$ ``monomers'' that are \emph{replicates} of the $N_p$-particle system, one replicate per trajectory timeslice, where the only interactions are ``bonds'' between corresponding atoms in sequential timeslices.
The interaction energies for these ``bonds'' are harmonic in the noise variables $\xi_{ti}$ corresponding to the atomic displacement between sequential timeslices.

Surprisingly, the path Hamiltonians and generalized path heat capacities are independent of the choice of collision rate $\gamma_i$ for Langevin and Brownian dynamics, except when $\gamma_i$ is identically zero, at which point the particles with zero associated collision rates no longer contribute to the path Hamiltonian or heat capacity since they evolve by deterministic dynamics.
This is analogous to the fact that the heat capacity of an ideal polymer is independent of the spring constant and particle masses.
Only in the case of Andersen dynamics will the collision rate $\nu$ modulate the heat capacity through determining the number of particle collisions observed in a trajectory of fixed length $N_t$.

It should be noted that there is often additional information available that could be incorporated into the reweighting scheme through the use of control variates~\cite{tan:j-am-stat-assoc:2004:extended-bridge-sampling}.
For example, for the model of Hamiltonian trajectories with randomized momenta, the normalization constants $\hat{Z}_i$ that are inferred during the procedure actually contain the temperature-dependent momentum partition function, which can be computed analytically.
Incorporation of this as a constraint, or other integrals of the dynamics that are known exactly, can further reduce the variance in the estimated properties of interest, though the degree to which this may occur will be problem-dependent.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACKNOWLEDGMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgments}

The authors would like to thank Hans C.~Andersen (Stanford University), Nicolae-Viorel Buchete (University College Dublin), Ken Dill (University of California, San Francisco), Peter Eastman (Stanford), Phillip L. Geissler (University of California, Berkeley), Gerhard Hummer (National Institutes of Health), Jes\'{u}s Izaguirre (University of Notre Dame), Antonia S.~J.~S.~Mey (University of Nottingham), David D.~L.~Minh (Argonne National Laboratory), and Jed W.~Pitera (IBM Almaden Research Center) for stimulating discussions on this topic.
JDC gratefully acknowledges support from the HHMI and IBM predoctoral fellowship programs, NIH grant GM34993 through Ken A. Dill (UCSF), and NSF grant for Cyberinfrastructure (NSF CHE-0535616) through Vijay S. Pande (Stanford), and a QB3-Berkeley Distinguished Postdoctoral Fellowship at various points throughout this work.
JHP and FN acknowledge funding from the German Science Foundation (DFG) through the Research Center Matheon, the International Graduate College 710, and grant NO 825/2.
VSP acknowledges NIH grant U54 GM072970.
JDC thanks {\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 associate {\sc PyOpenMM} Python wrappers.
This research was supported in part by the National Science Foundation through TeraGrid resources provided by NCSA Lincoln under grant number TG-MCB100015. 
We specifically acknowledge the assistance of John Estabrook of NCSA Consulting Services in aiding in the use of these resources.

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

%\bibliography{dynamical-reweighting}
%merlin.mbs 2010-03-15 4.21a (PWD, AO, DPC)
%Control: key (0)
%Control: author (8) initials jnrlst
%Control: editor formatted (1) identically to author
%Control: production of article title (-1) disabled
%Control: page (0) single
%Control: year (1) truncated
%Control: production of eprint (0) enabled
\begin{thebibliography}{80}%
\makeatletter
\providecommand \@ifxundefined [1]{%
 \@ifx{#1\undefined}
}%
\providecommand \@ifnum [1]{%
 \ifnum #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \@ifx [1]{%
 \ifx #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \natexlab [1]{#1}%
\providecommand \enquote  [1]{``#1''}%
\providecommand \bibnamefont  [1]{#1}%
\providecommand \bibfnamefont [1]{#1}%
\providecommand \citenamefont [1]{#1}%
\providecommand \href@noop [0]{\@secondoftwo}%
\providecommand \href [0]{\begingroup \@sanitize@url \@href}%
\providecommand \@href[1]{\@@startlink{#1}\@@href}%
\providecommand \@@href[1]{\endgroup#1\@@endlink}%
\providecommand \@sanitize@url [0]{\catcode `\\12\catcode `\$12\catcode
  `\&12\catcode `\#12\catcode `\^12\catcode `\_12\catcode `\%12\relax}%
\providecommand \@@startlink[1]{}%
\providecommand \@@endlink[0]{}%
\providecommand \url  [0]{\begingroup\@sanitize@url \@url }%
\providecommand \@url [1]{\endgroup\@href {#1}{\urlprefix }}%
\providecommand \urlprefix  [0]{URL }%
\providecommand \Eprint [0]{\href }%
\@ifxundefined \urlstyle {%
  \providecommand \doi  [0]{\begingroup \@sanitize@url \@doi}%
  \providecommand \@doi [1]{\endgroup \@@startlink {\doibase
  #1}doi:\discretionary {}{}{}#1\@@endlink }%
}{%
  \providecommand \doi  [0]{doi:\discretionary{}{}{}\begingroup
  \urlstyle{rm}\Url }%
}%
\providecommand \doibase [0]{http://dx.doi.org/}%
\providecommand \Doi [0]{\begingroup \@sanitize@url \@Doi }%
\providecommand \@Doi  [1]{\endgroup\@@startlink{\doibase#1}\@@Doi}%
\providecommand \@@Doi [1]{#1\@@endlink}%
\providecommand \selectlanguage [0]{\@gobble}%
\providecommand \bibinfo  [0]{\@secondoftwo}%
\providecommand \bibfield  [0]{\@secondoftwo}%
\providecommand \translation [1]{[#1]}%
\providecommand \BibitemOpen [0]{}%
\providecommand \bibitemStop [0]{}%
\providecommand \bibitemNoStop [0]{.\EOS\space}%
\providecommand \EOS [0]{\spacefactor3000\relax}%
\providecommand \BibitemShut  [1]{\csname bibitem#1\endcsname}%
%</preamble>
\bibitem [{\citenamefont {Marinari}\ and\ \citenamefont
  {Parisi}(1992)}]{marinari:europhys-lett:1992:simulated-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Marinari}}\ and\ \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Parisi}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal}
  {Europhys. Lett.},\ }\textbf {\bibinfo {volume} {19}},\ \bibinfo {pages}
  {451} (\bibinfo {year} {1992})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Lyubartsev}\ \emph {et~al.}(1992)\citenamefont
  {Lyubartsev}, \citenamefont {Martsinovski}, \citenamefont {Shevkunov},\ and\
  \citenamefont
  {Vorontsov-Velyaminov}}]{lyubartsev:jcp:1992:expanded-ensembles}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~P.}\ \bibnamefont
  {Lyubartsev}}, \bibinfo {author} {\bibfnamefont {A.~A.}\ \bibnamefont
  {Martsinovski}}, \bibinfo {author} {\bibfnamefont {S.~V.}\ \bibnamefont
  {Shevkunov}}, \ and\ \bibinfo {author} {\bibfnamefont {P.~N.}\ \bibnamefont
  {Vorontsov-Velyaminov}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {J. Chem. Phys.},\ }\textbf {\bibinfo {volume} {96}},\ \bibinfo
  {pages} {1776} (\bibinfo {year} {1992})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Mitsutake}\ and\ \citenamefont
  {Okamoto}(2000)}]{mitsutake:2000a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Mitsutake}}\ and\ \bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Okamoto}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Chem.
  Phys. Lett.},\ }\textbf {\bibinfo {volume} {332}},\ \bibinfo {pages} {131}
  (\bibinfo {year} {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Geyer}\ and\ \citenamefont
  {Thompson}(1995)}]{geyer-thompson:j-am-stat-assoc:1995:parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.~J.}\ \bibnamefont
  {Geyer}}\ and\ \bibinfo {author} {\bibfnamefont {E.~A.}\ \bibnamefont
  {Thompson}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Am. Stat. Assoc.},\ }\textbf {\bibinfo {volume} {90}},\ \bibinfo {pages}
  {909} (\bibinfo {year} {1995})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hukushima}\ and\ \citenamefont
  {Nemoto}(1996){\natexlab{a}}}]{hukushima-nemoto:j-phys-soc-jpn:1996:parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Hukushima}}\ and\ \bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Nemoto}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Phys. Soc. Jpn.},\ }\textbf {\bibinfo {volume} {65}},\ \bibinfo {pages}
  {1604} (\bibinfo {year} {1996}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Hansmann}(1997)}]{hansmann:chem-phys-lett:1997:parallel-tempering-monte-carlo}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {U.~H.~E.}\
  \bibnamefont {Hansmann}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {Chem. Phys. Lett.},\ }\textbf {\bibinfo {volume} {281}},\ \bibinfo
  {pages} {140} (\bibinfo {year} {1997})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Sugita}\ and\ \citenamefont
  {Okamoto}(1999)}]{sugita:chem-phys-lett:1999:REMD}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Sugita}}\ and\ \bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Okamoto}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Chem.
  Phys. Lett.},\ }\textbf {\bibinfo {volume} {314}},\ \bibinfo {pages} {141}
  (\bibinfo {year} {1999})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Liu}(2002)}]{jun-s-liu:mcmc}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~S.}\ \bibnamefont
  {Liu}},\ }\href@noop {} {\emph {\bibinfo {title} {{Monte Carlo} strategies in
  scientific computing}}},\ \bibinfo {edition} {2nd}\ ed.\ (\bibinfo
  {publisher} {Springer-Verlag},\ \bibinfo {address} {New York},\ \bibinfo
  {year} {2002})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Trebst}\ \emph
  {et~al.}(2006){\natexlab{a}}\citenamefont {Trebst}, \citenamefont {Troyer},\
  and\ \citenamefont
  {Hansmann}}]{hansmann:jcp:2006:optimized-parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Trebst}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Troyer}}, \
  and\ \bibinfo {author} {\bibfnamefont {U.~H.~E.}\ \bibnamefont {Hansmann}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {124}},\ \bibinfo {pages} {174903} (\bibinfo
  {year} {2006}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Machta}(2009)}]{machta:pre:2009:strengths-and-weaknesses-of-parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Machta}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. E},\ }\textbf {\bibinfo {volume} {80}},\ \bibinfo {pages} {056706}
  (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Rosta}\ and\ \citenamefont
  {Hummer}(2009)}]{hummer:jcp:2009:error-and-efficiency-parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Rosta}}\ and\ \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Hummer}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {2009}},\ \bibinfo {pages} {165102} (\bibinfo
  {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ferrenberg}\ and\ \citenamefont
  {Swendsen}(1989)}]{ferrenberg:1989a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~M.}\ \bibnamefont
  {Ferrenberg}}\ and\ \bibinfo {author} {\bibfnamefont {R.~H.}\ \bibnamefont
  {Swendsen}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal}
  {Phys. Rev. Lett.},\ }\textbf {\bibinfo {volume} {63}},\ \bibinfo {pages}
  {1195} (\bibinfo {year} {1989})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kumar}\ \emph {et~al.}(1992)\citenamefont {Kumar},
  \citenamefont {Bouzida}, \citenamefont {Swendsen}, \citenamefont {Kollman},\
  and\ \citenamefont {Rosenberg}}]{kumar:1992a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Kumar}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Bouzida}},
  \bibinfo {author} {\bibfnamefont {R.~H.}\ \bibnamefont {Swendsen}}, \bibinfo
  {author} {\bibfnamefont {P.~A.}\ \bibnamefont {Kollman}}, \ and\ \bibinfo
  {author} {\bibfnamefont {J.~M.}\ \bibnamefont {Rosenberg}},\ }\href@noop {}
  {\bibfield  {journal} {\bibinfo  {journal} {J. Comput. Chem.},\ }\textbf
  {\bibinfo {volume} {13}},\ \bibinfo {pages} {1011} (\bibinfo {year}
  {1992})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bartels}\ and\ \citenamefont
  {Karplus}(1997)}]{bartels:1997a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.}~\bibnamefont
  {Bartels}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Karplus}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Comput. Chem.},\ }\textbf {\bibinfo {volume} {18}},\ \bibinfo {pages} {1450}
  (\bibinfo {year} {1997})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gallicchio}\ \emph {et~al.}(2005)\citenamefont
  {Gallicchio}, \citenamefont {Andrec}, \citenamefont {Felts},\ and\
  \citenamefont {Levy}}]{gallicchio:2005a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Gallicchio}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Andrec}},
  \bibinfo {author} {\bibfnamefont {A.~K.}\ \bibnamefont {Felts}}, \ and\
  \bibinfo {author} {\bibfnamefont {R.~M.}\ \bibnamefont {Levy}},\ }\href@noop
  {} {\bibfield  {journal} {\bibinfo  {journal} {J. Phys. Chem. B},\ }\textbf
  {\bibinfo {volume} {109}},\ \bibinfo {pages} {6722} (\bibinfo {year}
  {2005})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bartels}(2000)}]{bartels:cpl:2000:wham}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.}~\bibnamefont
  {Bartels}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Chem.
  Phys. Lett.},\ }\textbf {\bibinfo {volume} {331}},\ \bibinfo {pages} {446}
  (\bibinfo {year} {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Shirts}\ and\ \citenamefont
  {Chodera}(2008)}]{shirts-chodera:jcp:2008:mbar}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.~R.}\ \bibnamefont
  {Shirts}}\ and\ \bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont
  {Chodera}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {129}},\ \bibinfo {pages} {124105}
  (\bibinfo {year} {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Chodera}\ \emph
  {et~al.}(2007){\natexlab{a}}\citenamefont {Chodera}, \citenamefont {Swope},
  \citenamefont {Pitera}, \citenamefont {Seok},\ and\ \citenamefont
  {Dill}}]{chodera:jctc:2007:parallel-tempering-wham}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont
  {Chodera}}, \bibinfo {author} {\bibfnamefont {W.~C.}\ \bibnamefont {Swope}},
  \bibinfo {author} {\bibfnamefont {J.~W.}\ \bibnamefont {Pitera}}, \bibinfo
  {author} {\bibfnamefont {C.}~\bibnamefont {Seok}}, \ and\ \bibinfo {author}
  {\bibfnamefont {K.~A.}\ \bibnamefont {Dill}},\ }\href@noop {} {\bibfield
  {journal} {\bibinfo  {journal} {J. Chem. Theor. Comput.},\ }\textbf {\bibinfo
  {volume} {3}},\ \bibinfo {pages} {26} (\bibinfo {year}
  {2007}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Buchete}\ and\ \citenamefont
  {Hummer}(2008){\natexlab{a}}}]{hummer:pre:2008:kinetics-from-replica-exchange}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {N.-V.}\ \bibnamefont
  {Buchete}}\ and\ \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Hummer}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. E},\ }\textbf {\bibinfo {volume} {7}},\ \bibinfo {pages} {030903(R)}
  (\bibinfo {year} {2008}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Skeel}(2009)}]{skeel:siam-j-sci-comput:2008:what-makes-md-work}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~D.}\ \bibnamefont
  {Skeel}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {SIAM J.
  Sci. Comput.},\ }\textbf {\bibinfo {volume} {31}},\ \bibinfo {pages} {1363}
  (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Minh}\ and\ \citenamefont
  {Chodera}(2009)}]{minh-chodera:asymptotic-variance}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~D.~L.}\
  \bibnamefont {Minh}}\ and\ \bibinfo {author} {\bibfnamefont {J.~D.}\
  \bibnamefont {Chodera}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {J. Chem. Phys.},\ }\textbf {\bibinfo {volume} {131}},\ \bibinfo
  {pages} {134110} (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Crooks}(1999)}]{crooks:pre:1999:fluctuation-theorem}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.~E.}\ \bibnamefont
  {Crooks}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. E},\ }\textbf {\bibinfo {volume} {60}},\ \bibinfo {pages} {2721}
  (\bibinfo {year} {1999})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Crooks}(2000)}]{crooks:2000a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.~E.}\ \bibnamefont
  {Crooks}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. E},\ }\textbf {\bibinfo {volume} {61}},\ \bibinfo {pages} {2361}
  (\bibinfo {year} {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Meng}\ and\ \citenamefont
  {Wong}(1996)}]{meng-wong:1996:statistica-sinica:estimating-normalizing-constants}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {X.-L.}\ \bibnamefont
  {Meng}}\ and\ \bibinfo {author} {\bibfnamefont {W.~H.}\ \bibnamefont
  {Wong}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal}
  {Statistica Sinica},\ }\textbf {\bibinfo {volume} {6}},\ \bibinfo {pages}
  {831} (\bibinfo {year} {1996})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kong}\ \emph {et~al.}(2003)\citenamefont {Kong},
  \citenamefont {McCullagh}, \citenamefont {Meng}, \citenamefont {Nicolae},\
  and\ \citenamefont {Tan}}]{kong:j-royal-stat-soc-b:2003:bridge-sampling}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Kong}}, \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {McCullagh}},
  \bibinfo {author} {\bibfnamefont {X.-L.}\ \bibnamefont {Meng}}, \bibinfo
  {author} {\bibfnamefont {D.}~\bibnamefont {Nicolae}}, \ and\ \bibinfo
  {author} {\bibfnamefont {Z.}~\bibnamefont {Tan}},\ }\href@noop {} {\bibfield
  {journal} {\bibinfo  {journal} {J. Royal Stat. Soc. B.},\ }\textbf {\bibinfo
  {volume} {65}},\ \bibinfo {pages} {585} (\bibinfo {year} {2003})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont
  {Tan}(2004)}]{tan:j-am-stat-assoc:2004:extended-bridge-sampling}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Z.}~\bibnamefont
  {Tan}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Am.
  Stat. Assoc.},\ }\textbf {\bibinfo {volume} {99}},\ \bibinfo {pages} {1027}
  (\bibinfo {year} {2004})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bennett}(1977)}]{bennett:1977a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.~H.}\ \bibnamefont
  {Bennett}},\ }\enquote {\bibinfo {title} {Algorithms for chemical
  computations},}\ \ (\bibinfo  {publisher} {Am. Chem. Soc.},\ \bibinfo
  {address} {Washington, D.C.},\ \bibinfo {year} {1977})\ Chap.~\bibinfo
  {chapter} {4}, pp.\ \bibinfo {pages} {63--97}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Chandler}(1998)}]{chandler:1998:cqdcps:tps}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Chandler}},\ }\enquote {\bibinfo {title} {Classical and quantum dynamics in
  condensed phase simulations},}\ \ (\bibinfo  {publisher} {World Scientific},\
  \bibinfo {year} {1998})\ Chap.\ \bibinfo {chapter} {Finding transition
  pathways: throwing ropes over rough mountain passes, in the dark}, pp.\
  \bibinfo {pages} {51--66}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Andersen}(1980)}]{andersen:jcp:1980:andersen-thermostat}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.~C.}\ \bibnamefont
  {Andersen}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {72}},\ \bibinfo {pages} {2384}
  (\bibinfo {year} {1980})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hoover}(1985)}]{hoover:pra:1985:nose-hoover}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.~G.}\ \bibnamefont
  {Hoover}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. A},\ }\textbf {\bibinfo {volume} {31}},\ \bibinfo {pages} {1695}
  (\bibinfo {year} {1985})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Evans}\ and\ \citenamefont
  {Holian}(1985)}]{evans:jcp:1985:nose-hoover}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~J.}\ \bibnamefont
  {Evans}}\ and\ \bibinfo {author} {\bibfnamefont {B.~L.}\ \bibnamefont
  {Holian}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {83}},\ \bibinfo {pages} {4069}
  (\bibinfo {year} {1985})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Berendsen}\ \emph
  {et~al.}(1984){\natexlab{a}}\citenamefont {Berendsen}, \citenamefont
  {Postma}, \citenamefont {van Gunsteren}, \citenamefont {DiNola},\ and\
  \citenamefont {Haak}}]{berendsen:1984:jcp:weak-coupling-thermostat}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.~J.~C.}\
  \bibnamefont {Berendsen}}, \bibinfo {author} {\bibfnamefont {J.~P.~M.}\
  \bibnamefont {Postma}}, \bibinfo {author} {\bibfnamefont {W.~F.}\
  \bibnamefont {van Gunsteren}}, \bibinfo {author} {\bibfnamefont
  {A.}~\bibnamefont {DiNola}}, \ and\ \bibinfo {author} {\bibfnamefont {J.~R.}\
  \bibnamefont {Haak}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {J. Chem. Phys.},\ }\textbf {\bibinfo {volume} {81}},\ \bibinfo
  {pages} {3684} (\bibinfo {year} {1984}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Morishita}(2000)}]{morishita:2000:jcp:berendsen-analysis}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Morishita}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {113}},\ \bibinfo {pages} {2976}
  (\bibinfo {year} {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Cooke}\ and\ \citenamefont
  {Schmidler}(2008)}]{schmidler:2008:jcp:boltzmann-in-replica-exchange}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {B.}~\bibnamefont
  {Cooke}}\ and\ \bibinfo {author} {\bibfnamefont {S.~C.}\ \bibnamefont
  {Schmidler}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {129}},\ \bibinfo {pages} {164112}
  (\bibinfo {year} {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Rosta}\ \emph {et~al.}(2009)\citenamefont {Rosta},
  \citenamefont {Buchete},\ and\ \citenamefont
  {Hummer}}]{hummer:2009:jctc:thermostat-artifacts}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Rosta}}, \bibinfo {author} {\bibfnamefont {N.-V.}\ \bibnamefont {Buchete}},
  \ and\ \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Hummer}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Theor.
  Comput.},\ }\textbf {\bibinfo {volume} {5}},\ \bibinfo {pages} {1393}
  (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Hastings}(1970)}]{hastings:1970:biometrika:metropolis-hastings}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.~K.}\ \bibnamefont
  {Hastings}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal}
  {Biometrika},\ }\textbf {\bibinfo {volume} {57}},\ \bibinfo {pages} {97}
  (\bibinfo {year} {1970})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Duane}\ \emph {et~al.}(1987)\citenamefont {Duane},
  \citenamefont {Kennedy}, \citenamefont {Pendleton},\ and\ \citenamefont
  {Roweth}}]{duane:1987:phys-lett-b:hmc}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Duane}}, \bibinfo {author} {\bibfnamefont {A.~D.}\ \bibnamefont {Kennedy}},
  \bibinfo {author} {\bibfnamefont {B.~J.}\ \bibnamefont {Pendleton}}, \ and\
  \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Roweth}},\ }\href@noop {}
  {\bibfield  {journal} {\bibinfo  {journal} {Phys. Lett. B},\ }\textbf
  {\bibinfo {volume} {195}},\ \bibinfo {pages} {216} (\bibinfo {year}
  {1987})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {E}\ and\ \citenamefont
  {Li}(2008)}]{e:commun-pure-appl-math:2008:andersen-thermostat}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {E}}\
  and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Li}},\ }\href@noop
  {} {\bibfield  {journal} {\bibinfo  {journal} {Commun. Pure Appl. Math},\
  }\textbf {\bibinfo {volume} {61}},\ \bibinfo {pages} {0096} (\bibinfo {year}
  {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Juraszek}\ and\ \citenamefont
  {Bolhuis}(2006)}]{bolhuis:pnas:2006:trp-cage}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Juraszek}}\ and\ \bibinfo {author} {\bibfnamefont {P.~G.}\ \bibnamefont
  {Bolhuis}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Proc.
  Nat. Acad. Sci. USA},\ }\textbf {\bibinfo {volume} {103}},\ \bibinfo {pages}
  {15859} (\bibinfo {year} {2006})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Andersen}(1983)}]{andersen:1983:j-comput-phys:rattle}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.~C.}\ \bibnamefont
  {Andersen}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Comput. Phys.},\ }\textbf {\bibinfo {volume} {52}},\ \bibinfo {pages} {24}
  (\bibinfo {year} {1983})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Izaguirre}\ \emph {et~al.}(2010)\citenamefont
  {Izaguirre}, \citenamefont {Sweet},\ and\ \citenamefont
  {Pande}}]{izaguirre:2010:psb:normal-mode-langevin}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~A.}\ \bibnamefont
  {Izaguirre}}, \bibinfo {author} {\bibfnamefont {C.~R.}\ \bibnamefont
  {Sweet}}, \ and\ \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont
  {Pande}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Pac.
  Symp. Biocomp.},\ }\textbf {\bibinfo {volume} {15}},\ \bibinfo {pages} {240}
  (\bibinfo {year} {2010})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ermak}\ and\ \citenamefont
  {Yeh}(1974)}]{ermak-yeh:cpl:1974:brownian-dynamics}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~L.}\ \bibnamefont
  {Ermak}}\ and\ \bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont {Yeh}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Chem. Phys.
  Lett.},\ }\textbf {\bibinfo {volume} {24}},\ \bibinfo {pages} {243} (\bibinfo
  {year} {1974})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ermak}(1975)}]{ermak:jcp:1975:brownian-dynamics}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~L.}\ \bibnamefont
  {Ermak}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {62}},\ \bibinfo {pages} {4189}
  (\bibinfo {year} {1975})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Buchete}\ and\ \citenamefont
  {Hummer}(2008){\natexlab{b}}}]{hummer:jpcb:2008:coarse-master-equations}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {N.-V.}\ \bibnamefont
  {Buchete}}\ and\ \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Hummer}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Phys. Chem. B},\ }\textbf {\bibinfo {volume} {112}},\ \bibinfo {pages} {6057}
  (\bibinfo {year} {2008}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Berendsen}\ \emph
  {et~al.}(1984){\natexlab{b}}\citenamefont {Berendsen}, \citenamefont
  {Postma}, \citenamefont {van Gunsteren}, \citenamefont {DiNola},\ and\
  \citenamefont {Haak}}]{berendsen:1984a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.~J.~C.}\
  \bibnamefont {Berendsen}}, \bibinfo {author} {\bibfnamefont {J.~P.~M.}\
  \bibnamefont {Postma}}, \bibinfo {author} {\bibfnamefont {W.~F.}\
  \bibnamefont {van Gunsteren}}, \bibinfo {author} {\bibfnamefont
  {A.}~\bibnamefont {DiNola}}, \ and\ \bibinfo {author} {\bibfnamefont {J.~R.}\
  \bibnamefont {Haak}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {J. Chem. Phys.},\ }\textbf {\bibinfo {volume} {81}},\ \bibinfo
  {pages} {3684} (\bibinfo {year} {1984}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Sindhikara}\ \emph {et~al.}(2008)\citenamefont
  {Sindhikara}, \citenamefont {Meng},\ and\ \citenamefont
  {Roitberg}}]{roitberg:jcp:2008:exchange-frequency}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Sindhikara}}, \bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont {Meng}}, \
  and\ \bibinfo {author} {\bibfnamefont {A.~E.}\ \bibnamefont {Roitberg}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {128}},\ \bibinfo {pages} {024103} (\bibinfo
  {year} {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Vlugt}\ and\ \citenamefont
  {Smit}(2001)}]{vlugt-smit:2001:phys-chem-comm:parallel-tempering-tps}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~J.~H.}\
  \bibnamefont {Vlugt}}\ and\ \bibinfo {author} {\bibfnamefont
  {B.}~\bibnamefont {Smit}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Chem. Comm.},\ }\textbf {\bibinfo {volume} {2}},\ \bibinfo
  {pages} {1} (\bibinfo {year} {2001})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bolhuis}\ \emph {et~al.}(2002)\citenamefont
  {Bolhuis}, \citenamefont {Chandler}, \citenamefont {Dellago},\ and\
  \citenamefont
  {Geissler}}]{bolhuis:annu-rev-phys-chem:2002:transition-path-sampling}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.~G.}\ \bibnamefont
  {Bolhuis}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Chandler}},
  \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Dellago}}, \ and\
  \bibinfo {author} {\bibfnamefont {P.~L.}\ \bibnamefont {Geissler}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Annu. Rev. Phys.
  Chem.},\ }\textbf {\bibinfo {volume} {53}},\ \bibinfo {pages} {291} (\bibinfo
  {year} {2002})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Geyer}(191)}]{geyer:1991:parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.~J.}\ \bibnamefont
  {Geyer}},\ }in\ \href@noop {} {\emph {\bibinfo {booktitle} {Computing Science
  and Statistics: The 23rd symposium on the interface}}},\ \bibinfo {editor}
  {edited by\ \bibinfo {editor} {\bibfnamefont {E.}~\bibnamefont {Keramigas}}}\
  (\bibinfo  {publisher} {Interface Foundation},\ \bibinfo {year} {191})\ pp.\
  \bibinfo {pages} {156--163}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hukushima}\ and\ \citenamefont
  {Nemoto}(1996){\natexlab{b}}}]{hukushima-nemoto:1996:j-phys-soc-japan:parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Hukushima}}\ and\ \bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Nemoto}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Phys. Soc. Japan},\ }\textbf {\bibinfo {volume} {65}},\ \bibinfo {pages}
  {1604} (\bibinfo {year} {1996}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Pitera}\ \emph {et~al.}(2006)\citenamefont {Pitera},
  \citenamefont {Haque},\ and\ \citenamefont {Swope}}]{pitera:2006a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~W.}\ \bibnamefont
  {Pitera}}, \bibinfo {author} {\bibfnamefont {I.}~\bibnamefont {Haque}}, \
  and\ \bibinfo {author} {\bibfnamefont {W.~C.}\ \bibnamefont {Swope}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {124}},\ \bibinfo {pages} {141102} (\bibinfo
  {year} {2006})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Swope}\ \emph
  {et~al.}(1982){\natexlab{a}}\citenamefont {Swope}, \citenamefont {Andersen},
  \citenamefont {Berens},\ and\ \citenamefont {Wilson}}]{swope:1982a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.~C.}\ \bibnamefont
  {Swope}}, \bibinfo {author} {\bibfnamefont {H.~C.}\ \bibnamefont {Andersen}},
  \bibinfo {author} {\bibfnamefont {P.~H.}\ \bibnamefont {Berens}}, \ and\
  \bibinfo {author} {\bibfnamefont {K.~R.}\ \bibnamefont {Wilson}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {76}},\ \bibinfo {pages} {637} (\bibinfo {year}
  {1982}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Flyvbjerg}\ and\ \citenamefont
  {Petersen}(1989)}]{flyvbjerg:1989a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
  {Flyvbjerg}}\ and\ \bibinfo {author} {\bibfnamefont {H.~G.}\ \bibnamefont
  {Petersen}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {91}},\ \bibinfo {pages} {461}
  (\bibinfo {year} {1989})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Janke}(2002)}]{janke:2002a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
  {Janke}},\ }in\ \href@noop {} {\emph {\bibinfo {booktitle} {Quantum
  Simulations of Complex Many-Body Systems: From Theory to Algorithms}}},\
  Vol.~\bibinfo {volume} {10},\ \bibinfo {editor} {edited by\ \bibinfo {editor}
  {\bibfnamefont {J.}~\bibnamefont {Grotendorst}}, \bibinfo {editor}
  {\bibfnamefont {D.}~\bibnamefont {Marx}}, \ and\ \bibinfo {editor}
  {\bibfnamefont {A.}~\bibnamefont {Murmatsu}}}\ (\bibinfo  {publisher} {John
  von Neumann Institute for Computing},\ \bibinfo {year} {2002})\ pp.\ \bibinfo
  {pages} {423--445}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Chodera}\ \emph
  {et~al.}(2007){\natexlab{b}}\citenamefont {Chodera}, \citenamefont {Swope},
  \citenamefont {Pitera}, \citenamefont {Seok},\ and\ \citenamefont
  {Dill}}]{chodera:jctc:2007}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont
  {Chodera}}, \bibinfo {author} {\bibfnamefont {W.~C.}\ \bibnamefont {Swope}},
  \bibinfo {author} {\bibfnamefont {J.~W.}\ \bibnamefont {Pitera}}, \bibinfo
  {author} {\bibfnamefont {C.}~\bibnamefont {Seok}}, \ and\ \bibinfo {author}
  {\bibfnamefont {K.~A.}\ \bibnamefont {Dill}},\ }\href@noop {} {\bibfield
  {journal} {\bibinfo  {journal} {J. Chem. Theor. Comput.},\ }\textbf {\bibinfo
  {volume} {3}},\ \bibinfo {pages} {26} (\bibinfo {year}
  {2007}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Kofke}(2002)}]{kofke:2002:jcp:acceptance-probability}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont
  {Kofke}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {117}},\ \bibinfo {pages} {6911}
  (\bibinfo {year} {2002})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Katzgraber}\ \emph {et~al.}(2006)\citenamefont
  {Katzgraber}, \citenamefont {Trebst}, \citenamefont {Huse},\ and\
  \citenamefont
  {Troyer}}]{katzberger-trebst-huse-troyer:j-stat-mech:2006:feedback-optimized-parallel-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.~G.}\ \bibnamefont
  {Katzgraber}}, \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont {Trebst}},
  \bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont {Huse}}, \ and\
  \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Troyer}},\ }\href@noop {}
  {\bibfield  {journal} {\bibinfo  {journal} {J. Stat. Mech.},\ }\textbf
  {\bibinfo {volume} {2006}},\ \bibinfo {pages} {03018} (\bibinfo {year}
  {2006})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Trebst}\ \emph
  {et~al.}(2006){\natexlab{b}}\citenamefont {Trebst}, \citenamefont {Troyer},\
  and\ \citenamefont
  {Hansmann}}]{trebst-troyer-hansmann:jcp:2006:optimized-replica-selection}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Trebst}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Troyer}}, \
  and\ \bibinfo {author} {\bibfnamefont {U.~H.~E.}\ \bibnamefont {Hansmann}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {124}},\ \bibinfo {pages} {174903} (\bibinfo
  {year} {2006}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gront}\ and\ \citenamefont
  {Kolinski}(2007)}]{gront-kolinski:j-phys-condens-matter:2007:optimized-replica-selection}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Gront}}\ and\ \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Kolinski}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Phys.: Condens. Matter},\ }\textbf {\bibinfo {volume} {19}},\ \bibinfo
  {pages} {036225} (\bibinfo {year} {2007})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Park}\ and\ \citenamefont
  {Pande}(2007)}]{park-pande:pre:2007:choosing-weights-simulated-tempering}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Park}}\ and\ \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont
  {Pande}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. E},\ }\textbf {\bibinfo {volume} {76}},\ \bibinfo {pages} {016703}
  (\bibinfo {year} {2007})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Chandler}(1978)}]{chandler:jcp:1978:reactive-flux-theory}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Chandler}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {68}},\ \bibinfo {pages} {2959}
  (\bibinfo {year} {1978})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Jr.}\ \emph {et~al.}(1979)\citenamefont {Jr.},
  \citenamefont {Chandler},\ and\ \citenamefont
  {Berne}}]{chandler-berne:1979:jcp:reactive-flux-theory}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~A.~M.}\
  \bibnamefont {Jr.}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Chandler}}, \ and\ \bibinfo {author} {\bibfnamefont {B.~J.}\ \bibnamefont
  {Berne}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.},\ }\textbf {\bibinfo {volume} {70}},\ \bibinfo {pages} {4056}
  (\bibinfo {year} {1979})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Prinz}\ \emph {et~al.}(2010)\citenamefont {Prinz},
  \citenamefont {Chodera}, \citenamefont {Pande}, \citenamefont {Swope},
  \citenamefont {Smith},\ and\ \citenamefont
  {No\'{e}}}]{prinz:2010:temperature-dependent-markov-models}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.-H.}\ \bibnamefont
  {Prinz}}, \bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont {Chodera}},
  \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont {Pande}}, \bibinfo
  {author} {\bibfnamefont {W.~C.}\ \bibnamefont {Swope}}, \bibinfo {author}
  {\bibfnamefont {J.~C.}\ \bibnamefont {Smith}}, \ and\ \bibinfo {author}
  {\bibfnamefont {F.}~\bibnamefont {No\'{e}}},\ }\href@noop {} {,\ \bibinfo
  {pages} {Submitted to JCP.} (\bibinfo {year} {2010})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Zhang}\ \emph {et~al.}(2009)\citenamefont {Zhang},
  \citenamefont {Hou}, \citenamefont {Schafmeister}, \citenamefont {Ross},\
  and\ \citenamefont {Case}}]{leap-ambertools-12}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
  {Zhang}}, \bibinfo {author} {\bibfnamefont {T.}~\bibnamefont {Hou}}, \bibinfo
  {author} {\bibfnamefont {C.}~\bibnamefont {Schafmeister}}, \bibinfo {author}
  {\bibfnamefont {W.~S.}\ \bibnamefont {Ross}}, \ and\ \bibinfo {author}
  {\bibfnamefont {D.~A.}\ \bibnamefont {Case}},\ }\href
  {http://ambermd.org/#AmberTools} {\enquote {\bibinfo {title} {Leap from
  ambertools 1.2},}\ } (\bibinfo {year} {2009})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kollman}(1996)}]{kollman:acc-chem-res:1996:amber96}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.~A.}\ \bibnamefont
  {Kollman}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Acc.
  Chem. Res.},\ }\textbf {\bibinfo {volume} {29}},\ \bibinfo {pages} {461}
  (\bibinfo {year} {1996})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Jorgensen}\ \emph {et~al.}(1983)\citenamefont
  {Jorgensen}, \citenamefont {Chandrasekhar}, \citenamefont {Madura},
  \citenamefont {Impey},\ and\ \citenamefont {Klein}}]{jorgensen:1983a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.~L.}\ \bibnamefont
  {Jorgensen}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Chandrasekhar}}, \bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont
  {Madura}}, \bibinfo {author} {\bibfnamefont {R.~W.}\ \bibnamefont {Impey}}, \
  and\ \bibinfo {author} {\bibfnamefont {M.~L.}\ \bibnamefont {Klein}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {79}},\ \bibinfo {pages} {926} (\bibinfo {year}
  {1983})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Liu}\ and\ \citenamefont
  {Nocedal}(1989)}]{liu-nocedal:1989:math-prog:L-BFGS}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~C.}\ \bibnamefont
  {Liu}}\ and\ \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Nocedal}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Mathematical
  Programming},\ }\textbf {\bibinfo {volume} {45}},\ \bibinfo {pages} {503}
  (\bibinfo {year} {1989})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Nocedal}(1980)}]{nocedal:1980:math-comput:L-BFGS}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Nocedal}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal}
  {Mathematics of Computation},\ }\textbf {\bibinfo {volume} {35}},\ \bibinfo
  {pages} {773} (\bibinfo {year} {1980})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Beauchamp}\ \emph {et~al.}(2010)\citenamefont
  {Beauchamp}, \citenamefont {Bruns}, \citenamefont {Eastman}, \citenamefont
  {Friedrichs}, \citenamefont {Ku}, \citenamefont {Pande}, \citenamefont
  {Radmer},\ and\ \citenamefont {Sherman}}]{openmm-theory-guide}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Beauchamp}}, \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Bruns}},
  \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Eastman}}, \bibinfo
  {author} {\bibfnamefont {M.}~\bibnamefont {Friedrichs}}, \bibinfo {author}
  {\bibfnamefont {J.~P.}\ \bibnamefont {Ku}}, \bibinfo {author} {\bibfnamefont
  {V.}~\bibnamefont {Pande}}, \bibinfo {author} {\bibfnamefont
  {R.}~\bibnamefont {Radmer}}, \ and\ \bibinfo {author} {\bibfnamefont
  {M.}~\bibnamefont {Sherman}},\ }\href {http://simtk.org/home/openmm}
  {\enquote {\bibinfo {title} {Openmm user manual and theory guide, release
  2.0},}\ } (\bibinfo {year} {2010})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Eastman}\ and\ \citenamefont
  {Pande}(2010){\natexlab{a}}}]{eastman:2010:j-comput-chem:ccma}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Eastman}}\ and\ \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont
  {Pande}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Theor. Comput.},\ }\textbf {\bibinfo {volume} {6}},\ \bibinfo {pages}
  {434} (\bibinfo {year} {2010}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Miyamoto}\ and\ \citenamefont
  {Kollman}(1992)}]{kollman:1992:j-comput-chem:settle}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Miyamoto}}\ and\ \bibinfo {author} {\bibfnamefont {P.~A.}\ \bibnamefont
  {Kollman}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Comput. Chem.},\ }\textbf {\bibinfo {volume} {13}},\ \bibinfo {pages} {952}
  (\bibinfo {year} {1992})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Darden}\ \emph {et~al.}(1993)\citenamefont {Darden},
  \citenamefont {York},\ and\ \citenamefont {Pedersen}}]{darden:1993a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~A.}\ \bibnamefont
  {Darden}}, \bibinfo {author} {\bibfnamefont {D.~M.}\ \bibnamefont {York}}, \
  and\ \bibinfo {author} {\bibfnamefont {L.~G.}\ \bibnamefont {Pedersen}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {98}},\ \bibinfo {pages} {10089} (\bibinfo {year}
  {1993})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Allen}\ and\ \citenamefont
  {Tildesley}(1991)}]{allen-tildesley}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.~P.}\ \bibnamefont
  {Allen}}\ and\ \bibinfo {author} {\bibfnamefont {D.~J.}\ \bibnamefont
  {Tildesley}},\ }\href@noop {} {\emph {\bibinfo {title} {Computer simulation
  of liquids}}}\ (\bibinfo  {publisher} {Clarendon Press},\ \bibinfo {address}
  {Oxford},\ \bibinfo {year} {1991})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Friedrichs}\ \emph {et~al.}(2009)\citenamefont
  {Friedrichs}, \citenamefont {Eastman}, \citenamefont {Vaidyanathan},
  \citenamefont {Houston}, \citenamefont {LeGrand}, \citenamefont {Beberg},
  \citenamefont {Ensign}, \citenamefont {Bruns},\ and\ \citenamefont
  {Pande}}]{friedrichs:2009:j-comput-chem:openmm}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.~S.}\ \bibnamefont
  {Friedrichs}}, \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Eastman}},
  \bibinfo {author} {\bibfnamefont {V.}~\bibnamefont {Vaidyanathan}}, \bibinfo
  {author} {\bibfnamefont {M.}~\bibnamefont {Houston}}, \bibinfo {author}
  {\bibfnamefont {S.}~\bibnamefont {LeGrand}}, \bibinfo {author} {\bibfnamefont
  {A.~L.}\ \bibnamefont {Beberg}}, \bibinfo {author} {\bibfnamefont {D.~L.}\
  \bibnamefont {Ensign}}, \bibinfo {author} {\bibfnamefont {C.~M.}\
  \bibnamefont {Bruns}}, \ and\ \bibinfo {author} {\bibfnamefont {V.~S.}\
  \bibnamefont {Pande}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo
  {journal} {J. Comput. Chem.},\ }\textbf {\bibinfo {volume} {30}},\ \bibinfo
  {pages} {864} (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Eastman}\ and\ \citenamefont
  {Pande}(2010){\natexlab{b}}}]{eastman:2010:comp-sci-eng:openmm}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Eastman}}\ and\ \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont
  {Pande}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal}
  {Computing in Science and Engineering},\ }\textbf {\bibinfo {volume} {12}},\
  \bibinfo {pages} {34} (\bibinfo {year} {2010}{\natexlab{b}})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Eastman}\ and\ \citenamefont
  {Pande}(2010){\natexlab{c}}}]{eastman:2010:j-comput-chem:openmm}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Eastman}}\ and\ \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont
  {Pande}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Comput. Chem.},\ }\textbf {\bibinfo {volume} {31}},\ \bibinfo {pages} {1268}
  (\bibinfo {year} {2010}{\natexlab{c}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bruns}\ \emph {et~al.}(2010)\citenamefont {Bruns},
  \citenamefont {Radmer}, \citenamefont {Chodera},\ and\ \citenamefont
  {Pande}}]{pyopenmm}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.~M.}\ \bibnamefont
  {Bruns}}, \bibinfo {author} {\bibfnamefont {R.~A.}\ \bibnamefont {Radmer}},
  \bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont {Chodera}}, \ and\
  \bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont {Pande}},\ }\href
  {http://simtk.org/home/pyopenmm} {\enquote {\bibinfo {title} {Pyopenmm},}\ }
  (\bibinfo {year} {2010})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Swope}\ \emph
  {et~al.}(1982){\natexlab{b}}\citenamefont {Swope}, \citenamefont {Andersen},
  \citenamefont {Berens},\ and\ \citenamefont
  {Wilson}}]{swope:1982:jcp:velocity-verlet}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.~C.}\ \bibnamefont
  {Swope}}, \bibinfo {author} {\bibfnamefont {H.~C.}\ \bibnamefont {Andersen}},
  \bibinfo {author} {\bibfnamefont {P.~H.}\ \bibnamefont {Berens}}, \ and\
  \bibinfo {author} {\bibfnamefont {K.~R.}\ \bibnamefont {Wilson}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. Chem. Phys.},\
  }\textbf {\bibinfo {volume} {76}},\ \bibinfo {pages} {637} (\bibinfo {year}
  {1982}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Verlet}(1967)}]{verlet:1967:phys-rev:verlet-integrator-1}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
  {Verlet}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev.},\ }\textbf {\bibinfo {volume} {159}},\ \bibinfo {pages} {98} (\bibinfo
  {year} {1967})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont
  {Verlet}(1968)}]{verlet:1967:phys-rev:verlet-integrator-2}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
  {Verlet}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev.},\ }\textbf {\bibinfo {volume} {165}},\ \bibinfo {pages} {201} (\bibinfo
  {year} {1968})}\BibitemShut {NoStop}%
\end{thebibliography}%

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

\appendix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIMULATION DETAILS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Additional simulation details for alanine dipeptide}
\label{appendix:simulation-details}

Using the {\tt LEaP} program from the AmberTools 1.2 molecular mechanics package~\cite{leap-ambertools-12}, a terminally-blocked alanine peptide (sequence ACE-ALA-NME, see Figure~\ref{figure:alanine-dipeptide-psi-pmf}) was generated in the extended conformation, with peptide force field parameters taken from the AMBER {\tt parm96} parameter set \cite{kollman:acc-chem-res:1996:amber96}.
The system was subsequently solvated with 749 TIP3P water molecules \cite{jorgensen:1983a} in a cubic simulation box with dimensions chosen to ensure all box boundaries were at least 9 \AA\ from any atom of the extended peptide.
The system was subjected to energy minimization using L-BFGS~\cite{liu-nocedal:1989:math-prog:L-BFGS,nocedal:1980:math-comput:L-BFGS} to reduce the root-mean-square force to less than 1 kJ/mol/nm, and then equilibrated for 1 ns with a leapfrog Langevin integrator~\cite{izaguirre:2010:psb:normal-mode-langevin} at a control temperature of 300 K, using a 2 fs timestep and collision rate of 5/ps.
During equilibration, an isotropic Monte Carlo barostat with a control pressure of 1 atm was applied, with volume moves attempted every 25 timesteps, and the range of volume change proposals automatically adjusted to give an approximate acceptance probability of approximately 50\%~\cite{openmm-theory-guide}. 
All bonds to hydrogen, and hydrogen-hydrogen distances in waters, were constrained by the CCMA~\cite{eastman:2010:j-comput-chem:ccma} and SETTLE~\cite{kollman:1992:j-comput-chem:settle} algorithms, as appropriate.
The particle-mesh Ewald (PME) method~\cite{darden:1993a} was used to treat electrostatics, using a real-space cutoff of 9 \AA.
PME parameters were automatically selected by an algorithm that attempts to find the most efficient set of parameters for which an error bound is less than a specified error tolerance~\cite{openmm-theory-guide}, which was set to the default of $5\cdot 10^{-4}$.
Lennard-Jones interactions were truncated at 9 \AA\ without switching, and a homogeneous analytical long-range dispersion correction~\cite{allen-tildesley} was employed to compensate for dispersion interactions outside the cutoff.
The resulting system had a box volume of 23 073.7 \AA$^3$, and the box volume was fixed in the subsequent parallel tempering simulation.

A custom Python code making use of the GPU-accelerated {\sc OpenMM} package~\cite{friedrichs:2009:j-comput-chem:openmm,eastman:2010:comp-sci-eng:openmm,eastman:2010:j-comput-chem:openmm} and the {\sc PyOpenMM} Python wrapper~\cite{pyopenmm} was used to conduct the simulations.
Because {\sc OpenMM} lacks a velocity Verlet integrator, a hybrid velocity Verlet~\cite{swope:1982:jcp:velocity-verlet} / leapfrog~\cite{verlet:1967:phys-rev:verlet-integrator-1,verlet:1967:phys-rev:verlet-integrator-2} integration scheme was used for integration of the equations of motion, implemented as follows.
At the beginning of the dynamical propagation phase of parallel tempering iteration, on-step velocities were first generated from the Maxwell-Boltzmann distribution at the current replica temperature.
Components of the velocity along constrained bonds (all bonds to hydrogen and the hydrogen-hydrogen bonds of TIP3P water) were removed using the RATTLE algorithm~\cite{andersen:1983:j-comput-phys:rattle} with $10^{-16}$ relative tolerance, and a 1 fs ``half-kick backwards'' applied to the modified velocities using the force evaluated at the current positions.
The {\sc OpenMM} Verlet integrator was then used to evolve the positions by 0.5 ps (250 steps with a 2 fs timestep) on the GPU, and a 1 fs ``half-kick forwards'' was applied to the velocities to synchronize them with the positions after again applying the RATTLE algorithm.
During leapfrog integration, the CCMA and SETTLE algorithms were used to constrain bonds, as appropriate.



\end{document}