シミュレーションプログラムの作成

シミュレーションプログラムの流れ

分子動力学シミュレーションの大まかな流れを図示すると、 以下のようになります。

分子動力学シミュレーションの流れ

「位置の更新、速度の更新」、「力の計算」、「温度、圧力などの計算」 の部分をメインループと呼びます。 このメインループを実行するとΔt だけ時間が進みます。 Δt はシミュレーションの時間刻み幅で、 通常は0.5〜2.0(fs)程度の値を使用します。 このメインループを数万回〜数百万回実行することで、 数十(ps)〜数(ns)のシミュレーションを実行することが出来ます。

速度ベルレの方法

ここでは方法としてベルレ使するえますΔt  したときのベルレ
                               2
r(t + Δt ) = r(t) + Δtv (t) + Δt-F (t)
                             2
(1)

v(t + Δt) = v (t) + Δt- {F (t) + F (t + Δt )}
                    2
(2)

となりますここではによりm  = 1  としていますメインル 」、「ベルレってくと

      v′(t + Δt ) = v(t) + Δ2tF (t)       速 度 の 更V (1)
   r (t + Δt ) = r(t) + Δtv ′(t + Δt )  位 置 の 更V
                   ↓
    r (t + Δt )か  ら F (t + Δt )を 計Z     力 の 計Z                   ↓
              ′           Δt
v (t + Δt ) = v (t + Δt ) + 2-F (t + Δt ) 速 度 の 更V(2)

となります

粒子の位置を更新した場合の周期境界処理

速度が一定、力が常にゼロの状態で粒子が動くプログラムは以下の通りです。

main()
{
  int i,j,k;
  int step,total_step;
  double side = 9;
  double cd[N*3],vl[N*3];
  double dt = 3.0;

  cd[0] = 4.5; cd[1]  = 3.0; cd[2]  = 3.0;
  cd[3] = 4.5; cd[4]  = 3.0; cd[5]  = 6.0;
  cd[6] = 4.5; cd[7]  = 6.0; cd[8]  = 3.0;
  cd[9] = 4.5; cd[10] = 6.0; cd[11] = 6.0;
  srand(1);
  for(i = 0; i < N*3; i++){
    vl[i] = ((double)rand()/RAND_MAX-0.5)*.1;
  }
  total_step = 100;
  for(step = 1; step <= total_step; step++){
    for(i = 0; i < N*3; i++){
      cd[i] += vl[i]*dt;
    }
    output_file(cd,side);
  }
}

プログラムの全体はmd_program6_1.cです。 このプログラムでは4つの粒子が一定の速度で移動し、 各ステップの粒子の位置情報を cdview用のファイ ルに出力します。ファイル名はcd0000.dat〜cd0099.datです。

> gcc md_program6_1.c -o md_program6_1
> md_program6_1
> cdview3 cd0000.dat

上記のようにプログラムをコンパイル、 実行した後にcdviewを実行し、'c'キーを押すと粒子が動く様子がわかります ('z'で戻る)。 cdviewについては初期配置の設定 も参考にしてください。

演習

上記のプログラムでは周期境界条件処理を行っていないため、 粒子はシミュレーションセルの外側にはみ出る。 if文 を使用して 周期境界条件 処理を行い、粒子がはみ出さないようにプログラムを変更せよ。

MDプログラムの作成

以上のことを踏まえて、 分子動力学シミュレーションのプログラムをC言語風に書くと。

main()
{
  (変数の初期化)
  (位置の初期化)
  (速度の初期化)
  (温度の補正)

  for(メインループの回数){

   (速度の更新(1))
   (位置の更新)

   (周期境界条件処理)

    for(i = 0; i < 粒子数; i++){
      for(j = i+1; j < 粒子数; j++){
        (力の計算)
      }
    }

   (速度の更新(2))
   (温度の計算)
   (温度の補正)

  }
}

となります。

(変数の初期化)は変数の無次元化
(位置の初期化)は初期位置の設定
(速度の初期化)は初期速度の設定
(力の計算)は力の計算
(位置の更新)、(速度の更新)は数値積分の実行
に詳しく書いてあります。

(周期境界処理)は周期境界条件 を実現するため処理です。

温度の補正

上記のようなMDプログラムではエネルギーが保存されるため、 大抵の場合において温度は設定温度からずれていきます。 そのため、設定温度でのシミュレーションを行うためには、 適当なステップ間隔で温度を補正する必要があります。

温度を補正するためには、現在の温度を設定温度に合わせる処理を行います。 具体的には全粒子の速度に現在の温度と設定温度から計算される定数をかけることになります。 詳しくは 「初期速度の設定」-「温度の補正」 を参照して下さい。

温度の補正はシミュレーションのステップが500以下の場合のみ、 50ステップに一回の割合で行うようにして下さい。 50ステップに一回、 というような処理を行う場合は割り算の余りを求める%演算子を使用します。

初期力

力の初期値は全てゼロにします。 また、力は各ステップで常にゼロから計算するので、 力の計算の前に力の変数全てにゼロを代入する必要があります。 ポテンシャルエネルギーを計算する場合も同様に各ステップでゼロにする必要があります。

レポート課題

これまでの結果を元に、 アルゴンの分子動力学シミュレーションプログラムを完成させ、 以下の事項についてのレポートを作成せよ。 レポートにはシミュレーションのプログラムも印刷して、添付すること。

  1. 速度ベルレの式を導出せよ。
  2. 各ステップで運動エネルギーとポテンシャルエネルギーを計算し、 それらの和からトータルエネルギーの時間変化を求め図示せよ。
  3. トータルエネルギーが500ステップ以降で一定となることを確認し、 その標準偏差を求めよ。
  4. エネルギーの時間変化の様子を図から確認し、 なぜそのように変化するか考察せよ。

参考

粒子数が256の場合、アルゴンのMD開始直後での温度、全エネルギー、 運動エネルギー、ポテンシャルエネルギーの参考値は以下の通り。

温度 0.72197
全エネルギー -1460.3
運動エネルギー 276.15
ポテンシャルエネルギー-1736.4