分子動力学（Molecular Dynamics）シミュレーション

Lennard-Jonesモデル流体（粒子数 N＝１０８個）

粒子間ポテンシャル関数：　f(r) = 4e[(s/r)12-(s/r)6]
質量：  m

１．Verocity Verlet 法による運動方程式の差分化（これもまたVerlet法の別表現、講 義ノート参照）

Newtonの運動方程式：

  1. ri(t+dt)  = ri(t) + vi(t)dt + [fi(t)/m]dt2/2
  2. vi(t+dt)  = vi(t) + [{fi(t+dt)+fi(t)}/m]dt/2

Nose-Hooverの運動方程式：　[参考文献： "Understanding Molecular Simulation" by Frenkel and Smit (Academic Press)*]

  1. ri(t+dt)  = ri(t) + vi(t)dt + [fi(t)/m-z(t)vi(t)]dt2/2
  2. vi(t+dt)  = vi(t) + [fi(t+dt)/m-z(t+dt)vi(t+dt)+fi(t)-z(t)vi(t)]dt/2
  3. ln s(t+dt) = ln s(t) + z(t)dt + [Smvi(t)2 - gT]dt2/2Q
  4. z(t+dt) = z(t) + [{Smvi(t+dt)2 - gT} + {Smvi(t)2 - gT}]dt/2Q

    ２、４式において右辺にz(t+dt)vi(t+dt)、Smvi(t+dt)2の 項があるために陽解法では解けない。さらに線形で はないので、行列演算やGauss-Seidelの反復法も使えない。→非線形でも使えるNewton- Raphson法を使えばよい。[参考文献： "数値計算" 赤坂隆（コロナ社）*] 　どのみちfi(t)の計算が律則となるため、２、４式の反復にかかる時間は全体の速度に大して影響しない ことに注意する。

    （ここに挙げた参考文献は２冊とも大変役に立つよい本です）

２．ミクロカノニカル（NVE一定）分子動力学シミュレーション

まず最も簡単な場合として分子動力学シミュレーションでミクロカノニカルアンサンブルを構成することを考えよう。この場合はN個の粒子系についてNewton の運動方程式をコンピューター上で数値的に解けばよい。ここではVelocity Verlet法を用いたプログラムを考える。

サンプルプログラム[ fortran, C ] ←C code added (11/07/2002)
サンプル入力データ[ nve.dat ]

> gfortran -O -o mdnve mdnve.f
or
> gcc -O -o mdnve mdnve.c -lm
> ./mdnve < nve.dat

MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS
ENTER NUMBER OF STEPS
ENTER NUMBER OF STEPS BETWEEN OUTPUT LINES
ENTER THE FOLLOWING IN LENNARD-JONES UNITS
ENTER DENSITY
ENTER POTENTIAL CUTOFF DISTANCE
ENTER THE CORE RADIUS TO AVOID OVERLAP
ENTER TIME STEP
NUMBER OF ATOMS = 108
NUMBER OF STEPS = 15000
OUTPUT FREQUENCY = 100
POTENTIAL CUTOFF = 2.5650
DENSITY = 0.8000
TIME STEP = 0.005000
INITIAL V/N = -2.6597
INITIAL W/N = 45.1070
INITIAL P = 12.0285


**** START OF DYNAMICS ****


TIMESTEP ..ENERGY.. CUTENERGY. ..KINETIC. ..POTENT.. .PRESSURE. ..TEMPER..

100 -2.6670 -1.8815 2.4706 -5.1376 3.2947 1.6470
200 -2.6659 -1.8808 2.4179 -5.0838 3.4000 1.6120
...
...
14900 -2.6711 -1.8805 2.2953 -4.9664 3.6735 1.5302
15000 -2.6678 -1.8831 2.4406 -5.1084 3.6081 1.6270

**** END OF DYNAMICS ****


AVERAGES -2.66762 -1.88130 2.38414 -5.05176 3.63882 1.58943
FLUCTS 0.00398 0.01108 0.11368 0.11406 0.45524 0.07578
　

入力ファイル[nve.dat]の説 明

15000    トータルMDステップ数
100      中間データの出力間隔
0.8      密度（単位：数／σ3）
3.       相互作用の打ち切り距離（単位：σ）
0.9      初期座標作成のための粒子のコア直径（単位：σ）
0.005    数値積分のための時間刻み（単位：(mσ2/e)1/2）
　

３．カノニカル（NVT一定）分子動力学シミュレーション

次に分子動力学シミュレーションでカノニカル（温度一定）アンサンブルを構成することを考えよう。講義で説明したように、この場合はN個の粒子＋粒子速度をス ケールするパラメータ(s)の拡張系についてNose-Hooverの運動方程式を解けばよい。上のNVE同様に、ここではVelocity Verlet法を用いたプログラムを考える。

サンプルプログラム[ fortran, C ] ←C code added (11/07/2002)
サンプル入力データ[ nvｔ.dat ]

> gfortran -O -o mdnvt mdnvt.f
or
> gcc -O -o mdnvt mdnvt.c -lm
> ./mdnvt < nvt.dat

MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS
ENTER NUMBER OF STEPS
ENTER NUMBER OF CYCLES FOR EQUILIBRATION
ENTER NUMBER OF STEPS BETWEEN OUTPUT LINES
ENTER STEPS FOR CALCULATING TIME CORRELATIONS
ENTER THE FOLLOWING IN LENNARD-JONES UNITS
ENTER DENSITY
ENTER TEMPERATURE
ENTER Q
ENTER POTENTIAL CUTOFF DISTANCE
ENTER THE CORE RADIUS TO AVOID OVERLAP
ENTER TIME STEP
NUMBER OF ATOMS = 108
NUMBER OF STEPS = 15000
NUMBER OF EQ STEPS 5000
OUTPUT FREQUENCY = 100
POTENTIAL CUTOFF = 2.5650
DENSITY = 0.8000
TIME STEP = 0.005000
INITIAL V/N = -2.6597
INITIAL W/N = 45.1070
INITIAL P = 12.0285


**** START OF DYNAMICS ****


TIMESTEP ..H_NOSE.. ..H_0..... ..KINETIC. ..POTENT.. .PRESSURE. ..TEMPER..

100 -1.8828 -4.4094 1.3342 -5.7436 -0.0242 0.8895
200 -1.8830 -4.5447 1.2178 -5.7625 -0.1501 0.8119
...
...
14900 -1.8839 -3.9266 1.5484 -5.4750 1.3874 1.0323
15000 -1.8840 -3.8765 1.7704 -5.6469 0.7188 1.1803

**** END OF DYNAMICS ****


AVERAGES -1.88368 -4.03045 1.49995 -5.53041 1.00765 0.99996
FLUCTS 0.01187 0.16732 0.13027 0.09425 0.46261 0.08685
Velocity auto correlation function and mean square dis.:
Number of samples : 10000
Number of initial times : 1000
Timestep between samples: 0.00500
Timestep between t=0 : 0.05000
Diffusion coef. : 0.06603
　

入力ファイル[nvt.dat]の説 明

15000    トータルMDステップ数
5000     平衡化のためのステップ数（始めのこの部分はデータの平均操作に含まない）
100      中間データの出力間隔ステップ数
1        速度自己相関関数を計算する際の間隔ステップ数
0.8      密度（単位：数／σ3）
1.0      温度（単位：ε）
4.0      パラメータsの有効質量（単位：mσ2）
3.       相互作用の打ち切り距離（単位：σ）
0.85     初期座標作成のための粒子のコア直径（単位：σ）
0.005    数値積分のための時間刻み（単位：(mσ2/e)1/2）


出力ファイル[fort.10]

0. 3.0006547 0. 0.00500109093
0.00499999989 2.98975682 7.51437256E-05 0.00998401921
0.00999999978 2.95835733 0.000300056068 0.0149146141
0.0149999997 2.90666604 0.000673147908 0.0197590571
...
...

t        VAF(t)       MSD(t)      ∫0t [VAF(t')/3] dt'
時間  速度自己相関関数  平均２乗変位  VAC(t)の積分値/3


> gnuplot
...
...
gnuplot> set data style lines
gnuplot> plot 'fort.10' u 1:2
gnuplot> plot 'fort.10' u 1:3
　

《演習》

  1. NVE一定の分子動力学シミュレーションのサンプルプログラムをよく読み、内容を理解せよ。CやC＋＋言語に慣れている人はFortranのサンプルプログラムを見て書き 換えて見るとよい（うまく書けた人はCのサンプルにしますのでソースを送ってください）。
    　
  2. Nose-Hoover運動方程式によるNVT一定の分子動力学シミュレーションのサンプルプログラムをよく読み、内容を理解せよ。CやC＋＋言語に慣れている人は Fortranのサンプルプログラムを見て書き換えて見るとよい（うまく書けた人はCのサンプルにしますのでソースを送ってください）。
    　
  3. 前章のモンテカルロシミュレーションの演習と同様に、温度（単位：ε）を４種類[0.5, 1.0, 2.0, 4.0]、密度（単位：数／σ3）を５種類[0.1, 0.2, 0.4, 0.6, 0.8, 1.0]の全組み合わせ（合計２０組）の状態についてNVT一定の条件でシミュレーションを行い、圧力のカノニカル平均値を求めて密度（あるいはその逆数∝体積）−圧力の 等温線（４本）をグラフ化せよ。
    　
  4. モンテカルロシミュレーションと違い、分子動力学シミュレーションでは輸送係数も計算できる。サンプルプログラムでは久保公式にしたがって速度自己相関関数の時間積分から 粒子の自己拡散係数を求めている。前問で行ったシミュレーションそれぞれの温度について、拡散係数を密度に対してプロットし、密度依存性、温度依 存性を議論せよ。
    　
  5. 出力ファイルfort.10には速度自己相関関数VAC(t)、平均２乗変位MSD(t)などのデータが書き込まれている。平均２乗変位を両対数のグラフにプロットし、 Einstein則（MSD(t)=6Dt）が成立していることを確認せよ。またMSD(t)の短時間と長時間における振る舞い（時間依存性） を調べ、その物理的な意味を考察せよ。
    > gnuplot
    ...
    ...
    gnuplot> set data style lines
    gnuplot> set logscale xy
    gnuplot> plot 'fort.10' u 1:3