%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DOCUMENT CLASS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Single-spaced, two-column with PRL look and style (easy on the eyes)
%\documentclass[aps,pre,twocolumn,superscriptaddress]{revtex4-1}
%\documentclass[aps,pre,twocolumn,superscriptaddress]{revtex4}

% Double-spaced, one-column style (for submission/review/editing)
\documentclass[aps,preprint,prl,superscriptaddress,showpacs]{revtex4}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PREAMBLE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

%\bibliographystyle{apsrevlong}
\bibliographystyle{apsrev}

% 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{\der}[3]{\frac{d^{#3} #1}{d #2^{#3}}}
\newcommand{\dlam}{\der{}{\lambda}{}}
\newcommand{\dudlam}{\der{u}{\lambda}{}}
\newcommand{\daudlam}{\der{\expect{u}}{\lambda}{}}
\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{\Var}[1] {\mathrm{var}\left( #1 \right)}
\newcommand{\Ept}[1] {{\mathrm E}\left[ #1 \right]}
\newcommand{\Eptk}[2] {{\mathrm E}_{#1}\left[ #2\right]}
\newcommand{\T}{\mathrm{T}}                                % T used in matrix transpose

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TITLE AND AUTHORS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Estimators for Entropy and Enthalpy of Molecular Transformations from Multistate Sampling}

 \author{Michael R. Shirts}
  \email{michael.shirts@virgina.edu}
 \affiliation{Department of Chemical Engineering, University of Virginia, Charlottesville, VA 22094-0471}

 \author{Jacob I. Monroe}
 \affiliation{Department of Chemical Engineering, University of Virginia, Charlottesville, VA 22094-0471}

 \author{John D. Chodera}
 \affiliation{Memorial Sloan-Kettering Cancer Center, New York, USA}
 
\date{\today}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BODY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

People are often interested in calculating the entropy difference
between two different chemical states.

Examples: the change in entropy between solvation, or the change in
entropy on binding.

There are a number of ways to calculate it.

\section{Theory}

The entropy can be defined most generally by the Gibbs equation for entropy:
\begin{eqnarray}
S = k_B \sum_m p_m \ln p_m
\end{eqnarray}
It can also be defined thermodynamically by: 
\[S = -\der{A}{T}{V,N_i}\]
in the canonical ensemble, or:
\[S = -\der{G}{T}{P,N_i}\] in the isobaric-isothermal ensemble.
To more generalize these approaches, we 

We define the \emph{reduced potential function} $u_i(\x)$ for defined statistical mechanical state $i$ to be
\begin{eqnarray}
u_i(\x) &=& \beta_i [ U_i(\x) + p_i V(\x) + \bfv{\mu}_i^\T \bfv{n}(\x) ] \label{equation:reduced-potential}
\end{eqnarray}

where $\x \in \bfm{\Gamma}$ denotes the configuration of the system
within a configuration space $\bfm{\Gamma}$, with volume $V(\x)$ (in
the case of a constant pressure ensemble) and $\bfv{n}(\x)$ the number
of molecules of each of $M$ components of the system (in the case of a
(semi)grand ensemble).  For each state $i$, $\beta_i$ denotes the
inverse temperature, $U_i(\x)$ the potential energy function (which
may include biasing weights), $p_i$ the external pressure, and
$\bfv{\mu}_i$ the vector of chemical potentials of the $M$ system
components. For each state, we define reduced free energy 
\begin{eqnarray}
f_i -\ln int_V e^{-u_i(\vec{x})} d\vec{x}
\end{eqnarray}
In the canonical ensemble $f_i = \beta A_i$, in the
isothermal-isobaric ensemble, $f_i = \beta G_i$, and in the grand
canonical ensemble $f_i = \beta P_i V$. We will generally use $F$ to
represent each of these quantities in the appropriate ensemble.  We
also define an reduced entropy $s = S/k_B$, so that $f_i = u_i-s_i$.

We then obtain the following identities.  We assume that derivatives
assuming that $P$, $N_i$, or other variables other than $S$ and $T$
are held constant, and using $U$ to represent the nonreduced
potential, which may be equal to $H$ or some other quantity depending
on the ensemble.

\begin{eqnarray}
\der{F}{T}{} &=& S  \\
\end{eqnarray}
\begin{eqnarray}
\der{f}{T}{} &=& \der{F/k_B T}{T}{}  \\
           &=& k_B^{-1}\left[-F T^{-2} + T^{-1}\der{F}{T}{}\right]  \\
           &=& -f/T - s/T  \\
           &=& -u/T
\end{eqnarray}

%F = H - TS
%F = H - dS/dH dS
%S = H/T - F/T 
%\beta F = \beta U - \beta TS
%[d(beta)/dT]^-1
%= [k_B^-1 d(T^-1)/dT]^-1
%= [-k_B^-1 T^-2)]^-1
%= -k_B T^2
%= -k_B T^2
%= -k_B^-1 beta^{-2}
%[d(beta)/dT]
%= k_B beta^{2}

%dG/dbeta = dG/dT[dT/dB] dG/dT(dB/dT)^-1
%dG/dbeta = -S*( -k_BT^2)
%dG/dbeta = k_BST^2
%dG/dbeta = S/k_B beta^-2 
%         = s/beta^2

% in summary dG/dT = -S
%            dG/dbeta = s/beta^2 
%            df/dT = h/T
%            df/dbeta = H
%Where $X$ stands for the variables that are held constant while the
%temperature is varied, and $F$ is either $G$, $A$, or $PV$ depending
%on the choice of $X$.  

We see that the temperature derivative of the free energy is the
entropy, whereas the temperature derivative of the reduced free energy
is actually the reduced energy (or enthalpy) divided by the
temperature.  Thus, one can find either the average energy or the
entropy can be written as the derivative of a free energy.  We can
find equivalent calculations with respect to $\beta$ 

\begin{eqnarray}
\der{f}{\beta} &=& \der{\beta F/}{\beta}  \\
           &=& F+\beta \der{F}{\beta} \\
           &=& F+\beta\der{F}{T}\der{\beta}{T}^{-1} \\
           &=& F-\beta S(-k_B T^2)  \\
           &=& F+TS  \\
           &=& U  \\
           &=& \frac{u}{\beta}
\end{eqnarray}
\der{F}{\beta}{} &=& \der{F}{T}{}\der{T}{\beta}{} \\
                 &=& S(-k_B^-1 beta^{-2}) \\
                 &=& -\frac{s}{\beta^2} 
\end{eqnarray}
Again, with the reduced free energy derivatives corresponding to
average potentials, and the free energy derivatives corresponding to
entropies.

We note that the free energy difference between two systems at
different temperatures is not the simple ratio of partition functions,
although the difference in reduced free energy can always be expressed
as a ratio of partition functions

\begin{eqnarray}
F_2 - F_1 = k_B \left[ T_2 ln Q_2 - T_1 Q_1 \right]
f_2 - f_1 = ln Q_2 - ln Q_1 = \frac{Q_2}{Q_1}
\end{eqnarray}

%F_2 - F_1 = T_2 f_2 - T_1 f_1
%
%F_2 = F_1 = \int_{T=T_0}^{T_1} \der{F}{T} dT 
%          = \int \der{Tf}{T} dT = \int \left[f + T\df{f}{T}] dT 
%          = \int \left[f + u\right\right] 
%
% 
%If we choose a reference state for which u=0, then we have 

The reduced free energy of the system can be directly computed from
statistical mechanics as $f = \int e^(-u) d\vec{x}$ where here the
$\vec{x}$ is an integral of all degrees of freedom of the system,
which include the values of $N$ and $V$ as well as the coordinates of
particle.

There are a number of different ways to calculate the free energy,
entropy, and energy from a simulation or series of simulations of a
system.  All three of these quantites, at least classically, are only
defined up to a constant, so we are generally interested in estimating
the difference of free energies, entropies, and enthalpies given two
simulations or a series of simulations.

In this paper, we compare a number of ways to calculate the entropy
difference between two systems.  We specifically look at:

\begin{enumerate}
\item Subtraction of endpoint energies
\item Finite differences in temperature
\item Thermodynamic integration
\item Energy-entropy decomposition using MBAR
\item ``virtual'' Finite difference using reweighted ensembles
\item Finite difference using multiple temperature ensembles
\end{enumerate}
%using higher order derivatives with O(dT^4) estimator (varying delta T)
%\item Finite difference by using higher order derivatives reweighting O(dT^4) estimator
%- TI
%- TI with MBAR reweighting 

We test these estimators on three systems:
\begin{enumerate}
\item a set of harmonic oscillators for which the answers are known analytically
\item the solvation of a small molecule 
\item the binding free energy of a host-guest system
\end{enumerate}


\subsection{Endpoint Calculations}

A straightforward way to calculate the reduced entropy is to first
calculate the difference in reduced energies and the initial and final
states, and calculate the free energy difference between the initial
and final states.  The reduced entropy difference can then calculated
from $\Delta s_i = \Delta h_i - \Delta f_i$, and divide by $k_B$ to
obtain the entropy. Assuming both initial and final states are at the
same temperature, then then this is equivalent to $\Delta S_i =
\frac{1}{T} (U_i - F_i)$.  We can use MBAR to calculate $f_i$, and use
either the direct average to calculate $u_i$ and $u_j$ from the data
at states $i$ and $j$ respectively, or use MBAR to calculate $u_i$ and
$u_j$ using data from all intermediate states.

\subsection{MBAR Estimator} 

Once the solve for the normalization constants $f_k$, then we can
compute the expectation value of any observable from simulations using
Eq. X of Chodera and Shirts.
\begin{eqnarray}
\expect{u}_i = \frac{\int u_i(x) \exp(-u_i(x)) dx}{\sum_{i=1}^K\exp(f_k-u_k(x) )} \\ 
\end{eqnarray}
%FIX THIS

\subsection{Thermodynamic Integration}

If $u(x)$ is a function of some coupling parameter $\lambda$ that
connects the two states of interest, such that $u(x,\lambda=0) =
u_0(x)$ and $u(x,\lambda=1) = u_1(x)$, then we will have 

\begin{eqnarray*}
\der{f}{\lambda}{} &=& -\dlam \ln \sum_x \exp(-u(x)) \\
                    &=& -\frac{\dlam \sum_x \exp(-u(x))}{\sum_x \exp(-u(x))} \\
                    &=& \frac{\sum_x \dudlam \exp(-u(x))}{\sum_x \exp(-u(x))} \\ 
                    &=& \expect{\dudlam} \\ 
\end{eqnarray*}
We also note that 
\begin{eqnarray*}
\dlam Q^{-1} &=& -Q^{-2} \dlam Q \\
             &=& -Q^{-2} \dlam \sum_x \exp(-u(\lambda,x)) \\
             &=& -Q^{-2} \sum_x -\dudlam \exp(-u(\lambda,x)) \\
             &=&  Q^{-1} \sum_x \frac{\dudlam \exp(-u(\lambda,x))}{Q} \\
             &=&  Q^{-1} \expect{\dudlam}
\end{eqnarray*}

We then look at $\lambda$ derivatives of expectatons that are not functions of $u$.
\begin{eqnarray*}
\der{\expect{A}}{\lambda}{}    &=& \dlam \left[\sum_x A \exp(-u(x))\right] Q^{-1} \\ 
                             &=& \left[\sum_x \dlam A \exp(-u(x))\right] Q^{-1} + \\
                             & &   \left[\sum_x A \exp(-u(x))\right] \dlam Q^{-1}  \\
                             &=& -\left[\sum_x A \dudlam|_x \exp(-u(x)) \right] Q^{-1} + \\
                             & &    \left[\sum_x A \exp(-u(x))\right] Q^{-1} \expect{\dudlam} \\
                             &=& -\expect{A\dudlam} + \expect{A}\expect{\dudlam} \\ 
                             &=& -\Cov{A,\dudlam}         
\end{eqnarray*}

Next, we look at the derivative with respect to $u$.
\begin{eqnarray*}
\daudlam  &=& \dlam \left[\sum_x u \exp(-u(x))\right] Q^{-1} \\ 
         &=& \left[\sum_x \dlam u \exp(-u(x))\right] Q^{-1} + \\
         & &   \left[\sum_x u \exp(-u(x)) \right] \dlam Q^{-1}  \\
         &=& \left[\sum_x \dudlam|_x \exp(-u(x)) - u\dudlam|_x \exp(-u(x))\right] Q^{-1} + \\
         & &    \left[\sum_x u \exp(-u(x))\right]  Q^{-1} \expect{\dudlam} \\
         &=& \expect{\dudlam} - \expect{u\dudlam}+ \expect{u}\expect{\dudlam} \\
         &=& \expect{\dudlam} - \Cov{u,\dudlam}
\end{eqnarray*}

The reduced entropy $s = \beta T S = S/k$.  From $F = U-TS$, we have $f =
u-s$, so the derivative $s$ with respect to $\lambda$ is:
\begin{eqnarray}
\der{s}{\lambda}{} &=& \daudlam - \der{f}{\lambda}{} \\
\der{s}{\lambda}{} &=& \daudlam - \expect{\dudlam} \\
                    &=& - \Cov{u,\dudlam}
\end{eqnarray}

We can compute the covariance along a path from the MBAR estimator of
the covariance estimate.

If we assume a path linear in spring constant $K$ and linear in offset
$O$, such that $K(\lambda) = c_k \left[k_0 + \lambda(k_1 - k_0)\right]$ and $O(\lambda) =
c_o\lambda$, and using a reduced spring constant $k(\lambda) = \beta K(\lambda)$,
\begin{eqnarray}
u(x,\lambda) &=& \beta U(x,\lambda) = \frac{\beta c_k \left[k_0 + \lambda(k_1 - k_0)\right]}{2}(x-c_o\lambda)^2 \nonumber \\
\frac{du}{d\lambda} &=& \frac{\beta c_k(k_1-k_0)}{2}(x-c_o\lambda)^2 + c_k c_o \left[k_0 + \lambda(k_1 - k_0)\right](x-c_o\lambda) \nonumber \\
                   &=& \frac{k-\beta c_k k_o}{2\lambda}(x-O)^2 + \frac{kO}{\lambda}(x-O)  \nonumber\\
\expect{\frac{du}{d\lambda}} &=& \int_{-\infty}^{\infty}\frac{du(x,\lambda)}{d\lambda}\exp(f(\lambda)-u(x,\lambda))dx  \nonumber\\
                           &=& \int_{-\infty}^{\infty}\left[\frac{k-\beta c_k k_o}{2\lambda}(x-O)^2 + \frac{kO}{\lambda}(x-O)\right] \sqrt{k/2\pi} \exp\left(- \frac{k}{2}(x-O)^2\right) dx \nonumber \\
                           &=& \sqrt{k/2\pi} \left[\int_{-\infty}^{\infty}\frac{k-\beta c_k k_o}{2\lambda}(x-O)^2 \exp\left(-\frac{k}{2}(x-O)^2\right) dx\right] \nonumber \\ 
                           &+& \sqrt{k/2\pi} \left[\int_{-\infty}^{\infty}\frac{kO}{\lambda}(x-O)\exp\left(-\frac{k}{2}(x-O)^2\right) dx\right] \nonumber\\ 
                           &=& \frac{k-\beta c_k k_0}{2k\lambda}  \nonumber\\
                           &=& \frac{1}{2}\frac{1}{\lambda + \frac{k_0}{k_1 - k_0}}  
\end{eqnarray}
Note that $\expect{du/d\lambda}$ is independent of both $c_k$ and $o_k$. 
\subsubsection{Single state}

\subsection{Finite differences}
We can perform our series of simulations at two different
temperatures, find the free energy at these temperatures, and perform
a finite difference calculation. We perform a central difference
approximation, as this has error $\mathcal{O}(\Delta T^2)$ rather than
$\mathcal{O}(\Delta T)$.
\begin{equation}
S(T) \approx \frac{F(T+\Delta T) - F(T-\Delta T)}{2T}
\end{equation}

\subsection{Finite Difference using MBAR from the expectations}
Using MBAR, we can use reweighting techniques to find the free energy
at temperatures $T+\Delta T$ and $T-\Delta T$, and subtract directly.

Because MBAR gives us the covariance matrix, if we analyze the finite
difference equations using MBAR, we can determine the covariance, and
compute the uncertainty in $\Var{F_2 - F_1}$ directly from the
simulations, including covariance.

\subsection{Finite Difference derivative from MBAR and a single simulation}

Using MBAR, we can take a single simulation estimate at both $T-\Delta
T$ and $T+\Delta T$ using reweighting from a simulation at $T$.

\subsection{Finite Difference derivative from MBAR}

\subsection{Alternate estimators}

\subsubsection{Analytical derivative with resepct to the MBAR estimator} 

Pappu and coworkers proposed a new estimator for entropy differences
between thermodynamic states~\cite{pappu:jpcb:2010:entropy-estimator}
based on the multistate Bennett acceptance
ratio~\cite{shirts-chodera:jcp:2008:mbar}.  The identity that forms
the basis of their proposed estimator, given in Eq.~33 of
\cite{pappu:jpcb:2010:entropy-estimator}, is written in reduced
(dimensionless) units as
\begin{eqnarray}
\lefteqn{s_i = - f_i + \sum_{j=1}^K N_j \expect{\mu_i u_i}_j} \nonumber \\
&& \mbox{} + \left[ \sum_{j=1}^K N_j \frac{T_i}{T_j} (\expect{\mu_i u_j}_j - \expect{\mu_i}_j \expect{u_j}_j)  \right. \nonumber \\
&& \left. + \sum_{j=1}^K \sum_{k=1}^K N_j N_k \frac{T_i}{T_k} \left[ (f_k + s_k) \expect{\mu_i \mu_k}_j - \expect{\mu_i \mu_k u_k}_j \right] \right] \nonumber
\end{eqnarray}
where $s_i$ denotes the entropic contribution to the free energy in dimensionless units (having subsumed the temperature $T_i$ written in Eq.~33 of \cite{pappu:jpcb:2010:entropy-estimator}), $f_i$ is the dimensionless free energy, $N_i$ denotes the number of samples from state $i$, $u_i(x)$ is the reduced potential energy of state $i$, and $\mu_i(x)$ is the weighting function
\begin{eqnarray}
w_i(x) &=& \frac{\exp[f_i - u_i(x)]}{\sum\limits_{k=1}^K \exp[f_k - u_k(x)]} \nonumber
\end{eqnarray}
Replacing expectations $\expect{A}_j$ with empirical estimators $N_j^{-1} \sum_{n=1}^{N_j} A(x_{jn})$, we obtain the estimator for $s_i$:
\begin{eqnarray}
\lefteqn{\hat{s}_i = - \hat{f}_i + \sum_{j=1}^K \sum_{n=1}^{N_j} \mu_i(x_{jn}) v_i(x_{jn})} \label{equation:rohit-entropy} \\
&& \mbox{} + \left[ \sum_{j=1}^K  \frac{T_i}{T_j} \sum_{n=1}^{N_j} \mu_i(x_{jn}) v_j(x_{jn}) \right. \nonumber \\
&& \mbox{} - \sum_{j=1}^K \frac{T_i}{T_j} N_j^{-1} \left( \sum_{n=1}^{N_j} \mu_i(x_{jn}) \right) \left( \sum_{n=1}^{N_j} v_j(x_{jn}) \right)  \nonumber \\
&& \mbox{} + \sum_{j=1}^K \sum_{k=1}^K N_k \frac{T_i}{T_k} (\hat{f}_k + \hat{s}_k) \sum_{n=1}^{N_j} \mu_i(x_{jn}) \mu_k(x_{jn}) \nonumber \\
&& \left. - \sum_{j=1}^K \sum_{k=1}^K N_k \frac{T_i}{T_k} \sum_{n=1}^{N_j} \mu_i(x_{jn}) \mu_k(x_{jn}) v_k(x_{jn}) \right] \nonumber
\end{eqnarray}
Because $\hat{s}_i$ appears on both the left-hand side (LHS) and right-hand side (RHS), these equations must be solved self-consistently for the $\hat{s}_i$; this is done by fixing the $\hat{f}_i$ determined from MBAR and iteratively updating the estimate for $\hat{s}_i$ on the LHS~\cite{pappu:jpcb:2010:entropy-estimator}.
Note that only relative entropies, differences of the form $\Delta s_{ij} \equiv s_j - s_i$, are meaningful, just as only relative free energy differences $\Delta f_{ij} \equiv f_j - f_i$ are can be estimated because MBAR provides estimates of $f_i$ only up to an arbitrary additive constant.

The reduced free energy $f_i$, average reduced energy $u_i = \expect{v_i}_i$, and reduced entropy $s_i$ obey the Helmholtz relationship
\begin{eqnarray}
f_i &=& u_i - s_i \nonumber
\end{eqnarray}
An alternative estimator for the reduced entropy $s_i$ is therefore given by
\begin{eqnarray}
\hat{s}_i &=& - \hat{f}_i + \hat{u}_i = -\hat{f}_i + \sum_{n=1}^N w_{ni} u_i(x_n) \label{equation:simple-entropy-estimator}
\end{eqnarray}
The free energy $f_i$ (up to an additive constant across all states) can easily be estimated with MBAR using Eq.~11 of \cite{shirts-chodera:jcp:2008:mbar}, and $u_i$ can be estimated using Eq.~15 of the same reference.  

This then provides an alternative, simpler approach to the estimation of the reduced entropy $s_i$ over Eq.~\ref{equation:rohit-entropy}.
In fact, we recognize that the estimator on the RHS of Eq.~\ref{equation:simple-entropy-estimator} corresponds to the first two terms on the RHS of Eq.~\ref{equation:rohit-entropy}, as $\hat{u}_i \equiv \sum_{j=1}^K \sum_{n=1}^{N_j} w_i(x_{jn}) v_i(x_{jn})$.
The remaining terms (those enclosed in brackets), might be surmised to correspond to a quantity that, once the iterative estimation of these equations has converged, is identically or nearly zero.

To test both entropy estimators, we consider a system consisting of a simple one-dimensional harmonic oscillator, for which all thermodynamic properties can be computed exactly.
The system has reduced potential
\begin{eqnarray}
v(x) &=& \frac{\beta K}{2} (x - x_0)^2 \nonumber
\end{eqnarray}
with $\beta = 1$ and $K = 0.2$.
For each experiment, 40 independent and identically distributed (\emph{i.i.d.}) samples were collected from each of $K = 7$ thermodynamic states, differing only in x0, where x0 was one of $\{0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8\}$.
1000 independent trials of this experiment were run, and statistics were collected on the error in the estimated entropy difference between state 1 ($x_0 = 0$) and state 7 ($x_0 = 1.8$), where the true thermodynamic differences are identically zero ($\Delta s_{17} = \Delta u_{17} = \Delta f_{17} = 0$).  

The estimator in Eq.~\ref{equation:rohit-entropy} (Eq.~30 of Ref.~\cite{pappu:jpcb:2010:entropy-estimator} had a standard deviation of $0.071 \pm 0.001$ and bias of $0.001 \pm 0.002$, while the estimator in Eq.~\ref{equation:simple-entropy-estimator} had a standard deviation of $0.070 \pm 0.001$ and a bias of $0.000 \pm 0.002$.
(Uncertainty estimates over replicates were obtained from expressions in Sections II.A and II.B of \cite{ahn-fessler:2003:standard-errors}.)
The difference is statistically indistinguishable and suggests that the simple estimator of Eq.~\ref{equation:simple-entropy-estimator} is as efficient as the more complex estimator (Eq.~\ref{equation:rohit-entropy}), as would be expected from the fact that the MBAR estimates for $f_i$ and $u_i$ are asymptotically efficient~\cite{shirts-chodera:jcp:2008:mbar}.

The bracketed terms in Eq.~\ref{equation:rohit-entropy} have a mean of $-0.0004 \pm 0.0003$, suggesting that while the expectation of the quantity in brackets is zero.
While it would be difficult to prove that the bracketed terms increase the asymptotic variance of the estimator over that of the simple estimator in Eq.~\ref{equation:simple-entropy-estimator}, the terms have a standard deviation of $0.0108 \pm 0.0002$, over the replicates due to statistical or numerical noise, suggesting this contributes some small amount to an increase in variance of the estimator in Eq.~\ref{equation:rohit-entropy}.

We also compared the two entropy estimators in the estimation of the entropy change in transferring 3-methylindole from vacuum to water in explicit solvent.
A set of 11 alchemical intermediates were used in which the interactions of the solute and solvent were decoupled in stages, and 50 ns trajectories in the isothermal-isobaric ensemble were generated, with potential energies at all alchemical intermediates stored every 1 ps.
MBAR was applied to the entire dataset to estimate the \emph{target} entropy difference between the vacuum and solvated endstates.
The dataset was then divided into 100 sequential blocks of 500 snapshots each, and both entropy estimators were applied to each block.
The standard deviation and bias of each estimator were computed from the deviation from the \emph{target} entropy difference over the blocks.
Here, the bias of the estimator in Eq.~\ref{equation:rohit-entropy} was significant ($-1.6 \pm 0.4$) compared to that of the simpler estimator of Eq.~\ref{equation:simple-entropy-estimator} ($0.0 \pm 0.3$).  
The standard deviation of the estimator of Eq.~\ref{equation:rohit-entropy} ($-4.0 \pm 0.3$) was also larger than that of the simpler estimator of Eq.~\ref{equation:simple-entropy-estimator} ($3.1 \pm 0.2$). 
%MRS: this doesn't make sense -- you have a negative standard deviation.

Together with the extra complexity of including these terms, this
suggests that users of MBAR interested in entropy estimates should
prefer the simpler estimator of
Eq.~\ref{equation:simple-entropy-estimator} over the more complex
estimator of Eq.~\ref{equation:rohit-entropy}.

We also note that the statistical error in entropy differences
estimated through the simple expression in
Eq.~\ref{equation:simple-entropy-estimator} can be straightforwardly
obtained from the covariance estimation equation given in Eq.~10 of
\cite{shirts-chodera:jcp:2008:mbar}.  As described in Section IV of
\cite{shirts-chodera:jcp:2008:mbar}, we first augment the $N \times K$
weight matrix by another $K$ columns corresponding to the reduced
potential $v(x)$ as an observable, such that the corresponding
$q_{K+i}(x) = v_i(x) q_i(x)$ for column $K+i$.  The corresponding
weight column is given by
\begin{eqnarray}
W_{n,K+i} &=& \hat{c}_{K+i}^{-1} \frac{u_i(x_n) \exp[-u_i(x_n)]}{\sum\limits_{k=1}^K N_k \exp[\hat{f}_k - u_k(x_n)]}
\end{eqnarray}

Making note of the fact that, in terms of the log normalization
constants $\theta_i = \ln c_i$ defined in
\cite{shirts-chodera:jcp:2008:mbar}, we have
\begin{eqnarray}
\Delta f_{ij} &=& - \ln c_j + \ln c_i = e^{\theta_i - \theta_j} \nonumber \\
\Delta u_{ij} &=& \frac{c_{K+j}}{c_j} - \frac{c_{K+i}}{c_i} = e^{\theta_{K+j} - \theta_j} - e^{\theta_{K+i} - \theta_i} \nonumber \\
\Delta s_{ij} &=& \ln \frac{c_j}{c_i} + \frac{c_{K+j}}{c_j} - \frac{c_{K+i}}{c_i} = \theta_j - \theta_i + e^{\theta_{K+j} - \theta_j} - e^{\theta_{K+i} - \theta_i} \nonumber 
\end{eqnarray}
\begin{eqnarray}
\delta^2 \Delta f_{ij} &=& \hat{\Theta}_{i,i} + \hat{\Theta}_{j,j} - 2 \hat{\Theta}_{i,j} \nonumber \\
\delta^2 \Delta u_{ij} &=& \left[\begin{array}{rrrr} \frac{\partial \Delta u_{ij}}{\partial \theta_i} & \frac{\partial \Delta u_{ij}}{\partial \theta_j} & \frac{\partial \Delta u_{ij}}{\partial \theta_{K+i}} & \frac{\partial \Delta u_{ij}}{\partial \theta_{K+j}}\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r} \frac{\partial \Delta u_{ij}}{\partial \theta_i} \\ \frac{\partial \Delta u_{ij}}{\partial \theta_j}\\ \frac{\partial \Delta u_{ij}}{\partial \theta_{K+i}}\\ \frac{\partial \Delta u_{ij}}{\partial \theta_{K+j}}\end{array}\right] \nonumber \\
&=& \left[\begin{array}{rrrr}u_i & -u_j & -u_i & u_j\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r}u_i\\-u_j\\-u_i\\u_j\end{array}\right] \nonumber \\
\delta^2 \Delta s_{ij} &=& \left[\begin{array}{rrrr} \frac{\partial \Delta s_{ij}}{\partial \theta_i} & \frac{\partial \Delta s_{ij}}{\partial \theta_j} & \frac{\partial \Delta s_{ij}}{\partial \theta_{K+i}} & \frac{\partial \Delta s_{ij}}{\partial \theta_{K+j}}\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r} \frac{\partial \Delta s_{ij}}{\partial \theta_i} \\ \frac{\partial \Delta s_{ij}}{\partial \theta_j}\\ \frac{\partial \Delta s_{ij}}{\partial \theta_{K+i}}\\ \frac{\partial \Delta s_{ij}}{\partial \theta_{K+j}}\end{array}\right] \nonumber \\
&=& \left[\begin{array}{rrrr}(-1+u_i) & (1-u_j) & -u_i & u_j\end{array}\right] \tilde{\bf \Theta} \left[\begin{array}{r}(-1+u_i)\\(1-u_j)\\-u_i\\u_j\end{array}\right] \nonumber 
\end{eqnarray}
where 
\begin{eqnarray}
\tilde{\bf \Theta} &\equiv& \left[\begin{array}{llll}
\hat{\Theta}_{i,i} & \hat{\Theta}_{i,j} & \hat{\Theta}_{i,K+i} & \hat{\Theta}_{i,K+j} \\
\hat{\Theta}_{j,i} & \hat{\Theta}_{j,j} & \hat{\Theta}_{j,K+i} & \hat{\Theta}_{j,K+j} \\
\hat{\Theta}_{K+i,i} & \hat{\Theta}_{K+i,j} & \hat{\Theta}_{K+i,K+i} & \hat{\Theta}_{K+i,K+j} \\
\hat{\Theta}_{K+j,i} & \hat{\Theta}_{K+j,j} & \hat{\Theta}_{K+j,K+i} & \hat{\Theta}_{K+j,K+j}
\end{array}\right] \nonumber
\end{eqnarray}
and the asymptotic covariance matrix estimate $\hat{\bf \Theta}$ is computed according to Eq.~8 of \cite{shirts-chodera:jcp:2008:mbar}.

Estimators for dimensionless free energy differences $\Delta
\hat{f}_{ij}$, entropy differences $\Delta \hat{s}_{ij}$, and energy
differences $\Delta \hat{u}_{ij}$ using MBAR have been available in
the freely-provided GPL Python software package pyMBAR, available at
\url{http://simtk.org/home/pymbar} since 2008.

\section{appendix}

If $u(x)$ is a function of $\lambda$, then:
\begin{eqnarray*}
\der{f}{\lambda}{} &=& -\dlam \ln \sum_x \exp(-u(\lambda,x)) \\
                    &=& -\frac{\dlam \sum_x \exp(-u(\lambda,x))}{\sum_x \exp(-u(\lambda,x))} \\
                    &=& \frac{\sum_x \dudlam \exp(-u(\lambda,x))}{\sum_x \exp(-u(\lambda,x))} \\ 
                    &=& \expect{\dudlam} \\ 
\end{eqnarray*}
We also note that 
\begin{eqnarray*}
\dlam Q^{-1} &=& -Q^{-2} \dlam Q \\
             &=& -Q^{-2} \dlam \sum_x \exp(-u(\lambda,x)) \\
             &=& -Q^{-2} \sum_x -\dudlam \exp(-u(\lambda,x)) \\
             &=&  Q^{-1} \sum_x \frac{\dudlam \exp(-u(\lambda,x))}{Q} \\
             &=&  Q^{-1} \expect{\dudlam}
\end{eqnarray*}

Next, we look at the derivative with respect to $u$.
\begin{eqnarray*}
\daudlam  &=& \dlam \left[\sum_x u(\lambda,x) \exp(-u(\lambda,x))\right] Q^{-1} \\ 
         &=& \left[\sum_x \dlam u(\lambda,x) \exp(-u(\lambda,x))\right] Q^{-1} + \\
         & &   \left[\sum_x u(\lambda,x) \exp(-u(\lambda,x)) \right] \dlam Q^{-1}  \\
         &=& \left[\sum_x \left.\dudlam\right|_x \exp(-u(\lambda,x)) - u\left.\dudlam\right|_x \exp(-u(\lambda,x))\right] Q^{-1} + \\
         & &    \left[\sum_x u(\lambda,x) \exp(-\lambda,u(x))\right]  Q^{-1} \expect{\dudlam} \\
         &=& \expect{\dudlam} - \expect{u\dudlam}+ \expect{u}\expect{\dudlam} \\
         &=& \expect{\dudlam} - \Cov{u,\dudlam}
\end{eqnarray*}

The reduced entropy $s = \beta T S = S/k$.  From $F = U-TS$, we have $f =
u-s$, so the derivative $s$ with respect to $\lambda$ is:
\begin{eqnarray}
\der{s}{\lambda}{} &=& \daudlam - \der{f}{\lambda}{} \nonumber \\
\der{s}{\lambda}{} &=& \daudlam - \expect{\dudlam} \nonumber \\
                   &=& -\Cov{u,\dudlam}
\end{eqnarray}

\begin{table}
\begin{tabular}{|l|ccc|}
\hline
                           & $\Delta f (k_BT)$   & $\Delta s (k_BT)$   &  $\Delta u (k_BT)$   \\ 
\hline
Method                        &                   &         &                \\
\hline    
Analytical                    &  1.151292         & -1.151292            &   0                 \\          
MBAR Expectations             &  1.14 $\pm$ 0.11  & -1.14 $\pm$ 0.11     &  -0.003 $\pm$ 0.023 \\
Endpoint Expectations         &  1.14 $\pm$ 0.11  & -1.14 $\pm$ 0.11     &   0.039 $\pm$ 0.032 \\
Thermodynamic Integration     &  1.20 $\pm$ 0.12  & -1.19 $\pm$ 0.20     &   0.09 $\pm$ 0.16 \\
Finite Diff. ($s = -\beta \der{F}{T}{}$) &  1.14 $\pm$ 0.11  & -1.00 $\pm$ 1.45     &    0.14  $\pm$ 1.45 \\ 
Finite Diff. ($u = T\der{f}{T}{}$) &  1.14 $\pm$ 0.11  & -0.98 $\pm$ 1.45     &    0.16  $\pm$ 1.45 \\ 
Virtual Finite Diff.          &  1.14 $\pm$ 0.11  & -1.14 $\pm$ 0.11     &  -0.003 $\pm$ 0.023 \\
Finite Diff. (MBAR)           &  
\end{tabular}
\caption{\label{tab:harmonic} Entropy and enthalpy for harmonic oscillators.}
\end{table}


\begin{table}
\begin{tabular}{|l|ccc|}
                              & $\Delta G (kJ/mol)$   & $\Delta S (J/mol/K)$   &  $\Delta H (kJ/mol)$   \\ 
\hline
Method                        &                       &                &                  \\
\hline    
MBAR Expectations             &   9.033  $\pm$ 0.029  & -35 $\pm$ 3    &   -1.4 $\pm$ 1.0 \\
Endpoint Expectations         &   9.033  $\pm$ 0.029  & -43 $\pm$ 6    &   -3.8 $\pm$ 1.7 \\
Thermodynamic Integration     &   9.211  $\pm$ 0.031  & -33 $\pm$10    &   -0.7 $\pm$ 2.9 \\
Finite Diff. ($10 K$)         &   9.033  $\pm$ 0.029  & -40 $\pm$ 4    &   -2.9 $\pm$ 1.0 \\ 
Finite Diff. (MBAR, $10 K$)   &   9.033  $\pm$ 0.029  & -40 $\pm$ 4    &   -2.9 $\pm$ 1.0 \\ 
Virtual Finite Diff.          &   9.033  $\pm$ 0.029  & -35 $\pm$ 3    &   -1.4 $\pm$ 1.0 \\
\end{tabular}
\caption{\label{tab:methane} Entropy and enthalpy for methane solvation.}
\end{table}


\begin{table}
\begin{tabular}{|l|c|c|}
Finite Diff.        &   $1 K$   & -60 $\pm$ 39  \\
Finite Diff.        &  $10 K$   & -40 $\pm$ 4   \\
Finite Diff.        &  $20 K$   & -37 $\pm$ 2   \\
Finite Diff. (MBAR) &   $1 K$   & -60 $\pm$ 39  \\
Finite Diff. (MBAR) &  $10 K$   & -40 $\pm$ 4   \\
Finite Diff. (MBAR) &  $20 K$   & -37 $\pm$ 2   \\
Virtual Finite Diff & $0.1 K$   & -35 $\pm$ 3   \\
Virtual Finite Diff &   $1 K$   & -35 $\pm$ 3   \\
Virtual Finite Diff &   $5 K$   & -19 $\pm$ 6   \\
Virtual Finite Diff &  $10 K$   &  66 $\pm$ 16  \\
\end{tabular}
\caption{Finite difference (virtual and actual). Turn into chart}
\end{table}


\begin{table}
\begin{tabular}{|l|c|c|c|}
                              & $\Delta S_{solvent}$   & $\Delta S_{complex}$   &  $\Delta S_{binding}$   \\ 
\hline
Method                        &                       &                &                  \\
\hline    
MBAR Expectations             &   180 $\pm$  40  & -20 $\pm$ 34   &  -200 $\pm$ 53 \\
Endpoint Expectations         &   110 $\pm$  34  & -32 $\pm$ 55   &   -78 $\pm$ 65 \\
Thermodynamic Integration     &   310 $\pm$ 150  & -30 $\pm$ 150  &  -340 $\pm$ 210 \\
%Finite Difference             &  
%Finite Difference (MBAR)      &  
\end{tabular}
\caption{Host/Guest, all in J/mol*K}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIBLIOGRAPHY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{prsty} 
\bibliography{comment-on-entropy-estimator}

\end{document}
