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

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

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

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

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

速度ベルレの方法

ここでは、数値積分の方法として速度ベルレの式を使用する場合を考えます。時間刻み幅を $\Delta t$ としたときの速度ベルレの式は

$\begin{equation*} \boldsymbol{r}(t+\Delta t) = \boldsymbol{r}(t) + \Delta t\boldsymbol{v}(t) + \frac{\Delta t^2}{2}\boldsymbol{F}(t) \end{equation*}$ (1)

$\begin{equation*} \boldsymbol{v}(t +\Delta t) = \boldsymbol{v}(t)+\frac{\Delta t}{2}\left\{\boldsymbol{F}(t)+\boldsymbol{F}(t+\Delta t)\right\} \end{equation*}$ (2)

となります。ここでは無次元化により $m=1$ としています。

実際の計算では、 $\boldsymbol{F}(t)$ と $\boldsymbol{F}(t+\Delta t)$ を入れる変数は同時に保持していないため、速度の更新の式(2)を以下の2つに分割します。

$\begin{equation*} \boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle) = \boldsymbol{v}(t)+\frac{\Delta t}{2}\boldsymbol{F}(t) \end{equation*}$ (3)

$\begin{equation*} \boldsymbol{v}(t+\Delta t) = \boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle)+\frac{\Delta t}{2}\boldsymbol{F}(t+\Delta t) \end{equation*}$ (4)

メインループ内の「位置の更新、速度の更新」、「力の計算」を速度ベルレの式に従って書くと

$\boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle) = \boldsymbol{v}(t)+\frac{\Delta t}{2}\boldsymbol{F}(t)$	速度の更新[1], 式(3)
$\boldsymbol{r}(t+\Delta t) = \boldsymbol{r}(t)+\Delta t\boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle)$	位置の更新, 式(1)
$\boldsymbol{r}(t+\Delta t)\;\;\;\rightarrow\;\;\;\boldsymbol{F}(t+\Delta t)$	力の計算
$\boldsymbol{v}(t+\Delta t) = \boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle)+\frac{\Delta t}{2}\boldsymbol{F}(t+\Delta t)$	速度の更新[2], 式(4)

となります。

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

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

C言語

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

enum {
    N = 4
};

void output_file(int np, double *cd, double side)
{
  （略）
}

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

  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 < np*3; i++){
    vl[i] = ((double)rand()/RAND_MAX-.5)*.1;
  }
  total_step = 100;
  for(step = 1; step <= total_step; step++){
    for(i = 0; i < np*3; i++){
      cd[i] += vl[i]*dt;
    }
    output_file(np,cd,side);
    for(i = 0; i < np*3; i += 3){
      printf("%d %d %f %f %f\n",step,i/3,cd[i],cd[i+1],cd[i+2]);
    }
  }
  return 0;
}

Python

import numpy as np

def output_file(npa, side, cd, step):
  （略）

npa = 4
side = 9.0
dt = 3.0
cd = np.zeros(npa*3)
vl = np.zeros(npa*3)

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

np.random.seed(1)
for i in range(npa*3):
  vl[i] = (np.random.random() - .5) *.1

total_step = 100
for step in range(total_step):
  for i in range(npa*3):
    cd[i] += vl[i]*dt

  output_file(npa,side,cd,step)

  for i in range(0,npa*3,3):
    print(f'{step} {i//3} {cd[i]:.6f} {cd[i+1]:.6f} {cd[i+2]:.6f}')

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

プログラムを実行すると、cd0000.cdv 〜 cd0099.cdv のファイルが出力されます。マウスで cd0000.cdv を cdview.exe にドラッグ＆ドロップすると、 cd0000.cdv での粒子が表示さます。

'c' キーを押すと cd0001.cdv, cd0002.cdv, cd0003.cdv, … と表示するファイルが切り替わるので、子が動く様子がわかります。ファイルを逆順に切り替えるときは 'z'キーを使います。また、最初のファイルを表示するときは'a'キー、最後のファイルを表示するときは's'キーを使い、ファイルを連続で切り替えるときは SHIFT + 'c'キーや SHIFT + 'z' キーを使います。

演習

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

MDプログラムの作成

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

C言語

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

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

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

   （周期境界条件処理）

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

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

  }
}

となり、 Python 風に書くと

Python

（変数の初期化）
（位置の初期化）
（速度の初期化）
（温度の補正）

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

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

  （周期境界条件処理）

  for i in (粒子数):
    for j in (j+1,粒子数):
      （力の計算）

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

となります。

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

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

温度の補正

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

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

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

初期力

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

このページではプログラミング言語として C 言語と Python を使用していますが、 Python で作ったプログラムの動作速度は C 言語に比べて遅いので、 MD シミュレーションのプログラムは C 言語のものを使うことを想定しています。

Python のプログラムは C 言語のプログラムを理解することを目的に作られており、なるべく C 言語に近い書き方になっています。そのため、 Python のプログラムでは、効率の悪い書き方になっている箇所もありますので注意してください。

Koishi's Page