分子動力学（MD）の数値解法

分子動力学は、連続時間の運動方程式を離散時間へ写像し、その誤差と安定性を制御しながら、原子運動の統計性とダイナミクスを再現する計算法である。数値解法の選択は、再現したいアンサンブル（NVE/NVT/NPTなど）と、必要な時間・長さスケールによって決まる。

参考ドキュメント

1. L. Verlet, Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules (1967) https://link.aps.org/pdf/10.1103/PhysRev.159.98
2. 吉井範行, 分子動力学法 I（応用物理 75巻6号, 2006） https://www.jstage.jst.go.jp/article/oubutsu/75/6/75_718/_article/-char/ja/
3. 理研数理, 国産の分子動力学ソフトウェア「GENESIS」の産業界における利用をサポートするサービスの提供開始（2021） https://www.riken-suuri.jp/news-release/20211012.html

1. 連続時間の方程式

1.1 ニュートン方程式（古典MD）

粒子 $i$ の位置 $r_i$、速度 $v_i$ に対して

$$m_i{\ddot{\mathbf{r}}}_i={\mathbf{F}}_i({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)=-{\mathrm{\nabla }}_{{\mathbf{r}}_i}U({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)$$

である。ハミルトニアン形式では

$$H(\mathbf{p},\mathbf{q})=\sum _i\frac{{\mathbf{p}}_i^2}{2m_i}+U(\mathbf{q})$$

であり、数値積分では「このハミルトニアン流れを、有限の刻み幅で近似的に追跡する」ことになる。

1.2 確率微分方程式（Langevinなど）

熱浴を確率過程で表すと

$$d\mathbf{r}=\mathbf{v}dt,md\mathbf{v}=\mathbf{F}dt-\gamma m\mathbf{v}dt+\sqrt{2\gamma m{k}_{B}T}d{\mathbf{W}}_t$$

となる。ここで $d{W}_t$ はウィーナー過程である。

2. 数値解法に求められる性質

MDでは短時間の「精密な軌道」より、長時間での「統計量の正しさ」と「保存則の破綻のしにくさ」が要点になりやすい。したがって以下が重要である。

• 安定性：高速振動（最短周期）に対して発散しないこと
• 長時間挙動：エネルギーや運動量が系統的にドリフトしにくいこと
• 時間反転対称性：速度反転で同じ軌道を戻れる構造を保つこと -（ハミルトニアン系では）シンプレクティック性：位相空間体積を保つ離散写像であること

これらは「時間刻みを小さくすれば常に解決」ではなく、同じ刻み幅でも積分器の構造で長時間誤差が大きく変わる。

3. 基本の時間積分：Verlet系

3.1 Verlet（位置Verlet）

$${\mathbf{r}}_{n+1}=2{\mathbf{r}}_n-{\mathbf{r}}_{n-1}+\frac{\mathrm{\Delta }{t}^2}{m}\mathbf{F}({\mathbf{r}}_n)$$

であり、速度を陽に持たない。2次精度であり、時間反転性を持つ。

3.2 速度Verlet（velocity Verlet）

力評価を 1 回/step に保ちつつ速度も更新できる形で、最も広く使われる。

$${\mathbf{r}}_{n+1}={\mathbf{r}}_n+{\mathbf{v}}_n\mathrm{\Delta }t+\frac{\mathrm{\Delta }{t}^2}{2m}\mathbf{F}({\mathbf{r}}_n)$$$${\mathbf{v}}_{n+1}={\mathbf{v}}_n+\frac{\mathrm{\Delta }t}{2m}[\mathbf{F}({\mathbf{r}}_n)+\mathbf{F}({\mathbf{r}}_{n+1})]$$

3.3 Leapfrog（スタッガード速度）

速度を半ステップずらして更新する：

$${\mathbf{v}}_{n+\frac{1}{2}}={\mathbf{v}}_{n-\frac{1}{2}}+\frac{\mathrm{\Delta }t}{m}\mathbf{F}({\mathbf{r}}_n),{\mathbf{r}}_{n+1}={\mathbf{r}}_n+\mathrm{\Delta }t{\mathbf{v}}_{n+\frac{1}{2}}$$

Verlet系と等価な見方ができる。

3.4 代表的な積分器の比較

積分器	次数（大域）	時間反転	シンプレクティック（NVE）	力評価/step	特徴	注意点
前進Euler	1	×	×	1	最小構成	エネルギー発散しやすい
Verlet / velocity Verlet / leapfrog	2	○	○	1	長時間安定・標準	刻み幅は高速振動で制限
Runge–Kutta（典型RK4）	4	×	×	4	短時間精度は高い	長時間で保存則が崩れやすい
予測子修正子（Gear系）	高次	△	×	複数	摩擦・速度依存力で使われた歴史	ハミルトニアン系ではドリフトが出やすい
分割（Trotter/Strang）	設計次第	設計次第	設計次第	設計次第	熱浴・圧力制御と相性が良い	分割順序が統計誤差を左右

4. 時間刻み $\Delta t$ の決め方

4.1 線形安定性（調和振動子の目安）

調和振動子の角振動数 $\omega$ に対して、leapfrog/Verlet系は概ね

$$\mathrm{\Delta }t<\frac{2}{\omega }$$

が安定性の目安となる。

4.2 物理的な「最速運動」に合わせる

• 軽元素（特にH）を含む結合伸縮は高速であり、これが $\Delta t$ を厳しくする
• 拘束（SHAKE/RATTLE）で X–H 伸縮を固定すると、$\Delta t$ を大きくできる場合がある
• 電子状態計算を含む第一原理MD（AIMD）は力評価が高価であり、$\Delta t$の選択は計算コストと物理精度の折衷になる

5. 拘束条件：SHAKE / RATTLE

5.1 何を拘束するか

結合長などの拘束を

$$g_k(\mathbf{r})=0(k=1,\dots ,M)$$

で課すと、運動方程式はラグランジュ未定乗数 ${\lambda }_k$ を用いて

$$m_i{\ddot{\mathbf{r}}}_i={\mathbf{F}}_i+\sum _k{\lambda }_k{\mathrm{\nabla }}_{{\mathbf{r}}_i}{g}_k$$

となる。ここでの数値解法の要点は「毎ステップ、拘束を許容誤差内で満たすように補正する」ことである。

5.2 SHAKE（位置拘束）とRATTLE（速度まで整合）

• SHAKEは位置更新後、拘束違反を反復的に補正する
• RATTLEは位置だけでなく速度にも拘束整合を課し、温度や運動量に関わる量の破綻を抑える

拘束の反復許容誤差が緩すぎると、見かけの温度やエネルギーのドリフトを誘発する。

6. 温度制御（NVT）と確率積分

6.1 Langevin熱浴：分割積分（BAOABなど）

Langevin方程式はハミルトニアン流れ＋減衰＋ランダムキックに分けられるため、分割（operator splitting）が自然である。分割順序の違いは、同じ$\Delta t$でも配置分布の誤差に大きく影響し得る。

6.2 Nosé–Hoover（連鎖）と拡張系の可逆積分

Nosé–Hoover系は、追加自由度（熱浴変数）を含む拡張ハミルトニアン（拡張力学系）として扱い、可逆な分割積分で時間発展させるのが典型である。連鎖（chains）を用いることで、単一熱浴の非エルゴード性を緩和する設計がなされる。

温度制御では、時定数（coupling time）を短くしすぎると物理ダイナミクスを歪め、長すぎると熱平衡化が遅くなるため、観測量に対して妥当な帯域で設定する。

7. 圧力制御（NPT）とセル自由度

等方体積変化やセル行列の自由度を導入して

• 等方NPT（体積のみ）
• Parrinello–Rahman型（セル形状も） などを扱う。ここでも拡張系の分割積分が基本であり、温度・圧力の同時制御では「正しいアンサンブルを与える積分器（例：MTK系）」が選択される。

8. 多重時間刻み（MTS）：r-RESPA

相互作用を「速い力」と「遅い力」に分割できるとき

$$\mathbf{F}={\mathbf{F}}_{\text{fast}}+{\mathbf{F}}_{\text{slow}}$$

内側刻み $\delta t$ で fast を、外側刻み $\Delta t$（$=n\delta t$）で slow を評価することで、計算量を削減できる。RESPA系は可逆性と安定性を保つように設計されるが、時間スケール分離が不十分な場合、共鳴（resonance）により不安定化することがある。

9. 長距離相互作用：EwaldとPME

周期境界条件下のクーロン相互作用は単純なカットオフで誤差が大きくなりやすい。Ewald和は実空間と逆空間に分割して収束させ、Particle Mesh Ewald（PME）は逆空間計算をFFTで高速化して

• 系サイズ $N$ に対して概ね $NlogN$ スケール の計算を実現する。導電・誘電応答や欠陥電荷の扱いでは、この選択がエネルギーや力の系統誤差に直結する。

10. 数値誤差の診断

10.1 NVE（保存則の検査）

• 全エネルギー $H(t)$ のドリフト：平均的な傾きと揺らぎを分けて評価する
• 線運動量・角運動量：境界条件や拘束・熱浴で破れていないか確認する

10.2 NVT/NPT（分布の検査）

• 速度分布がMaxwell分布に近いか
• 温度・圧力の自己相関時間が設定と整合しているか
• RDF、MSD、時間相関関数などの統計量が、Δtやサーモスタット設定に対して不変か（感度解析）

10.3 連続系の不連続が生む数値人工物

• カットオフで力が不連続になるとエネルギー保存が破れやすい
• 近接リスト更新間隔が長すぎると力の取りこぼしが出る
• 拘束反復の収束判定が緩いと温度が系統的にずれる

11. 国内動向の一例：大規模計算とアルゴリズム選択

超並列環境では、単に計算機が速いだけでなく、PMEや拘束、分割積分、近接探索などのアルゴリズムと実装の整合が性能と再現性の両方を左右する。国産MDコードの超並列展開（例：富岳での大規模計算）に関する公開情報もあり、計算規模がアルゴリズム選択に影響することが示唆される。

まとめ

MDの数値解法は、運動方程式の時間積分だけでなく、拘束条件、熱浴・圧力制御、長距離相互作用、MTSといった周辺要素を含む「離散時間力学系の設計」である。目的の物性（統計量・ダイナミクス）に対して、保存則・アンサンブル・安定性・計算量のバランスが取れる積分器と設定を選ぶことが肝要である。