Statistical Mechanics: Theory and Molecular Simulation Chapter 5 - The isobaric ensemble
Statistical Mechanics: Theory and Molecular SimulationはMark Tuckermanによって著された分子シミュレーションに関連した熱・統計物理学，及び量子力学の参考書になります。
現時点では（恐らく）和訳がないこともあり，日本での知名度はあまりないような気がしますが，被引用数が1300を超えていることを考えると海外では高く認知されている参考書みたいです。
本記事では，Chapter 5（The isobaric ensemble）の章末問題の和訳とその解答例を紹介します。解答例に間違いが見受けられた場合はお知らせいただけると助かります。
Problem 5.1
等温等圧(NPT)アンサンブルの分布関数は熱浴と機械的ピストンの両方と接続している系のミクロカノニカル描像から導けることを示せ。
系の粒子数，体積，ハミルトニアンをそれぞれ$${N_1,V_1,\mathcal{H}_1(\mathbf{x}_1)}$$，熱浴の粒子数，体積，ハミルトニアン，温度をそれぞれ$${N_2,V_2,\mathcal{H}_2(\mathbf{x}_2),T}$$とする。
熱浴と系は共に圧力$${P}$$の機械的ピストンと接続しており，体積$${V_1, V_2}$$は可変であるとする。
熱浴は系に対して十分に巨視的である（$${N_1\ll N_2, V_1\ll V_2}$$）とする。
このとき，系＋熱浴に対してNPHアンサンブル（$${N=N_1+N_2}$$）の分配関数$${\Gamma(N,P,H)}$$を考えると，
$$
\begin{align*}
\Gamma(N,P,H)&=\mathcal{M}\int_{0}^{\infty}{\rm d}V\int_{D(V_1)}{\rm d}\mathbf{x}_1\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+\mathcal{H}_2(\mathbf{x}_2)\right\}+P\left\{V_1+V_2\right\}-H\right)\\
&=\mathcal{M}\int_{0}^{\infty}{\rm d}V\int_{D(V_1)}{\rm d}\mathbf{x}_1\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}+\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\\
&=\mathcal{M}\int_{0}^{\infty}{\rm d}V_2\int_{D(V-V_2)}{\rm d}\mathbf{x}_1\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}+\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\\
\end{align*}
$$
となる。
系の分布関数$${f(\mathbf{x}_1)}$$は
$$
\begin{align*}
f(\mathbf{x}_1)&\propto\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}+\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\\
\end{align*}
$$
となる。
両辺の対数を取り，$${\mathcal{H}_1+PV_1\ll \mathcal{H}_2+PV_2}$$が成立することを利用して$${\mathcal{H}_1+PV_1}$$の１次の項までで近似すると，
$$
\begin{align*}
\ln f(\mathbf{x}_1)&=\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}+\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]+{\rm const.}\\
&\simeq\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]+\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}\left.\frac{\partial}{\partial (\mathcal{H}_1+PV_1)}\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}+\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]\right|_{\mathcal{H}_1(\mathbf{x}_1)+PV_1=0}+{\rm const.}\\
&=\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]-\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}\left.\frac{\partial}{\partial H}\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}+\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]\right|_{\mathcal{H}_1(\mathbf{x}_1)+PV_1=0}+{\rm const.}\\
&=\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]-\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}\frac{\partial}{\partial H}\ln\left[\int_{0}^{\infty}{\rm d}V_2\int_{D(V_2)}{\rm d}\mathbf{x}_2\delta\left(\left\{\mathcal{H}_2(\mathbf{x}_2)+PV_2\right\}-H\right)\right]+{\rm const.}\\
&=\ln\Gamma_2(N_2,P,H_2)-\left\{\mathcal{H}_1(\mathbf{x}_1)+PV_1\right\}\frac{\partial}{\partial H}\ln\Gamma_2(N_2,P,H_2)+{\rm const.}\\
&=\frac{S_2(N_2,P,H_2)}{k_{\rm B}}-\frac{\mathcal{H}_1(\mathbf{x}_1)+PV_1}{k_{\rm B}T}+{\rm const.}\\
\end{align*}
$$
$$
\begin{align*}
\therefore f(\mathbf{x}_1)&\propto\exp\left(-\frac{\mathcal{H}_1(\mathbf{x}_1)+PV_1}{k_{\rm B}T}\right)\\
\end{align*}
$$
Problem 5.2
等温等圧アンサンブルにおける等温圧縮率$${\kappa}$$を体積揺らぎ$${\Delta V}$$を用いて表し，$${\Delta V/\langle V\rangle\sim1/\sqrt{N}}$$となる，つまり熱力学極限において体積揺らぎは無視できることを示せ。
$$
\begin{align*}
\Delta V&=\sqrt{\left\langle V^2\right\rangle-\left\langle V\right\rangle^2}\\
\kappa&=-\frac{1}{\left\langle V\right\rangle}\left(\frac{\partial\left\langle V\right\rangle}{\partial P}\right)_{N,T}
\end{align*}
$$
$$
\begin{align*}
\left(\frac{\partial\left\langle V\right\rangle}{\partial P}\right)_{N,T}&=\frac{\partial}{\partial P}\left\{\frac{\int_0^{\infty}{\rm d}VV\int_{D(V)}\mathbf{x}{\rm e}^{-\frac{\mathcal{H}(\mathbf{x})+PV
}{k_{\rm B}T}}}{\int_0^{\infty}{\rm d}V\int_{D(V)}\mathbf{x}{\rm e}^{-\frac{\mathcal{H}(\mathbf{x})+PV
}{k_{\rm B}T}}}\right\}\\
&=-\frac{\Delta V^2}{k_{\rm B}T}\\
\therefore \kappa&=-\frac{1}{\left\langle V\right\rangle}\left(\frac{\partial\left\langle V\right\rangle}{\partial P}\right)_{N,T}\\
&=\frac{\Delta V^2}{k_{\rm B}T\left\langle V\right\rangle}
\end{align*}
$$
$${\kappa}$$は示強性のため，$${\Delta V^2\sim N}$$となる。
以上より，$${\Delta V/\langle V\rangle\sim \sqrt{N}/N=1/\sqrt{N}}$$となることが示される。
Problem 5.3
仕事に関するビリアル定理のテンソル版が式(5.6.15)で与えられることを示せ。
$$
\begin{align*}
\left\langle P_{\alpha\beta}^{\rm(int)}{\rm det}(\mathbf{h})\right\rangle+k_{\rm B}T\delta_{\alpha\beta}&=P\left\langle{\rm det}(\mathbf{h})\right\rangle\delta_{\alpha\beta}
\end{align*}\tag{5.6.15}
$$
参考文献１の式(5.6.10)の被積分関数に$${{\rm det}(\mathbf{h})}$$をかけて式展開すると，
$$
\begin{align*}
\left\langle P_{\alpha\beta}^{\rm(int)}{\rm det}(\mathbf{h})\right\rangle&=-\frac{k_{\rm B}T}{V_0\Delta(N,P,T)}\int{\rm d}\mathbf{h}\frac{\partial}{\partial h_{\alpha\gamma}}\left\{[{\rm det}(\mathbf{h})]^{-2}{\rm e}^{-\beta P {\rm det}(\mathbf{h})}\sum_{\gamma=1}^3h_{\beta\gamma}\right\}Q(N,\mathbf{h},T)\\
&=-\frac{k_{\rm B}T}{V_0\Delta(N,P,T)}\int{\rm d}\mathbf{h}\sum_{\gamma=1}^3\left\{-2[{\rm det}(\mathbf{h})]^{-3}\frac{\partial {\rm det}(\mathbf{h})}{\partial h_{\alpha\gamma}}h_{\beta\gamma}-\beta P[{\rm det}(\mathbf{h})]^{-2}\frac{\partial {\rm det}(\mathbf{h})}{\partial h_{\alpha\gamma}}h_{\beta\gamma}+[{\rm det}(\mathbf{h})]^{-2}\frac{\partial h_{\beta\gamma}}{\partial h_{\alpha\gamma}}\right\}{\rm e}^{-\beta P {\rm det}(\mathbf{h})}Q(N,\mathbf{h},T)\\
&=-\frac{k_{\rm B}T}{V_0\Delta(N,P,T)}\int{\rm d}\mathbf{h}[{\rm det}(\mathbf{h})]^{-2}\left\{-2-\beta P{\rm det}(\mathbf{h})+3\right\}\delta_{\alpha\beta}{\rm e}^{-\beta P {\rm det}(\mathbf{h})}Q(N,\mathbf{h},T)\\
&=P\left\langle{\rm det}(\mathbf{h})\right\rangle\delta_{\alpha\beta}-k_{\rm B}T\delta_{\alpha\beta}\\
\therefore \left\langle P_{\alpha\beta}^{\rm(int)}{\rm det}(\mathbf{h})\right\rangle+k_{\rm B}T\delta_{\alpha\beta}&=P\left\langle{\rm det}(\mathbf{h})\right\rangle\delta_{\alpha\beta}
\end{align*}
$$
Problem 5.4
a. 4章のProblem 4.6の理想気体に対して，シリンダーの長さのみが可変である場合の等温等圧条件における分配関数を計算せよ。
ヒント：必要に応じて二項定理を利用せよ。
b. シリンダーの長さの平均を算出せよ。
Problem 4.6の解答例については以下を参照のこと。
a. $${L_0}$$を分配関数の参照長に選ぶことにする。
このとき，等温等圧条件における分配関数$${\Delta(N,P,T)}$$は，
$$
\begin{align*}
\Delta(N,P,T)&=\frac{1}{L_0}\int_{0}^{\infty}{\rm d}L{\rm e}^{-\beta P\pi a^2 L}Q(N,V,T)\\
&=\frac{1}{L_0N!}\int_{0}^{\infty}{\rm d}L{\rm e}^{-\beta P\pi a^2 L}\left[\frac{1}{h^{3}}\int{\rm d}\mathbf{p}\int_{D(V)}{\rm d}\mathbf{r}{\rm e}^{-\beta h(\mathbf{r},\mathbf{p})}\right]^N\\
&=\frac{1}{L_0N!}\left[\frac{(2\pi)^{5/2}}{h^3\beta^{7/2}m^{1/2}w^2g}\left\{{\rm e}^{\beta m \omega^2a^2}-1\right\}\right]^N\int_{0}^{\infty}{\rm d}L{\rm e}^{-\beta P\pi a^2 L}\left(1-{\rm e}^{-\beta mgL}\right)^N\\
&=\frac{1}{L_0N!}\left[\frac{(2\pi)^{5/2}}{h^3\beta^{7/2}m^{1/2}w^2g}\left\{{\rm e}^{\beta m \omega^2a^2}-1\right\}\right]^N\sum_{n=0}^N\binom{N}{n}\int_{0}^{\infty}{\rm d}L{\rm e}^{-\beta( P\pi a^2+nmg) L}\\
&=\frac{1}{L_0N!\beta}\left[\frac{(2\pi)^{5/2}}{h^3\beta^{7/2}m^{1/2}w^2g}\left\{{\rm e}^{\beta m \omega^2a^2}-1\right\}\right]^N\sum_{n=0}^N\binom{N}{n}\frac{1}{P\pi a^2+nmg}\\
\end{align*}
$$
b.
$$
\begin{align*}
\langle L\rangle&=\frac{\sum_{n=0}^N\binom{N}{n}\int_{0}^{\infty}{\rm d}LL{\rm e}^{-\beta( P\pi a^2+nmg) L}}{\sum_{n=0}^N\binom{N}{n}\int_{0}^{\infty}{\rm d}L{\rm e}^{-\beta( P\pi a^2+nmg) L}}\\
&=\frac{\sum_{n=0}^N\binom{N}{n}(P\pi a^2+nmg)^{-2}}{\beta\sum_{n=0}^N\binom{N}{n}(P\pi a^2+nmg)^{-1}}
\end{align*}
$$
Problem 5.5
以下の各ケースに対して，節4.9の技法を用いて等方的NPT用の運動方程式である式(5.9.5)が正しいアンサンブル分布を生成することを示せ。
a. $${\sum_{i=1}^N\mathbf{F}_i\neq 0}$$
b. $${\sum_{i=1}^N\mathbf{F}_i= 0}$$，この場合は以下の量も保存することになる
$$
\begin{align*}
\mathbf{K}&=\mathbf{P}\exp\left[\left(1+\frac{1}{N}\right)\epsilon+\eta_1\right]
\end{align*}
$$
ここで，$${\mathbf{P}=\sum_{i=1}^N\mathbf{p}_i}$$は重心運動量である。
c. "Massive"熱浴を課した場合
$$
\begin{align*}
\dot{\mathbf{r}}_i&=\frac{\dot{\mathbf{p}}_i}{m_i}+\frac{p_{\epsilon}}{W}\mathbf{r}_i\\
\dot{\mathbf{p}}_i&=\mathbf{F}_i-\left(1+\frac{1}{N}\right)\frac{p_{\epsilon}}{W}\mathbf{p}_i-\frac{p_{\eta_1}}{Q_1}\mathbf{p}_i\\
\dot{V}&=\frac{dVp_{\epsilon}}{W}\\
\dot{p}_{\epsilon}&=dV(\mathcal{P}^{\rm (int)}-P)+\frac{1}{N}\sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i}-\frac{p_{\xi_1}}{Q_1'}p_{\epsilon}\\
\dot{\eta}_j&=\frac{p_{\eta_j}}{Q_j}\\
\dot{\xi}_j&=\frac{p_{\xi_j}}{Q_j'}\\
\dot{p}_{\eta_j}&=G_j-\frac{p_{\eta_{j+1}}}{Q_{j+1}}p_{\eta_j}\\
\dot{p}_{\eta_M}&=G_M\\
\dot{p}_{\xi_j}&=G_j'-\frac{p_{\xi_{j+1}}}{Q_{j+1}'}p_{\xi_j}\\
\dot{p}_{\xi_M}&=G_M'
\end{align*}\tag{5.9.5}
$$
a. 位相空間ベクトルを
$$
\begin{align*}
\mathbf{x}&=\begin{bmatrix}\mathbf{r}_1^{\rm T}&\cdots&\mathbf{r}_N^{\rm T}&\mathbf{p}_1^{\rm T}&\cdots&\mathbf{p}_N^{\rm T}&\boldsymbol\eta^{\rm T}&\mathbf{p}_{\eta}^{\rm T}&\boldsymbol\xi^{\rm T}&\mathbf{p}_{\xi}^{\rm T}&V&p_{\epsilon}\end{bmatrix}^{\rm T}
\end{align*}
$$
と定義する。
このとき，
$$
\begin{align*}
\nabla_{\mathbf{x}}\cdot\dot{\mathbf{x}}&=\sum_{i=1}^N\left\{\frac{\partial \dot{\mathbf{r}}_i}{\partial \mathbf{r}_i}+\frac{\partial \dot{\mathbf{p}}_i}{\partial \mathbf{p}_i}\right\}+\sum_{j=1}^M\left\{\frac{\partial \dot{\eta}_j}{\partial \eta_j}+\frac{\partial \dot{p}_{\eta_j}}{\partial p_{\eta_j}}+\frac{\partial \dot{\xi}_j}{\partial \xi_j}+\frac{\partial \dot{p}_{\xi_j}}{\partial p_{\xi_j}}\right\}++\frac{\partial \dot{V}}{\partial V}++\frac{\partial \dot{p}_{\epsilon}}{\partial p_{\epsilon}}\\
&=dN\left\{\frac{p_{\epsilon}}{W}-\left(1+\frac{1}{N}\right)\frac{p_{\epsilon}}{W}-\frac{p_{\eta_1}}{Q_1}\right\}-\sum_{j=1}^{M-1}\left\{\frac{p_{\eta_{j+1}}}{Q_{j+1}}+\frac{p_{\xi_{j+1}}}{Q_{j+1}'}\right\}+\frac{dp_{\epsilon}}{W}-\frac{p_{\xi_1}}{Q_1'}\\
&=-dN\frac{p_{\eta_1}}{Q_1}-\sum_{j=2}^{M}\frac{p_{\eta_{j}}}{Q_{j}}-\sum_{j=1}^{M}\frac{p_{\xi_{j}}}{Q_{j}'}\\
&=-\frac{{\rm d}}{{\rm d}t}\left\{dN\eta_1+\sum_{j=2}^{M}\eta_{j}+\sum_{j=1}^{M}\xi_{j}\right\}\\
\therefore \sqrt{g}&=\exp\left\{dN\eta_1+\sum_{j=2}^{M}\eta_{j}+\sum_{j=1}^{M}\xi_{j}\right\}\\
\end{align*}
$$
保存量が参考文献1の式(5.9.7)で与えられる$${\mathcal{H}'}$$のみであるため，分配関数$${\mathcal{Z}}$$は
$$
\begin{align*}
\mathcal{Z}&=\int{\rm d}\mathbf{x}\exp\left\{dN\eta_1+\sum_{j=2}^{M}\eta_{j}+\sum_{j=1}^{M}\xi_{j}\right\}\delta\left(\mathcal{H}(\mathbf{r},\mathbf{p})+\frac{p_{\epsilon}^2}{2W}+PV+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}+k_{\rm B}T\xi_j\right]+dNk_{\rm B}T\eta_1+k_{\rm B}T\sum_{j=2}^M\eta_j-\mathcal{H}'\right)\\
&=\frac{1}{dNk_{\rm B}T}\int{\rm d}V\int{\rm d}^N\mathbf{r}\int{\rm d}^N\mathbf{p}\int{\rm d}\boldsymbol\eta_{-1}\int{\rm d}\mathbf{p}_{\eta}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\int{\rm d}p_{\epsilon}\exp\left\{dN\frac{\mathcal{H}'-\mathcal{H}(\mathbf{r},\mathbf{p})-\frac{p_{\epsilon}^2}{2W}-PV-\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}+k_{\rm B}T\xi_j\right]-k_{\rm B}T\sum_{j=2}^M\eta_j}{dNk_{\rm B}T}+\sum_{j=2}^{M}\eta_{j}+\sum_{j=1}^{M}\xi_{j}\right\}\\
&=\frac{1}{dNk_{\rm B}T}\int{\rm d}V\int{\rm d}^N\mathbf{r}\int{\rm d}^N\mathbf{p}\int{\rm d}\boldsymbol\eta_{-1}\int{\rm d}\mathbf{p}_{\eta}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\int{\rm d}p_{\epsilon}\exp\left\{\frac{\mathcal{H}'-\mathcal{H}(\mathbf{r},\mathbf{p})-\frac{p_{\epsilon}^2}{2W}-PV-\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]}{k_{\rm B}T}\right\}\\
&=\frac{\exp\left(\frac{\mathcal{
H}'}{k_{\rm B}T}\right)}{dNk_{\rm B}T}\int{\rm d}V\int{\rm d}\boldsymbol\eta_{-1}\int{\rm d}\mathbf{p}_{\eta}\exp\left\{-\frac{\sum_{j=1}^M\frac{p_{\eta_j}^2}{2Q_j}}{k_{\rm B}T}\right\}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\exp\left\{-\frac{\sum_{j=1}^M\frac{p_{\xi_j}^2}{2Q_j}}{k_{\rm B}T}\right\}\int{\rm d}p_{\epsilon}\exp\left\{-\frac{p_{\epsilon}^2}{2Wk_{\rm B}T}\right\}\int{\rm d}^N\mathbf{r}\int{\rm d}^N\mathbf{p}\exp\left\{-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right\}\\
&\propto\int{\rm d}V\int{\rm d}^N\mathbf{r}\int{\rm d}^N\mathbf{p}\exp\left\{-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right\}\\
\end{align*}
$$
となる。
b. Problem 4.3 (g)と同じ変数変換を考える。
２つの保存量があることを考慮すると，分配関数$${\mathcal{Z}}$$は
$$
\begin{align*}
\mathcal{Z}&=\int{\rm d}V\int{\rm d}^{N-1}\mathbf{r}'\int{\rm d}^{N-1}\mathbf{p}'\int{\rm d}\mathbf{R}\int{\rm d}\mathbf{P}\int{\rm d}\boldsymbol\eta\int{\rm d}\mathbf{p}_{\eta}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\int{\rm d}p_{\epsilon}\exp\left\{dN\eta_1+\sum_{j=2}^{M}\eta_{j}+\sum_{j=1}^{M}\xi_{j}\right\}\delta\left(\mathcal{H}(\mathbf{r}',\mathbf{p}')+\frac{1}{2M}\mathbf{P}^2+\frac{p_{\epsilon}^2}{2W}+PV+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}+k_{\rm B}T\xi_j\right]+dNk_{\rm B}T\eta_1+k_{\rm B}T\sum_{j=2}^M\eta_j-\mathcal{H}'\right)\delta\left(\mathbf{P}\exp\left[\left(1+\frac{1}{N}\right)\epsilon+\eta_1\right]-\mathbf{K}\right)\\
&=\int{\rm d}V\int{\rm d}^{N}\mathbf{r}\int{\rm d}^{N-1}\mathbf{p}'\int{\rm d}\mathbf{P}\int{\rm d}\boldsymbol\eta\int{\rm d}\mathbf{p}_{\eta}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\int{\rm d}p_{\epsilon}\exp\left\{dN\eta_1+\sum_{j=2}^{M}\eta_{j}+\sum_{j=1}^{M}\xi_{j}\right\}\delta\left(\sum_{i=1}^{N-1}\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{1}{2M}\mathbf{P}^2+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\frac{p_{\epsilon}^2}{2W}+PV+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}+k_{\rm B}T\xi_j\right]+dNk_{\rm B}T\eta_1+k_{\rm B}T\sum_{j=2}^M\eta_j-\mathcal{H}'\right)\delta\left(\mathbf{P}\exp\left[\left(1+\frac{1}{N}\right)\epsilon+\eta_1\right]-\mathbf{K}\right)\\
&=\frac{{\rm e}^{\frac{\mathcal{H}'}{k_{\rm B}T}}}{k_{\rm B}T}\int{\rm d}V\int{\rm d}^{N}\mathbf{r}\int{\rm d}^{N-1}\mathbf{p}'\int{\rm d}\mathbf{P}\int{\rm d}\boldsymbol\eta_{-M}\int{\rm d}\mathbf{p}_{\eta}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\int{\rm d}p_{\epsilon}\exp\left\{-\frac{\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{1}{2M}\mathbf{P}^2+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\frac{p_{\epsilon}^2}{2W}+PV+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]}{k_{\rm B}T}\right\}\delta\left(\mathbf{P}\exp\left[\left(1+\frac{1}{N}\right)\epsilon+\eta_1\right]-\mathbf{K}\right)\\
&=\frac{{\rm e}^{\frac{\mathcal{H}'}{k_{\rm B}T}}}{k_{\rm B}T}\int{\rm d}V\int{\rm d}^N\mathbf{r}\int{\rm d}^{N-1}\mathbf{p}'\int{\rm d}\boldsymbol\eta_{-M}\int{\rm d}\mathbf{p}_{\eta}\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}\int{\rm d}p_{\epsilon}\exp\left\{-\frac{\frac{(\mathbf{p}_i')^2}{2\mu_i}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\frac{\mathbf{K}^2}{2M}\exp\left\{-2\left(1+\frac{1}{N}\right)\epsilon-2\eta_1\right\}+\frac{p_{\epsilon}^2}{2W}+PV+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]}{k_{\rm B}T}\right\}\exp\left\{-d\left(1+\frac{1}{N}\right)\epsilon-d\eta_1\right\}\\
&=\frac{{\rm e}^{\frac{\mathcal{H}'}{k_{\rm B}T}}}{k_{\rm B}T}\int{\rm d}V\int{\rm d}^N\mathbf{r}{\rm e}^{-\frac{U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+PV}{k_{\rm B}T}}\int{\rm d}^{N-1}\mathbf{p}'{\rm e}^{-\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}}{k_{\rm B}T}}\int{\rm d}\eta_1\exp\left\{-\frac{\mathbf{K}^2\exp\left\{-2\left(1+\frac{1}{N}\right)\epsilon-2\eta_1\right\}}{2Mk_{\rm B}T}\right\}\exp\left\{-d\left(1+\frac{1}{N}\right)\epsilon-d\eta_1\right\}\left(\prod_{j=2}^{M-1}\int{\rm d}\eta_j\right)\int{\rm d}\boldsymbol\xi\int{\rm d}\mathbf{p}_{\xi}{\rm e}^{-\frac{\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]}{k_{\rm B}T}}\int{\rm d}p_{\epsilon}{\rm e}^{-\frac{p_{\epsilon}^2}{2Wk_{\rm B}T}}
\end{align*}
$$
となる。
ここで，$${\eta_1}$$の積分に着目する。
$${\exp\{-(1+1/N)\epsilon-\eta_1\}=x}$$と変数変換すると，
$$
\begin{align*}
\int_{-\infty}^{\infty}{\rm d}\eta_1\exp\left\{-\frac{\mathbf{K}^2\exp\left\{-2\left(1+\frac{1}{N}\right)\epsilon-2\eta_1\right\}}{2Mk_{\rm B}T}\right\}\exp\left\{-d\left(1+\frac{1}{N}\right)\epsilon-d\eta_1\right\}&=\int_{0}^{\infty}{\rm d}xx^{d-1}\exp\left(-\frac{\mathbf{K}^2x^2}{2Mk_{\rm B}T}\right)
\end{align*}
$$
となり，積分結果は$${\epsilon}$$に依存しない（＝$${V}$$に依存しない）。
以上より，
$$
\begin{align*}
\mathcal{Z}&\propto\int{\rm d}V\int{\rm d}^N\mathbf{r}{\rm e}^{-\frac{U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+PV}{k_{\rm B}T}}\int{\rm d}^{N-1}\mathbf{p}'{\rm e}^{-\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}}{k_{\rm B}T}}
\end{align*}
$$
となり，配置と体積について正しいアンサンブルが得られ，また相対運動量がMaxwell速度分布に従う。
c. "Massive"熱浴の場合，重心運動量に関する保存則が破綻するため，運動量も含めて正しいアンサンブルが得られる。
Problem 5.6
$${\sum_{i=1}^N\mathbf{F}_i=\mathbf{0}}$$の場合，式(5.10.2)の運動方程式が正しいアンサンブル分布を生成することを示せ。
$$
\begin{align*}
\dot{\mathbf{r}}_i&=\frac{\mathbf{p}_i}{m_i}+\frac{\mathbf{p}_g}{W_g}\mathbf{r}_i\\
\dot{\mathbf{p}}_i&=\mathbf{F}_i-\frac{\mathbf{p}_g}{W_g}\mathbf{p}_i-\frac{1}{dN}\frac{{\rm Tr}[\mathbf{p}_g]}{W_g}\mathbf{p}_i-\frac{p_{\eta_1}}{Q_1}\mathbf{p}_i\\
\dot{\mathbf{h}}&=\frac{\mathbf{p}_g\mathbf{h}}{W_g}\\
\dot{\mathbf{p}}_g&={\rm det}[\mathbf{h}]\left(\mathbf{P}^{\rm(int)}-\mathbf{I}P\right)+\frac{1}{dN}\sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i}\mathbf{I}-\frac{p_{\xi_1}}{Q_1'}\mathbf{p}_g\\
\dot{\eta}_j&=\frac{p_{\eta_j}}{Q_j}\\
\dot{\xi}_j&=\frac{p_{\xi_j}}{Q_j'}\\
\dot{p}_{\eta_M}&=G_M\\
\dot{p}_{\xi_M}&=G_j'-\frac{p_{\xi_{j+1}}}{Q_{j+1}'}p_{\xi_j}\\
\dot{p}_{\xi_M}&=G_M'\\
\end{align*}\tag{5.10.2}
$$
Problem 4.3 (g)と同じ変数変換を考える。
２つの保存量があることを考慮すると，分配関数$${\mathcal{Z}}$$は
$$
\begin{align*}
\mathcal{Z}&=\int{\rm d}^{N-1}\mathbf{p}'{\rm d}\mathbf{P}{\rm d}^{N}\mathbf{r}{\rm d}\mathbf{h}{\rm d}\mathbf{p}_g{\rm d}\eta_1{\rm d}\eta_c{\rm d}\xi_1{\rm d}\xi_c{\rm d}^Mp_{\eta}{\rm d}^Mp_{\xi}[{\rm det}(\mathbf{h})]^{1-d}{\rm e}^{dN\eta_1+\eta_c}{\rm e}^{d^2\xi_1+\xi_c}\\
&\times\delta\left(\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\mathbf{P}^2}{2M}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]+dNk_{\rm B}T\eta_1+d^2k_{\rm B}T\xi_1+k_{\rm B}T[\eta_c+\xi_c]+\frac{{\rm Tr}\left[\mathbf{p}_g\mathbf{p}_g^{\rm T}\right]}{2W_g}+P{\rm det}[\mathbf{h}]-H\right)\\
&\times\delta\left(\mathbf{K}-\mathbf{h}\mathbf{P}\{{\rm det}[\mathbf{h}]\}^{\frac{1}{dN}}{\rm e}^{\eta_1}\right)\\
&=\int{\rm d}^{N-1}\mathbf{p}'{\rm d}\mathbf{P}'{\rm d}^{N}\mathbf{r}{\rm d}\mathbf{h}{\rm d}\mathbf{p}_g{\rm d}\eta_1{\rm d}\eta_c{\rm d}\xi_1{\rm d}\xi_c{\rm d}^Mp_{\eta}{\rm d}^Mp_{\xi}[{\rm det}(\mathbf{h})]^{-d}{\rm e}^{dN\eta_1+\eta_c}{\rm e}^{d^2\xi_1+\xi_c}\\
&\times\delta\left(\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\left(\mathbf{h}^{-1}\mathbf{P}'\right)^2}{2M}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]+dNk_{\rm B}T\eta_1+d^2k_{\rm B}T\xi_1+k_{\rm B}T[\eta_c+\xi_c]+\frac{{\rm Tr}\left[\mathbf{p}_g\mathbf{p}_g^{\rm T}\right]}{2W_g}+P{\rm det}[\mathbf{h}]-H\right)\\
&\times\delta\left(\mathbf{K}-\mathbf{P}'\{{\rm det}[\mathbf{h}]\}^{\frac{1}{dN}}{\rm e}^{\eta_1}\right)\\
&=\int{\rm d}^{N-1}\mathbf{p}'{\rm d}\mathbf{P}'{\rm d}^{N}\mathbf{r}{\rm d}\mathbf{h}{\rm d}\mathbf{p}_g{\rm d}\eta_1{\rm d}\eta_c{\rm d}\xi_1{\rm d}^Mp_{\eta}{\rm d}^Mp_{\xi}[{\rm det}(\mathbf{h})]^{-d}\\
&\times\exp\left\{-\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\left(\mathbf{h}^{-1}\mathbf{P}'\right)^2}{2M}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]+\frac{{\rm Tr}\left[\mathbf{p}_g\mathbf{p}_g^{\rm T}\right]}{2W_g}+P{\rm det}[\mathbf{h}]-H}{k_{\rm B}T}\right\}\\
&\times\delta\left(\mathbf{K}-\mathbf{P}'\{{\rm det}[\mathbf{h}]\}^{\frac{1}{dN}}{\rm e}^{\eta_1}\right)\\
&=\int{\rm d}^{N-1}\mathbf{p}'{\rm d}^{N}\mathbf{r}{\rm d}\mathbf{h}{\rm d}\mathbf{p}_g{\rm d}\eta_1{\rm d}\eta_c{\rm d}\xi_1{\rm d}^Mp_{\eta}{\rm d}^Mp_{\xi}[{\rm det}(\mathbf{h})]^{-d-\frac{1}{N}}{\rm e}^{-d\eta_1}\\
&\times\exp\left\{-\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\{{\rm det}[\mathbf{h}]\}^{-\frac{2}{dN}}{\rm e}^{-2\eta_1}\left(\mathbf{h}^{-1}\mathbf{K}\right)^2}{2M}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\sum_{j=1}^M\left[\frac{p_{\eta_j}^2}{2Q_j}+\frac{p_{\xi_j}^2}{2Q_j'}\right]+\frac{{\rm Tr}\left[\mathbf{p}_g\mathbf{p}_g^{\rm T}\right]}{2W_g}+P{\rm det}[\mathbf{h}]-H}{k_{\rm B}T}\right\}\\
&\propto\int{\rm d}^{N-1}\mathbf{p}'{\rm d}^{N}\mathbf{r}{\rm d}\mathbf{h}{\rm d}\eta_1[{\rm det}(\mathbf{h})]^{-d-\frac{1}{N}}{\rm e}^{-d\eta_1}\\
&\times\exp\left\{-\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\{{\rm det}[\mathbf{h}]\}^{-\frac{2}{dN}}{\rm e}^{-2\eta_1}\left(\mathbf{h}^{-1}\mathbf{K}\right)^2}{2M}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+P{\rm det}[\mathbf{h}]}{k_{\rm B}T}\right\}\\
\end{align*}
$$
となる。
ここで，$${\{{\rm det}[\mathbf{h}]\}^{-\frac{1}{dN}}{\rm e}^{-\eta_1}(\mathbf{h}^{-1}\mathbf{K})=x}$$とおくと，$${\eta_1}$$に関する積分は
$$
\begin{align*}
\int_{-\infty}^{\infty}{\rm d}\eta_1[{\rm det}(\mathbf{h})]^{-d-\frac{1}{N}}{\rm e}^{-d\eta_1}\exp\left\{-\frac{\{{\rm det}[\mathbf{h}]\}^{-\frac{2}{dN}}{\rm e}^{-2\eta_1}\left(\mathbf{h}^{-1}\mathbf{K}\right)^2}{2Mk_{\rm B}T}\right\}&=[{\rm det}(\mathbf{h})]^{-d}\left(\frac{1}{\mathbf{h}^{-1}\mathbf{K}}\right)^d\int_{0}^{\infty}{\rm d}xx^{d-1}\exp\left(-\frac{x^2}{2Mk_{\rm B}T}\right)\\
&\propto\left[{\rm det}(\mathbf{h})\mathbf{h}^{-1}\mathbf{K}\right]^{-d}
\end{align*}
$$
$$
\begin{align*}
\therefore \mathcal{Z}
&\propto\int{\rm d}^{N-1}\mathbf{p}'{\rm d}^{N}\mathbf{r}{\rm d}\mathbf{h}\left[{\rm det}(\mathbf{h})\mathbf{h}^{-1}\mathbf{K}\right]^{-d}\exp\left\{-\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+P{\rm det}[\mathbf{h}]}{k_{\rm B}T}\right\}\\
\end{align*}
$$
Problem 5.7
以下の運動方程式に従う分子動力学によって等方的NPTアンサンブルを生成するアルゴリズムが1985年にHooverによって提案された，
$$
\begin{align*}
\dot{\mathbf{r}}_i&=\frac{\mathbf{p}_i}{m_i}+\frac{p_{\epsilon}}{W}\mathbf{r}_i\\
\dot{\mathbf{p}}_i&=-\frac{\partial U}{\partial\mathbf{r}_i}-\frac{p_{\epsilon}}{W}\mathbf{p}_i-\frac{p_{\eta}}{Q}\mathbf{p}_i\\
\dot{V}&=\frac{dVp_{\epsilon}}{W}\\
\dot{p}_{\epsilon}&=dV(\mathcal{P}^{\rm (int)}-P)-\frac{p_{\eta}}{Q}p_{\epsilon}\\
\dot{\eta}&=\frac{p_{\eta}}{Q}\\
\dot{p}_{\eta}&=\sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i}+\frac{p_{\epsilon}^2}{W}-(dN+1)k_{\rm B}T
\end{align*}
$$
ここで，$${\mathcal{P}^{\rm (int)}}$$は式(5.7.28)（等方性なので，式(4.6.57)に等しい）で与えられる圧力の瞬時値である。
$$
\begin{align*}
\mathcal{P}^{\rm (int)}&=\mathcal{P}(\mathbf{r},\mathbf{p})\\
&=\frac{1}{3V}\sum_{i=1}^N\left[\frac{\mathbf{p}_i}{m_i}+\mathbf{r}_i\cdot\mathbf{F}_i\right]
\end{align*}\tag{4.6.57}
$$
これらの運動方程式では以下のエネルギーが保存される。
$$
\begin{align*}
\mathcal{H}'&=\mathcal{H}(\mathbf{r},\mathbf{p})+\frac{p_{\epsilon}^2}{2W}+\frac{p_{\eta}^2}{2Q}+(dN+1)k_{\rm B}T\eta+PV
\end{align*}
$$
以下のそれぞれの場合におけるアンサンブルの分布関数$${f(\mathbf{r},\mathbf{p},V)}$$を決定せよ。
a. $${\sum_{i=1}^N\mathbf{F}_i\neq\mathbf{0}}$$の場合。分布関数は，熱力学極限において正しい等温押圧アンサンブルの分布関数に収束するであろうか？
b. $${\sum_{i=1}^N\mathbf{F}_i=\mathbf{0}}$$の場合。この場合には以下の保存量が追加で存在する，
$$
\begin{align*}
\mathbf{K}&=\mathbf{P}{\rm e}^{\epsilon+\eta}
\end{align*}
$$
ここで，$${\mathbf{P}}$$は重心運動量である。全ての非物理的な変数に対する積分を確実に実行せよ。
a.
$$
\begin{align*}
\nabla_{\mathbf{x}}\cdot\dot{\mathbf{x}}&=\frac{dp_{\epsilon}}{W}-\frac{(dN+1)p_{\eta}}{Q}\\
&=-\frac{\rm d}{{\rm d}t}\left[-\ln V + (dN+1)\eta\right]\\
&=\dot{w}\\
\therefore \sqrt{g}&={\rm e}^{-w}\\
&={\rm e}^{-\ln V + (dN+1)\eta}\\
&=V^{-1}{\rm e}^{(dN+1)\eta}\\
\end{align*}
$$
保存量が$${\mathcal{H}'}$$のみの場合，分配関数$${\mathcal{Z}}$$は
$$
\begin{align*}
\mathcal{Z}&=\int{\rm d}^N\mathbf{r}{\rm d}^N\mathbf{p}{\rm d}V{\rm d}\eta{\rm d}p_{\eta}{\rm d}p_{\epsilon}V^{-1}{\rm e}^{(dN+1)\eta}\delta\left(\mathcal{H}(\mathbf{r},\mathbf{p})+\frac{p_{\epsilon}^2}{2W}+\frac{p_{\eta}^2}{2Q}+(dN+1)k_{\rm B}T\eta+PV-\mathcal{H}'\right)\\
&=\int{\rm d}^N\mathbf{r}{\rm d}^N\mathbf{p}{\rm d}V{\rm d}p_{\eta}{\rm d}p_{\epsilon}V^{-1}\exp\left(-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+\frac{p_{\epsilon}^2}{2W}+\frac{p_{\eta}^2}{2Q}+PV-\mathcal{H}'}{k_{\rm B}T}\right)\\
&\propto\int{\rm d}^N\mathbf{r}{\rm d}^N\mathbf{p}{\rm d}VV^{-1}\exp\left(-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right)\\
\end{align*}
$$
熱力学極限において，$${\Delta V\rightarrow 0}$$となるため，
$$
\begin{align*}
\int{\rm d}^N\mathbf{r}{\rm d}^N\mathbf{p}{\rm d}VV^{-1}\exp\left(-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right)&\simeq\frac{1}{\langle V\rangle}\int{\rm d}^N\mathbf{r}{\rm d}^N\mathbf{p}{\rm d}V\exp\left(-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right)\\
&\propto\int{\rm d}^N\mathbf{r}{\rm d}^N\mathbf{p}{\rm d}V\exp\left(-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right)
\end{align*}
$$
$$
\begin{align*}
\therefore f(\mathbf{r},\mathbf{p},V)&\propto \exp\left(-\frac{\mathcal{H}(\mathbf{r},\mathbf{p})+PV}{k_{\rm B}T}\right)\ \ \ \ ({\rm for}\ N\gg 1)
\end{align*}
$$
b. 保存量が$${\mathcal{H}', \mathbf{K}}$$の2つある場合，分配関数$${\mathcal{Z}}$$は
$$
\begin{align*}
\mathcal{Z}&=\int{\rm d}^N\mathbf{r}{\rm d}^{N-1}\mathbf{p}'{\rm d}\mathbf{P}{\rm d}V{\rm d}\eta{\rm d}p_{\eta}{\rm d}p_{\epsilon}V^{-1}{\rm e}^{(dN+1)\eta}\delta\left(\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\mathbf{P}^2}{2M}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\frac{p_{\epsilon}^2}{2W}+\frac{p_{\eta}^2}{2Q}+(dN+1)k_{\rm B}T\eta+PV-\mathcal{H}'\right)\delta\left(\mathbf{P}{\rm e}^{\epsilon+\eta}-\mathbf{K}\right)\\
&=\int{\rm d}^N\mathbf{r}{\rm d}^{N-1}\mathbf{p}'{\rm d}V{\rm d}\eta{\rm d}p_{\eta}{\rm d}p_{\epsilon}V^{-1}{\rm e}^{(dN+1-d)\eta}{\rm e}^{-d\epsilon}\delta\left(\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\mathbf{K}^2}{2M}{\rm e}^{-2(\epsilon+\eta)}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+\frac{p_{\epsilon}^2}{2W}+\frac{p_{\eta}^2}{2Q}+(dN+1)k_{\rm B}T\eta+PV-\mathcal{H}'\right)\\
&=4\pi\sqrt{WQ}V_0\int{\rm d}^N\mathbf{r}{\rm d}^{N-1}\mathbf{p}'{\rm d}V{\rm d}\eta{\rm d}p_{r}p_{r}V^{-2}{\rm e}^{(dN+1-d)\eta}\delta\left(\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\mathbf{K}^2}{2M}\left(\frac{V}{V_0}\right)^{-2/d}{\rm e}^{-2\eta}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+p_r^2+(dN+1)k_{\rm B}T\eta+PV-\mathcal{H}'\right)\\
&=2\pi\sqrt{WQ}V_0\int{\rm d}^N\mathbf{r}{\rm d}^{N-1}\mathbf{p}'{\rm d}V{\rm d}\eta{\rm d}yV^{-2}{\rm e}^{(dN+1-d)\eta}\delta\left(\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+\frac{\mathbf{K}^2}{2M}\left(\frac{V}{V_0}\right)^{-2/d}{\rm e}^{-2\eta}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+y+(dN+1)k_{\rm B}T\eta+PV-\mathcal{H}'\right)\\
&=\frac{2\pi\sqrt{WQ}V_0}{(dN+1)k_{\rm B}T}\int{\rm d}^N\mathbf{r}{\rm d}^{N-1}\mathbf{p}'{\rm d}V{\rm d}yV^{-2}\left\{W\left(-\frac{\mathbf{K}^2V_0^{1/d}\exp\left(2\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+y+PV-\mathcal{H}'}{(dN+1)k_{\rm B}T}\right)}{2M(dN+1)k_{\rm B}TV^{1/d}}\right)+1\right\}^{-1}\left\{-\frac{MV^{1/d}(dN+1)k_{\rm B}T}{\mathbf{K}^2V^{1/d}}W\left(-\frac{\mathbf{K}^2V_0^{1/d}\exp\left(2\frac{\sum_{i=2}^N\frac{(\mathbf{p}_i')^2}{2\mu_i}+U(\mathbf{r}_1,\cdots,\mathbf{r}_N)+y+PV-\mathcal{H}'}{(dN+1)k_{\rm B}T}\right)}{2M(dN+1)k_{\rm B}TV^{1/d}}\right)\right\}^{-\frac{dN+1-d}{2}}\\
\end{align*}
$$
Problem 5.8
圧力を分子ビリアルを用いて決定した場合において，式(5.11.6)が正しい等圧等エンタルピーアンサンブルの分布を生成することを示せ。
$$
\begin{align*}
\dot{\mathbf{r}}_{i,\alpha}&=\frac{\mathbf{p}_{i,\alpha}}{m_{i,\alpha}}+\frac{p_{\epsilon}}{W}\mathbf{R}_i\\
\dot{\mathbf{p}}_{i,\alpha}&=\mathbf{F}_{i,\alpha}-\left(1+\frac{1}{N}\right)\frac{p_{\epsilon}}{W}\frac{m_{i,\alpha}}{M_i}\mathbf{P}_i\\
\dot{V}&=\frac{dVp_{\epsilon}}{W}\\
\dot{p}_{\epsilon}&=dV(\mathcal{P}_{\rm mol}-P)+\frac{1}{N}\sum_{i=1}\frac{\mathbf{P}_i^2}{M_i}\\
\end{align*}\tag{5.11.6}
$$
分子圧力の瞬時値$${\mathcal{P}_{\rm mol}(\mathbf{R},\mathbf{P})}$$は以下の式で与えられる。
$$
\begin{align*}
\mathcal{P}_{\rm mol}(\mathbf{R},\mathbf{P})&=\frac{1}{dV}\sum_{i=1}^N\left[\frac{\mathbf{P}_i^2}{M_i}+\mathbf{R}_i\cdot\left(\sum_{\alpha=1}^n\mathbf{F}_{i,\alpha}\right)\right]
\end{align*}
$$
また，以下のエネルギー$${\mathcal{H}'}$$が保存する。
$$
\begin{align*}
\mathcal{H}'&=\sum_{i,\alpha}\frac{\mathbf{p}_{i,\alpha}^2}{2m_{i,\alpha}}+U(\mathbf{r})+\frac{p_{\epsilon}^2}{2W}+PV
\end{align*}
$$
$$
\begin{align*}
\nabla_{\mathbf{x}}\cdot\dot{\mathbf{x}}&=0
\end{align*}
$$
より，分配関数$${\mathcal{Z}}$$は
$$
\begin{align*}
\mathcal{Z}&\propto\int{\rm d}V{\rm d}^{nN}\mathbf{r}{\rm d}^{nN}\mathbf{p}{\rm d}p_{\epsilon}\delta\left(\sum_{i,\alpha}\frac{\mathbf{p}_{i,\alpha}^2}{2m_{i,\alpha}}+U(\mathbf{r})+\frac{p_{\epsilon}^2}{2W}+PV-\mathcal{H}'\right)
\end{align*}
$$
となり，分配関数は式(5.8.7)と同様の表式に帰着する。
式(5.8.7)と同様に真のエンタルピーと$${\frac{p_{\epsilon}^2}{2W}}$$だけ異なるが，$${N\gg 1}$$の状況では$${\frac{p_{\epsilon}^2}{2W}}$$の寄与は無視できる程度に小さい。
Problem 5.9
ナノワイヤを通る粒子の運動の単純モデルとして，周期ポテンシャル中を運動するN粒子の１次元理想気体を考える。座標$${q}$$，運動量$${p}$$としときの１粒子のハミルトニアンは以下のようになる。
$$
\begin{align*}
h(q,p)&=\frac{p^2}{2m}+\frac{kL^2}{4\pi^2}\left[1-\cos\left(\frac{2\pi q}{L}\right)\right]
\end{align*}
$$
ここで，$${m}$$は粒子の質量，$${k}$$は定数，$${K}$$は１次元箱（もしくは単位胞）の長さである。
a. 単位胞の長さを$${L_1}$$から$${L_2}$$に変化させた際の1粒子あたりのヘルムホルツ自由エネルギー変化を計算せよ。0次の変形ベッセル関数を用いて解答せよ。
$$
\begin{align*}
I_0(x)&=\frac{1}{\pi}\int_0^{\pi}{\rm d}\theta{\rm e}^{\pm x\cos\theta}
\end{align*}
$$
b. 一次元圧力"$${P}$$"を求めることによって状態方程式を計算せよ。理想気体の状態方程式が得られるであろうか？得られた状態方程式について考察せよ。必要に応じて，変形ベッセル関数の以下の性質を利用せよ。
$$
\begin{align*}
\frac{{\rm d}I_{\nu}(x)}{{\rm d}x}&=\frac{1}{2}\left[I_{\nu+1}(x)+I_{\nu-1}(x)\right]\\
I_{\nu}(x)&=I_{-\nu}(x)
\end{align*}
$$
c. 等温等圧アンサンブルにおける位置と長さの分布関数を積分表示で書き下せ。
a. 1粒子あたりのNVTアンサンブルの分配関数を$${q(L,T)}$$とおくと，
$$
\begin{align*}
q(L,T)&=\frac{1}{h}\int_{0}^{L}{\rm d}q\int_{-\infty}^{\infty}{\rm d}p\exp\left(-\frac{h(q,p)}{k_{\rm B}T}\right)\\
&=\frac{1}{h}\int_{0}^{L}{\rm d}q\exp\left(-\frac{kL^2}{4\pi^2 k_{\rm B}T}\left[1-\cos\left(\frac{2\pi q}{L}\right)\right]\right)\int_{-\infty}^{\infty}{\rm d}p\exp\left(-\frac{p^2}{2mk_{\rm B}T}\right)\\
&=\frac{\sqrt{2\pi mk_{\rm B}T}}{h}\exp\left(-\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)\int_{0}^{L}{\rm d}q\exp\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\cos\left(\frac{2\pi q}{L}\right)\right)\\
&=\frac{L\sqrt{2\pi mk_{\rm B}T}}{2\pi h}\exp\left(-\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)\int_{0}^{2\pi}{\rm d}\theta\exp\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\cos\theta\right)\\
&=\frac{L\sqrt{2\pi mk_{\rm B}T}}{h}\exp\left(-\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)I_{0}\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)\\
\therefore Q(N,L,T)&=\frac{1}{N!}\left[q(L,T)\right]^N\\
&=\frac{1}{N!}\left[\frac{L\sqrt{2\pi mk_{\rm B}T}}{ h}\exp\left(-\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)I_{0}\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)\right]^N\\
\end{align*}
$$
$${Q(N,L,T)}$$を用いて，$${L_1\rightarrow L_2}$$の自由エネルギー変化$${\Delta F_{12}}$$を計算すると，
$$
\begin{align*}
\Delta F_{12}&=F(N,L_2,T)-F(N,L_1,T)\\
&=-k_{\rm B}T\ln\left\{\frac{Q(N,L_2,T)}{Q(N,L_1,T)}\right\}\\
&=N\left[k_{\rm B}T\ln\left\{\frac{L_1I_{0}\left(\frac{kL_1^2}{4\pi^2 k_{\rm B}T}\right)}{L_2I_{0}\left(\frac{kL_2^2}{4\pi^2 k_{\rm B}T}\right)}\right\}+\frac{k}{4\pi^2}\left(L_2^2-L_1^2\right)\right]
\end{align*}
$$
b.
$$
\begin{align*}
P&=-\left(\frac{\partial F(N,L,T)}{\partial L}\right)_{N,T}\\
&=N\left[\frac{k_{\rm B}T}{L}+\frac{kL}{2\pi^2}\left\{\frac{I_{1}\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)}{I_{0}\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\right)}-1\right\}\right]
\end{align*}
$$
ポテンシャルによって粒子の分布はが非均一化する関係で，$${\rho=N/L}$$を固定して$${L}$$を大きくすると，圧力は一定ではなく増加する。
c.
$$
\begin{align*}
\langle q\rangle&=\frac{\int_0^{\infty}{\rm d}L\exp\left(-\frac{kL^2+4\pi^2 PL}{4\pi^2k_{\rm B}T}\right)\int_0^{L}{\rm d}qq\exp\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\cos\left(\frac{2\pi q}{L}\right)\right)}{\int_0^{\infty}{\rm d}L\exp\left(-\frac{kL^2+4\pi^2 PL}{4\pi^2k_{\rm B}T}\right)\int_0^{L}{\rm d}q\exp\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\cos\left(\frac{2\pi q}{L}\right)\right)}\\
\langle L\rangle&=\frac{\int_0^{\infty}{\rm d}LL\exp\left(-\frac{kL^2+4\pi^2 PL}{4\pi^2k_{\rm B}T}\right)\int_0^{L}{\rm d}q\exp\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\cos\left(\frac{2\pi q}{L}\right)\right)}{\int_0^{\infty}{\rm d}L\exp\left(-\frac{kL^2+4\pi^2 PL}{4\pi^2k_{\rm B}T}\right)\int_0^{L}{\rm d}q\exp\left(\frac{kL^2}{4\pi^2 k_{\rm B}T}\cos\left(\frac{2\pi q}{L}\right)\right)}\\
\end{align*}
$$
Problem 5.10
式(5.12.4)の数値積分を用いることにより，式(5.9.8)の一次元周期ポテンシャルに対して式(5.9.5)の等方的NPT運動方程式を積分するプログラムを書き，図 5.3にある分布を生成できることを確認せよ。
$$
\begin{align*}
\dot{q}&=\frac{p}{m}+\frac{p_{\epsilon}}{W}q\\
\dot{p}&=-\frac{\partial U}{\partial q}-\frac{2p_{\epsilon}}{W}p-\frac{p_{\eta_1}}{Q_1}p\\
&=-\frac{m\omega^2L}{2\pi}\sin\left(\frac{2\pi q}{L}\right)-\frac{2p_{\epsilon}}{W}p-\frac{p_{\eta_1}}{Q_1}p\\
\dot{L}&=\frac{Lp_{\epsilon}}{W}\\
\dot{p}_{\epsilon}&=L\left(\mathcal{P}^{\rm(int)}-P\right)+\frac{p^2}{m}-\frac{p_{\xi_1}}{Q_1'}p_{\epsilon}\\
&=L\left(\frac{1}{L}\left[\frac{p^2}{m}-q\frac{\partial U}{\partial q}\right]-P\right)+\frac{p^2}{m}-\frac{p_{\xi_1}}{Q_1'}p_{\epsilon}\\
&=\frac{2p^2}{m}-\frac{mq\omega^2L^2}{2\pi}\sin\left(\frac{2\pi q}{L}\right)-\frac{p_{\xi_1}}{Q_1'}p_{\epsilon}-P\\
\dot{\eta}_j&=\frac{p_{\eta_j}}{Q_j}\\
\dot{\xi}_j&=\frac{p_{\xi_j}}{Q_j'}\\
\dot{p}_{\eta_M}&=G_M\\
\dot{p}_{\xi_j}&=G_j'-\frac{p_{\xi_{j+1}}}{Q_{j+1}'}p_{\xi_j}\\
\dot{p}_{\xi_M}&=G_M'\\
\end{align*}\tag{5.9.5}
$$
$$
\begin{align*}
U(q, L)&=\frac{m\omega^2L^2}{4\pi^2}\left[1-\cos\left(\frac{2\pi q}{L}\right)\right]
\end{align*}\tag{5.9.8}
$$
$$
\begin{align*}
{\rm e}^{{\rm i}L\Delta t}&={\rm e}^{{\rm i}L_{\rm NHC-baro}\frac{\Delta t}{2}}{\rm e}^{{\rm i}L_{\rm NHC-part}\frac{\Delta t}{2}}\\
&\times {\rm e}^{{\rm i}L_{\epsilon,2}\frac{\Delta t}{2}}{\rm e}^{{\rm i}L_{2}\frac{\Delta t}{2}}\\
&\times {\rm e}^{{\rm i}L_{\epsilon,1}\Delta t}{\rm e}^{{\rm i}L_{1}\Delta t}\\
&\times {\rm e}^{{\rm i}L_{2}\frac{\Delta t}{2}}{\rm e}^{{\rm i}L_{\epsilon,2}\frac{\Delta t}{2}}\\
&\times{\rm e}^{{\rm i}L_{\rm NHC-part}\frac{\Delta t}{2}}{\rm e}^{{\rm i}L_{\rm NHC-baro}\frac{\Delta t}{2}}+\mathcal{O}\left(\Delta t^3\right)\\
\end{align*}\tag{5.12.4}
$$
Pythonによる実装例を以下に示す。
# Import libraries
from fractions import Fraction
import numpy as np
import itertools
import matplotlib.pyplot as plt
import sympy
### Input parameters ###
N_c = 4 # The number of decompotiosion for the operator exp(i * L_NHC * dt / 2)
dt = 0.05 # time step
n_steps = 5000 # The number of iteration
M = 2 # The number of chains
Q = 1 # The time scale parameter for all chains
W = 5
P = 1
# Seven weights of the Yodhida scheme for a sixth-order scheme
N_sy = 7
w = [0] * N_sy
w[0] = w[6] = 0.784513610477560
w[1] = w[5] = 0.235573213359357
w[2] = w[4] = -1.17767998417887
w[3] = 1 - 2 * np.sum(w[:3])
# Initial values of dynamical variables
q = [0.5]
p = [1.0]
L = [1.0]
p_epsilon = 0.1
eta = [0.1,] * M
p_eta = [0.1,] * M
xi = [0.1,] * M
p_xi = [0.1,] * M
# Define functions
def integrateNHC(p, eta, p_eta, weight):
for n_c, n_sy in itertools.product(range(N_c), range(N_sy)):
scale = 1.0
KE = p ** 2
delta_j = w[n_sy] * dt / N_c
p_eta[M - 1] += 0.25 * delta_j * (p_eta[M - 2] ** 2 / Q - 1)
for m in range(M - 2, -1, -1):
p_eta[m] *= np.exp(-0.125 * delta_j * p_eta[m + 1])
if m == 0:
p_eta[m] += 0.25 * delta_j * (KE - 1)
else:
p_eta[m] += 0.25 * delta_j * (p_eta[m - 1] ** 2 / Q - 1)
p_eta[m] *= np.exp(-0.125 * delta_j * p_eta[m + 1])
scale *= np.exp(-0.5 * delta_j * p_eta[0] / Q)
KE *= np.exp(-1.0 * delta_j * p_eta[0] / Q)
for m in range(M):
eta[m] += 0.5 * delta_j * p_eta[m] / Q
for m in range(0, M - 1, 1):
p_eta[m] *= np.exp(-0.125 * delta_j * p_eta[m + 1])
if m == 0:
p_eta[m] += 0.25 * delta_j * (KE - 1)
else:
p_eta[m] += 0.25 * delta_j * (p_eta[m - 1] ** 2 / Q - 1)
p_eta[m] *= np.exp(-0.125 * delta_j * p_eta[m + 1])
p_eta[M - 1] += 0.25 * delta_j * (p_eta[M - 2] ** 2 / Q - 1)
p *= scale
return p, eta, p_eta
def calculateForce(q, L, variable = "q"):
theta = 2 * np.pi * q / L
if variable == "q":
return -0.5 * L * np.sin(theta) / np.pi
else:
return 0.5 * q * np.sin(theta) / np.pi - 0.5 * L * (1 - np.cos(theta)) / (np.pi ** 2)
def integrateEpsilon1(q, epsilon, L, L_0, p_epsilon):
epsilon += p_epsilon * dt / W
q *= np.exp(epsilon)
L_new = L_0 * np.exp(epsilon)
return q, epsilon, L_new
def integrateEpsilon2(p, q, L, p_epsilon):
F_q = calculateForce(q, L, "q")
F_L = calculateForce(q, L, "L")
G_epsilon = 2 * p ** 2 + q * F_q + L * F_L - P * L
return p_epsilon + 0.5 * G_epsilon * dt
def integrate1(q, L, p, p_epsilon):
index = p_epsilon * dt / W
q *= np.exp(index)
# Use sympy to evaluate sinh(x) / x analytically
sympy.var('x')
sinh_per_x = sympy.limit(sympy.sinh(x) / x, x, 0.5 * index).evalf()
q += dt * p * np.exp(0.5 * index) * float(sinh_per_x)
if q > L:
while q > L:
q -= L
elif q < 0:
while q < 0:
q += L
return q
def integrate2(q, L, p, p_epsilon):
index = 2 * 0.25 * p_epsilon * dt / W
p *= np.exp(-2 * index)
F = calculateForce(q, L)
# Use sympy to evaluate sinh(x) / x analytically
sympy.var('x')
sinh_per_x = sympy.limit(sympy.sinh(x) / x, x, index).evalf()
p += 0.5 * dt * F * np.exp(-1 * index) * float(sinh_per_x)
return p
###---- The below is the code ----###
# The numerical integration
q_now = q[-1]
p_now = p[-1]
L_now = L[-1]
epsilon_now = 0
p_epsilon_now = p_epsilon
# NHC integration for equilibration
for _ in range(1000):
p_now, eta, p_eta = integrateNHC(p_now, eta, p_eta, 1)
p_now = integrate2(q_now, L_now, p_now, p_epsilon_now)
q_now = integrate1(q_now, L_now, p_now, p_epsilon_now)
p_now = integrate2(q_now, L_now, p_now, p_epsilon_now)
p_now, eta, p_eta = integrateNHC(p_now, eta, p_eta, 1)
# MTK integration
for n_step in range(n_steps):
if n_step % 25 == 0:
print(f"{n_step},{q_now},{p_now},{L_now}")
p_epsilon_now, xi, p_xi = integrateNHC(p_epsilon_now, xi, p_xi, W)
p_now, eta, p_eta = integrateNHC(p_now, eta, p_eta, 1)
p_epsilon_now = integrateEpsilon2(p_now, q_now, L_now, p_epsilon_now)
p_now = integrate2(q_now, L_now, p_now, p_epsilon_now)
q_now = integrate1(q_now, L_now, p_now, p_epsilon_now)
q_now, epsilon_now, L_now = integrateEpsilon1(q_now, epsilon_now, L_now, L[0], p_epsilon_now)
p_now = integrate2(q_now, L_now, p_now, p_epsilon_now)
p_epsilon_now = integrateEpsilon2(p_now, q_now, L_now, p_epsilon_now)
p_now, eta, p_eta = integrateNHC(p_now, eta, p_eta, 1)
p_epsilon_now, xi, p_xi = integrateNHC(p_epsilon_now, xi, p_xi, W)
q.append(q_now)
p.append(p_now)
L.append(L_now)
# Export time history of physical dynamical variables to csv file
result = np.array([range(n_steps + 1), q, p, L]).T
header = 'n_step, q, p, L'
with open('time_history.csv', 'wb') as f:
np.savetxt(f, result, delimiter = ",", header = header, fmt = "%.5e")
# Plot frequencies of q
plt.hist(q, bins = 25, density = True, label = 'numerical')
plt.legend()
plt.xlabel('q')
plt.ylabel('density')
plt.savefig('q2.png')
plt.figure()
# Plot frequencies of q
plt.hist(L, bins = 25, density = True, label = 'numerical')
plt.legend()
plt.xlabel('L')
plt.ylabel('density')
plt.savefig('L2.png')
plt.figure()
Problem 5.11
4.13節における動径分布関数を計算するためのアルゴリズムをNPTアンサンブル用にどのように修正すべきであろうか？
NPTアンサンブル下における体積$${V}$$の発生確率を$${p(V|N,P,T)}$$とおくと，NPTアンサンブル用に修正された式(4.13.26)は
$$
\begin{align*}
g_{ab}(r_i)&=\int_{0}^{\infty}{\rm d}Vp(V|N,P,T)\frac{h_{ab}(i)}{4\pi\rho_br_i^2\Delta rN_{\rm conf}N_a}\\
&=\frac{h_{ab}(i)}{4\pi N_br_i^2\Delta rN_{\rm conf}N_a}\int_{0}^{\infty}{\rm d}Vp(V|N,P,T)V\\
&=\frac{h_{ab}(i)\langle V\rangle}{4\pi N_br_i^2\Delta rN_{\rm conf}N_a}\\
&\simeq\frac{h_{ab}(i)}{4\pi\langle\rho_b\rangle r_i^2\Delta rN_{\rm conf}N_a}
\end{align*}
$$
となる。
つまり，$${\rho_b}$$をMD計算で得られた平均値$${\langle\rho_b\rangle}$$に置き換えればよい。
Problem 5.12
5.13節のROLLアルゴリズムをシミュレーションボックスの揺らぎが非等方的な場合に拡張せよ。
流石に式展開を追うのはしんど過ぎたので省略。。。
興味ある方は↓を参照のこと。
https://doi.org/10.1016/j.chemphys.2010.02.014
Problem 5.13
a. $${\alpha > \beta}$$に対して$${h_{\alpha\beta}=0}$$となる拘束条件に対してラグランジュの未定乗数法を用いることにより，単に上三角のセル行列を扱う問題となり，その結果として式(5.10.2)が持つセル全体の回転運動を消去できることを示せ。
b. $${\mathbf{p}_g-\mathbf{p}_g^{\rm T}=\mathbf{0}}$$となる拘束条件に対してラグランジュの未定乗数法を用いることにより，圧力テンソル$${P_{\alpha\beta}^{\rm (int)}}$$が陽的に対称化され，その結果として式(5.10.2)が持つセル全体の回転運動を消去できることを示せ。なぜこのスキームの方が5.13節のROLLアルゴリズムの中に含めることが容易なのであろうか？
a. $${\alpha > \beta}$$に対して$${h_{\alpha\beta}=0}$$のとき，ベクトル$${\mathbf{a}}$$は常に$${x}$$軸に平行であるため，セルの回転運動は生じない。
b. $${\mathbf{p}_g}$$が対称行列の場合，セルの角運動量が保存するため，初期条件で角運動量を$${\mathbf{0}}$$に選んでおけばセルの回転運動は生じない。
$${\mathbf{p}_g-\mathbf{p}_g^{\rm T}=\mathbf{0}}$$の条件は圧力テンソルを対称化することによって自動的に満たされる。対称化は圧力テンソルの計算時に追加で$${\tilde{P}_{\alpha\beta}^{\rm (int)}=(P_{\alpha\beta}^{\rm (int)}+P_{\beta\alpha}^{\rm (int)})/2}$$を計算すればよいだけなので実装は簡単である。
参考文献
Mark Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford Graduate Texts)