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

タイトルのとおりです。ソースコード（md_water_ewald.cpp）は
https://github.com/poinco-gogo/md_water
で配布しています。学生時代に作ったものなので色々とツッコミどころ満載ですが、動くことは動きます（重要）。
力場は TIP3P で、 NVE/NVT のみの対応となります。
位置 Verlet 法を採用[^1]し、強制的に SHAKE します。[^2]
温度制御は Langevin 方程式に従っています（[BBK integrator](http://localscf.com/localscf.com/LangevinDynamics.aspx.html)）。[^3]
式を explicit に書くなら

```math
r(t+\Delta{t})=\frac{1}{B}\left\{2r(t)-Ar(t-\Delta{t})+\frac{\Delta t^2}{m}\left[F(t)+R(t)\right]\right\}\\
A=1-\gamma\frac{\Delta t}{2}\\
B=1+\gamma\frac{\Delta t}{2}\\
<R(0)R(t)>=2k_{B}T\gamma m \delta(t)
```
$R(t)$ はガウスランダム過程で、上式の最終行を満たします。 $T$ は熱浴の温度、 $\delta(t)$ はディラックのデルタ関数です。摩擦係数 $\gamma=0$ とすると通常の位置 Verlet の式に戻ることに注意してください。インプットの config ファイルで `gamma 0`とすることで NVE のシミュレーションをすることができます。

以下、コードの解説をしていきたいと思います。
（commit: 8b01725 に基づく）

# 冒頭
```C++
#include "vector3.h"
```
このプログラムでは 3 次元ベクトルの計算に
http://cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/SRC/vector3.h
を使用します（ただし、 abs を vabs に変更して下さい）。今なら [Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page) を使う方があとで行列演算をしたくなった時に便利でしょう。

# クラス定義（Atom）
```C++
class Atom
{       
        public: 
        Vector3 position;
        Vector3 velocity;
        Vector3 fold, fnew, vnew, rold, rnew;
        string  atomname;
        double  charge, mass, R; 
        Atom(double x, double y, double z)
        {       
                position.x = x;
                position.y = y;
                position.z = z;
        }
};
```
原子を記述するクラスです。R は温度制御の時に使用する定数です（下記参照）。

# クラス定義（System）
```C++
class System
{
        public:
        int nstep, natom, kmax;
        double dt, T, gamma, boxsize, cutoff, volume, ewcoeff, lj, es;
        Vector3 origin, a1, a2, a3, g1, g2, g3, b1, b2, b3;

        System()
        {
                dt = 9999;
                T = 9999;
                gamma = 9999;
                nstep = 0;
                natom = 0;
                boxsize = 1;
                origin.x = 9999; origin.y = 9999; origin.z = 9999;
                cutoff = 1.;
                lj = 0;
                es = 0;
        }

        void setup_box()
        {
                volume = boxsize * boxsize * boxsize;
                a1.x = boxsize; a1.y = 0; a1.z = 0;
                a2.x = 0; a2.y = boxsize; a2.z = 0;
                a3.x = 0; a3.y = 0; a3.z = boxsize;
                double rvol2pi = 2. * M_PI / volume;
                g1 = (a2 % a3) * rvol2pi;
                g2 = (a3 % a1) * rvol2pi;
                g3 = (a1 % a2) * rvol2pi;
        }
};
```
系の情報を扱うクラスです。gamma は Langevin 方程式の friction coefficient というやつです。

# main 関数
```C++
if (argc < 4)
{
        cout << "usage: ./a.out trjout psf pdb config\n";
        return 1;
}
```
引数が足りないときに使い方を表示します（return 1 でいいのか？）。

| 引数       |       意味 |
|:----------------:|:-----------------|
| trjout           | トラジェクトリを保存するアウトプットファイル名。トラジェクトリは XYZ 形式で保存します。            |
| psf           |系に存在するすべての分子の結合（bond, angle,...）の情報を記述するファイル（psf 形式）|
| pdb             |初期座標（PDB 形式）|
| config               |MD のパラメータを記述するファイル|
[Github](https://github.com/poinco-gogo/md_water) に、 psf, pdb, config ファイルのサンプルがあります。

```C++
srand( (unsigned int)time(NULL) );
```
現在時刻をもとに乱数の発生系列を定義します。あとで rand() によって乱数を発生させていますが、 C++11 ではメルセンヌ・ツイスターが標準ライブラリに実装されているので、そちらを使う方がよいでしょう（[参考](https://cpprefjp.github.io/reference/random/mt19937.html)）。

```C++
const int nfree = natom * 3 - nwat * 3;
```
系の自由度を定義しています。強制的に SHAKE しているので、水の結合の自由度を全体から引いておきます。

```C++
const double ps2asu = 20.455;
```
本プログラムでは原子の質量を原子質量単位（atomic mass unit）、エネルギーの単位を kcal/mol、長さの単位をオングストロームで表現しています。この場合、時間の単位はピコ秒に ps2asu をかけたものになります。

```C++
const double dt_ps = system.dt * 0.001; // in ps
```
インプットの config ファイルに `dt 2` と書いた場合、それはタイムステップが 2 フェムト秒であることを意味します。なのでいったんピコ秒に変換します。

```C++
const double A = 1. - gamma * dt * 0.5;
const double B = 1. + gamma * dt * 0.5;
const double inB = 1. / B;
```
冒頭の式の $A$ と $B$ を計算しています。

```C++
for (int i = 0; i < natom; i++)
{      
        Atom& at = atomVector[i];
        at.R = sqrt(2. * kb * T * gamma * at.mass / dt);
}
```
各原子ごとに上式の $R$ の平均的な大きさを計算します。これに正規分布 $N(0, 1)$ に従う乱数を乗じたものが、ランダム力になります。

```C++
const double rOH  = 0.9572;
const double aHOH = 104.52 / 180. * acos(-1.0);
```
TIP3P のジオメトリを定義しています。

```C++
// make reciprocal vectors
vector<Vector3> g;
int kmax = system.kmax;
int sqkmax = kmax - 1;
sqkmax *= sqkmax;
double dum = 0;
for (int k = 0; k < kmax; k++)
for (int i = -kmax + 1; i < kmax; i++)
for (int j = -kmax + 1; j < kmax; j++)
{
        if (k * k + j * j + i * i > sqkmax or (i==0 && j==0 && k==0))
                continue;
        Vector3 vtmp(i * system.g1.x, j * system.g2.y, k * system.g3.z);
        if (vabs(vtmp) > dum)
                dum = vabs(vtmp);
        g.push_back(vtmp);
}
```
逆格子ベクトル $\vec{G}=2 \pi \left(\frac{k_{x}}{L_{x}},\frac{k_{y}}{L_{y}},\frac{k_{z}}{L_{z}} \right)$ を定義しています。 $k_{x}$, $k_{y}$ については、 config ファイルで指定した `kmax` の値をもとに -kmax から +kmax までとりますが、 $k_{z}$ については0 から +kmax までしかとりません。この条件の下で結果を 2 倍すれば求めるべき値を得ることができます。[^4]
`kmax` は、静電相互作用エネルギーが収束するまで充分大きくとりましょう。
（エワルド係数 `ewcoeff` との兼ね合いで設定するべき `kmax` の値が変わることに注意してください。）

```C++
double ew_self = 0;
for (int i = 0; i < atomVector.size(); i++)
{
        ew_self += atomVector[i].charge * atomVector[i].charge;
}
ew_self *= -system.ewcoeff / sqrt(M_PI) * 332.0636;
```
エワルドの自己相互作用エネルギー項です。332.0636 は、$e^2/4 \pi \epsilon_{0}$ だったような気がしますが、記憶が定かではありません。

```C++
for (int i = 0; i < natom; i++)
{
        Atom& at = atomVector[i];
        at.velocity.x = gauss() * sqrt(kb * T / at.mass);
        at.velocity.y = gauss() * sqrt(kb * T / at.mass);
        at.velocity.z = gauss() * sqrt(kb * T / at.mass);
}
```
初速度をボルツマン分布からサンプリングしています。

```C++
for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.fold = at.fnew;
        at.rold = at.position;
        at.position = at.rold + dt * at.velocity + dtdt_div2 / at.mass * at.fold;
}
```

冒頭の式で、 $r(t+\Delta t)$ を計算するために $r(t)$ と $r(t-\Delta t)$ が必要です。ここではインプットの座標が $r(t-\Delta t)$ だと思うことにして、 $r(t)$ をテイラー展開の 2 次までとることで求めます。

```math
r(t) = r(0)+v(0)\Delta t + \frac{\Delta t^2}{2m}F(0)
```

```C++
for (int istep = 1; istep <= nstep; istep++)
{
```
この中が MD のメインループです。

```C++
calc_frc(atomVector, lj_pair_list, el_pair_list, system, g);
```
まず `position` に保存された座標に基づいて force を計算します。計算された force は `fnew` に保存されます。

```C++
for (int i = 0; i < atomVector.size(); i++)
{
        Vector3 noise( gauss(), gauss(), gauss());
        Atom& at = atomVector[i];
        at.rnew = 2. * at.position - at.rold + gamma * dt_div2 * at.rold + dt * dt / at.mass * (at.fnew + at.R * noise);
        at.rnew = at.rnew * inB;
}
```
冒頭の式に基づいて $r(t+\Delta t)$ を計算します（ `rnew` に保存される）。

```C++
for (int ishake = 0; ishake < 100; ishake++)
{
        if (shake(atomVector, shake_list))
                break;
        if (ishake == 99)
        {
                cerr << "error: shake does not converged.\n";
                return 1;
        }
}
```
SHAKE ループです。

```C++
for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.vnew = 0.5 / dt * (at.rnew - at.rold);
}
```
$v(t+\Delta t)=\frac{r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}$ で、更新された速度ベクトルを計算します。結果は `vnew` に保存されます。

```C++
double K = calc_kin(atomVector);
calc_pot(atomVector, lj_pair_list, el_pair_list, system, g);
system.es += ew_self;
```
運動エネルギー、ポテンシャルエネルギーを計算します。
運動エネルギーは `vnew` に基づいて計算しています。
ポテンシャルエネルギーは `position` に基づいて計算しているみたいです。これは変ですね。 `rnew` に基づく計算であるべきではないだろうか？

```C++
for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.rold = at.position;
        at.position = at.rnew;
        at.velocity = at.vnew;
}
```
次のステップに備えて変数の値を入れ替えています。

以上で MD のループは終わりです。
ちょっと長くなってきたので、続きは別にします。
続き→[水（TIP3P）の古典 MD プログラム part 2](https://qiita.com/poinco_gogo/items/34c925f7b66cbe4bf576)

[^1]: 拘束の導入が簡単なので学習目的にはよいでしょう。しかし速度 Verlet の方が velocity のエラーが少ないので良いかもしれません。有名な MD プログラムは大抵速度 Verlet か leap frog なのではないでしょうか。
[^2]: 時間の刻み幅を大きくしたいからですね。
[^3]: BBK integrator の速度 velocity バージョンは Bussi, G.; Parrinello, M. "Accurate sampling using Langevin dynamics", Phys. Rev. E, **2007**, *75*, 056707. で提案された integrator の、時間の刻み幅が小さい極限で一致する、気がする。
[^4]: Tuckerman, M. E. "Statistical Mechanics: Theory and Molecular Simulation", Oxford Graduate Texts, pp. 661