Verlet法の数値計算例（C言語）

プログラムソースコードの解説

今回行ったシミュレーションは、Verlet 法を用いて、質量１の粒子をランダムに多数ばら撒き、その後の動きをシミュレーションするものである。粒子間には、互いの距離の２乗に反比例する反発力が働くとする。

プログラムは、一定時間おきに粒子の座標を記録するが、それだけでは、系が安定状態になったのかわかりにくいので、運動エネルギーの値を見て、系が安定状態になったか判断する。そのために、プログラムでは、速度の２乗値から運動エネルギーを計算し、その値も記録する。

ただし、本シミュレーションでは、平衡化という作業を行なっていないので不十分である事に注意する必要がある。つまり、系の設定温度になるように速度を調整する事と、運動量がゼロになるように粒子の速度を調整する必要がある。しかし、今回は、大まかな分子動力学法の流れを見るために、平衡化は行わないとする。

ソースコードはC言語で書かれているが、verlet1.cをダウンロードし、適切なソフトを使えば見られる。

gnuplotを用いた出力結果

プログラムを実行し、処理が終了すると、幾つかのファイルが生成される。ファイルから、gnuplot を用いてグラフを出力すれば、シミュレーション結果が見られる。

上図は、energy.dat をプロットしたものであり、横軸が time 値、縦軸が k 値である。つまり、運動エネルギーが、時間に対してどのように変化するか見るためのグラフである。グラフを見るとわかる通り、横軸の値が約 1000 のところまで、運動エネルギーの値は急上昇している。これは、初速度 0 で、ランダムに配置された粒子が、他の粒子との相互作用力により動き出し、速度を持ち始めたからである。その後は、大体一定値をとるようになっており、系は安定な状態へ向かっているのがわかる。