チュートリアル03:water finger の解析

　本章ではFreeFlexによるwater finger の解析について説明する。系としては水-DCM系における Na⁺ の輸送を取り扱う。注意：このチュートリアルは古い &bias を用いており、現行FreeFlexに対応していない。water fingerの計算は、新しい章チュートリアル08:water finger の解析と2次元自由エネルギー面の計算を参照のこと。

初期構造の作成

　water finger を解析するための初期構造を作るには、イオンが界面油相側に存在する方が扱いやすい。そこでまず、Na⁺を水重心から 20 Å 離れた位置に置き、その周りにバイアス・ポテンシャル

$U^{bias}={\frac {k_{B}T}{2\sigma ^{2}}}(z-z^{0})^{2}$

を置いた状態で平衡化をする。

Na⁺ の追加はadd.exe を下記のように用いることで容易に実現できる。

=== add_Na.nml =================================================================================
&system
  output = "system.ms"
  size = 25.0 25.0 85.0
  Tinit = 298.15
/

# z_center =   1.213686084985836E-009
&add
  file="water_DCM.ms"
  n=1
  area="fix 0.0 0.0 -12.13686084985836"
  velocity ="maxwell"
  direction ="hold"
/

&add
  file="Na.ms"
  n=1
  area="fix 0.0 0.0 20.0"
  velocity ="fix 0.0 0.0 0.0"
/
================================================================================================

（water_DCM.ms は平衡化済みの水-DCM系、Na.ms は Na^+ イオンのmsファイル）

> ./add.exe add_Na.nml system.ms
> ./unlap.exe system.ms system2.ms

一方、バイアス・ポテンシャルを用いたMD計算の実行は、FreeFlex の入力ファイルに

=== input3.nml =================================================================================
  .
  .
  .
&bias
  calculate_bias = T,
  bias_type = "WF2D"
  bias_param_file = "bias_param3.dat"
  bias_output_file = "bias_data3.dat"
/
================================================================================================

のように &bias ネームリストを加えることで実現できる（&biasは使い勝手が悪かった為廃止予定のシステム。実装の都合上本計算では &bias を用いるが、将来的には改良版の &biasmk2か更にそれを改善したシステムに移り変わる予定であり、&bias についての説明は最小限に留める。）。ここで bias_param3.dat はバイアスのパラメータを記載したファイルであり、

=== bias_param3.dat ============================================================================
&bias
 ndim = 2
 nparam = 5
/
---
20.0d-10 0.5d-10 0.00d-10 0.0d-10 0.0d0
---
877
1 2 3 ...
================================================================================================

のようになっている。 6行目の 20.0d-10 が z^0、0.5d-10 が σ に対応する。 8行目の877 はNa⁺のグループ番号（msファイルの3番目の値）、9行目以下の1 2 3 … は水分子のグループ番号である。

> ./FreeFlex.exe input3.nml &

を実行することで平衡化された初期構造のmsファイル system3.ms が得られる（result/ 以下に上記計算を実行して得られたsystem3.msが置かれているので、計算時間を省略したい場合はこちらを使用すること。）。

US計算

　system3.ms が得られたら、次にREUS計算の準備としてアンブレラサンプリングによる平衡化＆自由エネルギー計算を行う。今回のREUS計算では系に、

$U_{i}^{bias}(z,w)={\frac {k_{B}T}{2}}\left[\left({\frac {z-z_{i}^{0}}{\sigma _{i}^{z}}}\right)^{2}+\left({\frac {w-w_{i}^{0}}{\sigma _{i}^{w}}}\right)^{2}\right]-U^{0}(z,w)$

というバイアス・ポテンシャルをかけた状態でシミュレーションを行う。ここで i はバイアス・ポテンシャルの番号であり、今回は全部で160種類のバイアス・ポテンシャルを用いる。また、 $w$ はwater finger の座標である。 U⁰ はサンプリング効率を上げる為に後から決定するポテンシャルであり、まずは第一項のバイアス・ポテンシャルをかけた状態でUS計算を行う（2015/07/30現在、REUS計算ではリスタートできないため）。このようなバイアス・ポテンシャルは、

=== bias_param4.dat ============================================================================
&bias
 ndim = 2
 nparam = 5
/
---
 9.5d-10 0.5d-10  3.0d-10 0.0d-10 0.0d0
10.5d-10 0.5d-10  3.0d-10 0.0d-10 0.0d0
11.5d-10 0.5d-10  3.0d-10 0.0d-10 0.0d0
  .
  .
  .
---
877
1 2 3 ...
================================================================================================

のような形式のバイアスパラメータファイルを用いれば実現できる。一つ目と二つ目の ”- - -“ の間の数字は各行がバイアス・ポテンシャルの番号 $i$ を意味しており、一列目が $z_{i}^{0}$ 、二列目が $\sigma _{i}^{z}$ 、三列目が $w_{i}^{0}$ 、四列目が $\sigma _{i}^{w}$ を意味する。第五列目は $w$ についてのバイアス・ポテンシャルを入れるか入れないかを指定するパラメータであり、0.0d0 のときは入れず、1.0d0 のとき入れた状態で計算を行う（一部のバイアス・ポテンシャルでは $w$ についてのバイアス・ポテンシャルを除いた状態で計算を実行するため。）。

　water finger の座標については計算の設定の一部がソースコードに直接書き込まれている。 user_bias/watfing2d.F 中の subroutine calculate_xbias において、zmin（考慮する水分子位置の下限）、 zmax（考慮する水分子位置の上限）、wcorr（水分子とイオンの大きさの差の補正、今回は+0.4d-10とする。）、zbulk（バルクと見なす水分子位置の上限）がそれにあたるので適宜訂正すること。

　上記により計算の準備が整ったので、

> export OMP_NUM_THREADS=1; nohup mpiexec –machinefile machinefile4 –n 160 ./FreeFlex.exe input4.nml > output4 &

を実行することにより平衡化を実行できる。なお、input4.nml や machinefile4 はパラメータを一部変更したものであり、内容は実際のソースファイルを参照すること。

　上記プログラムを実行すると作業ディレクトリ中には0000/ ~ 0159/ までのディレクトリが作られる（result/ 以下に計算により得られたsystem4.ms およびbias_data4.datがあるのでMD計算の実行を省略したい場合はこちらを使用すること。）。その中の bias_data4.dat を見てみると、

=== bias_data4.dat ============================================================================
      0.21009730E-08      0.70341072E-09     878     450       2
      0.20860834E-08      0.67332344E-09     878     450       2
      0.20798348E-08      0.65120214E-09     878     450       2
.
.
.
================================================================================================

のようなデータが記載されている。それぞれ第一項目は $z$ の瞬間値、第二項目は $w$ の瞬間値、第三、第四項目は $w$ の両端に対応する重心の番号、第五項目は水和核中に水分子の数に対応する。これらのデータを用いてWHAM計算を行うにはwham.exeを用いる。 wham.exeをコンパイル後、

> ./wham.exe -s 40000 –t 300.0d0 -p bias_param4.dat -o free_energy4.dat -h histogram4.dat 0*/bias_data4.dat

のようなコマンドを実行する。

　ここで -s 40000 は40000個のデータをスキップする（bias_data4.datの出力間隔はデフォルトで10ステップに一回。50万ステップのシミュレーションを実行したので全部で50000個のデータがある。）ことを意味する。-t は温度（K）を設定する（何も指定しないと298.15 Kとなる）。-p bias_param4.datと0*/bias_data4.dat（0000/bias_data4.dat ~ 0159/bias_data4.dat に展開される。）は入力ファイル、-o free_energy4.dat と -h histogram4.dat は出力ファイルを指定するオプションになっている（詳細はwham.F内の説明を参照）。上記コマンドを実行することにより、

ファイル:Free energy.jpg

のような自由エネルギー面を得ることができる。

※gnuplot のインストールとX window システムの設定が必要になる。

> gnuplot
gnuplot> plot sin(x)

によりグラフの表示が可能なことを確認すること。

U⁰の決定と再度の平衡化

　0節で得られた自由エネルギー曲面を U⁰ として利用することでサンプリング効率が増大する。そこで本節では U⁰ を追加したシミュレーションについて説明する。得られた自由エネルギー面は離散データであるが、U⁰ は連続関数である必要があるため、

$U^{0}(z,w)=\sum _{i}^{100}a_{i}exp\left[-{\frac {1}{2}}\left({\frac {(z-\mu _{z,i})^{2}}{s_{z,i}^{2}}}+{\frac {(w-\mu _{w,i})^{2}}{s_{w,i}^{2}}}\right)\right]$

のように自由エネルギー面を100個のガウス関数の重ね合わせでフィッティングする。フィッティング計算は gauss_fit.exe により実行できる。 gauss_fit.exe は2015/07/29 現在、いくつかの直打ちパラメータを含んでいる。入出力ファイル名、alpha、ngauss、nloop、sigma（数式の $s$ に対応。）の初期値及び最小値等がそれにあたるため適宜修正したのちコンパイルする。また、読み込ませる自由エネルギー面のファイル（free_energy4.dat）の一番上の行に “#plot 2 1 40000” という文を追加する必要がある。初めの2は座標の次元、1は座標上のデータの数、40000=200×200 はデータ数に対応する。実行の準備が整ったら、

> ./gauss_fit.exe

を実行する。実行の結果出力ファイルとして、

=== gauss4.dat =================================================================================
#gauss        100         2
   -0.1576857363E+02    0.3024592293E-10    0.2501689972E-08    0.1915664544E-08    0.1676712562E-08
   -0.8884478250E+01    0.1025523130E-08    0.2854676391E-11    0.6924484395E-09    0.1534887776E-09
   -0.4044577445E+01    0.7416068230E-09    0.2602252801E-09    0.1160021206E-09    0.6628152500E-10
.
.
.
================================================================================================

のようなファイルが得られる。はじめの行の #gauss 以下の100はガウス関数の個数、2は座標の次元を意味し、2行目以下はそれぞれ $a_{i},\mu _{(}z,i),\mu _{(}w,i),s_{(}z,i),s_{(}w,i)$ を意味する。

　上記出力ファイルをFreeFlexで利用するには、バイアスパラメータファイルに出力ファイルの結果をコピーすればよい。つまり、

=== bias_param5.dat ============================================================================
  .
  .
  .
---
877
1 2 3 ...

#gauss        100         2
   -0.1576857363E+02    0.3024592293E-10    0.2501689972E-08    0.1915664544E-08    0.1676712562E-08
   -0.8884478250E+01    0.1025523130E-08    0.2854676391E-11    0.6924484395E-09    0.1534887776E-09
   -0.4044577445E+01    0.7416068230E-09    0.2602252801E-09    0.1160021206E-09    0.6628152500E-10
.
.
.
================================================================================================

のようなファイルを作成することになる。後は上記ファイルを用いてFreeFlexを実行することにより平衡化を行う。 wham.exe により自由エネルギーを計算すると0節では穴になってしまっていた領域の曲面を得られていることが分かる。

REUS計算

　前節までの作業により作業ディレクトリ中には0000/system5.ms ~ 0159/system5.ms までの160個のそれぞれ別のバイアス・ポテンシャルのかかった系の最終構造が得られた（これまでの作業を省略したい場合はresult/ 以下のファイルを作業ディレクトリにコピーすること。）。 FreeFlex で REUS計算を行うには、ネームリストファイルに下記のように

=== input6.nml =================================================================================
  .
  .
  .
&reus
  calculate_reus = T
/
================================================================================================

&reus ネームリストを追加すればよい。また、バイアスパラメータファイルに

=== bias_param6.dat ============================================================================
&bias
 ndim = 2
 nparam = 5
 calctype = "REUS"
/
---
 9.5d-10 0.5d-10  3.0d-10 0.0d-10 0.0d0
10.5d-10 0.5d-10  3.0d-10 0.0d-10 0.0d0
.
.
.
---
877
1 2 3 ...

#gauss        100         2
   -0.1576857363E+02    0.3024592293E-10    0.2501689972E-08    0.1915664544E-08    0.1676712562E-08
   -0.8884478250E+01    0.1025523130E-08    0.2854676391E-11    0.6924484395E-09    0.1534887776E-09
.
.
.

NEIGHBOR_LIST
                  2
              1   3
              2   4
  	.
	.
	.
================================================================================================

のように、ネームリスト中に calctype = "REUS" の指定を追加すること、そして、 NEIGHBOR_LIST キーワードとともにレプリカ交換の相手の番号のリストを追加することが必要となる。ここでNEIGHBOR_LIST キーワードの次の行はバイアス番号1（---以下の1行目のバイアス）とバイアス番号2が交換可能なことを意味し、その次の行は、バイアス番号2とバイアス番号1、および、3が交換可能であることを意味する。以下の行も同様でありバイアス・ポテンシャルの数と同じ行数のリストを追加することになる。　以上により計算の準備は完了したので、後は

> export OMP_NUM_THREADS=1; nohup mpiexec –machinefile machinefile6 –n 160 ./FreeFlex.exe input6.nml > output6 &

によりFreeFlexを実行すればよい。

上記の二次元自由エネルギー曲面のデータはfree_energy4.datに格納されており、例えば、等高線を表示するには contour4.plt のようなファイル（詳細は実際のファイルを参照。）を用いて、

> gnuplot contour4.plt

を実行すればよい。実行結果として

ファイル:Contour4plt.jpg

のようなグラフが得られる。また、二次元の自由エネルギー曲面を

$G^{(1)}(z)=-k_{B}Tln\int dwexp[-\beta G^{(2)}(z,w)]$

により一次元に縮約する場合、tools/fep_2d_to_1d.exe が使用できる。ソース中に入出力する自由エネルギー面のファイル名、WHAM計算で使用したデータの分割数（nx1, nx2。デフォルトは200）、および温度（T）を指定する箇所があるので適宜それらを書き換えてコンパイルする。

> ./fep_2d_to_1d.exe

を実行すると作業ディレクトリ中に出力ファイルが作られるので、

> gnuplot
gnuplot>

gnuplot等のグラフソフトで結果を表示することにより、

チュートリアル03:water finger の解析

目次

初期構造の作成

US計算

U⁰の決定と再度の平衡化

REUS計算

案内メニュー

検索

チュートリアル03:water finger の解析

初期構造の作成

US計算

U0の決定と再度の平衡化

REUS計算

案内メニュー

検索

U⁰の決定と再度の平衡化