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

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

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

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

サンプルプログラム[Fortran, C]　←Cのサンプルを追加(12/06/2002)
サンプル入力データ[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での圧力（単位：σ³／ε）

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       密度（単位：数／σ³）
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]
   　
   • N個の粒子の中からランダムに１個( i )を選択
      [Line: 102]
     　
   • 粒子 i の現在の座標について他の粒子との相互作用エネルギー、ビリアルを計算
      [Line: 103 - 108]
     　
   • 粒子 i をその近傍の空間に均一の確率で移動させる
      [Line: 110 - 115]
     　
   • 粒子 i の新しい座標について他の粒子との相互作用エネルギー、ビリアルを計算
      [Line: 117]
     　
   • メトロポリスの確率にしたがって
     
     新しい座標のエネルギーのほうが小さければ確率１で新しい座標を採用する
      [Line: 123 - 129]
     
     新しい座標のエネルギーのほうが大きければ、新しい一様乱数を振って確率exp(-bDU)で新しい座標を採用する
      [Line: 130 - 137]
     　
   • 現在の瞬間的な物理量を平均操作に反映する
      [Line: 138 - 146]
     　
7. 指定した周期で行う処理を実行
   　
   • 途中のデータの出力
      [Line: 149 - 157]
   • 選択した粒子の最大移動距離の見直し
      [Line: 159 - 161]
     　
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. 上の説明を参考にしながらサンプルプログラムをよく読み、内容を理解せよ。 サンプルでは乱数を発生させる関数にFortranやＣに付属（＝いい加減）のものを用いています。そこが気になる人はNumerical Reciepesなどを参考にして改善を試みてください。実際モンテカルロシミュレーションでは良い乱数を用いることが本質的に重要です。
   　
2. 温度（単位：ε）を４種類[0.5, 1.0, 2.0, 4.0]、密度（単位：数／σ³）を５種類[0.1, 0.2, 0.4, 0.6, 0.8, 1.0]の全組み合わせ（合計２０組）の状態についてモンテカルロシミュレーションを行い、圧力のカノニカル平均値を求めて密度（あるいはその逆数∝体積）−圧力の等温線（４本）をグラフ化せよ。
   　
3. 演習２の等温線の特徴を考察せよ。データが足りなければシミュレーションを追加して適宜補うこと。
   
   考察ポイントの例：
   　　・低温と高温で圧力の密度依存性にどのような違いがあるか
   　　・低温での圧力等温線の挙動は物理的にどのような現象を反映しているか
   　　・なぜ大きな負の圧力が実現されるのか（現実の実験ではどうなるべきか？）
   　　・理想気体やファンデルワールスの状態方程式（＝平均場理論）との比較
   　　・マクスウェルの等面積則（体積−圧力のプロット）と気体液体の共存
   　　・臨界点のおおざっぱな見積り
   　

  《レポート》

1. 演習２，３をメールで以下のアドレスに送ること。
   
   提出先：　ryoichi@r02.mbox.media.kyoto-u.ac.jp
   締切り：　　１２月　１４日（火）
   
   氏名、学籍番号を忘れずに。授業に関する要望・感想があれば遠慮せずに書いてください。
   
   提出者：　
   
   0500-16-3529
   
   0500-15-0173
   0500-15-0389
   0500-15-0404
   0500-15-0638
   0500-15-0656
   0500-15-0807
   0500-15-0852
   0500-15-1279
   0500-15-1528
   0500-15-1537
   0500-15-1804
   0500-15-1921
   0500-15-1958
   0500-15-1976
   0500-15-2033
   0500-15-2991
   0500-15-3077
   0500-15-3246
   0500-15-3596
   0500-15-3676
   0500-15-3845
   0500-15-3999
   0500-15-4252
   0500-15-4271
   0500-15-4539
   0500-15-4557
   0500-15-5142
   0500-15-6211
   0500-15-6220
   0500-15-6337
   0500-15-6991
   0500-15-7236
   0500-15-7281
   0500-15-7451
   0500-15-7782
   0500-15-8091
   0500-15-8233
   0500-15-8289
   0500-15-8313
   0500-15-8725
   0500-15-8805
   0500-15-8823
   0500-15-9179
   0500-15-9464
   0500-15-9473
   0500-15-9571
   0500-15-9581
   0500-15-9689
   0500-15-9894
   
   0500-14-0721
   0500-14-1030
   0500-14-1746
   0500-14-2350
   0500-14-2476
   0500-14-2529
   0500-14-2645
   0500-14-2734
   0500-14-2841
   
   0500-13-0466
   0500-13-0420
   0500-13-2326
   
   0500-11-2154
   0500-11-3071
   
   6100-15-2893
   
   6100-14-1218
   
   [2004.12.16 23:35]
   
   ＊受け取りの確認の意味で名前を載せています。掲載を希望しない人があれば対応しますので、レポート毎にメールに希望を記載してください。