12 データ解析

熱平衡状態の分子動力学(MD)を考える。 力学量（温度、内部エネルギー、圧力など）がそれぞれの平均値の回りで 揺らぎ始めたとき、熱平衡に達したと考える。

（１）熱平衡状態

通常、$\Delta t$ をなるべく大きくとるように設定するため一部の粒子の 運動が激しい初期段階では数値積分の精度が悪い。

どのような熱平衡状態実現させるかによってMD法の計算方法が異なっている。

（２）NVE一定のMD法

粒子数$N$ を一定にして、粒子間のポテンシャルを設定し、初期条件を与えて 運動方程式を数値積分するとこの力学系は保存系である。

力学エネルギー$E$ (運動エネルギー＋ポテンシャルエネルギー)は一定である。 （ただし、差分化による誤差、まるめ誤差がなければ）

その他に
NVT一定（粒子数$N$ 、 体積$V$ 、温度$T$ を一定にしたシミュレーション）
NPT一定（粒子数$N$ 、 体積$P$ 、温度$T$ を一定にしたシミュレーション）
といったMD法がよく使われる。これらは実際のシミュレーションで使われる 発展した方法である。発展的な方法なので、説明はこの程度にどどめる）

（３）時間平均

力学量の平均（これが実験の測定値と比較する量）

$$\langle Q \rangle = \lim_{t\rightarrow\infty} \frac{1}{t} \int_{0}^{t} Q(s) ds$$ (163)

時刻$t_{1}$ と$t_{2}$ との間の平均

$$\langle Q \rangle = \frac{1}{t_{1}-t_{2}} \int_{t_{1}}^{t_{2}} Q(s) ds$$ (164)

$$= \frac{1}{N_{s}N} \sum_{n=0}^{N_{s}-1} \sum_{i=1}^{N} Q_{i}(t_{1}+n \Delta t)$$ (165)

$N$ : 粒子数、 $N_{s}$ :サンプル数、 $\Delta t$ : サンプルを採取するインターバル時間

サンプル数を増すほど統計精度は良くなる。 インターバル$\Delta t$ を大きくとるほど、 位相空間( $\vec{r}_{i}, \vec{p}_{i}$ または $\vec{r}_{i}, \vec{v}_{i}$ )の 広い範囲を考慮したことになる。

力学量$Q$ の分散$\sigma(Q)$

$$\sigma(Q) = \langle (Q- \langle Q \rangle)^{2} \rangle = \ \langle Q^{2} \rangle - \langle Q \rangle^{2}$$ (166)

$$\langle Q^{2} \rangle = \frac{1}{N_{s}N} \sum_{n=0}^{N_{s}-1} \sum_{i}^{N} Q_{i}^{2} (t_{i}+ n\Delta t)$$ (167)

分散を計算しておくと統計誤差を知ることができる。 時間平均と粒子平均を区別なく記号 $\langle \cdot\cdot\cdot \rangle$ を 用いる場合がある。

(4)温度

系の粒子運動の激しさを表す。

$$\langle １粒子の運動エネルギー \rangle = \frac{m}{2} \langle (\vec{v})^{2} \rangle = 3 \times \frac{1}{2} k_{\rm B} T$$ (168)

$$T = \frac{m}{3 k_{\rm B}} \langle (\vec{v})^{2} \rangle = \frac{2}{3 N k_{\rm B}} \langle K \rangle$$ (169)

$$K \equiv \frac{1}{2} \sum_{i} m_{i} (\vec{v}_{i})^{2} K:全運動エネルギー$$ (170)

$\langle$ 1粒子当たりのポテンシャルエネルギー $\rangle$

$$= \frac{1}{N} \langle U \rangle$$ (171)

$$= \frac{1}{N} \sum_{i=1}^{N}\sum_{j>i}^{N} \langle \phi(r_{ij}) \rangle$$ (172)

MD法では、温度や系の安定性を知るためにもこのようなデータを計算しておく 必要がある。

$$力学エネルギー E = K + U$$ (173)

$$平均\langle E \rangle = \langle K \rangle + \langle U \rangle$$ (174)

その他に、力学エネルギーの保存具合を知るために、運動エネルギーおよび ポテンシャルエネルギーを必ず計算する必要がある。 $\longrightarrow$ モニターリング

平均 $\langle E \rangle$ は、熱力学では内部エネルギーと呼んでいる量に対応する。

（５）圧力

運動方程式

$$m_{i} \frac{d^{2} \vec{r}_{i}}{d t^{2}} = \vec{G}_{i} = \vec{F}_{i} + \vec{W}_{i}$$ (175)

$\vec{r}_{i}$ をかけて時間についての平均を計算する

$$\langle m_{i} \vec{r}_{i} \cdot \frac{d^{2} \vec{r}_{i}}{d t^{2}} \rangle = \langle \vec{r}_{i} \cdot \vec{G}_{i} \rangle$$ (176)

$$(左辺) = \frac{1}{t_{12}} \int_{t_{1}}^{t_{2}} m_{i} \vec{r}_{i} \cdot \frac{d^{2} \vec{r}_{i}}{d t^{2}} ds$$ (177)

$$= \frac{1}{t_{12}} \left\{ \left[ m_{i} \vec{r}_{i} \cdot \frac{d \... ...2}} m_{i} \frac{d \vec{r}_{i}}{d t} \cdot \frac{d \vec{r}_{i}}{d t} ds \right\}$$ (178)

$$= \frac{1}{t_{12}} \left\{ \left[ \left. m_{i} \vec{r}_{i} \cdot \f... ... \vec{r}_{i} \cdot \frac{d \vec{r}_{i}}{d t} \right\vert _{t_1} \right] \right.$$ (179)

$$\ ... ...2}} m_{i} \frac{d \vec{r}_{i}}{d t} \cdot \frac{d \vec{r}_{i}}{d t} ds \right\}$$ (180)

$$= \ ... ...2}} m_{i} \frac{d \vec{r}_{i}}{d t} \cdot \frac{d \vec{r}_{i}}{d t} ds \right\}$$ (181)

$$= -2 \left\langle \frac{1}{2} m_{i} v_{i}^{2} \right\rangle$$ (182)

ビリアルの定理

$$-2 \left\langle \frac{1}{2} m_{i} v_{i}^{2} \right\rangle = \langle \vec{r}_{i} \cdot \vec{G}_{i} \rangle$$ (183)

$$\left\langle \sum_{i} \vec{r}_{i} \cdot \vec{G}_{i} \right\rangle = \left\langle \sum_{i} \vec{r}_{i} \cdot \vec{F}_{i} \right\rangle + \left\langle \sum_{i} \vec{r}_{i} \cdot \vec{W}_{i} \right\rangle$$ (184)

$$\left\langle \sum_{i} \vec{r}_{i} \cdot \vec{W}_{i} \right\rangle = \left\langle \int \sum_{i} \delta (\vec{r}-\vec{r}_{i}) \vec{r}_{i} \cdot \left( - \vec{n}(\vec{r}_{i}) p d S \right) \right\rangle$$ (185)

$$= - \int_{\rm surface} \vec{r} \cdot \vec{n}(\vec{r}) p dS$$ (186)

$$= - p \int_{\rm surface} \vec{r} \cdot \vec{n} dS$$ (187)

$$= - p \int_{\rm volum} {\rm div} \vec{r} dV$$ (188)

$$= - 3 p \int_{\rm volum} dV$$ (189)

$$= - 3 p \Omega$$ (190)

したがって、

$$-2 \left\langle \frac{1}{2} m_{i} v_{i}^{2} \right\rangle = \left\langle \sum_{i} \vec{r}_{i} \cdot \vec{F}_{i} \right\rangle - 3 p \Omega$$ (191)

$$p = \frac{1}{3 \Omega} \left\{ 2 \left\langle K \right\rangle + \left\langle \sum_{i} \vec{r}_{i} \cdot \vec{F}_{i} \right\rangle \right\}$$ (192)

$$= \frac{2}{3\Omega} \left\langle K \right\rangle + \frac{1}{3\Omega} \left\langle \sum_{i} \vec{r}_{i} \cdot \vec{F}_{i} \right\rangle$$ (193)

$$= \frac{Nk_{\rm B} T}{\Omega} + \frac{1}{3\Omega} \left\langle \sum_{i} \vec{r}_{i} \cdot \vec{F}_{i} \right\rangle$$ (194)

内力 $\vec{F}_{i}$ を考えないと、 $p\Omega=Nk_{\rm B} T$ (理想気体の 状態方程式)となる。

２体間ポテンシャルの場合

$$\vec{F}_{i} = \sum_{j\ne i} \vec{f}_{ij} = \ \sum_{j\ne i} (-1) \left. \frac{d \phi}{d r} \right\vert _{r=r_{ij}} \frac{\vec{r}_{ij}}{r_{ij}}$$ (195)

$$\left\langle \sum_{i} \vec{r}_{i} \cdot \vec{F}_{i} \right\rangle = ........... = - \frac{1}{2} \left\langle \sum_{ij}^{i\ne j} r_{ij} \left. \frac{d \phi}{d r} \right\vert _{r=r_{ij}} \right\rangle$$ (196)

$$p = \frac{2}{3\Omega} \left\langle K \right\rangle - \frac{1}{6\Omega... ...{i\ne j} r_{ij} \left. \frac{d \phi}{d r} \right\vert _{r=r_{ij}} \right\rangle$$ (197)

問：式(196)を導出せよ。

状態方程式(まとめ)

理想気体の状態方程式

$$p\Omega = Nk_{\rm B} T$$ (198)

ビリアル展開式(カマリン-オネスの状態方程式)

$$p\Omega = Nk_{\rm B} T \left(1+ \frac{A_{2}}{\Omega} + \frac{A_{3}}{\Omega^{2}} + .... \right)$$ (199)

         $A_{2}, A_{3}$ ： 第２,第３ビリアル係数。温度のみの関数。 これらは統計力学の手法を用いて導出することができる。

ファンデルワールス(van del Waals)の実在気体の方程式

$$\left(p+\frac{a}{\Omega^{2}}\right) \left(\Omega-b \right) = Nk_{\rm B} T$$ (200)

（６）局所構造

各粒子の配置は時々刻々変化して行くが、その変化の仕方は無規則なものでは ない。温度、圧力といったマクロな(巨視的な)物理量に依存して、 ミクロな(微視的な)各粒子のお互いの配置に一定の性質が存在する。 そのような性質を局所構造と呼ぶ。ここでは代表的な物理量だけを学習する。

動径分布関数
1つの粒子のまわりの距離$r$ の位置に存在するもう一つの粒子の平均的な 分布を表す関数

$$g(r) = \frac{\Omega}{N} \frac{n(r)}{4\pi r^{2}}$$ (201)

ここで、 $\frac{\Omega}{N}$ は1粒子当たりの体積、 $r$ と$r+dr$ で 作られる球殻に存在する粒子数が$n(r)dr$ となるように、$n(r)$ が 決められる。

$r\rightarrow\infty$ のとき

$$n(r)dr \longrightarrow \frac{N}{\Omega} 4\pi r^{2} dr$$ (202)

$$g(r) \longrightarrow \frac{\Omega}{N} \frac{1}{4\pi r^{2}} \frac{N}{\Omega} 4\pi r^{2} dr = 1$$ (203)

$g(r)$ の計算法 $r$ の刻幅を$\Delta r$ と決めて、粒子間距離$r_{ij}$ のサンプルデータ に対して

$$m = \left[ \frac{r_{ij}}{\Delta r} \right] + 1$$ (204)

$$N(m) = N(m) + 1$$ (205)

を計算する。その後、十分なサンプル数のもとで、

$$g(m) = \frac{\Omega}{N} \frac{1}{4\pi r^{2}_{m}} \frac{\langle N(m) \rangle }{N} \frac{1}{\Delta r}$$ (206)

から$g(r)$ を得る。

結合角分布関数

構造因子