metadynamicsにおける自由エネルギーの時間依存定数項を決定する方法

well-tempered metadynamicsが十分に収束すると，時間依存バイアスポテンシャル$${V(s,t)}$$と自由エネルギープロファイル（平均力ポテンシャル）$${F(s)}$$には

$$
\begin{align*}
F(s) &=-\left[1+\frac{T}{\Delta T}\right]V(s,t)+C(t)
\end{align*}\tag{1}
$$

の関係が成立します。ここで，$${\Delta T}$$はwell-tempered metadynamicsにおけるスケーリングを調整するパラメータです。
本記事では，時間変化はするが，$${s}$$には依存しない定数項$${C(t)}$$が(2)式で表すことができることを解説します。

$$
\begin{align*}
C(t) &=k_{\rm B}T\ln\left(\int{\rm d}s\exp\left[\left(1+\frac{T}{\Delta T}\right)\frac{V(s,t)}{k_{\rm B}T}\right]\right)
\end{align*}\tag{2}
$$

ここで，$${k_{\rm B}}$$はボルツマン定数，$${T}$$は系の温度を表します。
(1)-(2)式が意味するところは，時間依存する$${V(s,t)}$$から時間依存しない$${F(s)}$$を完全に決定できることです。

Appendix Aにmetadynamicsの簡単な概説を記しましたので，そもそもmetadynamicsとは何ぞや？という方はよろしければAppendix Aの内容も併せてご参照ください。

Quasi-stationary limit

目的の導出にあたり，重要な仮定としてquasi-stationary limitがあります。quasi-stationary limitとは，バイアスポテンシャルの時間変化が系の平衡化時間よりも十分に遅くなった状態となります。well-tempered metadynamicsは付加するガウス関数の高さを徐々に0へとスケーリングしていきますので，十分に時間が経過すればいずれquasi-stationary limitに到達します。
quasi-stationary limitに到達すると，系の自由度の出現確率$${P(s,t)}$$がボルツマン分布に従うと仮定できます。

$$
\begin{align*}
P(s,t) &=\frac{\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}
\end{align*}\tag{3}
$$

バイアスポテンシャルの時間変化

(3)式の両辺を$${t}$$で微分すると，

$$
\begin{align*}
\dot{P}(s,t) &=-\frac{1}{k_{\rm B}T}\left\{\dot{V}(s,t)-\left\langle\dot{V}\right\rangle_{V_t}\right\}P(s,t)
\end{align*}\tag{4}
$$

となります。ここで，$${\langle\cdots\rangle_{V_t}}$$は$${P(s,t)}$$を用いた期待値を表します。
$${V(s,t)}$$の時間微分は(A1)式より，

$$
\begin{align*}
\dot{V}(s,t)&\simeq\frac{V(s,t)-V(s,t-\Delta t)}{\Delta t}\\
&=\frac{W(t)}{\Delta t}\exp\left[-\frac{\left\{s-\hat{s}(q(t))\right\}^2}{2\sigma^2}\right]
\end{align*}\tag{5}
$$

となります。
ガウス関数はデルタ関数に近似できる程度に$${\sigma}$$が小さいと仮定します。

$$
\begin{align*}
\dot{V}(s,t)&\simeq\frac{W(t)}{\Delta t}\sqrt{2\pi}\sigma\delta(s-\hat{s}(q(t)))
\end{align*}\tag{6}
$$

(6)式を(4)式の期待値の部分に挿入することにより，

$$
\begin{align*}
\dot{P}(s,t) &\simeq-\frac{1}{k_{\rm B}T}\left\{\dot{V}(s,t)-\frac{W(t)}{\Delta t}\sqrt{2\pi}\sigma\left\langle\delta(s-\hat{s}(q(t)))\right\rangle_{V_t}\right\}P(s,t)\\
&=-\frac{1}{k_{\rm B}T}\left\{\dot{V}(s,t)-\frac{W(t)}{\Delta t}\sqrt{2\pi}\sigma P(\hat{s}(q(t)),t)\right\}P(s,t)\\
\end{align*}\tag{7}
$$

となります。
quasi-stationary limitでは，$${P(s,t)}$$の時間変化も小さいことが期待されます。
任意の$${s}$$に対して$${\dot{P}(s,t)/P(s,t)\ll 1}$$を仮定できるためには，$${s=\hat{s}(q(t))}$$に対して

$$
\begin{align*}
\dot{V}(s=\hat{s}(q(t)),t) &\simeq\frac{W(t)}{\Delta t}\sqrt{2\pi}\sigma P(\hat{s}(q(t)),t)\\
&=\frac{\sqrt{2\pi}\sigma W_0}{\Delta t}\exp\left[-\frac{V(\hat{s}(q(t)), t)}{k_{\rm B}\Delta T}\right]P(\hat{s}(q(t)),t)
\end{align*}\tag{8}
$$

が成立する必要があります。
(8)式の$${\hat{s}(q(t))}$$を改めて$${s}$$に書き換えると，

$$
\begin{align*}
\dot{V}(s,t) &=\frac{\sqrt{2\pi}\sigma W_0}{\Delta t}\exp\left[-\frac{V(s, t)}{k_{\rm B}\Delta T}\right]P(s,t)\\
&=\frac{\sqrt{2\pi}\sigma W_0}{\Delta t}\exp\left[-\frac{V(s, t)}{k_{\rm B}\Delta T}\right]\frac{\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}\\
\dot{V}(s,t)\exp\left[\left(1+\frac{\Delta T}{T}\right)\frac{V(s,t)}{k_{\rm B}\Delta T}\right]&=\frac{\sqrt{2\pi}\sigma W_0}{\Delta t}\frac{\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}\\
&=\frac{\sqrt{2\pi}\sigma W_0}{\Delta t}P_0(s)\frac{\int{\rm d}s\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}\\
\therefore \frac{{\rm d}}{{\rm d}t}\exp\left[\left(1+\frac{\Delta T}{T}\right)\frac{V(s,t)}{k_{\rm B}\Delta T}\right]&=\left(1+\frac{\Delta T}{T}\right)\frac{\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}P_0(s)\frac{\int{\rm d}s\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}\\
\end{align*}\tag{9}
$$

が得られます。ここで，$${P_0(s)}$$はバイアスポテンシャルなしの場合の自由度$${s}$$の確率分布です。
(9)式は見た目が複雑なので，

$$
\begin{align*}
\gamma&:=1+\frac{\Delta T}{T}\\
{\rm e}^{\frac{c(t)}{k_{\rm B}T}}&:=\frac{\int{\rm d}s\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}
\end{align*}\tag{10}
$$

と定義し，

$$
\begin{align*}
\frac{{\rm d}}{{\rm d}t}\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]&=\frac{\gamma\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}P_0(s){\rm e}^{\frac{c(t)}{k_{\rm B}T}}\\
\end{align*}\tag{11}
$$

と表すことにします。

平均力ポテンシャルの表式の導出

(11)式を$${t}$$と$${s}$$に関して積分することより，$${V(s,t)}$$を用いた$${F(s)}$$の表式が得られます。
まず，(11)式を$${t}$$に関して積分すると，

$$
\begin{align*}
\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]-1&=\frac{\gamma\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}P_0(s)\int_{0}^t{\rm d}t'{\rm e}^{\frac{c(t')}{k_{\rm B}T}}\\
&=\frac{\gamma\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}P_0(s)\tau\\
\tau&:=\int_{0}^t{\rm d}t'{\rm e}^{\frac{c(t')}{k_{\rm B}T}}
\end{align*}\tag{12}
$$

となります。
(12)式の左辺の1は計算がある程度進むと無視できる大きさと見なせるため，

$$
\begin{align*}
\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]&\simeq\frac{\gamma\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}P_0(s)\tau
\end{align*}\tag{13}
$$

が成立します。
(13)式を$${s}$$に関して積分すると，

$$
\begin{align*}
\int{\rm d}s\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]&=\frac{\gamma\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}\tau\int{\rm d}sP_0(s)\\
&=\frac{\gamma\sqrt{2\pi}\sigma W_0}{\Delta t k_{\rm B}\Delta T}\tau\\
\therefore \tau&=\frac{\Delta t k_{\rm B}\Delta T}{\gamma\sqrt{2\pi}\sigma W_0}\int{\rm d}s\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]
\end{align*}\tag{14}
$$

と$${\tau}$$を$${V(s,t)}$$を積分する形で表現できます。
(14)式を(13)式に代入すると，

$$
\begin{align*}
\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]&=P_0(s)\int{\rm d}s\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]\\
P_0(s)&=\frac{\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]}{\int{\rm d}s\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]}\\
\frac{\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}&=\frac{\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]}{\int{\rm d}s\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]}\\
\therefore F(s)&=-\frac{T\gamma V(s,t)}{\Delta T}+k_{\rm B}T\int{\rm d}s\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]\\
&\ \ \ \ -k_{\rm B}T\int{\rm d}s\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]
\end{align*}\tag{15}
$$

(15)式の右辺の第三項は$${s}$$にも$${t}$$にも依存しない定数なので，無視することができます。さらに，$${\gamma=1+\Delta T/T}$$であることを思い出すと，(15)式は(1)-(2)式に帰着します。
以上より，目的の式を導出することができました。

今回の定式化の利点

参考文献１では，今回の定式化が単に理論的な興味を満たすだけでなく，実用的にも利点があることを主張しています。

１つめの利点は，metadynamics収束の解析の際に異なるmetadynamicsシミュレーションで得られた平均力ポテンシャルを直接比較できることが挙げられます。
$${C(t)}$$が未定のままでも適当にalignmentすることによって，異なるシミュレーション結果を比較することは原理的に可能です。ただし，その場合，何を以って「適当にalignment」とするか，という問題が生じます。逆にいうと，今回の定式化は「理論的に正当なalignment」の基準を与えると解釈できます。

もう一つの利点は，定式化に至るまでの式変形において，配置$${q}$$に関するre-weighting factorを見積もるための式の新しい定式化もできることが挙げられます。
metadynamicsで得られたtrajectoryはバイアスポテンシャルの影響があるため，そのままでは配置$${q}$$に関する解析に利用することはできません。
ですが，実は(10)式で定義した$${c(t)}$$（$${C(t)}$$ではないことに注意）とバイアスポテンシャルを用いて，$${\exp[\frac{V(\hat{s}(q(t)),t)-c(t)}{k_{\rm B}T}]}$$でre-weightingすることによって，通常のMD trajectoryと同じように配置$${q}$$に関する解析（＝つまり，metadynamicsで対象とした自由度$${s}$$以外の自由度も解析対象とできる）ができるようになります。なぜ$${\exp[\frac{V(\hat{s}(q(t)),t)-c(t)}{k_{\rm B}T}]}$$がre-weighting factorとなるかはAppendix Bで解説しています。

$${{\rm e}^{\frac{c(t)}{k_{\rm B}T}}}$$の定式は，(14)式の両辺を$${t}$$で微分することで得られます。

$$
\begin{align*}
{\rm e}^{\frac{c(t)}{k_{\rm B}T}}&=\frac{{\rm d}\tau}{{\rm d}t}\\
&=\frac{\Delta t}{\sqrt{2\pi}\sigma W_0}\int{\rm d}s\dot{V}(s,t)\exp\left[\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]\\
&=\frac{1}{\sqrt{2\pi}\sigma}\int{\rm d}s\exp\left[-\frac{(s-\hat{s}(q(t)))^2}{2\sigma^2}+\frac{\gamma V(s,t)}{k_{\rm B}\Delta T}\right]\\
\end{align*}\tag{16}
$$

参考文献

1. Pratyush Tiwary and Michele Parrinello, J. Phys. Chem. B 2015, 119, 736−742

Appendix A: metadynamicsの概説

metadynamicsは，系の自由度に対してバイアスポテンシャルをon-the-flyで付加していくことによって，その自由度に対する自由エネルギープロファイル（平均力ポテンシャル）を計算するための手法です。原子間距離，二面角，RMSDといった位置座標の関数として表される自由度（複数の座標で定義される自由度なので集団座標，Collective Variableと称することが一般的です）がmetadynamicsの対象となります。

on-the-flyで付加されるバイアスポテンシャルの関数形はガウス関数を使用することが一般的であり，(A1)式にようになります。

$$
\begin{align*}
V(s,t) &=\sum_{k\Delta t\leq t}W(k\Delta t)\exp\left[-\frac{\left\{s-\hat{s}(q(k\Delta t))\right\}^2}{2\sigma^2}\right]
\end{align*}\tag{A1}
$$

ここで，$${s}$$はmetadynamicsの対象となる自由度，$${q(k\Delta t)}$$は$${t=k\Delta t}$$における系の位置座標（例えば，3次元空間に$${N}$$粒子ある場合は$${q=\begin{bmatrix}\mathbf{r}_1&\cdots&\mathbf{r}_N\end{bmatrix}, \mathbf{r}_i=\begin{bmatrix}x_i&y_i&z_i\end{bmatrix}^{\rm T}}$$を表します），$${\Delta t}$$はバイアスポテンシャルを更新する時間間隔，$${\sigma^2}$$はガウス分布の分散，及び$${W(k\Delta t)}$$はガウス関数の高さになります。通常のmetadynamicsでは$${W(k\Delta t)=W_0}$$（＝時間に依存しない定数）ですが，well-tempered metadynamicsと称される手法では滑らかに収束させる目的で(A2)式を採用します。

$$
\begin{align*}
W(k\Delta t)&=W_0\exp\left[-\frac{V(\hat{s}(q(k\Delta t)), k\Delta t)}{k_{\rm B}\Delta T}\right]
\end{align*}\tag{A2}
$$

$${\Delta T\rightarrow \infty}$$で通常のmetadynamicsに収束します。

(A1)式では対象となる自由度が１つしかない場合になりますが，原理的には複数の自由度を設定することが可能です。ただし，自由度空間にくまなくガウス関数を付加することを鑑みると，現実的には３自由度以上に対して収束させることは殆どの場合において困難です。そのため，１つのmetadynamicsで取り扱う自由度数は1～２に調整されることが多い印象があります。

well-tempered metadynamicsが十分に収束した後には，自由度の平均力ポテンシャル（自由エネルギープロファイル）$${F(s)}$$とバイアスポテンシャルには(A3)式の関係式が成立します。

$$
\begin{align*}
F(s)&=-\left(1+\frac{T}{\Delta T}\right)V(s,t)+C(t)
\end{align*}\tag{A3}
$$

Appendix B: バイアポテンシャル有無の配置qの確率分布の関係

系のポテンシャルを$${U(q)}$$とすると，quasi-stationary limitに達したmetadynamicsにおける配置$${q}$$の確率分布$${P(q,t)}$$は$${U(q)+V(\hat{s}(q),t)}$$のボルツマン分布となります。

$$
\begin{align*}
P(q,t) &=\frac{\exp\left[-\frac{U(q)+V(\hat{s}(q),t)}{k_{\rm B}T}\right]}{\int{\rm d}q\exp\left[-\frac{U(q)+V(\hat{s}(q),t)}{k_{\rm B}T}\right]}
\end{align*}\tag{B1}
$$

(B1)式右辺の分母にデルタ関数を用いて$${s}$$に関する積分を挿入すると，

$$
\begin{align*}
&\ \int{\rm d}q\exp\left[-\frac{U(q)+V(\hat{s}(q),t)}{k_{\rm B}T}\right]\\
&=\int{\rm d}s\int{\rm d}q\delta(\hat{s}(q)-s)\exp\left[-\frac{U(q)+V(\hat{s}(q),t)}{k_{\rm B}T}\right]\\
&=\int{\rm d}s\exp\left[-\frac{V(s,t)}{k_{\rm B}T}\right]\int{\rm d}q\delta(\hat{s}(q)-s)\exp\left[-\frac{U(q)}{k_{\rm B}T}\right]\\
&=\int{\rm d}s\exp\left[-\frac{V(s,t)}{k_{\rm B}T}\right]\frac{\int{\rm d}q\delta(\hat{s}(q)-s)\exp\left[-\frac{U(q)}{k_{\rm B}T}\right]}{\int{\rm d}q\exp\left[-\frac{U(q)}{k_{\rm B}T}\right]}\int{\rm d}q\exp\left[-\frac{U(q)}{k_{\rm B}T}\right]\\
&=\int{\rm d}s\exp\left[-\frac{V(s,t)}{k_{\rm B}T}\right]P_0(s)\int{\rm d}q\exp\left[-\frac{U(q)}{k_{\rm B}T}\right]\\
\end{align*}\tag{B2}
$$

と式変形できます。
(B2)式を(B1)式に代入すると，

$$
\begin{align*}
{\rm (r.h.s)} &=\frac{\exp\left[-\frac{U(q)+V(\hat{s}(q),t)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{V(s,t)}{k_{\rm B}T}\right]P_0(s)\int{\rm d}q\exp\left[-\frac{U(q)}{k_{\rm B}T}\right]}\\
&=\frac{P_0(q)\exp\left[-\frac{V(\hat{s}(q),t)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{V(s,t)}{k_{\rm B}T}\right]P_0(s)}\\
\end{align*}\tag{B3}
$$

となります。ここで，$${P_0(q)}$$はバイアスポテンシャルがない場合の配置$${q}$$の確率分布です。
(B3)式を$${P_0(q)}$$についてまとめると，最終的に

$$
\begin{align*}
P_0(q)&={\exp\left[\frac{V(\hat{s}(q),t)}{k_{\rm B}T}\right]}\int{\rm d}s\exp\left[-\frac{V(s,t)}{k_{\rm B}T}\right]P_0(s)P(q,t)\\
&=\exp\left[\frac{V(\hat{s}(q),t)}{k_{\rm B}T}\right]\frac{\int{\rm d}s\exp\left[-\frac{F(s)+V(s,t)}{k_{\rm B}T}\right]}{\int{\rm d}s\exp\left[-\frac{F(s)}{k_{\rm B}T}\right]}P(q,t)\\
&=\exp\left[\frac{V(\hat{s}(q),t)-c(t)}{k_{\rm B}T}\right]P(q,t)\\
\end{align*}\tag{B4}
$$

が得られます。
$${P(q,t)}$$の情報はquasi-stationary limitに達した後にmetadynamics trajectoryをサンプリングすることで得ることができます。(B4)式はサンプリングされた各スナップショット（フレーム）に$${\exp\left[\frac{V(\hat{s}(q),t)-c(t)}{k_{\rm B}T}\right]}$$をかけることにより，$${P_0(q)}$$の情報に変換できることを意味しています。