分子動力学シミュレーション（MD）は化学のシミュレーションの中でもかなりメジャーな方法です。タンパク質のダイナミックな動きをシミュレーションしたMDを見た方も多いのではないでしょうか？
MDは古典力学に則るだけと思われがちですが、精度よくシミュレーションしたり、温度が一定になるように計算するには工夫が必要など、奥が深い分野です。
このページではMDでよく用いられる速度ベルレ法（NVE一定）と能勢・フーバー法（NVT一定）のうち、速度ベルレ法のアルゴリズムを確認し、実際に単原子分子に対して分子動力学シミュレーションを行った結果をご紹介します。
- NVEが一定となる各粒子の運動方程式
- NVEが一定の分子動力学シミュレーションアルゴリズム（速度ベルレ法のアルゴリズム）
- NVEが一定のシミュレーション結果
内部エネルギー$E$ではなく、温度$T$を一定とする能勢・フーバーのアルゴリズムは↓のリンクで解説しています。併せてご覧ください。
NVE一定のシミュレーション
MDで最も単純なアルゴリズムがNVE（粒子数、体積、内部エネルギー）が一定のシミュレーションです。結論を述べると、速度ベルレ法を使うことで精度よく内部エネルギーを保存させながらシミュレーションを行うことができます。
皆さんも高校生の時にばねによる振動では運動エネルギーとポテンシャルエネルギーの和は保存されることを学んだのではないでしょうか？これと同じように、内部エネルギー$E$が一定のシミュレーションは単純に$ma=F$を解くことで得ることができます。
内部エネルギー保存はどんなときでも成り立つの？
粒子に働く力$F$がポテンシャル$V(r)$の微分で得られるとき（分子間力、静電引力など）のときに内部エネルギーは保存されるよ。
NVE一定の運動方程式
つまり、NVT一定のMDを行うには、単純に運動方程式$ma=F$を各粒子ごと（$i=1,2,\cdots,N$）に解けばよいことがわかります。
$\frac{\text{d}r_i}{\text{d}t}=v_i$（位置の変化量=速度×$\Delta{t}$）
$\frac{\text{d}p_i}{\text{d}t}=F_i$（運動方程式 $ma=F$で運動量$p=mv$とする）
微分使われてるし、運動量$p$を使われるとあまりなじみないな。。
位置の変化量の式と$ma=F$で運動量$p=mv$を利用しているだけで、複雑ではありませんよ。
この方程式では内部エネルギー$E=\sum_i\frac{p_i^2}{2m_i}+\sum_{i<j}V(r_{ij})$が保存量として保存されます。
この運動方程式を100個や1000個の粒子に対して一つ一つ解いていくことで分子動力学シミュレーションを行うことができるのです。
単原子状分子のMDをグラフィック化したものをご紹介しますね。
分子動力学シミュレーションをグラフィック化したものです。（VMDで可視化しています）
アルゴリズム
NVE一定の運動方程式を実現する最も単純なアルゴリズムは次ではないでしょうか？
- $p_i(t+\Delta{t})=p_i(t)+F_i*\Delta{t}$
- $x_i(t+\Delta{t})=x_i(t)+\frac{p_i(t+\Delta{t})}{m_i}*\Delta{t}$
①によって、運動量を時間$t$から$t+\Delta{t}$へ発展させ、②で位置を時間$t$から$t+\Delta{t}$へ発展させます。この①と②をループさせることで大量の粒子や分子を時間発展させることができます。
力$F_i$は粒子$i$がほかの粒子から感じる力$f_i$を足し合わせたものだね。
このアルゴリズムでもプログラムとしては問題なく動くのですが、実は精度がよくありません。
そこで、もっとよい精度を出すアルゴリズムが速度ベルレ法です。次のようにひと手間加えます。
- $p_i(t+\frac{\Delta{t}}{2})=p_i(t)+F_i*\frac{\Delta{t}}{2}$
- $x_i(t+\Delta{t})=x_i(t)+\frac{p_i(t+\Delta{t}/2)}{m_i}*\Delta{t}$
- $p_i(t+\Delta{t})=p_i(t+\frac{\Delta{t}}{2})+F_i*\frac{\Delta{t}}{2}$
運動量$p_i$の発展を二回にわけているんだね。
理論的な説明はこのページでは省きますが、単純なアルゴリズムに比べ、速度ベルレ法では保存量である内部エネルギーの保存性が格段に上がります。
実際のシミュレーション結果(NVE一定)
実際に1000粒子系について、NVE一定のシミュレーションを行い保存量$E$の変位を調べました（タイムステップ$\Delta{t}=0.005$）。
保存量のグラフなので、より一定に近いほうがよいシミュレーションということになります。
赤い線がオイラー法で、青い線が速度ベルレ法です。
見ての通りですが、オイラー法では内部エネルギー$E$が$0.1$程度の振れ幅で増減してしまっていますが、ベルレ法ではそれよりも振れ幅が小さいことがわかります
速度の時間変化アルゴリズムにひと手間加えただけで、青い線はずっと一定で直線に見えるね！
もちろん、シミュレーションである以上、速度ベルレ法でも内部エネルギー$E$は厳密には一定にはならず、おおよそ$0.005$程度の幅で増減を繰り返していることがわかります。
同じタイムステップ$\Delta{t}=0.005$でもオイラー法では$0.1$の振れ幅、ベルレ法では$0.005$の振れ幅とかなり振れ幅が異なってくるのです。
ベルレ法は簡便ながら、内部エネルギーの保存を高いレベルで実現できる方法です。
まとめ
このページでは次のことを確認しました。
- NVEが一定となる各粒子の運動方程式
- NVEが一定の分子動力学シミュレーションアルゴリズム
- NVEが一定のシミュレーション結果
NVE一定の分子動力学シミュレーションは分子動力学シミュレーションのなかでもかなり基本的でよくつかわれる手法なので、みなさんも一度試してみてください。
