太陽系の惑星運動シミュレーションと分子動力学計算のスケールを比較

ふとしたきっかけで、地球と火星が数十億年後に衝突する可能性があるよー、というシミュレーション研究が数年前に話題になっていたことを思い出し、元記事を見つけました。

35億年後、地球が火星と衝突する？パリ天文台が試算

論文で主張している結果については上記記事を見てほしいですが、簡単にまとめておくと以下のような感じです。

• 今後50億年間の太陽系惑星（水星、金星、地球、火星、木星、土星、天王星、海王星、冥王星）の軌道をコンピュータシミュレーションで解析
• 2501本のシミュレーションを実行したところ、地球と火星が実質的に衝突したシミュレーションが１本得られた
• この地球と火星が衝突したシミュレーションの途中から火星の位置を微妙にずらして200本のシミュレーションを行ったところ、地球が水星、金星、火星と衝突するケースが出現した

多体系のカオス理論（わずかでも初期条件のが異なっていれば、時間の経過とともにその差は拡大する）を考慮すれば、数十億年後の太陽系惑星の運動を正確にシミュレートすることは不可能であるため、このように多数のシミュレーションを実行して結果を解析する必要があるわけです。

（実際に地球と他の惑星が衝突するとしても数十億年後の話ですので、特に私たちの生活に影響がある話ではありませんね）

惑星運動シミュレーションのスケール

私は１個１個の原子が見えるレベルでのシミュレーション（分子動力学シミュレーション）を本業としていますので、このような惑星規模のシミュレーション論文を読んで非常に刺激を受けました。

まず面白いと感じたのは使っている単位の違いです。

この論文中の図の横軸は「Time (Myr)」とあり Fig.1 では 0-5000 までの値をとっているのですが、この Myr はなんだろう？としばらく考え込んでしまいました。
しばらく考えたところ、50億年のシミュレーションを実行したという記述から Myr = Mega Year(100万年）のことだと気づき、改めて分子レベルとは時間スケールが違うことを実感しました。

数値積分の積分刻みは 2.5×10^-2 years (つまり 0.025年 ＝ 9.13日 ＝ 79万秒）であるというのも時間スケールの違いを感じさせられます。

50億年のシミュレーションということは、数値積分を2000億ステップ行う必要がありますが、この実行には Xeon 5472 CPU で 2500 時間〜３か月半かかったようです。（Supporting Information参照）

なお積分刻みが9.13日ということは、およそ１か月に１度（正確には27日7時間43分）地球の周りを公転している月の運動はまともに積分することは不可能になります（月の公転一周期を積分刻み３回で数値積分するのは誤差が大きすぎる）。
そのため月の影響は一般相対論を考慮した上で平均場的に取り込んでいるようです。

また多数のシミュレーション実行を行うための初期状態として、地球の初期座標を3.8cmずつずらして実行しているというのも興味深く感じました。

分子シミュレーションの場合

分子シミュレーションであれば、水素原子を含む結合の伸縮運動を凍結させるSHAKE法を適用する場合、積分刻みは 1-2 fs（フェムト秒、10^-15秒）程度を使います。

もし水素原子を含む結合の振動も考慮するのであれば積分刻みは 0.2 fs 程度とより細かくとる必要があります。
これは最も周期が短い水素原子を含む結合の振動周期が 10 fs 程度であり、この振動周期の運動を正しく数値積分する必要が出てくるためです。
詳細は下記の本のＦ１５節に解説されていますので、興味のあるかたは読んでみてください。

すぐできる分子シミュレーションビギナーズマニュアル (KS化学専門書)

水素原子を含む結合が振動している様子は下の動画で見ることができます。
動画の右下にシミュレーション時間を ps (ピコセカンド、1 ps = 1000 fs) で表記しています。（このシミュレーションでは積分刻みは 0.2 fs で実行しました）

特に目的がない限り、通常はSHAKE法を適用して水素原子を含む振動を凍結し、2 fs程度の積分刻みを適用します。

最近の標準的な分子シミュレーション研究では ns （ナノ秒, 10^-9秒）〜 μs （マイクロ秒, 10^-6秒）のシミュレーションを実行しますので、 2 fs の積分刻みであれば 100万〜10億 step の数値積分を計算することになります。

惑星軌道シミュレーションと分子シミュレーションの単位を比較

惑星軌道シミュレーションと分子シミュレーションの単位を比べてみました。

.	惑星軌道シミュレーション	分子シミュレーション
数値積分刻み	10^6 sec （数日）	10^-15 sec（フェムト秒）
解析対象とする現象の時間スケール	10^16 sec（10億年）	10^-9 〜 10^-6 sec（ナノ秒〜マイクロ秒）
高速化のために無視する運動	月の公転運動	水素原子を含む振動

積分刻みは21桁異なり、対象とする現象のスケールも 億年：ナノ〜マイクロ秒 と大きく異なっています。

これほど対象の時間スケールは異なっていますが、分子動力学計算のSHAKE法と同様に積分刻みを長くとるための工夫（月の公転の平均場近似）がされているという点は極めて興味深く感じました。