チュートリアル07:REUSを用いた1次元自由エネルギー面の計算のソースを表示

　本章ではREUS(レプリカ交換アンブレラサンプリング)を用いた1次元自由エネルギー面の計算を説明する。 

== 計算する系の作成 ==
　[[チュートリアル01：水-DCM 界面のシミュレーション#各分子モデルを作成する|チュートリアル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：水-DCM 界面のシミュレーション#分子の重なりを除去する|チュートリアル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_system|&BIAS]]を加えればよい。
 === input.nml =================================================================================
 
   (略)
  
 &BIAS
   calculate_bias = .TRUE.
   nmax = 10
 /
 ===============================================================================================
バイアスの設定は、ネームリストの末尾に&biasdataを加えることで行う。&biasdata は[[Hamiltonian Replica Exchange Umbrella Sampling (HREUS)#Generating &biasdata for 2D REUS using biasdata_gen.exe|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  16
 <s>PARABOLIC 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</s>
 &end
 ===============================================================================================

== 計算の実行①（系の平衡化） ==
　[[チュートリアル02：並列実行|チュートリアル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を行うためにはネームリストに[[Hamiltonian Replica Exchange Umbrella Sampling (HREUS)|&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 System|&MONITOR]]を加えればよい。
 === input.nml =================================================================================
 
   (略)
  
 &MONITOR
   interval = 1
   file = "z.output"
   condition = "REUS_qU 1"
   separator = "$"
 /
 ===============================================================================================

== 計算の実行②（REUS計算） ==
　[[チュートリアル02：並列実行|チュートリアル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を[[Hamiltonian Replica Exchange Umbrella Sampling (HREUS)#Using WHAM to analyze data|&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 &

== 計算の結果 ==
　上記の計算を実行することによって以下のような一次元の自由エネルギー面が得られる。
以下の図において界面の位置は <small><math>z=0</math></small> である。

[[ファイル:T7_free_energy_1d.png|500px]]