チュートリアル01：水-DCM 界面のシミュレーション
本章ではFreeFlexを用いたMD計算の実行までを説明する。
FreeFlex のビルド
FreeFlex には 2026/03/16 現在、 FreeFlex/Makefile、FreeFlex/tools/Makefile, FreeFlex/fftw3/Makefile の4 つのMakefile が存在する。 fftw については初回に
> cd fftw3/ > make
を実行しておけば以後は気にする必要はない。その他2つのMakefileについては、作成したいプログラムに対応するソースファイルが存在するディレクトリに移動し、そのディレクトリ内のMakefileの始めの方にある “target: ” で始まる文にプログラム名を指定し make コマンドを実行することで作成される。（例えば、Freeflex.exeを作成したい場合、Makefile 中の“target: ” で始まる文に「Freeflex.exe」と書き込んだ後、 makeコマンドを実行すればよい。）また、再度初めからビルドをやり直したい場合は、
> make clean
デバッグモードでコンパイルしたい場合（通常は書き換えたファイルのみデバッグモードでコンパイルされる。全てのファイルをデバッグしたい場合は先にmake clean を行う）は、
> make debug
を実行することで自動処理を行ってくれる。 その他、Makefile 中にはMPI を使用するかどうかなどの設定項目があるが現在機能していないものも多いため、使用する場合は実際に何を行っているか確認することをお勧めする。 また、Makefile が煩雑になってきているので一度整備する必要があるかもしれない。
デバッグモードでコンパイルしたFreeFlex.exeは計算を行う際かなり遅くなるので、実際に計算を流す際はmake clean後にmakeし直す必要がある。
分子科学研究所の計算科学研究センターでのコンパイル時には、makeの代わりに
> make rccs
を用いることで適切なコンパイルオプションを選択できる。
東北大学金属材料研究所の計算材料学センターでのコンパイル時には、コンパイラがcrayコンパイラになっているので、
> module switch PrgEnv-cray PrgEnv-intel > module module load cray-fftw > module module unload cray-libsci
上記のコマンドでintelコンパイラに変更する。その後、makeの代わりに
> make imr
を用いることで適切なコンパイルオプションを選択できる。
各分子モデルを作成する
まず、作業を行う前に/sample/tutorial/01をFreeFlex外部にコピーをし、そこで作業することを推奨する。
FreeFlexでは分子の情報は msファイルで管理される。 そのため初めに各分子の msファイルを作成することになる。 以下に水分子と DCM分子の msファイルを示す(これらのファイルは sample/ms_0.3_sample/ に用意してあるものと同じ)。なお、msファイル ver0.3未満についてはリンク先を参照すること。
=== water.ms ===================================================================================== &system version = 0.3 n = 3 / # version: 読み取り形式のバージョン（必須） # n: 全サイト数（必須） # cf. Caldwell,J.W.; Kollman,P.A. J.Phys.Chem. 1995, 99, 6208 # water=POL3 SITE_DATA # サイトデータ開始点 # 通し番号 サイト名 グループ番号 質量 電荷 分極率 LJ_sigma LJ_epsilon 座標 速度 力 1 HW 1 1.0 0.365 0.170 0.0 0.0 9.625597 6.787278 12.673000 0.0 0.0 0.0 0.0 0.0 0.0 2 HW 1 1.0 0.365 0.170 0.0 0.0 9.625597 8.420323 12.673000 0.0 0.0 0.0 0.0 0.0 0.0 3 OW 1 16.0 -0.730 0.528 3.204 0.156 10.203012 7.603800 12.673000 0.0 0.0 0.0 0.0 0.0 0.0 BOND_DATA # 結合データ開始点 # サイト１の番号 サイト２の番号 keyword パラメータのセット 1 2 CONSTRAINT 1.63328 1 3 CONSTRAINT 1.0 2 3 CONSTRAINT 1.0 BOND_DATA_END # 結合データ終了点 ANGLE_DATA # 結合角データ開始点 # サイト１の番号 サイト２の番号 サイト３の番号 keyword パラメータのセット ANGLE_DATA_END # 結合角データ終了点 DIHEDRAL_DATA # 二面角データ開始点 # サイト１の番号 サイト２の番号 サイト３の番号 サイト４の番号 keyword パラメータのセット DIHEDRAL_DATA_END # 二面角データ終了点 OUTOFPLANE_DATA # 面外角データ開始点 # サイト１の番号 サイト２の番号 サイト３の番号 サイト４の番号 keyword パラメータのセット OUTOFPLANE_DATA_END # 面外角データ終了点 PAIRTYPE_DATA # ユーザ定義の pairtype データ開始点 # サイト１の番号 サイト２の番号 pairtypeパラメータ PAIRTYPE_DATA_END # ユーザ定義の pairtype データ終了点 ==============================================================================================
=== DCM.ms ===================================================================================== &system version = 0.3 n = 5 / # version: 読み取り形式のバージョン（必須） # n: 全サイト数（必須） # cf. Dang, L.X. J.Chem.Phys. 1999, 110, 10113 SITE_DATA # サイトデータ開始点 # 通し番号 サイト名 グループ番号 質量 電荷 分極率 LJ_sigma LJ_eps 座標xyz 速度xyz 力xyz 1 C 1 12.0 -0.2720 0.878 3.410 0.137 0.0000 0.0000 0.7761 0.0 0.0 0.0 0.0 0.0 0.0 2 Cl 1 35.5 -0.0537 1.910 3.450 0.280 0.0000 1.4675 -0.2178 0.0 0.0 0.0 0.0 0.0 0.0 3 Cl 1 35.5 -0.0537 1.910 3.450 0.280 0.0000 -1.4675 -0.2178 0.0 0.0 0.0 0.0 0.0 0.0 4 H 1 1.0 0.1897 0.135 2.400 0.040 -0.8854 0.0000 1.3734 0.0 0.0 0.0 0.0 0.0 0.0 5 H 1 1.0 0.1897 0.135 2.400 0.040 0.8854 0.0000 1.3734 0.0 0.0 0.0 0.0 0.0 0.0 BOND_DATA # 結合データ開始点 # サイト１の番号 サイト２の番号 keyword パラメータのセット 1 2 CONSTRAINT 1.772 1 3 CONSTRAINT 1.772 1 4 CONSTRAINT 1.070 1 5 CONSTRAINT 1.070 2 3 CONSTRAINT 2.92420 2 4 CONSTRAINT 2.33057 2 5 CONSTRAINT 2.33057 3 4 CONSTRAINT 2.33057 3 5 CONSTRAINT 2.33057 BOND_DATA_END # 結合データ終了点 ANGLE_DATA # 結合角データ開始点 # サイト１の番号 サイト２の番号 サイト３の番号 keyword パラメータのセット ANGLE_DATA_END # 結合角データ終了点 DIHEDRAL_DATA # 二面角データ開始点 # サイト１の番号 サイト２の番号 サイト３の番号 サイト４の番号 keyword パラメータのセット DIHEDRAL_DATA_END # 二面角データ終了点 OUTOFPLANE_DATA # 面外角データ開始点 # サイト１の番号 サイト２の番号 サイト３の番号 サイト４の番号 keyword パラメータのセット OUTOFPLANE_DATA_END # 面外角データ終了点 PAIRTYPE_DATA # ユーザ定義の pairtype データ開始点 # サイト１の番号 サイト２の番号 pairtypeパラメータ 4 5 10 PAIRTYPE_DATA_END # ユーザ定義の pairtype データ終了点 ==============================================================================================
各ファイルのはじめに記載してあるネームリストでは、msファイルのバージョンとサイト数が指定されており、以下に続く#以下はコメント行となる（以下同様）。
SITE_DATA等のキーワードについては msファイルの仕様を参照すること。
分子を配置する
系全体の情報も msファイルで管理される。 系全体の msファイルは add.exe で作成する。 具体的には、
> ./add2.exe add.nml system1.ms
の形で使用する。 (ここで、記載したコマンドはあくまで例であり、実際にコマンドを実行する際には、ファイルの存在する位置を適宜読み替える必要がある。以下のコマンドの記載も同様である。) ここで add.nmlは add.exe用のネームリストで今回は下記のものを用いる。
=== add.nml ===================================================================================＝= &system size = 25.0 25.0 85.0 Tinit = 298.15 / # output: 出力先ファイル名 # size: 系のx,y,z方向のサイズ（単位はÅ）。はみ出た場合周期境界として扱う。 # Tinit: 系の温度（単位はK）。初期速度の作成に用いる。 &add file="water.ms" n=523 area="cuboid 0.0 0.0 0.0 25.0 25.0 25.0" velocity ="maxwell" direction ="random" / # file: 分子情報が記載されたmsファイル名 # n: 設置する分子数 # area, velocity, direction: 設置の仕方 &add file="DCM.ms" n=353 area="cuboid 0.0 0.0 25.0 25.0 25.0 85.0" / ==============================================================================================
はじめの &system以下では周期境界条件、系の初期温度を指定する。 その後の &add以下ではどの分子を何個、どのように設置するかを指定する。 今回は初めに水分子を523個、続いてDCM分子を353個、右図のような形に設置している。出力として下記が得られ、ディレクトリ内に system1.msが作成される。
（以下で各化学種のファイル "water.ms, DCM.ms" を読み込むときに、「WARNING: &system の変数 size が指定されていない。...」と出ているが、ここでは問題ない。sizeパラメータは、add.nml 中ですでに指定している。）
=== add.exe の実行結果 ============================================================================== Input File: add.nml &SYSTEM SIZE = 2*2.500000000000000E-009 , 8.500000000000000E-009, TINIT = 298.150000000000 / water.ms から 523個追加 system.F(393): WARNING: バージョン0.3以上の ms ファイルを読み込みます。 system.F(400): WARNING: &system の変数 size が指定されていない。デフォルト値:[100.0d0,100.0d0,100.0d0] を使用する。 DCM.ms から 353個追加 system.F(393): WARNING: バージョン0.3以上の ms ファイルを読み込みます。 system.F(400): WARNING: &system の変数 size が指定されていない。デフォルト値:[100.0d0,100.0d0,100.0d0] を使用する。 system1.ms に出力 ==============================================================================================
作成した構造を確認する
作成したmsファイルはms2xyz.exeでxyzファイルに変換できる。
> ./ms2xyz.exe system1.ms system1.xyz
上記コマンドによりディレクトリ内にsystem1.xyzが作られる。後はvmdのようなビュアーで構造を確認すれば、下記のような構造が得られる。
分子の重なりを除去する
作成した構造には分子間の重なりが存在し、このままMDを実行するとエネルギーが発散してしまう。 そのため重なりを除去するため Packmol を用いて 原子配置を作り直して unlap する。PACKMOLウェブサイトより Packmol をダウンロードし、User guide に従ってコンパイルすること。
まず 以下のような xyzファイルを用意する。 これらは、それぞれ
> ./ms2xyz.exe water.ms water.xyz > ./ms2xyz.exe DCM.ms DCM.xyz
で作成できる。
=== water.xyz ====================================================================================
3
made from ms2xyz
HW 9.62560 6.78728 12.67300
HW 9.62560 8.42032 12.67300
OW 10.20301 7.60380 12.67300
==============================================================================================
=== DCM.xyz =====================================================================================
5
made from ms2xyz
C 0.00000 0.00000 0.77610
Cl 0.00000 1.46750 -0.21780
Cl 0.00000 -1.46750 -0.21780
H -0.88540 0.00000 1.37340
H 0.88540 0.00000 1.37340
==============================================================================================
次に作成する系の情報を記した system.inp を用意する。Packmol は周期境界条件に非対応であるため、本来のボックスサイズから 1Å ずつ内側にずらして設定する。これにより生じる空間はシミュレーション初期の平衡化によって消失するため問題ない。
=== system.inp ================================================================================== tolerance 2.0 filetype xyz output system.xyz structure water.xyz number 523 inside box -11.5 -11.5 -24.0 11.5 11.5 0.0 end structure structure DCM.xyz number 353 inside box -11.5 -11.5 0.0 11.5 11.5 59.0 end structure ==============================================================================================
以下のコマンドにより、重なりが除去された構造が記録された system.xyz が作成される。
> ./packmol < system.inp
system.xyz の構造は下図のようになり、重なりが除去されていることが確かめられる。（分かりにくいかもしれないが、ビュアー上で動かせば重なりがないことが分かる。）
最後に、先ほど作成された system1.ms の座標部分を削除し、Packmol で作成した system.xyz の座標部分で置換することで unlapされた msファイル system1.ms を生成する。この操作は、以下のコマンドで実行できる。（ここでは、xyz2ms.exe の出力ファイル system2.ms で元の system1.ms を上書きしている。）
> ./xyz2ms.exe system.xyz system1.ms system2.ms > mv system2.ms system1.ms
MDシミュレーションを行う
初期構造の作成が終了したので本題のMD計算に入る。MD計算はFreeFlex.exeにて行う。 使い方は以下の通り。
> ./FreeFlex.exe input.nml &
input.nmlはFreeFlex.exe用のネームリストで今回は下記を用いる。 各設定の意味については&FreeFlexを参照すること。
=== input.nml ====================================================================================== &FreeFlex ms_init_file = "system1.ms", xyz_traj_file = "traj.xyz", log_file = "log", ms_last_file = "system2.ms", nml_output_file = "output.nml" ms_restart_file = "restart.ms" T = 298.15 nstep = 10000 dt = 1.0d-15 / ================================================================================================
上記コマンドを実行するとMD計算が開始される。 大量の警告メッセージが出力されるが、それほど気にしなくてよい。 これらのメッセージは、初期構造の不安定に伴うものだからである。数千ステップほど経過すると、通常は警告メッセージが表示されなくなる。 表示され続ける場合は何らかの問題を含んでいるので注意すること。
結果を見る
MDシミュレーションの結果はネームリスト（input.nml）で指定した各ファイルに出力される。 log_file にはMD計算のログとして下記に示す通り初期設定や各ステップの熱力学量等に関する情報が出力される。
=== log ========================================================================================
Start FreeFlex at 2015年 6月 1日 月曜日 17:24:49 JST
Read namelist file "input.nml"
random_seed = 1483 2237
Read ms init file.
Setting NVT simulation with Nose-Hoover Chain.
Chain number set to: 10
### cutoff parameters ###
rcut = 1.150000000000000E-009
rskin= 1.000000000000000E-010
rswitch= 6.000000000000000E-010
rswitch_0= 1.000000000000000E-010
### Ewald parameters ###
sigma = 2.187758465334345E-010
kappa = 3232106251.17377
hmax = 10 10 33
.
.
.
### Periodical Output ###############################################
step = 4000
H = -2.73810111E-17, T = 316.242
U = -5.09978E-17, K = 1.14678E-17
ULJ = 3.083E-19, UCoulomb =-5.131E-17
UCC =-3.968E-17
UCD =-1.162E-17, UDC =-1.162E-17, UDD =-2.912E-18, Uind = 1.453E-17
SHAKE loop: 25, RATTLE loop: 22
dipole loop: 5
total momentum = 7.88738E-36 1.91973E-35 -9.58818E-36
#####################################################################
### Periodical Output ###############################################
step = 4050
H = -2.73805442E-17, T = 307.091
U = -5.07263E-17, K = 1.11360E-17
ULJ = 5.668E-19, UCoulomb =-5.129E-17
UCC =-3.962E-17
UCD =-1.168E-17, UDC =-1.168E-17, UDD =-2.952E-18, Uind = 1.463E-17
SHAKE loop: 24, RATTLE loop: 21
dipole loop: 5
total momentum = 8.44042E-36 1.81884E-35 -5.98198E-36
#####################################################################
.
.
.
total: 5500.316 s
init: 2.501 s
move: 90.051 s
output: 1661.640 s
interact: 3745.717 s
distance: 138.325 s
book: 24.263 s
distance: 125.661 s
LJ: 38.211 s
coulomb: 3555.451 s
tensor: 300.099 s
induced: 2683.162 s
barrier: 679.601 s
sum: 543.062 s
ewald: 2432.536 s
real: 427.172 s
tensor: 62.810 s
charge_field: 49.282 s
dipole_field: 266.958 s
recip: 1797.643 s
expikx: 235.461 s
charge_field: 244.774 s
dipole_field: 1281.753 s
mask: 70.571 s
Finish FreeFlex at 2015年 6月 1日 月曜日 18:56:30 JST
================================================================================================
初期設定（自動決定されるパラメータ）に不備がないか？、保存量（Periodic Output の H）が保存しているか？（今の条件では4桁保存する模様）、エネルギーやループ回数に異常がないか？、初期緩和の時以外にWARNNINGが出ていないか？、計算時間に異常がないか？ 等を確かめることでMDが正常に動いているかを確認すること。 また、各種エネルギーの変化を視覚化するツールとして、tools/show_energy.exe がある。
> ./show_energy.exe log 200
のように実行することで、
のようなグラフが得られる。 保存量 H は中央付近の赤線でほぼ一定になっていることが分かる（一番上の赤線はUindで紛らわしいので注意）。
xyz_traj_file にはMDのトラジェクトリが保存される。VMD等のビュア―で分布や動きに異常がないか確かめること。
ms_last_file には最終構造のmsファイルが出力される。 ネームリストのms_init_file にこのファイル名を指定することで最終構造の続きからMDを実行できる。 また、ms_restart_file には定期的にmsファイルの構造が出力されるので異常終了等の場合はこちらから再開すること。