NPT条件下で時間依存の物理量を解析する方法

本記事は2022/8/17に修正されました。
  １．数式を画像からtexに差し替え
  ２．NPT/NVTアンサンブル間の期待値の関係における式展開の誤りを修正

NPTアンサンブル（粒子数，圧力，温度一定）やNVTアンサンブル（粒子数，体積，温度一定）等の平衡状態でのMDシミュレーションでは，配置変化の時刻歴（MD trajectory）が形式上得られるものの，それが時刻歴として物理的意味を成しているかは定かではありません（平衡状態のMDシミュレーションは，あくまでボルツマン分布等の確率分布のサンプリングです）。例えば，Langevin thermostatを用いた場合は，タンパクやポリマー等の高分子のダイナミクスがMDシミュレーション上で著しく遅く場合があることが知られています。
NVEアンサンブル（粒子数，体積，エネルギー一定）はハミルトンダイナミクスとなり，MD trajectoryはサンプリングという意味のみならず，時刻歴としても物理的意味を持ちます。そのため，MD trajectoryを時刻歴データとして取り扱いたい場合は，NVEアンサンブルを適用する必要があります。
ですが，300K/1atmといったように特定の温度/圧力条件下での物理量を求めたいことが多いと思います。その物理量が拡散係数等の時刻歴の解析を必要とするような場合，NVEアンサンブルでのMDシミュレーションをNPTアンサンブルと紐づける必要があります。

異なるアンサンブル間での期待値の関係

まずはNPTアンサンブルとNVTアンサンブルにおける，物理量$${A}$$の期待値がどのような関係式で結ばれるのかを考えます。

$$
\begin{align*}
\langle A\rangle_{NPT}&=\frac{\int{\rm d}V\int{\rm d}\mathbf{x}A(\mathbf{x}){\rm e}^{-\beta(\mathcal{H}(\mathbf{x})+PV)}}{\int{\rm d}V\int{\rm d}\mathbf{x}{\rm e}^{-\beta(\mathcal{H}(\mathbf{x})+PV)}}\\
&=\int{\rm d}V\left\{\frac{{\rm e}^{-\beta PV}\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}{\int{\rm d}V{\rm e}^{-\beta PV}\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\right\}\left\{\frac{\int{\rm d}\mathbf{x}A(\mathbf{x}){\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\right\}\\
&=\int{\rm d}V\left\{\frac{{\rm e}^{-\beta PV}\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}{\int{\rm d}V{\rm e}^{-\beta PV}\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\right\}\langle A\rangle_{NVT}\\
&=\int{\rm d}Vp(V|N,P,T)\langle A\rangle_{NVT}
\end{align*}
$$

$${p(V|N,P,T)}$$はNPTアンサンブルにおける体積$${V}$$の発生確率です。つまり，NPTアンサンブルでの期待値は，NVTアンサンブルでの期待値を$${p(V|N,P,T)}$$で重みづけし，足し合わせることで得られることになります。
同様にして，NVTアンサンブルとNVEアンサンブルにおける物理量$${A}$$の期待値は以下のような関係式で結ばれます。

$$
\begin{align*}
\langle A\rangle_{NVT}&=\frac{\int{\rm d}\mathbf{x}A(\mathbf{x}){\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\\
&=\frac{\int{\rm d}E\int{\rm d}\mathbf{x}\delta(\mathcal{H}(\mathbf{x})-E)A(\mathbf{x}){\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\\
&=\frac{\int{\rm d}E{\rm e}^{-\beta E}\int{\rm d}\mathbf{x}\delta(\mathcal{H}(\mathbf{x})-E)A(\mathbf{x})}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\\
&=\frac{\int{\rm d}E{\rm e}^{-\beta E}\left\{\langle A\rangle_{NVE}\int{\rm d}\mathbf{x}\delta(\mathcal{H}(\mathbf{x})-E)\right\}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\\
&=\frac{\int{\rm d}E\left\{{\rm e}^{-\beta E}\int{\rm d}\mathbf{x}\delta(\mathcal{H}(\mathbf{x})-E)\right\}\langle A\rangle_{NVE}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\\
&=\frac{\int{\rm d}E\left\{\int{\rm d}\mathbf{x}\delta(\mathcal{H}(\mathbf{x})-E){\rm e}^{-\beta\mathcal{H}(\mathbf{x})}\right\}\langle A\rangle_{NVE}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\\
&=\int{\rm d}E\left\{\frac{\int{\rm d}\mathbf{x}\delta(\mathcal{H}(\mathbf{x})-E){\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}{\int{\rm d}\mathbf{x}{\rm e}^{-\beta\mathcal{H}(\mathbf{x})}}\right\}\langle A\rangle_{NVE}\\
&=\int{\rm d}Ep(E|N,V,T)\langle A\rangle_{NVE}\\
\end{align*}
$$

以上より，NPTアンサンブルにおける期待値はNVEアンサンブルを用いて

$$
\begin{align*}
\langle A\rangle_{NPT}&=\int{\rm d}Vp(V|N,P,T)\int{\rm d}Ep(E|N,V,T)\langle A\rangle_{NVE}\\
\end{align*}
$$

と表すことができます。

計算手順

上式をまともに評価する場合，以下のような手順が必要になります。

1. NPTアンサンブルのサンプリングを実施し，$${p(V|N,P,T)}$$を評価する

2. $${p(V|N,P,T)}$$の主要な値を占める$${V}$$を幾つかピックアップし，それぞれの$${V}$$に対してNVTアンサンブルのサンプリングを実施し，$${p(E|N,V,T)}$$を評価する

3. $${p(E|N,V,T)}$$の主要な値を占める$${E}$$を幾つかピックアップし，幾つかの$${(V,E)}$$の組み合わせに対してNVEアンサンブルのサンプリングをそれぞれ実施する

4. 各NVEアンサンブルにて$${\langle A\rangle_{NVE}}$$を評価し，$${p(V|N,P,T)}$$と$${p(E|N,V,T)}$$を用いた重みづけ平均を行う

これだと，とても手間がかかってしまいますね。。。
$${V}$$と$${E}$$がどちらも期待値からあまり揺らがない場合，

$$
\begin{align*}
\langle A\rangle_{NPT}&=\int{\rm d}Vp(V|N,P,T)\int{\rm d}Ep(E|N,V,T)\langle A\rangle_{NVE}\\
&\simeq \int{\rm d}V\delta(V-\langle V\rangle_{NPT})\int{\rm d}Ep(E|N,V,T)\langle A\rangle_{NVE}\\
&=\int{\rm d}Ep(E|N,\langle V\rangle_{NPT},T)\langle A\rangle_{N\langle V\rangle_{NPT}E}\\
&\simeq\int{\rm d}E\delta(E-\langle E\rangle_{N\langle V\rangle_{NPT}T})\langle A\rangle_{N\langle V\rangle_{NPT}E}\\
&=\langle A\rangle_{N\langle V\rangle_{NPT}\langle E\rangle_{N\langle V\rangle_{NPT}T}}\\
&\simeq\langle A\rangle_{N\langle V\rangle_{NPT}\langle E\rangle_{NPT}}
\end{align*}
$$

と近似することができます。この式に従うと手順を簡略化することができます。

1. NPTアンサンブルのサンプリングを実施し，$${\langle V\rangle_{NPT}，\langle E\rangle_{NPT}}$$を評価する

2. $${V=\langle V\rangle_{NPT}，E=\langle E\rangle_{NPT}}$$でのNVEアンサンブルのサンプリングを実施し，$${\langle A\rangle_{NVE}}$$を計算する

原理的には，NPTアンサンブルとNVEアンサンブルのサンプリングを1回ずつ行えば目的の期待値が得られることになります。