クルックスの揺動定理からのベネット受容比法の導出

本記事では参考文献１を参考にして，クルックスの揺動定理（Crooks fluctuation theorem）と最尤法を用いてベネット受容比法（Bennett acceptance ratio method）を導出する手順を記します。

クルックスの揺動定理とは

クルックスの揺動定理は熱力学状態を変化させるために必要な仕事量と状態間の自由エネルギー差に関する関係式を与えます。
今，２つの熱力学状態$${\rm A,\ B}$$を考えることにします。系は温度$${T}$$の熱浴に接しており，初期状態は平衡状態にあるとします。このとき，$${{\rm  A \rightarrow B}}$$の変化を生じさせるために必要な仕事量が$${W}$$となる確率を$${P_{\rm A\rightarrow B}(\mathcal{W})}$$，逆過程である$${\rm B \rightarrow A}$$に必要となる仕事量が$${\mathcal{W}}$$となる確率を$${P_{\rm B \rightarrow A}(\mathcal{W})}$$とすると，２つの確率の比に関して以下の等式が成立します。

$$
\ln\left[\frac{P_{\rm A \rightarrow B}(\mathcal{W})}{P_{\rm B \rightarrow A}(-\mathcal{W})}\right]=\beta(\mathcal{W}-\Delta F_{\rm AB})
$$

ここで，$${\beta}$$は逆温度（$${=1/T}$$），$${\Delta F_{\rm AB}}$$は状態$${\rm A, B}$$間の自由エネルギー差です。つまり，$${P_{\rm A\rightarrow B}(\mathcal{W})，P_{\rm A\rightarrow B}(-\mathcal{W})}$$が何かしらの方法で分かれば$${\Delta F_{\rm AB}}$$が求まることになります。

尤度の導出

最尤法を適用するため，クルックスの揺動定理から尤度の表式を導出します。
仕事量$${\mathcal{W}}$$は生成分布$${P(\mathcal{W})}$$からサンプリングされるものとみなし，また$${\rm A \rightarrow B}$$の際には系にする仕事，逆過程である$${\rm B \rightarrow A}$$の際は系がする仕事を表すものとします。
上記設定の下，$${P_{\rm A\rightarrow B}(\mathcal{W})}$$，$${P_{\rm B\rightarrow A}(-\mathcal{W})}$$を事後確率の表記に変更します。

$$
\begin{array}{}P_{\rm A\rightarrow B}(\mathcal{W})&\Longrightarrow& P(\mathcal{W}|{\rm A \rightarrow B})\\P_{\rm B\rightarrow A}(-\mathcal{W})&\Longrightarrow& P(\mathcal{W}|{\rm B \rightarrow A})\end{array}
$$

今は，$${\rm A\rightarrow B\ or\ B\rightarrow A}$$の試行のみが生じる状況を考えているため，

$$
P(\mathcal{W})=P(\mathcal{W},{\rm A \rightarrow B})+P(\mathcal{W},{\rm B\rightarrow A})
$$

となります。これを式変形して，

$$
\begin{array}{}1&=&\frac{P(\mathcal{W},A\rightarrow B)}{P(\mathcal{W})}+\frac{P(\mathcal{W},B\rightarrow A)}{P(\mathcal{W})}\\1&=&P({\rm A\rightarrow B}|\mathcal{W})+P({\rm B\rightarrow A}|\mathcal{W})\end{array}
$$

が得られます。
ベイズの定理に従って，クルックスの揺動定理の式を事後確率から尤度での表記に変更します。

$$
\begin{array}{}\frac{P(\mathcal{W}|{\rm A\rightarrow B})}{P(\mathcal{W}|{\rm B\rightarrow A})}&=&\left\{\frac{P(\mathcal{W}, {\rm A\rightarrow B})}{P({\rm A\rightarrow B})}\right\}\left\{\frac{P({\rm B\rightarrow A})}{P(\mathcal{W}, {\rm B\rightarrow A})}\right\}\\&=&\left\{\frac{P({\rm A\rightarrow B}|\mathcal{W})P(\mathcal{W})}{P({\rm A\rightarrow B})}\right\}\left\{\frac{P({\rm B\rightarrow A})}{P({\rm B\rightarrow A}|\mathcal{W})P(\mathcal{W})}\right\}\\&=&\left\{\frac{P({\rm B\rightarrow A})}{P({\rm A\rightarrow B})}\right\}\left\{\frac{P({\rm A\rightarrow B}|\mathcal{W})}{P({\rm B\rightarrow A}|\mathcal{W})}\right\}\\&=&{\rm e}^{-\beta M}\left\{\frac{P({\rm A\rightarrow B}|\mathcal{W})}{P({\rm B\rightarrow A}|\mathcal{W})}\right\}\ \ \left({\rm e}^{-\beta M}:=\frac{P({\rm B\rightarrow A})}{P({\rm A\rightarrow B})}\right)\\&=&{\rm e}^{\beta(\mathcal{W}-\Delta F_{\rm AB})}\end{array}
$$

$${P({\rm A\rightarrow B)}，P({\rm A\rightarrow B)}}$$はそれぞれの過程が生じる確率を表します。
２つの尤度の関係式を用いて片方を消去すると

$$
\begin{array}{}{\rm e}^{-\beta M}\left\{\frac{P({\rm A\rightarrow B|\mathcal{W})}}{1-P({\rm A\rightarrow B|\mathcal{W})}}\right\}&=&{\rm e}^{\beta(\mathcal{W}-\Delta F_{\rm AB})}\\P({\rm A\rightarrow B|\mathcal{W})}&=&\frac{1}{1+{\rm e}^{-\beta(M+\mathcal{W}-\Delta F_{\rm AB})}}\end{array}
$$

が得られます。逆過程については，

$$
\begin{array}{}P({\rm B\rightarrow A|\mathcal{W})}&=&\frac{1}{1+{\rm e}^{\beta(M+\mathcal{W}-\Delta F_{\rm AB})}}\end{array}
$$

となります。

最尤法の適用

現段階で自由エネルギー差を変数とした確率の表式を得ることができました。これらの表式と最尤法を用いて自由エネルギー差を最尤値として求めることを考えます。
今，$${\rm A\rightarrow B}$$過程を$${n_{\rm F}}$$回，$${\rm B\rightarrow A}$$過程を$${n_{\rm R}}$$回，測定orシミュレーションした結果があるとし，尤度関数を

$$
L(\Delta F_{\rm AB})=\left\{\prod_{i=1}^{n_{\rm F}}P({\rm A\rightarrow B}|\mathcal{W}_i)\right\}\left\{\prod_{j=1}^{n_{\rm R}}P({\rm B\rightarrow A}|\mathcal{W}_j)\right\}
$$

と定義します。両辺の対数をとると，

$$
\begin{array}{}\ln \{L(\Delta F_{\rm AB})\}&=&\ln\left[\left\{\prod_{i=1}^{n_{\rm F}}P({\rm A\rightarrow B}|\mathcal{W}_i)\right\}\left\{\prod_{j=1}^{n_{\rm R}}P({\rm B\rightarrow A}|\mathcal{W}_j)\right\}\right]\\&=&\sum_{i=1}^{n_{\rm F}}\ln\{P({\rm A\rightarrow B}|\mathcal{W}_i)\}+\sum_{j=1}^{n_{\rm R}}\ln\{P({\rm B\rightarrow A}|\mathcal{W}_j)\}\end{array}
$$

と変形されます。$${\Delta F_{\rm AB}}$$に関する微分が0になる条件は，

$$
\begin{array}{}\frac{\partial\ln\left\{L(\Delta F_{\rm AB})\right\}}{\partial\Delta F_{\rm AB}}&=&\sum_{i=1}^{n_{\rm F}}\frac{\partial\ln\{P({\rm A\rightarrow B}|\mathcal{W}_i)\}}{\partial\Delta F_{\rm AB}}+\sum_{j=1}^{n_{\rm R}}\frac{\partial\ln\{P({\rm B\rightarrow A}|\mathcal{W}_j)\}}{\partial\Delta F_{\rm AB}}\\&=&-\beta\left[\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta(M+\mathcal{W}_i-\Delta F_{\rm AB})}}-\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{-\beta(M+\mathcal{W}_j-\Delta F_{\rm AB})}}\right]\\&=&0\end{array}
$$

となります。右辺の$${[\ ]}$$の内部に着目すると，$${\Delta F_{\rm AB}}$$に対して単調増加する関数形となっていますので，微分が0となる$${\Delta F_{\rm AB}}$$の存在は１つです。
$${P({\rm A\rightarrow B)}，P({\rm A\rightarrow B)}}$$をそれぞれの過程のデータ数で表すと，

$$
{\rm e}^{-\beta M}=\frac{P({\rm B\rightarrow A})}{P({\rm A\rightarrow B})}=\frac{n_{\rm R}}{n_{\rm F}}
$$

となります。
さらに，$${\rm B \rightarrow A}$$の仕事量を符号を反転させて系にする仕事に戻すと，尤度関数の微分が０となる条件式は

$$
\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta\mathcal{W}_i-(\beta\Delta F_{\rm AB}-\ln(n_{\rm F}/n_{\rm R}))}}-\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{\beta\mathcal{W}_j+(\beta\Delta F_{\rm AB}-\ln(n_{\rm F}/n_{\rm R}))}}=0
$$

となります。
ここで，

$$
C=\Delta F_{\rm AB}-\frac{1}{\beta}\ln\left(\frac{n_{\rm F}}{n_{\rm R}}\right) 
$$

を導入すると，

$$
\begin{array}{}\beta \Delta F_{\rm AB}&=&C-\frac{1}{\beta}\ln\left(\frac{n_{\rm R}}{n_{\rm F}}\right)\\&=&0+C-\frac{1}{\beta}\ln\left(\frac{n_{\rm R}}{n_{\rm F}}\right)\\&=&\ln\left[\frac{\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_j+C)}}}{\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_i-C)}}}\right]+C-\frac{1}{\beta}\ln\left(\frac{n_{\rm R}}{n_{\rm F}}\right)\\&=&\ln\left[\frac{\sum_{j=1}^{n_{\rm R}}f(\beta(\mathcal{W}_j+C))}{\sum_{i=1}^{n_{\rm F}}f(\beta(\mathcal{W}_i-C))}\right]+C-\frac{1}{\beta}\ln\left(\frac{n_{\rm R}}{n_{\rm F}}\right)\ \ \left(f(x):=\frac{1}{1+{\rm e}^{x}}\right)\end{array}
$$

となり，参考文献２の(12a), (12b)式が得られます。
以上より，目的のベネット受容比法の表式を導出することができました。

ジャルジンスキー等式との関係

データ数が$${\rm A \rightarrow B}$$過程，もしくは$${\rm B \rightarrow A}$$過程に偏っている場合を考えます。
$${n_{\rm F} \gg n_{\rm R}(-C \gg 1)}$$の場合，

$$
\begin{array}{}\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_i-C)}}&\simeq&\sum_{i=1}^{n_{\rm F}}{\rm e}^{-\beta(\mathcal{W}_i-C)}\\&=&n_{\rm F}{\rm e}^{\beta C}\langle{\rm e}^{-\beta\mathcal{W}}\rangle_{\rm F}\\\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_j+C)}}&\simeq&\sum_{j=1}^{n_{\rm R}}\frac{1}{1+0}\\&=&n_{\rm R}\\\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_i-C)}}-\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_j+C)}}&\simeq&n_{\rm F}{\rm e}^{\beta C}\langle{\rm e}^{-\beta\mathcal{W}}\rangle_{\rm F}-n_{\rm R}\\&=&0\\\langle{\rm e}^{-\beta\mathcal{W}}\rangle_{\rm F}&=&{\rm e}^{-\beta C}\frac{n_{\rm R}}{n_{\rm F}}\\&=&{\rm e}^{-\beta\Delta F_{\rm AB}}\end{array}
$$

となり，ジャルジンスキー等式（Jarzynski Equality）が得られます。
一方，$${n_{\rm R} \gg n_{\rm F}(C \gg 1)}$$の場合，

$$
\begin{array}{}\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_i-C)}}&\simeq&\sum_{i=1}^{n_{\rm F}}\frac{1}{1+0}\\&=&n_{\rm F}\\\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_j+C)}}&\simeq&\sum_{j=1}^{n_{\rm R}}{\rm e}^{-\beta(\mathcal{W}_j+C)}\\&=&n_{\rm R}{\rm e}^{-\beta C}\langle{\rm e}^{-\beta \mathcal{W}}\rangle_{\rm R}\\\sum_{i=1}^{n_{\rm F}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_i-C)}}-\sum_{j=1}^{n_{\rm R}}\frac{1}{1+{\rm e}^{\beta(\mathcal{W}_j+C)}}&\simeq&n_{\rm F}-n_{\rm R}{\rm e}^{-\beta C}\langle{\rm e}^{-\beta \mathcal{W}}\rangle_{\rm R}\\&=&0\\\langle{\rm e}^{-\beta\mathcal{W}}\rangle_{\rm R}&=&{\rm e}^{\beta C}\frac{n_{\rm F}}{n_{\rm R}}\\&=&{\rm e}^{\beta\Delta F_{\rm AB}}\\&=&{\rm e}^{-\beta\Delta F_{\rm BA}}\end{array}
$$

となり，やはりジャルジンスキー等式が得られます。
つまり，ジャルジンスキー等式は片方の過程にデータが偏った場合のベネット受容比法であると解釈することができます。

参考文献

1. Phys. Rev. Lett. 91, 140601 (2003).

2. J. Comput. Phys., 22, 245-268.