＊）分子動力学/テキスト目次

最終更新日時: 2025年08月25日 12:57

TODOs

• nothing

目次

• 第1章 分子動力学とは何か
  • 1.1 分子動力学(MD)の目的と意義
  • 1.2 MDシミュレーションの基本的な流れ
  • 1.3 古典近似の前提
• 第2章 理論的背景の再確認
  • 2.1 古典力学の再確認
  • 2.2 統計力学・熱力学の再確認
  • 2.3 分子動力学における熱力学の位置づけ
• 第3章 原子間ポテンシャルと力場
  • 3.1 原子間相互作用の基本
  • 3.2 代表的な力場(Force Field)
  • 3.3 応用への注意点
• 第4章 分子動力学シミュレーションの実装
  • 4.1 数値積分法の実際（Verlet法の使い方）
  • 4.2 境界条件と系のセッティング
  • 4.3 長距離相互作用の扱い
  • 4.4 拘束手法
  • 4.5 温度・圧力制御
  • 4.6 小規模系での練習例
• 第5章 シミュレーションデータの解析・可視化
  • 5.1 構造の評価
  • 5.2 物性・動力学的評価
  • 5.3 可視化ツール
• 第6章 応用と高度な話題
  • 6.1 フリーエネルギー計算
  • 6.2 高度なサンプリング手法
  • 6.3 大規模計算・並列化
  • 6.4 分野別応用例

第1章 分子動力学とは何か

本章では分子動力学が解くことを目指す代表的な項目、それらをどういう流れ・どういう枠組みでそれらを解くのか。 また、古典力学、量子力学と比してどういう前提やどういう考え方でアプローチするのかを説明する

1.1 分子動力学(MD)の目的と意義

分子動力学(MD)は、多数の原子や分子が互いに作用し合いながら時間発展する過程を、古典力学に基づいて数値的に解く手法である。粒子同士のポテンシャルエネルギーを定義し、運動方程式を繰り返し解くことで、微視的な構造変化や物性変化を追跡できる。

質量$m_i$をもつ粒子$i$の位置ベクトルを$\mathbf{r}_i(t)$、相互作用ポテンシャルを$V(\mathbf{r}_1,\dots,\mathbf{r}_N)$とおくと、ニュートンの運動方程式は下記の通り。

$$m_i \frac{d^2 \mathbf{r}_i}{dt^2} = - \nabla_{\mathbf{r}_i} V(\mathbf{r}_1,\dots,\mathbf{r}_N)$$

これを離散時間ステップごとに数値解を求めることで、 粒子の位置と速度を逐次更新する。系の時間発展を追い続けると、 構造ゆらぎや運動エネルギーの時系列データが得られ、そこから拡散や相転移などの動的現象を評価することが目的である。

1.1.1 応用範囲と実例

1. 材料物性の評価
   例えば合金や高分子材料における拡散係数や弾性係数を推定する。局所的な欠陥生成や粒界挙動のシミュレーションも可能。

2. バイオ分野(タンパク質・核酸)
   タンパク質の立体構造変化や酵素反応における基質の結合モードを解析する。リガンドドッキングや分子認識のダイナミクスを評価できる。

3. 化学反応の経路探索
   触媒表面上での吸着・脱離過程や、液相中での分子衝突頻度を追うことで、反応速度論の微視的起源を検討する。

4. ナノ材料の設計
   グラフェンやカーボンナノチューブなどでの機械的強度や層間相互作用を評価し、新規材料設計の指針を得る。

5. 液体・溶液の熱物性
   粘度、熱伝導率、相変化(凝固・融解)などのマクロ物性を、粒子同士の相互作用と速度分布の観点から推定する。

1.1.2 古典近似

MDは量子力学的な電子状態を明示的に扱わず、原子全体を質点として扱うのが古典近似の基本となる。 電子構造は明示的に考慮せず、相互作用ポテンシャルでまとめて表す。

室温付近かつ重原子が多い系では、量子効果による誤差は比較的小さい。

量子力学的効果を無視する背景には、熱的な波動長の評価がある。室温付近で重い原子の場合、粒子の熱的ドブロイ波長

$$\Lambda \sim \sqrt{\frac{h^2}{2\pi m k_B T}}$$

が原子間距離より十分小さくなるため、量子重ね合わせや不確定性の影響が小さい。例えばアルゴン原子(質量約$40\,\mathrm{u}$)を室温($\sim 300\,\mathrm{K}$)で扱うと、$\Lambda$はオーダー$10^{-11}\,\mathrm{m}$以下になり、原子間距離($\sim 3\times 10^{-10}\,\mathrm{m}$)より格段に短い。

高周波振動(例: 水素結合の伸縮)などは必要に応じて拘束や修正ポテンシャルを導入することで近似的に扱える。

高周波振動を含む系では、原子間の振動数が$10^{13}\,\mathrm{Hz}$以上(振動エネルギーが数百$\mathrm{meV}$相当)になることが多い。特に水素結合やO-H伸縮振動は$3000\,\mathrm{cm}^{-1}$前後(およそ$9\times 10^{13}\,\mathrm{Hz}$)の領域にあり、時間刻みを極端に小さくしないと正確な再現が難しくなる。そこで、拘束(例: SHAKE/RATTLE)や修正ポテンシャルを導入して振動自由度をある程度固定し、より長い時間刻みでシミュレーションを進める方法が選ばれる。こうすると高周波成分を省略または近似しつつ、低周波の主要挙動(例: 分子間緩和、拡散運動)を効率よく追える。

1.1.3 エルゴード仮説とマクロ物性予測

エルゴード仮説とは、系の時間平均とアンサンブル平均が一致するという仮説である。 これにより、長時間のシミュレーションから得られる物理量の時間平均を 統計力学的なアンサンブル平均とみなし、平均的な物性を計算できるようになる。

統計力学で言うアンサンブル平均とは、 系のすべての可能な状態(位相空間上の座標と運動量の組)を、 一定の確率分布に従って重み付けして求める平均のことである。

カノニカルアンサンブル(NVT)の確率密度$\rho(\Gamma)$は、 $\Gamma$は全粒子の座標・運動量をまとめた位相空間点、 $H(\Gamma)$はハミルトニアン、 $\beta = 1/(k_B T)$とすると、下記のように表現できる。

$$\rho(\Gamma) \propto \exp\bigl(-\beta H(\Gamma)\bigr)$$

アンサンブル平均とは系がとりうるあらゆる状態を統計力学的な重み付けで平均化することである。 ある物理量$A(\Gamma)$の期待値を$\langle A \rangle$とすると下記のように 位相空間全体にわたる積分(あるいは総和)によって表すことができる。

$$\langle A \rangle = \int A(\Gamma)\,\rho(\Gamma)\, d\Gamma$$

時間平均$\overline{A}$は下記のように表現できる。

$$\overline{A} \;=\; \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} A\bigl(\mathbf{r}(t), \mathbf{p}(t)\bigr) \, dt$$

時間平均$\overline{A}$が、アンサンブル平均$\langle A \rangle \;$に一致するというのがエルゴード仮説であり、 系がこれを満たす時、「系がエルゴード的」であるという。

1.1.4 時間相関を考慮する、より進んだマクロ物性予測

シミュレーションを十分に長時間実行すれば、エルゴード仮説によって時間平均はアンサンブル平均に対応づけられる。しかし、粘度や拡散係数のような輸送係数は、流体の流れや分子の拡散における「時間相関」が鍵となるため、単純にある物理量の時間平均だけでは不十分になる。

時間相関とは、 ある物理量が「時刻0」と「それ以降の時刻t」でどの程度 相関があるかという概念である。

時刻が経過すると、その相関がどのように消えていくか(減衰速度)を見ることで、 粒子が元の運動の状態をどのくらい保持し続けるかを知ることができる。

拡散係数の場合、はじめに分子が特定方向に高速で運動していたとき、その速度ベクトルをどれだけ保持し続けるか(速度自己相関)が拡散過程を支配する。完全ランダムに散らばるなら相関はすぐにゼロに近づき、ある程度の慣性や相互作用で同じ方向を維持すれば相関がゆっくり減衰する。

粘度は「せん断応力が空間的・時間的にどのように伝わり、緩和していくか」 を反映する量なので、 運動量フラックス(応力テンソル)がどれくらい過去と似た振る舞いを保ち続けるかを調べる必要がある。

要するに、時間相関は「現在の状態が過去の状態をどの程度引きずっているか」 を示す指標で、その緩和パターンから粒子間相互作用の影響度合いを知ることができる。

グリーン-クボ(Green-Kubo)の公式は、運動量フラックスや速度自己相関関数の時間積分から輸送係数を求める枠組みを与える。粘度$\eta$の一例としては、

$$\eta = \frac{V}{k_B T} \int_0^\infty \langle P_{xy}(0)\,P_{xy}(t)\rangle \, dt$$

ここで

• $P_{xy}$は応力テンソルの$xy$成分
• $\langle \cdot \rangle$は統計力学的平均を示す
• $V$は系の体積
• $k_B$はボルツマン定数、$T$は温度

また拡散係数$D$は速度自己相関関数を積分する形で、

$$D = \int_0^\infty \langle \mathbf{v}(0)\cdot\mathbf{v}(t)\rangle \,dt$$

表現できる(定数係数は省略)。

これらはいずれも「時間相関関数を時間積分する」手法であり、 単に物理量そのものを時間平均するだけでは得られない。 輸送係数の本質は流体や分子の“動き方”にあり、 その動きの履歴(相関)をしっかり評価するための仕組みが必要となる。

こうした理論的背景に基づき、粘度や拡散係数を定量的に取り扱えるのが分子動力学シミュレーションの強みでもある。

1.1.5 ドブロイの関係式

ドブロイ波長を求め、粒子が古典的にふるまうのか、波動性が顕著になるのかを理解するためのドブロイ波長を求めるための説明の端緒となるドブロイの関係式の理解から始める。

ドブロイの関係式 (波長 $\lambda$ と運動量 $p$ の関係 $\lambda = \frac{h}{p}$) は、もともと光子に対するエネルギーと運動量の式から着想を得て、それを物質粒子にも拡張したものとして導かれる。以下に初学者向けに簡単な流れを示す。

1.1.5.1. 光子のエネルギーと運動量

電磁波(光)は粒子(光子)としての性質を持ち、プランク=アインシュタインの関係式がある

$$E = h \nu$$

ここで $h$ はプランク定数、$\nu$ は電磁波の振動数ある。

光子のエネルギーは運動量 $p$ と光速 $c$ を用いて

$$E = pc$$

とも書ける。すると

$$pc = h \nu$$

となる。

1.1.5.2. 光の振動数と波長の関係

波としての光に対しては、位相速度が $c$ であるため

$$\nu \,\lambda = c$$

が成り立つ。よって

$$\nu = \frac{c}{\lambda}$$

を用いると、(1) の $pc = h \nu$ に代入して

$$p c = h \frac{c}{\lambda} \quad \Rightarrow\quad p = \frac{h}{\lambda}$$

となる。これが光(質量0粒子)に対する運動量と波長を関連付ける式となる。

1.1.5.3. ドブロイの仮説

ドブロイは、光子だけでなく質量をもつ物質粒子についても、同様に波動性を持つと提案した。すなわち「粒子の運動量 $p$ と対応する波(物質波)の波長 $\lambda$ は

$$\lambda = \frac{h}{p}$$

で結びつく」というのがドブロイの仮説である。歴史的には電子回折(デヴィッソン=ガーマーの実験など)で実証され、後にあらゆる粒子に普遍的に適用可能だと理解された。

1.1.5.4. 物質粒子への適用

質量 $m$ の粒子が速度 $v$ で運動しているとき、運動量は $p = mv$ (相対論的効果が小さいと仮定)。すると

$$\lambda = \frac{h}{mv}$$

という式が得られる。これは粒子がもつ運動量に対応した波長(物質波長)とみなすことができる。

1.1.5.5. 古典系や量子系への示唆

• 粒子が古典的に振る舞うか、波動性が顕著になるかは、この $\lambda$ と系の典型的長さスケール(例えば粒子間距離や散乱断面積)を比較して判断できる
• 波長が十分小さいとき、粒子は古典的な質点モデルで扱っても誤差が小さい
• 波長が大きい場合、干渉やトンネル効果など量子的振る舞いを無視できなくなる

1.1.6 熱的ドブロイ波長

熱的ドブロイ波長とは、 粒子の運動が主に古典的に振る舞うか、 あるいは量子的な波動性が顕著になるかを評価する目安の長さスケールを示す。 以下では、その導出手順をなるべく初学者向けに整理する。

1.1.6.1. 自由粒子の運動量とドブロイ波長

量子力学で粒子に割り当てられる波動性は、ドブロイの関係式

$$\lambda = \frac{h}{p}$$

によって定義される。ここで

• $h$ はプランク定数
• $p$ は粒子の運動量

1.1.6.2. 熱運動による典型的運動量

温度$T$の古典的系を考えると、粒子の運動エネルギー(運動量)の大きさは、マクスウェル=ボルツマン分布によって決まる。典型的な運動量のスケールを見積もるためには、

$$\frac{p^2}{2m} \sim k_B T$$

というオーダー評価を用いる。よって

$$p \sim \sqrt{2 m k_B T}$$

が「熱運動によって得られる典型的な運動量」となる。

1.1.6.3. 典型的運動量をドブロイ波長に当てはめる

この運動量をドブロイの式に代入すると、

$$\lambda_\mathrm{th} \,\sim\, \frac{h}{\sqrt{2 m k_B T}}$$

というスケールが得られる。これが量子波動性の大きさを示す、いわば「熱的ドブロイ波長」の概念的出発点になる。

1.1.6.4. 熱的ドブロイ波長のより正確な定義

ただし、より厳密には統計力学の手順を踏むと、$2\pi$ の因子などが導入される。最終的にしばしば用いられる定義は

$$\Lambda_\mathrm{th} \,=\, \sqrt{\frac{h^2}{2 \pi m k_B T}}$$

となる場合が多い。これは1粒子の運動量積分を量子化定数$h$で区切りながら計算し、相対的な誤差を含めて導入される補正因子(主に$2\pi$)から得られる。

1. マクスウェル=ボルツマン分布を量子統計へ繋げる際に、相空間の1セルを$(h)^{3}$(3次元の場合)の大きさで区切り、運動量積分を評価する
2. その結果、$2\pi$の補正がかかった$\Lambda_\mathrm{th}$が自然に定義される

1.1.6.5. 物理的意味

• 粒子同士の平均的な距離より$\Lambda_\mathrm{th}$が十分に小さいとき、粒子は古典粒子のように近似できる
• $\Lambda_\mathrm{th}$が相対的に大きい場合、波動重ね合わせなど量子力学的効果を考慮する必要が高まる

具体的には、室温付近($T \approx 300\,\mathrm{K}$) かつ重原子(質量$m$大きめ)の系では、$\sqrt{2\pi m k_B T}$が比較的大きくなるため、$\Lambda_\mathrm{th}$は原子間距離(数オングストローム)よりも小さくなる傾向がある。これにより量子的な波動性よりは古典的な粒子像が優勢になり、分子動力学の枠組み(古典近似)での解析が妥当になる。

ボルツマン定数と

以下では、単一粒子の平衡状態を統計力学的に扱う際の「運動量空間の取り扱い」と「量子論的な位相空間区画」のアイデアを組み合わせ、結果的に

$$\Lambda_\mathrm{th} = \sqrt{\frac{h^2}{2\pi m k_B T}}$$

という熱的ドブロイ波長の定義が導かれる流れを示す。最終的には

$$\frac{p^2}{2m} \sim k_B T$$

という関係も理解できるようにする。

---

1. ボルツマン定数とは何か

古典熱力学における理想気体の状態方程式

$$$$

を分子数 $N$ を用いる形

$$PV = N_{A}k_B T$$

で書き直すと

$$k_B = \frac{R}{N_A}$$

と定義される。

ここで

• $R$ は気体定数 ($\approx 8.314\,\mathrm{J/(mol\cdot K)}$)
• $N_A$ はアボガドロ数 ($\approx 6.022\times 10^{23}\,\mathrm{mol}^{-1}$)

要するに、$1\,\mathrm{K}$ の温度変化がどれだけのエネルギー変化に相当するかを示す“単粒子版の気体定数”とみなせ、 ボルツマン定数 $k_B$ は、熱力学的温度 $T$ とエネルギーの次元を結びつける定数となる。

---

2. 古典統計力学: 単一粒子の運動量分布

質量 $m$ の粒子が $3$次元空間を自由に運動しているとき、カノニカル分布(温度 $T$)に従う確率密度は、運動量ベクトル $\mathbf{p} = (p_x,\,p_y,\,p_z)$ の空間で

$$f(\mathbf{p}) \,\propto\, \exp\!\Bigl(-\,\beta\, \frac{p^2}{2m}\Bigr)$$

となる。ここで

• $p^2 = p_x^2 + p_y^2 + p_z^2$
• $\beta = 1/(k_B T)$

この指数因子が「運動エネルギー $\frac{p^2}{2m}$ が高い状態ほど確率が低い」ことを示している。

2.1 平均運動エネルギー

運動エネルギー $\frac{p^2}{2m}$ の平均は、等分配則により

$$\bigl\langle \tfrac{p^2}{2m} \bigr\rangle = \tfrac{3}{2} k_B T$$

となる。1自由度あたり $\frac{1}{2} k_B T$ が配分されるためで、3次元なら $\frac{3}{2} k_B T$。これをざっくり

$$\frac{p^2}{2m} \sim k_B T$$

と見なせば、熱運動のオーダー評価が得られる。

---

3. 量子論的相空間区画: $h^3$ のセル

古典的には運動量空間の体積要素 $d^3 p = dp_x\,dp_y\,dp_z$ は連続無限小として扱う。しかし、量子力学の枠では、1粒子あたりの実空間と運動量空間を合わせた「相空間」を

$$\Delta x\,\Delta p \sim h$$

程度の大きさで区切る考え方(パウリによる位相空間の1セル)を導入する。この $h$ はプランク定数であり、3次元であれば1セルの大きさは $h^3$ となる。

• 位相空間とは (座標, 運動量) の6次元空間
• 古典的には連続だが、量子効果で “1つの量子状態” が占める体積を $h^3$ とみなす

---

4. 単一粒子の平衡分配関数: 例として $Z_\text{trans}$

3次元の自由粒子が体積 $V$ 中を運動するとき、平衡分配関数の運動量積分部分 (翻訳運動の寄与) を

$$Z_\mathrm{trans} = \frac{1}{h^3} \int \exp\bigl(-\beta\,\tfrac{p^2}{2m}\bigr)\, d^3 p \;\times\; \text{(実空間積分)}$$

のように書く。実空間積分は単に体積 $V$ になるので、

$$Z_\mathrm{trans} = \frac{V}{h^3} \int_0^\infty \!\!\int \exp\bigl(-\beta\,\tfrac{p^2}{2m}\bigr)\, 4\pi\, p^2\, dp$$

(角度成分の積分で $4\pi$ が出る)。ここで上式の運動量積分を行うと、

$$\int_0^\infty \exp\!\Bigl(-\beta\, \tfrac{p^2}{2m}\Bigr)\,4\pi\,p^2\, dp = \bigl(2\pi m k_B T\bigr)^{3/2}$$

という結果が得られる。よって

$$Z_\mathrm{trans} = \frac{V}{h^3} \,\bigl(2\pi m k_B T\bigr)^{3/2}$$

とまとまる。

---

5. 熱的ドブロイ波長 $\Lambda_\mathrm{th}$ の定義

上式を

$$Z_\mathrm{trans} = \frac{V}{\Lambda_\mathrm{th}^3}$$

という形で書く場合が多い。ここで

$$\Lambda_\mathrm{th} \;=\; \frac{h}{\sqrt{2\pi m k_B T}}$$

を「熱的ドブロイ波長」と定義する。すなわち、

$$= \sqrt{\frac{h^2}{2\pi m k_B T}}.$$

これは、量子力学的に見た “1セル当たりの運動量分布” を考え、古典的なカノニカル分布とのつながりを確立したときに自然に導入される波長スケールである。

---

6. 意味するところ

• 粒子間距離が $\Lambda_\mathrm{th}$ より十分大きい系では、重なり(干渉)の影響が小さく、古典近似で扱える
• $\Lambda_\mathrm{th}$ が大きくなる低温や軽粒子(例: 電子)の系では量子効果が顕著になる
• 結果的に $$\frac{p^2}{2m} \sim k_B T$$ が統計力学的に理解できる: カノニカル分布から導出される運動量積分で、運動エネルギーの平均スケールが $k_B T$ となる

---

まとめ

1. ボルツマン定数 $k_B$ が温度とエネルギーを関連づける
2. カノニカル分布における運動量空間の積分に量子論的なセル $h^3$ を導入すると、翻訳運動の分配関数から $\Lambda_\mathrm{th} = \sqrt{\frac{h^2}{2\pi m k_B T}}$ が定義される
3. この $\Lambda_\mathrm{th}$ が粒子の平均自由度と量子効果を比較する指標となり、古典力学近似の成立範囲を見積もる
4. $\tfrac{p^2}{2m} \sim k_B T$ は、カノニカル分布でのエネルギー等分配則とも対応している

このようにして「熱的ドブロイ波長」は量子力学的視点から見た粒子1個あたりの「有効波動パッケージ」の広がりを示す量になり、そこから $\tfrac{p^2}{2m} \sim k_B T$ という古典的なエネルギー評価も自然に浮かび上がる。

補足

以下では「気体定数」「ボルツマン定数」「エネルギー」の相互関係を理解するために、理想気体の内部エネルギーまで踏み込み、なぜ $pV = nRT$ や $pV = Nk_B T$ の式が“エネルギー”と結びつくのかを示す。

1. 理想気体の状態方程式

理想気体では、圧力 $p$ と体積 $V$、物質量 $n$、温度 $T$ が

$$pV = nR\,T$$

という関係に従う。$R \approx 8.314\,\mathrm{J/(mol\cdot K)}$ は気体定数。さらに、1モルあたりの粒子数をアボガドロ数 $N_A$ とすると、

$$n = \frac{N}{N_A} \quad\Longrightarrow\quad pV = N k_B T \quad (\text{ただし } k_B = \frac{R}{N_A})$$

として書ける。ここでは $N$ は粒子数で、$k_B$ がボルツマン定数になる。

なぜこれがエネルギーを示唆するのか

上式だけを見ると、$pV$ が「圧力 $\times$ 体積」であり、物理量として“仕事”の次元と関係することはわかる。ただしこの式自体に「内部エネルギー」という言葉は出てこないため、次に理想気体の内部エネルギーを併せて考える。

---

2. 理想気体の内部エネルギー

1原子分子(単原子分子)の理想気体の場合、内部エネルギー $U$ はすべて粒子の平動運動エネルギーに対応し、

$$U = \frac{3}{2} n R T$$

となる。ここで $3/2$ は3自由度(3次元空間)における等分配則の結果。この式に $n = N / N_A$ を代入すると、

$$U = \frac{3}{2} \,\frac{N}{N_A} R\,T = \frac{3}{2} N k_B T.$$

つまり1粒子あたりの平均運動エネルギーは $\tfrac{3}{2} k_B T$ で、温度が1 K 上がるごとに粒子1個当たり $\tfrac{3}{2} k_B$ だけのエネルギーが増える計算になる。

---

3. $pV$ と内部エネルギーの関係

上記の内部エネルギーと状態方程式 $pV = N k_B T$ とを組み合わせると、

$$U = \frac{3}{2} N k_B T = \frac{3}{2} pV.$$

(単原子理想気体の場合) すなわち圧力と体積の積 $pV$ は、内部エネルギーの2/3に相当する。これは「温度 $T$ が上がることは、粒子の運動エネルギーが上がること」を意味し、$pV = Nk_B T$ の背後にエネルギーが潜んでいるとわかる。

• $pV$ が上昇するとき、粒子の運動エネルギー(内部エネルギー)が増えている
• 結果的に $T$ も上昇している

---

4. ボルツマン定数の役割

ボルツマン定数 $k_B$ は、温度 1 K あたりのエネルギー尺度を示す。1粒子単位で考えると

• $T$ が 1 K 上昇すると、$\Delta E \approx k_B$ (オーダー評価で、3自由度の場合はさらに3/2を掛ける)

というふうに「温度」と「エネルギー」を直接結びつけている。これは分子1つ1つに着目したスケール。モル単位では $R$ (気体定数) が同様の役割を果たす。結局、

$$k_B = \frac{R}{N_A}$$

であり、「1モル単位か、1粒子単位か」の違いのみで本質は同じ。つまり、

1. $pV = n R T$ はモル単位での物質量と温度を結びつけ
2. $pV = N k_B T$ は粒子数ベースで同じことを表し
3. 内部エネルギー $U$ の式 (例えば $\tfrac{3}{2} n R T$) と組み合わせることで、温度変化とエネルギー変化の関係が可視化される

---

5. 最終的な解釈

• $pV = n R T$ は「圧力と体積と温度」を結ぶ式だが、理想気体の内部エネルギーが運動エネルギーのみで表されることを考えると、これが“エネルギー論”とつながる
• 温度を 1 K 上げると、1粒子あたり約 $k_B$ のエネルギー(厳密には次元や自由度による係数がつく)が増えるという目安が得られる
• したがって、1 K の温度変化に対応するエネルギー量は $k_B$ で測れるため、「ボルツマン定数は熱力学的温度とエネルギーの次元を結びつける定数」という理解になる

補足２

以下では、単原子理想気体の内部エネルギー $U$ が

$$U = \tfrac{3}{2} n R T$$

となる理由を、統計力学と等分配則を使って導出する。

---

1. 単原子分子(理想気体)の運動エネルギー

1原子分子は、平動以外の自由度(振動、回転)を持たないと考える。したがって、粒子の運動エネルギーは平動運動のみで決まる。1粒子の平動エネルギーは

$$E_\mathrm{trans} = \tfrac{1}{2} m (v_x^2 + v_y^2 + v_z^2)$$

である。3つの速度成分 $(v_x,\,v_y,\,v_z)$ によってエネルギーが分担されている。

---

2. 等分配則(Equipartition Theorem)

古典統計力学の等分配則によれば、温度 $T$ の熱平衡状態では「1自由度あたり $\tfrac{1}{2} k_B T$ の平均エネルギーを持つ」。ここで自由度とは、速度成分や位置成分など運動エネルギーに現れる独立な2次形式の数を指す。

• 単原子分子の平動自由度は3 (x, y, z 方向)
• よって1粒子の平均平動エネルギーは $\tfrac{3}{2} k_B T$

すなわち、

$$\langle E_\mathrm{trans} \rangle = \tfrac{3}{2} k_B T.$$

---

3. 多粒子系の内部エネルギー

内部エネルギー $U$ は、系全体の粒子の総エネルギー(相互作用がない理想気体ならば平動エネルギーの総和)とみなせる。粒子数が $N$ 個あるとすると、

$$U = N \,\times\, \bigl(\tfrac{3}{2} k_B T\bigr).$$

モル単位 $n$ (モル数) で表すには、$N = n \, N_A$ として、

$$U = n\, N_A \,\times\, \bigl(\tfrac{3}{2} k_B T\bigr).$$

ところで $R = N_A k_B$ (気体定数とボルツマン定数の関係) を用いると、

$$U = \tfrac{3}{2}\,n\,N_A\,k_B\,T = \tfrac{3}{2} \,n\,R\,T.$$

これが単原子分子(理想気体)における内部エネルギーの基本式

$$U = \tfrac{3}{2} n R T.$$

---

4. 3自由度の起源

なぜ3自由度なのかは、1原子分子が「x, y, z」の3方向の平動しかないから。回転運動や振動運動などを考えない(あるいはそれらが量子力学的に凍結されている)とき、3つの平動自由度でエネルギーを分担する。結果、運動エネルギーの期待値は $\frac{3}{2} k_B T$ / 1粒子 となる。

---

5. 意味するところ

1. 単原子分子は、いわゆる「核と電子殻からなる1つの球」として振る舞い、他の内部自由度が表に出てこない
2. 温度1 K 上昇あたり、1粒子の平均エネルギーは $\tfrac{3}{2} k_B$ 増加
3. 1モル(約 $6.02\times 10^{23}$ 個)あたりでは $\tfrac{3}{2} R$ [J/K] だけエネルギーが増える

これにより、気体の内部エネルギーが “温度しか依存しない” (体積や圧力には直接依存しない) 理由もわかる。単原子分子の理想気体では分子間相互作用を無視できるため、内部エネルギーは純粋に熱運動(平動)だけで決まるからである。

補足３単原子理想気体の内部エネルギー導出

等分配則の導出

カノニカル分布と分配関数

古典統計力学では、熱平衡状態にある系のミクロ状態の確率分布はカノニカル分布に従う。系のハミルトニアンを $H(\mathbf{p}, \mathbf{q})$、逆温度を $\beta = \frac{1}{k_B T}$ とすると、位相空間内の点 $(\mathbf{p}, \mathbf{q})$ に系が存在する確率密度は：

$$\rho(\mathbf{p}, \mathbf{q}) = \frac{1}{Z} e^{-\beta H(\mathbf{p}, \mathbf{q})}$$

ここで $Z$ は分配関数であり：

$$Z = \int e^{-\beta H(\mathbf{p}, \mathbf{q})} d\mathbf{p} d\mathbf{q}$$

二次形式のエネルギーと等分配則

ハミルトニアンが、ある座標（または運動量）$x$ に関して二次形式 $H = \frac{1}{2} \alpha x^2 + \ldots$ ($\alpha$ は定数) の項を含む場合、その座標に関する平均エネルギーを計算する。

平均エネルギーは期待値として：

$$\langle \frac{1}{2} \alpha x^2 \rangle = \frac{1}{Z} \int \frac{1}{2} \alpha x^2 e^{-\beta H(\mathbf{p}, \mathbf{q})} d\mathbf{p} d\mathbf{q}$$

$x$ に関する部分だけを取り出すと：

$$\langle \frac{1}{2} \alpha x^2 \rangle = \frac{1}{Z_x} \int \frac{1}{2} \alpha x^2 e^{-\beta \frac{1}{2} \alpha x^2} dx$$

ここで $Z_x = \int e^{-\beta \frac{1}{2} \alpha x^2} dx$ は $x$ に関する分配関数の一部である。

この積分は次のように計算できる：

$$\begin{aligned} \langle \frac{1}{2} \alpha x^2 \rangle &= \frac{1}{Z_x} \frac{1}{2} \alpha \int x^2 e^{-\beta \frac{1}{2} \alpha x^2} dx \\ &= \frac{1}{Z_x} \frac{1}{2} \alpha \cdot \frac{1}{\beta \alpha} \int \beta \alpha x^2 e^{-\beta \frac{1}{2} \alpha x^2} dx \end{aligned}$$

部分積分を用いると：

$$\int \beta \alpha x^2 e^{-\beta \frac{1}{2} \alpha x^2} dx = -2 \int x \frac{d}{dx}(e^{-\beta \frac{1}{2} \alpha x^2}) dx = 2 \int e^{-\beta \frac{1}{2} \alpha x^2} dx = 2Z_x$$

したがって：

$$\langle \frac{1}{2} \alpha x^2 \rangle = \frac{1}{Z_x} \frac{1}{2} \alpha \cdot \frac{1}{\beta \alpha} \cdot 2Z_x = \frac{1}{\beta} = k_B T$$

この結果より、二次形式の各自由度は平均エネルギー $\frac{1}{2}k_B T$ を持つことが導かれる。なぜなら、二次形式の項は通常 $\frac{1}{2} \alpha x^2$ の形をしており、これが $\frac{1}{2}k_B T$ の平均エネルギーを持つためである。

等分配則と自由度

以上の導出から、古典統計力学の等分配則によれば、温度 $T$ の熱平衡状態において、系の各二次の自由度は平均エネルギー $\frac{1}{2}k_B T$ を持つ。ここで $k_B$ はボルツマン定数である。

単原子気体の自由度

単原子気体の場合、各原子は3次元空間内を運動するため、並進運動に関する自由度は3つある。すなわち、$x$, $y$, $z$ 方向の運動エネルギーがそれぞれ独立した自由度となる。

単原子分子の運動エネルギーは次のように表される：

$$E = \frac{1}{2}mv_x^2 + \frac{1}{2}mv_y^2 + \frac{1}{2}mv_z^2$$

ここで、$m$ は分子の質量、$v_x$, $v_y$, $v_z$ はそれぞれの方向の速度成分である。

内部エネルギーの計算

等分配則によれば、各自由度は平均 $\frac{1}{2}k_B T$ のエネルギーを持つため、1つの分子あたりの平均エネルギーは：

$$\langle E \rangle = 3 \times \frac{1}{2}k_B T = \frac{3}{2}k_B T$$

$N$ 個の分子からなる系の全内部エネルギー $U$ は：

$$U = N \times \frac{3}{2}k_B T = \frac{3}{2}N k_B T$$

モル数を使った表現

アボガドロ数 $N_A$ を用いて、$n$ モルの気体に含まれる分子数は $N = n N_A$ である。また、気体定数 $R$ とボルツマン定数 $k_B$ の間には $R = N_A k_B$ の関係がある。

これらを用いると：

$$U = \frac{3}{2}N k_B T = \frac{3}{2}n N_A k_B T = \frac{3}{2}n R T$$

よって、$n$ モルの単原子理想気体の内部エネルギーは $U = \frac{3}{2}n R T$ となる。

カノニカル分布の導出と意味

統計力学の基本的考え方

統計力学では、多数の粒子からなるマクロな系の熱力学的性質を、ミクロな粒子の運動から理解することを目指す。系の状態は位相空間内の点 $(\mathbf{p}, \mathbf{q})$ で表される。ここで $\mathbf{p}$ は全粒子の運動量ベクトル、$\mathbf{q}$ は全粒子の位置ベクトルである。

カノニカル分布の導出

孤立系と小正準集団

まず、全体系（宇宙）を考える。これは孤立系であり、エネルギー $E_{tot}$ は保存される。この系の可能な微視的状態（ミクロ状態）の数を $\Omega(E_{tot})$ とする。

ここで、我々が関心を持つ小さな系（部分系）と、それを取り巻く大きな環境（熱浴）に分ける：

$$E_{tot} = E + E_{bath}$$

ここで、$E$ は系のエネルギー、$E_{bath}$ は熱浴のエネルギーである。

エントロピーと確率

熱浴のエントロピー $S_{bath}$ は、熱浴の可能な微視的状態の数 $\Omega_{bath}$ と以下の関係がある：

$$S_{bath} = k_B \ln \Omega_{bath}(E_{bath}) = k_B \ln \Omega_{bath}(E_{tot} - E)$$

系がエネルギー $E$ を持つ確率 $P(E)$ は、対応する熱浴の微視的状態の数に比例する：

$$P(E) \propto \Omega_{bath}(E_{tot} - E)$$

熱浴は系よりも非常に大きいため、$E_{bath}$ を $E$ で展開できる：

$$S_{bath}(E_{tot} - E) \approx S_{bath}(E_{tot}) - \frac{\partial S_{bath}}{\partial E_{bath}} E$$

熱力学の関係式より、$\frac{\partial S_{bath}}{\partial E_{bath}} = \frac{1}{T}$ となる。ここで $T$ は温度である。

したがって：

$$\begin{aligned} P(E) &\propto \Omega_{bath}(E_{tot} - E) \\ &= \exp\left(\frac{S_{bath}(E_{tot} - E)}{k_B}\right) \\ &\approx \exp\left(\frac{S_{bath}(E_{tot})}{k_B} - \frac{E}{k_B T}\right) \\ &\propto \exp\left(-\frac{E}{k_B T}\right) \end{aligned}$$

これにより、系のエネルギー $E$ に対する確率分布が指数関数形式で表されることがわかる。

ミクロカノニカル分布からカノニカル分布へ

系のハミルトニアン $H(\mathbf{p}, \mathbf{q})$ は、系の位相空間上の点 $(\mathbf{p}, \mathbf{q})$ におけるエネルギーを与える。したがって、系が位相空間内の特定の点 $(\mathbf{p}, \mathbf{q})$ に存在する確率密度は：

$$\rho(\mathbf{p}, \mathbf{q}) = \frac{1}{Z} \exp\left(-\frac{H(\mathbf{p}, \mathbf{q})}{k_B T}\right) = \frac{1}{Z} e^{-\beta H(\mathbf{p}, \mathbf{q})}$$

ここで、$\beta = \frac{1}{k_B T}$ は逆温度、$Z$ は分配関数と呼ばれる規格化定数である：

$$Z = \int e^{-\beta H(\mathbf{p}, \mathbf{q})} d\mathbf{p} d\mathbf{q}$$

この分布は、系が熱浴と接触して熱平衡にある場合の確率分布であり、カノニカル分布（正準分布）と呼ばれる。

カノニカル分布の物理的意味

カノニカル分布は、熱平衡状態における系の確率分布を表す。この分布に従うことで、系は以下の特性を持つ：

1. 低エネルギー状態ほど出現確率が高い
2. 高温になるほど、高エネルギー状態の出現確率が増加する
3. 分配関数 $Z$ から系の熱力学量（自由エネルギー、エントロピー、内部エネルギーなど）を計算できる

具体的には、以下の関係式が成り立つ：

• 自由エネルギー: $F = -k_B T \ln Z$
• 内部エネルギー: $U = -\frac{\partial \ln Z}{\partial \beta} = \langle H \rangle$
• エントロピー: $S = k_B \beta \langle H \rangle + k_B \ln Z$

これらの関係式を用いることで、ミクロな量子力学・統計力学からマクロな熱力学的性質を導出できる。

補足 統計熱力学の基本概念: エントロピーと確率

1. 熱浴のエントロピーと微視的状態数の関係

1.1 ボルツマンの関係式

統計力学において、エントロピー $S$ は系の微視的状態数 $\Omega$ と以下の関係がある:

$$S = k_B \ln \Omega$$

ここで $k_B$ はボルツマン定数（$1.380649 \times 10^{-23}$ J/K）である。

この関係式はボルツマンの関係式と呼ばれ、統計力学と熱力学をつなぐ重要な式である。エントロピーという巨視的な熱力学量が、微視的な状態数と対数関係にあることを示している。

1.2 ボルツマンの関係式の導出

この関係式をボルツマンの関係式と呼び、 熱力学的エントロピーの性質と微視的状態数の性質を結びつけることで以下に導出する。

エントロピー $S$ と微視的状態数 $\Omega$ の間に何らかの関数関係

$$S = f(\Omega)$$

があると仮定し、この関数 $f$ を解くことを考える。

熱力学において、独立した二つの系 A と B があるとき、全体のエントロピーは各部分のエントロピーの和になると考えられる。

$$S_{A+B} = S_A + S_B$$

一方、独立した二つの系 A と B の微視的状態数をそれぞれ $\Omega_A$ と $\Omega_B$ とすると、組み合わせの原理より、複合系の微視的状態数は積になると考えられる。

$$\Omega_{A+B} = \Omega_A \times \Omega_B$$

上記の性質を組み合わせると、関数 $f$ は次の関数方程式を満たす必要がある:

$$f(\Omega_A \times \Omega_B) = f(\Omega_A) + f(\Omega_B)$$

この関数方程式を満たす連続関数を求める。

$$\left\{ \begin{align*} \Omega_A = e^x\\ \Omega_B = e^y \end{align*} \right.$$

とおくと、

$$\Omega_A \times \Omega_B = e^x \times e^y = e^{x+y}$$

となる。

関数 $g(x) = f(e^x)$ を定義すると、元の関数方程式は:

$$g(x+y) = g(x) + g(y)$$

となる。これはコーシーの関数方程式であり、その連続解は:

$$g(x) = Cx$$

ここで $C$ は定数である。元の関数に戻すと:

$$f(e^x) = Cx$$

$\Omega = e^x$ より $x = \ln \Omega$ なので:

$$f(\Omega) = C \ln \Omega$$

熱力学との整合性から定数 $C$ はボルツマン定数 $k_B$ に一致する必要があり:

$$S = f(\Omega) = k_B \ln \Omega$$

これがボルツマンの関係式である。

1.3 物理的解釈

この式の物理的解釈:

• 微視的状態数 $\Omega$ が大きいほど、系は多くの異なる微視的配置をとれる
• 多くの配置が可能であるほど、系はより乱雑（無秩序）である
• エントロピーは無秩序の尺度として解釈できる
• 対数関数の使用により、非常に大きな数（典型的には $10^{23}$ オーダー）を扱いやすいスケールに変換できる
• 熱力学第二法則（孤立系ではエントロピーは増加する傾向がある）は、系が最も確率の高い状態（＝微視的状態数が最大の状態）へ向かうという統計的傾向として理解できる

熱浴についても同様に、そのエントロピー $S_{bath}$ は熱浴の微視的状態数 $\Omega_{bath}$ の対数に比例する:

$$S_{bath} = k_B \ln \Omega_{bath}(E_{bath})$$

ここで $E_{bath}$ は熱浴のエネルギーである。全系のエネルギー $E_{tot}$ が保存されるとすると、系がエネルギー $E$ を持つとき、熱浴のエネルギーは $E_{bath} = E_{tot} - E$ となる。従って:

$$S_{bath} = k_B \ln \Omega_{bath}(E_{tot} - E)$$

2. 系のエネルギー確率分布

系がエネルギー $E$ を持つ確率 $P(E)$ は、次の基本原理から導かれる:

等確率の原理: 孤立系において、すべての可能な微視的状態は等しい確率で出現する。

系と熱浴の合計系は孤立系と考えられるため、系がエネルギー $E$ を持つ確率は、そのときに実現可能な微視的状態の総数に比例する。系がエネルギー $E$ を持つとき、熱浴のエネルギーは $E_{bath} = E_{tot} - E$ となり、熱浴の取りうる微視的状態数は $\Omega_{bath}(E_{tot} - E)$ である。

したがって:

$$P(E) \propto \Omega_{bath}(E_{tot} - E)$$

エントロピーの定義を使って書き換えると:

$$P(E) \propto \exp\left(\frac{S_{bath}(E_{tot} - E)}{k_B}\right)$$

3. 熱力学的関係式 $\frac{\partial S}{\partial E} = \frac{1}{T}$ の導出

この関係式は熱力学の基本関係から導出できる。

まず、熱力学第一法則を考える:

$$dE = TdS - PdV + \mu dN + ...$$

系の体積 $V$ と粒子数 $N$ が一定の場合、この式は簡略化される:

$$dE = TdS$$

これを変形すると:

$$\frac{\partial S}{\partial E} = \frac{1}{T}$$

この関係式は、エネルギーの微小変化に対するエントロピーの変化率が、温度の逆数に等しいことを示している。

物理的な意味:

• 高温では、エネルギーの増加に対するエントロピーの増加は小さい
• 低温では、わずかなエネルギー増加でもエントロピーは大きく増加する

熱浴についても同様に:

$$\frac{\partial S_{bath}}{\partial E_{bath}} = \frac{1}{T}$$

この関係を使うと、熱浴のエントロピーを系のエネルギー $E$ の関数として展開できる:

$$S_{bath}(E_{tot} - E) \approx S_{bath}(E_{tot}) - \frac{\partial S_{bath}}{\partial E_{bath}} E$$ $$S_{bath}(E_{tot} - E) \approx S_{bath}(E_{tot}) - \frac{E}{T}$$

これを確率分布の式に代入すると、ボルツマン分布が導出される:

$$P(E) \propto \exp\left(\frac{S_{bath}(E_{tot}) - E/T}{k_B}\right) \propto \exp\left(-\frac{E}{k_B T}\right)$$

この式は、平衡状態における系のエネルギー分布を記述する基本的な関係式である。

演習問題

Q1: MDの概要 (基本問題)

分子動力学が「原子・分子の運動を時間発展で追う手法」であると言われる理由を1〜2行で述べよ

A1: 解答/解説

• 粒子の運動方程式を数値的に解き、時々刻々の位置と速度の変化から現象を再現する
• 時間発展を追うことで、構造・物性・相転移など動的な性質を調べることができる

Q2: 応用分野 (基本問題)

「これを学ぶと何ができるようになるのか？」という実例を5つ以上挙げよ

A2: 解答/解説

• 1. 材料物性の評価 (拡散係数、機械的強度)
• 2. バイオ分野 (タンパク質の立体構造変化、酵素反応部位のダイナミクス)
• 3. 化学反応の経路探索 (触媒表面上の反応)
• 4. ナノ材料の設計 (グラフェンやカーボンナノチューブの力学特性)
• 5. 液体・溶液の粘度や相変化 (凝固・融解など)

Q3: 古典近似 (基本問題)

古典近似では量子効果を一部無視するが、どのような状況でそれが妥当と考えられるか1つ例を挙げよ

A3: 解答/解説

• 室温付近で重い原子が主体の系 (例: アルゴン液体) では粒子の熱的波長が非常に短く、量子効果は小さい
• 高温・高自由度でブロックされない振動モードが多い固体や液体も古典扱いが成立する場合が多い

Q4: マクロ物性予測 (応用問題)

粘度や拡散係数などのマクロ物性をMDで評価する意義を説明し、実験との比較でどのような検証が可能か述べよ

A4: 解答/解説

• MDでは粒子の運動と相互作用から物性を理論的に求められる
• 実験値と一致すれば、モデル化や力場の妥当性が検証できる
• 一致しない場合は力場のパラメータ再検討や相互作用の不足が示唆される

Q5: 熱的ドブロイ波長の評価 (基本問題)

アルゴン原子 (質量約 40 u) を室温 (300 K) で扱うとき、熱的ドブロイ波長

$$\Lambda \sim \sqrt{\frac{h^2}{2\pi m k_B T}}$$

を数値的に求めよ。単位系をSIに変換して計算し、結果をアルゴン間の代表的な原子間距離 (約 3.4 A) と比較して、古典近似がどの程度妥当か考察せよ。

なお、

• プランク定数 $h \approx 6.626\times 10^{-34}\,\mathrm{J\cdot s}$
• ボルツマン定数 $k_B \approx 1.381\times 10^{-23}\,\mathrm{J/K}$
• アルゴンの質量 $m \approx 6.63\times 10^{-26}\,\mathrm{kg}$ (40 u に相当)
• 室温 $T = 300\,\mathrm{K}$

とする。

A5: 解答/解説のポイント

• 計算した $\Lambda$ (単位 m) を Aに換算し、3.4 A と比較
• $\Lambda$ が原子間距離より十分に小さければ量子効果は小さいとみなせる

---

Q6: 拡散係数の数値積分 (応用問題)

速度自己相関関数が指数的減衰

$$\langle \mathbf{v}(0)\cdot\mathbf{v}(t)\rangle = \langle v^2 \rangle \exp\bigl(-t/\tau\bigr)$$

で表されると仮定する。拡散係数 $D$ が

$$D = \int_0^\infty \langle \mathbf{v}(0)\cdot\mathbf{v}(t)\rangle \, dt$$

で定義されるとき、$\tau$ と $\langle v^2 \rangle$ を用いて $D$ を求めよ。

A6: 解答/解説のポイント

• 速度自己相関関数を $C(t) = \langle v^2\rangle \exp(-t/\tau)$ とおく

$$D = \int_0^\infty C(t)\,dt = \langle v^2 \rangle \int_0^\infty \exp(-t/\tau)\, dt$$

• 上式の指数積分結果は $\tau \langle v^2\rangle$
• 単位系(例: m^2/s, 1/s など)を適切に考慮すれば数値例にも対応可能

1.2 MDシミュレーションの基本的な流れ

• ワークフローの説明
  • 系の設定 → 数値積分による時間発展 → 結果解析・可視化 → （必要に応じて）大規模化・高度手法へ
• 実際の研究・産業応用での典型的ワークフロー
  • 基本との違いを学ぶ

Q1: 基本ワークフロー (基本問題)

系の設定 → 数値積分による時間発展 → 結果解析・可視化 → 大規模化・高度手法へ、という流れが一般的な理由を1〜2行で述べよ

A1: 解答/解説

• 粒子配置と力場が決まらなければシミュレーションを始められない
• 時間発展結果を解析し、必要ならスケールアップや追加手法を適用するのが自然な手順

Q2: 研究・産業応用のワークフロー (基本問題)

研究や産業用途でのワークフローは、基本流れとどう異なるか1つ挙げよ

A2: 解答/解説

• 材料設計やドラッグデザインの場合、実験との反復比較や大量のパラメータスキャンを行うので、自動化された大規模パイプライン(ハイスループット計算)を用いることが多い

Q3: 緩和と平衡 (基本問題)

数値積分による時間発展では、初期状態から目的とする熱力学的平衡までに緩和が必要となる。なぜか1行で述べよ

A3: 解答/解説

• 初期状態が不自然な配置や速度分布の場合が多く、系が平衡分布に到達するまである程度のステップを要する

Q4: シミュレーションの成功要因 (応用問題)

MDが信頼できる結果を得るために注意すべき点を2つ挙げ、それぞれ理由を簡潔に説明せよ

A4: 解答/解説

• 時間刻みの適切な選択
  • 大きすぎるとエネルギー暴走が起き、小さすぎると計算コストが肥大化する
• 力場の妥当性
  • 不適切な力場パラメータでは、平衡構造や物性が実系と大きくずれる

1.3 古典近似の前提

• 量子化学計算・モンテカルロ法との位置づけ比較
• モンテカルロ法の問題点と発展

Q1: 古典MDと量子化学計算 (基本問題)

分子軌道法や密度汎関数法(DFT)と比べて、古典MDが簡易である理由を1行で述べよ

A1: 解答/解説

• 電子レベルの波動関数を扱わず、原子核を点粒子(古典質点)として相互作用ポテンシャルのみを用いるため計算負荷が軽い

Q2: モンテカルロ法の違い (基本問題)

MDとモンテカルロ法(MC)では同じアンサンブルを生成できることがあるが、その主要な違いを1つ挙げよ

A2: 解答/解説

• MCは系のエネルギー変化を受容・拒否する確率過程によりサンプリングを行い、系の「状態遷移」しか扱わない
• MDは運動方程式を逐次的に解き、実際の「時間発展」を追う

Q3: MCの問題点と発展 (基本問題)

単純なMC法が大規模系や複雑なポテンシャルで苦戦する理由を簡潔に述べよ

A3: 解答/解説

• 大きな相互作用エネルギー勾配を乗り越える経路が探索しにくく、受容率が低下して収束が遅くなる場合が多い

Q4: 半古典的/ハイブリッド法 (応用問題)

古典と量子化学の中間的アプローチ(QM/MMなど)が提案されている。どのような系に有効か述べ、例を1つ挙げよ

A4: 解答/解説

• タンパク質内の活性部位だけを量子力学的に扱い、それ以外の大部分を古典力場で近似するなど
• 例: 酵素反応の活性中心 (QM) とタンパク質全体 (MM) を同時に扱うときに計算コストを大幅削減できる

第2章 理論的背景の再確認

本章では理論的背景の復習をして分子動力学の基礎を学びなおすことを目的とする

2.1 古典力学の再確認

• ハミルトニアン形式、ラグランジュ形式を再度把握

  • $F=ma$の運動方程式からのハミルトニアン形式、ラグランジュ形式の導出と方程式を理解
  • ハミルトニアン、ラグランジュ形式の運動方程式の解法
  • 振動子を初めとした典型的な例を二つ取り上げ、解析に解く方法を提示

• キーワード

• 「ハミルトニアン」: 運動エネルギー＋ポテンシャルエネルギー

$$H = T + V$$

• 「運動方程式」:

$$\dot{p_i} = -\frac{\partial H}{\partial q_i}, \; \dot{q_i} = \frac{\partial H}{\partial p_i}$$

Q1: 運動方程式からハミルトニアン形式への変換 (基本問題)

$F=ma$ からハミルトニアン形式への変換過程で登場する一般化運動量$p_i$をどのように定義するか述べよ

A1: 解答/解説

• $p_i = \frac{\partial L}{\partial \dot{q_i}}$
• ただし $L = T - V$ はラグランジュ関数で、$q_i$ は一般化座標

Q2: 1次元調和振動子のハミルトニアン (基本問題)

質量$m$、ばね定数$k$の1次元調和振動子のハミルトニアンを示せ

A2: 解答/解説

$$H = \frac{p^2}{2m} + \frac{1}{2} k x^2$$

Q3: ラグランジュ形式の運動方程式 (基本問題)

ラグランジュ形式ではオイラー=ラグランジュ方程式で運動を記述する。その式を1行で書け

A3: 解答/解説

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q_i}}\right) - \frac{\partial L}{\partial q_i} = 0$$

Q4: 振動子の解析解 (応用問題)

調和振動子のハミルトニアンから運動方程式を導出し、周期$T$がどのように与えられるか述べよ

A4: 解答/解説

• 運動方程式は $m \ddot{x} + k x = 0$
• 一般解は $x(t) = A \cos(\omega t) + B \sin(\omega t)$
• 角振動数 $\omega = \sqrt{\frac{k}{m}}$
• 周期 $T = \frac{2\pi}{\omega} = 2\pi \sqrt{\frac{m}{k}}$

2.2 統計力学・熱力学の再確認

• 温度・エネルギー・エンタルピー・エントロピー・自由エネルギーの関係を導出しその意味を学習
• ボルツマン分布、カノニカル分布等の位置づけを理解し、数式を導出しその意味を学習
• 分配関数を導出し数式の意味を理解
• アンサンブルとは何かを数式の上から理解
• キーワード
  • 「ボルツマン分布」
  • 「自由エネルギー」($F = U - TS$)、エンタルピー($H = U + PV$)
  • 「アンサンブル」(NVE, NVT, NPT)

Q1: ボルツマン分布 (基本問題)

ボルツマン分布の形式を示し、その物理的意味を1行で述べよ

A1: 解答/解説

$$P(\text{状態} \propto e^{-\beta E})$$

• エネルギーが低い状態ほど確率が高く、$\beta = 1/(k_B T)$ で温度依存が決まる

Q2: 自由エネルギーとエンタルピー (基本問題)

$F = U - TS$ と $H = U + PV$ の違いを簡潔に述べよ

A2: 解答/解説

• $F$ (ヘルムホルツ自由エネルギー) は一定温度・体積の系を扱うとき有用
• $H$ (エンタルピー) は一定圧力の系でのエネルギーバランスを扱うとき用いられる

Q3: カノニカル分布 (基本問題)

カノニカルアンサンブルにおける分配関数$Z$の定義を示せ

A3: 解答/解説

$$Z = \sum_{\text{全状態}} e^{-\beta E_i} \quad \text{(離散系の場合)}$$

あるいは積分系なら

$$Z = \int e^{-\beta E(q,p)} \, dq \, dp$$

Q4: 分配関数から物理量を導く (応用問題)

分配関数$Z$から内部エネルギー$\langle U \rangle$や圧力$P$を導く代表的な微分関係を2つ挙げよ

A4: 解答/解説

• $\langle U \rangle = -\frac{\partial \ln Z}{\partial \beta}$
• $P = k_B T \frac{\partial \ln Z}{\partial V}$ (あるいは $P = -\frac{\partial F}{\partial V}$)

2.3 分子動力学における熱力学の位置づけ

• 時間平均とアンサンブル平均の等価性（エルゴード仮説）を導出しこの前提から得られるメリットを説明
• 大規模系での「統計的揺らぎ」の説明
• その評価方法を説明し理解する

Q1: 時間平均とアンサンブル平均 (基本問題)

時間平均とアンサンブル平均の等価性を前提とするエルゴード仮説が何を意味するか1行で述べよ

A1: 解答/解説

• 長時間の観測(時間平均)が、多数の初期状態を無作為にとったときの平均(アンサンブル平均)に等しいとみなす

Q2: 大規模系での統計的揺らぎ (基本問題)

粒子数$N$が非常に大きい系ほど、マクロ物性の揺らぎはどのように振る舞うか1行で述べよ

A2: 解答/解説

• 相対的揺らぎは $1/\sqrt{N}$ 程度に減少し、系が大きくなるほど揺らぎは小さくなる

Q3: エルゴード仮説のメリット (基本問題)

エルゴード仮説を採用することでMDではどのように物性を計算できるか、1〜2行で述べよ

A3: 解答/解説

• 時系列データを十分に長く取れば、それを平均するだけでアンサンブル平均を近似できる

Q4: 系のサイズ拡大と平衡 (応用問題)

MDで巨大系を扱うとき、平衡化やエルゴード性に関してどんな問題が生じるか2つ挙げ、その理由を簡単に説明せよ

A4: 解答/解説

• 初期緩和に時間がかかる
  • 大規模系だと局所相互作用の伝播が遅く、系全体が平衡になるまで時間を要する
• すべての位相空間を網羅しにくい
  • ポテンシャル面が非常に複雑になり、エネルギー障壁を越えてサンプリングが困難になる場合がある

第3章 原子間ポテンシャルと力場

分子動力学で必要となるいくつかの考え方を導入する。 2章の理論的背景とこの章の内容を理解して初めて分子動力学の基本的なシミュレーションの理論的背景を構成する。

3.1 原子間相互作用の基本

• 分子動力学の“物理モデル”を形づくる力場（Force Field）の数式を示しその意味を学習

• 物理モデルとは何か？物理モデル以外にどのようなモデルがあるか？の説明

• ForceFieldとは何か？なぜこの考え方が必要であるか。なぜその仮説が成立すると考えられるのかの説明

• キーワード

  • 「Lennard-Jonesポテンシャル」

$$V_{\mathrm{LJ}}(r) = 4\epsilon \bigl[(\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6 \bigr]$$

• 「クーロン相互作用」:

$$\propto \frac{1}{r}$$

Q1: 力場の意義 (基本問題)

分子動力学において力場(ポテンシャル関数)を導入する理由を1行で述べよ

A1: 解答/解説

• 原子・分子間の相互作用エネルギーを数式で表し、運動方程式へ力として反映させるため

Q2: Lennard-Jonesポテンシャル (基本問題)

$V_{\mathrm{LJ}}(r) = 4 \epsilon \bigl[ (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6 \bigr]$ がどのような相互作用を表すか2行で述べよ

A2: 解答/解説

• 短距離での強い斥力 (12乗項) と、中距離での引力 (6乗項) をモデル化した相互作用
• 不活性分子や貴ガスなどの単純分子間力をよく近似できる

Q3: クーロン相互作用 (基本問題)

$\propto \frac{1}{r}$ で表されるクーロン相互作用は、どのような系で重要か1つ例を挙げよ

A3: 解答/解説

• イオン結晶やタンパク質内の電荷をもつ側鎖など、電気的に帯電した粒子同士が存在する系

Q4: 物理モデルと他のモデル (応用問題)

「物理モデル」とは何かを定義し、分子軌道法や統計的学習モデル(機械学習ポテンシャル)との違いを1点説明せよ

A4: 解答/解説

• 物理モデルは、観測結果を物理法則(相互作用ポテンシャル)で直接に説明する枠組み
• 機械学習ポテンシャルは多くのデータからブラックボックス関数を作るため、明示的な物理定義が薄い場合がある

3.2 代表的な力場(Force Field)

• AMBER, CHARMM, OPLS, GROMOSなどの概要

  • 数学的に各力場を理解し特性を知る

• 結合項（Stretch, Bend, Torsion）と非結合項（LJ, Coulomb）の考え方

  • 数学的に結合項は何であるかを理解しその特性を理解する
  • 数学的に非結合項は何であるかを理解しその特性を理解する

Q1: AMBER, CHARMM, OPLS, GROMOS (基本問題)

これらの力場が主に対象とする分子系を1〜2行でまとめよ

A1: 解答/解説

• AMBER, CHARMM, OPLSはタンパク質や核酸などのバイオ分子を主な対象とした力場
• GROMOSは汎用性が高く、生体高分子や有機分子にも対応可能

Q2: 結合項 (基本問題)

Stretch, Bend, Torsionといった結合項はどのような数式形式が一般的か、1つ例を挙げよ

A2: 解答/解説

• 例: 伸縮 (Stretch) の場合

$$E_\text{bond} = \frac{1}{2} k_b (r - r_0)^2$$

Q3: 非結合項 (基本問題)

LJやCoulomb以外に非結合相互作用を表すモデルがある場合、1例挙げよ

A3: 解答/解説

• 水素結合を明示的に取り込む補正項 (HBond項) が導入される力場もある

Q4: 力場のパラメータ化 (応用問題)

力場パラメータは実験データや量子化学計算によって決定されるが、その過程に生じる不確かさをどう扱うか簡潔に述べよ

A4: 解答/解説

• パラメータフィッティングでは複数のターゲット(結合長、スペクトル、エネルギー)を同時に最小二乗などで最適化
• 不確かさは誤差推定やクロスバリデーションで把握し、パラメータの選択範囲を検討する

3.3 応用への注意点

• 生体・材料・高分子での適用力場の違い
• ポラリザブル力場、金属間ポテンシャルなど

Q1: 生体・材料・高分子での力場の違い (基本問題)

生体高分子と金属材料などで使用する力場が異なる理由を1行で述べよ

A1: 解答/解説

• 取り扱う結合様式や電子の自由度、相互作用(メタリックボンドなど)が異なり、ポテンシャル形式も異なるため

Q2: ポラリザブル力場 (基本問題)

ポラリザブル力場が従来の力場に比べてどのような改良点をもたらすか2行で述べよ

A2: 解答/解説

• 分子間の電荷再分配をモデルに取り込み、環境や相手分子に応じて誘起双極子を変化させる
• 静電相互作用の精度が向上し、特に水や高誘電系での再現性が高まる

Q3: 金属間ポテンシャル (基本問題)

EAM(Embedded Atom Method)などが用いられるが、どのような特徴があるか1行で述べよ

A3: 解答/解説

• 原子の電子密度に基づき、埋め込みエネルギーを導入することで金属結合の集団的性質を再現する

Q4: 力場選択の指針 (応用問題)

バイオ系、金属系、高分子系などで力場を選ぶとき、何を基準に選定すべきか2つ挙げよ

A4: 解答/解説

• 対象分子の結合様式や電荷状態を適切に扱えるか
• 実績(公開された検証データ)やサポートされる分子種の範囲を確認し、研究目的に合った力場を用いる

---

第4章 分子動力学シミュレーションの実装

実際にシミュレーションを実装するにあたり、 理論と実装との間にあるギャップをこの章で埋める

4.1 数値積分法の実際（Verlet法の使い方）

• 運動方程式を数値的に解く際の実装要点を押さえる
• いくつかの数値積分手法とその特性を理解する
• なぜVerlet/Velocity-Verlet方が適しているかの説明とその意義を学習する
• キーワード
  • 「Verlet/Velocity-Verlet」: 実装が単純でエネルギー保存に優れ、MDでよく用いられる
  • 時間刻み$\Delta t$の選択指針

Q1: Velocity-Verletのアルゴリズムフロー(基本問題)

1次元系で、位置$x$、速度$v$、加速度$a$を用いるVelocity-Verlet法の更新式(1ステップ分)を列挙せよ。

A1: 解答/解説

• Velocity-Verletでは、まず速度の半ステップ更新 → 位置の更新 → 力(加速度)計算 → 残りの半ステップ分の速度更新 という流れになります。
• 具体的には次の式で表されます。

1. $v\left(t + \frac{\Delta t}{2}\right) = v(t) + \frac{1}{2} a(t)\,\Delta t$
2. $x(t + \Delta t) = x(t) + v\left(t + \frac{\Delta t}{2}\right) \Delta t$
3. $a(t + \Delta t)$ を新しい位置 $x(t + \Delta t)$ から計算する
4. $v(t + \Delta t) = v\left(t + \frac{\Delta t}{2}\right) + \frac{1}{2} a\left(t + \Delta t\right)\Delta t$

Q2: エネルギー保存性の確認(基本問題)

小さな系(例: 2粒子のLennard-Jonesポテンシャル系)で、オイラー法とVelocity-Verlet法を比較したときのエネルギー推移の違いを、概念的に述べよ。

A2: 解答/解説

• オイラー法は逐次的な誤差の蓄積が大きく、MDシミュレーションでは全エネルギー(運動エネルギー+ポテンシャルエネルギー)が発散・収束する傾向が高い。
• 一方、Velocity-Verletはシンプレクティック積分の一種であり、長時間シミュレーションでもエネルギー保存が良好。
• 実際の実装では、時間刻みを大きくするとどちらも誤差は増えるが、比較するとVelocity-Verletがエネルギー保存に優れる。

Q3: 時間刻み$\Delta t$の選択基準(基本問題)

分子動力学で扱う系の中に高周波振動(例: 水素結合の伸縮モード)が存在するとき、$\Delta t$ はどう設定すべきか？理由も述べよ。

A3: 解答/解説

• 一般に、系が持つ最も速い振動(高周波モード)を十分にサンプリングできるように、周期の10分の1〜20分の1程度の$\Delta t$を設定するのが目安。
• 水素結合の伸縮など強い結合は振動周期が短いので、$\Delta t$を小さくしないと数値的不安定・エネルギー暴走が起きやすい。

Q4: Verlet積分の拡張(応用問題)

Leapfrog-Verlet法、Gear predictor-corrector法など、Velocity-Verlet以外の手法が使われる場合がある。 (1) MDでVelocity-Verletが選択される主な利点を2つ挙げよ。 (2) Gear法が好まれるケースはどのような状況か考察せよ。

A4: 解答/解説

• (1) Velocity-Verletの利点

1. 実装が簡単である。
2. エネルギー保存特性が良好（シンプレクティック積分）。速度と位置を同時に扱いやすい。

• (2) Gear法が好まれるケース
• 位置・速度・加速度などの高次の微分を扱い、予測・補正のステップを挟むため、高次精度を狙う場合や、制約条件が多い系で特定の安定性が必要な場合に用いられる。
• ただし、単純なMDシミュレーションだとVelocity-Verletでも十分な精度が得られるため、計算コストとの兼ね合いで選択される。

---

4.2 境界条件と系のセッティング

• 周期境界条件をはじめどのような境界条件設定が存在するか
• それらの境界条件の間の特性の違いを理解し適用性を理解する
• PBCとは何か、どんな特性があるかを理解する
• クリスタル構造、液体、気体、など各種状態の置ける初期状態設定の考え方を学ぶ
• 相転移を考える際の初期状態設定の考え方を学ぶ
• キーワード
• 「Periodic Boundary Condition (PBC)」
• クリスタル構造・液体などに適用する際の初期配置と平衡化
• 演習問題

Q1: 周期境界条件 (PBC) の概念(基本問題)

PBC(Periodic Boundary Condition)を導入する主な目的を1?2行で述べよ。

A1: 解答/解説

• 有限サイズのシミュレーションセルでも「無限系」に近い挙動を模擬するため。
• 系の境界面での人工的な効果(表面効果)を最小化するために用いられる。

Q2: イメージセルの取り扱い(基本問題)

x方向に長さ$L_x$のシミュレーションボックスを想定する。位置座標$x_i$がボックス外に出た場合($x_i > L_x$)、どのようにしてPBC下で座標を再計算するか、式で示せ。

A2: 解答/解説

• 単純には「モジュロ演算」を用いる。

$$x_i^\prime = x_i \mod L_x$$

ただしプログラム実装では「floor」や「fmod」関数を使う場合が多い。

Q3: 初期構造の構築方法(基本問題)

密度が既知の液体をシミュレーションしたいとき、初期配置をどう生成するか？ 1つの具体例を挙げ、メリット・デメリットを述べよ。

A3: 解答/解説

• 一例として「単純立方格子（SC格子）に分子を整列して配置し、緩やかに加熱・緩和して液体状態へ持っていく」方法がある。
• メリット: 初期配置の作り方が単純で自動生成が容易。
• デメリット: 人工的な初期配列なので、序盤の緩和(Equilibration)に時間がかかることが多い。

Q4: 固液共存系の境界条件(応用問題)

固体と液体の界面を同時にシミュレーションしたい場合、PBCをどのように設定すれば良いか？ (1) 界面方向にPBCを適用しない(自由表面)パターンと、(2) すべての方向にPBCをかけるパターンの違いを考察せよ。

A4: 解答/解説

• (1) 界面方向にPBCを適用しない
• 実際の自由表面を再現できる反面、表面効果や表面張力に関する挙動は観測できる。
• ただし、有限サイズ効果を大きく受けやすい。
• (2) すべての方向にPBCをかける
• 固-液-固 と周期的に並ぶことになり、本当の自由表面は消失する。
• 代わりに「無限に反復した界面構造」のシミュレーションとなるため、界面近傍の影響が重複しうる。
• 実際の系と違う挙動を示す可能性があるが、単位表面積あたりの計算量は少なく済む。

---

4.3 長距離相互作用の扱い

• クーロン相互作用にカットオフだけを使うと精度が低下することを、数式を使って理解する
• Ewald法について、数学的に導出し、その意味を理解する
• 従来モデルに対してどう修正し、どう適用し精度が向上するか、数学的に物理的に理解する
• Particle Mesh Ewaldについて、数学的に導出し、その意味を理解する
• 従来モデルに対してどう修正し、どう適用し精度が向上するか数学的に物理的に理解する
• Reaction Field 法
• 従来モデルに対してどう修正し、どう適用し精度が向上するか数学的に物理的に理解する
• Multipole Expansion
• 従来モデルに対してどう修正し、どう適用し精度が向上するか数学的に物理的に理解する

Q1: Coulomb相互作用の収束性(基本問題)

なぜ単純に$\frac{1}{r}$で無限に広がるクーロン相互作用をリアル空間で切り取る(カットオフ)だけでは不十分な場合があるのか？

A1: 解答/解説

• クーロン相互作用は遠距離でも減衰が遅く、カットオフを設けるとエネルギーや力の計算に大きな誤差が出る。
• カットオフ付近で不連続が生じ、系全体のエネルギーが不安定になる場合もある。
• よって、Fourier空間(またはその他の手法)を使って長距離和を取り扱うEwaldサミングなどが必要。

Q2: Ewaldサミングの概念(基本問題)

Ewaldサミングで、「リアル空間項」と「逆空間(Fourier)項」に分割して和を取るのはなぜか？ 簡単に理由を述べよ。

A2: 解答/解説

• $\frac{1}{r}$の収束を高速化・安定化するために、$\mathrm{erf}$関数等を使って相互作用を2つの領域(近接はリアル空間で、遠隔はFourier空間で)に分けて評価する。
• こうすることで、両側の和がそれぞれ比較的早く収束し、全系としての無限和を効率的に計算できる。

Q3: Particle Mesh Ewald (PME)(基本問題)

PME法ではなぜ「ファストフーリエ変換(FFT)」が重要な役割を果たすのか？

A3: 解答/解説

• 逆空間での電荷分布を格子(メッシュ)に投影し、FFTを用いて畳み込み計算を高速に行うため。
• FFTを使うことで、$O(N \log N)$程度の計算量でCoulomb相互作用の長距離項を取り扱える。

Q4: Ewaldパラメータの最適化(応用問題)

Ewaldサミングでは「$\alpha$パラメータ」(分割関数に登場する遮蔽パラメータ)や「リアル空間カットオフ$r_c$」の最適値を設定する必要がある。 これらを適切に調整しないと、計算コストが増えたり誤差が大きくなったりするが、 (1) $\alpha$と$r_c$の関係を概念的に述べ、(2) 大まかな選び方を説明せよ。

A4: 解答/解説

• (1) $\alpha$と$r_c$の関係:
• $\alpha$が大きいほどリアル空間側で相互作用が急激に減衰し、$r_c$を小さく設定できるが、逆空間での分解能は上がり計算が増える。
• $\alpha$が小さいほどリアル空間での減衰が弱くなり、$r_c$を大きく取る必要があるが、逆空間は楽になる。
• (2) 大まかな選び方:
• 系のサイズや電荷分布に応じて、リアル空間と逆空間での計算コストがバランスするように$\alpha$と$r_c$を合わせる。
• 実際には経験的・数値実験的に最適値を見つけるケースが多い。

---

4.4 拘束手法

• 水分子のO-H結合長固定に使われるSHAKE法・RATTLE法などを導入する際の力学的背景を説明する
• バイオ分子や水分子モデル（例: TIP3P水）などを扱うときになぜこのテクニックが必要になるかを説明
• 結合長や角度を固定する拘束法の意味について説明
• SHAKEについて数学的に説明し、その数学的に物理的意味をに理解する
• RATTLEについて数学的に説明し、その数学的に物理的意味をに理解する

Q1: 拘束手法の導入目的 (基本問題)

分子動力学シミュレーションで拘束手法 (SHAKE, RATTLE など) を導入する主な目的を簡潔に述べよ

A1: 解答/解説

• 高周波振動 (例: O-H結合) を固定し、より大きな時間刻みでも数値的に安定して解くため
• 特定の結合長や角度を固定することで自由度を削減し、計算コストを下げる場合がある
• TIP3P水などでよく用いられる (O-H結合の拘束)

Q2: SHAKEとRATTLE (基本問題)

SHAKEとRATTLEの主な違いを1〜2行で述べよ

A2: 解答/解説

• SHAKEは位置座標の更新時に拘束条件を補正する方法
• RATTLEは速度 (加速度) の更新ステップでも拘束条件を考慮する拡張版

Q3: 水分子拘束の利点と注意点 (基本問題)

TIP3P水モデルなどでO-H結合を拘束するメリットと、拘束によって失われる可能性のある現象や自由度を1つずつ挙げよ

A3: 解答/解説

• メリット
• 高周波振動を省略できるため、時間刻みを大きく取っても数値的安定性が高い
• 失われる自由度
• O-H結合の振動モード (実際の分子振動) は再現できない

Q4: 大規模タンパク質系への応用 (応用問題)

大規模タンパク質をシミュレーションする際、主鎖に部分的拘束をかけることがある (1) なぜ部分的拘束を行うか (2) その利点と欠点をそれぞれ1つ挙げよ

A4: 解答/解説

• (1) 理由
• タンパク質の大きな変形を抑え、特定部位 (活性部位など) の挙動を重点的に解析する
• 利点
• 不要な揺らぎを抑えて計算負荷を下げる
• 欠点
• タンパク質本来の柔軟性を制限し、生理学的挙動を一部失う可能性がある

---

温度・圧力制御

• カノニカル(NVT)とは何か、その実現方法は何かを数学的に理解する
• 等圧等温(NPT)アンサンブルとは何か、その実現方法は何かを数学的に理解する
• 実際の制御手法を数学的に理解する
• 練習問題
• キーワード
• Thermostat: Nose-Hoover, Langevin, Andersen
• Barostat: Berendsen, Parrinello-Rahman

Q1: Thermostatの目的(基本問題)

NVTアンサンブルを実現するために、Thermostatが必要となる理由を簡潔に述べよ。

A1: 解答/解説

• 分子動力学の基本形(NVE)ではエネルギーが保存され、温度が一定に保たれるわけではない。
• 実際の系では外部との熱交換が可能と考えるため、Thermostatを導入し、温度を制御(一定)に保つ必要がある。

Q2: Nose-Hoover Thermostat(基本問題)

Nose-Hoover Thermostatの導入によって、運動方程式に付加される「仮想質量」や「熱浴変数」の役割を概念的に説明せよ。

A2: 解答/解説

• Nose-Hoover法では、系に「熱浴変数」(しばしば$\xi$など)を追加し、その変数が系の運動方程式に連動する形で温度を制御する。
• 仮想質量($Q$)は、この$\xi$の慣性に相当し、熱交換の緩慢さ(カップリングの強さ)を調整するパラメータ。
• 「$\dot{\xi}$」がゼロ付近を安定的に保つように温度が目標値に近づいていく。

Q3: Barostatの役割(基本問題)

NPTアンサンブル下で、密度や体積が変化しながら系が平衡化する挙動を再現するために、Barostatが必要となる理由を説明せよ。

A3: 解答/解説

• NPT系では圧力が一定になるように体積が変化できる状態を想定する。
• 分子動力学においては、シミュレーションボックスのサイズを動的に変化させる必要があり、この操作をBarostatが司る。
• Barostatを用いないとボックスサイズが固定され、等積条件(NVTやNVE)になってしまう。

Q4: 圧力制御の実装(応用問題)

Parrinello-Rahman Barostatを用いると、単にボックスの大きさだけでなく形状(斜方変形)も変わる可能性がある。 これが有効になる典型的なケース(物質系)を1つ挙げ、その理由を述べよ。

A4: 解答/解説

• 例: 結晶構造が高圧力下で相転移するような材料系(例: Si, Ge などの半導体結晶)。
• 剛体的な拡大・縮小だけでは再現できない格子定数の異方的変化や変形(歪み)をシミュレーションするために、ボックス形状そのものを変化させるBarostatが有効。
• Parrinello-Rahman法はボックスの形状行列を変化させるので、斜方晶・三斜晶への変形を自然に取り込む。

4.6 小規模系での練習例

• プログラムはTypescript/vite/three.js/等の方法を基本とする。python/c++での実装との違いを説明する
• アルゴン気体・Lennard-Jones系を例にした簡易MDプログラム
• どういうワークフロー、どういうプログラム構成にすべきかを提示する
• EquilibrationとProduction Runの流れ
• どういうワークフロー、どういうプログラム構成にすべきかを提示する

Q1: アルゴン気体の初期速度分布(基本問題)

アルゴン気体(N個のアルゴン原子)をシミュレーションする際、初期速度分布をどう設定するのが一般的か？理由も含めて述べよ。

A1: 解答/解説

• 一般的には、マクスウェル・ボルツマン分布(対象とする温度Tに対応)に従う乱数で速度を初期化する。
• 理由: シミュレーション開始直後から系が所望の温度周辺にいるため、平衡化時間が短縮される。

Q2: Equilibrationフェーズの目的(基本問題)

小規模のAr系で初期状態からMDを走らせたとき、なぜ「Equilibration」と「Production Run」を分けるのか？

A2: 解答/解説

• Equilibrationフェーズでは、系が熱平衡状態に達するまでの一過性(初期緩和)を含む。
• Production Runでは、平衡化した後の状態をサンプリングして、物性値の統計を取る。
• もし両者を分けない場合、初期緩和が統計に混ざり正しい平均を阻害する可能性がある。

Q3: シミュレーションステップ数(基本問題)

練習例としてAr原子100個を10psほどシミュレーションする場合、$\Delta t = 1$ fsとしたとき約何ステップ計算すればよいか？

A3: 解答/解説

• 1 fs = $1 \times 10^{-15}$ s
• 10 ps = $10 \times 10^{-12}$ s
• ステップ数 = $10\,\text{ps} / 1\,\text{fs} = 10{,}000$ステップ程度。

Q4: プログラム構成の設計(応用問題)

Web系(Typescript + three.js等)でMDを可視化する場合と、C++やPythonでファイル出力重視のオフライン計算をする場合、それぞれメリットとデメリットを1つずつ挙げよ。

A4: 解答/解説

• Web可視化のメリット:
  • リアルタイムにブラウザ上でインタラクティブ表示でき、教育用デモに最適。
• Web可視化のデメリット:
  • 大規模計算に向かない、ブラウザのパフォーマンス制限がある。
• C++/Pythonオフライン計算のメリット:
  • 並列化や最適化がしやすく、大規模計算が可能。
• C++/Pythonオフライン計算のデメリット:
  • 開発・実行環境が専用的になり、可視化には別ツールを組み合わせる必要があることが多い。

---

第5章 シミュレーションデータの解析・可視化

実際にシミュレーションを実装したときにその結果を どの様に可視化をして、どの様に評価をするかその考え方と手法を学ぶ

5.1 構造の評価

• 分子配置の特徴を解析し、物質構造を理解する
• キーワード
• 「Radial Distribution Function (RDF)」
• 配向秩序パラメータ
• 演習問題

5.2 物性・動力学的評価

• エネルギー成分解析（ポテンシャル・運動エネルギー）
• 自己拡散係数, Mean Square Displacement(MSD)
• 時系列データ解析（自己相関関数など）
• 演習問題

5.3 可視化ツール

• VMD, OVITOなどによる視覚化
• Python等でのスクリプト可視化例
• 演習問題

第6章 応用と高度な話題

5章までの基本的な考え方を拠り所とし、 さらに分子動力学を深化するための考え方を学びさらなる応用に進む。

6.1 フリーエネルギー計算

• 物性評価や反応経路解析において不可欠な“自由エネルギー”の定量的評価手法
• キーワード
• 「熱力学的積分 (Thermodynamic Integration; TI)」
• 「自由エネルギー摂動法 (FEP)」
• 「Umbrella Sampling / WHAM解析」
• 演習問題

6.2 高度なサンプリング手法

• Replica Exchange MD (温度交換法)
• Metadynamicsなど
• 演習問題

6.3 大規模計算・並列化

• 並列化の基本(MPI, GPU)
• LAMMPS, GROMACS, NAMD, AMBERなど主要MDソフトの並列手法の概観
• 演習問題

6.4 分野別応用例

• バイオ系（タンパク質-リガンド相互作用）
• 材料系（結晶、アモルファス、ナノ材料）
• 演習問題

付録

A. 学習用の参考文献リスト

• Daan Frenkel & Berend Smit, “Understanding Molecular Simulation”
• 分子シミュレーション全般（MDとMC法）を広範にカバーする定番。
• M. P. Allen & D. J. Tildesley, “Computer Simulation of Liquids”
• MDの理論背景をしっかり固めたい人向けの古典的名著。
• Andrew R. Leach, “Molecular Modelling: Principles and Applications”
• バイオ系を含む力場や分子モデリングの基礎を網羅。
• Rapaport, “The Art of Molecular Dynamics Simulation”
• コード例付きで実装を学びたい場合に。
• Mark E. Tuckerman, “Statistical Mechanics: Theory and Molecular Simulation”
• 統計力学を深く学びたい人向け。
• 日本語の書籍例
• 『分子シミュレーション入門』（化学同人）
• 『分子動力学シミュレーションの基礎』（東京大学出版会）

B. MDプログラムの簡易サンプル

• C++/Python実装例：
• Velocity-Verletアルゴリズム (Lennard-Jones系)
• Thermostat例: シンプルなLangevin温度制御

C. 関連ツールとWebリソース

• GROMACS, LAMMPS, NAMD, AMBERなど公式ドキュメントへのリンク
• Pythonライブラリ（MDAnalysis, MDtraj等）での解析例