# 8.6: The "blue moon'' Ensemble Approach

• • Mark Tuckerman
• New York University
$$\newcommand{\vecs}{\overset { \rightharpoonup} {\mathbf{#1}} }$$ $$\newcommand{\vecd}{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$$$\newcommand{\AA}{\unicode[.8,0]{x212B}}$$

The term "blue moon'' in the present context describes rare events, i.e. events that happen once in a blue moon. The blue moon ensemble approach was introduced by Ciccotti and coworkers as a technique for computing the free energy profile along a reaction coordinate direction characterized by one or more barriers high enough that they would not likely be crossed in a normal thermostatted molecular dynamics calculation.

Suppose a process of interest can be monitored by a single reaction coordinate $$q_1=f_1({\bf r}_1,...,{\bf r}_N)$$ so that eqns. (29) and (30) reduce to

\begin{align} P(s) &= { {C_N \over Q(N,V,T)}\int\;d^N{\bf p}\;d^N{\bf r}e^{-\beta H({\bf p},{\bf r})}\delta(f_1({\bf r}_1,...,{\bf r}_N)-s) } \\[4pt] &= { {1 \over N!\lambda^{3N} Q(N,V,T)}\int\;d^N{\bf r}e^{-\beta U({\bf r})}\delta(f_1({\bf r}_1,...,{\bf r}_N)-s) }\\[4pt] A (s) &=-kT\ln P(s) \label{31} \end{align}

The "1'' subscript on the value $$s$$ of $${q_1}$$ is superfluous and will be dropped throughout this discussion. In the second line, the integration over the momenta has been performed giving the thermal prefactor factor $$\lambda ^{3N}$$. In the blue moon ensemble approach, a holonomic constraint $$\sigma({\bf r}_1,...,{\bf r}_N) =f_1({\bf r}_1,...,{\bf r}_N)-s$$ is introduced in a molecular dynamics calculation as a means of "driving'' the reaction coordinate from an initial value $${s_i}$$ to a final value $${s_f}$$ via a set of intermediate points $${s_1,...,s_n }$$ between $${s_i}$$ and $${s_f}$$. Unfortunately, the introduction of a holonomic, constraint does not yield the single $$\delta$$-function condition $$\delta(\sigma({\bf r}) =\delta(f_1({\bf r})-s)$$, where $${{\bf r}\equiv {\bf r}_1,...,{\bf r}_N }$$ required by Equation \ref{31} but rather the product of $$\delta$$-functions $$\delta(\sigma({\bf r}))\delta(\dot{\sigma}({\bf r},{\bf p}))$$, since both the constraint and its first time derivative are imposed in a constrained dynamics calculation. We will return to this point a bit later in this section. In addition to this, the blue moon ensemble approach does not yield $$A(s)$$ directly but rather the derivative

${dA \over ds} = -{kT \over P(s)}{dP \over ds} \label{32}$

from which the free energy profile $$A (q)$$ along the reaction coordinate and the free energy difference $$\Delta A = A(s_f)-A(s_i)$$ are given by the integrals

$A(q) = A(s_i) + \int_{s_i}^q ds {dA \over ds}\;\;\;\;\;\;\;\;\;\;\Delta A = \int_{s_i}^{s_f} ds {dA \over ds} \label{33}$

In the free-energy profile expression $$A (s_i)$$ is just an additive constant that can be left off. The values $${s_1,...,s_n }$$ at which the reaction coordinate is constrained can be chosen at equally-spaced intervals between $${s_i}$$ and $${s_f}$$, in which a standard numerical quadrature can $$q=f_1({\bf r})$$ be applied for evaluating the integrals in Equation \ref{33}, or they can be chosen according to a more sophisticated quadrature scheme.

We next turn to the evaluation of the derivative in Equation \ref{32}. Noting that $$P(s) = \langle \delta(f_1({\bf r})-s)\rangle$$, the derivative can be written as

${1 \over P(s)}{dP \over ds} ={C_N \over Q(N,V,T)}{\int\;d^N {\bf p} d^N {\bf r} e^{-\beta H_{(p,r)}} {\partial \over \partial s}\delta(f_1({\bf r})-s)\over \langle \delta(f_1({\bf r})-s)\rangle} \label{34}$

In order to avoid evaluating the derivative of the $$\delta$$-function, an integration by parts can be used. First, we introduce a complete set of $$3N$$ generalized coordinates:

$q_{\alpha} = f_{\alpha}({\bf r}_1,...,{\bf r}_N) \label{35}$

and their conjugate momenta $${p_{\alpha} }$$. Such a transformation has a unit Jacobian so that $$d^N{\bf p}\;d^N{\bf r}= d^{3N}p\;d^{3N}q$$. Denoting the transformed Hamiltonian as $${\tilde{H}(p,q) }$$, Equation \ref{34} becomes

${1 \over P(s)}{dP \over ds} ={C_N \over Q(N,V,T)}{\int\;d^{3N} P d^{3N} q e^{-\beta \tilde {H} (p, q)}{\partial \over \partial s}\delta(q_1-s)\over \langle \delta(q_1-s)\rangle} \label{36}$

Changing the derivative in front of the $$\delta$$-function from $$\partial/\partial s$$ to $$\partial/\partial q_1$$, which introduces an overall minus sign, and then integrating by parts yields

\begin{align} {1 \over P(s)}{dP \over ds} &= {C_N \over Q(N,V,T)}{\int\;d^{3N}p\;d^{3N}q\;\left[{\partial \over \partial q_1}e^{-\beta \tilde{H}(p,q)}\right]\delta(q_1-s)\over \langle \delta(q_1-s)\rangle} \\[4pt] &= -{\beta C_N \over Q(N,V,T)}{\int\;d^{3N}p\;d^{3N}q\;{\partial\tilde {H} \over \partial q_1}e^{-\beta \tilde{H}(p,q)}\delta(q_1-s)\over \langle \delta(q_1-s)\rangle} \\[4pt] & = -\beta {\left<\left({\partial \tilde{H} \over \partial q_1}\right)\delta(q_1-s)\right>\over \langle \delta(q_1-s)\rangle} \label{37} \end{align}

The last line defines a new ensemble average, specifically an average subject to the condition (not constraint) that the coordinate $${q_1}$$ have the particular value $$s$$. This average will be denoted $$\langle\cdots\rangle^{\rm cond}_{s}$$. Thus, the derivative becomes

${1 \over P(s)}{dP \over ds} =-\beta\left<{\partial \tilde{H} \over \partial q_1}\right>^{\rm cond}_s \label{38}$

Substituting Equation \ref{38} yields a free energy profile of the form

$A(q) = A(s_i) + \int_{s_i}^q\;ds\;\left<{\partial \tilde{H} \over \partial q_1}\right>^{\rm cond}_s \label{39}$

from which $$\Delta A$$ can be computed by letting $${q = s_f}$$. Given that - $${\langle \partial \tilde{H}/\partial q_1\rangle^{\rm cond}_s }$$ is the expression for the average of the generalized force on $${q_1}$$ when $${q_1 = s}$$, the integral represents the work done on the system, i.e. the negative of the work done by the system, in moving from $${s_i}$$ to an arbitrary final point $${q}$$. Since the conditional average implies a full simulation at each fixed value of $${q_1}$$, the thermodynamic transformation is certainly carried out reversibly, so that Equation \ref{39} is consistent with the Clausius inequality.

Although Equation \ref{39} provides a very useful insight into the underlying statistical mechanical expression for the free energy, technically, the need for a full canonical transformation of both coordinates and momenta is inconvenient since, from the chain rule

${\partial \tilde{H} \over \partial q_1} =\sum_{i=1}^N\left [ {\partial H \over \partial {\bf p} _i } \cdot {\partial {\bf p}_i \over \partial q_1} + {\partial H \over \partial {\bf r}_i }\cdot{\partial {\bf r}_i \over \partial q_1}\right] \label{40}$

A more useful expression results if the momenta integrations are performed before introducing the transformation to generalized coordinates. Starting again with Equation \ref{34}, we carry out the momentum integrations, yielding

${1 \over P(s)}{dP \over ds} ={1 \over N!\lambda^{3N}Q(N,V,T )} {\int d^N {\bf r} e^{-\beta U (r)} {\partial \over \partial s} \delta(f_1({\bf r})-s)\over \langle \delta(f_1({\bf r})-s)\rangle} \label{41}$

Now, we introduce only the transformation of the coordinates to generalized coordinates $$q_{\alpha} = f_{\alpha}({\bf r}_1,...,{\bf r}_N)$$. However, because there is no corresponding momentum transformation, the Jacobian of the transformation is not unity. Let $$J(q)\equiv J(q_1,...,q_{3N}) =\partial ({\bf r}_1,...,{\bf r}_N)/\partial (q_1,...,q_{3N})$$ denote the Jacobian of the transformation. Then, Equation \ref{41} becomes

\begin{align} {1 \over P(s)}{dP \over ds} &= { {1 \over N!\lambda^{3N}Q(N,V,T)}{\int\;d^{3N}q\;J(q)e^{-\beta {\tilde U} (q)}{\partial \over \partial s}\delta(q_1-s)\over \langle \delta(q_1-s)\rangle}} \\[4pt] &= {1 \over N!\lambda^{3N}Q(N,V,T)}{\int\;d^{3N}q\;e^{-\beta \left(U(q) - kT J(q) \right )}{\partial \over \partial s}\delta(q_1-s)\over \langle \delta(q_1-s)\rangle} \label{42} \end{align}

where, in the last line, the Jacobian has been exponentiated. Changing the derivative $$\partial/\partial s$$ to $$\partial/\partial q_1$$ and performing the integration by parts as was done in Equation \ref{37}, we obtain

\begin{align} {1 \over P(s)}{dP \over ds} &= {1 \over N!\lambda^{3N}Q(N,V,T)}{\int\;d^{3N}q\;{\partial \over \partial q_1}e^{-\beta \left (\tilde{U}(q)-kT\ln J(q)\right)}\delta(q_1-s)\over \langle \delta(q_1-s)\rangle} \\[4pt] &= -{\beta \over N!\lambda^{3N}Q(N,V,T)}{\int\;d^{3N}q\;\left[{\partial \tilde {U} \over partial q_1} -KT{\partial \over \partial q_1} \ln J (q) \right]e^{-\beta \left (\tilde{U}(q)-kT\ln J(q)\right)}\delta(q_1-s)\over \langle \delta(q_1-s)\rangle} \\[4pt] &= {-\beta\left<\left[{\partial\tilde{U} \over \partial q_1}-kT{\partial \over \partial q_1}\ln J(q)\right]\right>^{\rm cond}_s } \label{43} \end{align}

Therefore, the free energy profile becomes

$A(q) = A(s_i) +\int_{s_i}^q\;ds\;\left<\left[{\partial\tilde {U} \over \partial q_1} - KT {\partial \over \partial q_1}\ln J(q)\right]\right>^{\rm cond}_s \label{44}$

Again, the derivative of $${\tilde U}$$, the transformed potential, can be computed form the untransformed potential via the chain rule

${\partial \tilde{U} \over \partial q_1} =\sum_{i=1}^N {\partial U \over \partial {\bf r}_i}\cdot {\partial {\bf r}_i \over \partial q_1} \label{45}$

Equation \ref{44} is useful for simple reaction coordinates in which the full transformation to generalized coordinates is known. We will see shortly how the expression for $$A (q)$$ can be further simplified in a way that does not require knowledge of the transformation at all. First, however, we must tackle the problem alluded to earlier of computing the conditional ensemble averages from the constrained dynamics employed by the blue moon ensemble method.

This page titled 8.6: The "blue moon'' Ensemble Approach is shared under a CC BY-NC-SA 4.0 license and was authored, remixed, and/or curated by Mark Tuckerman.