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

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

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

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

「位置の更新、速度の更新」、「力の計算」、「温度、圧力などの計算」の部分 をメインループと呼びます。このメインループを実行するとΔ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-.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
> cv109 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ステップに一 回の割合で行うようにして下さい。

レポート課題

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

  1. 各ステップで運動エネルギーとポテンシャルエネルギーを計算し、それら の和からトータルエネルギーの時間変化を求めよ。
  2. トータルエネルギーが時間変化に対して保存することを確認し、その精度 を求めよ。
  3. エネルギーの時間変化の様子を調べ、なぜそのように変化するか考察せよ。

参考

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

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