メトロポリスの方法によるモンテカルロシミュレーション

Lennard-Jonesモデル流体（カノニカルアンサンブル、粒子数＝１０８個）

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

[...]　の中の第一項は近距離での粒子の反発を表し、第二項は粒子間の引力を表す。このLennard-Jonesモデルはアルゴンなど希ガスの性質を比較的よく記述し、与えられた熱力学状態によって気体・液体・固体の３相（あるいはそれらの共存状態）をとることが出来る。

ここではメトロポリスの方法によって、カノニカルアンサンブルに基づいたサンプルを次々に発生し、与えられた熱力学状態[NVT]のもとで系のポテンシャルエネルギーと圧力の平均値を与えることのできるプログラムを考える。

サンプルプログラム[Fortran, C]
サンプル入力データ[in.dat]

> f77 -O -o mcnvt mcnvt.f (or gcc -O -o mcnvt mcnvt.c -lm)
> ./mcnvt < in.dat
...
...
AVERAGES

 <V/N> = -4.885490   与えたNVTでの１粒子当たりのポテンシャルエネルギー（単位：ε）
<P> = 0.009233     与えたNVTでの圧力（単位：σ3／ε）

FLUCTUATIONS

FLUCTUATION IN <V/N> = 0.000000
FLUCTUATION IN <P> = 0.364621

END OF SIMULATION
　

[in.datの説明]

550000    トータルモンテカルロステップ数
50000     平衡化のためのステップ数（始めのこの部分は平均操作に含まない）
10000     中間データの出力間隔
10000     最大変位距離の見直し間隔（新しい座標のアクセプト率が0.5以上なら小さく、以下なら大きくする）
0.7       密度（単位：数／σ3）
1.0       温度（単位：ε）
3.        相互作用の打ち切り距離（単位：σ）
0.85      初期座標作成のための粒子のコア直径（単位：σ）
0.1       最大変位距離の初期値（単位：σ）


[mcnvt.fの説明]

  1. 変数の定義
     [Line: 2 - 14]
    　
  2. 入力パラメータの読み込みと確認
     [Line: 20 - 50]
    　
  3. 粒子の初期座標の設定
     [Line: 51 - 67]
    　
  4. エネルギー、圧力に対する打ち切り補正の見積もり
     [Line: 68 - 77]
    　
  5. 初期座標に対するエネルギー、圧力の計算
     [Line: 79 - 85]
    　
  6. メトロポリス法によるモンテカルロループ
     [Line: 86 - 166]
    　
  7. 指定した周期で行う処理を実行
    　
  8. エネルギーと圧力のカノニカル平均値の計算と出力
     [Line: 175 - 186]
    　
  9. サブルーチンの説明

    SUBROUTINE SUMUP ( RCUT, L, V, W )

    このサブルーチンは打ち切り距離RCUT以内にあるすべての粒子対のポテンシャルエネルギーとビリアルを計算する。

    RCUT (INPUT): 粒子間相互作用の打ち切り距離（<0.5L）
    L (INPUT):    シミュレーションセルの１辺の長さ
    V (OUTPUT):   RCUT以内にある粒子対の全ポテンシャルエネルギー
    W (OUTPUT):   RCUT以内にある粒子対の全ビリアル

    粒子の座標RX, RY, RZはLで規格化してあるので[-0.5:0.5]の範囲にある。
    粒子Iから見てJとの[XYZ]座標の差が0.5を超えるか-0.5以下の場合、周期境界条件によってさらに近い距離に別のセルに属するJの虚像が存在する。

    ...
          DO I = 1, N - 1
             RXI = RX(I)
             RYI = RY(I)
             RZI = RZ(I)
             DO J = I + 1, N
                RXIJ  = RXI - RX(J)     粒子IとJのX座標の差
                RYIJ  = RYI - RY(J)     粒子IとJのY座標の差
                RZIJ  = RZI - RZ(J)     粒子IとJのZ座標の差
    C230** MINIMUM IMAGE THE PAIR SEPARATIONS **    以下で周期境界条件を考慮して
                RXIJ  = RXIJ - ANINT ( RXIJ )    Iの最も近くにいるJとのX座標の差
                RYIJ  = RYIJ - ANINT ( RYIJ )    Iの最も近くにいるJとのY座標の差
                RZIJ  = RZIJ - ANINT ( RZIJ )    Iの最も近くにいるJとのZ座標の差
                RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ    粒子IとJの距離^2
                IF ( RIJSQ .LT. RCUTSQ ) THEN    距離^2がRCUT^2より小さい場合だけ計算
                   SR2 = SIGSQ / RIJSQ
                   SR6 = SR2 * SR2 * SR2
                   VIJ = SR6 * ( SR6 - 1.0 )     IJ間のポテンシャルエネルギーの計算
                   WIJ = SR6 * ( SR6 - 0.5 )     IJ間のビリアルの計算
     240           V   = V + VIJ                 全ポテンシャルエネルギーの積算
                   W   = W + WIJ                 全ビリアルの積算
                ENDIF
             ENDDO
          ENDDO
    ...


    SUBROUTINE ENERGY ( RXI, RYI, RZI, I, RCUT, L, V, W )

    このサブルーチンは座標[RXI, RYI, RZI]にある粒子Iと打ち切り距離RCUT以内にあるその他の粒子のポテンシャルエネルギーとビリアルを計算する。

    RXI, RYI, RZI (INPUT): 移動させる粒子の座標（X, Y, Z）
    I (INPUT):    移動させる粒子の番号
    RCUT (INPUT): 粒子間相互作用の打ち切り距離（<0.5L）
    L (INPUT):    シミュレーションセルの１辺の長さ
    V (OUTPUT):   粒子IとIからRCUT以内にある粒子とのポテンシャルエネルギー
    W (OUTPUT):   粒子IとIからRCUT以内にある粒子とのビリアル
    　

　

《演習》

  1. 上の説明を参考にしながらサンプルプログラムをよく読み、内容を理解せよ。CやC＋＋言語に慣れている人は、上のFortranのサンプルプログラムを見て書き換えて見るとよい（うまく書けた人はCのサンプルにしますのでソースを送ってください）。 サンプルでは乱数を発生させる関数にFortranに付属（＝いい加減）のものを用いています。そこが気になる人はNumerical Reciepesなどを参考にして改善を試みてください。実際モンテカルロシミュレーションでは良い乱数を用いることが本質的に重要です。
    　
  2. 温度（単位：ε）を４種類[0.5, 1.0, 2.0, 4.0]、密度（単位：数／σ3）を５種類[0.1, 0.2, 0.4, 0.6, 0.8, 1.0]の全組み合わせ（合計２０組）の状態についてモンテカルロシミュレーションを行い、圧力のカノニカル平均値を求めて密度（あるいはその逆数∝体積）−圧力の等温線（４本）をグラフ化せよ。
    　
  3. 演習２の等温線の特徴を考察せよ。データが足りなければシミュレーションを追加して適宜補うこと。

    考察ポイントの例：
    　　・低温と高温で圧力の密度依存性にどのような違いがあるか
    　　・低温での圧力等温線の挙動は物理的にどのような現象を反映しているか
    　　・なぜ大きな負の圧力が実現されるのか（現実の実験ではどうなるべきか？）
    　　・理想気体やファンデルワールスの状態方程式（＝平均場理論）との比較
    　　・マクスウェルの等面積則（体積−圧力のプロット）と気体液体の共存
    　　・臨界点のおおざっぱな見積り