ハイブリッドモンテカルロ法とタンパク質フォールディング問題

概要

本稿では運動方程式を用いたモンテカルロ法である一般化ハイブリッドモンテカルロ法をとりあげ，方法論の基礎について比較的詳しく述べた．また，マルチカノニカル分布を生成する一般化ハイブリッドモンテカルロ法を粗視化したモデルタンパク質に適用し，その熱力学的性質を調べた例について解説した．

1　 はじめに

分子シミュレーションの主要な方法がモンテカルロ法と分子動力学法であることは言を俟たない．モンテカルロ法は確率過程の定常状態としてカノニカル分布のような目的とする分布を実現する．分子動力学法は対象とする系の力学を数値的に解くことにより，その性質を調べる．統計力学的な観点からはミクロカノニカル分布が生成されている．どちらの方法も長所・短所があり，どちらが優れた方法なのか一口に判定することは難しい．研究開発のツールとして使う場合にはシミュレーションの汎用アプリで利用できる分子動力学法を選ぶことが多いかもしれない．モンテカルロ法は効率を左右するパラメータの最適化や場合によっては対象とする系に応じたプログラミングも必要なので少しハードルが高いと感じるかもしれない．

さて，本稿はハイブリッドモンテカルロ（HMC）法［1–3］について解説する．この方法はモンテカルロ法の一種であり，運動方程式を用いて試行配位を生成することにより，効率的に大域的な状態更新を可能とするところに特色がある．HMC法は，元々は経路積分モンテカルロ計算を分子動力学法を使って実現するというアイデアに端を発する［4–6］．化学物理のコミュニティでは，経路積分量子化された系を古典的な“高分子”の系と解釈し［7,8］，方法論が構築されている．その高分子内のモード間でエネルギー移動がスムーズに進まず，単純に分子動力学法を適用するだけではエルゴード性に関する問題が生じる場合があることが知られている［6,9］．当初は分子動力学計算の中で定期的に運動量のリサンプリングを行うなど［10］アドホックな対処がなされていた．この方向で定式化されたのがHMC法であり，運動方程式の離散化誤差を伴わないカノニカル分布を厳密に生成する手法である［1］．分子動力学法ベースで能勢-Hoover鎖のような非線形性の強いサーモスタットを用いて経路積分計算に伴うエルゴード問題を回避する手法も確立している［3,11］．本稿ではHMC法の基礎から比較的詳しく解説したので，興味を持った読者がモンテカルロ法を選択するハードルが少しでも下がれば幸いである．

2　 方法論

この節ではハイブリッドモンテカルロ（HMC）法の方法論的な基礎について詳述する．標準的なHMC法を少し一般化した方法（一般化ハイブリッドモンテカルロ法）をその対象とする［12–15］．標準的なHMC法は運動量をマクスウェル分布から選び，その運動量と現在の系の配位から出発して数ステップMD計算を行い，その結果のハミルトニアン誤差から試行配位の採択・不採択を決定する方法である．一般化ハイブリッドモンテカルロ（GHMC）法は最初に運動量を用意する際に，前のステップでの運動量を部分的に再利用できるようよう一般化されている．この節では主にカノニカル分布の生成について考え，2.4節でマルチカノニカル分布の生成について考察する．

ここで共通する表記をまとめておこう．対象とする系の座標と運動量を

  

$$\mathit{q}=(q_1,\dots ,q_f),\mathit{p}=(p_1,\dots ,p_f)$$

と集合的に表す．ここで$f$は自由度の数である．系は逆温度$\beta$のカノニカル分布に従っている．分布関数はハミルトニアン$H(\mathit{q},\mathit{p})$を用いて

  

$$\begin{array}{c}\hfill P_{\mathrm{c}}(\mathit{q},\mathit{p})\propto e^{-\beta H(\mathit{q},\mathit{p})}\end{array}$$(1)

と書くことができる．ハミルトニアンは慣性を表すパラメータ$m$を用いて，

  

$$\begin{array}{c}\hfill H(\mathit{q},\mathit{p})=\frac{{\mathit{p}}^2}{2m}+U\left(\mathit{q}\right)\end{array}$$(2)

で与えられる．$U\left(\mathit{q}\right)$は相互作用ポテンシャルである．以下，運動量を用意する部分と運動方程式により試行配位を生成する部分を2.1と2.2節で個別に説明する．

2.1　 ランダム運動量の混合

この小節では表記の煩雑さを避けるため，1次元で考えるが，多次元への一般化は容易である．逆温度$\beta$での運動量分布は

  

$$\begin{array}{c}\hfill P_{\mathrm{G}}(p)=\sqrt{\frac{\beta }{2\pi m}}{e}^{-\frac{\beta p^2}{2m}}\end{array}$$(3)

と表されるマクスウェル分布（平均0，分散$m/\beta$のガウス分布）である．次のように現在の運動量に同じ分布に従うランダム運動量$\xi$を混合する試行を考える．

  

$$\begin{array}{c}\hfill \left(\genfrac{}{}{0ex}{}{p^{\prime }}{{\xi }^{\prime }}\right)=\left(\begin{array}{cc}\hfill \cos \varphi \hfill & \hfill \sin \varphi \hfill \\ \hfill -\sin \varphi \hfill & \hfill \cos \varphi \hfill \end{array}\right)\left(\genfrac{}{}{0ex}{}{p}{\xi }\right)\end{array}$$(4)

ここで$0<\varphi \le \pi /2$である．混合後の運動量$p^{\prime }$は，

  

$$\begin{array}{c}\hfill p^{\prime }=p\cos \varphi +\xi \sin \varphi \end{array}$$(5)

で与えられる．運動量$p^{\prime }$の分布は，

  

$$\begin{array}{ccc}P\left(p^{\prime }\right)\hfill & =\hfill & {\int }_{-\infty }^{\infty }d\xi {\int }_{-\infty }^{\infty }dp\hfill \\ & & \delta (p^{\prime }-(p\cos \varphi +\xi \sin \varphi ))P_{\mathrm{G}}(p)P_{\mathrm{G}}(\xi )\\ & =\hfill & \sqrt{\frac{\beta }{2\pi m}}{e}^{-\frac{\beta {p^{\prime }}^2}{2m}}=P_{\mathrm{G}}\left(p^{\prime }\right)\hfill \end{array}$$(6)

であるので，同じマクスウェル分布に従っていることがわかる．そのため，この試行は確率1で採択される．

2.2　 運動方程式による試行配位の生成

まずは位相空間内での状態更新について一般的に考察する．状態の更新はマルコフ過程を用いる．状態更新により，位相空間の“配位”$(\mathit{q},\mathit{p})$は$({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })$に移るとする．この遷移確率を

  

$$\begin{array}{c}\hfill P_{\mathrm{M}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\end{array}$$(7)

と書く．この遷移により目的とする分布が乱されないためには，詳細釣り合い

  

$$\begin{array}{c}{P}_{\mathrm{c}}(\mathit{q},\mathit{p})P_{\mathrm{M}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\hfill \\ =P_{\mathrm{c}}({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })P_{\mathrm{M}}(({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })\to (\mathit{q},\mathit{p}))\hfill \end{array}$$(8)

が満たされていれば十分である．ここで遷移確率を配位の提案確率$P_{\mathrm{S}}$とその試行配位の採択確率$P_{\mathrm{A}}$の積で表す．

  

$$\begin{array}{c}{P}_{\mathrm{M}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\hfill \\ =P_{\mathrm{S}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))P_{\mathrm{A}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\hfill \end{array}$$(9)

詳細釣り合い(8)を満たすメトロポリス型の採択確率は以下のようになる．

  

$$\begin{array}{c}{P}_{\mathrm{A}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\hfill \\ =\min \{1,\frac{P_{\mathrm{c}}({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })P_{\mathrm{S}}(({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })\to (\mathit{q},\mathit{p}))}{P_{\mathrm{c}}(\mathit{q},\mathit{p})P_{\mathrm{S}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))}\}\hfill \end{array}$$(10)

よく見かけるタイプのものは提案確率が定数であり，その部分が少し一般化されている¹［16］．さて，GHMC法では試行配位は二段階を経て生成する．まず運動方程式により配位$(\mathit{q},\mathit{p})=(\mathit{q}(0),\mathit{p}(0))$を時間$\tau$だけ進め，

  

$$\begin{array}{c}\hfill (\overline{\mathit{q}},\overline{\mathit{p}})=(\mathit{q}(\tau ),\mathit{p}(\tau ))\end{array}$$(11)

その後，運動量を反転する．

  

$$\begin{array}{c}\hfill ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })=(\overline{\mathit{q}},\mathbf{-}\overline{\mathit{p}})\end{array}$$(12)

運動方程式による試行配位の生成は決定論的に行われるため，配位$(\overline{\mathit{q}},\overline{\mathit{p}})$を生成する確率は

  

$$\begin{array}{c}\hfill P_{\mathrm{H}}((\mathit{q},\mathit{p})\to (\overline{\mathit{q}},\overline{\mathit{p}}))=\delta ((\overline{\mathit{q}},\overline{\mathit{p}})-(\mathit{q}\left(\mathit{\tau }\right),\mathit{p}\left(\mathit{\tau }\right)))\end{array}$$(13)

と与えられる．また，運動量反転については

  

$$\begin{array}{c}\hfill P_{\mathrm{F}}((\overline{\mathit{q}},\overline{\mathit{p}})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))=\delta (({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })-(\overline{\mathit{q}},\mathbf{-}\overline{\mathit{p}}))\end{array}$$(14)

である．提案確率は，この二つを用いて，

  

$$\begin{array}{c}{P}_{\mathrm{S}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\hfill \\ ={\iint }_{}^{}d\overline{\mathit{q}}d\overline{\mathit{p}}{P}_{\mathrm{H}}((\mathit{q},\mathit{p})\to (\overline{\mathit{q}},\overline{\mathit{p}}))P_{\mathrm{F}}((\overline{\mathit{q}},\overline{\mathit{p}})\to ({\mathit{q}}^{\mathbf{\prime }},{\mathit{p}}^{\mathbf{\prime }}))\hfill \end{array}$$(15)

と表される．変数$(\overline{\mathit{q}},\overline{\mathit{p}})$に関する積分を実行すると

  

$$\begin{array}{c}\hfill P_{\mathrm{S}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))=\delta (({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })-(\mathit{q}(\tau ),-\mathit{p}(\tau )))\end{array}$$(16)

が得られる．同様に逆方向の提案確率は

  

$$\begin{array}{c}\hfill P_{\mathrm{S}}(({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })\to (\mathit{q},\mathit{p}))=\delta ((\mathit{q},\mathit{p})-(\mathit{q}(0),\mathit{p}(0)))\end{array}$$(17)

となり，力学の時間反転対称性より両者は一致する．

  

$$\begin{array}{c}\hfill P_{\mathrm{S}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))=P_{\mathrm{S}}(({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })\to (\mathit{q},\mathit{p}))\end{array}$$(18)

この関係により，採択確率(10)は以下のようになる．

  

$$\begin{array}{c}\hfill P_{\mathrm{A}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))=\min \{1,e^{-\beta \mathrm{\Delta }H}\}\end{array}$$(19)

ここで

  

$$\begin{array}{c}\hfill \mathrm{\Delta }H=H({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })-H(\mathit{q},\mathit{p})\end{array}$$(20)

である．詳細釣り合いについては，次のように確認することができる．関係式

  

$$\begin{array}{ccc}{e}^{-\beta H(\mathit{q},\mathit{p})}\min \{1,e^{-\beta \mathrm{\Delta }H}\}\hfill & =& \min \{e^{-\beta H(\mathit{q},\mathit{p})},e^{-\beta H({\mathit{q}}^{\mathbf{\prime }},{\mathit{p}}^{\mathbf{\prime }})}\}\hfill \\ & =& e^{-\beta H({\mathit{q}}^{\mathbf{\prime }},{\mathit{p}}^{\mathbf{\prime }})}\min \{e^{\beta \mathrm{\Delta }H},1\}\hfill \end{array}$$(21)

を用いると，

  

$$\begin{array}{c}{P}_{\mathrm{c}}(\mathit{q},\mathit{p})P_{\mathrm{A}}((\mathit{q},\mathit{p})\to ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime }))\hfill \\ =P_{\mathrm{c}}({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })P_{\mathrm{A}}(({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })\to (\mathit{q},\mathit{p}))\hfill \end{array}$$(22)

を直ちに示すことができる．式(9), (18)と(22)より，詳細釣り合い（式(8)）が満たされていることがわかる．採択された試行配位は，$P_{\mathrm{H}}$で時間発展した後に，運動量の符号が反転している．実用上は高い採択率が想定されているので，試行が棄却された場合に反転するほうが好ましい．そのためにもう一度最後に運動量反転の操作を行う．ハミルトニアンはこの操作に対して不変なので，反転操作は確率1で採択される．

最後に試行配位を生成する運動方程式についてコメントしておこう．ハミルトニアン式(2)に対する正準方程式は

  

$$\begin{array}{c}\frac{d\mathit{q}}{dt}=\frac{\partial H}{\partial \mathit{p}}=\frac{\mathit{p}}{m}\hfill \\ \frac{d\mathit{p}}{dt}=-\frac{\partial H}{\partial \mathit{q}}=-\frac{\partial U}{\partial \mathit{q}}\hfill \end{array}$$(23)

である．試行配位を生成する段階では目的とする分布のポテンシャル$U$でなく別のものを用いてもよく，計算を効率化できる可能性があることに注意しておく［1］．この運動方程式により生成された配位$(\overline{\mathit{q}},\overline{\mathit{p}})$は時間反転に対して対称であり，時間発展を変数変換とみなしたヤコビ行列式は

  

$$\begin{array}{c}\hfill \frac{\partial (\overline{\mathit{q}},\overline{\mathit{p}})}{\partial (\mathit{q},\mathit{p})}=1\end{array}$$(24)

である．拡張系の方法で使われる運動方程式は，式(24)のヤコビ行列式が1でないこともあるので［17］，その際は採択確率式(19)に修正が必要である［18］．また，運動方程式は有限の時間刻み$\mathrm{\Delta }t$の差分式を用いて数値的に解くことになる．その際には時間反転対称性や式(24)の性質を持っている手法が望ましい．シンプレクティック数値解法はその良い例である［19,20］．

2.3　 アルゴリズムの要約

ここでGHMC法のアルゴリズムを簡単にまとめる．初期配位は$(\mathit{q},\mathit{p})$である．まず逆温度$\beta$のマクスウェル乱数と現在の運動量を式(5)により混合する．その後，分子動力学計算を$n_{\mathrm{MD}}$ステップ行い，時間を$n_{\mathrm{MD}}\times \mathrm{\Delta }t$だけ進める．これに伴うハミルトニアン誤差$\mathrm{\Delta }H$を用いて確率$\min (1,e^{-\beta \mathrm{\Delta }H})$で試行配位を採択する．棄却された場合は初期配位の運動量を反転させる．

  

$$\begin{array}{c}\hfill ({\mathit{q}}^{\prime },{\mathit{p}}^{\prime })=(\mathit{q},\mathbf{-}\mathit{p})\end{array}$$(25)

運動量の反転がこの段階で行われるのは，2.2節での2回の反転操作を繰り込んだ結果である．

効率化を図るために調節することができるGHMCパラメータは，$\varphi ,n_{\mathrm{MD}},\mathrm{\Delta }t$の三つである．ランダム運動量の混合割合を定める角$\varphi$はGHMCの計算のあいだ数値を固定しても良いし，ランダム変数としても良い．$\varphi =\pi /2$の場合は，最初に運動量はマクスウェル乱数から選ばれるので標準的なHMC法に帰着する．棄却された時の式(25)の運動量反転操作は，次のステップでマクスウェル乱数により全て上書きされるので標準的なHMC計算の場合は必要ない．

2.4　 マルチカノニカル集合

ここまではカノニカル分布の生成方法について概観してきた．この小節では拡張アンサンブル法［21］のひとつであるマルチカノニカル分布［22］を取り上げ，そのGHMC法による生成について考えよう［23］．マルチカノニカル集合でのエネルギー分布関数は

  

$$\begin{array}{c}\hfill P_{\mathrm{mc}}(U)\propto \Omega (U)e^{-W(U)}\end{array}$$(26)

で与えられる．$\Omega (U)$はポテンシャルエネルギーに関する状態密度である．$W(U)$はマルチカノニカル重み関数であり，確率密度$P_{\mathrm{mc}}$がポテンシャルエネルギーに関して一定値になるよう導入される．また，配位空間ではマルチカノニカル分布関数は$e^{-W\left(U\left(\mathit{q}\right)\right)}$に比例する．ここで以下のような“マルチカノニカルポテンシャル”

  

$$\begin{array}{c}\hfill U_{\mathrm{mc}}\left(U\left(\mathit{q}\right)\right)=W\left(U\left(\mathit{q}\right)\right)/{\beta }_0\end{array}$$(27)

を定義する．これにより配位空間でのマルチカノニカル分布は逆温度${\beta }_0$でのカノニカル分布と解釈することができる．この読み替えにより，カノニカル分布を生成するGHMC法を用いてマルチカノニカル分布を生成することができる．運動を生成するハミルトニアンは

  

$$\begin{array}{c}\hfill H(\mathit{q},\mathit{p})=\frac{{\mathit{p}}^2}{2m}+U_{\mathrm{mc}}\left(\mathit{q}\right)\end{array}$$(28)

である．対応する運動方程式は［24,25］

  

$$\begin{array}{c}\frac{d\mathit{q}}{dt}=\frac{\partial H}{\partial \mathit{p}}=\frac{\mathit{p}}{m}\hfill \\ \frac{d\mathit{p}}{dt}=-\frac{\partial H}{\partial \mathit{q}}=-\frac{\partial U_{\mathrm{mc}}\left(\mathit{q}\right)}{\partial \mathit{q}}=-\frac{\partial U_{\mathrm{mc}}}{\partial U}\frac{\partial U}{\partial \mathit{q}}\hfill \end{array}$$(29)

となる．パラメータ${\beta }_0$は任意であるが，平坦化したポテンシャルエネルギー分布の最も高いエネルギー領域に対応する温度にとることが多いようである［23］．

3　 タンパク質の熱力学的フォールディング転移

この節ではごくシンプルなタンパク質のモデルにマルチカノニカル分布を生成するGHMC法を適用した例について概観する［23,26］．このモデルではアミノ酸残基を一つのビーズとして粗視化しており，モデルの提案者名の頭文字をとってHTモデルと呼ぶことにする［27］．また，HTモデルとエネルギー最安定構造は同一であるが，その構造を他の極小構造と比べて安定化させたモデルをGo-like HTモデル［28］と呼ぶことする．この名称は郷の整合性原理［29］にインスパイアされたモデルであることに由来する．本稿ではHTモデルは一般的なヘテロ高分子のモデルとして，Go-like HTモデルはタンパク質のモデルとして考える．

本稿では46残基からなるモデルを考える．具体的なビーズ間の相互作用については原論文等を参照して欲しい．このHTモデルのエネルギー最小構造は4本のβストランドからなるβバレル構造である．このモデルはエネルギー最小構造と同程度のエネルギーをもつ複数のエネルギー極小構造があることが知られている．次にエネルギー最小構造が他の極小構造より安定化されているモデルを構築しよう．これはエネルギー最小構造でコンタクトしていない残基間の引力相互作用を取り除くことにより実現する．このモデルがGo-like HTモデルである．また，本稿ではエネルギー最小構造を“天然構造”と定義することにし，この構造が支配的な状態を“天然状態”と呼ぶことにする．

最初にHTモデルにマルチカノニカルGHMC法を適用した場合，その計算効率のGHMCパラメータ依存性について示そう．図1は標準的なHMC（$\varphi =\pi /2$）に対して，相関時間を時間刻み$\mathrm{\Delta }t$の関数として示している．相関時間は相関しているモンテカルロステップ数を計算するために必要なCPU時間である．ここでは$n_{\mathrm{MD}}=10$と固定されている．相関時間は$\mathrm{\Delta }t=0.025$のところで極小になっており，対応する採択率は63%である．この時間刻みに固定して$n_{\mathrm{MD}}$を変え相関時間を調べた結果が図2である．$n_{\mathrm{MD}}=40-150$において相関時間にプラトー領域が見て取れる．高密度レナードジョーンズ（LJ）液体にマルチカノニカルGHMC法を適用した場合$n_{\mathrm{MD}}=10$が最適であったが［23］，タンパク質のモデルの場合，最適な$n_{\mathrm{MD}}$はLJ液体の場合と比べて大きい．このような傾向は孤立分子系や中低密度流体のカノニカルHMC計算でも見られる［30–32］．大きな$n_{\mathrm{MD}}$値を採用する場合，その中間での構造を利用する方法が提案されている［30］．しかしながら，複数のタイプの試行配位生成を組み合わせて利用する場合は，ロードバランスの観点から$n_{\mathrm{MD}}$は大きすぎないほうが望ましい．図3にランダム運動量の混合割合を変えた場合（$\varphi =\pi /8$）の相関時間を示す．この場合は$n_{\mathrm{MD}}=20$が極小となり，運動量を再利用することにより小さな$n_{\mathrm{MD}}$を用いても十分高い効率を得ることができることがわかる．

図1  時間刻み$\mathrm{\Delta }t$の関数としての相関時間$\tau$（丸印）および採択率（四角印）．GHMCパラメータ$\varphi ＝\pi /2,n_{\mathrm{MD}}=10$は固定．文献［23］より転載．

図2  $n_{\mathrm{MD}}$の関数としての相関時間$\tau$（丸印）および採択率（四角印）．GHMCパラメータ$\varphi ＝\pi /2,\mathrm{\Delta }t=0.025$は固定．文献［23］より転載．

図3  $n_{\mathrm{MD}}$の関数としての相関時間$\tau$（丸印）および採択率（四角印）．GHMCパラメータ$\varphi ＝\pi /8,\mathrm{\Delta }t=0.025$は固定．文献［23］より転載．

次に粗視化モデルの熱力学的な性質を見てみよう．図4にマルチカノニカル計算から得られたポテンシャルエネルギー分布を示す．HTモデル，Go-like HTモデルともに広いエネルギー領域で一定値であり，効率的なサンプリングが実現されていることがわかる．そのために，得られた分布からからカノニカル分布での期待値をReweighting法により精度良く求めることができる．任意の物理量$A$の逆温度$\beta$におけるカノニカル平均${\mathcal{⟨}A\mathcal{⟩}}_{\mathrm{c}}$は次のように求められる．

図4  マルチカノニカルGHMC計算により生成されたポテンシャルエネルギー分布$H(U)$．上のパネルはHTモデル，下のモデルはGo-like HTモデルの分布である．文献［26］より転載．

  

$$\begin{array}{c}\hfill {\mathcal{⟨}A\mathcal{⟩}}_{\mathrm{c}}=\frac{{\mathcal{⟨}A{e}^{W(U)-\beta U}\mathcal{⟩}}_{\mathrm{mc}}}{{\mathcal{⟨}{e}^{W(U)-\beta U}\mathcal{⟩}}_{\mathrm{mc}}}\end{array}$$(30)

ここで${\mathcal{⟨}\cdots \mathcal{⟩}}_{\mathrm{mc}}$はマルチカノニカル分布による期待値を示す．カノニカル分布による内部エネルギー$E$，

  

$$\begin{array}{c}\hfill E=\frac{3N}{2\beta }+{\mathcal{⟨}U\mathcal{⟩}}_{\mathit{c}}\end{array}$$(31)

および，熱容量$C$は

  

$$\begin{array}{c}\hfill \frac{C}{k_{\mathrm{B}}}=\frac{3N}{2}+{\beta }^2\{{\mathcal{⟨}{U}^2\mathcal{⟩}}_c-{\mathcal{⟨}U\mathcal{⟩}}_c^2\}\end{array}$$(32)

と与えられる．ここで$k_B$はボルツマン定数である．一方，ミクロカノニカル分布に基づき熱力学量を計算することもできる．内部エネルギーに関する状態密度$D(E)$は次のように求められる．

  

$$\begin{array}{c}\hfill D(E)=C_N{\int }_{U_0}^E(E-U{)}^{\frac{3N}{2}-1}\mathrm{\Omega }(U)dU\end{array}$$(33)

$C_N$はエネルギーに依らない定数である．エントロピーは次のとおり状態密度により与えられる．

  

$$\begin{array}{c}\hfill \mathrm{S}(E)=k_{\mathrm{B}}\ln D(E)\end{array}$$(34)

このエントロピーから温度$T$は

  

$$\begin{array}{c}\hfill \frac{1}{T(E)}=\frac{\partial \mathrm{S}(E)}{\partial E}\end{array}$$(35)

と求めることができる．図5にカノニカル分布から計算した内部エネルギー$E$，熱容量$C$，慣性半径$R_g$を温度の関数として示す．HTモデル，Go-like HTモデルともに内部エネルギーは温度の低下と共に減少する．Go-like HTモデルは，温度$T=0.6$あたりで急激な内部エネルギーの減少が見られる．その大きな変化率は熱容量のピークとなり，相転移様の振る舞いを示している．慣性半径をみると，熱容量にピークを与える温度近傍で，急激に減少していることが見てとれる．つまりこの温度で空間的に広がっている状態からコンパクトな構造が支配的な状態に転移している．HTモデルもGo-like HTモデルほどシャープではないが，同様の振る舞いを示している．この転移をコラプス転移と呼ぶことにしよう．後述のようにHTモデルでは，この転移は天然状態にうつるものではない．一方，Go-like HTモデルでは，天然構造の4本のβストランドのうち，3本が天然構造と同一であり，部分的に天然構造にフォールドしている状態である．また図にはミクロカノニカル分布により求めた温度も示してある．HTモデルは統計集合による違いは見られないが，Go-like HTモデルはコラプス転移温度付近に顕著な統計集合依存性が見られる．ミクロカノニカル分布による結果に見られるファンデルワールスループ様の振る舞い（内部エネルギーのS字型温度依存性）は，Go-like HTモデルのコラプス転移，つまり部分的なフォールド状態への転移は，一次相転移的な性質を持っていることを示している．上述の温度－エネルギー曲線の統計集合依存性はGo-like HTモデルのコラプス転移が一次相転移的であり，かつ系が有限であることに由来すると考えられる．

図5  カノニカル分布により計算したHTモデルおよびGo-like HTモデルの内部エネルギー$E$，熱容量$C$，慣性半径$R_g$．横軸は温度$T$．青の実線はHTモデル，赤の実線はGo-like HTモデルの結果である．内部エネルギーの破線はミクロカノニカル分布により計算したもの．文献［26］より転載．

天然構造との“距離”を適切に導入することにより，秩序変数$Q$を定義する．図6に秩序変数の温度依存性を示す．Go-like HTモデルでは$T=0.6$あたりで急激に増大し，さらなる温度の低下と共に1へ近づいていく．またそのゆらぎは$T=0.57$でピークを持っており，この温度で天然状態へ転移していることがわかる．これをフォールディング転移と呼ぶことにしよう．図5の熱容量は，この転移温度あたりに裾野をひいており，転移の影響が見てとれる．Go-like HTモデルは，温度の低下と共に一次相転移様のコラプス転移により部分的に天然構造に折りたたまれた状態になり，引き続きわずかに温度を低下させると最終的にフォールディング転移を経て天然状態となる．一方，HTモデルは温度の低下とともに秩序変数は徐々に増大し，その分散は$T=0.19$にピークがあり，コラプス転移とフォールデング転移の温度差は大きい．この温度差が小さい場合は，速く天然構造に折りたたまることが経験的に知られており［33］，Go-like HTモデルをタンパク質のモデルとしていることと整合している．詳しくは文献［23,26,34］を参照して欲しい．

図6  カノニカル分布により計算した秩序変数の期待値（上のパネル）およびその分散（下のパネル）．青の実線はHTモデル，赤の実線はGo-like HTモデルの結果である．文献［26］より転載．

4　 まとめ

本稿では運動方程式にガイドされたモンテカルロ法であるハイブリッドモンテカルロ（HMC）法について解説した．筆者のグループではHMC法をベースにした方法の開発をいくつか行っている．手前味噌で恐縮だが，興味ある読者の役に立つことを祈りつつ，この場を借りて紹介したい．量子モンテカルロ法の枠組みでは，ボーズ粒子系の経路積分ハイブリッドモンテカルロ（PIHMC）法［31,35］，基底状態の厳密計算手法である変分経路積分ハイブリッドモンテカルロ（VPIHMC）法［32］が挙げられる．PIHMC法では量子化された剛体も取り扱えるよう方法を拡張し，分子をドープした超流動ヘリウムクラスターの性質を調べた［36,37］．VPIHMC法を液体ヘリウム に適用した例は文献［38］である．VPIHMC法とレプリカ交換法と結合して，パラ水素クラスターの基底状態を探索した例は文献［39］を参照して欲しい．古典系の手法としては，マルチカノニカル分布を生成する一般化ハイブリッドモンテカルロ（GHMC）法を手掛けた［23］．この手法を粗視化タンパク質モデルに適用した例は本稿の3節で紹介した．全原子モデルによるシニョリンのフォールディング転移の研究にも利用されている［40］．また，HMC法を用いたマルチカノニカル分布生成の効率については文献［41］を，生体分子に対するGHMC法の計算効率については文献［42］を参照して欲しい．

注

1  この種の一般化は異なるコミュニティで独立に発見されている．例えば，文献［16］．

参考文献

• [1]   S.  Duane,  A. D.  Kennedy,  B. J.  Pendleton and  D.  Roweth, Phys. Lett. B, 195, 216 (1987).
• [2]   B.  Mehlig,  D. W.  Heermann and  B. M.  Forrest, Phys. Rev. B, 45, 679 (1992).
• [3]   M.  Tuckerman,  B. J.  Berne,  G. J.  Martyna and  M. L.  Klein, J. Chem. Phys., 99, 2796 (1993).
• [4]   D.  Callaway and  A.  Rahman, Phys. Rev. Lett., 49, 613 (1982).
• [5]   M.  Parrinello and  A.  Rahman, J. Chem. Phys., 80, 860 (1984).
• [6]   R.  Hall and  B.  Berne, J. Chem. Phys., 81, 3641 (1984).
• [7]   J. A.  Baker and  D.  Henderson, Rev. Mod. Phys., 48, 587 (1976).
• [8]   D.  Chander and  P.  Wolynes, J. Chem. Phys., 74, 4078 (1981).
• [9]   J. B.  Kogut, Eur. Phys. J. H, 49, 25 (2024).
• [10]   J.  Schnitker and  P.  Rossky, J. Chem. Phys., 86, 3471 (1987).
• [11]   三浦 伸一，アンサンブル，28, 6 (2004).
• [12]   A. M.  Horowitz, Phys. Lett. B, 268, 247 (1991).
• [13]   A. D.  Kennedy and  B.  Pendleton, Nucl. Phys. B, 607, 456 (2001).
• [14]   E.  Akhmatskaya,  N.  Bou-Rabee and  S.  Reich, J. Comput. Phys., 228, 2256 (2009).
• [15]   E.  Akhmatskaya,  N.  Bou-Rabee and  S.  Reich, J. Comput. Phys., 228, 7492 (2009).
• [16]  M. H. Kalos and P. A. Whitlock, Monte Carlo Methods, Volume I: Basics, Wiley, New York (1986).
• [17]   D.  Evans and  G.  Morriss, Statistical Mechanics of Nonequilibrium Liquids, Cambridge University Press, Cambridge (2008).
• [18]   W.  Shinoda,  M.  Shiga and  M.  Mikami, Phys. Rev. B, 69, 134103 (2004).
• [19]  大貫義郎，吉田春夫，現代の物理学1 力学，岩波書店 (1994).
• [20]   E.  Hairer,  C.  Lubich and  G.  Wanner, Geometric Numerical Integration 2nd ed., Springer, Heidelberg (2006).
• [21]  岡崎進，岡本祐幸編，生体系のコンピュータ・シミュレーション，化学同人 (2002).
• [22]   B. A.  Berg and  T.  Neuhaus, Phys. Rev. Lett., 68, 9 (1992).
• [23]   N.  Mukuta and  S.  Miura, J. Chem. Phys., 149, 072322 (2018).
• [24]   U.  Hansmann,  Y.  Okamoto and  F.  Eisenmenger, Chem. Phys. Lett., 59, 321 (1996).
• [25]   N.  Nakajima,  H.  Nakamura and  A.  Kidera, J. Phys. Chem. B, 101, 817 (1997).
• [26]   N.  Mukuta and  S.  Miura, Biophys. Physicobiol., 17, 14 (2020).
• [27]   J. F.  Honeycutt and  D.  Thirumalai, Proc. Natl. Acad. Sci. USA, 87, 3526 (1990).
• [28]   J.  Kim and  T.  Keyes, J. Phys. Chem. B, 112, 954 (2008).
• [29]   N.  Go, Annu. Rev. Biophys. Bioeng., 12, 183 (1983).
• [30]   N.  Matubayasi and  M.  Nakahara, J. Chem. Phys., 110, 3291 (1999).
• [31]   S.  Miura and  J.  Tanaka, J. Chem. Phys., 120, 2160 (2004).
• [32]   Y.  Kamibayashi and  S.  Miura, J. Chem. Phys., 145, 074114 (2016).
• [33]   D.  Klimov and  D.  Thirumalai, Phys. Rev. Lett., 76, 4070 (1996).
• [34]  三浦伸一，技術情報協会編，タンパク質の構造解析手法とIn silicoスクリーニングへの応用事例，第5章第5節，技術情報協会 (2023).
• [35]   S.  Miura and  S.  Okazaki, Chem. Phys. Lett., 308, 115 (1999).
• [36]   S.  Miura, J. Chem. Phys., 126, 114308 (2007).
• [37]   S.  Miura, J. Chem. Phys., 126, 114309 (2007).
• [38]  S. Miura, in: S. Tanaka, S. Rothstein, W. Lester (Eds.), Advances in Quantum Monte Carlo, the ACS Symposium Series 1094, Chapter 14, American Chemical Society, Washington DC (2012).
• [39]   S.  Miura, J. Chem. Phys., 148, 102333 (2018).
• [40]   Y.  Yamada,  D.  Suzuki and  S.  Miura, J. Phys. Soc. Jpn., 93, 084804 (2024).
• [41]   T.  Tagawa,  T.  Kaneko and  S.  Miura, Mol. Simul., 43, 13 (2017).
• [42]   D.  Suzuki,  T.  Hori and  S.  Miura, J. Phys.: Conf. Ser., 2207, 012018 (2022).

著者紹介

三浦 伸一（博士（理学））

〔経歴〕1995年京都大学理学研究科博士課程単位取得満期退学，その後，ペンシルバニア大学博士研究員，東京工業大学総合理工学研究科助手，分子科学研究所助手を経て，2007年から金沢大所属．〔専門〕液体物性理論，量子モンテカルロ．〔趣味〕読書，海外旅行，茶の湯．