チュートリアル10:溶液界面のグランドカノニカルMD

溶液の界面の系を適用対象として液相のバルクにおける溶質密度を一定とするグランドカノニカルMDの実行方法を解説する。理論的な説明はこちらを参照。

目的

今回実行するシミュレーションの目的は、バルク溶質密度をgivenとする条件下における界面に垂直な一次元密度分布の計算である。粒子数を一定とする溶液界面のMD計算では、スラブ厚さが有限であることにより界面の密度構造に依存してバルク部分の溶質密度が変動する。これに対して今回のシミュレーション手法は、スラブ中央部分における溶質密度を固定することで望みの溶質濃度・望みの界面厚さを仮定したときの界面密度構造を求めることを可能とする。

計算の流れ

本手法では密度分布のサンプリングに先立ち、溶質分子のON-OFF切り替えの自由エネルギー差を打ち消すバイアスポテンシャルを決定する必要がある。計算の流れは以下の通りである。

1. 粒子数一定のMD計算による初期バイアスポテンシャルの生成。

2. グランドカノニカルMDによるバイアスポテンシャルの最適化。

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

バイアスポテンシャルは、自由エネルギー面の全体に-1を掛けたものを解析関数の和でfittingすることで求められる。自由エネルギー面の計算手法として、溶質座標軸上の各点各点における平均力を数値積分する熱力学的積分法を採用する。 1.では1つの溶質分子にON-OFFスイッチングを適用してその溶質座標を拘束し、拘束位置の異なるトラジェクトリを並列計算して各点上の平均力を求める。 2.では全ての溶質分子にON-OFFスイッチングを適用し、溶質座標をビンで区切ることで平均力のヒストグラムを求める。

バイアスポテンシャルの最適化では、バイアスを適用したMD計算により自由エネルギー面を求め、それをfittingして得られた新たなバイアスをもとのバイアスに加える、という操作を繰り返す。最終的な自由エネルギー面の平坦さの精度は、密度分布に系統的な誤差を与える。

使用プログラム

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

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

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

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

初期バイアスポテンシャルの生成

以下、グランドカノニカルMDの最も単純な適用例として、希ガス混合流体の気液界面の計算を示す。溶媒は液体アルゴン、溶質はキセノンとする。

はじめに初期構造生成用スクリプトを実行する（gitで公開されているチュートリアル用スクリプトは各種プログラムのパス部分のみ修正が必要）。実行すると指定のMPI並列数と同数のディレクトリ0000/, 0001/, ...が生成され、その下には初期構造のmsファイル1つとFreeFlex.exeの入力ネームリストファイル2つが置かれる。

2つのネームリストファイルはそれぞれ平衡化用とサンプリング用であり、内容は同一で、拘束位置が記入されたネームリスト&addconstrのみが記載されている。以下、スクリプトを実行したディレクトリをワーキングディレクトリ、生成された0000/などのディレクトリをMPIディレクトリと呼ぶ。

$ ls
wd1/ wd2/
$ cd ./wd1/
$ vi ./gc_gen_tdi_input.sh
$ sh ./gc_gen_tdi_input.sh
$ ls ./0000/
input_1.nml input_2.nml system1.ms

$ cat ./0000/input_1.nml
&addconstr
condition = "SPCOORD 1 -6.35d-10  1.0d-10 0.01d-10"
/

次にワーキングディレクトリにFreeFlex.exe入力ネームリストファイルを作成し、粒子数一定のMD計算を実行する。現在は下のコマンドをそのまま実行するのではなく、PBSproを利用することが推奨される。

拘束される溶質分子は初期構造で拘束位置付近に置かれるため、平衡化の最初から拘束をかけることができる。拘束位置はMPIディレクトリのネームリストファイルから読み込まれる。グランドカノニカルMD用spcoordは自動生成されるが、拘束用およびmonitor出力用のspcoordは別途明示的に定義する必要がある。平衡化でははじめ気液界面に壁を設置して分子が気相に弾き出されることを防ぎ、サンプリング計算では拘束された溶質分子が感じる力をmonitor出力する。希ガス混合流体の計算ではcalculate_coulomb=.FALSE.を指定することで静電相互作用計算を省略し計算の高速化が見込める。

$ mpiexec –machinefile machinefile –n 128 ~/FreeFlex/FreeFlex.exe ./input_1.nml
$ mpiexec –machinefile machinefile –n 128 ~/FreeFlex/FreeFlex.exe ./input_2.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 = 83.79
  nstep = 100000
  dt = 1.0d-15
  omp_number = 1
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "gc"
  output_energy_interval = 1000
  output_traj_interval = 1000
  output_restart_interval = 10000
  set_wall = .TRUE.
  z_wall = "-55.13d-10 55.13d-10"
  nstep_wall = 50000
  v_correct = .TRUE.
/
&addCOM
  condition = "Ar"
/
&spcoord
  condition = "Z 5431 5432"
/
&constraint2
  tolerance = 1.0d-16
  maxiter = 200
  threshold = 0.2d0
/
&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .FALSE.
  COMid_solvent = 5432
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  tapering_delta = "3.0d-10"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 3.0d-10
  c_corr = 100 
  sigma_corr = 5.0d-11
/

$ 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 = 83.79
  nstep = 500000
  dt = 1.0d-15
  omp_number = 1
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "gc"
  output_energy_interval = 1000
  output_traj_interval = 1000
  output_restart_interval = 10000
  set_wall = .FALSE.
  v_correct = .FALSE.
/
&addCOM
  condition = "Ar"
/
&spcoord
  condition = "Z 5431 5432"
/
&constraint2
  tolerance = 1.0d-16
  maxiter = 200
  threshold = 0.2d0
/
&MONITOR
  interval = 1
  file = "spcoord_f_2.dat"
  condition = "SPCOORD_F 1"
  separator = "$"
/
&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .FALSE.
  COMid_solvent = 5432
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 3.0d-10
  c_corr = 100 
  sigma_corr = 5.0d-11
/

計算終了後、ワーキングディレクトリで解析プログラムtherm_int.exeを実行する。 MPIディレクトリのmonitor出力ファイルを読み込んで各拘束位置における平均力を計算し、熱力学的積分法を適用して一次元自由エネルギー面を求める。

$ ~/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 = "1000000"
  mpi_number = 128
  start = -6.35d-10
  interval = 0.1d-10
/

出力されるキセノンの自由エネルギー面の例を以下に示す。横軸は溶質座標で、z=0が仮想的な半透膜、右側がON、左側がOFFであり、縦軸はOFFを基準とした自由エネルギーである。 tapering関数の値が変化する範囲で自由エネルギーも大きく変化し、溶媒和を感じるONが安定化されていることが分かる。

得られた自由エネルギー面を解析関数でfittingし、グランドカノニカルMDの初期バイアスポテンシャルを決定する。第一にON-OFF間の自由エネルギー差を相殺するため、 tapstep_fit.exeを実行して自由エネルギー面を1つのTAPERED_STEPバイアスでfittingする。第二に半透膜近傍に残った凹凸を打ち消すため、 gauss_fit.exeを実行してfitting残差を任意数のGAUSS_1Dバイアスの和でfittingする。

$ ~/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 = "bias1.dat"
  plotfile = "plot1.dat"
  restfile = "gfinput.dat"
  center = 0.0d-10
  delta = 2.2d-10
/

$ cat ./gauss_fit.nml
&gauss_fit
  input_file_name = "gfinput.dat"
  param_file_name = ""
  bias_file_name = "bias2.dat"
  plot_file_name = "plot2.dat"
  log_file_name = ""
  spcoord_ID = ""
  ngauss_max = 50
  niter_max = 1000
  thres_finish = 1.5d-1
  thres_unchange = 1.0d-1
  lambda_1 = 1.0d-3
  lambda_2 = 1.0d-4
  lambda_3 = 1.0d-15
  ncolumn_data = 1
  ncolumn_weight = 0
  omp_number = 1
  plot_result = .TRUE.
/

TAPERED_STEPの中心centerと幅delta, gauss_fit.exeのパラメータであるthresholdやlambdaにはそれほど拘らなくてもfittingの成否に支障はないと思われるが、少ない数のバイアスできれいにfittingしたい場合は試行錯誤的な調整が必要になることがある。

Gauss関数によるfittingでは、入力データにおける離散点の存在領域より外側は速やかに一定値に収束することが望ましい。しかしgauss_fit.exeは離散点上の残差のみを計算し、より外側のバイアス形状については関知しない。離散点の存在領域よりも外側にピークのあるバイアスを採用した場合、その後のバイアス更新では修正の難しいアーティファクトとなる危険性がある。これを防ぐため、biasplot.exeを用いてより外側のバイアス形状を確認しておくとよい（使用方法はソースコード冒頭のコメントを参照）。離散点の存在領域外のバイアス形状に問題が見られた場合、gauss_fit.exeの入力パラメータ等を試行錯誤的に調整する。

最後にtapstep_fit.exeとgauss_fit.exeの結果を合わせ、グランドカノニカルMDの初期バイアスポテンシャルを得る。使用する際は&gcbiasの仕様に合わせて適宜修正してからFreeFlex.exe入力ネームリストファイルに貼り付ける。

$ head -n 2 ./bias1.dat > ./bias3.dat
$ tail -n +2 ./bias2.dat >> ./bias3.dat

$ cat ./bias3.dat
&biasdata
TAPERED_STEP 1 1 3  1.3160346281E-020  0.0000000000E+000  2.2000000000E-010 0
GAUSS_1D 1 1 3  0.3882291961E-21 -0.3922052640E-10  0.3000026325E-10 0
GAUSS_1D 1 1 3 -0.2180851492E-21  0.5253120809E-10  0.3358765741E-10 0
GAUSS_1D 1 1 3 -0.1803022811E-21 -0.1019720958E-09  0.5803369593E-10 0
GAUSS_1D 1 1 3  0.1272394604E-21  0.1996464698E-09  0.1593226247E-09 0
GAUSS_1D 1 1 3 -0.2005171690E-21  0.1746743029E-09  0.4626665197E-10 0
GAUSS_1D 1 1 3  0.8201276248E-22 -0.1748582373E-09  0.3123850723E-10 0
&end

バイアスポテンシャルの更新

グランドカノニカルMDを実行して自由エネルギー面を求め、バイアスポテンシャルの改良を行う。今回はサンプル数を増やすためにMPI並列計算として実行する例を示す。

初期平衡化は2段階で行う。第1段階として溶媒とONだけの系の初期構造を生成し、粒子数一定のMD計算で緩和させる。初期構造の生成には専用のスクリプト1を用いる（gitで公開されているチュートリアル用スクリプトは各種プログラムのパス部分のみ修正が必要）。ただし、このスクリプトはMPI並列計算用のディレクトリを生成しないため、必要であれば事前に作成しておく。

$ cd ../wd2/
$ for dir in `seq -f %04g 0 15`
> do
> mkdir $dir
> done
$ vi gc_gen_system_1.sh
$ sh gc_gen_system_1.sh
$ ls ./0000/
system1.ms
$ mpiexec –machinefile machinefile –n 16 ~/FreeFlex/FreeFlex.exe ./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_on.ms",
  nml_output_file = "output_1.nml"
  ms_restart_file = "restart_1.ms"
  T = 83.79
  nstep = 10000
  dt = 1.0d-15
  omp_number = 1
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "default"
  output_energy_interval = 100
  output_traj_interval = 100
  output_restart_interval = 1000
  set_wall = .TRUE.
  z_wall = "-55.13d-10 55.13d-10"
  nstep_wall = 100000
  v_correct = .TRUE.
/
&gc
  perform_gc = .FALSE.
/

第2段階として緩和した溶媒とONの系にOFFとDummyを組み合わせて全系の初期構造を生成し、グランドカノニカルMDによる平衡化を行う。初期構造の生成には専用のスクリプト2を用いる（gitで公開されているチュートリアル用スクリプトは各種プログラムのパス部分のみ修正が必要）。前節で作成した初期バイアスポテンシャルを全ての溶質分子に適用する。平均力サンプリングではネームリスト&gc_freeを用いてデータの出力を行う。

$ vi ./gc_gen_system_2.sh
$ sh ./gc_gen_system_2.sh
$ ls ./0000/
... system1.ms system2.ms system2_on.ms ...
$ mpiexec –machinefile machinefile –n 16 ~/FreeFlex/FreeFlex.exe ./input_2.nml
$ mpiexec –machinefile machinefile –n 16 ~/FreeFlex/FreeFlex.exe ./input_3.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 = 83.79
  nstep = 200000
  dt = 1.0d-15
  omp_number = 4
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "gc"
  output_energy_interval = 1000
  output_traj_interval = 1000
  output_restart_interval = 10000
  set_wall = .TRUE.
  z_wall = "-55.13d-10 55.13d-10"
  nstep_wall = 100000
  v_correct = .FALSE.
/
&addCOM
  condition = "Ar"
/
&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .TRUE.
  COMid_solvent = 6174
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 3.0d-10
  c_corr = 100 
  sigma_corr = 5.0d-11
/
&gc_bias
  calculate_bias = .TRUE.
  nmax = 10
/
&gcbias
comment 1
t1 s1
TAPERED_STEP 1 1 3  1.3160346281E-020  0.0000000000E+000  2.2000000000E-010 0
GAUSS_1D 1 1 3  0.3882291961E-21 -0.3922052640E-10  0.3000026325E-10 0
GAUSS_1D 1 1 3 -0.2180851492E-21  0.5253120809E-10  0.3358765741E-10 0
GAUSS_1D 1 1 3 -0.1803022811E-21 -0.1019720958E-09  0.5803369593E-10 0
GAUSS_1D 1 1 3  0.1272394604E-21  0.1996464698E-09  0.1593226247E-09 0
GAUSS_1D 1 1 3 -0.2005171690E-21  0.1746743029E-09  0.4626665197E-10 0
GAUSS_1D 1 1 3  0.8201276248E-22 -0.1748582373E-09  0.3123850723E-10 0
&end

$ cat ./input_3.nml
&FreeFlex
  ms_init_file = "system3.ms"
  xyz_traj_file = "traj_3.xyz",
  log_file = "log_3",
  ms_last_file = "system4.ms",
  nml_output_file = "output_3.nml"
  ms_restart_file = "restart_3.ms"
  T = 83.79
  nstep = 3000000
  dt = 1.0d-15
  omp_number = 4
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "gc"
  output_energy_interval = 1000
  output_traj_interval = 1000
  output_restart_interval = 10000
  set_wall = .FALSE.
  v_correct = .FALSE.
/
&addCOM
  condition = "Ar"
/
&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .TRUE.
  COMid_solvent = 6174
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 3.0d-10
  c_corr = 100 
  sigma_corr = 5.0d-11
/
&gc_bias
  calculate_bias = .TRUE.
  nmax = 10
/
&gc_free
  calculate_free = .TRUE.
  is_formatted = .TRUE.
  output_filename = "gc_free_3.dat"
  output_interval = 10000
  lbin = 0.2d-10
/
&gcbias
comment 1
t1 s1
TAPERED_STEP 1 1 3  1.3160346281E-020  0.0000000000E+000  2.2000000000E-010 0
GAUSS_1D 1 1 3  0.3882291961E-21 -0.3922052640E-10  0.3000026325E-10 0
GAUSS_1D 1 1 3 -0.2180851492E-21  0.5253120809E-10  0.3358765741E-10 0
GAUSS_1D 1 1 3 -0.1803022811E-21 -0.1019720958E-09  0.5803369593E-10 0
GAUSS_1D 1 1 3  0.1272394604E-21  0.1996464698E-09  0.1593226247E-09 0
GAUSS_1D 1 1 3 -0.2005171690E-21  0.1746743029E-09  0.4626665197E-10 0
GAUSS_1D 1 1 3  0.8201276248E-22 -0.1748582373E-09  0.3123850723E-10 0
&end

計算終了後、gc_analyzer.exeを実行して一次元自由エネルギー面を求める。 MPI並列化された全トラジェクトリの出力データを総合して1つの自由エネルギー面を計算する。なお、gc_analyzer.exeは半透膜の位置q=0を前提として解析を行うため、他の座標を指定してグランドカノニカルMDを実行した場合は正しい自由エネルギー面を得ることができない。

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

$ cat ./analyze_1.nml
&free
  is_formatted = .TRUE.
  inputfile = "gc_free_3.dat"
  outputfile = "free3.dat"
  max_nread = "300"
  dir_list = "./0000 ./0001 ./0002 ./0003 ./0004 ./0005 ./0006 ./0007 ./0008 ./0009 ./0010 ./0011 ./0012 ./0013 ./0014 ./0015"
  ndivide = "4"
  nskip = 0
  nread = 300
  ntype = 1
  nswitch = 1
  lbin = 0.2d-10
  zbox = 300.0d-10
/

出力ファイルの1列目が溶質座標、2列目が溶質の自由エネルギー面、3列目が溶質の標準誤差（エラーバー）である。出力自由エネルギー面は溶質座標の全範囲に渡るため、バイアスポテンシャルの決定に先立ち半透膜の近傍を抜き出す必要がある。以下に自由エネルギー面の例を示す。前節の結果と比べてON-OFF境界近傍がフラットになることが期待される。

続いてtapstep_fit.exeおよびgauss_fit.exeを実行し、バイアスポテンシャルを生成する。得られたバイアスを上で使用した初期バイアスに追記することでバイアスポテンシャルの改良が完了する。

$ ~/FreeFlex/tools/tapstep_fit.exe ./tapstep_fit.nml
$ ~/FreeFlex/gauss_fit.exe ./gauss_fit.nml
$ head -n 2 ./bias1.dat > ./bias3.dat
$ tail -n +2 ./bias2.dat >> ./bias3.dat

$ cat ./tapstep_fit.nml
&tapstepfit
  inputfile = "free3.dat"
  biasfile = "bias1.dat"
  plotfile = "plot1.dat"
  restfile = "gfinput.dat"
  center = 2.0d-10
  delta = 0.8d-10
/

$ cat ./gauss_fit.nml
&gauss_fit
  input_file_name = "gfinput.dat"
  param_file_name = ""
  bias_file_name = "bias2.dat"
  plot_file_name = "plot2.dat"
  log_file_name = ""
  spcoord_ID = ""
  ngauss_max = 50
  niter_max = 1000
  thres_finish = 1.5d-1
  thres_unchange = 1.0d-1
  lambda_1 = 1.0d-3
  lambda_2 = 1.0d-4
  lambda_3 = 1.0d-15
  ncolumn_data = 1
  ncolumn_weight = 0
  omp_number = 1
  plot_result = .TRUE.
/

$ cat ./bias3.dat
&biasdata
TAPERED_STEP 1 1 3 -8.4022715727E-023  2.0000000000E-010  8.0000000000E-011 0
GAUSS_1D 1 1 3 -0.5362199787E-22  0.4181974358E-09  0.5879943512E-10 0
GAUSS_1D 1 1 3 -0.1479540976E-20  0.1992593700E-09  0.8023420425E-10 0
GAUSS_1D 1 1 3  0.5886826408E-22 -0.2233793837E-10  0.2164258711E-10 0
GAUSS_1D 1 1 3  0.1489705494E-20  0.1985126857E-09  0.7506830971E-10 0
GAUSS_1D 1 1 3  0.1476489327E-22 -0.8139569836E-10  0.4428985243E-10 0
&end

密度分布の計算

最適化されたバイアスポテンシャルを用いてグランドカノニカルMDを実行し、界面の密度分布を計算する。平均力のサンプリングも同時に行い、自由エネルギー面が十分フラットであることを確認する。実用上は、前節で示した1回目のバイアス改良作業の後すぐに本節の計算を実行し、自由エネルギー面の凹凸が大きい場合はバイアスを更新して本節の計算を再実行、十分に平坦であればバイアスの最適化を終了してそのときの密度分布を採用する。

前節のMD計算の最終的な構造を初期構造として、グランドカノニカルMDによる平衡化とサンプリングを実行する。平衡化はバイアスポテンシャルの変更に伴うものである。サンプリング計算ではネームリスト&gc_freeと&gc_densityを同時に指定して平均力と密度のサンプリングを行う。なお、チュートリアル用の入力ネームリストには十分にフラットな自由エネルギー面を与えるバイアスを記載した。

$ mpiexec –machinefile machinefile –n 16 ~/FreeFlex/FreeFlex.exe ./input_4.nml
$ mpiexec –machinefile machinefile –n 16 ~/FreeFlex/FreeFlex.exe ./input_5.nml

$ cat ./input_4.nml
&FreeFlex
  ms_init_file = "system4.ms"
  xyz_traj_file = "traj_4.xyz",
  log_file = "log_4",
  ms_last_file = "system5.ms",
  nml_output_file = "output_4.nml"
  ms_restart_file = "restart_4.ms"
  T = 83.79
  nstep = 100000
  dt = 1.0d-15
  omp_number = 4
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "gc"
  output_energy_interval = 1000
  output_traj_interval = 1000
  output_restart_interval = 10000
  set_wall = .FALSE.
  v_correct = .FALSE.
/
&addCOM
  condition = "Ar"
/
&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .TRUE.
  COMid_solvent = 6174
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 3.0d-10
  c_corr = 100 
  sigma_corr = 5.0d-11
/
&gc_bias
  calculate_bias = .TRUE.
  nmax = 10
/
&gcbias
comment 1
t1 s1
TAPERED_STEP 1 1 3  1.3160346281E-020  0.0000000000E+000  2.2000000000E-010 0
GAUSS_1D 1 1 3  0.3882291961E-21 -0.3922052640E-10  0.3000026325E-10 0
GAUSS_1D 1 1 3 -0.2180851492E-21  0.5253120809E-10  0.3358765741E-10 0
GAUSS_1D 1 1 3 -0.1803022811E-21 -0.1019720958E-09  0.5803369593E-10 0
GAUSS_1D 1 1 3  0.1272394604E-21  0.1996464698E-09  0.1593226247E-09 0
GAUSS_1D 1 1 3 -0.2005171690E-21  0.1746743029E-09  0.4626665197E-10 0
GAUSS_1D 1 1 3  0.8201276248E-22 -0.1748582373E-09  0.3123850723E-10 0
TAPERED_STEP 1 1 3 -8.4022715727E-023  2.0000000000E-010  8.0000000000E-011 0
GAUSS_1D 1 1 3 -0.5362199787E-22  0.4181974358E-09  0.5879943512E-10 0
GAUSS_1D 1 1 3 -0.1479540976E-20  0.1992593700E-09  0.8023420425E-10 0
GAUSS_1D 1 1 3  0.5886826408E-22 -0.2233793837E-10  0.2164258711E-10 0
GAUSS_1D 1 1 3  0.1489705494E-20  0.1985126857E-09  0.7506830971E-10 0
GAUSS_1D 1 1 3  0.1476489327E-22 -0.8139569836E-10  0.4428985243E-10 0
TAPERED_STEP 1 1 3 -6.7351505906E-023  5.0000000000E-011  2.5000000000E-010 0
GAUSS_1D 1 1 3  0.7981367928E-23  0.2891055373E-09  0.7274529369E-10 0
GAUSS_1D 1 1 3  0.9202499884E-23  0.3056362429E-09  0.2259943582E-10 0
GAUSS_1D 1 1 3  0.9669170145E-23  0.5569579921E-10  0.2001062402E-10 0
GAUSS_1D 1 1 3  0.1040307521E-22  0.1399215532E-09  0.1977047796E-10 0
&end

$ cat ./input_5.nml
&FreeFlex
  ms_init_file = "system5.ms"
  xyz_traj_file = "traj_5.xyz",
  log_file = "log_5",
  ms_last_file = "system6.ms",
  nml_output_file = "output_5.nml"
  ms_restart_file = "restart_5.ms"
  T = 83.79
  nstep = 3000000
  dt = 1.0d-15
  omp_number = 4
  rcut = 11.5d-10
  calculate_coulomb = .FALSE.
  interact_type = "gc"
  output_energy_interval = 1000
  output_traj_interval = 1000
  output_restart_interval = 10000
  set_wall = .FALSE.
  v_correct = .FALSE.
/
&addCOM
  condition = "Ar"
/
&gc
  perform_gc = .TRUE.
  do_off_dummy_switch = .TRUE.
  COMid_solvent = 6174
  spcoord_type = "GC_Z_DQDX_CONST"
  type_tapering = "6"
  tapering_delta = "3.0d-10"
  delta_tapering = "3.0d-10"
  qswitch = "0.0d0"
  coulomb_sigma = 3.0d-10
  c_corr = 100 
  sigma_corr = 5.0d-11
/
&gc_bias
  calculate_bias = .TRUE.
  nmax = 10
/
&gc_density
  calculate_density = .TRUE.
  is_formatted = .TRUE.
  output_filename = "gc_density_5.dat"
  output_interval = 10000
  lbin = 0.5d-10
  COMid_append = ""
/
&gc_free
  calculate_free = .TRUE.
  is_formatted = .TRUE.
  output_filename = "gc_free_5.dat"
  output_interval = 10000
  lbin = 0.2d-10
/
&gcbias
comment 1
t1 s1
TAPERED_STEP 1 1 3  1.3160346281E-020  0.0000000000E+000  2.2000000000E-010 0
GAUSS_1D 1 1 3  0.3882291961E-21 -0.3922052640E-10  0.3000026325E-10 0
GAUSS_1D 1 1 3 -0.2180851492E-21  0.5253120809E-10  0.3358765741E-10 0
GAUSS_1D 1 1 3 -0.1803022811E-21 -0.1019720958E-09  0.5803369593E-10 0
GAUSS_1D 1 1 3  0.1272394604E-21  0.1996464698E-09  0.1593226247E-09 0
GAUSS_1D 1 1 3 -0.2005171690E-21  0.1746743029E-09  0.4626665197E-10 0
GAUSS_1D 1 1 3  0.8201276248E-22 -0.1748582373E-09  0.3123850723E-10 0
TAPERED_STEP 1 1 3 -8.4022715727E-023  2.0000000000E-010  8.0000000000E-011 0
GAUSS_1D 1 1 3 -0.5362199787E-22  0.4181974358E-09  0.5879943512E-10 0
GAUSS_1D 1 1 3 -0.1479540976E-20  0.1992593700E-09  0.8023420425E-10 0
GAUSS_1D 1 1 3  0.5886826408E-22 -0.2233793837E-10  0.2164258711E-10 0
GAUSS_1D 1 1 3  0.1489705494E-20  0.1985126857E-09  0.7506830971E-10 0
GAUSS_1D 1 1 3  0.1476489327E-22 -0.8139569836E-10  0.4428985243E-10 0
TAPERED_STEP 1 1 3 -6.7351505906E-023  5.0000000000E-011  2.5000000000E-010 0
GAUSS_1D 1 1 3  0.7981367928E-23  0.2891055373E-09  0.7274529369E-10 0
GAUSS_1D 1 1 3  0.9202499884E-23  0.3056362429E-09  0.2259943582E-10 0
GAUSS_1D 1 1 3  0.9669170145E-23  0.5569579921E-10  0.2001062402E-10 0
GAUSS_1D 1 1 3  0.1040307521E-22  0.1399215532E-09  0.1977047796E-10 0
&end

計算終了後、gc_analyzer.exeを実行して自由エネルギー面と密度分布を得る。ネームリスト&density, &freeを同時に与えると一度に両方の結果を出力する。

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

$ cat analyze_2.nml
&density
  is_formatted = .TRUE.
  inputfile = "gc_density_5.dat"
  outputfile = "dens5.dat"
  max_nread = "300"
  dir_list = "./0000 ./0001 ./0002 ./0003 ./0004 ./0005 ./0006 ./0007 ./0008 ./0009 ./0010 ./0011 ./0012 ./0013 ./0014 ./0015"
  ndivide = "4"
  nskip = 0
  nread = 300
  ntype = 1
  nswitch = 1
  lbin = 0.5d-10
  zbox = 300.0d-10
  xbox = 50.0d-10
  ybox = 50.0d-10
/
&free
  is_formatted = .TRUE.
  inputfile = "gc_free_5.dat"
  outputfile = "free5.dat"
  max_nread = "300"
  dir_list = "./0000 ./0001 ./0002 ./0003 ./0004 ./0005 ./0006 ./0007 ./0008 ./0009 ./0010 ./0011 ./0012 ./0013 ./0014 ./0015"
  ndivide = "4"
  nskip = 0
  nread = 300
  ntype = 1
  nswitch = 1
  lbin = 0.2d-10
  zbox = 300.0d-10
/

自由エネルギー面の例を以下に示す（半透膜近傍の抜粋）。出力ファイルの1列目が溶質座標、2列目が溶質の自由エネルギー面、3列目が溶質の標準誤差（エラーバー）である。この結果はON-OFF境界近傍において十分にフラットと判断できる。（バイアスの最適化が不十分とみられる場合、前節までと同様にtapstep_fit.exeとgauss_fit.exeを実行してバイアスを更新し、再度平衡化とサンプリングを行う）

続いて、溶媒および溶質の密度分布の例を示す。出力ファイルの1列目が溶質座標、2列目と3列目がそれぞれ溶質と溶媒の密度分布、 4列目と5列目がそれぞれ溶質と溶媒の標準誤差（エラーバー）である。溶質は溶質座標の正負を合算し、半透膜の位置をz=0としてz>0にON, z<0にOFFがプロットされる。溶媒は溶質座標原点を基準とした相対座標でサンプリングし、やはり正負を合算してz>0にプロットされる。溶質密度分布がON-OFF境界近傍でフラットであることを確認する。（不連続に見える等の場合はサンプリングの継続を検討する）

密度分布の計算結果を発表する際には、溶質と溶媒をそれぞれバルク部分における密度で規格化して1枚のグラフに描画するとよい。こうして作成したグラフの例を以下に示す。ここではz<0における溶媒の密度を表示しないため、0.0をNaNに置換してからプロットしている。

複数分子種・複数半透膜への拡張

ここまで希ガス混合流体の気液界面系を例にグランドカノニカルMDの実行手順を説明してきた。この系は溶質分子種が1種類のみ、かつ仮想的な半透膜が1枚だけの場合である。一方、電解質水溶液のグランドカノニカルMDでは複数分子種の取り扱いが必要であり、また液液界面の計算では複数半透膜を設置しなければならない。

そこで、複数分子種・複数半透膜のグランドカノニカルMDにおける主な留意点を以下に述べる。入力ネームリスト等の具体的な設定方法は各プログラムの解説ページやソースファイル冒頭のコメントを参照。

粒子数一定のMD計算
- （溶質分子種の数）×（半透膜の数）個のワーキングディレクトリを用意し、それぞれ着目する分子種と半透膜を設定する。計算終了後はワーキングディレクトリごとに自由エネルギー面を求めてfittingを行い、着目している分子種・半透膜を対象とした初期バイアスを得る。
グランドカノニカルMD
- 1本のトラジェクトリ（またはひとまとまりのMPI並列計算）で全分子種・全半透膜の自由エネルギー面や密度分布を計算できる。入力ファイルには（溶質分子種の数）×（半透膜の数）個の&gcbiasセクションを記述する。
- gc_analyzer.exeは半透膜3枚以上をサポートしない。半透膜2枚のときはその位置q=0, Lz/2（または-Lz/2）を前提とする。
- gc_analyzer.exeの出力ファイルでは各分子種がそれぞれ異なる列に書き出される。状態ごとに溶質座標の正負を合算して平均値を計算し、OFF1, ON, OFF2の順で一続きにプロットする。バイアスポテンシャルの生成では、出力自由エネルギー面から各半透膜の近傍領域を切り出し、それぞれ適宜平行移動してからfittingする必要がある。
- &gc_densityのネームリスト変数COMid_appendは、半透膜2枚の場合に溶質座標原点の定義に関わらない溶媒（油水界面系で溶質座標を水相の重心として定義したときの油相）の密度分布をサンプリングするためのオプションである。半透膜1枚のとき空文字列とし、2枚のときサンプリングしたい溶媒全分子の重心として別途定義した重心番号を1つだけ指定する。gc_analyzer.exeではこのように指定したことを前提としてファイルの読み込みと解析を行う。