Hirano: ページの作成:「　本章ではFreeFlexを用いたMD計算の実行までを説明する。 == FreeFlex のビルド == 　FreeFlex には 2026/03/16 現在、 FreeFlex/Makefile、Fre…」

2026-05-26T02:34:21Z

ページの作成:「　本章ではFreeFlexを用いたMD計算の実行までを説明する。 == FreeFlex のビルド == 　FreeFlex には 2026/03/16 現在、 FreeFlex/Makefile、Fre…」

新規ページ

　本章では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 が煩雑になってきているので一度整備する必要があるかもしれない。

<span style="color:red">デバッグモードでコンパイルしたFreeFlex.exeは計算を行う際かなり遅くなるので、実際に計算を流す際はmake clean後にmakeし直す必要がある。</span>

分子科学研究所の計算科学研究センターでのコンパイル時には、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未満については[[チュートリアル01：水-DCM 界面のシミュレーション(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"
/

==============================================================================================
[[ファイル:position.jpg|right|分子の配置]]
　はじめの &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のようなビュアーで構造を確認すれば、下記のような構造が得られる。
[[ファイル:system1xyz.jpg|center]]

== 分子の重なりを除去する ==
　作成した構造には分子間の重なりが存在し、このままMDを実行するとエネルギーが発散してしまう。そのため重なりを除去するため Packmol を用いて原子配置を作り直して unlap する。[http://www.ime.unicamp.br/~martinez/packmol/home.shtml 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 の構造は下図のようになり、重なりが除去されていることが確かめられる。（分かりにくいかもしれないが、ビュアー上で動かせば重なりがないことが分かる。）
[[ファイル:Unlap_system1.jpg|center]]

最後に、先ほど作成された 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

のように実行することで、
[[ファイル:show_energy.jpg|center]]
のようなグラフが得られる。保存量 H は中央付近の赤線でほぼ一定になっていることが分かる（一番上の赤線はUindで紛らわしいので注意）。

　xyz_traj_file にはMDのトラジェクトリが保存される。VMD等のビュア―で分布や動きに異常がないか確かめること。

　ms_last_file には最終構造のmsファイルが出力される。ネームリストのms_init_file にこのファイル名を指定することで最終構造の続きからMDを実行できる。また、ms_restart_file には定期的にmsファイルの構造が出力されるので異常終了等の場合はこちらから再開すること。

チュートリアル01：水-DCM 界面のシミュレーション - 版の履歴

Hirano: ページの作成:「 本章ではFreeFlexを用いたMD計算の実行までを説明する。 == FreeFlex のビルド == FreeFlex には 2026/03/16 現在、 FreeFlex/Makefile、Fre…」

Hirano: ページの作成:「　本章ではFreeFlexを用いたMD計算の実行までを説明する。 == FreeFlex のビルド == 　FreeFlex には 2026/03/16 現在、 FreeFlex/Makefile、Fre…」