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

Github にアップロードした水（TIP3P）の古典 MD プログラムですが、今までに 3 回に分けて解説を行ってきました。

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

このプログラム（コミット： 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 2.13 （シングルコア）で実行すると・・・

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

となり、自作プログラムの 134 倍速いことがわかります。
そこでもう少し使い勝手をよくするために、自作プログラムを整頓しました。プログラムは Github にアップロードしてあります。

以下の解説はコミット: 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 サイト で探して（toppar_c36_xxx.tgz の中に入っていると思う）、 sample ディレクトリに入れてください。

必要な外部ファイル - Eigen

このプログラムから、 C++ 線形代数ライブラリ Eigen を採用しました。適当な場所に 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 のソースファイルから以下の 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 を使用します。

最新バージョンでは 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 実装
• 並列化