チュートリアル04:自由エネルギー摂動計算

　ここでは、水分子1つがDCM液体に加わった/消えたときの自由エネルギー差として、溶媒和自由エネルギーを求める。実際の計算は摂動前後の状態の間を内挿する点をおいて、overlapping distribution methodを用いている。

計算の基礎

状態の内挿

系A（始状態）とB（終状態）の自由エネルギー差 $\Delta F=F_{B}-F_{A}$ を求める。摂動前後の状態を内挿するパラメータ λ を考え、λ=0 はA状態、λ=1 はB状態とする。それぞれの状態のポテンシャルエネルギーを U_A, U_B とすると、パラメータ λ の状態でのポテンシャルエネルギーは

$U(\lambda )=(1-\lambda )U_{A}+\lambda U_{B}$

となる。 λ=1 から 0 の間を n-1 分割して（等分でなくてよい）、

$1=\lambda _{1}>\lambda _{2}>\cdots >\lambda _{n}=0$

と中間状態を n-2 個設定して、隣接する2つの中間状態 λ_j と λ_j+1 間の自由エネルギー差 $\Delta F_{j}=F(\lambda _{j+1})-F(\lambda _{j})$ を求める。 状態A (λ=0) からB (λ=1) の差は、それらの隣接した自由エネルギー差の総和をとって求められる。 すなわち、$\Delta F=\sum _{j=1}^{n-1}\Delta F_{j}$

分割数 n は、通常自由エネルギー摂動計算のMPI並列数と一致するようにとられる。また、λ₁, ..., λ_n は、 配列変数 FEP_lambda で指定される。

overlapping distribution

隣接する2つの状態において、$\Delta U=U(\lambda _{j+1})-U(\lambda _{j})$ とする。状態 λ_j の熱平衡における ΔU の確率分布 $p_{0}(\Delta U)$、および状態 λ_j+1 の熱平衡における ΔU の確率分布 $p_{1}(\Delta U)$ をMD計算から求める。そのとき、自由エネルギー差 ΔF_j は overlapping distribution method により、

$\Delta F_{j}=\Delta F_{j}(\Delta U)=k_{\mathrm {B} }T\ln {\frac {p_{1}(\Delta U)}{p_{0}(\Delta U)}}+\Delta U$

として求められる。この左辺は、ΔU の値に依存しない定数となることに注意せよ。 実際の計算では、$\Delta F_{j}$ は $\Delta F_{j}(\Delta U)$ の重み付き平均で計算する。

$\Delta F_{j}=\int d(\Delta U)w(\Delta U)\Delta F_{j}(\Delta U)$

ここでの重み関数 $w(\Delta U)$ は、分布の重なりが大きいところに大きな重みがくるように、以下のようにとる。

$w(\Delta U)={\frac {1}{W}}{\frac {p_{0}(\Delta U)p_{1}(\Delta U)}{p_{0}(\Delta U)+p_{1}(\Delta U)}},\qquad \left(W=\int d(\Delta U){\frac {p_{0}(\Delta U)p_{1}(\Delta U)}{p_{0}(\Delta U)+p_{1}(\Delta U)}}\right)$

相互作用エネルギーのスケーリング（分子間のみ）

状態AからBへの中間状態 λ は、AとBで違いがある原子サイトをそれぞれ FEP_list_A, FEP_list_B のリストに挙げておくことで、 以下のように実装している。 前者はLennard-Jonesの$\varepsilon$およびクーロン力の電荷$q$に$(1-\lambda )$をかけたものを実効的な$\varepsilon$および$q$としてMD計算を行う。後者はそれらに$\lambda$をかけたものを実効的なパラメータとして用いる。

サイト$i$が FEP_list_A に属するとき、

${\begin{aligned}U_{ij}^{\mathrm {LJ} }(r_{ij})&\to U_{ij}^{\mathrm {LJ} }(r_{ij};1-\lambda )=4{\sqrt {(1-\lambda )\varepsilon _{i}\varepsilon _{j}}}{\biggl [}{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{12}-{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{6}{\biggr ]}\\U_{ij}^{\mathrm {Coulomb} }(r_{ij})&\to U_{ij}^{\mathrm {Coulomb} }(r_{ij};1-\lambda )={\frac {(1-\lambda )q_{i}q_{j}}{r_{ij}}}\end{aligned}}$

サイト$i$が FEP_list_B に属するとき、

${\begin{aligned}U_{ij}^{\mathrm {LJ} }(r_{ij})&\to U_{ij}^{\mathrm {LJ} }(r_{ij};\lambda )=4{\sqrt {\lambda \varepsilon _{i}\varepsilon _{j}}}{\biggl [}{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{12}-{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{6}{\biggr ]}\\U_{ij}^{\mathrm {Coulomb} }(r_{ij})&\to U_{ij}^{\mathrm {Coulomb} }(r_{ij};\lambda )={\frac {\lambda q_{i}q_{j}}{r_{ij}}}\end{aligned}}$

したがって中間状態 λ では、FEP_list_A に属するサイトと FEP_list_B に属するサイトの両者ともに中間的な相互作用を及ぼす形で存在している。両端のA (λ=0) と B (λ=1) では、それぞれ一方のみが相互作用を及ぼすことになる。

相互作用エネルギーのスケーリング（分子内も含めて）

FEP は指定したイオンあるいは分子の周りとの相互作用をスケールするものであり、分子内の相互作用はスケールされるべきではない。しかし、FreeFlex での相互作用計算は分子間と分子内を分けずに単にサイト同士の相互作用として計算されている。そのため分子内の相互作用は以下のように計算され、スケールされてしまう。（以下は FEP_list_B に指定された分子の分子内相互作用の例。FEP_list_A に指定された分子の分子内相互作用については $\lambda _{i/j}$ を $1-\lambda _{i/j}$ に置き換えたものである。）

$U_{ij}^{\mathrm {LJ} }(r_{ij})=4{\sqrt {\lambda _{i}\lambda _{j}\varepsilon _{i}\varepsilon _{j}}}{\biggl [}{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{12}-{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{6}{\biggr ]}$

$U_{ij}^{\mathrm {Coulomb} }(r_{ij})={\frac {\lambda _{i}\lambda _{j}q_{i}q_{j}}{r_{ij}}}$

ここで $\lambda _{i}=\lambda _{j}=\lambda$ であるが、それぞれの $\lambda$ が FEP でスケールされる別々のサイトにかかっていることを明示するためにこのように表記した。

これは本来

$U_{ij}^{\mathrm {LJ} }(r_{ij})=4{\sqrt {\varepsilon _{i}\varepsilon _{j}}}{\biggl [}{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{12}-{\biggl (}{\frac {\sigma _{ij}}{r_{ij}}}{\biggr )}^{6}{\biggr ]}$

$U_{ij}^{\mathrm {Coulomb} }(r_{ij})={\frac {q_{i}q_{j}}{r_{ij}}}$

として扱われるべきものである。

　そこで、LJ に関してはサイト間の LJ パラメータを決める際に FEP でスケールされる分子の分子内相互作用に関係する場合は $\varepsilon$ に $\lambda$ をかけないことで対応している。

　一方でクーロンポテンシャルについてはこれまでのスケールされたポテンシャルに対して補正用のポテンシャルを加えた以下のものを新しいクーロンポテンシャルとしている。

$U_{ij}^{\mathrm {Coulomb} }(r_{ij})={\frac {\lambda _{i}\lambda _{j}q_{i}q_{j}}{r_{ij}}}+{\frac {(1-\lambda _{i}\lambda _{j})q_{i}q_{j}}{r_{ij}}}$

これと同様に、電荷と誘起双極子 (CD)、誘起双極子同士 (DD) の相互作用についても以下のように補正している。

$U_{ij}^{\mathrm {CD} }(r_{ij})=\lambda _{i}\mu _{i}\mathbf {T} ^{(1)}\lambda _{j}q_{j}+(1-\lambda _{i}\lambda _{j})\mu _{i}\mathbf {T} ^{(1)}q_{j}$

$U_{ij}^{\mathrm {DD} }(r_{ij})=\lambda _{i}\mu _{i}\mathbf {T} ^{(2)}\lambda _{j}\mu _{j}+(1-\lambda _{i}\lambda _{j})\mu _{i}\mathbf {T} ^{(2)}\mu _{j}$

ここで、$\mathbf {T} ^{(1)}=\mathbf {\nabla } {\frac {1}{\mathbf {r} }}$, $\mathbf {T} ^{(2)}=\mathbf {\nabla } \mathbf {T} ^{(1)}$ である。

これらのポテンシャルの補正項により誘起双極子の形も変わる。

まず、先の式からすべての静電相互作用のポテンシャルは以下のように表せる。

ただし、表記のしやすさのために FEP で相互作用がスケールされるサイトに指定されなかったサイトについては $\lambda _{i}=1$, $\lambda _{j}=1$ と考えることにした。

$U={\frac {1}{2}}\Sigma _{i}\Sigma _{(j\neq i)}U_{ij}^{\mathrm {coulomb} }+\Sigma _{i}\Sigma _{(j\neq i)}U_{ij}^{\mathrm {CD} }-{\frac {1}{2}}\Sigma _{i}\Sigma _{(j\neq i)}U_{ij}^{\mathrm {DD} }+\Sigma _{i}{\frac {|\mathbf {\mu } _{i}|^{2}}{2\alpha _{i}}}$

すると、誘起双極子を決定する SCF 方程式は以下のようになる。

${\frac {\mathrm {d} U}{\mathrm {d} \mathbf {\mu } _{i}}}=0$

$\mathbf {\mu } _{i}=\lambda _{i}\alpha _{i}\Sigma _{(j\neq i)}{\biggl (}-\mathbf {T} ^{(1)}\lambda _{j}q_{j}+\mathbf {T} ^{(2)}\lambda _{j}\mu _{j}{\biggr )}+\alpha _{i}\Sigma _{(j\neq i)}{\biggl (}-\mathbf {T} ^{(1)}(1-\lambda _{i}\lambda _{j})q_{j}+\mathbf {T} ^{(2)}(1-\lambda _{i}\lambda _{j})\mu _{j}{\biggr )}$

第１項はスケールされた電場をスケールされた分極率で感じることを意味する。一方で、第２項はスケールされていない電場をスケールされていない分極率で感じたものを足し、スケールされた相互作用を打ち消すことを意味する。

これらの実装は、$U_{ij}^{\mathrm {coulomb} },U_{ij}^{\mathrm {CD} },U_{ij}^{\mathrm {DD} },\mathbf {\mu } _{i}$ の第１項をこれまでの実装と同様に計算したうえで、FEP でスケールされる分子の分子内相互作用がある部分についてのみ第２項を計算し、足す形になっている。

特に指定をしなければ、FreeFlex では 1-i(>=5) 以上の分子内相互作用を考える仕様になっているため、第２項に含まれるのは FEP でスケールされる分子内のこれに該当するペアである。 ただし、ネームリスト &FreeFlex の calcualte_1_4 という変数で 1-4 相互作用を計算する程度を指定することができるようになっており、 これによって 1-4 相互作用を考える場合には $U_{ij}^{\mathrm {coulomb} },U_{ij}^{\mathrm {CD} },U_{ij}^{\mathrm {DD} },\mathbf {\mu } _{i}$ の第１項、第２項それぞれに 1-4 相互作用を考える程度を示すスケールパラメータがかかることになる。（実装済み）

自由エネルギー計算の収束は、解析プログラム“tools/FEPSHOW.exe”の結果のoverlapping plotで確認する。

計算手続き

この計算は、以下の3つのステップA-C からなる。

A.初期設定

• ネームリスト&FEP をFreeFlex本体の入力ファイルに設定する。

項目	変数の種類	Default	Description
FEP_mode	LOGICAL	.FALSE.	自由エネルギー計算をするかどうか
FEP_lambda	実数配列	""	両端および中間状態の λ の一式。この個数は、mpi並列数と同じとする。 λは 1 から 0 の順に並べること。
FEP_list_A	CHARACTER	""	初期状態で存在するサイトのリスト(*)
FEP_list_B	CHARACTER	""	終状態で存在するサイトのリスト(*)
FEP_output_file	CHARACTER	""	FEP 出力ファイル ($j+1$状態での$\Delta U$).
FEP0_output_file	CHARACTER	""	FEP 出力ファイル ($j$状態での$\Delta U$).

(*) FEP_list_Aに指定したサイトは、初期状態 ($\lambda =\lambda _{n}=0$) で存在し終状態 ($\lambda =\lambda _{1}=1$)で消滅する。FEP_list_Bに指定したサイトは、初期状態 (λ=0) では存在せず終状態 (λ=1) で存在する。自由エネルギー摂動計算では、どちらのリストで指定したサイトもmsファイル中の&systemのSITE_DATA中に指定しておく必要がある。

• 初期状態の .ms ファイルを用意する。(チュートリアル01：水-DCM 界面のシミュレーション を参照）

B.MD計算の実行

• 上の&FEPを含む入力ファイルで、MDを実行する。（チュートリアル01：水-DCM 界面のシミュレーション 、チュートリアル02：並列実行 を参照)

C.結果の解析

• MD計算が終わると、作業用ディレクトリ（&FreeFlex内のdirで指定したもの。何も指定しないとカレントディレクトリ ”./”）の下に0000,0001,…,000N という名前のディレクトリができているはずである。N+1 は、前後と内挿を含めて扱った状態の数で、FEP_lamdba で指定された要素の数と同じである。
• MD計算がディレクトリj内で状態$\lambda _{j}$で走っている時に、状態$\lambda _{j-1}$および状態$\lambda _{j+1}$に変化させた時のエネルギー差がそれぞれディレクトリj内のFEP0_output_file, FEP_outputに出力されている。すなわち、分布としてoverlapしているのはディレクトリjのFEP_output_file(j→j+1)とディレクトリj+1のFEP0_output_file(j+1→j)であることに注意する。

解析プログラムtools/FEPSHOW.exe

　これは上で得られた$\Delta U$の生データを解析して、自由エネルギー差を計算し、収束を確かめるための重なりグラフを出力する。これは、独自の入力ファイルを必要とし、ネームリスト&INPUTの各項目を設定する。入力ファイルは、FEPSHOW.exe実行時の “-i” オプションで引き渡す。

　“tools/FEPSHOW.exe” の入力ネームリスト&INPUT

項目	変数の種類	Default	Description
input_file	CHARACTER	"FEP_output"	FEP出力ファイル名 ($j+1$状態での$\Delta U$).
input_file0	CHARACTER	"FEP0_output"	FEP 出力ファイル名 ($j$状態での$\Delta U$).
n_bin	INTEGER	100	ヒストグラム作成のビン数。
cutoff_range	DOUBLE PRECISION	1.0d-17	ヒストグラムのプロット範囲（横軸）を -cutoff_range ~ cutoff_range 以内に限る。 エネルギー差のプロットで、極端に大きな値を無視するため。単位: J.
n_dir	INTEGER	16	自由エネルギーデータのディレクトリ数 (0000, ...,000n_dir).
leap_line	INTEGER	0	自由エネルギー計算時に先頭から読み飛ばす行数。初期の平衡化の分を除くために用いる。
T	DOUBLE PRECISION	298.15d0	温度(K)
show	LOGICAL	.FALSE.	Δ U 依存性のグラフを出力するかどうか
show_mode	INTEGER	1	上でshow=.TRUE.のとき 1: $p_{0}(\Delta U),p_{1}(\Delta U),w(\Delta U)$ と自由エネルギー $\Delta F(\Delta U)$を出力 2: $p_{0}(\Delta U),p_{1}(\Delta U),w(\Delta U)$ のみ出力 3: 自由エネルギー$\Delta F(\Delta U)$のみ出力
output_graph	CHARACTER	"overlapping.dat"	重なりグラフのデータの出力先ファイル

実行例－ 水の溶媒和自由エネルギー

　sample/tutorial/04 ディレクトリには、ジクロロメタン)(DCM)溶液中の水の自由エネルギーを計算するため、以下の3つのファイルがある。

• input.nml: FreeFlex本体の入力ファイル。
• w1along.ms: FreeFlex本体のmsファイル。水1分子とDCM598分子の初期配置。
• nmlfep: tools/FEPSHOW.exeの入力ファイル。

実行例は以下の通りである。 この際、初期構造は全てのディレクトリに共通のものを用いる必要がある。 これはFEPを用いた後のrestartファイルや、outputのmsファイルはFEPによりスケーリングを受けた後のパラメータとなっているためである。

(1) 実行ディレクトリに移る。

> cd FreeFlex/sample/tutorial/04

(2) MD計算を実行する。

> nohup mpiexec -machinefile machinefile -n 8 ../../../FreeFlex.exe input.nml > logfile &

• MPI並列実行のノード数 (-n 8) は、λ で分割した状態数（ネームリスト &FEP 中の FEP_lambda の要素の数）と同じにする。

• logfile は FreeFlex用の標準出力ファイル。

(3) MD計算結果の確認

MD計算実行後に、実行ディレクトリの下にノード番号に対応するディレクトリ (0000, 0001, ..., 0007) が確認される。

この各ディレクトリ内に、以下の2つのファイルが出力されている。

* FEP0_output --- ネームリスト &FEP 中の FEP0_output_file で指定したファイル。
* FEP_output --- FEP_output_file で指定したファイル。

(4) tools/FEPSHOW.exe による解析を実行する。

> ../../../tools/FEPSHOW.exe -i nmlfep

標準出力のなかで、隣接する λ 状態間の自由エネルギー差と標準偏差、 および全体の自由エネルギー差と標準偏差が表示される。

********** STEP  1 **********
 (1)l0= 0.00E+00 -> l1= 1.00E-04
  delta_G (kcal/mol):    1.26004731950368      +-  0.226145296970347

********** STEP  2 **********
 (2)l0= 1.00E-04 -> l1= 1.00E-03
  delta_G (kcal/mol):   0.445784144692485      +-  0.113937233985755

********** STEP  3 **********
 (3)l0= 1.00E-03 -> l1= 1.00E-02
  delta_G (kcal/mol):   0.371481706915616      +-  0.117521396550999

********** STEP  4 **********
 (4)l0= 1.00E-02 -> l1= 1.00E-01
  delta_G (kcal/mol):  -7.107622766591130E-002 +-  0.127231527920189

********** STEP  5 **********
 (5)l0= 1.00E-01 -> l1= 3.33E-01
  delta_G (kcal/mol):  -0.755530734850856      +-  0.132509635469720

********** STEP  6 **********
 (6)l0= 3.33E-01 -> l1= 6.66E-01
  delta_G (kcal/mol):   -1.49635701253196      +-  0.239194458862602

********** STEP  7 **********
 (7)l0= 6.66E-01 -> l1= 1.00E+00
  delta_G (kcal/mol):   -2.16760102634211      +-  0.136389226627518

 === Calculated total delta_G (kcal/mol) ===
  -2.41325183027906      +-  0.433009633561746

入力ファイル (nmlfep) 中で show=.TRUE. を指定したときには、ディレクトリ 0001, ..., 0007 中に 結果の出力ファイル overlapping.dat (ネームリスト &FEP 中の output_graph で指定したファイル）ができる。 （ディレクトリ 0000 にはできないので注意。これは状態の隣接ペアの数は、状態数よりも1つ少なくなるため。）

エネルギー差 ΔU 、　状態0での分布 p0(ΔU)、　状態1での分布 p1(ΔU)、　重み関数 (p0 p1)/(p0+p1)、　自由エネルギー差 ΔF(ΔU)
  -9.8464711561E-21    0.0000000000E+00    0.0000000000E+00    0.0000000000E+00    0.0000000000E+00
   1.6838274198E-22    1.4139927387E-01    8.5537500000E-01    6.3815850960E-01    8.9127007844E-01
   1.0183236640E-20    1.1618094397E-01    1.3857500000E-01    3.3236640178E-01    1.2830625000E+00
   2.0198090538E-20    6.8753475289E-02    5.7500000000E-03    2.7906644622E-02    1.1742812522E+00
   3.0212944436E-20    5.0469368397E-02    3.0000000000E-04    1.5684439995E-03    1.0720740910E+00
   4.0227798334E-20    3.6535505184E-02    0.0000000000E+00    0.0000000000E+00    0.0000000000E+00
   ...

また、show=.TRUE. のときには、これらのグラフが画面にポップアップする。表示されるグラフの種類は、show_mode =1, 2, 3 によって変わる。