kumpeishiraishi/mc_dynamics_tips.md

## mc_dynamics_tips.md

      
    Raw
  

              mc_dynamics_tips.md
            
          
    モンテカルロ法の軌道の解析について私が知っている二、三の事柄

モンテカルロ法は、粒子系のシミュレーションをする際、とてもよく使われる計算手法です。
静的な性質（ポテンシャルなど）を調べる際の配置サンプリングに使われるのはもちろん、モンテカルロ法のダイナミクスそれ自体が計算されることもよくあります。
ここで言うダイナミクスとは、平均二乗変位や様々な時間相関関数のことです。
分子動力学法（MD）のダイナミクスの計算方法はよく解説を見かける一方、モンテカルロ法（MC）の解説はさほど多くないように思います。
ここでは、粒子系のMCシミュレーションをどう実行し、得られた軌道をどのように解析するのか、簡単にまとめます。
まずは、本稿で考える粒子モデルであるKob-Andersen（KA）モデルを導入します。
KAモデルとは、過冷却液体やガラスの最も基本的なモデルです。
粒子はLennard-Jonesポテンシャルで相互作用します。
$$
V(r_{ij}) = 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right].
$$
系内の粒子はA粒子とB粒子の2成分であり、上記の相互作用中のパラメータは以下のように取ります：
$\epsilon_{AA} = \epsilon$ 、 $\epsilon_{AB} = 1.5\epsilon$ 、
$\epsilon_{BB} = 0.5\epsilon$ 、 $\sigma_{AA} = \sigma$ 、
$\sigma_{AB} = 0.8\sigma$ 、 $\sigma_{BB} = 0.88\sigma$ 。
ポテンシャルはカットオフ距離 $r_c = 2.5\sigma_{ij}$ で連続になるようにシフトします。
$N$ 粒子が一辺の大きさ $L$ の箱の中に入っており、各方向で周期境界条件を採用します。
$L$ は数密度 $\rho = N/L^3 = 1.204$ より決まります。
長さ、エネルギーの単位はそれぞれ $\sigma$ 、 $\epsilon$ とし、時間の単位は1 MCステップ、つまり $N$ 回の試行です。
さて、ここで導入したKAモデルのダイナミクスとして、平均二乗変位（mean-squared displacements, MSD）を計算しましょう。
$$
\Delta^2(t) = \frac{1}{N} \sum_{i=1}^N | \vec{r}_i(t) - \vec{r}_i(0) |^2
$$
以下、

モンテカルロ法による粒子系の時間発展
シミュレーションで得られた軌道の解析

に分けて説明します。
適宜、C++によるサンプルコードを挿入しますが、その中で使用する疑似乱数生成器、一様実数分布生成器は
std::mt19937 mt(123);
std::uniform_real_distribution<double> udist(0.0, 1.0);

とします。 粒子位置は配列
double conf[N][deg];

に格納されているとします（ N は粒子数、 deg は空間次元）。
モンテカルロ法による粒子系の時間発展

MCの時間発展は主に3ステップに分けられます。

粒子をランダムに選択する
選択した粒子に小さい変位を与えて、次の状態を提案する
エネルギー差を評価して、メトロポリス条件で提案を採択するか決定する

順番に見ていきましょう。
粒子 $i$ の選択は、以下のように行えます。
const int i = N * udist(mt);

次に、選択した粒子の位置に対して、小さな変位ベクトルを与えます。
$$
\vec{r}_i = \vec{r}_i + \delta\vec{r}.
$$
変位ベクトルは、一辺の大きさ $\delta$ の箱の中のランダムなベクトルで、以下のように得られます。
const double dx = delta * (udist(mt) - 0.5);
const double dy = delta * (udist(mt) - 0.5);
const double dz = delta * (udist(mt) - 0.5);

$\delta$ の大きさは通常、時間相関関数の緩和が最も速い値を選びます。
KAモデルの場合は $\delta = 0.15$ が用いられます。
こうして発生させた変位ベクトルを、上で選択した粒子 $i$ の位置に与えましょう。
conf[i][X] += dx;
conf[i][Y] += dy;
conf[i][Z] += dz;

この変位ベクトルの加算の前後で、粒子 $i$ のポテンシャル差 $\Delta E_i$ を評価し、メトロポリス条件で採択を判断します。 コーディング的には
udist(mt) < std::exp(-beta*(Enew - Eold));

のとき採択、という if 文を書けば良いでしょう（ $\beta = 1/T$ ）。
粒子 $i$ のポテンシャル（ Eold と Enew ）は試行毎に計算します。
その際、MDでよく行うように、ベルレの隣接リストを使うと計算時間を短縮できます。
提案を採択する場合、周期境界条件には注意が必要です。
時々刻々の粒子位置はminimum-image conventionで配列 conf に保存していきますが、周期境界条件の何番目のセルに粒子が位置しているのかも、同時に記録しなければなりません。
その用途に、配列
int box[N][deg];

を用いることにします。 上記の提案が採択された場合は、
double xbox = std::floor(conf[i][X] / Lbox + 0.5);
double ybox = std::floor(conf[i][Y] / Lbox + 0.5);
double zbox = std::floor(conf[i][Z] / Lbox + 0.5);
box[i][X] += xbox;
box[i][Y] += ybox;
box[i][Z] += zbox;

として、周期境界条件のセル番号を記録します（ Lbox はシミュレーションの箱の大きさ）。
粒子の位置はminimum-image conventionに戻しておきましょう。
conf[i][X] -= Lbox * xbox;
conf[i][Y] -= Lbox * ybox;
conf[i][Z] -= Lbox * zbox;

これを $N$ 回繰り返し、1 MCステップとします。
そして何ステップも繰り返してシミュレーションを長時間行いますが、後の解析のために、対数間隔の時間間隔で conf と box を保存していきます。
対数間隔で取る理由は、今使っているKAモデルの緩和が低温において極めて遅くなるからです（過冷却液体のモデルだから）。
シミュレーションで得られた軌道の解析方法

さて、MCシミュレーションが終わりました。
保存していた軌道のデータを使って、MSDの解析を行いましょう。
MSDを計算するため、時刻 $0$ の配置 conf_0 と時刻 $t$ の配置 conf_t を読み込みましょう。 表式の通り、MSDは
double MSD() {
    double ans = 0.0;
    for (int i=0; i<N_A; i++) {
        double dx = conf_t[i][X] - conf_0[i][X];
        double dy = conf_t[i][Y] - conf_0[i][Y];
        double dz = conf_t[i][Z] - conf_0[i][Z];
        const double rij2 = dx*dx + dy*dy + dz*dz;
        ans += rij2;
    }
    return ans/N_A;
}

で計算できます。但し、KAモデルは2成分系なので、このコードでは $N_A$ 個のA粒子のみについてMSDを計算しています。
時刻 $t$ の配置を読み込む際には注意が必要です。
MCでは、 時間発展とともに系の重心位置が動きます 。
MDでは通常、シミュレーションの初期速度を与える際に重心速度を除去しますから、時間発展しても系の重心は動きません。
MCでは重心がドリフトするため、これをシミュレーション後の解析の段階で補正する必要があります。
さらに、時刻 $t$ での配置の重心位置は、 周期境界条件を解いた粒子位置で 計算しなければいけません。
言葉で説明するよりコードで書いた方が分かりやすいと思うので、周期境界条件を解いた粒子位置は
for (int i=0; i<N; i++) {
    conf_t[i][X] += Lbox * box_t[i][X];
    conf_t[i][Y] += Lbox * box_t[i][Y];
    conf_t[i][Z] += Lbox * box_t[i][Z];
}

で得られます。この処理をした後で系の重心位置を計算し、重心位置が原点に来るように全体を平行移動します。
double posX_g = 0.0, posY_g = 0.0, posZ_g = 0.0;
for (int i=0; i<N; i++) {
    posX_g += conf_t[i][X];
    posY_g += conf_t[i][Y];
    posZ_g += conf_t[i][Z];
}
posX_g /= N;
posY_g /= N;
posZ_g /= N;
for (int i=0; i<N; i++) {
    conf_t[i][X] -= posX_g;
    conf_t[i][Y] -= posY_g;
    conf_t[i][Z] -= posZ_g;
}

さて、以上の前処理を施してからMSDの計算をしましょう。結果をプロットすると以下のようになります。過冷却液体に特徴的な遅いダイナミクス（中間領域でのMSDのプラトー）が、低温できちんと観察されました。

ちなみに、よくある間違い（筆者が犯した間違い）として、

重心位置の補正を忘れる
重心計算の際に周期境界条件を考慮するのを忘れる

といったものがあります。これらのミスをすると、1 MCステップ進めただけでMSDが大きな値を持つ、という挙動になってしまいます。
時間相関関数である自己中間散乱関数 $F_s(q,t)$ で見ると、ある程度の長時間が経過した時点で相関が急激に減少する、という振る舞いになります。

  
## msd.png

      
    Raw
  

              msd.png