分子動力学法ステップ・バイ・ステップ その6

はじめに

Step 5で温度制御まで書いた。次は圧力制御、と行きたいが、その前に圧力の測定ルーチンを書いてみよう。

今回のコードはgithub/kaityo256/mdstep/step6に置いておく。

注意：結果が正しいかあまりちゃんと確認していません。

• その1　 $O(N^2)$のルーチンまで
• その2　 ペアリストの構築とBookkeeping法
• その3　 メッシュ探索
• その4 　少しだけ高速化
• その5 　温度制御法を三種類実装してみる
• その6 　圧力測定ルーチンの実装 ← イマココ

圧力の導出

分子動力学法において圧力(応力)の定義はさほど自明ではない。圧力には様々な導出があり、ビリアル定理を使うのが一般的であるが、ここでは個人的な趣味で、ハミルトニアンを内部エネルギーとみなして熱力学関係式から導く方法を採用する。

3次元$N$粒子系を考える。粒子番号と座標軸を混ぜたインデックスを取る。すなわち、系は一般化座標$q_i$と一般化運動量$p_i$で記述され、インデックス$i$は$1$から$3N$までの値を取る。ハミルトニアンは

$$
H(\{p_i\},\{q_i\}) = \sum_i^{3N} \frac{p_i^2}{2m} + U(\{q_i\})
$$

簡単のため質量を全て同じとした。さて、ハミルトニアンを内部エネルギーと同一視する。すると、圧力と内部エネルギーの関係から、

$$
p = - \left .\frac{\partial H}{\partial V} \right|_S
$$

であった。さて、いま系が一辺$L$の立方体であるとする。この系をスケール変数$s$を変化させることで一様に膨脹、収縮させよう。スケールされた体積は$V^* = s^3 L^3$である。このスケール変換で運動量、座標はそれぞれ

\begin{align}
p_i^* &= \frac{p_i}{s} \\
q_i^* &= s q_i 
\end{align}

とスケールされる¹。スケールされた変数で記述されたハミルトニアンを$V^*$で偏微分すると、

\begin{align}
P &= - \left. \frac{\partial H}{\partial V^*} \right|_{V^*=V} \\
&= - \left. \frac{d s}{d V^*}\right|_{s = 1}
\left.  \frac{\partial H }{\partial s} \right|_{s = 1}\\
&= - \frac{1}{3V} \sum_i^{3N} \left(
-  \frac{p_i^2}{m}
+ q_i \frac{\partial U}{\partial q_i}
\right)
\end{align}

ただし、途中で$L^3 = V$を用いた。

$k$をボルツマン定数として、分子動力学法における(運動)温度の定義から

\sum_i^{3N} \frac{p_i^2}{m} = 3 N k T

であるから、最終的に圧力の定義として

PV = N k T - \frac{1}{3} \sum_i^{3N} q_i \frac{\partial U}{\partial q_i}

を得る。右辺第一項が理想気体からの寄与分、第二項が相互作用からの寄与分である。

後の便利のため、圧力のポテンシャル由来の項を

\Phi = - \frac{1}{3V} \sum_i^{3N} q_i \frac{\partial U}{\partial q_i}

と表現しておこう。

Lennard-Jones粒子系の場合

先の公式は一般のポテンシャルに適用可能だが、二体相互作用の場合にはもっと簡単になる。表記が一貫せず申し訳ないが、以下では$i$や$j$のインデックスは粒子番号を指すとする。相互作用が二体の場合には、ポテンシャル関数$U$が、二体ポテンシャルエネルギー$\phi(r_{ij})$の和、

$$
U = \sum_{<ij>} \phi(r_{ij})
$$

でかける。ただし

$$
r_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2}
$$

であり、$<i,j>$は全てのペアについて和を取ることを意味する。ここで、$i$粒子の座標を$\vec{q}_i = (x_i, y_i, z_i)$と表記した²。

ここで、

$$
\frac{\partial U}{\partial x_i} =
\sum_{i \ne j} \frac{(x_i-x_j)}{r_{ij}}\phi'(r_{ij})
$$

等に注意して整理すると、

\Phi = - \frac{1}{3V} \sum_{<ij>} r_{ij} \phi'(r_{ij})

一般の6-12型Lennard-Jonesポテンシャル

\phi(r) = 4 \left(r^{-12} - r^{-6}  \right)

であれば、

r \phi'(r) = -48 r^{-12} + 24 r^{-6}

を全てのペアについて加えたものを$-1/3$倍して体積で割ったものが圧力のポテンシャル貢献分$\Phi$となる。

実装

分子動力学法においては、ボルツマン定数を1とする単位系を取ることが多い。従って、密度$\rho = N/V$を使って、

P = \rho T + \Phi

と計算できる。相互作用ペアリストができていれば、$\Phi$の計算は容易であろう。該当コードはこんな感じになる。

double
Observer::pressure(Variables *vars, std::vector<Pair> &pairs) {
  const double N = static_cast<double>(vars->number_of_atoms());
  const double V = L*L*L;
  const double T = temperature(vars);
  double phi = 0.0;
  const int pp = pairs.size();
  Atom *atoms = vars->atoms.data();
  for (int k = 0; k < pp; k++) {
    const int i = pairs[k].i;
    const int j = pairs[k].j;
    double dx = atoms[j].qx - atoms[i].qx;
    double dy = atoms[j].qy - atoms[i].qy;
    double dz = atoms[j].qz - atoms[i].qz;
    adjust_periodic(dx, dy, dz);
    double r2 = (dx * dx + dy * dy + dz * dz);
    if (r2 > CL2)continue;
    double r6 = r2 * r2 * r2;
    double r12 = r6 * r6;
    phi += 48.0/r12 - 24.0/r6;
  }
  phi = phi / 3.0 / V;
  const double density = N/V;
  return density*T + phi;
}

定義そのままの形で計算しているため、難しいところはないと思う。

計算結果

以下の条件で計算してみる。

• 温度 $T = 1.0$
• 密度 $\rho = 0.5$
• 時間刻み $dt = 0.005$
• 温度制御：Nose-Hoover
• システムサイズ: $50 \times 50 \times 50$
• 粒子数: $N=62500$
• カットオフ: 2.5

ちなみに格子定数2で面心立方格子(FCC)に組むと、ちょうど数密度が0.5になるため、密度0.5は都合がよろしい。また、この密度で温度$1.0$というのも圧力の計算のチェックには都合が良い。

まず、温度の時間発展はこんな感じになる。

だいたい$t=20$くらいで$T=1.0$に落ち着いてるみたいですね。

次に圧力はこんな感じ。

密度0.5、温度1.0、カットオフ2.5のLJ系の圧力はほぼゼロになる。運動エネルギーからの圧力の寄与は$0.5$であるから、運動エネルギーによる圧力と、相互作用からくる負圧がだいたいキャンセルしていることになる。

ちなみに$20<t$の領域で圧力を平均すると、$P \sim 0.02$程度になる。データを拡大してみるとわかるが、こういった短距離相互作用系の圧力ゆらぎはかなり大きく、意味のある値を得るには長時間平均する必要がある。

まとめ

分子動力学法における圧力の定義と測定ルーチンの実装を行った。僕にとっては分子動力学法における圧力は非自明な量で、ちゃんとわかってから教科書読み直すとちゃんと書いてあるのだが、わからない状態で読んでも実装できなかった。

ちなみに計算結果があってるかどうかちゃんと確認していない。NISTがLennard-Jones粒子系の熱力学量などを公開しているのだが、ここで使っている単純カットオフとは異なるカットオフ法(Linear-force Shift)を使っていたり、カットオフのlong range correctionsを施していたりと、直接比較できる数字が無いので、比較を諦めた³。気になる人は同じ条件のカットオフ形状を指定してLAMMPSとかで計算すれば良いと思う⁴。

1. ここで$s$によるスケールの仕方が$p$と$q$で反対なのは、一般化速度(ラグランジアン)から一般化運動量(ハミルトニアン)への変換による。$q \rightarrow s q$という変換で$\dot{q} \rightarrow s \dot{q}$となるが、$p \rightarrow p/s$となる(接空間での拡大は、余接空間での縮小になる)。 ↩

2. 座標として$q$使ったり$(x,y,z)$使ったりと揺らいで申し訳ないが、もう記号考えるの面倒になったので許して。 ↩

3. 一応自分が昔書いたコードと結果を比較したが、それでは思想バグは取れないので・・・ ↩

4. っていうか普通にMDやるならLAMMPSなりNAMDなり使った方が早いと思う。 ↩