水（TIP3P）の古典 MD プログラム part 2

これは 水（TIP3P）の古典 MD プログラム part 1 の続きです。
引き続きソースコードの解説をしていきます（自明な部分は適当に端折っています）。
part 1 では、クラスの定義と main 関数の説明を行いました。
part 2 以降は、各種関数の解説を行います。

SHAKE 関数

bool shake(vector<Atom>& atomVector, vector<int>& shake_list)
{
        static const double eps = 1e-6;
        static const double eps2 = eps * eps;
        static const double rOH = 0.9572;
        static const double aHOH = 104.52 / 180. * acos(-1.0);
        static const double dOH = rOH * rOH;
        static const double rHH = rOH * sin(aHOH/2.) * 2.;
        static const double dHH = rHH * rHH;

SHAKE 関数の冒頭です。各種定数を定義しています。
定数 eps は SHAKE エラーの閾値です。 1e-6 はちょっと緩めかもしれませんね。
dOH と dHH がそれぞれ OH と HH の固定距離（の 2 乗）です。

for (int i = 0; i < shake_list.size() / 2; i++)
{
        Atom& at1 = atomVector[shake_list[2 * i]];
        Atom& at2 = atomVector[shake_list[2 * i + 1]];
        double gamma;
        if (at1.atomname == "O" && at2.atomname == "H")
        {
                gamma = (dOH - norm(at1.rnew - at2.rnew)) /
                        (2.*(1./at1.mass+1./at2.mass)*((at1.position - at2.position) * (at1.rnew - at2.rnew)));
        }
        else
        {
                gamma = (dHH - norm(at1.rnew - at2.rnew)) /
                        (2.*(1./at1.mass+1./at2.mass)*((at1.position - at2.position) * (at1.rnew - at2.rnew)));

        }
        at1.rnew = at1.rnew + gamma * (at1.position - at2.position) / at1.mass;
        at2.rnew = at2.rnew + gamma * (at2.position - at1.position) / at2.mass;
}

SHAKE のメインループです。あらかじめ SHAKE すべき結合のペアを探しておいて、 shake_list に保存してあります。
$N_c$ 個の拘束条件を考えます。¹

\sigma_k (\vec{r}_1,\vec{r}_2,...,\vec{r}_N)=0\qquad k=1,2,\cdots,N_c

ここで $\vec{r}_i$ は $i$ 番目の粒子の座標です。
拘束が存在する条件での運動方程式は、

m_i \ddot{\vec{r}}_i =\vec{f}_i+\sum_{k=1}^{N_c}\lambda_k\frac{\partial\sigma_k}{\partial\vec{r}_i}\qquad i=1,2,\cdots,N

ここで $\lambda_k$ は未定乗数です。この定数を、シミュレーションの最中に繰り返し計算で決定する方法が SHAKE です。²
時間に関するテイラー展開の 2 次までとると時間発展は

\vec{r}_i(\Delta t)=\vec{r}_i(0)+\dot{\vec{r}}_i(0)\Delta t + \frac{\Delta t^2}{2m_i}\left[ \vec{f}_i(0)+\sum_{k=1}^{N_c}\lambda_k\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}\right]

ただし $\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}$ は、 偏微分を t = 0 で評価することを意味します。
ここで以下のように分離すると：¹

{\vec{r}_i}^\prime=\vec{r}_i(0)+\dot{\vec{r}}_i(0)\Delta t+\frac{\Delta t^2}{2m_i}\vec{f}_i(0)\\
\vec{r}_i(\Delta t)={\vec{r}_i}^\prime+\frac{\Delta t^2}{2m_i}\sum_{k=1}^{N_c}\lambda_k\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}

$ {\vec{r}_i}^\prime $ が拘束のないときの通常の 1 step で、拘束による補正を加えたものが $\vec{r}(\Delta t)$ だと考えることができます。

$\vec{r}_i(\Delta t)$ にも拘束条件を課すと、

\sigma_l\left( \vec{r}_1(\Delta t), \cdots,\vec{r}_N(\Delta t)\right)=\\
\sigma_l\left( {\vec{r}_1}^\prime+\frac{\Delta t^2}{2m_1}\sum_{k=1}^{N_c}\lambda_k\frac{\partial \sigma_k(0)}{\partial\vec{r}_1},
\cdots,
{\vec{r}_N}^\prime+\frac{\Delta t^2}{2m_N}\sum_{k=1}^{N_c}\lambda_k\frac{\partial \sigma_k(0)}{\partial\vec{r}_N}
 \right)=0,\qquad l=1,2,\cdots,N_c

この $N_c$ 個の方程式を直接解いて解（未定乗数）を得るのは難しいので、 SHAKE の出番というわけです。

（SHAKE）
まず、未定乗数の initial guess {$\lambda_k^{(1)} $} を与えます（本ソースコードでは 0 です）。
このとき座標を以下のようにアップデートします：

{\vec{r}_i}^{(1)}={\vec{r}_i}^{\prime}+\frac{\Delta t^2}{2m_i}\sum_{k=1}^{N_c}\lambda_k^{(1)}\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}

未定乗数は

\lambda_k=\lambda_k^{(1)}+\delta \lambda_k^{(1)}

と書けます。あとは $\delta \lambda_k^{(1)}$ が求まれば解決です。その近似的な表式の導出¹は割愛しますが結果は

\delta \lambda_k^{(1)}=-\frac{\sigma_k({\vec{r}_1}^{(1)},\cdots,{\vec{r}_N}^{(1)})}
{\sum_{i=1}^N\frac{\Delta t^2}{2m_i}\frac{\partial\sigma_k({\vec{r}_1}^{(1)},\cdots,{\vec{r}_N}^{(1)})}{\partial \vec{r}_i}\frac{\partial\sigma_k(0)}{\partial\vec{r}_i}}

となります。すると

\vec{r}_i(\Delta t)={\vec{r}_i}^{(1)}+\frac{\Delta t^2}{2m_i}\sum_{k=1}^{N_c}\delta\lambda_k^{(1)}\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}

によって、アップデートされた座標を得ることができます。 $ \delta\lambda_k^{(1)} $ は近似式なので実際には

{\vec{r}_i}^{(n+1)}={\vec{r}_i}^{(n)}+\frac{\Delta t^2}{2m_i}\sum_{k=1}^{N_c}\delta\lambda_k^{(n)}\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}

として、全ての拘束条件が満たされるまで SHAKE のループが続きます。
上式では、 $\sum_{k=1}^{N_c}$ を計算するために、まず、未定乗数 $\delta\lambda_k$ をあらかじめ全て決定しておきます。座標の更新はそのあとに初めて行えます。
これに対して上のソースコードでは、ある拘束条件 $\sigma_k$ に付随する未定乗数 $\delta\lambda_k$ を決定するたびに、その拘束に関与する原子の座標を

{\vec{r}_i}^{\text{new}}={\vec{r}_i}^{\text{old}}+\frac{\Delta t^2}{2m_i}\delta\lambda_k\frac{\partial \sigma_k(0)}{\partial\vec{r}_i}

と更新しています（ for ループの中身）。このループを一巡したところで、拘束条件が満たされているかチェックし（下記のソースコード）、満たされていなければ次のループに入ります。
どちらの方が効率が良いかはわかりません。

ここまでは拘束の具体的な形を特に定めていませんでした。今回のように結合長を固定したい場合、 $k$ 番目の拘束の式は

\sigma_k(\vec{r}_a,\vec{r}_b)=(\vec{r}_a-\vec{r}_b)^2-d_{ab}^2=0

となるでしょう。ここで a と b は $k$ 番目の結合長の拘束に関わる粒子のインデックスです。すると $\delta \lambda_k^{(1)}$ の具体的な表式は

\delta \lambda_k^{(1)}=-\frac{({\vec{r}_a}^{(1)}-{\vec{r}_b}^{(1)})^2-d_{ab}^2}
{2\Delta t^2 \left( \frac{1}{m_a}+\frac{1}{m_b} \right)\left( {\vec{r}_a}^{(1)}-{\vec{r}_b}^{(1)}  \right)\cdot \left( \vec{r}_a(0)-\vec{r}_b(0)  \right)}

となります。ただし、教科書によっては未定乗数の定義が異なるため、上の $\delta \lambda_k^{(1)}$ に相当する表式が微妙に異なることがあるので注意してください。上のソースコードは、岡崎・吉井の本³の表式に従っています。

上のソースコードで、未定乗数の値を特に保存していないことに注意してください。これは NVE/NVT しか想定していないためです。圧力計算のように、ビリアルへの拘束力の寄与を explicit に計算する必要があるときは、拘束力を保存しておく必要があります（いわゆる Monte Carlo barostat を使えばビリアルを計算することなく圧力を制御することが可能なようですね）。

for (int i = 0; i < shake_list.size() / 2; i++)
{
        Atom& at1 = atomVector[shake_list[2 * i]];
        Atom& at2 = atomVector[shake_list[2 * i + 1]];
        double r = vabs(at1.rnew - at2.rnew);
        double error;
        if (at1.atomname  == "O" && at2.atomname == "H")
        {
                error = abs(r - rOH);
        }
        else
        {
                error = abs(r - rHH);
        }
        if (error > eps) return false;
}

補正された座標が拘束条件を満たしているかをチェックしています。拘束条件を満たさない場合、 false を返して main 関数に戻ります。

正規分布に従う乱数の生成

double gauss()
{
        const double A1 = 3.949846138, A3 = 0.252408784;
        const double A5 = 0.076542912, A7 = 0.008355968;
        const double A9 = 0.029899776;
        double sum = 0.;
        for (int i = 0; i < 12; i++)
        {
                sum += rand() / static_cast<double>(RAND_MAX);
        }

        double R = (sum - 6.) / 4.;
        double R2 = R * R;
        return ((((A9 * R2 + A7) * R2 + A5) * R2 + A3) * R2 + A1) * R;
}

Allen と Tildesley の本（初版）⁴に載っている、一様分布から正規分布を生成する方法です。
12 個の一様乱数を足し合わせる（ for 文）と、その和は中心極限定理によって正規分布に従うようになるでしょう、というアイデアです。
その下の式は、 for 文によって生成した、 0 から 12 の間の正規分布を、中心 0 、分散 1 に修正する式のようですが、詳細は Knuth の本⁵に書いてあるみたいです（調べたことないです）。
ちなみに初版にはこの方法の他に ボックス＝ミュラー法 というアルゴリズムも紹介されています。この方法は、上記の 12 個（だけ）の乱数を使用する方法よりも安定なようですが、 log, sin, cos, sqrt を使用するので、初版が発行された当時は計算速度の観点から、使用をためらうには十分な理由になったようです。
しかし、 2017 年の世界ではさすがに問題ないと判断されたのか、 2017 年刊行の第二版ではボックス＝ミュラーだけが紹介されており、中心極限定理に依存する方法は削除されてしまいました。
なお、 C++11 では、正規分布を生成する標準ライブラリが利用できます（正規分布を生成する方法は、処理系によって異なるみたいですが、その辺りの品質保証はどうなっているのだろう？）。

以上で part 2 を終わります。次回はいよいよエネルギーと力の計算の解説に入ります。
続き→水（TIP3P）の古典 MD プログラム part 3

1. Tuckerman, M. E. (2010) "Statistical mechanics: Theory and molecular simulation", Oxford University Press. ↩ ↩² ↩³

2. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. (1977) "Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes", J. Comput. Phys., 23, 327-341. ↩

3. 岡崎進; 吉井範行 (2011) "コンピュータ・シミュレーションの基礎", 第 2 版, 化学同人. ↩

4. Allen, M. P.; Tildesley, D. J. (1987) "Computer simulation of liquids", Oxford University Press. ↩

5. Knuth, D. (1973) "The art of computer programming", 2nd edn, Addison-Wesley, Reading, MA. ↩