Weighted Histogram Analysis Method

1.Multiple Histogram Reweighting

カノニカルアンサンブルにおいてある座標$Q(\mathbf {r} ^{N})$ に関してバイアスポテンシャル$W_{i}(Q)$ を加えることを 考える。ここでi はトラジェクトリーの番号であり，それぞれのトラジェクトリーで違うバイアスを掛けて いる。この時，それぞれのトラジェクトリーから直接計算されるQ についてのバイアス付きの分布関数は，

$p_{i}(Q)={\frac {1}{Z_{i}}}\int d\mathbf {r} ^{N}exp\left[-\beta \left\{U(\mathbf {r} ^{N})+W_{i}(Q(\mathbf {r} ^{N}))\right\}\right]\delta (Q-Q(\mathbf {r} ^{N}))$

(1)

$={\frac {1}{Z_{i}}}exp\left[-\beta W_{i}(Q)\right]\int d\mathbf {r} ^{N}exp[-\beta U(\mathbf {r} ^{N})]\delta (Q-Q(\mathbf {r} ^{N}))$

(2)

$Z_{i}=\int d\mathbf {r} ^{N}exp[-\beta {U(\mathbf {r} ^{N})+W_{i}(Q(\mathbf {r} ^{N}))}]$

(3)

と表わせる。一方，バイアスなしの場合の分布関数は,

$p_{0}(Q)={\frac {1}{Z_{0}}}\int d\mathbf {r} ^{N}exp[-\beta U(\mathbf {r} ^{N})]\delta (Q-Q(\mathbf {r} ^{N}))$

(4)

$Z_{0}=\int d\mathbf {r} ^{N}exp[-\beta U(\mathbf {r} ^{N})]$

(5)

であるから，バイアス付きの分布関数とは

${\begin{aligned}p_{0}(Q)&={\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]{\frac {1}{Z_{i}}}exp[-\beta W_{i}(Q)]\int d\mathbf {r} ^{N}exp[-\beta U(\mathbf {r} ^{N})]\delta (Q-Q(\mathbf {r} ^{N}))\\&={\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]p_{i}(Q)\end{aligned}}$

(6)

の関係をもつ。

　式6 は理想的には$p_{i}(Q)$ からでも同じ$p_{0}(Q)$ を計算できることを意味している。勿論，有限のトラジェ クトリーでは$p_{i}(Q)$ がほぼ0 となる領域については正しくサンプリングできないため，実際にはそうはいか ない。そこでそれぞれのトラジェクトリーから実際に得られる分布関数

$p_{i}^{est}(Q)\delta Q={\frac {H_{i}(Q)}{M_{i}}}$

(7)

（$H_{i}(Q)$ はQ に関するヒストグラム，$\delta Q$はヒストグラムの幅，$M_{i}$ は全サンプル数）の線形結合

$p_{0}^{est}(Q)=\sum _{i}w_{i}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]p_{i}^{est}(Q)$

(8)

$\sum _{i}w_{i}(Q)=1$

(9)

を取り，最も真の分布p0(Q) に近くなるよう重み関数$w_{i}(Q)$ を決定することにする。

　式8 が最も真の分布に近くなるのは、 $p_{0}^{est}(Q)$の分散が最小になる時である。$p_{0}^{est}(Q)$の分散は

$\sigma ^{2}[p_{0}^{est}(Q)]=<p_{0}^{est}(Q)^{2}>-<p_{0}^{est}(Q)>^{2}$

(10)

$=\sum _{i}\sum _{j}w_{i}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]w_{j}(Q){\frac {Z_{j}}{Z_{0}}}exp[\beta W_{j}(Q)]\quad \times \quad \left(<p_{i}^{est}(Q)p_{j}^{est}(Q)>-<p_{i}^{est}(Q)><p_{j}^{est}(Q)>\right)$

(11)

$=\sum _{i}\left[w_{i}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]\right]^{2}\left(<p_{i}^{est}(Q)^{2}>-<p_{i}^{est}(Q)>^{2}\right)$

(12)

$=\sum _{i}\left[w_{i}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]\right]^{2}{\frac {<H_{i}(Q)^{2}>-<H_{i}(Q)>^{2}}{(M_{i}\delta Q)^{2}}}$

(13)

と表わされる。ヒストグラムを作るために使用したサンプルの間に相関がない場合，あるサンプルが$[Q,Q+\delta Q)$ に入る確率，つまり，$H_{i}(Q)$ としてカウントされる確率は$p_{i}(Q)\delta Q$となるから、$H_{i}(Q)=k$ と なる確率は二項分布

$_{M_{i}}C_{k}[p_{i}(Q)\delta Q]^{k}[1-p_{i}(Q)\delta Q]^{M_{i}-k}$

(14)

に従う。従ってその分散は

$<H_{i}(Q)^{2}>-<H_{i}(Q)>^{2}=M_{i}p_{i}(Q)\delta Q[1-p_{i}(Q)\delta Q]\fallingdotseq M_{i}p_{i}(Q)\delta Q$

(15)

となる。式15 を式13 に代入すると，

$\sigma ^{2}[p_{0}^{est}(Q)]=\sum _{i}\left[w_{i}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]\right]^{2}{\frac {p_{i}(Q)}{M_{i}\delta Q}}$

(16)

$=p_{0}(Q)\sum _{i}w_{i}^{2}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]{\frac {1}{M_{i}\delta Q}}$

(17)

が得られる。後はLagrange の未定乗数法により式9 の条件のもと，式17 を最小化する$w_{i}(Q)$ を求め ればよい。Lagrange の未定乗数法より

${\frac {\partial }{\partial w_{i}(Q)}}\left[\sigma ^{2}[p_{0}^{est}(Q)]-2\alpha \left\{\sum _{i}w_{i}(Q)-1\right\}\right]$

(18)

$=p_{0}(Q)2w_{i}(Q){\frac {Z_{i}}{Z_{0}}}exp[\beta W_{i}(Q)]{\frac {1}{M_{i}\delta Q}}-2\alpha =0$

(19)

従って，

$w_{i}(Q)=\alpha {\frac {1}{p_{0}(Q)}}{\frac {Z_{0}}{Z_{i}}}M_{i}\delta Qexp[-\beta W_{i}(Q)]$

(20)

式9 からα は規格化定数なので

$w_{i}(Q)={\frac {M_{i}exp[-\beta W_{i}(Q)]/Z_{i}}{\sum _{j}M_{j}exp[-\beta W_{j}(Q)]/Zj}}$

(21)

よって，求めるべき分布関数は式8 から

$p_{0}^{est}(Q)={\frac {1}{\delta Q}}{\frac {\sum _{i}H_{i}(Q)}{\sum _{j}M_{j}exp[-\beta W_{j}(Q)]Z_{0}/Zj}}$

(22)

となる。これで未知変数は${\frac {Z_{i}}{Z_{0}}}$ のみとなった。式3,4 から

${\frac {Z_{i}}{Z_{0}}}={\frac {1}{Z_{0}}}\int d\mathbf {r} ^{N}exp[-\beta \{U(\mathbf {r} ^{N})+W_{i}(Q(\mathbf {r} ^{N}))\}]$

(23)

$={\frac {1}{Z_{0}}}\int d\mathbf {r} ^{N}\int dQexp[-\beta U(\mathbf {r} ^{N})]\delta (Q-Q(\mathbf {r} ^{N}))exp[-\beta W_{i}(Q)]$

(24)

$=\int dQp_{0}(Q)exp[-\beta W_{i}(Q)]$

(25)

となるので，$p_{0}(Q)=p_{0}^{est}(Q)$として式22 と式25 をiterative に解くことで$p_{0}^{est}(Q)$を計算できる。自由エネルギー曲線は，

$F(Q)=-k_{B}Tlnp_{0}^{est}(Q)$

(26)

により得られる。

2. 鏡像電荷によるポテンシャル

　誘電体の境界近傍の電荷は鏡像電荷によるポテンシャルを感じる。本章では，そのポテンシャルの大きさ について議論する。

　誘電体$\varepsilon _{0}$,$\varepsilon _{1}$ の境界を$x=0$ にとり，$\varepsilon _{0}$側$x=a$ の位置に電荷q が存在する系を考える。鏡像電荷の方法によると，観測者にはそれぞれの誘電体$\varepsilon _{0}$,$\varepsilon _{1}$ 中で図のような電荷q, q', q" が存在しているように見える。

　そこでまずはq', q" をq を用いて表すことを目指す。

　$x>0$　の任意の点$\mathbf {r}$ における電場は左図より，

$\mathbf {E} _{0}(\mathbf {r} )={\frac {1}{4\pi \varepsilon _{0}}}\left[{\frac {q}{(\mathbf {r} -a{\hat {x}})^{2}}}{\frac {\mathbf {r} -a{\hat {x}}}{|\mathbf {r} -a{\hat {x}}|}}+{\frac {q'}{(\mathbf {r} +a{\hat {x}})^{2}}}{\frac {\mathbf {r} +a{\hat {x}}}{|\mathbf {r} +a{\hat {x}}|}}\right]$

(27)

と表わされる。一方$x<0$ 　における電場は右図より，

$\mathbf {E} _{1}(\mathbf {r} )={\frac {1}{4\pi \varepsilon _{1}}}{\frac {q''}{(\mathbf {r} -a{\hat {x}})^{2}}}{\frac {\mathbf {r} -a{\hat {x}}}{|\mathbf {r} -a{\hat {x}}|}}$

(28)

となる。

　電場の接続条件（$\nabla \times E=0$　より）から，境界面において水平方向の電場は等しくなる（$E_{0,y}(0,y,z)=E_{1,y}(0,y,z)$）ので

${\frac {1}{4\pi \varepsilon _{0}}}{\frac {(q+q')y}{(a^{2}+y^{2}+z^{2})^{\frac {3}{2}}}}={\frac {1}{4\pi \varepsilon _{1}}}{\frac {q''y}{(a^{2}+y^{2}+z^{2})^{\frac {3}{2}}}}$

(29)

${\frac {(q+q')}{\varepsilon _{0}}}={\frac {q''}{\varepsilon _{1}}}$

(30)

が得られる。また，電束密度の接続条件($\nabla \cdot D=0$　より）から，境界面において垂直方向の電束密度は等し くなる（$D_{0,x}(0,y,z)=D_{1,x}(0,y,z)$）ので，

${\frac {1}{4\pi }}{\frac {(-q+q')a}{(a^{2}+y^{2}+z^{2})^{\frac {3}{2}}}}={\frac {1}{4\pi }}{\frac {-q''a}{(a^{2}+y^{2}+z^{2})^{\frac {3}{2}}}}$

(31)

$-q+q'=-q''$

(32)

が得られる。式30 ，32 から求める電荷は，

$q'={\frac {\varepsilon _{0}-\varepsilon _{1}}{\varepsilon _{0}+\varepsilon _{1}}}q$

(33)

$q''={\frac {2\varepsilon _{1}}{\varepsilon _{0}+\varepsilon _{1}}}q$

(34)

となる。

　q', q" が求まったので電荷q にかかるポテンシャルを計算しよう。電荷q はq' から

$F_{x}(a)=-{\frac {1}{4\pi \varepsilon _{0}}}{\frac {\varepsilon _{0}-\varepsilon _{1}}{\varepsilon _{0}+\varepsilon _{1}}}{\frac {q^{2}}{4a^{2}}}$

(35)

の力を受ける。従ってq を無限遠から$x=a$まで持っていくために必要な仕事（＝ポテンシャルエネル ギー）は，

$U_{mirror}(a)=\int _{\infty }^{a}F_{x}(a')da'={\frac {1}{4\pi \varepsilon _{0}}}{\frac {\varepsilon _{0}-\varepsilon _{1}}{\varepsilon _{0}+\varepsilon _{1}}}{\frac {q^{2}}{4a^{2}}}$

(36)

となる。鏡像電荷の移動には仕事を伴わないので，電荷q とq' が距離2a だけ離れた電荷間に働くポテン シャルエネルギー

$U(a)={\frac {1}{4\pi \varepsilon _{0}}}{\frac {\varepsilon _{0}-\varepsilon _{1}}{\varepsilon _{0}+\varepsilon _{1}}}{\frac {q^{2}}{2a^{2}}}$

(37)

の半分になることに注意せよ。