Weighted Histogram Analysis Method - 版の履歴

Hirano: ページの作成:「== 1.Multiple Histogram Reweighting == カノニカルアンサンブルにおいてある座標 $Q(\mathbf{r}^N)$ に関してバイアスポテンシャル
2026-05-26T04:13:50Z

ページの作成:「== 1.Multiple Histogram Reweighting == カノニカルアンサンブルにおいてある座標<math>Q(\mathbf{r}^N)</math> に関してバイアスポテンシャル<ma…」

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

{{NumBlk|::|<math>
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))
</math>|{{EquationRef|1}}}}

{{NumBlk|::|<math>
= \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))
</math>|{{EquationRef|2}}}}

{{NumBlk|::|<math>
Z_i =\int d \mathbf{r}^N exp[-\beta{ U( \mathbf{r}^N) +W_i(Q( \mathbf{r}^N))}]
</math>|{{EquationRef|3}}}}

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

{{NumBlk|::|<math>
p_0(Q) = \frac{1}{Z_0}\int d \mathbf{r}^N exp[- \beta U(\mathbf{r}^N)] \delta(Q-Q(\mathbf{r}^N))
</math>|{{EquationRef|4}}}}

{{NumBlk|::|<math>
Z_0 =\int d\mathbf{r}^N exp[-\beta U(\mathbf{r}^N)]
</math>|{{EquationRef|5}}}}

であるから，バイアス付きの分布関数とは
{{NumBlk|::|<math>\begin{align}
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{align}</math>|{{EquationRef|6}}}}

の関係をもつ。

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

{{NumBlk|::|<math>
p^{est}_i(Q)\delta Q = \frac{H_i(Q)} {M_i}
</math>|{{EquationRef|7}}}}

（<math>H_i(Q)</math> はQ に関するヒストグラム，<math>\delta Q </math>はヒストグラムの幅，<math>M_i</math> は全サンプル数）の線形結合

{{NumBlk|::|<math>
p^{est}_0 (Q) =\sum_i w_i(Q)\frac{Z_i}{Z_0} exp [\beta W_i(Q)] p^{est}_i (Q)
</math>|{{EquationRef|8}}}}

{{NumBlk|::|<math>
\sum_i w_i(Q) = 1
</math>|{{EquationRef|9}}}}

を取り，最も真の分布p0(Q) に近くなるよう重み関数<math> w_i(Q)</math> を決定することにする。

　式{{EquationNote|8}} が最も真の分布に近くなるのは、 <math>p^{est}_0 (Q) </math>の分散が最小になる時である。<math>p^{est}_0 (Q) </math>の分散は

{{NumBlk|::|<math>
\sigma ^2[p^{est}_0 (Q)] = < p^{est}_0 (Q)^2>- < p^{est}_0 (Q)>^2
</math>|{{EquationRef|10}}}}

{{NumBlk|::|<math>
=\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^{est}_i (Q) p^{est}_j (Q)> - < p^{est}_i (Q)> < p^{est}_j (Q) > \right)
</math>|{{EquationRef|11}}}}

{{NumBlk|::|<math>
=\sum_i \left[ w_i(Q) \frac{Z_i}{Z_0} exp [\beta W_i(Q)] \right]^2 \left( < p^{est}_i (Q)^2> - < p^{est}_i (Q)>^2 \right)
</math>|{{EquationRef|12}}}}

{{NumBlk|::|<math>
=\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 }
</math>|{{EquationRef|13}}}}

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

{{NumBlk|::|<math>
_{M_i} C_k [p_i(Q) \delta Q]^k [1-p_i(Q)\delta Q]^{M_i-k}
</math>|{{EquationRef|14}}}}

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

{{NumBlk|::|<math>
< 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
</math>|{{EquationRef|15}}}}

となる。式{{EquationNote|15}} を式{{EquationNote|13}} に代入すると，

{{NumBlk|::|<math>
\sigma ^2[p^{est}_0 (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}
</math>|{{EquationRef|16}}}}

{{NumBlk|::|<math>
= 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}
</math>|{{EquationRef|17}}}}

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

{{NumBlk|::|<math>
\frac{\partial}{\partial w_i(Q)}\left[ \sigma ^2[p^{est}_0 (Q)] -2 \alpha \left\{\sum_i w_i(Q) - 1 \right\} \right]
</math>|{{EquationRef|18}}}}

{{NumBlk|::|<math>
= 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
</math>|{{EquationRef|19}}}}

従って，

{{NumBlk|::|<math>
w_i(Q) = \alpha \frac{1}{p_0(Q)} \frac{Z_0}{Z_i} M_i\delta Q exp [- \beta W_i(Q)]
</math>|{{EquationRef|20}}}}

式{{EquationNote|9}} からα は規格化定数なので

{{NumBlk|::|<math>
w_i(Q) = \frac{ M_i exp [- \beta W_i(Q)] /Z_i } {\sum_j M_j exp [ -\beta W_j(Q)] /Zj }
</math>|{{EquationRef|21}}}}

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

{{NumBlk|::|<math>
p^{est}_0 (Q) =\frac{1}{\delta Q} \frac{\sum_i H_i(Q)} {\sum_j M_j exp [-\beta W_j(Q)] Z_0/Zj }
</math>|{{EquationRef|22}}}}

となる。これで未知変数は<math>\frac{Z_i}{Z_0}</math> のみとなった。式{{EquationNote|3}},{{EquationNote|4}} から

{{NumBlk|::|<math>
\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)) \}]
</math>|{{EquationRef|23}}}}

{{NumBlk|::|<math>
=\frac{1}{Z_0}\int d \mathbf{r}^N \int dQ exp[-\beta U(\mathbf{r}^N)] \delta(Q - Q(\mathbf{r}^N)) exp [-\beta W_i(Q)]
</math>|{{EquationRef|24}}}}

{{NumBlk|::|<math>
=\int dQ p_0(Q) exp [-\beta W_i(Q)]
</math>|{{EquationRef|25}}}}

となるので，<math>p_0(Q) = p^{est}_0 (Q) </math>として式{{EquationNote|22}} と式{{EquationNote|25}} をiterative に解くことで<math>p^{est}_0 (Q)</math>を計算できる。自由エネルギー曲線は，

{{NumBlk|::|<math>
F(Q) = -k_B T ln p^{est}_0 (Q)
</math>|{{EquationRef|26}}}}

により得られる。

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

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

[[ファイル:Wham-2-electricfield.jpg]]

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

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

{{NumBlk|::|<math>
\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]
</math>|{{EquationRef|27}}}}

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

{{NumBlk|::|<math>
\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} |}
</math>|{{EquationRef|28}}}}

となる。

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

{{NumBlk|::|<math>
\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} } }
</math>|{{EquationRef|29}}}}

{{NumBlk|::|<math>
\frac{(q + q')} {\varepsilon_0} = \frac{q''} {\varepsilon_1}
</math>|{{EquationRef|30}}}}

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

{{NumBlk|::|<math>
\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} } }
</math>|{{EquationRef|31}}}}

{{NumBlk|::|<math>
-q + q' = -q''
</math>|{{EquationRef|32}}}}

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

{{NumBlk|::|<math>
q'=\frac{ \varepsilon_0 - \varepsilon_1 } {\varepsilon_0 + \varepsilon_1}q
</math>|{{EquationRef|33}}}}

{{NumBlk|::|<math>
q''=\frac{ 2 \varepsilon_1 } {\varepsilon_0 + \varepsilon_1}q
</math>|{{EquationRef|34}}}}

となる。

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

{{NumBlk|::|<math>
F_x(a) =-\frac{1}{4 \pi \varepsilon_0} \frac{ \varepsilon_0 - \varepsilon_1 } {\varepsilon_0 + \varepsilon_1} \frac{q^2}{4a^2}
</math>|{{EquationRef|35}}}}

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

{{NumBlk|::|<math>
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}
</math>|{{EquationRef|36}}}}

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

{{NumBlk|::|<math>
U(a) =\frac{1}{4 \pi \varepsilon_0} \frac{ \varepsilon_0 - \varepsilon_1 } {\varepsilon_0 + \varepsilon_1} \frac{q^2}{2a^2}
</math>|{{EquationRef|37}}}}

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