分子動力学法

数学者のための統計力学３：物理に寄り添った Monte Carlo 法

Nature
Computation
Author

司馬博文

Published

6/29/2024

概要
本質的に Metropolis 法がサンプリング法であるならば，MD 法は \(N\)-体問題に対する数値解法であると言える．しかし，Hamiltonian Monte Carlo は元々 Monte Carlo 法と MD 法との融合を目指したものであること，Event-Chain Monte Carlo 法も MD 法における古典的手法の輸入と理解できること，Langevin 動力学も正準集団に対する MD 法と捉えられることを考えると，尽きぬ計算テクニックの源泉であると言える．

1 概観

1.1 はじめに

分子動力学法 (MD: Molecular Dynamics simulation) は，系のミクロなダイナミクスを実験的に調べる代わりに，シミュレーションによって必要条件を絞り込んでいく構成的な・計算機集約的なアプローチである．

特に，多体系の運動方程式を計算機を用いて数値的に解くことを指す (Alder and Wainwright, 1959)．それ故，物理や化学はもちろん，生化学，物質科学，さらには工学分野にまで幅広く用いられるシミュレーション法になる．

the inescapable conclusion is that MD will — if it hasn’t already — become an indispensable part of the theorist’s toolbox. (Rapaport, 2004, p. ix)

従来は実験と理論が相補的な関係にあったが，まるで基盤モデルが強化学習のための効率的な世界モデルになるように，計算機シミュレーションと Monte Carlo 法が新時代の科学の第三の要素になるのかも知れない．

“What distinguishes computer simulation in general from other forms of computation, if such a distinction can be made, is the manner in which the computer is used: instead of merely performing a calculation, the computer becomes the virtual laboratory in which a system is studied - a numerical experiment.” (Rapaport, 2004, p. 3)．

Article Image

計算とは何か

数値実験と LLM とはいずれもシミュレーションに使えるが，用いる形式が違う（数字と文字）．これにより，物理的な用途と社会的な用途とに別れている．この形式の違いを超克するのが機械学習の悲願であるとするならば，計算とはなんだろうか？ Monte Carlo 法とはシミュレーションと計算を架橋する存在であるならば，今後どのような貢献ができるのであろうか？

1.2 歴史

アイデア自体は当然極めてナイーブであり，Laplace の悪魔や気体分子運動論の時代から大変に意識されたが，その最初の非自明な適用はデジタルコンピュータ完成のあとであった (Alder and Wainwright, 1958), (Alder and Wainwright, 1959)．

彼らの関心は，ポロニウムなどの剛体球 (hard-sphere) の系とみれる物質の温度に依存した相転移現象にあった (Alder and Wainwright, 1957)．

(Alder and Wainwright, 1959) では一般の撃力相互作用 (impulsive interaction) に拡張され，event-driven アプローチ (ED-MD scheme) が取られた．1

その他に，粒子に内部構造がない場合は time-driven な Newton 力学のシミュレーションで良いかもしれないが，剛体で体積を持つ場合は Euler 方程式，内部構造を持つ場合は Langrange 方程式のシミュレーションも伴うかもしれない (Rapaport, 2004, p. 4)．

通常の平衡状態の（エネルギー・体積・粒子数一定の）分子系は小正準集団に対応するが，特定の温度条件下での物性を考えたい場合は，運動方程式を修正してシミュレーションする方法もある．この方法の欠点は，物性を再現できても，個々の分子の軌道を正しく模倣しているわけではないという点であるが，そもそもそのカオス的な振る舞いから，正しい軌道を得ることは諦めることが多い．このこともあり，MD 法で用いられる数値積分法は比較的低次元で軽量なものでも十分なのである (Rapaport, 2004, p. 4)．2

セルオートマトンや格子ボルツマン法などの格子ベースの手法の方が計算量は安価であるが，表現力に劣ることになる．

Lanvegin 方程式に基づいた Brownian dynamics などのさらに連続時間ベースの手法もある．

1.3 タイムステップの問題

ここで，optimal scaling と並行な，タイムステップ選択の問題 がやはり中核となることが確認された．

how to choose a good step size has always been an art in the field. (Liu, 2004, p. 183)

タイムステップが系の特徴的な事象をよく捉えるように注意して取る必要がある．特に，系のハミルトニアンが変化しすぎないように離散化誤差を抑える必要がある (Liu, 2004, p. 184)．

特に，凝縮系において系を放置し過ぎると，すぐにハミルトニアンが保存されなくなってしまうため，MD を用いた Monte Carlo 法では極めてタイムステップを細かくすることが必須であることが一番の問題である．

a main problem with MD simulation is the stringent requirement of a small time-step size. (Liu, 2004, p. 189)3

かといって，MD から遊離した MCMC を実行して性能を上げるためには，何かしらの形で遷移核に分布の形状を通知したいところであるが，これは MD 法なしには難しいということになる．4

明らかに分布が特殊な形状をしているのに，ランダムウォークがその方向を見つけるまで待つのでは効率が悪い．

1.4 ポテンシャルの役割

1.4.1 モデル選択の見方

例え量子系のシミュレーションであろうとも，Born-Oppenheimer などの近似を重ねて，原子に働く有効ポテンシャルにのみ注目することにより，古典系のシミュレーションと同様の方法で実行することができ，多くのマクロ的な性質を再現することが出来ることが，(Griebel et al., 2007, p. 17) に詳細に説明されている．

実際の系での粒子間相互作用を決定する際，純粋に Schrödinger 方程式などから理論的に導出するだけでなく，解析的に記述されたポテンシャルをフィッティングしてみてそれが理論や実験に合致するかどうかを見る統計的な方法もとられる (Griebel et al., 2007, p. 27)．

このようなモデル選択の立場に立つのが良い (戸田盛和 et al., 2011, p. 34)．

1.4.2 ポテンシャルの推定

多くの場合ポテンシャルは，粒子間距離や角度，座標などの変数が入っており，これらのパラメータを推定することで探されるが，当然このステップは難しいものである：

The construction of good potentials is still a form of art and requires much skill, work, and intuition. (Griebel et al., 2007, p. 28)

Morse ポテンシャルや Lennard-Jones ポテンシャルがその代表例である．いわば，これらのポテンシャルも模型なのである．半導体のモデリングなどでは，さらに複雑なポテンシャルが必要になる (Griebel et al., 2007, p. 30)．

このように，多体問題をポテンシャル関数によって解く手法は 1980 年代からである (Griebel et al., 2007, p. 30)．

ポテンシャル関数の形によって、適した数値計算アルゴリズムが異なる。特に、長距離力の高速計算のためには、ポテンシャルの性質を巧みに利用したアルゴリズムが必要である。

1.4.3 システム同定としての物理学

What do we mean by “understanding” something? We can imagine that this complicated array of moving things which constitutes “the world” is something like a great chess game being played by the gods, and we are observers of the game. (Feynman et al., 1964)

今まで見てきたように，状態空間モデルは運動学習過程を統一的に説明するモデルであり，最適推定や最適制御の定式化に欠かせない．一方，状態空間モデルに含まれる行列の値は，運動方程式などといった物理的要請から決定できることもあるが，運動適応の学習係数などは未知のパラメータであることが多い．従って，与えられた実験データから状態空間モデルのパラメタを決定する必要がある．このパラメタ推定は制御理論において システム同定 と呼ばれる．冒頭の引用のように，ファインマンは数学や物理を自然界のシステム同定とみなした．(田中宏和, 2019, p. 173)

2 分子動力学法

2.1 剛体円板のシミュレーション

剛体円板のシミュレーションは，衝突イベントを先に予測し，その間を線型に補間すれば良いから，典型的には event-driven で行われる．

2.1.1 アルゴリズムの概略

  • 初期位置と初期速度 \(\{v_i\}_{i=1}^N\subset\mathbb{R}^2\) をランダムに定める．
  • 次のように系を発展させていく： \[ x_i(t)\gets x_i(t_0)+(t-t_0)v_i(t), \tag{1}\] \[ v_i(t)\gets v_i(t_0). \]
  • 衝突イベントの発生時刻 \(t^*\) は次の方程式によって予測できる： \[ d(x_i(t^*),x_j(t^*))=2\sigma. \] 衝突が起こると，衝突の角度に応じて，速度を反転させる： \[ v_i(t^*)\gets v_i(t_0)-\left(\frac{v_{ij}(t_0)^\top x_{ij}(t^*)}{4\sigma^2}\right)x_{ij}(t^*). \]

2.1.2 数値誤差の蓄積

式 (1) が正確な時間発展を記述できている以上，誤差は衝突イベントの処理の際の数値誤差にのみ起因する．

しかし，角度計算の小さな誤差は，大きな軌道の違いを生むため，この数値誤差の蓄積が典型的に起こる．5

2.1.3 アルゴリズムのエルゴード性

\(\mathbb{T}^2\) 上の剛体円板ダイナミクスのエルゴード性（Bernoulli 撹拌性）は (Simányi, 2003) によって示された．

特にエルゴード理論の観点からは，ビリヤード問題とも呼ばれる (Tabachnikov, 2005)．

2.2 滑らかなポテンシャルを持つ小正準集団のシミュレーション

剛体円板のポテンシャルは単関数であった．より現実的な系は，滑らかなポテンシャルによって記述される．

2.2.1 力学

ポテンシャル \(U\) に対して，Newton の運動方程式は \[ m_t\frac{d ^2x_i}{d t^2}=-D_iU(x) \] で表される．

これを一階系 \[ \dot{x}_i=v_i \] \[ \dot{v}_i=-\frac{D_iU(x)}{m_i} \] にすることで，数値的に可解になる．

この系，または Hamiltonian 系の数値計算法については，(Bou-Rabee and Sanz-Serna, 2018) などのレビューも参照．

2.2.2 速度に関する (Verlet, 1967) 法

Hamiltonian 系の数値解法において，最も一般的なアルゴリズムは次のようなものである：

ステップサイズ \(\epsilon>0\) を設定しておく．

  1. 運動量を半分だけ進める： \[ p_i\left(\frac{\epsilon}{2}\right)\gets p_i(0)+\frac{\epsilon}{2}F_i(x(0)). \]
  2. 次のように位置と運動量を更新する： \[ x_i(n\epsilon)\gets x_i((n-1)\epsilon)+\epsilon p_i\biggr(\frac{(n-1/2)\epsilon}{m_i}\biggl), \] \[ p_i\biggr((n+1/2)\epsilon\biggl)\gets p_i\biggr((n-1/2)\epsilon\biggl)+\epsilon F_i(x(n\epsilon)). \]
  3. 残りの運動量を進める．

２のダイナミクスは leapfrog dynamics とも呼ばれる．

運動量の更新を kick，位置の更新を drift と呼ぶこともある．

この手法は，イテレーション毎に \(F\) の評価が１回で済むこと，\(O(\epsilon^2)\) 誤差の性能が知られている．6

この手法は (Störmer, 1907) もすでに用いていたもので，Störmer-Verlet integration とも呼ばれる．

2.2.3 小正準分布のシミュレーション

上述のダイナミクスは，Hamiltonian の等位集合上のみを探索するという意味で，小正準集団のシミュレーションのみが可能である．

すなわち，Boltzmann-Gibbs 分布のシミュレーションが可能でない．\(\left\{U\le H_0\right\}\) 上しか探索できないためである．

単に剛体円板（第 2.1 節）においては \(U\) が単関数であったために，任意の \(H_0\) に対して集合 \(\left\{U\le H_0\right\}\) が変わらなかったため，問題が起こらなかったのである．

2.3 滑らかなポテンシャルを持つ正準集団のシミュレーション

2.3.1 Langevin 動力学

従って，系が外界とエネルギーを交換するような系（正準集団）を考えたい場合は，単に Hamilton 方程式を数値的に解くだけでは不十分である．

そこで，Hamilton 方程式の運動量の項に，\(\beta\) に依存した係数を持つランダムな項を加えて得る (underdamped) Langevin 動力学を考えることができる：7 \[ dx_i=m_i^{-1}p_idt, \] \[ dp_i=-\nabla U(x)dt-\gamma m_i^{-1} p_idt+\sqrt{2\gamma\beta^{-1}}dB_i. \]

\(\gamma\) は粘性係数，\(\beta\) は逆温度である．

2.3.2 シミュレーション法

次のようにして underdamped Langevin 動力学を数値的に解くことができる：8

\(\epsilon>0\) をタイムステップとする．

  1. 確率項を消去した Hamilton 系をシミュレーションし，\(\epsilon>0\) だけ進める．
  2. 残る Ornstein-Uhlenbeck 過程 \[ dp_i=-\gamma m_i^{-1} p_idt+\sqrt{2\gamma\beta^{-1}}dB_i \] を解き，\(\epsilon>0\) だけ進める．

以上を繰り返す．

2.3.3 エルゴード性

カップリングアプローチによるエルゴード性証明は (Eberle et al., 2019) により初めて行われた．

2.3.4 過減衰 Langevin 法

過減衰極限 \(\gamma\to\infty\) を取ることで，１本の方程式にまとまる：9 \[ dy_i=-\nabla U(y)dt+\sqrt{2\beta^{-1}}dB_i, \] \[ y(t):=x(\gamma t). \]

これを 過減衰 Langevin 動力学 または対称 Langevin 拡散という．統計力学では Brownian 動力学 とも呼ばれる (Rossky et al., 1978)．

これに対して，元の underdamped Langevin 拡散は，力学的 Langevin 動力学 (kinetic Langevin dynamics), ２階の Langevin 過程などとも呼ばれる．10

2.3.5 Langevin 力学によるサンプリング

統計学と機械学習では，特に (Roberts and Tweedie, 1996) 以来，上述の過減衰 Langevin 動力学の方がポピュラーである．

しかし，underdamped Langevin 動力学を用いた方が，サンプリング法の収束が速くなることが知られている (Dalalyan and Riou-Durand, 2020)．

3 Hybrid Monte Carlo

3.1 はじめに

MD 法はサンプリングに向いているわけではない．また，Monte Carlo 法も，提案核の設計によって効率は大きく変わる．

そこで，MD 法のように Hamiltonian からの知識を提案核に取り入れた Monte Carlo 法が HMC (hybrid Monte Carlo) (Duane et al., 1987) として提案された．11

これは混合モンテカルロ法ともいう：

混合モンテカルロ／分子動力学法：モンテカルロ法を用いた構造分布の生成において，試行に用いる構造をニュートンの運動方程式など決定論的手法に従って用意する手法．分子凝集系などの複雑な系に対して，効率的に構造空間探索を行えると期待される．(栗﨑 and 田中, 2022)

3.2 Hamiltonian Monte Carlo

後に，統計力学的な背景を持たないサンプリング問題についても，Hamilton 系から示唆された提案分布を用いる方が MCMC として性能が良いことが発見され，(Neal, 1993) 以降機械学習でも広く用いられるようになった．

今日ではこれらの手法は HMC (Hamiltonian Monte Carlo) 法 (Neal, 2011) として機械学習・統計学界隈に普及しており，双方の分野からアルゴリズム改善のアイデアが盛んに生まれている．12

3.3 MD に提案をもらう Monte Carlo 法

提案は MD によって示唆されたものの方が効率が良くなる場合がある：

The advantage of the MD proposal is that the resulting MCMC moves follow the dynamics of the target distribution more closely. (Liu, 2004, p. 184)

これは結局，背後の物理学的な本質を捉えているためとも言える：

A major advantage of molecular dynamics simulation in physical systems is its reliance on basic physics principles (e.g. Newton’s equation), which has been shown by nature to work well. (Liu, 2004, p. 189)

References

Alder, B. J., and Wainwright, T. E. (1957). Phase Transition for a Hard Sphere System. The Journal of Chemical Physics, 27(5), 1208–1209.
Alder, B. J., and Wainwright, T. E. (1958). Molecular dynamics by electronic computers. In I. Prigogine, editor, Proceedings of the international symposium on transport processes in statistical mechanics. Held in brussels, august 27-31, 1956, page 97.
Alder, B. J., and Wainwright, T. E. (1959). Studies in Molecular Dynamics. I. General Method. The Journal of Chemical Physics, 31(2), 459–466.
Bou-Rabee, N., and Sanz-Serna, J. M. (2018). Geometric integrators and the hamiltonian monte carlo method. Acta Numerica, 27, 113–206.
Bressloff, P. C. (2021). Stochastic processes in cell biology,Vol. 41. Springer Cham.
Dalalyan, A. S., and Riou-Durand, L. (2020). On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli, 26(3), 1956–1988.
Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987). Hybrid monte carlo. Physics Letters B, 195(2), 216–222.
Eberle, A., Guillin, A., and Zimmer, R. (2019). Couplings and quantitative contraction rates for Langevin dynamics. The Annals of Probability, 47(4), 1982–2010.
Faulkner, M. F., and Livingstone, S. (2024). Sampling Algorithms in Statistical Physics: A Guide for Statistics and Machine Learning. Statistical Science, 39(1), 137–164.
Feynman, R. P., Leighton, R. B., and Sands, M. (1964). The feynman lectures of physics. In,Vol. I. Addison-Wesley.
Griebel, M., Zumbusch, G., and Knapek, S. (2007). Numerical simulation in molecular dynamics: Numerics, algorithms, parallelization, applications,Vol. 5. Springer Belin, Heidelberg.
Leimkuhler, Ben, and Matthews, C. (2015). Molecular dynamics: With deterministic and stochastic numerical methods. Springer Cham.
Leimkuhler, Benedict, and Reich, S. (2005). Simulating hamiltonian dynamics. Cambridge University Press.
Liu, J. S. (2004). Monte carlo strategies in scientific computing. Springer New York.
Neal, R. M. (1993). Probabilistic Inference using Markov Chain Monte Carlo Methods. Department of Computer Science, University of Toronto. Retrieved from https://glizen.com/radfordneal/review.abstract.html
Neal, R. M. (2011). In S. Brooks, A. Gelman, G. Jones, and X.-L. Meng, editors, pages 113–162. Chapman; Hall/CRC.
Pavliotis, G. A. (2014). Stochastic processes and applications: Diffusion processes, the fokker-planck and langevin equations,Vol. 60. Springer New York.
Peters, E. A. J. F., and de With, G. (2012). Rejection-free monte carlo sampling for general potentials. Physical Review E, 85(2).
Rapaport, D. C. (2004). The art of molecular dynamics simulation. Cambridge University Press.
Roberts, G. O., and Tweedie, R. L. (1996). Exponential convergence of langevin distributions and their discrete approximations. Bernoulli, 2(4), 341–363.
Rossky, P. J., Doll, J. D., and Friedman, H. L. (1978). Brownian dynamics as smart Monte Carlo simulation. The Journal of Chemical Physics, 69(10), 4628–4633.
Simányi, N. (2003). Proof of the boltzmann-sinai ergodic hypothesis for typical hard disk systems. Inventiones Mathematicae, 154(1), 123–178.
Stoltz, G. (2021). Computational statistical physics and hypocoercivity.
Störmer, C. (1907). Sur les trajectoires des corpuscules électrisés dans l’espace. Applications à l’aurore boréale et aux perturbations magnétiques. Radium (Paris), 4(1), 2–5.
Tabachnikov, S. (2005). Geometry and billiards,Vol. 30. American Mathematical Society.
Verlet, L. (1967). Computer "experiments" on classical fluids. I. Thermodynamical properties of lennard-jones molecules. Phys. Rev., 159, 98–103.
戸田盛和, 斎藤信彦, 久保亮五, and 橋爪夏樹. (2011). 統計物理学,Vol. 5. 岩波書店.
栗﨑, and 田中. (2022). 用語解説. 生物物理, 62(4), 250–250.
田中宏和. (2019). 計算論的神経科学：脳の運動制御・感覚処理機構の理論的理解へ. 森北出版.

Footnotes

  1. とは (Peters and de With, 2012) でも言及されている．当時の計算機では，500個の粒子を扱うことが限界であったことが興味深い．↩︎

  2. 内部構造がある分子をシミュレーションする場合はよりやわらかい相互作用を考える必要があり，その際は積分は高次元になる上に，内部の運動が高速であるためにタイムステップも細かくする必要が出てくる．さらに拘束条件が存在する場合は，積分法よりも高い精度で取り扱う必要があり，特別な注意を要する．↩︎

  3. “For example, the protein folding process takes about \(10^{-3}\) seconds in nature. A proper MD simulation of such a process needs a step size of order \(10^{-12}\) and will take about \(10^6\) days using a current computer.”↩︎

  4. For example, if the system of interest consists of closely packed particles, a random proposal for moving a particle is most likely rejected because the proposed new position has been partially occupied by others. (Liu, 2004, p. 189)↩︎

  5. (Faulkner and Livingstone, 2024, p. 17) 4.3.1 節も参照．↩︎

  6. (Faulkner and Livingstone, 2024, p. 18) 4.3.2 節，(Bou-Rabee and Sanz-Serna, 2018) も参照．↩︎

  7. 詳しくは (Ben Leimkuhler and Matthews, 2015) 6.3.2 節，(Pavliotis, 2014, p. 181)，(Bressloff, 2021, p. 46) も参照．↩︎

  8. (Benedict Leimkuhler and Reich, 2005) に詳しい．↩︎

  9. この過程は，(Pavliotis, 2014) 第 6.5 節，そして (Stoltz, 2021) 参照．↩︎

  10. (Dalalyan and Riou-Durand, 2020) に倣った．↩︎

  11. 特に格子場の理論において，Fermion 自由度が存在する場合のシミュレーションを問題として扱っていた．↩︎

  12. (Faulkner and Livingstone, 2024, p. 19) 4.4 節も参照．↩︎