グランドカノニカルシミュレーションにおける化学ポテンシャルと気圧の関係性

グランドカノニカルアンサンブル（μVTアンサンブル）のシミュレーションでは温度 $T$ と化学ポテンシャル $\mu$ をパラメータとして与える必要があるものの、実験では化学ポテンシャルではなくリザーバーの気圧 $P$ が固定条件となる場合が多く、実験値とシミュレーション結果を比較するためにはリザーバーの温度・化学ポテンシャルと気圧の関係が必要となります。

# 分配関数（カノニカルアンサンブル）
カノニカルアンサンブル（NVTアンサンブル）における $N$ 個の原子を含む系の分配関数 $Q(N, V, T)$ は次の通り。
$$Q(N, V, T) = \frac{V^N}{\Lambda^{3N}N!}\int d\mathbf{s}^N \exp\left\\{-\beta U(\mathbf{s}^N)\right\\}\tag{1}$$
ここで、$\beta=\frac{1}{k_BT}$ は逆温度、$\Lambda=\frac{h}{\sqrt{2\pi mk_BT}}$ は熱的ドブロイ波長、$\mathbf{s}^N$は $N$ 粒子の座標です。


また、自由エネルギー $F$ は分配関数を用いて次のように書けます。
$$F=-\frac{1}{\beta}\ln Q(N, V, T)\tag{2}$$
これから、化学ポテンシャル $\mu$ は次のようになります。
$$\mu\equiv\frac{\partial F}{\partial N}\simeq-\frac{1}{\beta}\ln\left[\frac{Q(N+1, V, T)}{Q(N, V, T)}\right]\tag{3}$$

また、$M$ 原子分子を $N$ 個含む系の分配関数 $Q(N, M, V, T)$ は
$$Q(N, M, V, T) = \frac{q(T)^NV^N}{N!}\prod_{i=1}^M\int d\mathbf{s}_i^N \exp\left\\{-\beta U(\mathbf{s}_i^N)\right\\}\tag{4}$$
$q(T)$ は分子の分配関数の一部で、単原子分子の場合は $q(T)=\frac{1}{\Lambda^3}$ となります。また、$\mathbf{s}^N$ は分子内の原子のデカルト座標です。

# 理想気体の化学ポテンシャル
任意の系において $(密度)\to0$ の極限をとると理想気体として振舞うようになり、そのとき分子内相互作用のみが全ポテンシャルエネルギーに寄与します。
$$U\approx\sum_{i=1}^N U^{intra}(i)\tag{5}$$

### 単原子分子の場合
式$(1)$のように表されていた分配関数が、単原子分子のみからなる系では次のように短縮されます。
$$Q_{IG}(N, V, T)=\frac{V^N}{\Lambda^{3N}N!}\tag{6}$$
このような単原子分子の化学ポテンシャル $\mu_{ideal}$ は、分子の個数密度 $\rho$ を用いて
$$\mu_{ideal}=\mu_{ideal}^0+k_BT\ln\rho\tag{7}$$
と書くことができます。ここで、基準状態における化学ポテンシャル $\mu_{ideal}^0$ は次のように定義されます。
$$\mu_{ideal}^0\equiv k_BT\ln\Lambda^3\tag{8}$$

### 多原子分子の場合
式$(4)$で表されていた分配関数が、分子同士の相互作用がない系では次のように短縮されます。
$$Q_{ideal}(N, M, V, T)=\frac{q(T)^NV^N}{N!}\left\[\prod_{i=1}^M\int d\mathbf{s}_i\exp\\{-\beta U^{intra}(\mathbf{s}_i)\\}\right\]^N\tag{9}$$
これを式(3)に代入して、化学ポテンシャルは

$$\mu_{ideal}=\mu_{ideal}^0+k_BT\ln\rho\tag{10}$$
と書くことができます。ここで、基準状態における化学ポテンシャル $\mu_{ideal}^0$ は次のように定義されます。

```math
\begin{align}
\mu_{ideal}^0 &\equiv -k_BT\ln q(T)+\mu_{intra}^0\\
&=-k_BT\ln \left[q(T)\prod_{i=1}^M\int d\mathbf{s}_i\exp\left\{-\beta U^{intra}(\mathbf{s}_i)\right\}\right]\tag{11}
\end{align}
```
ここで、 $\mu_{ideal}^0$ は温度のみの関数となります。また、どんな温度を基準状態としても化学ポテンシャルは一定数シフトするだけで、系の観測可能な熱力学的性質には影響を及ぼしません。

# GCMC
グランドカノニカルモンテカルロ（GCMC）シミュレーションでは、以下の採択ルールに従います。
### 粒子の挿入
$\Delta U=U(N+1)-U(N)$ とすると、粒子の挿入の採択確率 $P_{accept}(N\to N+1)$ は、
$$P_{accept}(N\to N+1)=\min\left[1, \frac{Vq(T)}{N+1}\exp[-\beta\\{\Delta U-(\mu^B-\mu_{intra}^0)\\}]\right]\tag{12}$$
となります。
### 粒子の除去
$\Delta U=U(N-1)-U(N)$ とすると、粒子の除去の採択確率 $P_{accept}(N\to N-1)$ は、
$$P_{accept}(N\to N-1)=\min\left[1, \frac{N}{Vq(T)}\exp[-\beta\\{\Delta U+(\mu^B-\mu_{intra}^0)\\}]\right]\tag{13}$$
となります。
### 背景
これらの採択ルールは、気体分子が同じ分子を同じ化学ポテンシャルで含むリザーバーと交換されるという考え方に基づいています。ただし、リザーバーでは分子同士が相互作用しないという点が系内の状況とは異なります。
吸着現象など現実的な状況において、これは希薄な蒸気（リザーバー）と平衡状態にある高密度相に対応します。リザーバーの気圧は、先に述べた通りGCMCの条件設定において重要な量です。
### 気圧との関係性
リザーバー内の気圧と化学ポテンシャルの関係は次のようになります。

```math
\begin{align}
\mu^B&\equiv\mu_{ideal}^0+k_BT\ln\rho\\
&=\mu_{ideal}^0+k_BT\ln(\beta P_{ideal})\tag{14}
\end{align}
```
この関係式を採択ルールの式に代入すると、
$$P_{accept}(N\to N+1)=\min\left[1, \frac{\beta P_{ideal}V}{N+1}\exp[-\beta\Delta U]\right]\tag{15}$$
$$P_{accept}(N\to N-1)=\min\left[1, \frac{N}{\beta P_{ideal}V}\exp[-\beta\Delta U]\right]\tag{16}$$
つまり、系が理想気体と見なせるリザーバーと平衡状態にあるときはリザーバーの気圧のみが採択ルールの式に表れ、基準状態に関する情報は（当然ながら）入ってこないことがわかります。
### 実在気体の場合
リザーバーが理想気体と見なせないほどに高圧な場合は、次に示すリザーバーの化学ポテンシャルと気圧の関係式を用いる必要があります。
$$\mu_B=\mu_{ideal}^0+k_BT\ln(\beta P\phi)\tag{17}$$
ここで $\phi$ はリザーバー内流体のフガシティ係数で、これは温度と圧力の関数です。また、フガシティ係数はリザーバー内の蒸気の状態方程式から直接計算することができます。まとめると、非理想気体においては式(15), (16)の $P_{ideal}$ を $P\phi$ で置き換えればよいことになります。
# 参考文献
Daan Frenkel and Berend Smit. 1996. Understanding Molecular Simulation: From Algorithms to Applications (1st. ed.). Academic Press, Inc., USA.