チュートリアル07:REUSを用いた1次元自由エネルギー面の計算
本章ではREUS(レプリカ交換アンブレラサンプリング)を用いた1次元自由エネルギー面の計算を説明する。
計算する系の作成
チュートリアル01で用いた系にClイオンを加えた系を考える。このClイオンと水層の重心とのz方向の距離について1次元自由エネルギーを計算する。 Clイオンのmsファイルを以下に示す。
=== cl.ms ===================================================================================== &system version = 0.3 n = 1 / # version: 読み取り形式のバージョン # n: 全サイト数 # cf. Dang, L.X. J.Phys.Chem.B 2001,105,804-809 SITE_DATA # サイトデータ開始点 # 通し番号 サイト名 グループ番号 質量 電荷 分極率 LJ_sigma LJ_eps 座標xyz 速度xyz 力xyz 1 Cl 1 35.5 -1.0 3.690 4.410 0.100 0.0000 0.0000 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 BOND_DATA # 結合データ開始点 # サイト１の番号 サイト２の番号 keyword パラメータのセット 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 データ終了点 ================================================================================================
以下のadd.nmlを用いて、系を作成し、チュートリアル01と同様にしてunlapする。
=== add.nml ==================================================================================== （略） &add file="cl.ms" n=1 area="fix 0.0 0.0 20.0" velocity ="fix 0.0 0.0 0.0" / ================================================================================================
ネームリストの作成①（系の平衡化）
一般化座標zの設定
Clイオンと水層の重心とのz方向の距離を一般化座標zとして定義する。これをするには、ネームリストに&addCOMと&spcoordを加えればよい。
=== input.nml ================================================================================ &FreeFlex (略) / &addCOM condition = "OW HW" / &spcoord condition = "Z 877 878" / ===============================================================================================
ここで、"877"はClイオンの番号（今回の系には877個の分子が含まれており、Clイオンはその877番目にある）、"878"は&addCOM によって与えられた水層の重心番号である。
バイアスの設定
バイアスを用いるにはネームリストに&BIASを加えればよい。
=== input.nml ================================================================================= (略) &BIAS calculate_bias = .TRUE. nmax = 10 / ===============================================================================================
バイアスの設定は、ネームリストの末尾に&biasdataを加えることで行う。&biasdata はbiasdata_gen.exeを用いることで作成できる。
系の数に対応する&biasdataを指定する必要がある。（PBSの都合上、系の個数は16の倍数が好ましいため、17個目以降のバイアスは無視する。）
=== input.nml ================================================================================= (略) &biasdata PARABOLIC 1 1 2 1.0d-10 0.1000000000E-08 1 PARABOLIC 1 1 2 1.0d-10 0.1100000000E-08 2 PARABOLIC 1 1 2 1.0d-10 0.1200000000E-08 3 PARABOLIC 1 1 2 1.0d-10 0.1300000000E-08 4 PARABOLIC 1 1 2 1.0d-10 0.1400000000E-08 5 PARABOLIC 1 1 2 1.0d-10 0.1500000000E-08 6 PARABOLIC 1 1 2 1.0d-10 0.1600000000E-08 7 PARABOLIC 1 1 2 1.0d-10 0.1700000000E-08 8 PARABOLIC 1 1 2 1.0d-10 0.1800000000E-08 9 PARABOLIC 1 1 2 1.0d-10 0.1900000000E-08 10 PARABOLIC 1 1 2 1.0d-10 0.2000000000E-08 11 PARABOLIC 1 1 2 1.0d-10 0.2100000000E-08 12 PARABOLIC 1 1 2 1.0d-10 0.2200000000E-08 13 PARABOLIC 1 1 2 1.0d-10 0.2300000000E-08 14 PARABOLIC 1 1 2 1.0d-10 0.2400000000E-08 15 PARABOLIC 1 1 2 1.0d-10 0.2500000000E-08 16PARABOLIC 1 1 2 1.0d-10 0.2600000000E-08 17 PARABOLIC 1 1 2 1.0d-10 0.2700000000E-08 18 PARABOLIC 1 1 2 1.0d-10 0.2800000000E-08 19 PARABOLIC 1 1 2 1.0d-10 0.2900000000E-08 20 PARABOLIC 1 1 2 1.0d-10 0.3000000000E-08 21 PARABOLIC 1 1 2 1.0d-10 0.3100000000E-08 22 PARABOLIC 1 1 2 1.0d-10 0.3200000000E-08 23 PARABOLIC 1 1 2 1.0d-10 0.3300000000E-08 24&end ===============================================================================================
計算の実行①（系の平衡化）
チュートリアル02と同様にして並列計算を行う。100ps ほど平衡化を行えばよい。
> mpiexec –machinefile machinefile –n 16 ./FreeFlex.exe input.nml &
PBSのスクリプトを用いる際は
#PBS -l nselect=1:ncpus=16:mpiprocs=16:ompthread=1
とすればよい。
ネームリストの作成② （REUS計算）
REUSの設定
REUSを行うためにはネームリストに&REUSを加えればよい。
=== input.nml ================================================================================= (略) &REUS calculate_reus = .TRUE. dim = 1 max_iter = 100 exchange_interval = 100 bias_set_id = 1 dim_q = "1" pos_PARA = "2" spcoord_id = "1" / ===============================================================================================
データ出力の設定
後でWHAM (重み付きヒストグラム解析法)を行うためのデータを出力する。これを行うためには、ネームリストに&MONITORを加えればよい。
=== input.nml ================================================================================= (略) &MONITOR interval = 1 file = "z.output" condition = "REUS_qU 1" separator = "$" / ===============================================================================================
計算の実行②（REUS計算）
チュートリアル02と同様にして並列計算を行う。
> mpiexec –machinefile machinefile –n 16 ./FreeFlex.exe input.nml &
PBSのスクリプトを用いる際は
#PBS -l nselect=1:ncpus=16:mpiprocs=16:ompthread=1
とすればよい。
WHAMの実行
アンブレラサンプリングを行ったデータを解析するにはWHAMを用いる必要がある。WHAMを実行するためには、先のネームリスト"input.nml"の&MONITORを&WHAMSETTINGと&WHAM_INPUTに変えたネームリストを作らなければならない。
=== wham.nml ================================================================================== (略) &WHAMSETTING wham_mode = .TRUE. mpi_size = 16 / &WHAM_INPUT wham_mode = .TRUE. reus_mode = .TRUE. show = .TRUE. input_file_name = "z.output" free_energy_output_file_name = "free_energy.output" histogram_output_file_name = "histogram.out" nbin = 200 ndir = 16 ncore = 16 nprof_column = 1 data_column = "3" spcoord_ID = "1" bias_set_ID = 1 interval = 1 nskip = 14 nread = 1000000 / ===============================================================================================
また、この時ネームリスト&FreeFlexのms_init_fileで指定するファイルについて注意が必要である。ここではwham.nmlと同じディレクトリにあるsystemのmsファイルであれば何を使ってもよい。これは、whamを実行する際には&FreeFlexの情報は基本的に必要がないからである。ただし、FreeFlex.exeを実行するうえで必須なネームリストであるため、最低限、読み込める必要がある。
これを実行する。このときMPIは用いないことに注意する。
> ./FreeFlex.exe wham.nml &
計算の結果
上記の計算を実行することによって以下のような一次元の自由エネルギー面が得られる。 以下の図において界面の位置は である。
