水（TIP3P）の古典 MD プログラム 番外編

[Github](https://github.com/poinco-gogo/md_water) にアップロードした水（TIP3P）の古典 MD プログラムですが、今までに 3 回に分けて解説を行ってきました。

- [水（TIP3P）の古典 MD プログラム part 1](https://qiita.com/poinco_gogo/items/f7e6db1c85f51c739bfc)
- [水（TIP3P）の古典 MD プログラム part 2](https://qiita.com/poinco_gogo/items/34c925f7b66cbe4bf576)
- [水（TIP3P）の古典 MD プログラム part 3](https://qiita.com/poinco_gogo/items/0e18cd1ffeda4dcefb65)

このプログラム（コミット： 807c6db）をインテルコンパイラ（バージョン 18.0.1）で `-O3` オプションつきでコンパイルし、

```
$ icpc md_water_ewald.cpp -O3 -o md_water_ewald.x
```

MD のオプションを記述するファイル config の中で `numstep 10` とし、 10 ステップの MD を実行します。

```
$ time ./md_water_ewald.x out.xyz water.psf input.pdb config
```

上記の `time` はプログラムの実行時間を計測するコマンドですね。
プログラムの終了時に概ね以下のように表示されます。

```
real    0m21.418s
user    0m21.414s
sys     0m0.004s
```
今回は 10 ステップに 21 秒かかっています。
ちなみにプログラムを実行したマシンの CPU は Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10GHz です。

同等のシミュレーションを、代表的な生体高分子分子動力学シミュレーションプログラムである [NAMD](http://www.ks.uiuc.edu/Research/namd/) 2.13 （シングルコア）で実行すると・・・

```
real    0m0.160s
user    0m0.096s
sys     0m0.009s
```

となり、自作プログラムの **134 倍速い**ことがわかります。
そこでもう少し使い勝手をよくするために、自作プログラムを整頓しました。プログラムは [Github](https://github.com/poinco-gogo/md_alad) にアップロードしてあります。

以下の解説はコミット: ~~`7b19cd5`~~ ~~`6f05eeb`~~ `780a659` に基づいています。
（コミット `7b19cd5` では、最新の helPME ライブラリに対応していませんでした・・・）
（コミット `6f05eeb` では、 `g++ 9.1.0` でプログラムが途中で終了してしまうエラーを修正）


# 準備
Github からダウンロードした状態では以下のようなディレクトリ構成になっています：

```
$ tree ./
./
├── （略）
├── include
│   └── helpme_standalone.h
├── makefile_md_alad_ewald
├── md_alad_ewald.cpp
└── sample
    ├── config
    ├── input.pdb
    └── water.psf
```
`sample` ディレクトリの中にはサンプルインプットファイルを置いてあります。
準備として、

```
$ mkdir include/uiuc
```
としてください。するとディレクトリの状態は、

```
$ tree ./
./
├── （略）
├── include
│   ├── helpme_standalone.h
│   └── uiuc
├── makefile_md_alad_ewald
├── md_alad_ewald.cpp
└── sample
    ├── config
    ├── input.pdb
    └── water.psf
```
のようになります。

さらに、 CHARMM22 のパラメータファイル `par_all22_prot.inp` を [WEB サイト](http://mackerell.umaryland.edu/charmm_ff.shtml) で探して（toppar_c36_xxx.tgz の中に入っていると思う）、 `sample` ディレクトリに入れてください。

# 必要な外部ファイル - Eigen

このプログラムから、 C++ 線形代数ライブラリ [Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page) を採用しました。適当な場所に Eigen を展開し、上で作成した `include` ディレクトリの中に `eigen-eigen-5a0156e40feb/Eigen` （5a0...の文字列は適宜変更）へのリンクを作成してください。するとディレクトリ構造は、

```
$ tree ./
./
├── （略）
├── include
│   ├── Eigen -> /path/to/eigen-eigen-5a0156e40feb/Eigen/
│   ├── helpme_standalone.h
│   └── uiuc
├── makefile_md_alad_ewald
├── md_alad_ewald.cpp
└── sample
    ├── config
    ├── input.pdb
    └── water.psf
```
のようになります。
 

# 必要な外部ファイル - MDTraj
さらに、[MDTraj](https://github.com/mdtraj/mdtraj) のソースファイルから以下の 5 つを探して `include/uiuc` ディレクトリに入れてください。

```
dcdplugin.h   fastio.h      molfile_plugin.h
endianswap.h  largefiles.h  vmdplugin.h
```
また、 `dcdplugin.c` を `dcdplugin.cpp` という名前に変更してソースコードを展開したディレクトリの中に置いてください。するとディレクトリの状態は、

```
$ tree ./
./
├── （略）
├── dcdplugin.cpp
├── include
│   ├── Eigen -> /path/to/eigen-eigen-5a0156e40feb/Eigen/
│   ├── helpme_standalone.h
│   └── uiuc
│       ├── dcdplugin.h
│       ├── endianswap.h
│       ├── fastio.h
│       ├── largefiles.h
│       ├── molfile_plugin.h
│       └── vmdplugin.h
├── makefile_md_alad_ewald
├── md_alad_ewald.cpp
└── sample
    ├── config
    ├── input.pdb
    └── water.psf
```
のようになります。


# 必要な外部ファイル - helPME

particle mesh Ewald の計算ライブラリ [helPME](https://github.com/andysim/helpme) を使用します。

最新バージョンでは `include` ディレクトリの中にヘッダファイル（`helpme_standalone.h`）が同梱されているので特にすることはありません。

# コンパイル

```
$ make -f makefile_md_alad_ewald
```
Linux 用のインテルコンパイラ（バージョン 18.0.1）でコンパイルできることを確認しています。 MKL を使用する前提で makefile を書いています。
GNU コンパイラは `g++` バージョン 5.4.0 と 9.1.0 でコンパイルできることを確認しています。この場合、コンパイルオプションの `-qopenmp` を `-openmp` に変更し、

```
LIBS = -lfftw3
```

としてください（fftw のライブラリがインストールされている必要あり）。
C++11 を採用しているので、古い GNU コンパイラだとコンパイルできません。

# 実行

```
$ cd sample/
$ ../md_alad_ewald.x config
```

再び `time` コマンドで時間を測定すると・・・

```
real    0m0.309s
user    0m0.284s
sys     0m0.005s
```

となり、**NAMD の半分程度の速度まで改善できた**ことがわかります。

ちなみに PME を使わない（`config` ファイルの中で `usePME no` とする）と、

```
real    0m0.666s
user    0m0.662s
sys     0m0.005s
```

およそ 2 倍の時間がかかることがわかります。 PME は偉大ですね。

# 他の改善点

- 分子の内部自由度の運動を考慮（つまり SHAKE しない）（`rigidBonds no`）
- 境界条件なし（`boundaryType NOBC`）
- 速度 Verlet 法（`integrator VVER`）
- トラジェクトリの .dcd 出力（`DCDFreq 1`）
- 最後のスナップショットの NAMD バイナリ出力（.coor, .vel）

# 改善すべき点

- 分子の内部自由度のうち、二面角（dihedral, improper, CMAP）の考慮
- 今のところ、 NVE にするためには Langevin の damping coefficient を 0 にするしかない。これは NVE と NVT で同じ BBK integrator を共有しているため。NVE と NVT で、 integrator を分離すべき？
- NPT 実装
- 並列化