見出し画像

Parallel Bias Metadynamics

nori

本記事では，Metadynamicsの派生手法の一つであるParallel Bias Metadynamics（参考文献１）について解説します。

通常のMetadynamicsの問題点

分子シミュレーションの計算を実行する際に，系内の複数の自由度に興味があるのはよくあるシチュエーションかと思われます。
そのような場合に興味ある全自由度を通常のMetadynamicsで実行すると現実的には大きな困難にさいなまれます。
理由は，対象となる自由度を増やすことにより，それらの自由度で構成される探索空間の広さが指数関数的に増えてしまうからです（例えば，10自由度の場合，10次元空間を網羅的に探索する必要があります）。個人的には，３自由度の段階で相当に苦しくなる印象があります。
計算コストの困難を避けるため，指数関数的に計算コストを増大させないよう複数自由度を取扱う手法がいろいろと提案されてきました。Parallal Bias Metadynamicsもそのような背景の下，提案された手法になります。今回の主題とは外れますが，ほかにもBias Exchange Metadynamics（参考文献２）やSGOOP法（参考文献３）といった手法もこれまで提案されてきました。

数式に出てくる変数・パラメータの定義

Parallel Bias Metadynamicsの説明に移る前に，本記事で出てくる数式内の変数やパラメータについて，先にまとめて記しておきます。

$$
\begin{align*}
N&: 系内の粒子数\\
\mathbf{r}_i&:i番目の粒子の座標(=\begin{bmatrix}x_i&y_i&z_i\end{bmatrix}^{\rm T})\\
\mathbf{R}&:N粒子の座標を１表記にまとめたもの(=\begin{bmatrix}\mathbf{r}_1^{\rm T}&\cdots&\mathbf{r}_N^{\rm T}\end{bmatrix}^{\rm T})\\
\hat{\mathbf{s}}(\mathbf{R})&:\mathbf{R}の関数として定義された{\rm metadynamics}で対象となる複数自由度をベクトルとしてまとめたもの\\
\mathbf{s}&:{\rm metadynamics}で対象となる複数自由度をベクトルとしてまとめたもの\\
D&:{\rm metadynamics}で対象となる複数自由度の数\\
w_0&:ガウシアンカーネルの高さに関する量\\
\sigma_{\alpha}&:s_{\alpha}に課せられるガウシアンカーネルの幅\\
k_{\rm B}&:ボルツマン定数\\
T&:温度\\
\beta&:逆温度\left(=\frac{1}{k_{\rm B}T}\right)\\
t&:時間
\end{align*}
$$

Parallel Bias Metadynamicsのバイアスポテンシャル

Parallel Bias Metadynamicsでは内部ポテンシャル$${U(\mathbf{R})}$$の他に以下のバイアスポテンシャルを課していきます。

$$
\begin{align*}
V(\hat{\mathbf{s}}(\mathbf{R}),t)&=-\frac{1}{\beta}\ln\left\{\sum_{\alpha}{\rm e}^{-\beta V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t)}\right\}
\end{align*}\tag{1}
$$

ここで，$${V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t)}$$は

$$
\begin{align*}
V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t)&=\int_0^{t}{\rm d}\tau w_{\alpha}(\tau)\exp\left(-\frac{(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}(\tau))^2}{2\sigma_{\alpha}^2}\right)\\
w_{\alpha}(\tau)&=w_0\times\left[\frac{\exp\left(-\beta V_{\alpha}(s_{\alpha}(\tau))\right)}{\sum_{\alpha'}\exp\left(-\beta V_{\alpha'}(s_{\alpha'}(\tau))\right)}\right]
\end{align*}\tag{2}
$$

となり，１自由度に対するmetadynamicsのガウシアンカーネルに対応するようなバイアスポテンシャルが採用されていることになります。
ただし，通常のmetadynamicsと異なり，ガウシアンカーネルの高さは定数ではなく$${w_{\alpha}(\tau)}$$に従って時間によって異なる値を取ります。
各自由度に対するガウシアンカーネルの高さ$${w_{\alpha}(\tau)}$$の和が

$$
\begin{align*}
\sum_{\alpha}w_{\alpha}(\tau)&=w_0
\end{align*}\tag{3}
$$

となるよう制限されていることになります。
更にwell-tempered metadynamicsと組み合わせることも可能ですが，今回は考慮しないで議論を進めます。

バイアスポテンシャル有無の確率分布

バイアスポテンシャルが存在する場合の$${\mathbf{R}}$$の確率分布を$${P_{b}(\mathbf{R})}$$，バイアスポテンシャルが存在しない場合の$${\mathbf{R}}$$の確率分布を$${P_{0}(\mathbf{R})}$$とします。
metadynamicsがquasi-stationary limit（=バイアスポテンシャルの変化より系の平衡化が十分に早く，確率分布がボルツマン分布に従うと仮定できる状態。metadynamicsが収束したと表現することもあります）に達した場合，$${P_{b}(\mathbf{R},t)}$$は

$$
\begin{align*}
P_{b}(\mathbf{R},t)&=\frac{\exp\left(-\beta(U(\mathbf{R})+V(\hat{\mathbf{s}}(\mathbf{R}),t))\right)}{\int{\rm d}\mathbf{R}\exp\left(-\beta(U(\mathbf{R})+V(\hat{\mathbf{s}}(\mathbf{R}),t))\right)}
\end{align*}\tag{4}
$$

式(4)の右辺の分母に式(2)を代入して式変形すると，

$$
\begin{align*}
&\int{\rm d}\mathbf{R}{\rm e}^{-\beta(U(\mathbf{R})+V(\hat{\mathbf{s}}(\mathbf{R}),t))}\\
&=\sum_{\alpha}\int{\rm d}\mathbf{R}{\rm e}^{-\beta(U(\mathbf{R})+V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t))}\\
&=\sum_{\alpha}\int{\rm d}s_{\alpha}\int{\rm d}\mathbf{R}\delta(s_{\alpha}-\hat{s}_{\alpha}(\mathbf{R})){\rm e}^{-\beta(U(\mathbf{R})+V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t))}\\
&=\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta V_{\alpha}(s_{\alpha},t)}\int{\rm d}\mathbf{R}\delta(s_{\alpha}-\hat{s}_{\alpha}(\mathbf{R})){\rm e}^{-\beta U(\mathbf{R})}\\
&=\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta V_{\alpha}(s_{\alpha},t)}\int{\rm d}\mathbf{R}{\rm e}^{-\beta U(\mathbf{R})}\frac{\int{\rm d}\mathbf{R}\delta(s_{\alpha}-\hat{s}_{\alpha}(\mathbf{R})){\rm e}^{-\beta U(\mathbf{R})}}{\int{\rm d}\mathbf{R}{\rm e}^{-\beta U(\mathbf{R})}}\\
&=\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta V_{\alpha}(s_{\alpha},t)}\int{\rm d}\mathbf{R}{\rm e}^{-\beta U(\mathbf{R})}P_0(s_{\alpha})\\
\end{align*}\tag{5}
$$

と変形できます。ここで，$${P_0(s_{\alpha})}$$はバイアスポテンシャルがない場合の自由度$${s_{\alpha}}$$の確率分布であり，バイアスポテンシャルがない場合の自由度$${s_{\alpha}}$$の平均力ポテンシャル$${F(s_{\alpha})}$$と

$$
\begin{align*}
P_0(s_{\alpha})&=\frac{{\rm e}^{-\beta F(s_{\alpha})}}{\int{\rm d}s_{\alpha}{\rm e}^{-\beta F(s_{\alpha})}}
\end{align*}\tag{6}
$$

の関係があります。
式(5)，式(6)を式(4)に代入すると，

$$
\begin{align*}
P_{b}(\mathbf{R},t)&=\frac{{\rm e}^{-\beta(U(\mathbf{R})+V(\hat{\mathbf{s}}(\mathbf{R}),t))}}{\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta V_{\alpha}(s_{\alpha},t)}\int{\rm d}\mathbf{R}{\rm e}^{-\beta U(\mathbf{R})}P_0(s_{\alpha})}\\
&=\frac{{\rm e}^{-\beta U(\mathbf{R})}}{\int{\rm d}\mathbf{R}{\rm e}^{-\beta U(\mathbf{R})}}\frac{{\rm e}^{-\beta V(\hat{\mathbf{s}}(\mathbf{R}),t)}}{\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta(F(s_{\alpha})+V_{\alpha}(s_{\alpha},t))}/\int{\rm d}s_{\alpha}{\rm e}^{-\beta F(s_{\alpha})}}\\
&=P_{0}(\mathbf{R})\frac{{\rm e}^{-\beta V(\hat{\mathbf{s}}(\mathbf{R}),t)}}{\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta(F(s_{\alpha})+V_{\alpha}(s_{\alpha},t))}/\int{\rm d}s_{\alpha}{\rm e}^{-\beta F(s_{\alpha})}}\\
\end{align*}\tag{7}
$$

となり，バイアスポテンシャルが存在する場合の$${\mathbf{R}}$$の確率分布をバイアスポテンシャルが存在しない場合の$${\mathbf{R}}$$の確率分布との関係式に帰着させることができます。

各自由度の平均力ポテンシャル

Parallel Bias Metadynamicsを収束するまで実行することにより，各自由度$${s_{\alpha}}$$に対する平均力ポテンシャルを得ることができます。
自由度$${s_{\alpha}}$$のバイアスポテンシャルがある場合の確率分布を$${P_b(s_{\alpha},t)}$$とすると，式(7)を用いて

$$
\begin{align*}
P_b(s_{\alpha},t)&=\int{\rm d}\mathbf{R}\delta\left(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}\right)P_{b}(\mathbf{R},t)\\
&=\frac{\int{\rm d}\mathbf{R}\delta\left(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}\right)P_{0}(\mathbf{R}){\rm e}^{-\beta V(\hat{\mathbf{s}}(\mathbf{R}),t)}}{\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta(F(s_{\alpha})+V_{\alpha}(s_{\alpha},t))}/\int{\rm d}s_{\alpha}{\rm e}^{-\beta F(s_{\alpha})}}\\
&=\frac{\int{\rm d}\mathbf{R}\delta\left(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}\right)P_{0}(\mathbf{R})\sum_{\alpha'}{\rm e}^{-\beta V(\hat{s}_{\alpha'}(\mathbf{R}),t)}}{\sum_{\alpha}\int{\rm d}s_{\alpha}{\rm e}^{-\beta(F(s_{\alpha})+V_{\alpha}(s_{\alpha},t))}/\int{\rm d}s_{\alpha}{\rm e}^{-\beta F(s_{\alpha})}}\\
&\propto \int{\rm d}\mathbf{R}\delta\left(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}\right)P_{0}(\mathbf{R})\sum_{\alpha'}{\rm e}^{-\beta V_{\alpha}(\hat{s}_{\alpha'}(\mathbf{R}),t)}
\end{align*}\tag{8}
$$

と表すことができます。
ここで，$${w_{\alpha}(\tau)}$$が式(2)に従うことを思い出すと，十分に収束した場合には，

$$
\begin{align*}
V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t)&\simeq V_{\alpha'}(\hat{s}_{\alpha'}(\mathbf{R}),t)
\end{align*}\tag{9}
$$

が成立することが期待されます。
式(9)を考慮して式(8)を更に式変形すると，

$$
\begin{align*}
P_b(s_{\alpha},t)&\propto D\int{\rm d}\mathbf{R}\delta\left(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}\right)P_{0}(\mathbf{R}){\rm e}^{-\beta V_{\alpha}(\hat{s}_{\alpha}(\mathbf{R}),t)}\\
&\propto D{\rm e}^{-\beta V_{\alpha}(s_{\alpha},t)}\int{\rm d}\mathbf{R}\delta\left(\hat{s}_{\alpha}(\mathbf{R})-s_{\alpha}\right)P_{0}(\mathbf{R})\\
&\propto \frac{D{\rm e}^{-\beta (F(s_{\alpha})+V_{\alpha}(s_{\alpha},t))}}{\int{\rm d}s_{\alpha}{\rm e}^{-\beta F(s_{\alpha})}}\\
&\propto{\rm e}^{-\beta (F(s_{\alpha})+V_{\alpha}(s_{\alpha},t))}
\end{align*}\tag{10}
$$

となります。
自由度がdiffusive（$${P_b(s_{\alpha},t)\simeq s_{\alpha}に依存しない定数}$$）になるまでmetadynamicsを実行したとすると，

$$
\begin{align*}
F(s_{\alpha})+V_{\alpha}(s_{\alpha},t)&=c(t)\\
\therefore F(s_{\alpha})&=-V_{\alpha}(s_{\alpha},t)+c(t)
\end{align*}\tag{11}
$$

となり，$${V_{\alpha}(s_{\alpha},t)}$$から$${s_{\alpha}}$$の平均力ポテンシャルが得られるという結論に帰着します。

参考文献

  1. Jim Pfaendtner and Massimiliano Bonomi, J. Chem. Theory Comput. 2015, 11, 11, 5062–5067

  2. S. Piana and A. Laio, J. Phys. Chem. B 111 (17) (2007) 4553–9

  3. Pratyush Tiwarya and B. J. Berne, Proceedings of the National Academy of Sciences 113 (11), 2839-2844, 2016-02-29


この記事が気に入ったらサポートをしてみませんか？