(水溶液界面のグランドカノニカルMD)

溶液の界面の系を適用対象として液相のバルクにおける溶質密度を一定とするグランドカノニカルMDの実行方法、特にCoulomb相互作用を含む系を対象とする場合について解説する。 基本的な流れはチュートリアル10と同一。

計算の流れ

密度分布をサンプリングの前段階としてON-OFF間の溶質の自由エネルギー差を解消するBias Potentialを決定する必要がある。 Coulomb相互作用を含む系でのグランドカノニカルMDではこの自由エネルギー差の解消が重く、工夫が求められる計算である。 全体の計算の流れは次の3stepで構成される。

1. 粒子数一定のMD計算による初期Biasの作成

2. グランドカノニカルMDによるBiasの最適化

3. グランドカノニカルMDによる界面密度分布の計算

基本的に熱力学的積分法から得られる自由エネルギー面の正負を反転させたものを複数の解析関数の和でfittingを行いBias Potentialを決定する。 1.では系中の一つの溶質分子のみをGCで取り扱うこととして、溶質座標を拘束・平均力サンプリングを行う。tapering領域内での自由エネルギー面が必要であるため領域を覆えるだけの拘束位置が異なるトラジェクトリを並列計算を行い自由エネルギー面を作成する。溶質種が複数ある場合には溶質種毎にそれぞれ初期Biasを作成する必要がある。

2.ではすべての溶質分子をGCで取り扱う。溶質座標をビンで区切ることで平均力のヒストグラムを作成し、ヒストグラムを元に自由エネルギー面を得る。 最適化ではBiasを適用したGCMDから得られた自由エネルギー面を再度fitting、計算されたBiasを元のBiasに追加しMDを行う操作を繰り返す。 全ての溶質を対象に平均力を計算するため溶質種が複数ある場合でも全溶質種同時に自由エネルギー面を得ることができる。 本チュートリアルでは主に2.のBiasの最適化の方法について取り扱う。

3.ではBiasの最適化が十分に完了した状態で長時間計算を行い、密度分布を求める。3.における密度分布を元に表面過剰量や界面構造の評価を行う。

使用プログラム

グランドカノニカルMDで使用する主なプログラムおよびチュートリアルで使用するスクリプトの一覧。

名称	説明
FreeFlex.exe	MD計算のメインプログラム。
gauss_fit.exe	離散データをGauss関数の和でfittingするプログラム。
biasplot.exe	Gauss fit で作成したバイアスをプロットできる形に出力するプログラム。
tools/gc_analyzer.exe	グランドカノニカルMDの出力ファイルから自由エネルギー面と密度分布を計算するプログラム。
tools/tapstep_fit.exe	離散データをステップ状の連続関数でfittingするプログラム。
tools/therm_int.exe	溶質座標を拘束した平均力サンプリングの出力から自由エネルギー面を求めるプログラム。
gc_gen_tdi_cation.sh	粒子数一定のMD計算による平均力サンプリングの初期構造を生成するスクリプト。cationの初期Biasを作成する場合に使用。
gc_gen_tdi_anion.sh	粒子数一定のMD計算による平均力サンプリングの初期構造を生成するスクリプト。anionの初期Biasを作成する場合に使用。
gc_gen_system_1.sh	グランドカノニカルMDの初期構造（溶媒, ON）を生成するスクリプト。
gc_gen_system_2.sh	グランドカノニカルMDの初期構造（系全体）を生成するスクリプト。

以下のチュートリアルでは上記の他に add2.exe, ms2xyz.exe, xyz2ms.exe および packmol も補助的に使用するため、コンパイルしておくこと。

FreeFlex.exeのグランドカノニカルMD用入力ネームリストに関する 説明はこちらを参照。 また、gauss_fit.exe, gc_analyzer.exe, tapstep_fit.exe, therm_int.exeのネームリスト変数の一覧については各プログラムのページを参照。

初期構造生成用の各種スクリプトはファイル中のパラメータを直接編集して使用する。 add2.exe, ms2xyz.exe, xyz2ms.exeおよびpackmolがスクリプトから呼び出される。 スクリプトの機能と使用方法は以下のチュートリアルを、 パラメータ設定方法などの詳細はスクリプトファイル中のコメントを参照。

初期Biasの作成

Coulomb相互作用を含む系として、電解質水溶液を対象にグランドカノニカルMDを行う。 溶質はナトリウムイオン・ヨウ化物イオン、溶媒はD2Oである。水を重水素置換しているのはtime stepを広げる工夫である。 これらを含むスラブ系に対して、まず各溶質種の平均力ポテンシャルを、溶質座標を拘束した計算によってサンプリングする。

初期Bias作成に用いるファイルは aqueous/init_bias/init_bias_* の中に用意されている。 粒子数一定のMDを行うためのinputファイルを以下の平衡化を目的とするinput_0.nmlとinput_1.nml、サンプリングを目的とするinput_2.nmlのように用意する。 input_0.nmlでは分極計算を省いた計算をすることで、発散して系が崩れるのを防ぐ。 input_0.nmlはinput_1.nmlをコピーして元のファイルを作り、この内入出力するファイル名の番号を全て1つずつ下げ、nstepを10000に、calculate_polを.FALSE.にする。 input_1.nml を用いて平衡化計算を行う際には最初界面に対して壁を設置し、蒸発を防いだ状態で計算を行うためset_wallによって壁が指定されている。 ただし壁を置く時間はnstep_wallによって指定され、以下の例では壁付きの平衡化を60 ps, 次いで壁なしの平衡化を40 ps行う。

平均力サンプリングでは拘束位置の異なるトラジェクトリを複数用意し（ここでは128個）、それぞれで平均力サンプリングを行う。 初期構造の生成はaqueous/init_bias/init_bias_*においてgc_gen_tdi_cation.shまたはgc_gen_tdi_anion.shによって行う。実行前にファイルを開き、#出力ファイル名の部分で「system1.ms」を「system0.ms」に書き換え、さらに「inputnml_0="input_0.nml"」の行を追加する。 また、#一時ファイル名の部分に「inputnml__=$dir_"/"$inputnml_0」「cp $inputnml_ $inputnml__」の2行を追加する。（すぐ上に書いてある#入力ネームリストファイル2をコピペして書き換えるだけである） (スクリプト内のadd2、ms2xyz、xyz2msのpathは適宜修正が必要) スクリプト実行後には拘束位置の数の分だけ 0000/, 0001/, ...とディレクトリが生成され、その下に初期構造となるmsファイルとnmlファイルが二つ置かれる。 各トラジェクトリでは溶質の拘束位置が異なり、mpiディレクトリ下のinput_[1 or 2].nml 内の &addconstr によって指定される。 ただし平衡化とサンプリングで拘束位置は変化しないため、これらの内容は同一である。

input_1.nml

$ cat input_1.nml 
&FreeFlex
 ms_init_file = "system1.ms",
 xyz_traj_file = "traj_1.xyz",
 log_file = "log_1",
 ms_last_file = "system2.ms",
 nml_output_file = "output_1.nml"
 ms_restart_file = "restart_1.ms"
 T = 298.15
 nstep = 50000
 dt = 2.0d-15
 omp_number = 1
 ewald_mode = "spme"
 pme_mode = .TRUE.
 calculate_real= .TRUE.
 calculate_mask = .TRUE.
 calculate_recip = .TRUE.
 calculate_pol = .TRUE.
 calculate_damp = .FALSE.
 calculate_neut = .TRUE.
 calculate_surf = .FALSE.
 interact_type = "gc"
 output_energy_interval = 100
 output_traj_interval = 2500
 output_restart_interval = 10000
 pme_ngrid=32,32,192
 v_correct = .TRUE.
 rdf_calc = .TRUE.
 induced_Ethresh = 5.0d11
 rswitch = 0.0d0
 set_wall = .TRUE.
 z_wall = "-16.0d-10 16.0d-10"
 nstep_wall = 30000
/
&addCOM
 condition = "HW OW"
/
&spcoord
 condition = "Z 545 546"
/
&constraint
 tolerance = 1.0d-16
 maxiter = 200
 threshold = 0.2d0
/
&MONITOR
 interval = 1
 file = "spcoord_f_2.dat"
 condition = "SPCOORD_F 1"
 separator = "$"
/
&MONITOR
 interval = 1
 file = "spcoord_q_2.dat"
 condition = "SPCOORD_Q 1"
 separator = "$"
/
&gc
 perform_gc = .TRUE.
 do_off_dummy_switch = .FALSE.
 COMid_solvent = 546
 spcoord_type = "GC_Z_DQDX_CONST"
 type_tapering = "6"
 delta_tapering = "3.0d-10"
 qswitch = "0.0d0"
 coulomb_sigma = 5.0d-10
 c_corr = 100
 sigma_corr = 0.5d-10
/
&gcbias
 calculate_bias = .FALSE.
/

実行はPBS proを利用が望ましく、以下のようにPBSへscriptファイルを投入すればよい。 なお、FreeFlexのpath等は適宜修正すること。 また、実行ファイルをinput_0.nmlとした上で、その計算を実行した後にinput_1.nmlの計算を行う。

$qsub pbssub.sh

input_0.nmlの計算を行わなかった場合、次のことが起こり得る。 初期構造生成後の平衡化において分極が発散するなどして系の温度が異常に高温になり系が壊れることが10%程度の割合で起こる。これはGCMDにおいて近距離のswitchingであるrswitchを0にしており分極の発散が起きやすいことが原因であるとみられる。 その際には新たに系を作成、または系が壊れていない中から拘束位置が近いものを初期構造として採用し、再度平衡化計算を行うとよい。系の崩壊を防止する策としてはじめに分極を切って計算を行うことや拘束する溶質含め全分子通常のNVTで取り扱いrwitchをデフォルト値にした状態で緩和させることが考えられる。

系が壊れた際はlogファイルに以下のような警告メッセージが頻繁に出力される。このような場合は上述の対応をとり平衡化を行う。

$tail -n 15 0018/log_1 | grep -i warn
default.F(1759): WARNING: 誘起双極子が収束していません.   max(|mu-mupre|) =    0.1806849327E-34, mpirank:        18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =         14, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =         25, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =         47, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        151, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        191, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        298, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        362, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        397, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        473, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        517, mpirank:         18
move.F(400): WARNING: 力が大きすぎるため補正しました。 i =        518, mpirank:         18

平均力サンプリングを実行する際のinputファイルは以下の通り。pbssub.sh内の実行ファイルをinput_2.nmlと変更したうえで同様にPBSへ投入すればよい。

input_2.nml

$ cat input_2.nml 
&FreeFlex
 ms_init_file = "system2.ms",
 xyz_traj_file = "traj_2.xyz",
 log_file = "log_2",
 ms_last_file = "system3.ms",
 nml_output_file = "output_2.nml"
 ms_restart_file = "restart_2.ms"
 T = 298.15
 nstep = 50000
 dt = 2.0d-15
 omp_number = 1
 ewald_mode = "spme"
 pme_mode = .TRUE.
 calculate_real= .TRUE.
 calculate_mask = .TRUE.
 calculate_recip = .TRUE.
 calculate_pol = .TRUE.
 calculate_damp = .FALSE.
 calculate_neut = .TRUE.
 calculate_surf = .FALSE.
 interact_type = "gc"
 output_energy_interval = 100
 output_traj_interval = 2500
 output_restart_interval = 10000
 set_wall = .FALSE.
 v_correct = .TRUE.
 n_restart_file = 5
 rdf_calc = .TRUE.
 rswitch = 0.0d0
 rcut_coulomb=10.5d-10
 pme_ngrid=32,32,192
 induced_Ethresh = 5.0d11
 pme_sporder=8
/
&addCOM
 condition = "HW OW"
/
&spcoord
 condition = "Z 545 546"
/
&constraint
 tolerance = 1.0d-16
 maxiter = 500
 threshold = 0.2d0
/
&MONITOR
 interval = 1
 file = "spcoord_f_2.dat"
 condition = "SPCOORD_F 1"
 separator = "$"
/
&MONITOR
 interval = 1
 file = "spcoord_q_2.dat"
 condition = "SPCOORD_Q 1"
 separator = "$"
/
&gc
 perform_gc = .TRUE.
 do_off_dummy_switch = .FALSE.
 COMid_solvent = 546
 spcoord_type = "GC_Z_DQDX_CONST"
 type_tapering = "6"
 delta_tapering = "3.0d-10"
 qswitch = "0.0d0"
 coulomb_sigma = 5.0d-10
 c_corr = 100
 sigma_corr = 0.5d-10
/
&gc_bias
 calculate_bias = .FALSE.
/

$qsub pbssub.sh

サンプリング終了後therm_int.exeを実行し、一次元エネルギー面を得る。 拘束の始点と間隔、何step目まで読み飛ばすか与える必要がある。

$~/src/FreeFlex/tools/therm_int.exe ./therm_int.nml

$ cat therm_int.nml 
&thermint
 inputfile = "spcoord_f_2.dat"
 outputfile = "free_energy.dat"
 nskip = "14"
 nread = "330000"
 mpi_number = 128
 start = -3.36d-10
 interval = 0.05d-10
/

therm_int.exeを実行すると、free_energy.dat への出力としてそれぞれ以下のような自由エネルギー面が得られる。

次に得られた自由エネルギー面をfittingし、初期に適用するBiasを決定する。 fittingはtapered_step biasでON-OFF間の自由エネルギー差を大まかに解消し、次にその残差を複数のgauss_1d biasで埋める。

$~/FreeFlex/tools/tapstep_fit.exe ./tapstep_fit.nml
$ ~/FreeFlex/gauss_fit.exe ./gauss_fit.nml

$ cat tapstep_fit.nml 
&tapstepfit
 inputfile = "free_energy.dat"
 biasfile = "bias_na_tap.dat"
 plotfile = "plot_na_tap.dat"
 restfile = "gfinput2.dat"
 center = 0.2d-10
 delta = 1.9d-10
/

cat gauss_fit.nml 
&gauss_fit
 input_file_name = "gfinput2.dat"
 param_file_name = ""
 bias_file_name = "bias_na_gauss.dat"
 plot_file_name = "plot_na_gauss.dat"
 log_file_name = ""
 spcoord_ID = ""
 ngauss_max = 50
 niter_max = 1000
 thres_finish = 1.0d-2
 thres_unchange = 1.0d-2
 lambda_1 = 1.0d-5
 lambda_2 = 1.0d-5
 lambda_3 = 1.0d-4
 ncolumn_data = 1
 ncolumn_weight = 0
 omp_number = 1
 plot_result = .TRUE.
/

gauss_fit.exeを実行する際にはtapstep_fit.nmlで指定したrestfileの先頭の行に以下の情報を追加する必要がある。

$ #plot 1 1 n
  n: dataが存在する行数（この例では n = mpi_number + 1 = 129） 

tapered_stepによるfittingは中心の位置や幅は、Biasをかける領域からはみ出すことがなければさほど支障はないためパラメーターの調整は拘らなく良い。 gauss関数によるfittingはデータが存在しない領域の外、Biasをかける領域外では速やかに収束するように決める必要がある。 外側にピークがあるBiasを生成してしまった場合には試行錯誤的にパラメーターの調整を行い、アーティファクトが形成されないように修正を行う。 外側にピークができていないか確認するにはbiasplot.exeを実行するとよい。

$ ~/FreeFlex/biasplot.exe ./biasplot.nml

作成されたBias（bias_na_tap.dat および bias_na_gauss.dat に出力されている）を適用するにはinput.nml内のnml、&gcbiasに追加する。

$&gcbias
 Na 1
 t1 s1
 TAPERED_STEP 1 1 3   8.5665306443E-019  7.0000000000E-011  1.9000000000E-010 0
 GAUSS_1D   1 1   3  0.8643425869E-19 -0.1721676346E-10  0.6566984025E-10  0
 GAUSS_1D   1 1   3  0.2774969855E-19  0.1681943494E-09  0.3212891372E-10  0
 &end

同様の手順をanion(ヨウ化物イオン)に対して実行し、作成したBias（bias_io_tap.dat および bias_io_gauss.dat に出力されている）を&gcbiasに追加する。

Biasは溶質種毎に異なるため複数溶質種が存在する場合、溶質種毎に&gcbiasを用意する必要がある。 二種類目の溶質(type2)であれば適用されるBiasは&gcbias二行目t2以下で指定されるBiasとなる。

$&gcbias
 Io 1
 t2 s1　
 TAPERED_STEP 1 1 3   5.5975196375E-019  7.0000000000E-011  1.8000000000E-010 0
 GAUSS_1D   1 1   3  0.5717509241E-19 -0.2643199901E-10  0.5732894427E-10  0
 GAUSS_1D   1 1   3  0.2861966463E-19  0.1214287496E-09  0.6142124609E-10  0
 GAUSS_1D   1 1   3 -0.7556229278E-20  0.2131926188E-09  0.3008222160E-10  0
&end

Biasの更新・初期

グランドカノニカルMDの実行し自由エネルギー面の採取・Bias Potentialの更新を繰り返す過程。 特に、上述の初期Biasを適用した直後の自由エネルギー差が大きい状態での更新の例について説明する。 作業するディレクトリはbiasupdate_1であるため当該ディレクトリへ移動すること。

初期構造の生成は、系の破綻を防止するために二段階に分けて行う。初めに溶媒とONの溶質のみの系を作成し、粒子数一定のMD(非GCのMD)を行い緩和させる。 この計算においても、上述した手順で.shファイルを書き換えてsystem0.msやinput0.nmlを作成し、分極しない計算を10000ステップ程度行う。 このとき分極が発散する、警告メッセージが常に出続けるような系があれば該当する系を取り除き、系を作り直すか他ディレクトリの系を複製する [複製してはいけないのでは]。 gc_gen_system1.shの実行の前にmpi並列計算用のディレクトリを用意しなければならないことに注意。

$ for dir in `seq -f %04g 0 79`
 >do
 >mkdir $dir
 >done
 sh gc_gen_system_1.sh
 qsub pbssub.sh

input_1.nml

$ cat input1.nml 
&FreeFlex
 ms_init_file = "system1.ms",
 xyz_traj_file = "traj_1.xyz",
 log_file = "log_1",
 ms_last_file = "system2_on.ms",
 nml_output_file = "output_1.nml"
 ms_restart_file = "restart_1.ms"
 T = 298.15
 nstep = 50000
 dt = 2.0d-15
 omp_number = 1
 ewald_mode = "spme"
 pme_mode = .TRUE.
 calculate_pol = .TRUE.
 calculate_damp = .FALSE.
 calculate_neut = .TRUE.
 calculate_mask = .TRUE.
 calculate_surf = .FALSE.
 interact_type = "default"
 output_energy_interval = 100
 output_traj_interval = 1000
 output_restart_interval = 10000
 set_wall = .TRUE.
 z_wall = "-16.0d-10 16.0d-10"
 nstep_wall = 30000
 v_correct = .TRUE.
 n_restart_file = 5
 rcut_coulomb=11.87d-10
 pme_ngrid=32,32,192
 pme_kappa=0.3062d10
 induced_Ethresh =4.0d11
/
&gc
 perform_gc = .FALSE.
/

平衡化が終了した系に対してOFFとdummyの溶質を追加してグランドカノニカルMDで用いる初期構造を生成、グランドカノニカルMDによるサンプリングを行う。 gc_gen_system_2.shを実行することによってOFFやdummyの溶質の追加がなされる。 与えるinput file内に前節、初期Biasで作成したBiasをそれぞれの溶質に対して適用するよう記述すること。

$sh gc_gen_system_2.sh

Coulomb相互作用は初期Biasとのずれが大きく溶質の密度分布が崩れることがあるため、それを抑制するために&accel_dequilによって密度差解消のBiasを適用する。 &accel_dequilではtapering領域近傍を除く密度分布が平坦となることを期待する領域(count_area)、溶質に対して力をかける領域(force_area)、設定濃度(set_conc)、溶質の移動速度(shift_v)を設定する必要がある。 更新初期では長時間計算を行うと密度分布が崩れるため計算時間を短くし、大まかな自由エネルギー面の把握でよい。

input_2.nml

$ cat input2.nml 
&FreeFlex
 ms_init_file = "system2.ms",
 xyz_traj_file = "traj_2.xyz",
 log_file = "log_2",
 ms_last_file = "system3.ms",
 nml_output_file = "output_2.nml"
 ms_restart_file = "restart_2.ms"
 T = 298.15
 nstep = 30000
 dt = 2.0d-15
 omp_number = 1
 ewald_mode = "spme"
 pme_mode = .TRUE.
 calculate_real= .TRUE.
 calculate_recip = .TRUE.
 calculate_pol = .TRUE.
 calculate_damp = .FALSE.
 calculate_neut = .TRUE.
 calculate_mask = .TRUE.
 calculate_surf = .FALSE.
 interact_type = "gc"
 output_energy_interval = 100
 output_traj_interval = 500
 output_restart_interval = 10000
 v_correct = .TRUE.
 n_restart_file = 5
 rdf_calc=.TRUE.
 rswitch = 0.0d0
 rcut_coulomb=11.87d-10
 pme_ngrid=32,32,192
 pme_kappa=0.3062d10
 induced_Ethresh =3.0d11
/
&addCOM
 condition = "HW OW"
/
&gc
 perform_gc = .TRUE.
 do_off_dummy_switch = .TRUE.
 COMid_solvent = 882
 spcoord_type = "GC_Z_DQDX_CONST"
 type_tapering = "6"
 delta_tapering = "3.0d-10"
 qswitch = "0.0d0"
 coulomb_sigma = 5.0d-10
 c_corr = 100
 sigma_corr = 0.5d-10
/

&gc_bias
 calculate_bias = .TRUE.
 nmax = 10
/
&gc_density
 calculate_density = .TRUE.
 output_filename = "gc_density.dat"
 output_interval = 1000
 lbin = 0.05d-10
/
&gc_free
 calculate_free = .TRUE.
 output_filename = "gc_free.dat"
 output_interval = 1000
 lbin = 0.05d-10
/
&accel_dequil
 do_accel_dequil = .TRUE.
 ntype = 2 
 count_area = 4.5d-10,6.0d-10
 force_area = -3.0d-10,5.0d-10
 set_conc = 1.2d0
 shift_v = 5d1 
/
&gcbias
Na 1
t1 s1
TAPERED_STEP 1 1 3   8.5665306443E-019  7.0000000000E-011  1.9000000000E-010 0
GAUSS_1D   1 1   3  0.8643425869E-19 -0.1721676346E-10  0.6566984025E-10  0
GAUSS_1D   1 1   3  0.2774969855E-19  0.1681943494E-09  0.3212891372E-10  0
&end

&gcbias
Io 1
t2 s1
TAPERED_STEP 1 1 3   5.5975196375E-019  7.0000000000E-011  1.8000000000E-010 0
GAUSS_1D   1 1   3  0.5717509241E-19 -0.2643199901E-10  0.5732894427E-10  0
GAUSS_1D   1 1   3  0.2861966463E-19  0.1214287496E-09  0.6142124609E-10  0
GAUSS_1D   1 1   3 -0.7556229278E-20  0.2131926188E-09  0.3008222160E-10  0
&end

サンプリング終了後、gc_analyzer.exeを実行し自由エネルギー面を求める。自由エネルギー面作成の対象とするディレクトリはdir_listで指定すること。 &free内のnskipとinput_2.nml、&gc_free内のoutput_intevalとdtの積が平衡化の時間、nread*output_interval*dtがサンプリング時間である。

$~/FreeFlex/tools/gc_analyzer.exe ./analyze_free.nml 

$ cat analyze_free.nml 
&free
 is_formatted = .TRUE.
 inputfile = "gc_free.dat"
 outputfile = "free_.dat"
 max_nread = "10000000"
 dir_list ="./0000 ./0001 ./0002 ./0003 ./0004 ./0005 ./0006 ./0007 ./0008 ./0009 ./0010 ./0011 ./0012 ./0013 ./0014 ./0015 ./0016 ./0017 ./0018 ./0019 ./0020 ./0021 ./0022 ./0023 ./0024 ./0025 ./0026 ./0027 ./0028 ./0029 ./0030 ./0031 
./0032 ./0033 ./0034 ./0035 ./0036 ./0037 ./0038 ./0039 ./0040 ./0041  ./0042 ./0043 ./0044 ./0045 ./0046 ./0047  ./0048 ./0049 ./0050 ./0051 ./0052 ./0053  ./0054 ./0055 ./0056 ./0057 ./0058 ./0059  ./0060 ./0061 ./0062 ./0063 ./0064 
./0065 ./0066 ./0067 ./0068 ./0069 ./0070 ./0071 ./0072 ./0073 ./0074 ./0075 ./0076 ./0077  ./0078 ./0079   "
 nskip = 15
 nread = 15
 ndivide = "8"
 ntype = 2
 nswitch = 1
 lbin = 0.050d-10
 zbox = 150d-10
/

出力ファイルは以下の例のように1列目に溶質座標、2列目にtype1の自由エネルギー面、3列目にtype2の自由エネルギー面、4,5列目は標準偏差が出力される。

$ cat free.dat 
#     q              type1              type2              type1_sigma        type2_sigma
 -7.5000000000E-009  0.0000000000E+000  0.0000000000E+000  0.0000000000E+000  0.0000000000E+000
 -7.4950000000E-009 -2.0820471542E-010 -2.8209696892E-010  1.9919827051E-010  3.3805845880E-010
 -7.4900000000E-009 -2.3801176585E-010 -2.3670491143E-010  6.9396960310E-010  8.5913687822E-010
 -7.4850000000E-009 -2.7762470714E-010  3.6091948561E-010  1.0293821244E-009  1.5316345744E-009

初回のグランドカノニカルMDでは以下のような自由エネルギー面が得られる。

得られた自由エネルギー面のデータfree.datを溶質種毎に分割、Bias Potentialを適用する範囲を抽出する。これはつまり、溶質座標の範囲外にあるデータを削除したファイルを手作業で作成すると言うことである。 tapstep_fit.exeやgauss.fit.exeは一種類の溶質毎にしか実行できないため、チュートリアル内のsep.shなどを用いて1列目に溶質座標、2列目にfittingを行う溶質種の自由エネルギー面となるよう加工する。

$ sh sep.sh

次いで、Bias Potentialを適用し自由エネルギー面を平坦にする領域を決める。tapering領域を-3~3Åとした場合Bias Potential作成に用いる領域は-5~4.5Å(count_areaの左端)を範囲とするとよい。これらの範囲外にあるデータはファイルから削除する。 適用範囲端の自由エネルギーのピークはBias適用範囲外でBiasが速やかに減衰しない原因となりうるため、端のピークは無視しある程度マージンを確保したうえでfittingを行うことを推奨する。 自由エネルギーエネルギー面を加工の例を以下に示す。始点となるOFFの末端の自由エネルギーは0.00と設定しない場合、gauss関数による残差のfittingでアーティファクトの原因となる裾野が長い関数が採用されることがあるため例のように末端で自由エネルギーが0となるように加工する。

cat free_io.dat 
-4.5000000000E-010 0.0000000000E+000
-4.4000000000E-010 1.0748884251E-004
-4.3500000000E-010 1.0753750715E-004
.
.
.
4.0500000000E-010 1.2090785664E+001
4.1000000000E-010 1.2088002962E+001
4.1500000000E-010 1.2092397189E+001
4.2000000000E-010 1.2094084230E+001

データの加工終了後は、各溶質の自由エネルギー面のfittingをそれぞれ行いBiasを作成する。

$ ~/FreeFlex/tools/tapstep_fit.exe ./tapstep_cation.nml
$ ~/FreeFlex/gauss_fit.exe ./gauss_cation.nml
$ head -n 2 ./bias_na_tap.dat > ./bias_na.dat
$ tail -n +2 ./bias_na_gauss.dat >> ./bias_na.dat

$ ~/FreeFlex/tools/tapstep_fit.exe ./tapstep_anion.nml
$ ~/FreeFlex/gauss_fit.exe ./gauss_anion.nml
$ head -n 2 ./bias_io_tap.dat > ./bias_io.dat
$ tail -n +2 ./bias_io_gauss.dat >> ./bias_io.dat

作成されるbiasは以下のようになる。

$cat bias_na.dat
&biasdata
TAPERED_STEP 1 1 3  -1.4444799096E-019 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.6669429107E-19 -0.8575724076E-10  0.5489977544E-10  0
GAUSS_1D   1 1   3 -0.1726478574E-19  0.8329702944E-10  0.4764115363E-10  0
GAUSS_1D   1 1   3 -0.2106322823E-19  0.5362231666E-11  0.2120837581E-10  0
&end

$cat bias_io.dat 
&biasdata
TAPERED_STEP 1 1 3  -1.4444799096E-019 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.6669429107E-19 -0.8575724076E-10  0.5489977544E-10  0
GAUSS_1D   1 1   3 -0.1726478574E-19  0.8329702944E-10  0.4764115363E-10  0
GAUSS_1D   1 1   3 -0.2106322823E-19  0.5362231666E-11  0.2120837581E-10  0
&end

Biasの更新・中期

前節では更新初期の自由エネルギー差が大きく系の密度分布が目的構造から解離した状態でのBias更新を扱った。 本節では自由エネルギー差が1~2kcal/mol程度まで小さくなった状態を想定し、Biasの更新を行う。

作業する際はbias_update2にディレクトリ移動すること。

更新中期では得られた自由エネルギー面は前節、前々節と異なり、そのまま正負反転しfittingを行うのではなく、自由エネルギー面をscaleしたものを正負反転しfittingする。 これは、採取された自由エネルギー面は対となる溶質種との相互作用によって過剰に安定化・不安定化した状態のものであるためである。scaleせずにBiasを作成した場合凹凸の反転が発生しやすくなる。 過剰に安定化・不安定化を見積もることが原因であるためこれを解消するために0.5倍から0.8倍程度に自由エネルギー面のscaleを行う。

まず、input_4.nmlに従って自由エネルギー面の採取を行う。計算時間は緩和、サンプリングをまとめて400ps程度とる。計算を始める初期構造はmpiディレクトリの下にsystem4.msとして置いてあるためそれを利用すること。

$qsub pbssub.sh

input_4.nml

$cat input_4.nml 
&FreeFlex
  ms_init_file = "system_4.ms",
  xyz_traj_file = "traj_4.xyz",
  log_file = "log_4",
  ms_last_file = "system_5.ms",
  nml_output_file = "output_4.nml"
  ms_restart_file = "restart_4.ms"
  T = 298.15
  nstep = 200000
  dt = 2.0d-15
  omp_number = 1
  ewald_mode = "spme"
  pme_mode = .TRUE.
  calculate_real= .TRUE.
  calculate_recip = .TRUE.
  calculate_pol = .TRUE.
  calculate_damp = .FALSE.
  calculate_neut = .TRUE.
  calculate_mask = .TRUE.
  calculate_surf = .FALSE.
  interact_type = "gc"
  output_energy_interval = 100
  output_traj_interval = 2500
  output_restart_interval = 10000
  v_correct = .TRUE.
  n_restart_file = 5
  rdf_calc=.TRUE.
  rswitch = 0.0d0
  rcut_coulomb=11.87d-10
  pme_ngrid=32,32,192
  pme_kappa=0.3062d10
  induced_Ethresh =3.0d11
/

&addCOM
  condition = "HW OW"
/

&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .TRUE.
  COMid_solvent = 882
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 5.0d-10
  c_corr = 100
  sigma_corr = 0.5d-10
/

&gc_bias
  calculate_bias = .TRUE.
  nmax = 10
/



&gc_density
  calculate_density = .TRUE.
  output_filename = "gc_density.dat"
  output_interval = 2000
  lbin = 0.05d-10
/

&gc_free
  calculate_free = .TRUE.
  output_filename = "gc_free.dat"
  output_interval = 2000
  lbin = 0.05d-10
/

&accel_dequil
  do_accel_dequil = .TRUE.
  ntype = 2 
  count_area = 4.5d-10,6.0d-10
  force_area = -3.0d-10,5.0d-10
  nstep_dequil = 1000000
  set_conc = 1.2d0
  shift_v = 5d1
/

&gcbias
Na 1
t1 s1
TAPERED_STEP 1 1 3   8.5665306443E-019  7.0000000000E-011  1.9000000000E-010 0
GAUSS_1D   1 1   3  0.8643425869E-19 -0.1721676346E-10  0.6566984025E-10  0
GAUSS_1D   1 1   3  0.2774969855E-19  0.1681943494E-09  0.3212891372E-10  0
TAPERED_STEP 1 1 3  -1.4444799096E-019 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.6669429107E-19 -0.8575724076E-10  0.5489977544E-10  0
GAUSS_1D   1 1   3 -0.1726478574E-19  0.8329702944E-10  0.4764115363E-10  0
GAUSS_1D   1 1   3 -0.2106322823E-19  0.5362231666E-11  0.2120837581E-10  0
TAPERED_STEP 1 1 3   5.6367049694E-021 -3.5000000000E-011  1.4500000000E-010 0
GAUSS_1D   1 1   3  0.1746087643E-19  0.6905508635E-11  0.1965401420E-10  0
GAUSS_1D   1 1   3 -0.5531106888E-20  0.1723882708E-09  0.3694116545E-10  0
GAUSS_1D   1 1   3 -0.3679289242E-20 -0.1044577670E-09  0.4765584004E-10  0
GAUSS_1D   1 1   3 -0.6942386519E-20  0.1055004884E-09  0.1898341592E-10  0
GAUSS_1D   1 1   3  0.1246560465E-20 -0.1612260279E-09  0.1423549213E-10  0
GAUSS_1D   1 1   3  0.2096869336E-20  0.4704703871E-10  0.1538894498E-10  0
GAUSS_1D   1 1   3 -0.1054638856E-20  0.2407564803E-09  0.1411863110E-10  0
GAUSS_1D   1 1   3 -0.2008488555E-20 -0.2213342166E-09  0.2306952979E-10  0
TAPERED_STEP 1 1 3  -9.0513937394E-021 -3.5000000000E-011  1.4500000000E-010 0
GAUSS_1D   1 1   3  0.1867964227E-19  0.2844493693E-10  0.5633667847E-10  0
GAUSS_1D   1 1   3  0.1458957556E-19  0.9499252586E-10  0.2602991754E-10  0
GAUSS_1D   1 1   3 -0.2807485230E-20  0.2095724608E-09  0.5931251045E-10  0
TAPERED_STEP 1 1 3  -2.0161095525E-021 -3.5000000000E-011  1.4500000000E-010 0
GAUSS_1D   1 1   3  0.6640814517E-20  0.4889941131E-10  0.1318787347E-09  0
GAUSS_1D   1 1   3 -0.2681934535E-20  0.8700605321E-10  0.3940907930E-10  0
GAUSS_1D   1 1   3  0.3040795340E-20  0.1136563041E-09  0.1431926441E-10  0
GAUSS_1D   1 1   3 -0.1188108338E-20  0.1458639657E-10  0.1259017447E-10  0
GAUSS_1D   1 1   3  0.6936817359E-21 -0.7577943375E-10  0.1270915398E-10  0
GAUSS_1D   1 1   3 -0.1799509286E-20 -0.1354719190E-09  0.7065455500E-10  0
GAUSS_1D   1 1   3 -0.6490408596E-21  0.3266769284E-09  0.3437623905E-10  0
TAPERED_STEP 1 1 3   7.2737720846E-022 -1.0500000000E-010  1.4500000000E-010 0
GAUSS_1D   1 1   3 -0.3104946873E-20 -0.8028007779E-11  0.7150264513E-10  0
GAUSS_1D   1 1   3  0.3985709430E-20  0.9896383999E-10  0.3366050449E-10  0
GAUSS_1D   1 1   3  0.2235256436E-21  0.2471993044E-09  0.3488276262E-10  0
TAPERED_STEP 1 1 3   8.0290546345E-022 -1.0500000000E-010  1.4500000000E-010 0
GAUSS_1D   1 1   3  0.5544767035E-20  0.5512082241E-10  0.6120346181E-10  0
GAUSS_1D   1 1   3 -0.8277943433E-21  0.1690409452E-09  0.2109678370E-10  0
GAUSS_1D   1 1   3 -0.2703997221E-21 -0.1206809282E-09  0.3905443804E-10  0
GAUSS_1D   1 1   3 -0.2825926550E-21  0.2406317210E-09  0.3221819662E-10  0
TAPERED_STEP 1 1 3   1.2503126790E-021 -1.0500000000E-010  1.4500000000E-010 0
GAUSS_1D   1 1   3  0.3372855785E-20  0.1171527312E-09  0.3669376627E-10  0
GAUSS_1D   1 1   3 -0.8149085974E-21 -0.3902440824E-10  0.1530391074E-10  0
GAUSS_1D   1 1   3 -0.4001277698E-21  0.1762147429E-09  0.8534047577E-11  0
GAUSS_1D   1 1   3 -0.3655631251E-21  0.8385413308E-10  0.8835928664E-11  0
GAUSS_1D   1 1   3 -0.1162516238E-20  0.1546087623E-09  0.1133191277E-10  0
GAUSS_1D   1 1   3  0.5054645358E-21  0.2469386527E-09  0.5889168871E-10  0
GAUSS_1D   1 1   3 -0.6144094552E-21 -0.1152030682E-09  0.1241935710E-10  0
GAUSS_1D   1 1   3 -0.3193282321E-21 -0.1642222475E-09  0.1025797840E-10  0
GAUSS_1D   1 1   3  0.4681925782E-21  0.4024341889E-10  0.2311663018E-10  0
GAUSS_1D   1 1   3 -0.8945593293E-22 -0.3003724634E-09  0.3051087446E-10  0
GAUSS_1D   1 1   3 -0.2214580606E-21  0.2886091400E-09  0.1202662249E-10  0
GAUSS_1D   1 1   3 -0.3176615171E-21 -0.4191229025E-11  0.1123317252E-10  0
GAUSS_1D   1 1   3 -0.1228246496E-21 -0.1407703786E-09  0.9579335684E-11  0
GAUSS_1D   1 1   3 -0.6142466103E-21 -0.7590715324E-10  0.1461407018E-10  0
&end

&gcbias
Io 1
t2 s1
TAPERED_STEP 1 1 3   5.5975196375E-019  7.0000000000E-011  1.8000000000E-010 0
GAUSS_1D   1 1   3  0.5717509241E-19 -0.2643199901E-10  0.5732894427E-10  0
GAUSS_1D   1 1   3  0.2861966463E-19  0.1214287496E-09  0.6142124609E-10  0
GAUSS_1D   1 1   3 -0.7556229278E-20  0.2131926188E-09  0.3008222160E-10  0
TAPERED_STEP 1 1 3  -1.4444799096E-019 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.6669429107E-19 -0.8575724076E-10  0.5489977544E-10  0
GAUSS_1D   1 1   3 -0.1726478574E-19  0.8329702944E-10  0.4764115363E-10  0
GAUSS_1D   1 1   3 -0.2106322823E-19  0.5362231666E-11  0.2120837581E-10  0
TAPERED_STEP 1 1 3   9.6527894990E-021 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.1783590945E-19  0.6693420295E-10  0.4115871888E-10  0
GAUSS_1D   1 1   3 -0.7914073744E-20 -0.9992820655E-10  0.2995695109E-10  0
GAUSS_1D   1 1   3  0.2088972377E-19  0.5057548116E-11  0.2268126934E-10  0
GAUSS_1D   1 1   3  0.4178374621E-20  0.2087012452E-09  0.5493587958E-10  0
GAUSS_1D   1 1   3 -0.1885019211E-20 -0.2180675183E-09  0.2519511661E-10  0
TAPERED_STEP 1 1 3  -1.5372157299E-022 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.1752408350E-19  0.9621222313E-10  0.4578294564E-10  0
GAUSS_1D   1 1   3  0.1336596612E-19  0.9628978390E-11  0.3734069566E-10  0
GAUSS_1D   1 1   3 -0.9763620486E-21  0.2196621318E-09  0.4049786695E-10  0
TAPERED_STEP 1 1 3  -2.1093126368E-021 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.6381679068E-20 -0.1113437969E-10  0.6858410543E-10  0
GAUSS_1D   1 1   3  0.3012390267E-20  0.2079700873E-09  0.3994122373E-10  0
GAUSS_1D   1 1   3 -0.4475465732E-20  0.9283608262E-10  0.2923528631E-10  0
GAUSS_1D   1 1   3 -0.1627292168E-20  0.3354837852E-10  0.1595257691E-10  0
GAUSS_1D   1 1   3 -0.7333652889E-21 -0.6302201515E-10  0.1467032041E-10  0
GAUSS_1D   1 1   3  0.5281354769E-21  0.2953311713E-09  0.2269308957E-10  0
TAPERED_STEP 1 1 3   8.7839973552E-022 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3 -0.7068715200E-20  0.1043201921E-10  0.9445365973E-10  0
GAUSS_1D   1 1   3  0.3151319709E-20  0.1520062846E-09  0.6125479325E-10  0
GAUSS_1D   1 1   3  0.2591981037E-20  0.3274576916E-09  0.4142946617E-10  0
GAUSS_1D   1 1   3  0.1127481863E-20 -0.1068271723E-11  0.1202169924E-10  0
GAUSS_1D   1 1   3  0.7145393317E-21  0.2514314158E-09  0.2123354491E-10  0
TAPERED_STEP 1 1 3  -9.5389206780E-023 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3  0.3286452773E-20  0.3056853785E-09  0.6505478488E-10  0
GAUSS_1D   1 1   3  0.6201111547E-20  0.7523161444E-10  0.6303184745E-10  0
GAUSS_1D   1 1   3 -0.2991203328E-21 -0.9308836935E-10  0.4704561271E-10  0
GAUSS_1D   1 1   3  0.1045149122E-20  0.2159335116E-09  0.2887830983E-10  0
TAPERED_STEP 1 1 3   1.9578825718E-021 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3 -0.2535642578E-20 -0.1168692317E-10  0.6620591380E-10  0
GAUSS_1D   1 1   3  0.6117429793E-21  0.2430283357E-09  0.4486957722E-10  0
GAUSS_1D   1 1   3  0.7530810237E-21  0.2986520368E-10  0.1432540804E-10  0
GAUSS_1D   1 1   3  0.1478897486E-20  0.1581132726E-09  0.2810252773E-10  0
GAUSS_1D   1 1   3 -0.4184055158E-21 -0.1345513457E-09  0.2585004657E-10  0
GAUSS_1D   1 1   3  0.1394733825E-21 -0.2374826942E-09  0.2291198666E-10  0
&end

自由エネルギー面の計算は前節と同様にgc_analyzer.exeを実行すればよい。

$~/FreeFlex/tools/gc_analyzer.exe ./analyze_free.nml 

input_4.nmlを実行すると以下のような自由エネルギー面が得られる。

次に得られた自由エネルギー面を所望する倍率でscaleを行う。チュートリアルではscalefree.f90によって自由エネルギー面を元の0.6倍となるようデータの加工を行っている。 scaleの倍率を変更する際はscalefree.f90内のscale_paraの値を適宜変更する。 scaleした自由エネルギー面のデータファイルを前節と同様にsep.shで溶質種毎にファイルを分割する。

$ifort scalefree.f90 -o scalefree
$./scalefree
$sh sep.sh

その後、Bias適用の範囲外のデータを削除しfittingを行えるよう加工する。 ついでtapstep.fit.exe、gauss_fit.exeを実行し、scaleされた自由エネルギー面を反転させたものを解析関数でfittingする。

$ ~/FreeFlex/tools/tapstep_fit.exe ./tapstep_cation.nml
$ ~/FreeFlex/gauss_fit.exe ./gauss_cation.nml
$ head -n 2 ./bias_na_tap.dat > ./bias_na.dat
$ tail -n +2 ./bias_na_gauss.dat >> ./bias_na.dat

$ ~/FreeFlex/tools/tapstep_fit.exe ./tapstep_anion.nml
$ ~/FreeFlex/gauss_fit.exe ./gauss_anion.nml
$ head -n 2 ./bias_io_tap.dat > ./bias_io.dat
$ tail -n +2 ./bias_io_gauss.dat >> ./bias_io.dat

得られたBiasとそれをプロットした結果の一例を以下に示す。

$ cat bias_na.dat 
&biasdata
TAPERED_STEP 1 1 3  -7.1684784424E-022 -1.0500000000E-010  1.4500000000E-010 0
GAUSS_1D   1 1   3 -0.7033737394E-21 -0.1614014555E-10  0.2145777148E-10  0
GAUSS_1D   1 1   3  0.2803127863E-20  0.1068113514E-09  0.1349916449E-09  0
GAUSS_1D   1 1   3 -0.6466735617E-20  0.1152957934E-09  0.2039594985E-10  0
GAUSS_1D   1 1   3 -0.6444928768E-20  0.7674273103E-10  0.2198668193E-10  0
GAUSS_1D   1 1   3 -0.2900601376E-20  0.3489154148E-10  0.2073406799E-10  0
GAUSS_1D   1 1   3 -0.1605564529E-21  0.2558645614E-09  0.2526410878E-10  0
&end

$ cat bias_io.dat 
&biasdata
TAPERED_STEP 1 1 3   4.8057988661E-022 -1.2000000000E-010  1.3000000000E-010 0
GAUSS_1D   1 1   3 -0.6757111881E-20  0.8958268324E-10  0.7359367215E-10  0
GAUSS_1D   1 1   3  0.2411106512E-20 -0.6009709946E-10  0.7240444122E-10  0
GAUSS_1D   1 1   3 -0.1763642201E-20  0.2921879177E-09  0.4786810218E-10  0
&end

得られたBias Potentialを適用し、再度グランドカノニカルMDを実行しBiasの更新を繰り返していく。 更新を繰り返す際は、直前のMD計算の最後に出力されるsystem fileをinit fileに設定し計算を行うことで効率的に更新が進められる。 更新中期でのBias最適化はON-OFF間の自由エネルギー差が0.2kcal/mol以下となるまで、以下に示すような自由エネルギー面が得られるまで更新を行う。

Bias更新・後期

セルサイズ変更時ON-OFF間の0.2~0.3kcal/mol程度の自由エネルギー差が生じるため、小さな系でのBias最適化は0.2kcal/mol程度までで十分である。 表面過剰量の収束や界面の厚み探索を行う際には少なくとも自由エネルギー差を0.05kcal/mol以下にする必要がある。この0.2kcalから0.05kcal/mol以下まで自由エネルギー差を解消する過程を更新後期に位置付ける。

更新後期においても自由エネルギー差を解消する手順はこれまでと同一でありグランドカノニカルMDで採取した自由エネルギー面を反転させBiasを作成、Biasを適用しMDを行う作業を繰り返す。 自由エネルギー差が小さくなるためエラーバーが自由エネルギー差以下となるようにするため計算時間を長く取る必要があることに留意すること。