Metadynamicsの境界補正について
本記事では参考文献１で提案された，Metadynamicsの境界補正方法（McGovern−de Pablo boundary conditions）について解説します。
Metadynamics自体について知りたい方は，以下の記事のAppendix Aを併せてご参照ください。
問題設定
系の配置に関する自由度をベクトル$${\mathbf{q}=\begin{bmatrix}q_1&q_2&\dots\end{bmatrix}^{\rm T}}$$とし，$${N}$$個の集団座標$${\mathbf{s}=\begin{bmatrix}s_1&s_2&\dots&s_N\end{bmatrix}^{\rm T}}$$を対象としたmetadynamicsを実行することを考えます。
各$${s_i}$$が取る値の範囲はそれぞれ上下限が存在すると仮定します。
$$
\begin{align*}
L_i&\le s_i \le U_i
\end{align*}\tag{1}
$$
各時刻における$${\mathbf{s}(t)}$$の値は$${\mathbf{q}(t)}$$の関数として計算されます。$${\mathbf{q}(t)}$$の関数として扱う場合の集団座標を$${\hat{\mathbf{s}}(\mathbf{q}(t))}$$と表記することにします。
このとき，時刻$${t}$$における集団座標$${\mathbf{s}(t)}$$は
$$
\begin{align*}
\mathbf{s}(t)&=\hat{\mathbf{s}}(\mathbf{q}(t))
\end{align*}\tag{2}
$$
となります。
時間間隔$${\tau}$$毎にガウス型のバイアスポテンシャル$${V_{\rm G}(\mathbf{s})=W\exp(-\sum_{i=1}^N\frac{(s_i-\hat{s}_i(\mathbf{q}(t)))^2}{2\sigma_i^2})}$$を課していきます。
時刻$${t}$$におけるバイアスポテンシャル$${V(\mathbf{s},t)}$$は
$$
\begin{align*}
V(\mathbf{s},t)&=\sum_{k\tau<t}W\exp\left(-\sum_{i=1}^N\frac{(s_i-\hat{s}_i(\mathbf{q}(k\tau)))^2}{2\sigma_i^2}\right)
\end{align*}\tag{3}
$$
となります。
問題点
十分にmetadynamicsが進んだ結果，$${\mathbf{s}}$$の確率分布$${P(\mathbf{s})}$$が収束したとします。
この状態で更に$${V_{\rm G}}$$を追加した場合のバイアスポテンシャルの平均変化率を$${\overline{\Delta V}(\mathbf{s})}$$とすると，
$$
\begin{align*}
\overline{\Delta V}(\mathbf{s})&=W\int{\rm d}\mathbf{x}P(\mathbf{x})\exp\left(-\sum_{i=1}^N\frac{(s_i-x_i)^2}{2\sigma_i^2}\right)\\
\mathbf{x}&:=\mathbf{q}(t)
\end{align*}\tag{4}
$$
と計算されます。
(Non-tempered) metadynamicsの場合，$${P(\mathbf{s})}$$は一様分布（$${\mathbf{s}}$$に依存しない定数）に収束することが期待されます。$${\mathbf{s}}$$に境界条件が課せられていることを思い出すと，
$$
\begin{align*}
P(\mathbf{s})&=\prod_{i=1}^N\frac{1}{U_i-L_i}
\end{align*}\tag{5}
$$
となっているはずです。
(5)式を(4)式に代入することにより，$${\overline{\Delta V}(\mathbf{s})}$$をより具体的に評価することができます。
$$
\begin{align*}
\overline{\Delta V}(\mathbf{s})&=W\prod_{i=1}^N\frac{1}{U_i-L_i}\int_{L_i}^{U_i}{\rm d}x_i\exp\left(-\frac{(s_i-x_i)^2}{2\sigma_i^2}\right)\\
&=W\prod_{i=1}^N\frac{1}{U_i-L_i}\int_{\frac{L_i-s_i}{\sqrt{2}\sigma_i}}^{\frac{U_i-s_i}{\sqrt{2}\sigma_i}}{\rm d}x_i\exp\left(-x_i^2\right)\\
&=W\sqrt{\frac{\pi}{2}}\prod_{i=1}^N\frac{\sigma_i}{U_i-L_i}\left[{\rm erf}\left(\frac{U_i-s_i}{\sqrt{2}\sigma_i}\right)+{\rm erf}\left(\frac{s_i-L_i}{\sqrt{2}\sigma_i}\right)\right]\\
\end{align*}\tag{6}
$$
(6)式が意味するところは，収束後のバイアスポテンシャルの変化は$${\mathbf{s}}$$に依存するということです。$${\mathbf{s}}$$に依存するということはバイアスポテンシャルの形状が依然として変化することを表すのですが，これは$${P(\mathbf{s})}$$が一様分布に収束したという仮定と矛盾します。
つまり，上下限の境界条件を有する集団座標のmetadynamicsは厳密には一様分布に収束しないということになります。
一様分布に収束するという仮定は，自由エネルギー曲面（平均力ポテンシャル）を得るという目的に対して重要な意味を持つため，$${P(\mathbf{s})}$$が一様分布に収束するよう何かしらの補正が求められます。
補正方法
バイアスポテンシャルに追加していく$${V_{\rm G}}$$を
$$
\begin{align*}
\overline{\Delta V}(\mathbf{s})&=\int{\rm d}\mathbf{x}P(\mathbf{x})V_{\rm M}(\mathbf{s}, \mathbf{x})&=C\\
\end{align*}\tag{7}
$$
を満たす$${V_{\rm M}}$$に置き換えることを考えます。
参考文献１では，大別して掛け算補正と足し算補正の２種類が紹介されています。下記はあくまで一例であって，(7)式を満たす補正は他にも存在しえます。
掛け算補正
$$
\begin{align*}
V_{\rm M}(\mathbf{s}, \mathbf{x})&=C(\mathbf{s})V_{\rm G}(\mathbf{s}, \mathbf{x})\\
\end{align*}\tag{7}
$$
とすると，
$$
\begin{align*}
C(\mathbf{s})&=\frac{C}{\sqrt{\frac{\pi}{2}}\prod_{i=1}^N\frac{\sigma_i}{U_i-L_i}\left[{\rm erf}\left(\frac{U_i-s_i}{\sqrt{2}\sigma_i}\right)+{\rm erf}\left(\frac{s_i-L_i}{\sqrt{2}\sigma_i}\right)\right]}
\end{align*}\tag{8}
$$
とすると(7)式が満たされることが容易に示されます。
足し算補正
$$
\begin{align*}
V_{\rm M}(\mathbf{s}, \mathbf{x})&=V_{\rm G}(\mathbf{s}, \mathbf{x})+C(\mathbf{s})\\
\end{align*}\tag{9}
$$
とすると，
$$
\begin{align*}
C(\mathbf{s})&=-W\sqrt{\frac{\pi}{2}}\prod_{i=1}^N\frac{\sigma_i}{U_i-L_i}\left[{\rm erf}\left(\frac{U_i-s_i}{\sqrt{2}\sigma_i}\right)+{\rm erf}\left(\frac{s_i-L_i}{\sqrt{2}\sigma_i}\right)\right]+C
\end{align*}\tag{10}
$$
とすると(7)式が満たされることが容易に示されます。
元も子もない話。。。
例として，集団座標が1つで$${U=1,\ L=-1,\ \sigma=0.02, W=1}$$の場合を考えます。
このとき，(6)式を描画すると図1のようになります。
図1から分かる通り，境界付近を除けば期待通り一様分布が実現されることになります。
実践的には境界近傍の値に注目せざるを得ない状況というのは稀な気がしますので，そこまで頑張って補正しなくても問題ないかもしれません。
参考文献
この記事が気に入ったらサポートをしてみませんか？