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

\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{{fig/}}
\bibliographystyle{apsrev}
%

%% two-column format
%\documentclass[9pt,draft]{revtex4}   % to mimic appearance in journal

%% We will be using PDF figures and generating PDF output.
%\pdfoutput=1
%\pdfcompresslevel=9
%\usepackage[pdftex]{graphicx}
%\DeclareGraphicsExtensions{.pdf}

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

%% Import useful packages.
%\usepackage{times}                                 % nice fonts
%\usepackage{amsmath}                           % for 'case'
%\usepackage{amsthm}                             % for theorems and proofs
%\usepackage{mathrsfs}                            % script math font
%\usepackage{url}                                           % allow use of URLs
%\usepackage{latexsym}                            % some math symbols, like \Box
%\usepackage{subfigure}                           % allow use of subfigures

%\newcommand{\bfm}[1]{{\bf #1}}              % proper boldface for math environment
\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
\newcommand{\x}{\bfm{x}}
\newcommand{\y}{\bfm{y}}
\newcommand{\f}{\bfm{f}}
\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{\Ww}[1] {{\mathrm w}\left( #1 \right)}
\newcommand{\Ht}{\bfm{\tilde{H}}}
\newcommand{\gt}{\bfm{\tilde{g}}}
\newcommand{\Dt}{\bfm{\tilde{D}}}
\newcommand{\pt}{\bfm{\tilde{p}}}
\newcommand{\Nt}{\bfm{\tilde{N}}}
\newcommand{\ft}{\bfm{\tilde{f}}}

\DeclareSymbolFont{AMSb}{U}{msb}{m}{n}
\DeclareMathSymbol{\N}{\mathbin}{AMSb}{"4E}
\DeclareMathSymbol{\Z}{\mathbin}{AMSb}{"5A}
\DeclareMathSymbol{\R}{\mathbin}{AMSb}{"52}
\DeclareMathSymbol{\Q}{\mathbin}{AMSb}{"51}
\DeclareMathSymbol{\I}{\mathbin}{AMSb}{"49}
\DeclareMathSymbol{\C}{\mathbin}{AMSb}{"43}


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

%% TITLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Applying statistically optimal methods for analysis of multiple equilibrium simulations}


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

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

\date{\today}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}

We combine the idea of Wang-Landau simulations with the concept of
waste recycling to be able to converge simulations faster.  Fun ensues!

\emph{Keywords: Monte Carlo; Wang-Landau; Multistate Bennett acceptance ratio}

\end{abstract}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTRODUCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}


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

Before describing our suggested algorithmic modifications
(\emph{Algorithms}, Section~\ref{section:algorithms}), we first
present some theoretical tools we will use to analyze expanded
ensemble and replica exchange simulations in the context of Gibbs
sampling.

\subsection{Thermodynamic states and thermodynamic ensembles}

To be as general as possible, we describe the expanded ensemble and
replica exchange algorithms as sampling a mixture of $K$ thermodynamic
states.  Here, a \emph{thermodynamic state} is parameterized by a
vector of time-independent thermodynamic parameters $\lambda$.  For
notational convenience and to make what follows general, we define the
\emph{reduced potential}~\cite{shirts-chodera:jcp:2008:mbar} $u(x)$ of
a physical system,
\begin{eqnarray}
u(x) &=& \beta \left[ H(x) + p V(x) + \sum_i\mu_i n_i(x) + \cdots \right], \label{equation:reduced-potential}
\end{eqnarray}
corresponding to its thermodynamic state, where $x$ denotes the
configuration of the system specifying any physical variables allowed
to change, including the volume $V(x)$ (in the case of a constant
pressure ensemble) and $n_i(x)$ the number of molecules of each of
$i=1,\ldots,M$ components of the system, in the case of a (semi)grand
ensemble.

The reduced potential is a function of the Hamiltonian $H$, the
inverse temperature $\beta = (k_B T)^{-1}$, the pressure $p$, and the
vector of chemical potentials for each of $M$ components $\mu_i$.
Other thermodynamic parameters and their conjugate coordinates can be
included in a similar manner, or some of these can be omitted, as
required by the physics of the system.  We denote the set of all
thermodynamic parameters by $\lambda \equiv \{\beta, H, p, \vec{\mu},
\ldots\}$.

We next denote a configuration of the molecular system by $x \in
\Omega$, where $\Omega$ is allowed configuration space, which may be
continuous or discrete.  A choice of thermodynamic state gives rise to
set of configurations of the system that are sampled by a given
time-independent probability distribution at equilibrium.  So each $x$
will have associated unnormalized probability density $q(x)$, which is
a function of $\lambda$, where $q(x) > 0$ for all $x \in \Omega$. If
we define the normalization constant, or \emph{partition function},
$Z$ as:
\begin{eqnarray}
Z &\equiv& \int_\Omega dx \, q(x) 
\end{eqnarray}
we can define a normalized probability density
\begin{eqnarray}
\pi(x) &=& Z^{-1} \, q(x).
\end{eqnarray}

A physical system in equilibrium with its environment obeying
classical statistical mechanics will sample configurations distributed
according to the Boltzmann distribution,
\begin{eqnarray}
q(x) &\equiv& e^{-u(x)} .
\end{eqnarray}

In this paper, we consider a set of $K$ thermodynamic states defined
by their thermodynamic parameter vectors, $\lambda_k \equiv \{\beta_k,
H_k, p_k, \vec{\mu}_k, \ldots\}$, with $k = 1,\ldots,K$, where $H_k$
denotes any modifications of the Hamiltonian $H$ as a function of $k$,
including biasing potentials.  Each new choice of $k$ gives rise to a
reduced potential $u_k$, unnormalized and normalized probability
distributions $q_k(x)$ and $\pi(x,k)$, and a partition function $Z_k$.
To ensure that any configuration $x$ has finite, nonzero density in
all $K$ thermodynamic states, we additionally require that the same
thermodynamic parameters be specified for all thermodynamic states,
though their values may of course differ.

\subsection{Expanded ensembles}

In an expanded ensemble simulation~\cite{lyubartsev:jcp:1992:expanded-ensembles}, a single replica or ``walker'') samples pairs $(x,k)$ from a joint distribution of configurations $x \in \Gamma$ and state indices $k \in \{1,\ldots,K\}$ given by,
\begin{eqnarray}
\pi(x,k) &\propto& \exp[-u_k(x) + g_k] ,
\end{eqnarray}
where $g_k$ is an state-dependent weighting factor. 
This space is therefore a {\em mixed}, {\em generalized}, or {\em expanded} ensemble which samples from multiple thermodynamic ensembles simultaneously. 
{\color{red}[JDC: I really think the appropriate phrasing here should be: ``\ldots which samples from multiple thermodynamic states within the same thermodynamic ensemble.'']}
$g_k$ is chosen to give a specific weighting of each subensemble in the expanded ensemble, and is generally determined
through some iterative procedure~\cite{lyubartsev:jcp:1992:expanded-ensembles,marinari-parisi:europhys-lett:1992:simulated-tempering,wang-landau:prl:2001:wang-landau,park-ensign-pande:pre:2006:bayesian-weight-update,park-pande:pre:2007:choosing-weights-simulated-tempering,li-fajer-yang:jcp:2007:simulated-scaling,chelli:jctc:2010:optimal-weights-expanded-ensembles}. 
The set of $g_k$ is frequently chosen to give each thermodynamic ensemble equal probability (so-called ``perfect weights''~\cite{park-pande:pre:2007:choosing-weights-simulated-tempering}), in which case $g_k=-\ln Z_k$, but they can be set to arbitrary values as desired.

In the context of Gibbs sampling, an expanded ensemble simulation proceeds by alternating between sampling from the two conditional distributions,
\begin{eqnarray}
\pi(x | k) &=& \frac{q_k(x)}{\int_\Omega dx \, q_k(x)} = \frac{e^{-u_k(x)}}{\int_\Omega dx \, e^{-u_k(x)}} \\
\pi(k | x) &=& \frac{e^{g_k}q_k(x)}{\sum\limits_{k'=1}^K e^{g_{k'}}q_{k'}(x)} = \frac{e^{g_k - u_k(x)}}{\sum\limits_{k'=1}^K e^{g_{k'} - u_{k'}(x)}} .\label{equation:expanded-ensemble-gibbs-update}
\end{eqnarray}

In all but trivial cases, sampling from the conditional distribution
$\pi(x | k)$ must be done using some form of Markov chain Monte Carlo
sampler that generates correlated samples, due to the complex form of
$u_k(x)$ and the difficulty of computing the normalizing constant in
the denominator~\cite{jun-s-liu:mcmc}. Typically, Metropolis-Hastings
Monte
Carlo~\cite{metropolis:jcp:1953:metropolis-monte-carlo,hastings:biometrika:1970:metropolis-hastings}
or molecular dynamics is used~\cite{footnote3}, generating an updated
configuration $x^{(n+1)}$ that is correlated with the previous
configuration $x^{(n)}$.  

\subsection{Wang-Landau method}

In many cases, one wishes to compute the free energies by sampling from all states evenly. One method for this is the Wang-Landau algorithm.  
For each state, we keep a histogram with the counts of visiting each thermodynamic state, $h(i)$.  Each time we visit the state in the expanded ensemble simulation, we update the count by 1.
We also adjust the logarithmic weight of that state by $-\delta$, where $\delta$ is some small quantity that is generally reduced as the simulation proceeds.
\begin{eqnarray}
h_{\mathrm{new}}(i) = h_{\mathrm{old}}(i) + 1 \\
g_{i,\mathrm{new}} = g_{i,\mathrm{old}} - \delta
\label{eq:wang-landa}
\end{eqnarray}
This can be viewed as a transition matrix algorithm.  
{\color{red}[JDC: Cite an explanation of what this means?]}
We update the transition matrix {\color{red}[JDC: Which transition matrix?  Which element?  How is it updated? You can't mean that 1 is added to some element, since the rows (or columns) of a stochastic matrix have to \emph{sum} to one.]} by 1 whenever we accept a move to the state.

This algorithm is self-correcting; if a state is visited more frequently than it would through equal sampling at each state, its weight is decreased by the visits, and eventually, the additional weight accruing to that state falls back to the correct range.  
This proceeds until all states are visited with approximately even frequency.  
When we reach some rough level of flatness, we reduce $\delta$ by multiplying by some factor $k$, frequently chosen to be $1/2$, though this can lead to saturation of parts of the configuration space~\cite{Belardinelli2008}.
We continue this process until a desired level of flatness is achieved.  
At this point, the weights $g_k$ will be equal to the free energies $f_k$ within statistical noise.
{\color{red}[JDC: Is there any mathematical proof of this claim?]}

We seek a more robust calculation for the free energies.
{\color{red}[JDC: More robust than what?  Why is the previous algorithm not robust?]}  
If we are visiting all the states, nearby states sample similar phase space.
There must be some way to use information sampled at each state to gain information about the other states nearby.

We start with MBAR.  
MBAR is the provably minimum variance estimate of the free energy given $N_i$ samples from each of $i=1 \ldots K$ states. 
We can write an estimator for each free energy $f_i$ (up to an additive constant) using the self-consistent equation:
\begin{eqnarray}
\exp(-f_i) = \sum_{n=1}^N \frac{\exp(-u_i(x_n))}{\sum_k N_k \exp(f_k-u_k(x_n))}
\label{eq:MBAR}
\end{eqnarray}
{\color{red}[JDC: Might be good to use $\hat{f}_i$ notation, or $\hat{f}_i^{(n)}$ for estimate after $n$ data points.]}

Now assume that we have a decent estimate of the free energies $f_i$, and we wish to add one more sample, and then come up with an updated free energy estimate.  
Assuming we do not recalculate the free energies, we now have:
\begin{eqnarray}
%\exp(-f_i) &=& \sum_{n=1}^N \frac{\exp(-u_i(x_n))}{\sum_k N_k \exp(f_k-u_k(x_n))} \\ % duplicate of previous eqn
e^{-f_{i,\mathrm{new}}} &=& \sum_{n=1}^N \frac{\exp(-u_i(x_n))}{\sum\limits_k N_k \exp(f_{\mathrm{k,old}}-u_k(x_n))} + \nonumber \\
                          & & \frac{\exp(-u_i(x_{N+1}))}{\sum\limits_k N_k \exp(f_{\mathrm{k,old}}-u_k(x_{N+1}))} \\
e^{-f_{i,\mathrm{new}}} &=& e^{-f_{i,\mathrm{old}}} + \frac{\exp(-u_i(x_{N+1}))}{\sum_k N_k \exp(f_{\mathrm{k,old}}-u_k(x_{N+1}))} \\
&=& \exp(-f_{i,\mathrm{old}}) \left[1 + \frac{\exp(f_{i,\mathrm{old}}-u_i(x))}{\sum_k N_k \exp(f_{\mathrm{k,old}}-u_k(x))}\right] \nonumber \\
f_{i,\mathrm{new}} &=& f_{i,\mathrm{old}}- \ln\left[1 + \frac{\exp(f_{i,\mathrm{old}}-u_i(x))}{\sum_k N_k \exp(f_{\mathrm{k,old}}-u_k(x))}\right] \nonumber \\
                   &=& f_{i,\mathrm{old}}- \ln\left[1 + \frac{1}{N_i}\frac{N_i\exp(f_{i,\mathrm{old}}-u_i(x))}{\sum_k N_k \exp(f_{\mathrm{k,old}}-u_k(x))}\right] \nonumber \\ 
                   &\approx& f_{i,\mathrm{old}}- \frac{1}{N_i}\frac{N_i\exp(f_{i,\mathrm{old}}-u_i(x))}{\sum_k N_k \exp(f_{\mathrm{k,old}}-u_k(x))} 
\end{eqnarray} 
Where in the last step we expand the logarithm in the first order Taylor series, in the limit that $N_i$ is large, since the rightmost fraction must lie in the interval $[0,1]$.  
{\color{red}[JDC: You're also omitting the contribution from reoptimizing $f_k$ on the old data, the left-hand term involving the sum from $n=1$ to $N$ on the left-hand side of the first equation.  It's not clear to me that this contribution is not as large as the contribution from the new data point $N+1$.  I would be much happier if you could linearize \emph{both} the old sum from $n=1$ to $N$ \emph{and} the contribution from $N+1$ and solve for the updated $f_{i,new}$ using this linearized form.]}

This indicates that the appropriate accumulation term is proportional
to the number of samples collected so far, divided by the number of
states, rather than a fixed number, as suggested in the Wang-Landau
algorithm.  Interestingly, this is exactly the proposed accumulation
factor $1/t$ algorithm proposed by~\citet{Belardinelli2007}, which
states that the incrementor factor $\delta$ should asymptotically approach $< N/K$.

So at each step, we have a formula based on MBAR for how to update the
free energy in an adaptive manner given the information collected so
far.  To turn this into an all-state Wang-Landau-like algorithm,
however, we need to have some way to compute a histogram to determine
when the weights lead to sufficient flatness.  We could use a simple
discrete histogram update, as in the original Wang-Landau algorithm.
However, if we are computing the energies at all other states, we
could update all states by a weighted amount.  For example, we could
look at the conditional probability for $k$ give the sample $x$ given
by:
\begin{eqnarray}
\pi(i | x) &=& \frac{e^{g_k}q_i(x)}{\sum\limits_{k=1}^K e^{g_{k}}q_{k}(x)} = \frac{e^{g_i - u_i(x)}}{\sum\limits_{k=1}^K e^{g_{k} - u_{k}(x)}} .
\end{eqnarray}
This is the probability of a state $i$ given sampling from $x$.  Therefore, we can interpret it as amount we should increment each state $i$ when we perform a sample from any state. 

To make this more compact, we define the quantity
\begin{eqnarray}
\pi(i|N,x) = \frac{N_k\exp(f_i-u_i(x))}{\sum_k N_k \exp(f_k-u_k(x))}
\end{eqnarray}
Which closely resembles $\pi(i|x)$, and equals it in the limit where
all $N_k$ are equal, but takes into account the inequal weighting so
far.

This can be generalized to arbitrary weights $g_k$ that differ from
$g_k$ by $\Delta_k$.  This would be the case if we wish to sample
state $i$ by more than equal sampling proportional to
$\exp(\Delta_i)$.

\begin{eqnarray}
\pi(i|x,\vec{N},\vec{\Delta}) = \frac{N_k\exp((g_i-\Delta_i)-u_i(x))}{\sum_k N_k \exp((g_k-\Delta_k)-u_k(x))}
\end{eqnarray}

Where $\vec{N}$ is the list of all numbers of samples from each state
$N_k=1,\ldots,k$, and $\vec{\Delta}$ is the vector of all deviations
of the weights from equal sampling.

Given these definitions, we can define a new algorithm for iterative
updating of the weights and free energies by:
\begin{eqnarray}
N_{i,\mathrm{new}} = N_{i,\mathrm{old}}+\pi(i|x) \\
f_{i,\mathrm{new}} = f_{i,\mathrm{old}} - \delta \pi(i|x,\vec{N},\vec{\delta})
\end{eqnarray}
Where $\delta$ is a scaling factor that is decreased wherever the
flatness of the $\vec{N}$ reaches a desired level.  At the point when
$\delta < N/K$, we ignore the histogram, set $\delta = N/K$ at each
step, and simply update the free energies. 


\section{Applications} 
We compare different variants of the update scheme.

\begin{itemize}
\item Original Wang-Landau 
\item all-state Wang-Landau, weighted $N_k$ assumption, discrete histogram
\item all-state Wang-Landau, weighted $N_k$ assumption, fractional histogram
\end{itemize}

We analyze them by the time require to reach flatness, and the
convergence of energies over time.

We choose an initial $\delta=1$, and $k_{\mathrm{scale}} = 0.5$.  And
an initial flatness criteria of 0.7.  Once $\delta < N/K$, we stop
updating the histogram, and continue updating only the free energies using $\delta = N/K$.

We take 100 initial starting structures, and simulate 10 ns for each
of three molecule (UA methane, anthracene, and a dipole inversion).
We measure the variance across the 100 states at fixed numbers of samples.

%If we are currently at a state, we have visited it with probability 1.
%The probability of a neighboring state $B$. The probability of the
%next state over is $exp(\delta f - \delta u)$. The average of this
%foreign weighting is 

%\begin{eqnarray}
%\int P_B/P_A  P_A(x) dx \\
%\int \exp(\delta f -\delta u) \exp(f_A - u_A) \\ 
%\int \exp(f_B - u_B) = 1 
%\end{eqnarray}

%The variance is:
%\begin{eqnarray}
%& & \int \exp(\delta f -\delta u)**2 \exp(f_A - u_A) dx - (\int \exp(\delta f - \delta u) \exp(f_A - u_A))^2 \\
%&=& \int \exp(\delta f -\delta u)**2 \exp(f_A - u_A) dx - 1  \\
%&=& \int \exp(\delta f -\delta u) \exp(f_B - u_B) dx - 1  \\
%&=& \int \exp(f_C - u_C) dx - 1  = 0 \\
%\end{eqnarray}

%or 
%\begin{eqnarray}
%& & \int (P_B/P_A)^2 P_A dx - (\int (P_B/P_A) P_A dx )^2  \\
%&=& \int (P_B / P_A) P_B dx - 1  \\

%&=& \int \exp(\delta f -\delta u) \exp(f_B - u_B) dx - 1  \\
%&=& \int \exp(f_C - u_C) dx - 1  = 0 \\
%\end{eqnarray}

%Consider a system having a set $\Omega$ of possible
%microstates. Microstate $s(x,k) \in \Omega$ has Boltzmann probability
%$\pi_s(u_k(x),\vec{0}) = \exp(u_k(x)]/Z$, where $u_k$ is the reduced
%  potential of thermodynamic state $k$, and $Z = \sum_{s \in \Omega}
%  \exp(u_k(x)$. 

%We are interested in computing the ratio of normalization constants
%$Z_i/Z_j$, where $Z_i = \sum_{s=s(x|k=i)} \exp(u_k(x)$ and $Z_j =
%\sum_{s=s(x|k=j)} \exp(u_k(x)$.  If we are in equilibrium in the
%Markov chain, then

%\begin{eqnarray}
%\frac{T_{s(x|i)\rightarrow s(x|j)}{T_{s(x|j)\rightarrow s(x|i)} =  P(s(x|i)
%\end{eqnarray}


%the value of the Hamiltonian for microstate s, Z = exp[−\betaHs] is the partition function, s∈Ω
%and the notation ⃗0 will be explained shortly. Let S denote themacrosta te to which microstate
%s belongs, the relation S-to-s being one-to-many. The Boltzmann macrostate probabilities are


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END APPENDIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{weightedwl}

\end{document}

% LocalWords:  PRL Chodera QB VA  Monte nonparametric subensemble MBAR inequal
% LocalWords:  Belardinelli incrementor dx Microstate ôexp Hs Multistate UA
% LocalWords:  anthracene microstates
