Github にアップロードした水（TIP3P）の古典 MD プログラムですが、今までに 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 にアップロードしてあります。
以下の解説はコミット: 7b19cd56f05eeb780a659 に基づいています。
（コミット 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 実装
- 並列化