第16章力と応力・構造最適化・第一原理分子動力学

前章までで，Kohn–Sham方程式を実際に解くための道具立て——基底関数，Brillouin域のサンプリング(第12章)，擬ポテンシャル(第15章)——が出揃った．すなわち，原子核の配置 $\{\bm{R}_I\}$ を与えれば，全エネルギー $E(\{\bm{R}_I\})$ を計算できるようになった．しかし現実の問題では，原子配置そのものが未知であることが多い．分子はどんな形で安定するのか，結晶の格子定数はいくらか，化学反応はどんな経路をたどり，障壁は何eVか，有限温度で原子はどう動くのか——これらはすべて「エネルギーの原子座標に関する微分」，すなわち力と応力を計算できて初めて答えられる問いである．

本章ではまず，力の計算を支えるHellmann–Feynmanの定理を導出し，有限基底に由来するPulay補正を明らかにする．次に力の「格子版」である応力テンソルを定義し，構造最適化のアルゴリズム(最急降下法・Newton法・BFGS法)，反応経路探索のNEB法，そして第一原理分子動力学(Verlet積分・Nosé–Hoover熱浴・Car–Parrinello法)へと進む．静的な電子状態理論が，この章で初めて「動く」応用へとつながる．次章では電子の側を非平衡に駆動する量子輸送を扱う．

この章で学ぶこと

ポテンシャルエネルギー面(PES)の概念と，極小・鞍点の物理的意味
Hellmann–Feynmanの定理の完全な導出と，有限基底で現れるPulay力
応力テンソルのひずみ微分による定義と，ビリアル定理との整合性
最急降下法・Newton法・BFGS法(準Newton法)の導出と使い分け
調和近似・力の定数行列・動的行列の導出，フォノン分散と動的安定性の判定
NEB法による遷移状態探索:力の射影というアイデア
Verlet法の導出とBorn–Oppenheimer分子動力学
Nosé–Hoover熱浴が正準分布を生む仕組み，Car–Parrinello法の拡張ラグランジアン
メタダイナミクスとQM/MM法の概観

本章でもHartree原子単位系($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$，第1章)を用いる．エネルギーの単位は Hartree(1 Hartree $= 27.211$ eV)，長さの単位は Bohr(1 Bohr $= 0.5292$ Å)である．本章で新しく登場する組立単位として，力は Hartree/Bohr(1 Hartree/Bohr $= 51.42$ eV/Å)，応力・圧力は Hartree/Bohr$^3$(1 Hartree/Bohr$^3$ $= 2.942\times 10^4$ GPa)，時間は原子単位時間(1 a.u. $= \hbar/E_{\mathrm{h}} = 2.419\times 10^{-2}$ fs，すなわち 1 fs $\simeq 41.34$ a.u.)で測る．

16.1 ポテンシャルエネルギー面(PES)

本書では一貫してBorn–Oppenheimer近似を仮定してきた．原子核の配置 $\{\bm{R}_I\}$($I = 1,\dots,M$)を止めて電子系の基底状態を解くと，全エネルギー

$$ E(\bm{R}_1, \bm{R}_2, \dots, \bm{R}_M) = E_{\mathrm{elec}}(\{\bm{R}_I\}) + \sum_{I \lt J} \frac{Z_I Z_J}{\abs{\bm{R}_I - \bm{R}_J}} \label{eq:16-pes-def} $$

が得られる．ここで $E_{\mathrm{elec}}$ は電子系の基底状態エネルギー，第2項は原子核間の古典的Coulomb反発であり，$Z_I$ は原子核 $I$ の電荷である．式\eqref{eq:16-pes-def}は $3M$ 個の座標変数を持つ一つの関数，すなわち $3M$ 次元空間の上に張られた「面」を定義する．これをポテンシャルエネルギー面(potential energy surface, PES)と呼ぶ．原子核の立場から見れば，$E(\{\bm{R}_I\})$ は原子核が感じるポテンシャルエネルギーそのものであり，本章で扱うすべての手法——構造最適化・遷移状態探索・分子動力学——は「この面の上をどう歩くか」という一つの問いの変奏である．

PESの幾何学で特に重要なのは，勾配(1階微分)とHessian(2階微分)である．まず勾配の符号を反転したものが原子核 $I$ に働く力である:

$$ \bm{F}_I = -\frac{\partial E}{\partial \bm{R}_I}. \label{eq:16-force-def} $$

次に，勾配がゼロになる点($\bm{F}_I = \bm{0}$ がすべての $I$ で成立する点)を停留点と呼ぶ．停留点の性質は，2階微分係数を並べた $3M \times 3M$ の対称行列

$$ H_{I\alpha, J\beta} = \frac{\partial^2 E}{\partial R_{I\alpha}\, \partial R_{J\beta}} \qquad (\alpha, \beta = x, y, z) \label{eq:16-hessian-def} $$

——Hessian行列——の固有値で分類できる(並進・回転に対応するゼロ固有値は除く)．

固有値がすべて正の停留点は極小点である．どの方向に動かしてもエネルギーが上がるから，これは分子や結晶の安定構造(または準安定構造)に対応する．最も深い極小が最安定構造，それ以外は準安定構造である．
負の固有値がちょうど1個の停留点は1次の鞍点である．1方向にだけエネルギーが下がる峠であり，化学反応の遷移状態に対応する．峠の高さ(反応物側の極小からのエネルギー差)が活性化エネルギーである．
2つの極小を結び，途中で鞍点を通る「谷底の道」を最小エネルギー経路(minimum energy path, MEP)と呼ぶ．反応座標とはこの経路に沿った道のりのことである(16.5節)．

図16.1 ポテンシャルエネルギー面の模式図．等高線はエネルギーの等しい配置の集合．極小点は(準)安定構造，1次の鞍点は遷移状態に対応し，力は勾配の逆向き，すなわち「谷を下る」向きを向く．破線は2つの極小を結ぶ最小エネルギー経路(MEP)．

PESの次元の高さには注意が要る．原子100個の系なら $3M = 300$ 次元であり，1次元あたりたった10点でサンプルしても $10^{300}$ 点——面全体を「地図にする」ことは絶望的に不可能である．だからこそ，局所情報(その場での $E$，勾配 $-\bm{F}$，可能ならHessian)だけを頼りに，極小点や鞍点へ効率よくたどり着くアルゴリズムが本章の主役になる．そしてその大前提が，「勾配そのものを電子状態計算から正確に・安価に取り出せること」である．次節でこれを保証する定理を証明する．

物理的意味:第一原理計算が出力する3つの微分量

密度汎関数理論に基づく計算は，配置 $\{\bm{R}_I\}$ を入力すると (i) 全エネルギー $E$，(ii) 各原子に働く力 $\bm{F}_I = -\partial E/\partial \bm{R}_I$，(iii) 単位胞に働く応力テンソル $\sigma_{\alpha\beta}$(16.3節)を出力する．この3つ組があれば，分子の平衡構造，結晶の格子定数と弾性定数，表面吸着サイトの安定性，反応経路と活性化エネルギー，有限温度の動力学までを，実験パラメータなしに予測できる．本章はこの「微分量の使い方」の章である．

16.2 Hellmann–Feynmanの定理

力 $\bm{F}_I = -\partial E / \partial \bm{R}_I$ を計算する最も素朴な方法は，原子を少し動かして全エネルギーを2回計算し，差分で微分を近似することである．しかしこの方法では，$3M$ 成分の力を得るために $3M$ 回(中心差分なら $6M$ 回)ものSCF計算が必要になり，しかも差分幅の選び方に起因する誤差が避けられない．これから示すHellmann–Feynmanの定理は，1回のSCF計算で得た波動関数(密度)だけから，全成分の力が解析的に計算できることを保証する，構造最適化と分子動力学の礎石である．

16.2.1 定理の主張と導出

まず一般的な形で述べる．ハミルトニアン $\hat{H}(\lambda)$ が何らかのパラメータ $\lambda$(あとで原子核座標 $R_{I\alpha}$ と読み替える)に依存し，その規格化された固有状態 $\Psi(\lambda)$ とエネルギー $E(\lambda)$ が

$$ \hat{H}(\lambda)\, \Psi(\lambda) = E(\lambda)\, \Psi(\lambda), \qquad \braket{\Psi(\lambda) | \Psi(\lambda)} = 1 \label{eq:16-hf-setup} $$

を満たすとする．パラメータ $\lambda$ は何でもよい．本章では $\lambda$ を原子核座標に取るが，第11章11.2.3節では同じ定理を電子間相互作用の結合定数 $\lambda$ に適用して断熱接続公式 $E_{xc}=\int_0^1\dd\lambda\braket{\Psi_\lambda|\hat V_{ee}|\Psi_\lambda}-E_{\mathrm{H}}$ を導いている．以下の主張と証明は，そこで用いたもの(第11章の数学ノート)と同一である．

定理:Hellmann–Feynmanの定理

式\eqref{eq:16-hf-setup}のもとで，エネルギーのパラメータ微分は，ハミルトニアンの陽なパラメータ微分の期待値に等しい:

$$ \frac{\dd E}{\dd \lambda} = \left\langle \Psi(\lambda) \,\middle|\, \frac{\partial \hat{H}}{\partial \lambda} \,\middle|\, \Psi(\lambda) \right\rangle \equiv \int \Psi^*(\lambda)\, \frac{\partial \hat{H}}{\partial \lambda}\, \Psi(\lambda)\, \dd\tau. \label{eq:16-hf-theorem} $$

すなわち，波動関数が $\lambda$ とともに変化することによる寄与は，厳密に消える．ここで $\dd\tau$ は全電子座標(とスピン)についての積分を表す．

導出:波動関数微分の項が消えることの証明

エネルギーは期待値として書ける:

$$ E(\lambda) = \int \Psi^*(\lambda)\, \hat{H}(\lambda)\, \Psi(\lambda)\, \dd\tau. \label{eq:16-hf-e} $$

これを $\lambda$ で微分する．被積分関数は $\Psi^*$，$\hat{H}$，$\Psi$ の3つの因子の積だから，積の微分法則により3つの項が生じる(積分と微分の順序交換は，被積分関数が $\lambda$ について滑らかであれば許される):

\begin{align} \frac{\dd E}{\dd \lambda} &= \int \frac{\partial \Psi^*}{\partial \lambda}\, \hat{H}\, \Psi\, \dd\tau + \int \Psi^*\, \frac{\partial \hat{H}}{\partial \lambda}\, \Psi\, \dd\tau + \int \Psi^*\, \hat{H}\, \frac{\partial \Psi}{\partial \lambda}\, \dd\tau. \label{eq:16-hf-3terms} \end{align}

第1項と第3項(波動関数の変化に由来する項)を処理する．第1項では，$\Psi$ が固有状態であること $\hat{H}\Psi = E\Psi$ をそのまま代入する:

\begin{align} \int \frac{\partial \Psi^*}{\partial \lambda}\, \hat{H}\, \Psi\, \dd\tau = E \int \frac{\partial \Psi^*}{\partial \lambda}\, \Psi\, \dd\tau. \label{eq:16-hf-term1} \end{align}

第3項では，$\hat{H}$ がエルミートであることを使う．エルミート性とは任意の(境界で十分速く減衰する)関数 $f, g$ に対して $\int f^* \hat{H} g\, \dd\tau = \int (\hat{H} f)^* g\, \dd\tau$ が成り立つことだから，$f = \Psi$，$g = \partial\Psi/\partial\lambda$ と選んで

\begin{align} \int \Psi^*\, \hat{H}\, \frac{\partial \Psi}{\partial \lambda}\, \dd\tau = \int (\hat{H}\Psi)^*\, \frac{\partial \Psi}{\partial \lambda}\, \dd\tau = E \int \Psi^*\, \frac{\partial \Psi}{\partial \lambda}\, \dd\tau. \label{eq:16-hf-term3} \end{align}

最後の等号では $\hat{H}\Psi = E\Psi$ と，$E$ が実数であること($E^* = E$;エルミート演算子の固有値)を使った．式\eqref{eq:16-hf-term1}と\eqref{eq:16-hf-term3}を足すと，共通因子 $E$ でくくれて

\begin{align} E \int \left( \frac{\partial \Psi^*}{\partial \lambda}\, \Psi + \Psi^*\, \frac{\partial \Psi}{\partial \lambda} \right) \dd\tau = E\, \frac{\dd}{\dd \lambda} \int \Psi^* \Psi\, \dd\tau = E\, \frac{\dd}{\dd \lambda} (1) = 0. \label{eq:16-hf-cancel} \end{align}

1つ目の等号は積の微分法則を逆向きに使ったもの，2つ目の等号は規格化条件 $\int \Psi^*\Psi\, \dd\tau = 1$ がすべての $\lambda$ で成り立つ(定数関数である)ことによる．定数の微分はゼロである．したがって式\eqref{eq:16-hf-3terms}に残るのは第2項のみで，定理\eqref{eq:16-hf-theorem}が証明された．

得られたものを言葉で言い直すと:波動関数は $\lambda$ の変化に応じて変形するが，その変形はエネルギーの1階微分には寄与しない．これは変分原理の直接の帰結である．固有状態はエネルギー汎関数の停留点なのだから，波動関数を(規格化を保って)どの方向に微小変化させてもエネルギーは1次では変わらない．$\lambda$ の変化が引き起こす波動関数の変化も「その方向の一つ」にすぎない，というわけである．

16.2.2 力は静電気力である

次に $\lambda$ を原子核 $I$ の座標 $\bm{R}_I$ と読み替えて，力の具体形を導く．$N$ 電子系のハミルトニアン(第1章)は

$$ \hat{H} = \underbrace{-\frac{1}{2}\sum_{k=1}^{N} \nabla_k^2}_{\hat{T}_e} + \underbrace{\frac{1}{2}\sum_{k \neq l} \frac{1}{\abs{\rr_k - \rr_l}}}_{\hat{V}_{ee}} + \underbrace{\sum_{k=1}^{N} v_{\mathrm{ext}}(\rr_k)}_{\hat{V}_{\mathrm{ext}}} + \underbrace{\sum_{I \lt J} \frac{Z_I Z_J}{\abs{\bm{R}_I - \bm{R}_J}}}_{E_{nn}}, \qquad v_{\mathrm{ext}}(\rr) = -\sum_{I} \frac{Z_I}{\abs{\rr - \bm{R}_I}} \label{eq:16-full-h} $$

である．$\bm{R}_I$ に依存するのは $\hat{V}_{\mathrm{ext}}$ と $E_{nn}$ の2つだけであることに注意する．電子の運動エネルギー $\hat{T}_e$ と電子間反発 $\hat{V}_{ee}$ は電子座標だけの演算子だから，$\bm{R}_I$ 微分でゼロになる．したがって

$$ \frac{\partial \hat{H}}{\partial \bm{R}_I} = \sum_{k=1}^{N} \frac{\partial v_{\mathrm{ext}}(\rr_k)}{\partial \bm{R}_I} + \frac{\partial E_{nn}}{\partial \bm{R}_I}. \label{eq:16-dh-dr} $$

ここで必要になる微分を1つずつ実行しておく．$\bm{a}$ をベクトル変数として $1/\abs{\rr - \bm{a}}$ の $\bm{a}$ 微分を計算する．成分で書くと $\abs{\rr - \bm{a}} = [(x-a_x)^2 + (y-a_y)^2 + (z-a_z)^2]^{1/2}$ だから，$a_x$ 成分の微分は合成関数の微分法則により

\begin{align} \frac{\partial}{\partial a_x} \frac{1}{\abs{\rr - \bm{a}}} &= -\frac{1}{2} \left[ (x-a_x)^2 + (y-a_y)^2 + (z-a_z)^2 \right]^{-3/2} \times \left\{ 2 (x - a_x) \cdot (-1) \right\} = \frac{x - a_x}{\abs{\rr - \bm{a}}^3}. \label{eq:16-coulomb-grad-x} \end{align}

1つ目の等号は $u^{-1/2}$ 型の外側の微分($-\tfrac{1}{2}u^{-3/2}$)に内側 $u = \abs{\rr-\bm{a}}^2$ の $a_x$ 微分を掛けたもの，内側の微分で $(x-a_x)^2$ から因子 $-2(x-a_x)$ が出る．$y, z$ 成分も同様だから，ベクトルにまとめて

$$ \frac{\partial}{\partial \bm{a}} \frac{1}{\abs{\rr - \bm{a}}} = \frac{\rr - \bm{a}}{\abs{\rr - \bm{a}}^3}. \label{eq:16-coulomb-grad} $$

もう1つ，多体波動関数の期待値を密度で書き直す関係を思い出す(第4章)．1体演算子 $\sum_k f(\rr_k)$ の期待値は，$\abs{\Psi}^2$ が電子の入れ替えについて対称であること(反対称波動関数の絶対値の2乗)から，$N$ 個の項がすべて等しく，

\begin{align} \left\langle \Psi \,\middle|\, \sum_{k=1}^{N} f(\rr_k) \,\middle|\, \Psi \right\rangle = N \int f(\rr_1) \abs{\Psi(\rr_1, \dots, \rr_N)}^2 \dd\rr_1 \cdots \dd\rr_N = \int f(\rr)\, n(\rr)\, \dd^3 r \label{eq:16-onebody} \end{align}

となる．最後の等号では電子密度の定義 $n(\rr) = N \int \abs{\Psi(\rr, \rr_2, \dots, \rr_N)}^2 \dd\rr_2 \cdots \dd\rr_N$ を使った．

式\eqref{eq:16-hf-theorem}に式\eqref{eq:16-dh-dr}を代入し，式\eqref{eq:16-coulomb-grad}と\eqref{eq:16-onebody}を使って各項を評価する．まず外部ポテンシャル項:$v_{\mathrm{ext}}$ のうち $\bm{R}_I$ を含むのは $-Z_I/\abs{\rr - \bm{R}_I}$ だけだから

\begin{align} \left\langle \Psi \,\middle|\, \sum_k \frac{\partial v_{\mathrm{ext}}(\rr_k)}{\partial \bm{R}_I} \,\middle|\, \Psi \right\rangle = \int n(\rr)\, \frac{\partial}{\partial \bm{R}_I}\left( -\frac{Z_I}{\abs{\rr - \bm{R}_I}} \right) \dd^3 r = -Z_I \int n(\rr)\, \frac{\rr - \bm{R}_I}{\abs{\rr - \bm{R}_I}^3}\, \dd^3 r. \label{eq:16-force-elec-term} \end{align}

次に核間反発項:和のうち $\bm{R}_I$ を含む対だけが残り，式\eqref{eq:16-coulomb-grad}(今度は $\rr \to \bm{R}_J$，$\bm{a} \to \bm{R}_I$)により

\begin{align} \frac{\partial E_{nn}}{\partial \bm{R}_I} = \sum_{J \neq I} Z_I Z_J\, \frac{\partial}{\partial \bm{R}_I} \frac{1}{\abs{\bm{R}_I - \bm{R}_J}} = -\sum_{J \neq I} Z_I Z_J\, \frac{\bm{R}_I - \bm{R}_J}{\abs{\bm{R}_I - \bm{R}_J}^3}. \label{eq:16-force-nn-term} \end{align}

ここで符号に注意:式\eqref{eq:16-coulomb-grad}は「分母の第2引数」での微分だったから，$\abs{\bm{R}_I - \bm{R}_J}^{-1}$ を $\bm{R}_I$(第1引数)で微分すると符号が反転する．力 $\bm{F}_I = -\partial E/\partial \bm{R}_I$ は，以上の2項の符号を反転して

$$ \bm{F}_I = Z_I \int n(\rr)\, \frac{\rr - \bm{R}_I}{\abs{\rr - \bm{R}_I}^3}\, \dd^3 r \;-\; Z_I \sum_{J \neq I} Z_J\, \frac{\bm{R}_J - \bm{R}_I}{\abs{\bm{R}_J - \bm{R}_I}^3}. \label{eq:16-hf-force} $$

これがHellmann–Feynman力である．何が得られたかをよく見てほしい．第1項は「電荷密度 $-n(\rr)$ の電子雲が，電荷 $+Z_I$ の点電荷を引っぱる古典的な静電引力」，第2項は「他の原子核からの古典的な静電反発」である．すなわち:

物理的意味:静電定理

基底状態の電子密度 $n(\rr)$ さえ正確に分かってしまえば，原子核に働く力は純粋に古典静電気学だけで計算できる．量子力学(運動エネルギー・交換・相関のすべて)の役割は密度の形を決めることに尽き，力の表式そのものには現れない．この見方はFeynmanが1939年の論文で強調したもので，静電定理とも呼ばれる．化学結合を「原子核の間に密度がたまり，その負電荷が両側の原子核を内側へ引く」と描く直観的な結合描像は，この定理によって定量的に正当化される．

16.2.3 定理が成り立つ条件

導出を振り返ると，使った仮定は次の2つだけである:(i) $\Psi$ が $\hat{H}$ の厳密な固有状態であること(式\eqref{eq:16-hf-term1}, \eqref{eq:16-hf-term3}で使用)，(ii) 規格化が $\lambda$ によらず保たれること(式\eqref{eq:16-hf-cancel}で使用)．実際の数値計算では (i) が破れる——有限個の基底関数の張る空間の中でしか波動関数を最適化できないからである．ただし，厳密固有状態でなくても，使える変分パラメータのすべてについてエネルギーが停留していれば，パラメータ微分に由来する項は同じ理屈で消える．問題は，変分空間そのものが $\lambda$(原子核の位置)に依存する場合である．これが次項の主題，Pulay力である．

数学ノート:重なり行列と一般化固有値問題

互いに直交しない基底関数 $\{\phi_\mu\}_{\mu=1}^{n}$ で波動関数を $\psi = \sum_\mu c_\mu \phi_\mu$ と展開するとき，内積は重なり行列 $S_{\mu\nu} = \int \phi_\mu^* \phi_\nu\, \dd^3 r$ を使って $\braket{\psi|\psi} = \sum_{\mu\nu} c_\mu^* S_{\mu\nu} c_\nu = \bm{c}^\dagger S \bm{c}$ と書ける(平面波なら $S_{\mu\nu} = \delta_{\mu\nu}$ だが，原子中心の局在基底では一般に $S \neq I$)．ハミルトニアン行列を $H_{\mu\nu} = \int \phi_\mu^* \hat{H} \phi_\nu\, \dd^3 r$ とすると，拘束条件 $\bm{c}^\dagger S \bm{c} = 1$ のもとで $\bm{c}^\dagger H \bm{c}$ を最小化する条件は，Lagrange未定乗数 $\varepsilon$ を導入して $\partial/\partial c_\mu^* \left[ \bm{c}^\dagger H \bm{c} - \varepsilon (\bm{c}^\dagger S \bm{c} - 1) \right] = 0$，すなわち

$$ \sum_\nu H_{\mu\nu} c_\nu = \varepsilon \sum_\nu S_{\mu\nu} c_\nu \qquad \Longleftrightarrow \qquad H \bm{c} = \varepsilon S \bm{c} \label{eq:16-gep} $$

となる．これを一般化固有値問題と呼ぶ．$c_\mu$ と $c_\mu^*$ を独立変数として扱う操作は，実部と虚部を独立に動かすことと同値である(第4章の変分計算と同じ手続き)．未定乗数 $\varepsilon$ は軌道エネルギーの意味を持つ．

16.2.4 Pulay力:基底関数が原子と一緒に動く場合の補正

原子中心の基底関数(第14章の局在基底(LCPAO)や，量子化学のGauss基底)を使うと，基底関数そのものが $\phi_\mu(\rr - \bm{R}_{I_\mu})$ という形で原子核座標に依存する．原子が動けば「ものさし」も一緒に動くのである．このとき，Hellmann–Feynman力だけでは全エネルギーの微分にならない．これから，その不足分——Pulay力(Pulay, 1969)——を明示的に導出する．

見通しをよくするため，まず「与えられたエルミート演算子 $\hat{H}$ の期待値を有限基底の中で変分的に最小化した」問題を考える．1つの規格化軌道 $\psi = \sum_\mu c_\mu \phi_\mu$ のエネルギー

$$ \varepsilon(R) = \sum_{\mu\nu} c_\mu^* H_{\mu\nu}(R)\, c_\nu, \qquad \sum_{\mu\nu} c_\mu^* S_{\mu\nu}(R)\, c_\nu = 1 \label{eq:16-pulay-setup} $$

を考える．$R$ はある原子核座標の1成分で，行列要素 $H_{\mu\nu}, S_{\mu\nu}$ は2つの経路で $R$ に依存する:$\hat{H}$ の中のポテンシャルを通じて(陽な依存)，そして基底関数 $\phi_\mu(\rr; R)$ を通じて(基底の依存)．係数 $\bm{c}$ は各 $R$ で変分条件\eqref{eq:16-gep}を満たすように決められているとする．

導出:有限基底でのエネルギー微分とPulay項

Step 1:全微分を書き下す．式\eqref{eq:16-pulay-setup}を $R$ で微分する．依存性は $c_\mu^*$，$c_\nu$，$H_{\mu\nu}$ の3箇所にあるから，積の微分法則で

\begin{align} \frac{\dd \varepsilon}{\dd R} = \sum_{\mu\nu} \frac{\dd c_\mu^*}{\dd R} H_{\mu\nu} c_\nu + \sum_{\mu\nu} c_\mu^* H_{\mu\nu} \frac{\dd c_\nu}{\dd R} + \sum_{\mu\nu} c_\mu^* \frac{\dd H_{\mu\nu}}{\dd R} c_\nu. \label{eq:16-pulay-total} \end{align}

Step 2:係数微分の項を変分条件で書き換える．第1項に一般化固有値方程式 $\sum_\nu H_{\mu\nu} c_\nu = \varepsilon \sum_\nu S_{\mu\nu} c_\nu$(式\eqref{eq:16-gep})を代入し，第2項にはそのエルミート共役 $\sum_\mu c_\mu^* H_{\mu\nu} = \varepsilon \sum_\mu c_\mu^* S_{\mu\nu}$($H, S$ はエルミート，$\varepsilon$ は実数)を代入する:

\begin{align} \text{(第1項)} + \text{(第2項)} = \varepsilon \left[ \sum_{\mu\nu} \frac{\dd c_\mu^*}{\dd R} S_{\mu\nu} c_\nu + \sum_{\mu\nu} c_\mu^* S_{\mu\nu} \frac{\dd c_\nu}{\dd R} \right]. \label{eq:16-pulay-step2} \end{align}

Step 3:規格化条件の微分を使う．拘束 $\bm{c}^\dagger S \bm{c} = 1$ はすべての $R$ で成り立つから，その全微分はゼロである．積の微分法則で3項に展開すると

\begin{align} 0 = \frac{\dd}{\dd R}\left( \bm{c}^\dagger S \bm{c} \right) = \sum_{\mu\nu} \frac{\dd c_\mu^*}{\dd R} S_{\mu\nu} c_\nu + \sum_{\mu\nu} c_\mu^* S_{\mu\nu} \frac{\dd c_\nu}{\dd R} + \sum_{\mu\nu} c_\mu^* \frac{\partial S_{\mu\nu}}{\partial R} c_\nu. \label{eq:16-pulay-norm} \end{align}

すなわち式\eqref{eq:16-pulay-step2}の角括弧は $-\sum_{\mu\nu} c_\mu^* (\partial S_{\mu\nu}/\partial R) c_\nu$ に等しい．よって

\begin{align} \text{(第1項)} + \text{(第2項)} = -\varepsilon \sum_{\mu\nu} c_\mu^* \frac{\partial S_{\mu\nu}}{\partial R} c_\nu. \label{eq:16-pulay-step3} \end{align}

係数の微分 $\dd c/\dd R$ そのものを計算する必要が消えたことに注目してほしい．これがHellmann–Feynman定理の証明(式\eqref{eq:16-hf-cancel})の有限基底版である．ただし今回はゼロにはならず，重なり行列の微分が残る．

Step 4:行列要素の微分を展開する．$H_{\mu\nu} = \int \phi_\mu^* \hat{H} \phi_\nu\, \dd^3 r$ は $\phi_\mu^*$，$\hat{H}$，$\phi_\nu$ の3因子を含むから，

\begin{align} \frac{\dd H_{\mu\nu}}{\dd R} = \int \frac{\partial \phi_\mu^*}{\partial R} \hat{H} \phi_\nu\, \dd^3 r + \int \phi_\mu^* \frac{\partial \hat{H}}{\partial R} \phi_\nu\, \dd^3 r + \int \phi_\mu^* \hat{H} \frac{\partial \phi_\nu}{\partial R}\, \dd^3 r, \label{eq:16-pulay-dh} \\[2pt] \frac{\partial S_{\mu\nu}}{\partial R} = \int \frac{\partial \phi_\mu^*}{\partial R} \phi_\nu\, \dd^3 r + \int \phi_\mu^* \frac{\partial \phi_\nu}{\partial R}\, \dd^3 r. \label{eq:16-pulay-ds} \end{align}

Step 5:すべてを集めて整理する．式\eqref{eq:16-pulay-total}に式\eqref{eq:16-pulay-step3}, \eqref{eq:16-pulay-dh}, \eqref{eq:16-pulay-ds}を代入し，$\varepsilon \partial S/\partial R$ の項を $\dd H/\dd R$ の基底微分項と組にする:

\begin{align} \frac{\dd \varepsilon}{\dd R} &= \sum_{\mu\nu} c_\mu^* \left[ \int \phi_\mu^* \frac{\partial \hat{H}}{\partial R} \phi_\nu\, \dd^3 r \right] c_\nu \nonumber\\ &\quad + \sum_{\mu\nu} c_\mu^* c_\nu \int \frac{\partial \phi_\mu^*}{\partial R} \left( \hat{H} - \varepsilon \right) \phi_\nu\, \dd^3 r + \sum_{\mu\nu} c_\mu^* c_\nu \int \phi_\mu^* \left( \hat{H} - \varepsilon \right) \frac{\partial \phi_\nu}{\partial R}\, \dd^3 r. \label{eq:16-pulay-combined} \end{align}

ここで $\sum_\nu c_\nu \phi_\nu = \psi$，$\sum_\mu c_\mu^* \phi_\mu^* = \psi^*$ とまとめ，さらに第3項が第2項の複素共役であること($\hat{H} - \varepsilon$ のエルミート性から $\int \phi_\mu^* (\hat{H}-\varepsilon) \partial_R \phi_\nu\, \dd^3 r = [\int \partial_R \phi_\nu^* (\hat{H}-\varepsilon) \phi_\mu\, \dd^3 r]^*$)を使うと，最終形が得られる:

\begin{align} \frac{\dd \varepsilon}{\dd R} = \underbrace{\left\langle \psi \,\middle|\, \frac{\partial \hat{H}}{\partial R} \,\middle|\, \psi \right\rangle}_{\text{Hellmann–Feynman項}} + \underbrace{2\, \mathrm{Re} \sum_{\mu} c_\mu^* \left\langle \frac{\partial \phi_\mu}{\partial R} \,\middle|\, \hat{H} - \varepsilon \,\middle|\, \psi \right\rangle}_{\text{Pulay項}}. \label{eq:16-pulay-final} \end{align}

得られた式\eqref{eq:16-pulay-final}を吟味しよう．Pulay項が消えるのは次のいずれかの場合である．

基底が $R$ に依存しない場合:$\partial \phi_\mu / \partial R = 0$．平面波基底がまさにこれで，平面波 $e^{i\bm{G}\cdot\rr}$ は原子の位置を知らない．平面波基底ではPulay力はゼロであり，Hellmann–Feynman力(と擬ポテンシャル項の陽微分)だけで正確な力が得られる．これは平面波法の大きな実務的利点である．
基底が完全系の場合:$\psi$ が厳密固有状態になり $(\hat{H} - \varepsilon)\ket{\psi} = 0$，よって任意の左ベクトルとの内積が消える．

有限の原子中心基底ではどちらも成り立たない．変分条件\eqref{eq:16-gep}が保証するのは，残差ベクトル $(\hat{H} - \varepsilon)|\psi\rangle$ が基底の張る空間に直交することまでである．ところが $\partial \phi_\mu/\partial R$ は一般に基底の張る空間の外に突き出ているから，内積\eqref{eq:16-pulay-final}は消えない．「ものさしの動き」がエネルギー微分に混入する，これがPulay力の正体である．なお $|\partial \phi_\mu/\partial R\rangle$ がたまたま基底空間内に収まる特殊な場合(完全系に限らない)にもPulay項は消えるが，実用上は「原子中心基底なら必ずPulay補正を入れる」と覚えてよい．

多電子系への一般化は形式的には単純で，占有軌道 $\psi_i$($i = 1, \dots, N_{\mathrm{occ}}$，占有数 $f_i$)について同じ計算を繰り返せばよい．Kohn–Sham全エネルギーはバンドエネルギーの単純和ではない(二重数え補正がある;第10章10.4節)ため，厳密な表式では密度の変分停留性も併用して整理することになるが，結果の構造は同じであり，力は

$$ \bm{F}_I = \bm{F}_I^{\mathrm{HF}} + \bm{F}_I^{\mathrm{Pulay}}, \qquad \bm{F}_I^{\mathrm{Pulay}} = -2\, \mathrm{Re} \sum_{i}^{\mathrm{occ}} f_i \sum_{\mu \in I} \sum_{\nu} c_{\mu i}^* \left\langle \frac{\partial \phi_\mu}{\partial \bm{R}_I} \,\middle|\, \hat{H}_{\mathrm{KS}} - \varepsilon_i \,\middle|\, \phi_\nu \right\rangle c_{\nu i} \label{eq:16-pulay-force} $$

という「Hellmann–Feynman力+Pulay補正」の形にまとまる．ここで $\mu \in I$ は原子 $I$ を中心とする基底関数だけが和に入ることを表す($\partial \phi_\mu/\partial \bm{R}_I$ は $\phi_\mu$ が原子 $I$ 上にあるときだけ非零)．局在基底を用いる第一原理計算コードでは，全エネルギーを構成する各項(運動エネルギー，局所・非局所擬ポテンシャル，Hartree，交換相関，イオン間反発)の $\bm{R}_I$ 微分をこの形で解析的に実装しており，数値グリッド上の積分を使っていても，離散化された全エネルギーと厳密に整合した力が得られるように設計されている．

例:Pulay補正を忘れるとどうなるか

局在基底の計算でHellmann–Feynman項だけを使って構造最適化を行うと，力とエネルギーが整合しないため，「力はゼロなのにエネルギーはまだ下がる」「エネルギー極小と力ゼロの構造が食い違う」という病的な振る舞いが起きる．典型的な原子中心基底では，Pulay項の大きさはHellmann–Feynman項と同程度(結合長スケールで $10^{-2}$〜$10^{-1}$ Hartree/Bohr)にもなり，無視できる小さな補正では全くない．分子動力学では力の系統誤差が全エネルギーのドリフトとして現れるため，Pulay項の正確な実装は死活問題である．

16.3 応力テンソル:力の「格子版」

分子であれば，力\eqref{eq:16-hf-force}(と\eqref{eq:16-pulay-force})だけで構造を決められる．しかし周期系(結晶)では，原子の内部座標に加えて単位胞の形と大きさ——格子ベクトル $\bm{a}_1, \bm{a}_2, \bm{a}_3$ の9成分——も構造の自由度である．格子に対する「力」に相当する量が応力テンソル(stress tensor)であり，これがあって初めて格子定数の最適化や弾性定数の計算ができる．本節では応力を「一様ひずみに対するエネルギーの1階微分」として定義し，運動エネルギーとCoulombエネルギーのひずみ微分を実際に計算して，第2章のビリアル定理と正確につながることを確認する．

16.3.1 一様ひずみと応力の定義

定義:ひずみテンソルと応力テンソル

系のすべての位置ベクトル(格子ベクトル・原子座標)を，対称な $3\times 3$ 行列 $\varepsilon_{\alpha\beta}$(ひずみテンソル)によって一斉に線形変換する:

$$ \rr \;\longrightarrow\; \rr' = (\bm{1} + \bm{\varepsilon})\, \rr, \qquad r'_\alpha = r_\alpha + \sum_\beta \varepsilon_{\alpha\beta}\, r_\beta. \label{eq:16-strain-def} $$

対角成分 $\varepsilon_{xx}$ 等は $x$ 方向の伸び(正)・縮み(負)を，非対角成分 $\varepsilon_{xy}$ 等はせん断変形を表す．このひずみを受けた系の全エネルギーを $E(\bm{\varepsilon})$ と書くとき，応力テンソルを

$$ \sigma_{\alpha\beta} = \frac{1}{\Omega} \left. \frac{\partial E}{\partial \varepsilon_{\alpha\beta}} \right|_{\bm{\varepsilon} = 0} \label{eq:16-stress-def} $$

で定義する．$\Omega$ は単位胞の体積である(応力は「単位体積あたりのエネルギー微分」であり，圧力と同じ次元を持つ)．符号の規約はコードや教科書により逆のものもあるので，実際の計算結果を読むときは必ず確認すること．

図16.2 一様ひずみによる単位胞の変形．ひずみテンソル $\varepsilon_{\alpha\beta}$ は，格子ベクトルと胞内の原子位置をすべて同じ線形変換で動かす．エネルギーのひずみ微分が応力テンソルである．

まず，定義\eqref{eq:16-stress-def}から圧力との関係を確かめておく．等方的な微小ひずみ $\varepsilon_{\alpha\beta} = \eta\, \delta_{\alpha\beta}$ を掛けると，体積は $\Omega \to \Omega \det(\bm{1} + \bm{\varepsilon})$ と変わる．

数学ノート:$\det(\bm{1} + \bm{\varepsilon})$ の1次展開

$3 \times 3$ 行列式を定義どおり展開すると，非対角要素 $\varepsilon_{\alpha\beta}$($\alpha \neq \beta$)は必ず2個以上の積でしか現れない(行列式の各項は行・列から1個ずつ要素を取るため，非対角要素を1個使うと別の非対角要素も使わざるを得ない)．よって $\bm{\varepsilon}$ の1次までは対角積のみが効き，

$$ \det(\bm{1} + \bm{\varepsilon}) = (1 + \varepsilon_{xx})(1 + \varepsilon_{yy})(1 + \varepsilon_{zz}) + O(\varepsilon^2) = 1 + \mathrm{Tr}\, \bm{\varepsilon} + O(\varepsilon^2). \label{eq:16-det-expand} $$

したがって $\partial \det(\bm{1}+\bm{\varepsilon}) / \partial \varepsilon_{\alpha\beta} |_{\varepsilon=0} = \delta_{\alpha\beta}$，体積については $\partial \Omega(\bm{\varepsilon}) / \partial \varepsilon_{\alpha\beta}|_{\varepsilon=0} = \Omega\, \delta_{\alpha\beta}$ である．同様に逆行列の展開 $(\bm{1} + \bm{\varepsilon})^{-1} = \bm{1} - \bm{\varepsilon} + \bm{\varepsilon}^2 - \cdots$ も，両辺に $(\bm{1}+\bm{\varepsilon})$ を掛ければ確認できる(幾何級数の行列版)．以下では両方とも1次までで使う．

等方ひずみでは $E$ は体積だけの関数として変わるから，合成関数の微分により

\begin{align} \sum_\alpha \frac{\partial E}{\partial \varepsilon_{\alpha\alpha}} = \sum_\alpha \frac{\dd E}{\dd \Omega} \frac{\partial \Omega}{\partial \varepsilon_{\alpha\alpha}} = \frac{\dd E}{\dd \Omega} \cdot 3\Omega \qquad \Longrightarrow \qquad P = -\frac{\dd E}{\dd \Omega} = -\frac{1}{3} \sum_\alpha \sigma_{\alpha\alpha}. \label{eq:16-pressure} \end{align}

すなわち圧力は応力テンソルのトレースの $-1/3$ 倍である．応力の対角和がゼロになる格子定数が，外圧ゼロでの平衡格子定数にほかならない．

16.3.2 波動関数のスケーリングと変分原理

応力を計算するには $E(\bm{\varepsilon})$ の微分が要るが，ここで一つ概念的な問題がある．ひずみを掛けると単位胞(波動関数の定義域と境界条件)自体が変わってしまうので，「同じ波動関数のままハミルトニアンだけ微分する」というHellmann–Feynmanの定理の使い方が直接はできない．この問題は，波動関数をスケーリング変換して元の領域に引き戻すことで解決する．各軌道(あるいは多体波動関数の各電子座標)を

$$ \psi_{\bm{\varepsilon}}(\rr) = \left[ \det(\bm{1} + \bm{\varepsilon}) \right]^{-1/2} \psi\!\left( (\bm{1} + \bm{\varepsilon})^{-1} \rr \right) \label{eq:16-scaled-wf} $$

と変換する．前係数は規格化を保つためのものである．実際，$\bm{u} = (\bm{1}+\bm{\varepsilon})^{-1} \rr$ と置換すると体積要素は $\dd^3 r = \det(\bm{1}+\bm{\varepsilon})\, \dd^3 u$(線形変換のJacobian)だから，

\begin{align} \int \abs{\psi_{\bm{\varepsilon}}(\rr)}^2 \dd^3 r = \frac{1}{\det(\bm{1}+\bm{\varepsilon})} \int \abs{\psi(\bm{u})}^2 \det(\bm{1}+\bm{\varepsilon})\, \dd^3 u = \int \abs{\psi(\bm{u})}^2 \dd^3 u = 1 \label{eq:16-scaled-norm} \end{align}

となり，任意の $\bm{\varepsilon}$ で規格化(および軌道どうしの直交性)が保たれる．ひずんだ格子の周期境界条件も，$\psi$ が元の格子の境界条件を満たしていれば $\psi_{\bm{\varepsilon}}$ が自動的に満たす．

ここで変分原理が効いてくる．ひずんだ系の真の基底状態は $\psi_{\bm{\varepsilon}}$ そのものではない．しかし，$\psi_{\bm{\varepsilon}}$ は「$\bm{\varepsilon} = 0$ で真の基底状態を通る，規格直交条件を満たす試行関数の族」である．エネルギーの $\bm{\varepsilon}$ 微分のうち「波動関数が族からずれて緩和する」寄与は，基底状態がエネルギーの停留点であることから1次では消える(16.2節の議論と同一の論理)．したがって，応力の1階微分は，スケーリングした波動関数を代入した期待値 $E(\bm{\varepsilon}) = \langle \Psi_{\bm{\varepsilon}} | \hat{H}(\bm{\varepsilon}) | \Psi_{\bm{\varepsilon}} \rangle$ を $\bm{\varepsilon}$ で微分すれば正しく得られる．以下，これを運動エネルギーとCoulombエネルギーについて実行する．

16.3.3 運動エネルギーのひずみ微分

これから，1つの軌道の運動エネルギー $T(\bm{\varepsilon}) = \frac{1}{2} \int \abs{\nabla \psi_{\bm{\varepsilon}}}^2 \dd^3 r$ のひずみ微分を計算する(占有軌道すべての和は最後に取ればよい)．

導出:運動エネルギー項の応力

Step 1:置換積分で $\bm{\varepsilon}$ 依存性を露出させる．$\bm{A} \equiv (\bm{1}+\bm{\varepsilon})^{-1}$ と略記する．$\bm{u} = \bm{A}\rr$ と置換すると，$u_\beta = \sum_\gamma A_{\beta\gamma} r_\gamma$ だから $\partial u_\beta / \partial r_\alpha = A_{\beta\alpha}$ であり，合成関数の微分法則により

\begin{align} \frac{\partial}{\partial r_\alpha} \psi(\bm{A}\rr) = \sum_\beta \frac{\partial u_\beta}{\partial r_\alpha} \frac{\partial \psi}{\partial u_\beta} = \sum_\beta A_{\beta\alpha}\, \partial_\beta \psi(\bm{u}), \qquad \partial_\beta \equiv \frac{\partial}{\partial u_\beta}. \label{eq:16-kin-chain} \end{align}

これを $T(\bm{\varepsilon})$ に代入する．規格化因子の2乗 $\det^{-1}$ と体積要素の $\det$ が打ち消し合い($\dd^3 r = \det(\bm{1}+\bm{\varepsilon})\, \dd^3 u$)，

\begin{align} T(\bm{\varepsilon}) = \frac{1}{2} \sum_\alpha \int \left[ \sum_\beta A_{\beta\alpha}\, \partial_\beta \psi^*(\bm{u}) \right] \left[ \sum_\gamma A_{\gamma\alpha}\, \partial_\gamma \psi(\bm{u}) \right] \dd^3 u = \frac{1}{2} \sum_{\alpha\beta\gamma} A_{\beta\alpha} A_{\gamma\alpha} \int \partial_\beta \psi^*\, \partial_\gamma \psi\, \dd^3 u. \label{eq:16-kin-sub} \end{align}

ひずみ依存性が，積分の外の行列因子 $A_{\beta\alpha} A_{\gamma\alpha}$ に完全に押し出されたことに注目してほしい．積分自体は $\bm{\varepsilon}$ を含まない．

Step 2:1次まで展開する．$\bm{A} = \bm{1} - \bm{\varepsilon} + O(\varepsilon^2)$ を代入し，$\bm{\varepsilon}$ の1次までを拾う:

\begin{align} A_{\beta\alpha} A_{\gamma\alpha} = (\delta_{\beta\alpha} - \varepsilon_{\beta\alpha})(\delta_{\gamma\alpha} - \varepsilon_{\gamma\alpha}) + O(\varepsilon^2) = \delta_{\beta\alpha}\delta_{\gamma\alpha} - \varepsilon_{\beta\alpha}\, \delta_{\gamma\alpha} - \varepsilon_{\gamma\alpha}\, \delta_{\beta\alpha} + O(\varepsilon^2). \label{eq:16-kin-expand} \end{align}

第1項はひずみのない運動エネルギー $T(0)$ を再現する．第2・第3項を式\eqref{eq:16-kin-sub}に戻すと

\begin{align} T(\bm{\varepsilon}) = T(0) - \frac{1}{2} \sum_{\alpha\beta} \varepsilon_{\beta\alpha} \int \left( \partial_\beta \psi^*\, \partial_\alpha \psi + \partial_\alpha \psi^*\, \partial_\beta \psi \right) \dd^3 u + O(\varepsilon^2). \label{eq:16-kin-firstorder} \end{align}

ここで第2項ではクロネッカーのデルタで $\gamma$ の和を，第3項では $\beta$ の和を実行し，残った添字を付け替えた．

Step 3:微分を読み取る．$\varepsilon_{\alpha\beta}$ の係数がそのまま1階微分である:

\begin{align} \left. \frac{\partial T}{\partial \varepsilon_{\alpha\beta}} \right|_{\bm{\varepsilon}=0} = -\frac{1}{2} \int \left( \partial_\alpha \psi^*\, \partial_\beta \psi + \partial_\beta \psi^*\, \partial_\alpha \psi \right) \dd^3 r = -\,\mathrm{Re} \int \partial_\alpha \psi^*\, \partial_\beta \psi\, \dd^3 r. \label{eq:16-kin-stress} \end{align}

占有軌道について和を取れば，運動エネルギー項の応力への寄与 $\Omega\, \sigma^{\mathrm{kin}}_{\alpha\beta} = -\sum_i f_i\, \mathrm{Re} \int \partial_\alpha \psi_i^*\, \partial_\beta \psi_i\, \dd^3 r$ が得られる．

得られた式\eqref{eq:16-kin-stress}のトレースを取ると，$\sum_\alpha \partial_\alpha \psi^* \partial_\alpha \psi = \abs{\nabla\psi}^2$ だから

$$ \sum_\alpha \frac{\partial T}{\partial \varepsilon_{\alpha\alpha}} = -\int \abs{\nabla \psi}^2 \dd^3 r = -2\, T. \label{eq:16-kin-trace} $$

これは「長さを $\lambda$ 倍に伸ばすと運動エネルギーは $\lambda^{-2}$ 倍になる」というスケーリング則($T \propto \lambda^{-2}$ の対数微分が $-2$)の精密化である．符号を確かめておこう．式\eqref{eq:16-pressure}に代入すると $P^{\mathrm{kin}} = -\frac{1}{3\Omega}(-2T) = +2T/3\Omega \gt 0$，つまり運動エネルギーは常に正の圧力(系を押し広げる向き)に寄与する．量子力学的な「閉じ込めのコスト」が外向きに押す圧力になる——白色矮星を支える縮退圧と同じ物理である(第3章のフェルミ球の運動エネルギー密度を思い出そう)．

16.3.4 Coulombエネルギーのひずみ微分とビリアル定理

次に，Coulomb相互作用(電子間・電子核間・核間のすべて)のひずみ微分を計算する．どの項も「2点間距離の逆数」の和・積分だから，核となる計算は $1/\abs{(\bm{1}+\bm{\varepsilon})\bm{u}}$ の $\bm{\varepsilon}$ 微分である．

導出:Coulomb項の応力

Step 1:距離の2乗を展開する．電子座標をスケーリングした波動関数\eqref{eq:16-scaled-wf}で期待値を取り，すべての電子座標を $\rr_k = (\bm{1}+\bm{\varepsilon})\bm{u}_k$ と置換すると(核座標ももともと $(\bm{1}+\bm{\varepsilon})$ 倍されている)，Jacobianと規格化因子は運動エネルギーのときと同様に打ち消し合い，Coulomb核だけが $\bm{\varepsilon}$ 依存性を持つ:任意の2粒子の相対ベクトル $\bm{u}$ に対して $\abs{\bm{u}} \to \abs{(\bm{1}+\bm{\varepsilon})\bm{u}}$．その2乗を成分で書くと

\begin{align} \abs{(\bm{1}+\bm{\varepsilon})\bm{u}}^2 = \sum_\alpha \Bigl( u_\alpha + \sum_\beta \varepsilon_{\alpha\beta} u_\beta \Bigr)^2 = \abs{\bm{u}}^2 + 2 \sum_{\alpha\beta} \varepsilon_{\alpha\beta}\, u_\alpha u_\beta + O(\varepsilon^2). \label{eq:16-cou-dist2} \end{align}

2乗を展開して $\bm{\varepsilon}$ の1次までを残した．

Step 2:逆数を展開する．$x \equiv \sum_{\alpha\beta} \varepsilon_{\alpha\beta} u_\alpha u_\beta / \abs{\bm{u}}^2$ と置くと，式\eqref{eq:16-cou-dist2}より $\abs{(\bm{1}+\bm{\varepsilon})\bm{u}} = \abs{\bm{u}} (1 + 2x)^{1/2}$ であり，二項展開 $(1+2x)^{-1/2} = 1 - x + O(x^2)$ を使って

\begin{align} \frac{1}{\abs{(\bm{1}+\bm{\varepsilon})\bm{u}}} = \frac{1}{\abs{\bm{u}}} \left( 1 - \sum_{\alpha\beta} \varepsilon_{\alpha\beta} \frac{u_\alpha u_\beta}{\abs{\bm{u}}^2} \right) + O(\varepsilon^2) \qquad \Longrightarrow \qquad \left. \frac{\partial}{\partial \varepsilon_{\alpha\beta}} \frac{1}{\abs{(\bm{1}+\bm{\varepsilon})\bm{u}}} \right|_{\bm{\varepsilon}=0} = -\frac{u_\alpha u_\beta}{\abs{\bm{u}}^3}. \label{eq:16-cou-kernel} \end{align}

Step 3:全Coulombエネルギーに適用する．全Coulombエネルギー $E_{\mathrm{C}}$(電子間反発 $+$ 電子核引力 $+$ 核間反発)は，電荷対の相対ベクトル $\bm{u}$ を持つ $1/\abs{\bm{u}}$ 型の項の(符号付き)和・期待値だから，微分\eqref{eq:16-cou-kernel}を項ごとに適用すればよい．形式的にまとめて書けば

\begin{align} \frac{\partial E_{\mathrm{C}}}{\partial \varepsilon_{\alpha\beta}} = -\left\langle \sum_{\text{pairs}\, (i,j)} q_i q_j\, \frac{u_\alpha u_\beta}{\abs{\bm{u}}^3} \right\rangle, \qquad \bm{u} = \bm{x}_i - \bm{x}_j \label{eq:16-cou-stress} \end{align}

である($q_i$ は電子なら $-1$，核なら $+Z_I$;$\langle \cdot \rangle$ は電子座標についての期待値)．

Step 4:トレースを取る．$\sum_\alpha u_\alpha u_\alpha / \abs{\bm{u}}^3 = \abs{\bm{u}}^2/\abs{\bm{u}}^3 = 1/\abs{\bm{u}}$ だから，

\begin{align} \sum_\alpha \frac{\partial E_{\mathrm{C}}}{\partial \varepsilon_{\alpha\alpha}} = -\left\langle \sum_{\text{pairs}} \frac{q_i q_j}{\abs{\bm{u}}} \right\rangle = -E_{\mathrm{C}}. \label{eq:16-cou-trace} \end{align}

これは「Coulombエネルギーは長さの $-1$ 乗でスケールする」($V \propto \lambda^{-1}$)ことの表現である．

以上を合わせると，全エネルギー $E = T + E_{\mathrm{C}}$(占有軌道和を含めた全運動エネルギーを改めて $T$ と書く)のひずみトレースは

$$ \sum_\alpha \frac{\partial E}{\partial \varepsilon_{\alpha\alpha}} = -2T - E_{\mathrm{C}}. \label{eq:16-virial-trace} $$

孤立系(分子)の平衡構造では，あらゆる変形に対してエネルギーが停留するから，特にスケーリング変形に対しても $\sum_\alpha \partial E/\partial \varepsilon_{\alpha\alpha} = 0$．すなわち

$$ 2T = -E_{\mathrm{C}} \qquad (\text{平衡構造で}) \label{eq:16-virial} $$

——第2章2.3節で導いたビリアル定理($\bm{F}_n=\bm{0}$ での $2T+U=0$)が，応力のトレースとして再導出された．本章の全Coulombエネルギー $E_{\mathrm{C}}$ は第2章のポテンシャルエネルギー $U$ そのもの($U=E_{\mathrm{C}}$)であり，記号を読み替えれば両者は同じ式である．応力テンソルの定式化がビリアル定理の「異方的な一般化」になっていることが分かる．さらに外圧 $P$ の下では，エンタルピー $E + P\Omega$ の停留条件から $\sum_\alpha [\partial E/\partial\varepsilon_{\alpha\alpha} + P\Omega\,\delta_{\alpha\alpha}] = 0$，すなわち圧力付きビリアル定理 $2T + E_{\mathrm{C}} = 3P\Omega$ が得られる．これは第2章2.7節でGaussの定理から導いたバルク系のビリアル定理 $2T+U=3pV$ と，$U\to E_{\mathrm{C}}$，$p\to P$，$V\to\Omega$ の読み替えで完全に一致する(圧縮 $P\gt0$ で $2T+U\gt0$ という符号も同じ)．

導出:LDA交換相関エネルギーの応力(具体例)

Kohn–Sham理論での応力計算の感触をつかむため，LDAの交換相関項 $E_{xc} = \int e_{xc}(n(\rr))\, \dd^3 r$($e_{xc} = n\, \epsilon_{xc}(n)$ は単位体積あたりのxcエネルギー;第7章)のひずみ微分を計算しよう．スケーリングされた波動関数から作られる密度は，式\eqref{eq:16-scaled-wf}より

$$ n_{\bm{\varepsilon}}(\rr) = \left[ \det(\bm{1}+\bm{\varepsilon}) \right]^{-1} n\!\left( (\bm{1}+\bm{\varepsilon})^{-1}\rr \right) \label{eq:16-scaled-density} $$

である(密度は $\abs{\psi}^2$ の和だから前係数は $\det^{-1}$;全電子数は保存される)．これを代入し，$\bm{u} = (\bm{1}+\bm{\varepsilon})^{-1}\rr$ と置換する($d \equiv \det(\bm{1}+\bm{\varepsilon})$ と略記):

\begin{align} E_{xc}(\bm{\varepsilon}) = \int e_{xc}\!\left( d^{-1}\, n(\bm{u}) \right) d\, \dd^3 u. \label{eq:16-xc-sub} \end{align}

$\varepsilon_{\alpha\beta}$ で微分する．$d$ を通じた依存性が2箇所(体積要素と密度の引数)にあり，$\partial d/\partial \varepsilon_{\alpha\beta} = \delta_{\alpha\beta}$(式\eqref{eq:16-det-expand};$\bm{\varepsilon}=0$ で $d = 1$)を使うと，積の微分法則と合成関数の微分法則により

\begin{align} \left. \frac{\partial E_{xc}}{\partial \varepsilon_{\alpha\beta}} \right|_{\bm{\varepsilon}=0} &= \delta_{\alpha\beta} \int e_{xc}(n(\bm{u}))\, \dd^3 u + \int \frac{\dd e_{xc}}{\dd n} \cdot \left( -\delta_{\alpha\beta}\, n(\bm{u}) \right) \dd^3 u \nonumber\\ &= \delta_{\alpha\beta} \int \left[ e_{xc}(n) - n\, v_{xc}(n) \right] \dd^3 r, \qquad v_{xc} \equiv \frac{\dd e_{xc}}{\dd n}. \label{eq:16-xc-stress} \end{align}

LDAのxc応力は純粋に等方的(圧力のみでせん断応力なし)であることが分かる．これは局所汎関数が密度の値しか見ない(方向を知らない)ことの帰結である．GGAでは $\nabla n$ がひずみとともに変換するため，$\nabla_\alpha n\, \nabla_\beta n$ に比例する異方項が加わる(結果のみ述べる)．

物理的意味:応力で何ができるか

(1) 格子定数の最適化:応力がゼロ(あるいは目標の外圧に等しく)なるまでセルを変形させれば，平衡格子定数が得られる．原子内部座標の力と応力を同時にゼロにする「セル+内部座標の同時最適化」が実務の標準である．(2) 弾性定数:応力のひずみ微分 $C_{\alpha\beta\gamma\delta} = \partial \sigma_{\alpha\beta} / \partial \varepsilon_{\gamma\delta}$ はエネルギーの2階微分であり，微小ひずみを掛けて応力の変化を見れば弾性定数テンソルが計算できる．体積弾性率 $B = -\Omega\, \dd P/\dd \Omega$ もその特別な組合せである．(3) 状態方程式・相転移:圧力下のエンタルピー比較により，構造相転移の転移圧力を予測できる．

例:単位の換算と精度の感覚

応力の原子単位は Hartree/Bohr$^3 = 2.942 \times 10^4$ GPa である．典型的な固体の体積弾性率は $10^2$ GPa のオーダー(例えばSiで約 $98$ GPa，ダイヤモンドで約 $443$ GPa)だから，原子単位で見ると $10^{-3}$ 程度のごく小さな数になる．格子定数を $0.1\%$ 決めたければ応力を $\sim 0.1$ GPa($3 \times 10^{-6}$ a.u.)の精度で計算する必要があり，平面波基底ではカットオフエネルギー不足が「偽の圧力」(Pulay応力:セルが変わると基底関数の個数が変わることに起因する，応力版のPulay補正)として現れることに注意する．局在基底では逆に，基底が原子と一緒に変形するため，重なり行列のひずみ微分に由来するPulay型の応力項を解析的に組み込む．

16.4 構造最適化アルゴリズム

力と応力が計算できるようになった．残る問題は，それらを使ってどうやって極小点まで歩くかである．以下では原子座標をすべて1本のベクトルに並べて

$$ \bm{x} = (R_{1x}, R_{1y}, R_{1z}, R_{2x}, \dots, R_{Mz})^{\mathsf{T}} \in \mathbb{R}^{3M}, \qquad \bm{g}(\bm{x}) = \frac{\partial E}{\partial \bm{x}} = -\bm{F} \label{eq:16-opt-notation} $$

と書く($\bm{g}$ を勾配ベクトルと呼ぶ．力は勾配の符号を反転したものである)．第一原理計算が1回のSCFで与えてくれるのは，その点での $E$ と $\bm{g}$ だけである．2階微分($3M \times 3M$ のHessian $H$)は，力の数値微分なら $6M$ 回のSCF，摂動論(密度汎関数摂動論)なら専用の追加計算を要する高価な情報であり，「あるものとして使う」わけにはいかない．この「値と勾配は安い，2階微分は高い」という非対称性が，以下のアルゴリズム設計をすべて支配する．

数学ノート:多変数Taylor展開・2次形式・条件数

滑らかな関数 $E(\bm{x})$ を点 $\bm{x}_0$ のまわりで2次まで展開すると，変位 $\bm{\Delta} = \bm{x} - \bm{x}_0$ について

$$ E(\bm{x}_0 + \bm{\Delta}) = E_0 + \sum_{i=1}^{3M} \frac{\partial E}{\partial x_i}\bigg|_0 \Delta_i + \frac{1}{2} \sum_{i,j=1}^{3M} \frac{\partial^2 E}{\partial x_i \partial x_j}\bigg|_0 \Delta_i \Delta_j + O(\Delta^3) = E_0 + \bm{g}_0^{\mathsf{T}} \bm{\Delta} + \frac{1}{2} \bm{\Delta}^{\mathsf{T}} H_0 \bm{\Delta} + O(\Delta^3) \label{eq:16-taylor2} $$

である．$H$ は式\eqref{eq:16-hessian-def}のHessianで，偏微分の順序交換可能性から実対称行列である．実対称行列は直交する固有ベクトル $\bm{v}_k$ の完全系と実固有値 $h_k$ を持つから，$\bm{\Delta} = \sum_k c_k \bm{v}_k$ と展開すれば2次形式は

$$ \bm{\Delta}^{\mathsf{T}} H \bm{\Delta} = \sum_k h_k c_k^2 \label{eq:16-quadform} $$

と対角化される．すなわちPESは極小点の近傍では，$3M$ 本の独立な「放物線の谷」が直交して重なったものに見える．$h_k$ はモード $k$ の実効的なばね定数であり，質量で重み付けしたHessian $D_{I\alpha,J\beta} = H_{I\alpha,J\beta}/\sqrt{M_I M_J}$ の固有値は基準振動の角振動数の2乗 $\omega_k^2$ を与える(16.6節で時間刻みを決めるときに使う)．最後に，正定値なHessianの最大固有値と最小固有値の比

$$ \kappa = \frac{h_{\max}}{h_{\min}} \;(\ge 1) \label{eq:16-condnum} $$

を条件数と呼ぶ．$\kappa = 1$ なら等高線は真円，$\kappa \gg 1$ なら細長い谷である．以下で見るように，素朴なアルゴリズムの収束速度はこの $\kappa$ ただ一つでほぼ決まってしまう．

16.4.1 最急降下法とその遅さ

最も素朴な方法は，力の向きにそのまま進むことである:

$$ \bm{x}_{n+1} = \bm{x}_n - \alpha_n\, \bm{g}_n = \bm{x}_n + \alpha_n \bm{F}_n, \qquad \alpha_n > 0. \label{eq:16-sd-update} $$

これを最急降下法(steepest descent)と呼ぶ．$-\bm{g}$ はその点で $E$ が最も速く減る向きだから，少なくとも $\alpha$ を十分小さく取れば必ずエネルギーは下がる．歩幅 $\alpha_n$ の決め方として，その直線上でエネルギーを最小にする「厳密直線探索」を採用しよう．式\eqref{eq:16-taylor2}を2次で打ち切った模型では

\begin{align} E(\bm{x}_n - \alpha \bm{g}_n) &= E_n - \alpha\, \bm{g}_n^{\mathsf{T}} \bm{g}_n + \frac{\alpha^2}{2}\, \bm{g}_n^{\mathsf{T}} H \bm{g}_n, \label{eq:16-sd-line} \end{align}

となる(式\eqref{eq:16-taylor2}に $\bm{\Delta} = -\alpha\bm{g}_n$ を代入して展開しただけである)．$\alpha$ で微分してゼロと置くと $-\bm{g}_n^{\mathsf{T}}\bm{g}_n + \alpha\, \bm{g}_n^{\mathsf{T}} H \bm{g}_n = 0$，すなわち

$$ \alpha_n = \frac{\bm{g}_n^{\mathsf{T}} \bm{g}_n}{\bm{g}_n^{\mathsf{T}} H \bm{g}_n} \label{eq:16-sd-alpha} $$

である($H$ が正定値なら分母は正だから，これは確かに最小値を与える)．ここからが本題である．この「最良の歩幅」を使ってさえ，最急降下法は細長い谷で絶望的に遅い．それを2変数の2次形式で厳密に示す．

導出:最急降下法の収束率は $\left(\dfrac{\kappa-1}{\kappa+1}\right)^n$

Step 1:模型を設定する．Hessianの固有基底を座標に取り，2変数だけを残す(他の方向は同じ議論の繰り返しである):

$$ E(x_1, x_2) = \frac{1}{2}\left( \lambda_1 x_1^2 + \lambda_2 x_2^2 \right), \qquad \lambda_1 \ge \lambda_2 > 0, \qquad \kappa \equiv \frac{\lambda_1}{\lambda_2}. \label{eq:16-sd-model} $$

極小は原点，勾配は $\bm{g} = (\lambda_1 x_1, \lambda_2 x_2)^{\mathsf{T}}$ である．

Step 2:特別な出発点を選ぶ．$\bm{x}_0 = (1/\lambda_1,\; 1/\lambda_2)^{\mathsf{T}}$ から出発する．このとき勾配は $\bm{g}_0 = (1, 1)^{\mathsf{T}}$ である．式\eqref{eq:16-sd-alpha}に代入すると，$\bm{g}_0^{\mathsf{T}}\bm{g}_0 = 2$，$\bm{g}_0^{\mathsf{T}} H \bm{g}_0 = \lambda_1 \cdot 1^2 + \lambda_2 \cdot 1^2 = \lambda_1 + \lambda_2$ だから

$$ \alpha_0 = \frac{2}{\lambda_1 + \lambda_2}. \label{eq:16-sd-alpha0} $$

Step 3:1歩進める．各成分は $x_i \to x_i(1 - \alpha_0 \lambda_i)$ と変換される(式\eqref{eq:16-sd-update}に $g_i = \lambda_i x_i$ を代入した)．倍率を計算すると

\begin{align} 1 - \alpha_0 \lambda_1 &= 1 - \frac{2\lambda_1}{\lambda_1 + \lambda_2} = \frac{\lambda_2 - \lambda_1}{\lambda_1 + \lambda_2} = -\frac{\kappa - 1}{\kappa + 1}, \label{eq:16-sd-fac1}\\[2pt] 1 - \alpha_0 \lambda_2 &= 1 - \frac{2\lambda_2}{\lambda_1 + \lambda_2} = \frac{\lambda_1 - \lambda_2}{\lambda_1 + \lambda_2} = +\frac{\kappa - 1}{\kappa + 1}. \label{eq:16-sd-fac2} \end{align}

最後の等号では分子・分母を $\lambda_2$ で割り，$\lambda_1/\lambda_2 = \kappa$ を使った．以下 $r \equiv (\kappa-1)/(\kappa+1)$ と略記する．よって

$$ \bm{x}_1 = r \left( -\frac{1}{\lambda_1},\; +\frac{1}{\lambda_2} \right)^{\mathsf{T}}. \label{eq:16-sd-x1} $$

Step 4:次の1歩も同じであることを確認する．$\bm{x}_1$ での勾配は $\bm{g}_1 = r(-1, +1)^{\mathsf{T}}$ である．$\bm{g}_1^{\mathsf{T}}\bm{g}_1 = 2r^2$，$\bm{g}_1^{\mathsf{T}} H \bm{g}_1 = r^2(\lambda_1 + \lambda_2)$ だから，歩幅は $\alpha_1 = 2/(\lambda_1+\lambda_2) = \alpha_0$ で変わらない．したがって倍率も同じで，

$$ \bm{x}_2 = r^2 \left( \frac{1}{\lambda_1},\; \frac{1}{\lambda_2} \right)^{\mathsf{T}} = r^2\, \bm{x}_0. \label{eq:16-sd-x2} $$

2歩で出発点の $r^2$ 倍に相似縮小した．以降はこの相似変換の繰り返しだから，厳密に

$$ \abs{\bm{x}_n} \propto r^n = \left( \frac{\kappa - 1}{\kappa + 1} \right)^{n}, \qquad E(\bm{x}_n) = r^{2n}\, E(\bm{x}_0) \label{eq:16-sd-rate} $$

が成り立つ．これで主張が示された．$x_1$ 成分(硬い方向)の符号が毎回反転し，$x_2$ 成分(軟らかい方向)の符号は変わらないことに注意しよう．これが「谷を横切って何度も往復しながら，谷底に沿ってはなかなか進まない」というジグザグ運動の正体である．

この結果を数値で味わっておこう．誤差を $10^{-3}$ 倍に落とすのに必要な歩数 $n$ は $r^n = 10^{-3}$ から $n = 3\ln 10 / \ln(1/r)$ である．$\kappa \gg 1$ では $r = (1 - 1/\kappa)/(1 + 1/\kappa) \simeq 1 - 2/\kappa$，したがって $\ln(1/r) \simeq 2/\kappa$ であり

$$ n \simeq \frac{3 \ln 10}{2}\, \kappa \simeq 3.5\,\kappa. \label{eq:16-sd-steps} $$

歩数が条件数に比例してしまう．分子系での $\kappa$ を見積もると，最も硬いC–H伸縮のばね定数は $\sim 0.35$ Hartree/Bohr$^2$，最も軟らかい二面角のねじれに対する実効的なばね定数はその $10^{-2}$〜$10^{-3}$ 倍であるから $\kappa \sim 10^2$〜$10^3$，すなわち数百〜数千ステップ——1ステップが1回のSCF計算であることを思えば，実用にならない．

最急降下法のジグザグと準Newton法の比較．条件数16の2次形式の等高線の上に，最急降下法（赤の実線）が谷を横切って往復しながら極小点へ近づく経路と，準Newton法（緑の破線）が2歩で極小点に達する経路を描いた図． — 図16.3 細長い谷における最急降下法(実線・赤)と準Newton法(破線・緑)の比較．最急降下法の各ステップは等高線に垂直で，連続する2ステップは互いに直交する．その結果，谷を横切る往復運動になり，谷底方向にはわずかしか進まない．準Newton法は過去の勾配の変化から谷の曲率(Hessian)を推定し，谷底へ直接向かうステップを作る．

物理的意味:条件数を下げるという発想

式\eqref{eq:16-sd-steps}は「悪いのはアルゴリズムではなく座標系である」とも読める．実際，$\bm{x}' = H^{1/2}\bm{x}$ という座標変換を行えばHessianは単位行列になり($\kappa = 1$)，最急降下法は1歩で極小に到達する．もちろん $H$ を知っていればNewton法が使えるのだから循環論法だが，「おおよその $H$」で十分だという点が重要である．これを前処理(preconditioning)と呼ぶ．分子計算で結合長・結合角・二面角からなる内部座標を使うと収束が劇的に速くなるのは，内部座標がまさに近似的な前処理になっている(硬い自由度と軟らかい自由度が座標として分離される)からである．次項以降の準Newton法は，この前処理を計算の途中で自動的に学習する方法だと言ってよい．

16.4.2 Newton法:曲率を使えば2次収束する

2次までのTaylor展開\eqref{eq:16-taylor2}を「模型」とみなし，模型の極小へ一気に飛ぶことを考える．$\bm{\Delta}$ について微分してゼロと置くと

\begin{align} \frac{\partial}{\partial \bm{\Delta}} \left[ E_0 + \bm{g}^{\mathsf{T}}\bm{\Delta} + \frac{1}{2}\bm{\Delta}^{\mathsf{T}} H \bm{\Delta} \right] = \bm{g} + H \bm{\Delta} = \bm{0} \qquad \Longrightarrow \qquad H \bm{\Delta} = -\bm{g}. \label{eq:16-newton-eq} \end{align}

($H$ が対称であることから2次形式の微分は $H\bm{\Delta}$ になる．成分で書けば $\partial/\partial\Delta_k \left[\frac{1}{2}\sum_{ij}\Delta_i H_{ij}\Delta_j\right] = \frac{1}{2}\sum_j H_{kj}\Delta_j + \frac{1}{2}\sum_i \Delta_i H_{ik} = \sum_j H_{kj}\Delta_j$ である．)これを解いて

$$ \bm{x}_{n+1} = \bm{x}_n - H_n^{-1} \bm{g}_n \label{eq:16-newton-update} $$

を得る．これがNewton法(正確にはNewton–Raphson法)である．PESが厳密に2次形式なら1歩で厳密解に到達する．一般のPESでも，極小点に十分近ければ驚くほど速い．それを示そう．

導出:Newton法の2次収束

本質は1変数で尽きるので，$E(x)$，$g = E'$，$h = E''$ として1変数で示す．極小点を $x^*$($g(x^*) = 0$，$h(x^*) \neq 0$)，誤差を $e_n = x_n - x^*$ とする．更新式は $x_{n+1} = x_n - g(x_n)/h(x_n)$ である．

Step 1:分子と分母を $x^*$ のまわりで展開する．

\begin{align} g(x_n) &= g(x^*) + h(x^*) e_n + \frac{1}{2} g''(x^*) e_n^2 + O(e_n^3) = h^* e_n + \frac{1}{2} g''^{*} e_n^2 + O(e_n^3), \label{eq:16-newton-gexp}\\[2pt] h(x_n) &= h(x^*) + g''(x^*) e_n + O(e_n^2) = h^* + g''^{*} e_n + O(e_n^2). \label{eq:16-newton-hexp} \end{align}

1行目で $g(x^*) = 0$ を使った(これが極小点の定義である)．$g'' = E'''$ は3階微分である．

Step 2:誤差の更新式を作る．両辺から $x^*$ を引くと $e_{n+1} = e_n - g(x_n)/h(x_n)$ だから

\begin{align} e_{n+1} = e_n - \frac{h^* e_n + \frac{1}{2} g''^{*} e_n^2}{h^* + g''^{*} e_n} = e_n \cdot \frac{\left( h^* + g''^{*} e_n \right) - \left( h^* + \frac{1}{2} g''^{*} e_n \right)}{h^* + g''^{*} e_n} = e_n \cdot \frac{\frac{1}{2} g''^{*} e_n}{h^* + g''^{*} e_n}. \label{eq:16-newton-err} \end{align}

2つ目の等号では $e_n$ を括り出し，通分して分子をまとめた($e_n \cdot \frac{h^*+g''^*e_n}{h^*+g''^*e_n}$ から分数を引いた)．3つ目の等号は分子の引き算を実行しただけである．

Step 3:極限を取る．$e_n \to 0$ では分母は $h^*$ に収束するから

$$ e_{n+1} \simeq \frac{g'''(x^*)}{2\,E''(x^*)}\, e_n^2 \qquad\text{すなわち}\qquad \abs{e_{n+1}} \le C \abs{e_n}^2. \label{eq:16-newton-quad} $$

誤差が毎回2乗される．これを2次収束と呼ぶ．有効数字で言えば「正しい桁数が1ステップごとに倍になる」ということで，$10^{-2} \to 10^{-4} \to 10^{-8} \to 10^{-16}$ と，数ステップで機械精度に達する．多変数の場合も，$H$ が極小点で正定値かつLipschitz連続であれば同じ評価 $\abs{\bm{e}_{n+1}} \le C\abs{\bm{e}_n}^2$ が成り立つ．

ではなぜ実際の構造最適化でNewton法をそのまま使わないのか．理由は2つある．

(i) Hessianが高価である．$3M \times 3M$ の $H$ を数値微分で作るには $6M$ 回のSCF計算が必要で，原子50個なら300回——最適化そのものより高くつく．密度汎関数摂動論を使えば1回の追加計算で得られるが，実装は重い．

(ii) Hessianが不定値だと極小へ行かない．これは概念的に重要な点なので詳しく見る．$H$ を固有分解して $H = \sum_k h_k \bm{v}_k \bm{v}_k^{\mathsf{T}}$，勾配を $\bm{g} = \sum_k g_k \bm{v}_k$($g_k = \bm{v}_k^{\mathsf{T}}\bm{g}$)と展開すると，Newton変位は

$$ \bm{\Delta} = -H^{-1}\bm{g} = -\sum_k \frac{g_k}{h_k} \bm{v}_k, \qquad \Delta E \simeq \bm{g}^{\mathsf{T}}\bm{\Delta} + \frac{1}{2}\bm{\Delta}^{\mathsf{T}} H \bm{\Delta} = -\sum_k \frac{g_k^2}{h_k} + \frac{1}{2}\sum_k \frac{g_k^2}{h_k} = -\sum_k \frac{g_k^2}{2 h_k} \label{eq:16-newton-modes} $$

となる(2行目は式\eqref{eq:16-newton-modes}第1式を代入し，$\bm{v}_k$ の直交規格性 $\bm{v}_k^{\mathsf{T}}\bm{v}_l = \delta_{kl}$ を使って和を実行した)．$h_k > 0$ のモードでは $\Delta E < 0$ でエネルギーが下がるが，$h_k < 0$ のモードでは $\Delta E > 0$ でエネルギーが上がる．すなわちNewton法は「勾配がゼロの点」に向かうだけで，極小に向かうとは限らない——鞍点にも同じ勢いで収束する．出発構造が悪いと $H$ は容易に負の固有値を持つから，これは現実的な危険である．

対策は，負(または小さすぎる)固有値を持ち上げることである．実務では次のレベルシフト付きの方程式を解く:

$$ \left( H + \mu\, \bm{1} \right) \bm{\Delta} = -\bm{g}, \qquad \mu \ge 0 \;\text{は}\; H + \mu\bm{1} \succ 0 \;\text{となるように選ぶ}. \label{eq:16-levelshift} $$

$\mu$ を大きくすると $\bm{\Delta} \to -\bm{g}/\mu$ すなわち最急降下法に連続的に移行する．この一見場当たり的な処方には，次に述べる信頼領域という明快な意味がある．

導出:信頼半径の拘束がレベルシフトを生む

2次のTaylor模型が信用できるのは $\bm{x}_n$ の近傍だけである．そこで「模型を信じる半径」$R$(信頼半径)を決め，その内側でのみ模型を最小化する:

$$ \min_{\bm{\Delta}} \; m(\bm{\Delta}) = \bm{g}^{\mathsf{T}}\bm{\Delta} + \frac{1}{2}\bm{\Delta}^{\mathsf{T}} H \bm{\Delta} \qquad \text{s.t.} \qquad \abs{\bm{\Delta}}^2 \le R^2. \label{eq:16-tr-problem} $$

拘束が有効(等号で成立)な場合，Lagrangeの未定乗数法(第4章)を使う．未定乗数を $\mu$ として

$$ \mathcal{L}(\bm{\Delta}, \mu) = \bm{g}^{\mathsf{T}}\bm{\Delta} + \frac{1}{2}\bm{\Delta}^{\mathsf{T}} H \bm{\Delta} + \frac{\mu}{2}\left( \abs{\bm{\Delta}}^2 - R^2 \right) \label{eq:16-tr-lag} $$

を作り，$\bm{\Delta}$ で微分してゼロと置くと

\begin{align} \frac{\partial \mathcal{L}}{\partial \bm{\Delta}} = \bm{g} + H\bm{\Delta} + \mu \bm{\Delta} = \bm{0} \qquad \Longrightarrow \qquad (H + \mu \bm{1}) \bm{\Delta} = -\bm{g} \label{eq:16-tr-result} \end{align}

となり，式\eqref{eq:16-levelshift}が導かれた．$\mu$ は拘束条件 $\abs{\bm{\Delta}(\mu)} = R$ から決まる．$\abs{\bm{\Delta}(\mu)}$ は $\mu$ の単調減少関数だから，$R$ を小さくすれば $\mu$ が大きくなる——「模型を信用しないほど，最急降下法に近づく」という直観そのものである．実際のコードは，1ステップ進むごとに予測されたエネルギー低下 $m(\bm{\Delta})$ と実際の低下 $E(\bm{x}_n + \bm{\Delta}) - E(\bm{x}_n)$ の比を見て，比が1に近ければ $R$ を拡大，大きく食い違えば $R$ を縮小してステップをやり直す，という適応制御を行う．

なお，この枠組みで $\mu$ の選び方を変えると鞍点を狙うことすらできる．ある1つのモード $k_0$ についてだけ $h_{k_0} + \mu < 0$ となるように選べば，そのモードに沿ってはエネルギーを上り，他のすべてのモードでは下る——これが固有ベクトル追跡法(eigenvector following)であり，遷移状態探索の一手法である(16.5節のNEB法とは相補的に使われる)．

16.4.3 準Newton法とBFGS更新

Hessianは高価だが，実はただで手に入る情報がある．2点 $\bm{x}_k$，$\bm{x}_{k+1}$ で勾配を計算したなら，その差は曲率の情報を含んでいるからである．実際，$E$ が2次形式なら $\bm{g}(\bm{x}) = \bm{g}(\bm{x}^*) + H(\bm{x} - \bm{x}^*)$ だから，

$$ \bm{s}_k \equiv \bm{x}_{k+1} - \bm{x}_k, \qquad \bm{y}_k \equiv \bm{g}_{k+1} - \bm{g}_k \qquad \Longrightarrow \qquad \bm{y}_k = H \bm{s}_k \label{eq:16-secant} $$

が厳密に成り立つ．一般のPESでもこれは「$\bm{s}_k$ 方向の平均曲率」を与える．そこでHessianの近似行列 $B_k$ を逐次更新していき，少なくとも直前のステップについてはこの関係を再現するよう要求する:

$$ B_{k+1}\, \bm{s}_k = \bm{y}_k . \label{eq:16-secant-cond} $$

これをセカント条件と呼ぶ．これが準Newton法(quasi-Newton method)の出発点である．ただしセカント条件は $3M$ 本の方程式にすぎず，$B_{k+1}$ の $3M(3M+1)/2$ 個の自由度を決めきれない．残りは「対称性を保つ」「前の $B_k$ からの変更を最小にする」という2つの要請で埋める．前者は $B$ がHessianの近似である以上当然であり，後者は「新しく得た情報($\bm{s}_k, \bm{y}_k$ の1組)だけを反映し，それ以外は前の知識を保つ」という自然な原理である．実際にやってみよう．

導出:BFGS更新式

Step 1:変更の形を決める．「最小の変更」を具体化するため，更新を対称なランク2の行列に限る．ランク1($\bm{u}\bm{u}^{\mathsf{T}}$ の形，1個の方向しか変えない)ではセカント条件と正定値性を両立させにくく，ランク3以上は自由度が過剰だからである．すなわち

$$ B_{k+1} = B_k + a\, \bm{u}\bm{u}^{\mathsf{T}} + b\, \bm{v}\bm{v}^{\mathsf{T}} \label{eq:16-bfgs-ansatz} $$

と置く($\bm{u}\bm{u}^{\mathsf{T}}$ は成分で $(\bm{u}\bm{u}^{\mathsf{T}})_{ij} = u_i u_j$ の対称行列である)．以下，添字 $k$ を省いて $\bm{s}, \bm{y}, B$ と書く．

Step 2:セカント条件を課す．式\eqref{eq:16-bfgs-ansatz}を式\eqref{eq:16-secant-cond}に代入すると

$$ B\bm{s} + a\, \bm{u}\, (\bm{u}^{\mathsf{T}}\bm{s}) + b\, \bm{v}\, (\bm{v}^{\mathsf{T}}\bm{s}) = \bm{y}. \label{eq:16-bfgs-sec} $$

ここで $\bm{u}^{\mathsf{T}}\bm{s}$，$\bm{v}^{\mathsf{T}}\bm{s}$ はスカラーだから，左辺は「$B\bm{s}$，$\bm{u}$，$\bm{v}$ の線形結合」である．右辺 $\bm{y}$ を作り，かつ余計な $B\bm{s}$ を消すには，$\bm{u} = \bm{y}$，$\bm{v} = B\bm{s}$ と選ぶのが最も自然である(他の選び方では $B\bm{s}$ を打ち消せない)．代入すると

$$ B\bm{s} + a\, (\bm{y}^{\mathsf{T}}\bm{s})\, \bm{y} + b\, (\bm{s}^{\mathsf{T}}B\bm{s})\, B\bm{s} = \bm{y}. \label{eq:16-bfgs-sec2} $$

Step 3:係数を決める．$\bm{y}$ と $B\bm{s}$ は一般に独立なベクトルだから，両辺の係数を比較して

\begin{align} \bm{y}\;\text{の係数}: \quad & a\,(\bm{y}^{\mathsf{T}}\bm{s}) = 1 &&\Longrightarrow\quad a = \frac{1}{\bm{y}^{\mathsf{T}}\bm{s}}, \label{eq:16-bfgs-a}\\[2pt] B\bm{s}\;\text{の係数}: \quad & 1 + b\,(\bm{s}^{\mathsf{T}}B\bm{s}) = 0 &&\Longrightarrow\quad b = -\frac{1}{\bm{s}^{\mathsf{T}}B\bm{s}}. \label{eq:16-bfgs-b} \end{align}

Step 4:結果をまとめる．式\eqref{eq:16-bfgs-a}, \eqref{eq:16-bfgs-b}を式\eqref{eq:16-bfgs-ansatz}に戻して，添字を復活させる:

$$ B_{k+1} = B_k + \frac{\bm{y}_k \bm{y}_k^{\mathsf{T}}}{\bm{y}_k^{\mathsf{T}} \bm{s}_k} - \frac{B_k \bm{s}_k\, \bm{s}_k^{\mathsf{T}} B_k}{\bm{s}_k^{\mathsf{T}} B_k \bm{s}_k}. \label{eq:16-bfgs} $$

これがBFGS更新式である(Broyden, Fletcher, Goldfarb, Shannoが1970年に独立に到達した)．第2項が「新しい曲率情報 $\bm{y}$ を加える」働きを，第3項が「$\bm{s}$ 方向についての古い情報 $B_k\bm{s}_k$ を取り除く」働きをしている．$B_k$ が対称なら $B_{k+1}$ も明らかに対称である．

BFGS更新の決定的な長所は，正定値性を自動的に保つことである．これを証明しておこう．

定理:BFGS更新は正定値性を保存する

$B_k$ が正定値で，かつ曲率条件 $\bm{y}_k^{\mathsf{T}}\bm{s}_k > 0$ が成り立つならば，式\eqref{eq:16-bfgs}で作られる $B_{k+1}$ も正定値である．

導出:正定値性の証明

任意の $\bm{z} \neq \bm{0}$ に対して $\bm{z}^{\mathsf{T}} B_{k+1} \bm{z} > 0$ を示す．式\eqref{eq:16-bfgs}を $\bm{z}$ で挟むと(以下 $B = B_k$，$\bm{s}=\bm{s}_k$，$\bm{y}=\bm{y}_k$)

$$ \bm{z}^{\mathsf{T}} B_{k+1} \bm{z} = \underbrace{\left[ \bm{z}^{\mathsf{T}} B \bm{z} - \frac{(\bm{s}^{\mathsf{T}} B \bm{z})^2}{\bm{s}^{\mathsf{T}} B \bm{s}} \right]}_{\displaystyle \equiv A} + \underbrace{\frac{(\bm{y}^{\mathsf{T}}\bm{z})^2}{\bm{y}^{\mathsf{T}}\bm{s}}}_{\displaystyle \equiv C}. \label{eq:16-bfgs-pd1} $$

($\bm{z}^{\mathsf{T}}\bm{y}\bm{y}^{\mathsf{T}}\bm{z} = (\bm{y}^{\mathsf{T}}\bm{z})^2$，$\bm{z}^{\mathsf{T}} B\bm{s}\,\bm{s}^{\mathsf{T}}B\bm{z} = (\bm{s}^{\mathsf{T}}B\bm{z})^2$ を使った．$B$ の対称性が効いている．)

第1項 $A$ の符号．$B$ は正定値だから $\langle \bm{a}, \bm{b}\rangle_B \equiv \bm{a}^{\mathsf{T}} B \bm{b}$ は内積の公理を満たす．この内積についてのCauchy–Schwarzの不等式

$$ (\bm{s}^{\mathsf{T}} B \bm{z})^2 \le (\bm{s}^{\mathsf{T}} B \bm{s})(\bm{z}^{\mathsf{T}} B \bm{z}) \label{eq:16-bfgs-cs} $$

より $A \ge 0$ である．等号が成り立つのは $\bm{z}$ と $\bm{s}$ が平行なとき，すなわち $\bm{z} = c\,\bm{s}$($c \neq 0$)のときに限る．

第2項 $C$ の符号．曲率条件 $\bm{y}^{\mathsf{T}}\bm{s} > 0$ より分母は正，分子は2乗だから $C \ge 0$ である．

場合分けで結論する．$\bm{z}$ が $\bm{s}$ と平行でなければ $A > 0$ なので $A + C > 0$．$\bm{z} = c\bm{s}$ のときは $A = 0$ だが，このとき

$$ C = \frac{(\bm{y}^{\mathsf{T}} c\bm{s})^2}{\bm{y}^{\mathsf{T}}\bm{s}} = \frac{c^2 (\bm{y}^{\mathsf{T}}\bm{s})^2}{\bm{y}^{\mathsf{T}}\bm{s}} = c^2\, \bm{y}^{\mathsf{T}}\bm{s} > 0 \label{eq:16-bfgs-pd2} $$

となる．いずれの場合も $\bm{z}^{\mathsf{T}}B_{k+1}\bm{z} > 0$ であり，証明が完了した．

正定値性が保たれることの意義は大きい．$B_{k+1} \succ 0$ ならば探索方向 $\bm{d} = -B_{k+1}^{-1}\bm{g}$ は必ず降下方向である．実際 $\bm{g}^{\mathsf{T}}\bm{d} = -\bm{g}^{\mathsf{T}} B_{k+1}^{-1}\bm{g} < 0$(正定値行列の逆行列も正定値)だから，$\bm{d}$ 方向に十分小さく進めばエネルギーは必ず下がる．Newton法の「鞍点に落ちる」病がBFGSでは構造的に排除されているのである．

実装では $B^{-1}$ を毎回作り直すのは無駄なので，逆行列そのものを更新する．$G_k \equiv B_k^{-1}$，$\rho_k \equiv 1/(\bm{y}_k^{\mathsf{T}}\bm{s}_k)$ と置くと

$$ G_{k+1} = \left( \bm{1} - \rho_k \bm{s}_k \bm{y}_k^{\mathsf{T}} \right) G_k \left( \bm{1} - \rho_k \bm{y}_k \bm{s}_k^{\mathsf{T}} \right) + \rho_k\, \bm{s}_k \bm{s}_k^{\mathsf{T}}, \qquad \bm{\Delta}_{k+1} = -G_{k+1}\, \bm{g}_{k+1} \label{eq:16-bfgs-inv} $$

である．これがセカント条件の逆版 $G_{k+1}\bm{y}_k = \bm{s}_k$ を満たすことは直接確かめられる．$\rho_k \bm{s}_k^{\mathsf{T}} \bm{y}_k = 1$ に注意して右から $\bm{y}_k$ を掛けると

\begin{align} G_{k+1}\bm{y}_k &= \left( \bm{1} - \rho_k \bm{s}_k \bm{y}_k^{\mathsf{T}} \right) G_k \underbrace{\left( \bm{y}_k - \rho_k \bm{y}_k (\bm{s}_k^{\mathsf{T}}\bm{y}_k) \right)}_{= \bm{y}_k - \bm{y}_k = \bm{0}} + \rho_k \bm{s}_k (\bm{s}_k^{\mathsf{T}}\bm{y}_k) = \bm{0} + \bm{s}_k = \bm{s}_k \label{eq:16-bfgs-inv-check} \end{align}

となる(中括弧の中でスカラー $\rho_k(\bm{s}_k^{\mathsf{T}}\bm{y}_k) = 1$ を使った)．

実務上の要点を3つ挙げておく．

初期Hessian $B_0$ の質が収束を左右する．$B_0 = \bm{1}$ と取れば初手は最急降下法になり，その後 $\bm{s},\bm{y}$ の対を蓄積して徐々に曲率を学習する．しかし化学的な知識を入れた $B_0$——結合している原子対には伸縮のばね定数，結合角には曲げのばね定数を経験式で割り当てる(Schlegelの方法など)——を使えば，収束に要するステップ数はしばしば半分以下になる．「良い前処理を最初から与える」ことの効果である．
大規模系にはL-BFGS．$G_k$ を陽に持つと $O((3M)^2)$ の記憶容量が要る．式\eqref{eq:16-bfgs-inv}は $G_k$ に $(\bm{s},\bm{y})$ の対を作用させる形なので，直近 $m$ 組($m = 5$〜$20$)だけを保存して $G_k \bm{g}$ を再帰的に計算すれば，記憶容量は $O(m \cdot 3M)$ で済む．これを記憶制限BFGS(L-BFGS)と呼び，数千原子の系でも使える．
変種．過去の勾配の線形結合で残差を最小化するDIIS(直接反復部分空間法;第14章のSCF加速と同じ発想を座標空間に適用したもの)，レベルシフトを有理関数模型で自動決定するRF法，Hessianの固有値を監視して符号を制御するEF法など，多くの変種が実装されている．分子・固体の標準的なベンチマークでは，BFGSを土台にレベルシフトを組み合わせた手法が最も安定に $10^{-4}$ Hartree/Bohr まで収束することが知られている．

16.4.4 直線探索・曲率条件・収束判定

準Newton法が方向 $\bm{d}_k = -G_k \bm{g}_k$ を与えたあと，実際にどれだけ進むか($\bm{x}_{k+1} = \bm{x}_k + \alpha_k \bm{d}_k$)を決めるのが直線探索である．厳密に最小化する必要はなく，次の2条件(Wolfe条件)を満たせば十分であることが知られている:

\begin{align} &\text{(i) 十分な減少(Armijo条件)}: && E(\bm{x}_k + \alpha \bm{d}_k) \le E(\bm{x}_k) + c_1 \alpha\, \bm{g}_k^{\mathsf{T}} \bm{d}_k, \label{eq:16-wolfe1}\\[2pt] &\text{(ii) 曲率条件}: && \bm{g}(\bm{x}_k + \alpha \bm{d}_k)^{\mathsf{T}} \bm{d}_k \ge c_2\, \bm{g}_k^{\mathsf{T}} \bm{d}_k, \label{eq:16-wolfe2} \end{align}

ただし $0 < c_1 < c_2 < 1$(典型的には $c_1 = 10^{-4}$，$c_2 = 0.9$)である．(i) は「エネルギーが十分に下がった」ことを，(ii) は「勾配の傾きが十分に緩んだ(まだ急降下の途中で止まっていない)」ことを要求する．(ii) には見逃せない役割がある．$\bm{s}_k = \alpha_k \bm{d}_k$ に注意して式\eqref{eq:16-wolfe2}を書き直すと

\begin{align} \bm{y}_k^{\mathsf{T}} \bm{s}_k = \alpha_k \left[ \bm{g}_{k+1} - \bm{g}_k \right]^{\mathsf{T}} \bm{d}_k \ge \alpha_k \left( c_2 - 1 \right) \bm{g}_k^{\mathsf{T}} \bm{d}_k > 0 \label{eq:16-wolfe-curv} \end{align}

となる．最後の不等号は，$c_2 - 1 < 0$ と，$\bm{d}_k$ が降下方向であること($\bm{g}_k^{\mathsf{T}}\bm{d}_k < 0$)から，負×負で正になることによる．つまりWolfe条件を満たす直線探索は，BFGSの正定値性に必要な曲率条件 $\bm{y}^{\mathsf{T}}\bm{s} > 0$ を自動的に保証する．理論と実装がきれいに噛み合っている場面である．

最後に収束判定について述べる．実務では次の3つを併用する．

最大力:$\max_I \abs{\bm{F}_I} < F_{\mathrm{tol}}$．$F_{\mathrm{tol}}$ の典型値は $10^{-3}$〜$10^{-4}$ Hartree/Bohr である(それぞれ $0.051$，$0.0051$ eV/Åに相当)．
エネルギー変化:$\abs{E_{k+1} - E_k} < 10^{-6}$ Hartree 程度．
最大変位:$\max_I \abs{\Delta \bm{R}_I} < 10^{-3}$ Bohr 程度．

力の閾値を「どこまで厳しくすべきか」は，残留力が構造をどれだけずらすかで決まる．モード $k$(ばね定数 $h_k$)に残留力 $F$ があるとき，真の極小からのずれは $\delta = F/h_k$ である．硬い結合伸縮($h \simeq 0.3$ Hartree/Bohr$^2$)なら $F = 10^{-4}$ Hartree/Bohr に対して $\delta = 3\times10^{-4}$ Bohr $= 2\times10^{-4}$ Å と完全に無視できる．しかし軟らかいモード($h \simeq 5\times10^{-3}$ Hartree/Bohr$^2$，分子間の伸縮や層間距離など)では $\delta = 0.02$ Bohr $= 0.01$ Å になる．閾値の意味は系のばね定数によって全く違う——これは初学者が見落としやすい点である．

例:力の数値ノイズという床

第一原理計算が返す力には，SCFの収束打ち切り，実空間グリッドの離散化(グリッドに対する原子の相対位置でエネルギーがわずかに波打つ「eggbox誤差」)，Brillouin域サンプリングの有限性に由来するノイズが必ず乗っている．その大きさは典型的に $10^{-5}$〜$10^{-6}$ Hartree/Bohr であり，これより厳しい収束判定を課しても意味がない(むしろ最適化が振動して止まらなくなる)．したがって構造最適化の前に，SCF収束条件・カットオフ・$k$ 点数を「力が安定する」ところまで確認しておくのが正しい手順である．特に分子動力学(16.6節)では，力のノイズがエネルギーの系統的ドリフトとして蓄積するため，単発の構造最適化より1〜2桁厳しいSCF収束が要求される．

16.4.5 格子振動とフォノン:Hessianを周期系へ広げる

ここまでHessianは「極小へ速く歩くための道具」として登場した．しかしHessianはそれ自体が物理量である．極小点で評価したHessianは，原子が平衡位置のまわりでどう振動するかを完全に決めており，周期系に対してはフォノン(格子振動の量子)の分散関係を与える．フォノンは第一原理計算の主要な出力の一つであり(第1章)，比熱・熱膨張・熱伝導・赤外およびRamanスペクトル・超伝導転移温度，そして何より得られた構造が本当に安定かどうかの判定に用いられる．本項では調和近似から動的行列までを一段も飛ばさずに導き，1次元2原子鎖で分散関係を解析的に求め，最後に第一原理での2つの計算法を比較する．

(1) 調和近似．平衡配置を $\{\bm{R}_i^{(0)}\}$，そこからの変位を $\bm{u}_i = \bm{R}_i - \bm{R}_i^{(0)}$ と書く($i$ は原子の番号，$\alpha,\beta = x,y,z$)．全エネルギーを変位について2次までTaylor展開する(式\eqref{eq:16-taylor2}と同じ操作を，変数名だけ変えて書き下したものである):

\begin{align} E(\{\bm{R}_i^{(0)} + \bm{u}_i\}) = E_0 &+ \sum_{i\alpha} \left.\frac{\partial E}{\partial R_{i\alpha}}\right|_0 u_{i\alpha} + \frac{1}{2} \sum_{i\alpha}\sum_{j\beta} \left.\frac{\partial^2 E}{\partial R_{i\alpha}\, \partial R_{j\beta}}\right|_0 u_{i\alpha} u_{j\beta} + O(u^3). \label{eq:16-ph-taylor} \end{align}

ここで $|_0$ は平衡配置で評価することを表す．1次の項は消える．理由は定義そのものである．式\eqref{eq:16-force-def}より $\partial E/\partial R_{i\alpha} = -F_{i\alpha}$ であり，平衡配置とは「すべての原子に働く力がゼロの配置」だから $F_{i\alpha}|_0 = 0$，すなわち $\partial E/\partial R_{i\alpha}|_0 = 0$ である．したがって

$$ E = E_0 + \frac{1}{2} \sum_{i\alpha,\, j\beta} \Phi_{i\alpha,\, j\beta}\, u_{i\alpha} u_{j\beta}, \qquad \Phi_{i\alpha,\, j\beta} \equiv \left.\frac{\partial^2 E}{\partial R_{i\alpha}\, \partial R_{j\beta}}\right|_0 \label{eq:16-ph-harmonic} $$

となる．3次以上の項を捨てるこの扱いを調和近似と呼び，$\Phi_{i\alpha,j\beta}$ を力の定数行列(force constant matrix)と呼ぶ．$\Phi$ は式\eqref{eq:16-hessian-def}で定義したHessian $H_{I\alpha,J\beta}$ そのものである——構造最適化の文脈ではHessian，格子振動の文脈では力の定数行列と呼び分ける習慣があるだけで，中身は同じ2階微分係数の行列である．偏微分の順序が交換できることから，$\Phi$ は実対称行列である:

$$ \Phi_{i\alpha,\, j\beta} = \Phi_{j\beta,\, i\alpha}. \label{eq:16-ph-symm} $$

(2) 変位に比例する力．調和近似のもとで原子 $k$ に働く力を求める．式\eqref{eq:16-ph-harmonic}を $u_{k\gamma}$ で微分する．2重和の中で $u_{k\gamma}$ を含む項は「第1の添字が $k\gamma$ のもの」と「第2の添字が $k\gamma$ のもの」の2種類あるから，

\begin{align} \frac{\partial E}{\partial u_{k\gamma}} &= \frac{1}{2}\sum_{j\beta} \Phi_{k\gamma,\, j\beta}\, u_{j\beta} + \frac{1}{2}\sum_{i\alpha} \Phi_{i\alpha,\, k\gamma}\, u_{i\alpha} \label{eq:16-ph-force-step1}\\[2pt] &= \frac{1}{2}\sum_{j\beta} \Phi_{k\gamma,\, j\beta}\, u_{j\beta} + \frac{1}{2}\sum_{j\beta} \Phi_{k\gamma,\, j\beta}\, u_{j\beta} = \sum_{j\beta} \Phi_{k\gamma,\, j\beta}\, u_{j\beta} \label{eq:16-ph-force-step2} \end{align}

となる．2行目では第2項の対称性\eqref{eq:16-ph-symm}を使って添字を入れ替え，和の名前を $i\alpha \to j\beta$ と付け替えただけである．力は勾配の符号反転だから

$$ F_{k\gamma} = -\frac{\partial E}{\partial u_{k\gamma}} = -\sum_{j\beta} \Phi_{k\gamma,\, j\beta}\, u_{j\beta}. \label{eq:16-ph-force-lin} $$

すなわち調和近似とは「力が変位の線形関数である」——連成したばねの集まりである——という近似にほかならない．$\Phi_{k\gamma,j\beta}$ は「原子 $j$ を $\beta$ 方向に単位長さずらしたとき，原子 $k$ が $\gamma$ 方向に受ける力の符号を反転したもの」という直接的な意味を持つ．

(3) 並進不変性から従う総和則．力の定数行列は勝手な値を取れるわけではない．結晶(あるいは分子)全体を剛体的に平行移動しても，原子どうしの相対配置は何も変わらないから，エネルギーも各原子に働く力も一切変わらない．この当たり前の事実が $\Phi$ に厳しい拘束を課す．剛体並進を $u_{j\beta} = c_\beta$(すべての $j$ について同じ定数ベクトル)と置いて式\eqref{eq:16-ph-force-lin}に代入すると，力はゼロのままでなければならない:

\begin{align} 0 = F_{i\alpha} &= -\sum_{j\beta} \Phi_{i\alpha,\, j\beta}\, c_\beta \label{eq:16-ph-asr-step1}\\[2pt] &= -\sum_{\beta} c_\beta \left( \sum_{j} \Phi_{i\alpha,\, j\beta} \right). \label{eq:16-ph-asr-step2} \end{align}

2行目では，$c_\beta$ が $j$ に依存しないので $j$ についての和の外に出し，和の順序を入れ替えた．これが任意の $c_\beta$ について成り立たねばならないから，括弧の中が各 $\beta$ について独立にゼロである:

$$ \sum_{j} \Phi_{i\alpha,\, j\beta} = 0 \qquad (\text{すべての } i,\ \alpha,\ \beta \text{ について}). \label{eq:16-ph-asr} $$

これを音響総和則(acoustic sum rule, ASR)と呼ぶ．$j = i$ の項を分離すれば $\Phi_{i\alpha,i\beta} = -\sum_{j \neq i}\Phi_{i\alpha,j\beta}$ となり，「自分自身との力の定数(対角ブロック)は，他のすべての原子との力の定数の和で決まる」ことを意味する．後で見るように，この総和則が音響分枝の存在($\bm{q}\to\bm{0}$ で $\omega\to0$)を保証する．実際の数値計算では離散グリッドの誤差などで\eqref{eq:16-ph-asr}がわずかに破れるため，対角ブロックを総和則が厳密に成り立つよう事後的に補正するのが標準的な作法である．

(4) 運動方程式と動的行列．核を古典粒子として扱う(Born–Oppenheimer近似のもとでの核の運動，16.6節と同じ立場である)．質量 $M_i$ の原子 $i$ についてNewtonの運動方程式を書き，右辺に\eqref{eq:16-ph-force-lin}を代入すると

$$ M_i\, \ddot{u}_{i\alpha}(t) = F_{i\alpha} = -\sum_{j\beta} \Phi_{i\alpha,\, j\beta}\, u_{j\beta}(t) \label{eq:16-ph-newton} $$

である．これは連成振動の方程式であり，原子数が多ければ巨大な連立微分方程式になる．ところが周期系ではこれが波数 $\bm{q}$ ごとの小さな固有値問題に分解する．以下それを示す．

まず添字を周期系に合わせて付け替える．単位胞を格子ベクトル $\bm{R}_l$ で区別し($l$ は胞の番号)，胞内の原子を $\kappa = 1,\dots,n$ で区別する($n$ は単位胞あたりの原子数)．原子 $(l,\kappa)$ の平衡位置は $\bm{R}_l + \bm{\tau}_\kappa$，質量は $M_\kappa$(胞の番号によらない)，変位は $u_{l\kappa\alpha}$ である．結晶の並進対称性から，力の定数は2つの胞の差だけに依存する:

$$ \Phi_{l\kappa\alpha,\, l'\kappa'\beta} = \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R}_{l'} - \bm{R}_l). \label{eq:16-ph-periodic} $$

(胞 $l$ と胞 $l'$ を同時に同じ格子ベクトルだけずらしても結晶は自分自身に重なるから，$\Phi$ は差 $\bm{R}_{l'}-\bm{R}_l$ の関数である．)これを使うと運動方程式\eqref{eq:16-ph-newton}は

$$ M_\kappa\, \ddot{u}_{l\kappa\alpha} = -\sum_{l'}\sum_{\kappa'\beta} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R}_{l'} - \bm{R}_l)\, u_{l'\kappa'\beta} \label{eq:16-ph-eom-cell} $$

と書ける．

数学ノート:並進対称性・平面波型の解・エルミート行列

式\eqref{eq:16-ph-eom-cell}の右辺は「胞の差」だけに依存する係数と変位との畳み込みの形をしている．第12章でBlochの定理を学んだときと事情はまったく同じで，係数行列が離散並進に対して不変なら，解を平面波因子 $e^{i\bm{q}\cdot\bm{R}_l}$ で分類できる．電子の場合は波動関数がBloch形をとり，核の変位の場合は変位パターンが平面波型をとる．どちらも「並進対称な線形問題は波数で対角化される」という一つの事実の現れである．

もう一つ用意しておくべき道具がエルミート行列である．複素行列 $A$ が $A_{mn}^{*} = A_{nm}$(転置して複素共役を取ると自分に戻る)を満たすとき $A$ をエルミートと呼ぶ．エルミート行列の固有値は必ず実数であり，固有ベクトルは(縮退があっても取り方を工夫すれば)複素内積 $\braket{\bm{a},\bm{b}} = \sum_m a_m^{*} b_m$ について直交させられる．実対称行列は実数の範囲でのエルミート行列にほかならない．以下に現れる動的行列は複素行列だがエルミートであり，そのおかげで固有値 $\omega^2$ が実数であることが保証される．ただし実数であっても正とは限らない——この一点が(5)の動的安定性の議論の要になる．

そこで，質量の平方根で重み付けした平面波型の解を代入してみる:

$$ u_{l\kappa\alpha}(t) = \frac{1}{\sqrt{M_\kappa}}\, e_{\kappa\alpha}(\bm{q})\, e^{i(\bm{q}\cdot\bm{R}_l - \omega t)}. \label{eq:16-ph-ansatz} $$

$e_{\kappa\alpha}(\bm{q})$ は胞内の原子ごとの振幅と位相を表す複素数($3n$ 個)であり，これから決めるべき未知数である．$1/\sqrt{M_\kappa}$ という因子を付けておくのは，後で行列がエルミートになるようにするためのお膳立てである(付けなくても物理は同じだが，行列が非対称になって扱いにくい)．代入を1行ずつ実行する．

導出:動的行列と固有値問題

Step 1:左辺を計算する．式\eqref{eq:16-ph-ansatz}を時間で2回微分すると，時間依存は $e^{-i\omega t}$ だけだから $(-i\omega)^2 = -\omega^2$ が出る:

$$ M_\kappa\, \ddot{u}_{l\kappa\alpha} = M_\kappa \cdot \frac{1}{\sqrt{M_\kappa}} e_{\kappa\alpha}\, (-i\omega)^2 e^{i(\bm{q}\cdot\bm{R}_l - \omega t)} = -\omega^2 \sqrt{M_\kappa}\, e_{\kappa\alpha}\, e^{i(\bm{q}\cdot\bm{R}_l - \omega t)}. \label{eq:16-ph-lhs} $$

($M_\kappa/\sqrt{M_\kappa} = \sqrt{M_\kappa}$ を使った．)

Step 2:右辺に代入する．式\eqref{eq:16-ph-eom-cell}の右辺の $u_{l'\kappa'\beta}$ に\eqref{eq:16-ph-ansatz}を入れる:

$$ -\sum_{l'}\sum_{\kappa'\beta} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R}_{l'} - \bm{R}_l)\, \frac{1}{\sqrt{M_{\kappa'}}}\, e_{\kappa'\beta}\, e^{i(\bm{q}\cdot\bm{R}_{l'} - \omega t)}. \label{eq:16-ph-rhs} $$

Step 3:和の変数を差に取り替える．$\bm{R} \equiv \bm{R}_{l'} - \bm{R}_l$ と置く．$l$ を固定して $l'$ をすべての胞にわたって走らせることは，$\bm{R}$ をすべての格子ベクトルにわたって走らせることと同じである(無限結晶，あるいは周期境界条件のもとで一対一対応する)．指数を $\bm{q}\cdot\bm{R}_{l'} = \bm{q}\cdot\bm{R}_l + \bm{q}\cdot\bm{R}$ と分解して，$\bm{R}_l$ に依存する部分を和の外へ出す:

$$ \text{(右辺)} = -\,e^{i(\bm{q}\cdot\bm{R}_l - \omega t)} \sum_{\kappa'\beta} \left[ \sum_{\bm{R}} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R})\, e^{i\bm{q}\cdot\bm{R}} \right] \frac{e_{\kappa'\beta}}{\sqrt{M_{\kappa'}}}. \label{eq:16-ph-shift} $$

ここが計算の核心である．$\bm{R}_l$ 依存性が共通因子 $e^{i\bm{q}\cdot\bm{R}_l}$ に完全にくくり出された——無限個の胞についての連立方程式が，1つの胞についての方程式に化けたのである．これは並進対称性のおかげであり，上の数学ノートで述べたBlochの定理と同じ機構である．

Step 4:共通因子を払う．\eqref{eq:16-ph-lhs}と\eqref{eq:16-ph-shift}を等置すると，両辺に共通の $e^{i(\bm{q}\cdot\bm{R}_l - \omega t)}$ が消せる．さらに両辺を $\sqrt{M_\kappa}$ で割り，全体の符号を反転すると

$$ \omega^2\, e_{\kappa\alpha} = \sum_{\kappa'\beta} \left[ \frac{1}{\sqrt{M_\kappa M_{\kappa'}}} \sum_{\bm{R}} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R})\, e^{i\bm{q}\cdot\bm{R}} \right] e_{\kappa'\beta}. \label{eq:16-ph-cancel} $$

($1/\sqrt{M_\kappa}$ を角括弧の中に入れ，もともとあった $1/\sqrt{M_{\kappa'}}$ と合わせて $1/\sqrt{M_\kappa M_{\kappa'}}$ とした．)角括弧の中身に名前を付ければ導出は完了する．∎

すなわち，角括弧の中身を動的行列(dynamical matrix)

$$ D_{\alpha\beta}^{\kappa\kappa'}(\bm{q}) = \frac{1}{\sqrt{M_\kappa M_{\kappa'}}} \sum_{\bm{R}} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R})\, e^{i\bm{q}\cdot\bm{R}} \label{eq:16-ph-dyn-def} $$

と定義すると，運動方程式は $3n \times 3n$ 行列の固有値問題

$$ \sum_{\kappa'\beta} D_{\alpha\beta}^{\kappa\kappa'}(\bm{q})\, e_{\kappa'\beta}(\bm{q}) = \omega^2(\bm{q})\, e_{\kappa\alpha}(\bm{q}) \label{eq:16-ph-eig} $$

に帰着する．動的行列は力の定数行列を格子上でFourier変換し，質量で重み付けしたものである．16.4節冒頭の数学ノートに現れた質量重み付きHessian $H_{I\alpha,J\beta}/\sqrt{M_I M_J}$ の，周期系への一般化がまさにこれである．

2つの性質を確認しておく．第一に，$D(\bm{q})$ はエルミートである．式\eqref{eq:16-ph-symm}を周期系の記法に翻訳すると $\Phi_{\kappa\alpha,\kappa'\beta}(\bm{R}) = \Phi_{\kappa'\beta,\kappa\alpha}(-\bm{R})$ であるから

\begin{align} \left[ D_{\alpha\beta}^{\kappa\kappa'}(\bm{q}) \right]^{*} &= \frac{1}{\sqrt{M_\kappa M_{\kappa'}}} \sum_{\bm{R}} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R})\, e^{-i\bm{q}\cdot\bm{R}} \label{eq:16-ph-herm1}\\[2pt] &= \frac{1}{\sqrt{M_{\kappa'} M_{\kappa}}} \sum_{\bm{R}} \Phi_{\kappa\alpha,\, \kappa'\beta}(-\bm{R})\, e^{i\bm{q}\cdot\bm{R}} \label{eq:16-ph-herm2}\\[2pt] &= \frac{1}{\sqrt{M_{\kappa'} M_{\kappa}}} \sum_{\bm{R}} \Phi_{\kappa'\beta,\, \kappa\alpha}(\bm{R})\, e^{i\bm{q}\cdot\bm{R}} = D_{\beta\alpha}^{\kappa'\kappa}(\bm{q}) \label{eq:16-ph-herm3} \end{align}

となる．1行目は $\Phi$ が実数であることから複素共役が指数因子にしか効かないこと，2行目は和の変数を $\bm{R}\to-\bm{R}$ と置き換えたこと(格子ベクトルの集合は符号反転で不変)，3行目は上記の対称性を使ったことによる．したがって固有値 $\omega^2(\bm{q})$ は実数である．第二に，任意の逆格子ベクトル $\bm{G}$ について $e^{i\bm{G}\cdot\bm{R}} = 1$ だから $D(\bm{q}+\bm{G}) = D(\bm{q})$ であり，$\bm{q}$ は第一Brillouin域(第12章)の中だけ考えれば十分である．

(5) フォノン分散と分枝．各 $\bm{q}$ について $3n$ 個の固有値 $\omega_\nu^2(\bm{q})$($\nu = 1,\dots,3n$)が得られる．その平方根 $\omega_\nu(\bm{q})$ がフォノン振動数であり，$\bm{q}$ の関数として描いたものをフォノン分散と呼ぶ．$\nu$ で区別される $3n$ 本の曲線を分枝(branch)という．

このうち3本は必ず $\bm{q}\to\bm{0}$ で $\omega\to0$ になる．これは総和則\eqref{eq:16-ph-asr}の直接の帰結である．$\bm{q}=\bm{0}$ での動的行列は $D_{\alpha\beta}^{\kappa\kappa'}(\bm{0}) = (M_\kappa M_{\kappa'})^{-1/2}\sum_{\bm{R}}\Phi_{\kappa\alpha,\kappa'\beta}(\bm{R})$ である．ここに試しに「質量重み付きの一様並進」$e_{\kappa\alpha} = \sqrt{M_\kappa}\, c_\alpha$($c_\alpha$ は $\kappa$ によらない定数)を代入すると

\begin{align} \sum_{\kappa'\beta} D_{\alpha\beta}^{\kappa\kappa'}(\bm{0})\, \sqrt{M_{\kappa'}}\, c_\beta &= \sum_{\kappa'\beta} \frac{\sqrt{M_{\kappa'}}}{\sqrt{M_\kappa M_{\kappa'}}} \left[\sum_{\bm{R}}\Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R})\right] c_\beta \label{eq:16-ph-ac1}\\[2pt] &= \frac{1}{\sqrt{M_\kappa}} \sum_{\beta} c_\beta \left[ \sum_{\kappa'}\sum_{\bm{R}} \Phi_{\kappa\alpha,\, \kappa'\beta}(\bm{R}) \right] = 0 \label{eq:16-ph-ac2} \end{align}

となる．1行目は代入しただけ，2行目は $\sqrt{M_{\kappa'}}$ を約分して和の順序を整理したもので，最後の等号で角括弧の中が総和則\eqref{eq:16-ph-asr}そのもの($\kappa'$ と $\bm{R}$ の和は「他のすべての原子 $j$ についての和」にほかならない)であることを使った．すなわち $c_\alpha$ の3つの独立な取り方($x,y,z$)に対応して，$\bm{q}=\bm{0}$ で固有値 $\omega^2 = 0$ の固有ベクトルが3本存在する．これらを音響分枝と呼ぶ．物理的には「結晶全体を平行移動しても復元力が働かない」ことを言っているにすぎず，長波長極限で $\omega \simeq v_s |\bm{q}|$(音波)となるためこの名がある．残る $3n-3$ 本は $\bm{q}\to\bm{0}$ でも有限の振動数を持ち，光学分枝と呼ばれる．胞内の原子どうしが逆位相で動くモードであり，イオン結晶では電気双極子を作るので赤外光と結合する——これが「光学」の由来である．単原子単位胞($n=1$)なら光学分枝は存在しない．

(6) 具体例:1次元2原子鎖．抽象論を具体化するため，解析的に解ける最小の例を最後まで計算する．格子定数 $a$ の1次元鎖を考え，単位胞に質量 $M_1$ と $M_2$($M_1 \gt M_2$ とする)の2原子を置く．胞 $l$ の原子1は $x = la$，原子2は $x = la + a/2$ にあり，隣り合う原子はすべて同じ力の定数 $K$ のばねで結ばれているとする．変位は鎖に沿った方向(縦波)だけを考えるので，添字 $\alpha,\beta$ は1種類しかなく，以下では省略する．

(c) 1次元2原子鎖（M1=3M2）のフォノン分散関係．音響枝は原点から線形に立ち上がり，ゾーン境界で最大値に達する．光学枝はq=0で有限の値を取る．両枝の間に振動数ギャップ（黄色の帯）がある． — 図16.4 1次元2原子鎖のフォノン．(a) 単位胞に質量 $M_1 \gt M_2$ の2原子を持ち，隣接原子が力の定数 $K$ のばねで結ばれた鎖．(b) $q\to0$ での2つのモード．音響モードでは全原子が同じ向きに動き(結晶の並進)，光学モードでは2種類の原子が逆向きに動いて重心は静止する．振幅は質量に反比例するので，軽い原子のほうが大きく振れる．(c) 分散関係\eqref{eq:16-ph-chain-sol}($M_1 = 3M_2$)．下の分枝が音響枝で，原点で $\omega \simeq v_s\abs{q}$ と線形に立ち上がる(緑破線)．上の分枝が光学枝で，$q=0$ で $\sqrt{2K/\mu}$ の有限値を取る．両者の間には振動数ギャップ(黄色の帯)があり，その端はゾーン境界での値 $\sqrt{2K/M_1}$ と $\sqrt{2K/M_2}$ で与えられる．

導出:1次元2原子鎖の分散関係

Step 1:調和エネルギーと運動方程式．ばねの伸びは隣接原子の変位差だから，調和エネルギーは

$$ E = E_0 + \frac{K}{2}\sum_{l}\left[ \left(u_{l2} - u_{l1}\right)^2 + \left(u_{l+1,1} - u_{l2}\right)^2 \right] \label{eq:16-ph-chain-e} $$

である(1つ目の括弧が胞内のばね，2つ目が隣の胞へ渡るばねである)．$u_{l1}$ で微分して符号を反転すると，原子1に働く力が得られる．$u_{l1}$ を含む項は $\left(u_{l2}-u_{l1}\right)^2$ と $\left(u_{l1}-u_{l-1,2}\right)^2$ の2つだから

\begin{align} M_1 \ddot{u}_{l1} &= -\frac{\partial E}{\partial u_{l1}} = -\frac{K}{2}\Big[ 2\left(u_{l2}-u_{l1}\right)\cdot(-1) + 2\left(u_{l1}-u_{l-1,2}\right)\cdot(+1) \Big] \label{eq:16-ph-chain-eom1a}\\[2pt] &= K\left( u_{l2} + u_{l-1,2} - 2u_{l1} \right). \label{eq:16-ph-chain-eom1} \end{align}

1行目は連鎖律を機械的に適用しただけであり，括弧内の $(-1)$ と $(+1)$ は各括弧の中身を $u_{l1}$ で微分した値である．同様に原子2については

$$ M_2 \ddot{u}_{l2} = K\left( u_{l+1,1} + u_{l1} - 2u_{l2} \right). \label{eq:16-ph-chain-eom2} $$

Step 2:力の定数を読み取る．式\eqref{eq:16-ph-force-lin}と見比べると(力 $=-\sum\Phi u$ の形に揃えて係数を読む)，ゼロでない力の定数は

$$ \Phi_{11}(0) = 2K,\quad \Phi_{12}(0) = -K,\quad \Phi_{12}(-a) = -K,\quad \Phi_{22}(0) = 2K,\quad \Phi_{21}(0) = -K,\quad \Phi_{21}(+a) = -K \label{eq:16-ph-chain-fc} $$

である．ここで $\Phi_{\kappa\kappa'}(R)$ の $R$ は「相手の胞 $-$ 自分の胞」であり，$\Phi_{12}(-a)$ は「胞 $l$ の原子1が，1つ手前の胞 $l-1$ の原子2から受ける項」を表す．総和則\eqref{eq:16-ph-asr}の確認もしておこう:$\Phi_{11}(0)+\Phi_{12}(0)+\Phi_{12}(-a) = 2K - K - K = 0$ で，確かに成り立っている．

Step 3:動的行列を組み立てる．定義\eqref{eq:16-ph-dyn-def}に\eqref{eq:16-ph-chain-fc}を代入する．$\bm{q}\cdot\bm{R}$ は1次元では $qR$ である:

\begin{align} D_{11}(q) &= \frac{1}{M_1}\,\Phi_{11}(0)\,e^{i q \cdot 0} = \frac{2K}{M_1}, \label{eq:16-ph-chain-D11}\\[2pt] D_{22}(q) &= \frac{1}{M_2}\,\Phi_{22}(0) = \frac{2K}{M_2}, \label{eq:16-ph-chain-D22}\\[2pt] D_{12}(q) &= \frac{1}{\sqrt{M_1M_2}}\left[ \Phi_{12}(0)\,e^{i q\cdot 0} + \Phi_{12}(-a)\,e^{-iqa} \right] = \frac{-K\left( 1 + e^{-iqa} \right)}{\sqrt{M_1M_2}}, \label{eq:16-ph-chain-D12}\\[2pt] D_{21}(q) &= \frac{-K\left( 1 + e^{+iqa} \right)}{\sqrt{M_1M_2}} = \left[ D_{12}(q) \right]^{*}. \label{eq:16-ph-chain-D21} \end{align}

最後の行は，予告どおり $D$ がエルミートであることの確認になっている．

Step 4:固有値方程式を書く．$2\times2$ 行列の固有値は特性方程式 $\det\left[ D(q) - \omega^2 \bm{1}\right] = 0$ から決まる:

$$ \left( \frac{2K}{M_1} - \omega^2 \right)\left( \frac{2K}{M_2} - \omega^2 \right) - D_{12}(q)\, D_{21}(q) = 0. \label{eq:16-ph-chain-det} $$

Step 5:非対角項の絶対値2乗を実数化する．

\begin{align} D_{12} D_{21} = \abs{D_{12}}^2 &= \frac{K^2}{M_1M_2}\left( 1 + e^{-iqa} \right)\left( 1 + e^{+iqa} \right) \label{eq:16-ph-chain-abs1}\\[2pt] &= \frac{K^2}{M_1M_2}\left( 1 + e^{iqa} + e^{-iqa} + 1 \right) = \frac{K^2}{M_1M_2}\left( 2 + 2\cos qa \right) \label{eq:16-ph-chain-abs2}\\[2pt] &= \frac{4K^2}{M_1M_2}\cos^2\!\left( \frac{qa}{2} \right). \label{eq:16-ph-chain-abs3} \end{align}

2行目は展開して $e^{i\theta}+e^{-i\theta} = 2\cos\theta$ を使い，3行目は半角公式 $1 + \cos\theta = 2\cos^2(\theta/2)$ を $\theta = qa$ に適用した．

Step 6:$\omega^2$ の2次方程式に整理する．\eqref{eq:16-ph-chain-det}を展開して\eqref{eq:16-ph-chain-abs3}を代入する:

\begin{align} \omega^4 - 2K\left( \frac{1}{M_1} + \frac{1}{M_2} \right)\omega^2 &+ \frac{4K^2}{M_1M_2} - \frac{4K^2}{M_1M_2}\cos^2\!\left( \frac{qa}{2} \right) = 0 \label{eq:16-ph-chain-quad1}\\[2pt] \omega^4 - 2K\left( \frac{1}{M_1} + \frac{1}{M_2} \right)\omega^2 &+ \frac{4K^2}{M_1M_2}\sin^2\!\left( \frac{qa}{2} \right) = 0 . \label{eq:16-ph-chain-quad} \end{align}

1行目は括弧を展開して $\omega^2$ について整理しただけである．2行目では定数項を $1 - \cos^2 = \sin^2$ でまとめた．$q$ 依存性がこの1つの $\sin^2(qa/2)$ に集約されたことに注意しよう．

Step 7:解の公式を適用する．換算質量 $\mu$ を $\mu^{-1} = M_1^{-1} + M_2^{-1}$，すなわち $\mu = M_1M_2/(M_1+M_2)$ で定義すると，\eqref{eq:16-ph-chain-quad}は $\omega^4 - (2K/\mu)\omega^2 + (4K^2/M_1M_2)\sin^2(qa/2) = 0$ と書ける．$\omega^2$ についての2次方程式として解の公式を使うと

\begin{align} \omega^2 &= \frac{K}{\mu} \pm \sqrt{ \left(\frac{K}{\mu}\right)^2 - \frac{4K^2}{M_1M_2}\sin^2\!\left(\frac{qa}{2}\right) } \label{eq:16-ph-chain-sol1} \end{align}

となる(判別式の $1/4$ を平方根の中に入れて整理した)．平方根の中から $(K/\mu)^2$ をくくり出せば

$$ \omega_{\pm}^2(q) = \frac{K}{\mu}\left[ 1 \pm \sqrt{ 1 - \frac{4\mu^2}{M_1M_2}\sin^2\!\left( \frac{qa}{2} \right) } \right] \label{eq:16-ph-chain-sol} $$

を得る．$\omega_-$ が音響枝，$\omega_+$ が光学枝である．∎

この解析解から4つの極限を読み取る．いずれも1行ずつ確かめる．

(i) 長波長極限の音響枝．$qa \ll 1$ では $\sin(qa/2) \simeq qa/2$ だから，平方根の中の第2項は $\epsilon \equiv \dfrac{4\mu^2}{M_1M_2}\cdot\dfrac{q^2a^2}{4} = \dfrac{\mu^2 q^2 a^2}{M_1M_2}$ となる．$\sqrt{1-\epsilon} \simeq 1 - \epsilon/2$ を使うと

\begin{align} \omega_-^2 &\simeq \frac{K}{\mu}\left[ 1 - \left( 1 - \frac{\epsilon}{2} \right) \right] = \frac{K}{\mu}\cdot\frac{\epsilon}{2} = \frac{K}{\mu}\cdot\frac{\mu^2 q^2 a^2}{2M_1M_2} = \frac{K\mu\, q^2 a^2}{2M_1M_2} \label{eq:16-ph-chain-long1}\\[2pt] &= \frac{K a^2}{2(M_1+M_2)}\, q^2 . \label{eq:16-ph-chain-long} \end{align}

最後の等号では $\mu/(M_1M_2) = 1/(M_1+M_2)$(換算質量の定義から直ちに従う)を使った．したがって

$$ \omega_-(q) = v_s \abs{q}, \qquad v_s = a\sqrt{\frac{K}{2(M_1+M_2)}} \label{eq:16-ph-chain-vs} $$

である．振動数が波数に比例し，比例係数が音速になる——確かに音波である．分母の $M_1+M_2$ は単位胞の全質量であり，長波長では2つの原子が一体となって動くという描像と整合している．

(ii) $q\to0$ の光学枝．同じ展開で複号を $+$ に取ると，$\epsilon\to0$ の最低次で

$$ \omega_+^2(0) = \frac{K}{\mu}\left[ 1 + 1 \right] = \frac{2K}{\mu} = 2K\left( \frac{1}{M_1} + \frac{1}{M_2} \right) \label{eq:16-ph-chain-opt0} $$

である．これは「質量 $\mu$ の粒子がばね定数 $2K$ に繋がれた振動子」の振動数にほかならない．実際，固有ベクトルを調べると $\bm{e} \propto (\sqrt{M_2},\, -\sqrt{M_1})$，すなわち実際の変位は $u_1 \propto \sqrt{M_2/M_1}$，$u_2 \propto -\sqrt{M_1/M_2}$ であり，$M_1u_1 + M_2u_2 = 0$ ——重心が静止したまま2種類の原子が逆位相で振動する(図16.4(b)下段)．一方 $\omega_-$ 側の固有ベクトルは $\bm{e}\propto(\sqrt{M_1},\sqrt{M_2})$，実変位では $u_1 = u_2$ で，これは(5)で一般的に示した一様並進である．

(iii) ゾーン境界($q = \pi/a$)．$\sin^2(qa/2) = \sin^2(\pi/2) = 1$ だから平方根の中は

$$ 1 - \frac{4\mu^2}{M_1M_2} = 1 - \frac{4M_1M_2}{(M_1+M_2)^2} = \frac{(M_1+M_2)^2 - 4M_1M_2}{(M_1+M_2)^2} = \frac{(M_1-M_2)^2}{(M_1+M_2)^2} \label{eq:16-ph-chain-zb1} $$

となる(1つ目の等号で $\mu^2 = M_1^2M_2^2/(M_1+M_2)^2$ を代入，3つ目で $(A+B)^2 - 4AB = (A-B)^2$ を使った)．$M_1 \gt M_2$ としているので平方根は $(M_1-M_2)/(M_1+M_2)$ であり，$K/\mu = K(M_1+M_2)/(M_1M_2)$ と合わせると

\begin{align} \omega_{\pm}^2\!\left( \frac{\pi}{a} \right) &= \frac{K(M_1+M_2)}{M_1M_2}\cdot\frac{(M_1+M_2) \pm (M_1-M_2)}{M_1+M_2} \label{eq:16-ph-chain-zb2}\\[2pt] &= \frac{K}{M_1M_2}\Big[ (M_1+M_2) \pm (M_1-M_2) \Big] \qquad\Longrightarrow\qquad \omega_+^2 = \frac{2K}{M_2},\quad \omega_-^2 = \frac{2K}{M_1}. \label{eq:16-ph-chain-zb} \end{align}

1行目は $1 \pm r$ を通分しただけ，2行目は約分して複号をそれぞれ計算した．物理的な意味は明快である．$q=\pi/a$ では隣り合う胞が逆位相で動くので，片方の副格子が静止し，もう片方だけが「両隣を固定された1原子」として $2K$ のばねで振動する．したがって振動数は自分の質量だけで決まる．$M_1 \gt M_2$ なら重い原子のほうが低い振動数を持つ($\omega_-$)．

(iv) 振動数ギャップ．音響枝の最大値 $\sqrt{2K/M_1}$ と光学枝の最小値 $\sqrt{2K/M_2}$ の間には，どの $q$ でも実現されない振動数の帯がある(図16.4(c)の黄色の帯)．$M_1 \to M_2$ ではこのギャップが閉じる．当然である——2つの原子が同じ質量なら実質的に格子定数 $a/2$ の単原子鎖であり，単位胞を2倍に取って分散を折り返しただけのものになるからである．イオン結晶で赤外吸収が特定の振動数帯に現れるのは，この光学枝の存在による．

(7) 動的安定性の判定．ここまでの道具立ての，実務上おそらく最も重要な使い道がこれである．動的行列はエルミートだから固有値 $\omega^2$ は実数だが，正である保証はどこにもない．ある $\bm{q}$ で $\omega^2 \lt 0$ になったとき，$\omega = \pm i\abs{\omega}$ は純虚数になる．式\eqref{eq:16-ph-ansatz}の時間因子 $e^{-i\omega t}$ は $e^{\pm\abs{\omega}t}$ となり，振動ではなく指数関数的な発散を表す．これを虚振動数あるいはソフトモードと呼ぶ．

虚振動数が出たとき何が起きているのかを，エネルギーの立場から見ておこう．質量で重み付けした変位 $w_{i\alpha} \equiv \sqrt{M_i}\, u_{i\alpha}$ を導入すると，調和エネルギー\eqref{eq:16-ph-harmonic}は

\begin{align} \Delta E = \frac{1}{2}\sum_{i\alpha,j\beta} \Phi_{i\alpha,j\beta}\, u_{i\alpha}u_{j\beta} &= \frac{1}{2}\sum_{i\alpha,j\beta} \Phi_{i\alpha,j\beta}\, \frac{w_{i\alpha}}{\sqrt{M_i}}\frac{w_{j\beta}}{\sqrt{M_j}} \label{eq:16-ph-unstable1}\\[2pt] &= \frac{1}{2}\sum_{i\alpha,j\beta} \tilde{D}_{i\alpha,j\beta}\, w_{i\alpha}w_{j\beta}, \qquad \tilde{D}_{i\alpha,j\beta} \equiv \frac{\Phi_{i\alpha,j\beta}}{\sqrt{M_iM_j}} \label{eq:16-ph-unstable2} \end{align}

と書き直せる(単に $u = w/\sqrt{M}$ を代入して質量因子を $\Phi$ 側に押し付けただけである)．$\tilde{D}$ は実対称行列であり，上で行った平面波型への変換は，この巨大な行列を $\bm{q}$ ごとの $3n\times3n$ ブロックにブロック対角化する操作にほかならない．したがって $\tilde{D}$ の固有値の全体は $\{\omega_\nu^2(\bm{q})\}$ そのものである．いま $\bm{w}$ を固有値 $\omega^2$ の固有ベクトル方向に取れば，式\eqref{eq:16-quadform}と同じ計算で

$$ \Delta E = \frac{1}{2}\,\omega^2\, \abs{\bm{w}}^2 \label{eq:16-ph-unstable} $$

となる．$\omega^2 \lt 0$ ならば $\Delta E \lt 0$ ——その固有ベクトルの方向に原子を動かすとエネルギーが下がる．すなわち構造は極小点ではなく鞍点であり，動的に不安定である．16.1節の分類でいえば，Hessianに負の固有値がある停留点だったということである．3次以上の非調和項が最終的に発散を止めるので，実際にはエネルギー曲線は二重井戸になり，有限の振幅で新しい極小に落ち着く(図16.5)．

(b) ソフトモード方向のエネルギー．実線の真のエネルギーは二重井戸になり，±u0に極小がある．赤の破線はω²<0の調和近似で，下に開く放物線となる． — 図16.5 動的不安定性の判定．(a) ある分枝の $\omega^2(\bm{q})$ がBrillouin域のM点付近で負になっている．赤く塗った領域が $\omega^2 \lt 0$ の範囲で，そこでは $\omega$ が純虚数になる．(b) その固有ベクトルの方向に原子を変位させたときのエネルギー．調和近似(破線)は $\Delta E = \frac{1}{2}\omega^2 u^2 \lt 0$ で下に開く放物線となり，対称構造は極小ではなく停留点にすぎない．非調和項が加わると二重井戸になり，$\pm u_0$ の歪んだ構造が真の極小である．第18章18.9.3項の事例研究では，この判定によって候補構造が篩にかけられ，さらにソフトモードの波数からストライプ周期が予測される．

ソフトモードが $\bm{q}\neq\bm{0}$ に現れた場合，対応する歪みの周期は $2\pi/\abs{\bm{q}}$ であり，もとの単位胞では表現できない．安定な構造を得るには，その周期を含むように単位胞を拡大(スーパーセル化)してから構造最適化をやり直さねばならない．第18章18.9.3項で扱うシリセン薄膜の事例では，まさにM点のソフトモードが見つかり，そこからストライプ構造の周期が予言される．逆に，すべての $\bm{q}$，すべての分枝で $\omega^2 \gt 0$ であることが確認できて初めて，その構造は「動的に安定」と言える．全エネルギーが低いだけでは不十分なのである．

物理的意味:2種類の「安定」

構造の安定性には少なくとも2つの水準がある．第一は力学的(動的)安定性で，「すべてのフォノン振動数が実である」こと，すなわち微小変位に対して必ず復元力が働くことを意味する．これは局所的な性質であり，本項の判定法が答えるのはこれである．第二は熱力学的安定性で，「同じ組成の他のあらゆる構造(および分解生成物)より自由エネルギーが低い」ことを意味する．ダイヤモンドは常温常圧で動的には完全に安定だが，熱力学的にはグラファイトより高エネルギーな準安定相である．第一原理計算で新物質を提案するときは，この2つを別々に確認しなければならない．フォノン計算は前者の必須の検査であり，相図の計算(凸包解析)は後者に対応する．なお，有限温度では非調和効果でソフトモードが安定化することがあり(強誘電体の常誘電相や高温bcc相はその典型である)，0 Kのフォノン計算だけで「存在しない」と断ずるのは早計な場合もある．

(8) 第一原理での計算法．残るのは「$\Phi$ をどうやって求めるか」である．実用されている方法は大きく2つある．

(a) 有限変位法(フローズンフォノン法)．力の定数の定義に戻れば，$F_{j\beta} = -\partial E/\partial R_{j\beta}$ だから

$$ \Phi_{i\alpha,\, j\beta} = \frac{\partial^2 E}{\partial R_{i\alpha}\partial R_{j\beta}} = -\frac{\partial F_{j\beta}}{\partial R_{i\alpha}} \label{eq:16-ph-fd-def} $$

である．すなわち原子を1個だけ少しずらして，全原子に働く力を計算し，それを数値微分すればよい．中心差分を使えば

$$ \Phi_{i\alpha,\, j\beta} \simeq -\frac{F_{j\beta}(+h\,\hat{\bm{e}}_\alpha^{(i)}) - F_{j\beta}(-h\,\hat{\bm{e}}_\alpha^{(i)})}{2h} \label{eq:16-ph-fd} $$

となる($\hat{\bm{e}}_\alpha^{(i)}$ は原子 $i$ を $\alpha$ 方向に動かす単位変位，$h$ は変位量で $0.01$〜$0.02$ Å 程度)．中心差分を選ぶ理由は精度である．$F(\pm h) = F(0) \pm hF' + \frac{h^2}{2}F'' \pm \frac{h^3}{6}F''' + \cdots$ を引き算すると偶数次の項が消え

$$ \frac{F(+h) - F(-h)}{2h} = F'(0) + \frac{h^2}{6}F'''(0) + O(h^4) \label{eq:16-ph-fd-err} $$

となって誤差が $O(h^2)$ に落ちる．しかも消える $F''$ の項は非調和性(エネルギーの3階微分)そのものだから，中心差分は非調和汚染を1次で除去してくれるという副次的な利点も持つ．$h$ を小さくしすぎると力の数値ノイズ(16.4.4項)が $\delta F/h$ として増幅されるので，$h$ には最適値がある．

ここで本質的なのがスーパーセルである．1つの単位胞だけで計算すると，周期境界条件のせいで「原子 $i$ を動かす」ことが「すべての等価な原子を同時に動かす」ことになってしまい，得られるのは $\bm{q}=\bm{0}$ の情報だけである．そこで単位胞を $N_1\times N_2\times N_3$ 倍に拡大したスーパーセルを作り，その中で1原子だけを変位させる．すると，変位した原子から測ってスーパーセルの半分の大きさまでの範囲の力の定数 $\Phi(\bm{R})$ が正しく得られる．つまりスーパーセルの大きさが，取り込める原子間相互作用の到達距離を決める．$\Phi(\bm{R})$ が距離とともに十分減衰していれば打ち切ってよく，そうでなければスーパーセルを大きくするしかない．減衰が遅い代表例が，金属のFriedel振動($\cos(2k_Fr)/r^3$，第8章)に由来する長距離力の定数と，極性絶縁体の双極子–双極子相互作用($1/r^3$)である．後者は $\bm{q}\to\bm{0}$ で動的行列に非解析的な項を生み，縦光学モードと横光学モードの分裂(LO–TO分裂)を引き起こすので，Born有効電荷と誘電テンソルを別途計算して解析的に補ってやる必要がある．なお，スーパーセルの格子ベクトル $\bm{R}_s$ について $e^{i\bm{q}\cdot\bm{R}_s}=1$ となる $\bm{q}$(整合波数)では分散は厳密であり，それ以外の $\bm{q}$ は $\Phi(\bm{R})$ からのFourier補間で得られる．

(b) 密度汎関数摂動論(DFPT)．もう一つの道は，数値微分を一切使わず，変位に対する電子系の応答を摂動論で直接解くことである．原子の変位 $u$ は外部ポテンシャルの摂動 $\Delta v_{\mathrm{ext}}$ を生む．それに対する電子密度の1次応答 $\Delta n$ が求まれば，エネルギーの2階微分すなわち $\Phi$ が得られる．核心は，$\Delta n$ を求めるのに Kohn–Sham 軌道の1次変化 $\Delta\psi_i$ を満たす方程式

$$ \left( \hat{H}_{\mathrm{KS}} - \varepsilon_i \right)\ket{\Delta\psi_i} = -\left( \Delta \hat{V}_{\mathrm{SCF}} - \Delta\varepsilon_i \right)\ket{\psi_i}, \qquad \Delta n(\rr) = 4\,\mathrm{Re}\sum_{i}^{\mathrm{occ}} \psi_i^{*}(\rr)\, \Delta\psi_i(\rr) \label{eq:16-ph-sternheimer} $$

を解くことである(Sternheimer方程式)．ここで摂動ポテンシャルは

$$ \Delta V_{\mathrm{SCF}}(\rr) = \Delta v_{\mathrm{ext}}(\rr) + \int \frac{\Delta n(\rr')}{\abs{\rr - \rr'}}\,\dd\rr' + \left.\frac{\dd v_{xc}}{\dd n}\right|_{n_0} \Delta n(\rr) \label{eq:16-ph-dfpt-vscf} $$

であり，$\Delta n$ が $\Delta V_{\mathrm{SCF}}$ を作り，$\Delta V_{\mathrm{SCF}}$ が $\Delta n$ を作るという自己無撞着な構造をしている．これは第8章で扱った線形応答とまったく同じ精神である．第8章では一様電子ガスに点電荷を入れたときの遮蔽を $\chi$ で記述したが，ここでは非一様な結晶に原子変位という摂動を入れたときの応答を，同じ枠組みで自己無撞着に解いているにすぎない．

DFPTの決定的な利点は，波数 $\bm{q}$ を持つ摂動 $\Delta v_{\mathrm{ext}} \propto e^{i\bm{q}\cdot\rr}$ に対して，$\Delta\psi$ が波数 $\kk+\bm{q}$ のBloch関数になるだけで，すべての計算がもとの単位胞の中で閉じることである．したがってスーパーセルを一切作らずに，任意の $\bm{q}$ の動的行列が直接得られる．また，エネルギーの2階微分を得るのに波動関数の1階微分までしか要らない($2n+1$ 定理)ため，計算コストも1回のSCFと同程度で済む．

両者の長短をまとめる．

実装の負担．有限変位法は「力が出せるコード」さえあれば追加実装なしに使える．したがって，ハイブリッド汎関数・DFT+U・van der Waals補正・meta-GGAなど，摂動論の式を導くのが面倒な汎関数でもそのまま使える．DFPTは汎関数と擬ポテンシャルの2階微分まで実装する必要があり，対応していない汎関数では使えない．
任意の $\bm{q}$．DFPTは任意の $\bm{q}$ を単位胞計算だけで扱える．有限変位法はスーパーセルに整合する $\bm{q}$ しか厳密には扱えず，それ以外はFourier補間になる．整合しない波数に本質的な特徴(コーン異常や非整合なソフトモード)がある場合，有限変位法は見落としうる．
計算コスト．有限変位法は「スーパーセルのサイズ」×「対称性で減らした変位パターン数」で決まり，相互作用が長距離だと急激に高くなる．DFPTは「必要な $\bm{q}$ 点の数」×「摂動の数($3n$)」に比例し，$\bm{q}$ を細かく取りたいときは高くつく(実務では粗い $\bm{q}$ 格子で計算し，実空間の $\Phi(\bm{R})$ に戻して補間する)．
誤差の性質．有限変位法には $h$ の選択に伴う打ち切り誤差と力のノイズの兼ね合いがあるが，DFPTには数値微分がないのでその問題がない．一方，有限変位法は非調和項を計算するための出発点にもなる(変位を大きく取って3次・4次の力の定数を抽出する)．

例:振動数の単位と精度の感覚

フォノン振動数は分野によって単位が異なる．原子単位の角振動数 $\omega$ に対し，波数(分光でよく使う)は $\tilde{\nu} = \omega/(2\pi c)$，周波数は $\nu = \omega/2\pi$，エネルギーは $\hbar\omega$ である．換算の目安は $1$ THz $= 33.36$ cm$^{-1}$ $= 4.136$ meV である．典型的な値として，Siの光学フォノン($\Gamma$ 点)は $15.5$ THz $\simeq 517$ cm$^{-1}$，ダイヤモンドは $40$ THz $\simeq 1333$ cm$^{-1}$，金属のフォノンはおおむね $2$〜$10$ THz に収まる．LDA/GGAで計算した振動数の実験値からのずれは，良い条件で $2$〜$5$ % 程度である．ただしこの誤差は格子定数の誤差に敏感で，たとえば実験格子定数を固定して計算するか，計算で最適化した格子定数を使うかで，振動数が数 % 変わることは珍しくない．ソフトモードの有無を議論するときは，格子定数の取り方・$k$ 点数・カットオフ・スーパーセルサイズをすべて収束させたうえで判定しなければならない．虚振動数は，単に計算が収束していないだけでも簡単に現れるからである．

(9) 熱力学量への応用．フォノンが求まると，格子振動に由来する熱力学量が調和近似の範囲で完全に決まる．各モード $(\bm{q},\nu)$ を独立な調和振動子とみなすと，その分配関数から振動の自由エネルギーは

$$ F_{\mathrm{vib}}(T) = \sum_{\bm{q},\nu}\left[ \frac{\hbar\omega_{\bm{q}\nu}}{2} + k_{\mathrm{B}}T \ln\!\left( 1 - e^{-\hbar\omega_{\bm{q}\nu}/k_{\mathrm{B}}T} \right) \right] \label{eq:16-ph-free} $$

で与えられる．第1項は温度に依存しない零点振動エネルギー

$$ E_{\mathrm{ZP}} = \sum_{\bm{q},\nu} \frac{\hbar\omega_{\bm{q}\nu}}{2} \label{eq:16-ph-zpe} $$

である．これは絶対零度でも消えない量子効果であり，軽元素を含む系では無視できない大きさになる．水素結合系では格子定数や反応障壁を $0.1$ eV 規模で変えることもあり(16.5.6項で活性化エネルギーの零点補正に触れたのはこのためである)，同位体効果の起源でもある．定積比熱は $C_V = -T\,\partial^2 F_{\mathrm{vib}}/\partial T^2$ から

$$ C_V = k_{\mathrm{B}} \sum_{\bm{q},\nu} \left( \frac{\hbar\omega_{\bm{q}\nu}}{k_{\mathrm{B}}T} \right)^2 \frac{e^{\hbar\omega_{\bm{q}\nu}/k_{\mathrm{B}}T}}{\left( e^{\hbar\omega_{\bm{q}\nu}/k_{\mathrm{B}}T} - 1 \right)^2} \label{eq:16-ph-cv} $$

となり，高温では各モードが $k_{\mathrm{B}}$ ずつ寄与してDulong–Petitの法則($3nNk_{\mathrm{B}}$)に，低温では音響枝だけが励起されてDebyeの $T^3$ 則に漸近する．さらに，体積を変えながらフォノンを計算して $F(V,T) = E(V) + F_{\mathrm{vib}}(V,T)$ を各温度で最小化すれば，熱膨張係数が得られる．これを準調和近似(quasi-harmonic approximation)と呼ぶ．厳密に調和的な結晶は熱膨張しない($\omega$ が体積によらないなら $F_{\mathrm{vib}}$ の体積微分がゼロになる)から，熱膨張は「体積依存性を通じて非調和性を最低限だけ取り込んだ」効果だと理解できる．なお，虚振動数が1つでもあると式\eqref{eq:16-ph-free}の対数の中身が意味を失うので，熱力学量の計算に進む前に動的安定性を確認しておかねばならない．

16.5 反応経路探索:NEB法

構造最適化は「谷を下る」問題であった．化学反応の理解に必要なのは，その逆——2つの谷を隔てる峠を越える問題である．反応が起こる速さは峠の高さで決まるから，峠の位置とその高さを求めることは，触媒設計・拡散係数の予測・材料の熱的安定性の評価に直結する．しかし峠は極小ではないので，力をゼロにするだけでは求まらない(力がゼロの点は無数にある)．本節では，現在の標準手法であるNEB法(nudged elastic band method)を，その素朴な原型がなぜ失敗するかという物語とともに構成する．

16.5.1 遷移状態・反応座標・活性化エネルギー

定義:MEP・遷移状態・活性化エネルギー

反応物の極小 $\bm{x}_{\mathrm{R}}$ と生成物の極小 $\bm{x}_{\mathrm{P}}$ を結ぶ経路 $\bm{x}(\xi)$($\xi$ は経路に沿った弧長)のうち，経路上のすべての点で，経路に垂直な方向の力がゼロであるもの——

$$ \left. \bm{F}(\bm{x}(\xi)) \right|_{\perp} = \bm{F} - \left( \bm{F} \cdot \hat{\bm{\tau}} \right) \hat{\bm{\tau}} = \bm{0}, \qquad \hat{\bm{\tau}} = \frac{\dd \bm{x}/\dd\xi}{\abs{\dd \bm{x}/\dd\xi}} \label{eq:16-mep-def} $$

を最小エネルギー経路(MEP)と呼ぶ．直観的には「谷底に沿った道」である．この $\xi$ が反応座標にほかならない．MEPに沿ったエネルギー $E(\xi)$ の最大点が遷移状態(1次の鞍点)であり，

$$ E_a = E(\bm{x}^{\ddagger}) - E(\bm{x}_{\mathrm{R}}) \label{eq:16-ea-def} $$

を活性化エネルギー(順方向)と呼ぶ．逆反応の活性化エネルギーは $E(\bm{x}^{\ddagger}) - E(\bm{x}_{\mathrm{P}})$ であり，その差が反応エネルギー $\Delta E = E(\bm{x}_{\mathrm{P}}) - E(\bm{x}_{\mathrm{R}})$ に等しい．

式\eqref{eq:16-mep-def}が示すように，MEPの条件には「力の垂直成分」という射影が最初から現れている．この射影がNEB法の設計原理そのものになる．射影の道具立てを整理しておこう．

数学ノート:単位ベクトルによる平行・垂直分解

単位ベクトル $\hat{\bm{\tau}}$($\abs{\hat{\bm{\tau}}} = 1$)が与えられたとき，任意のベクトル $\bm{v}$ は

$$ \bm{v} = \bm{v}_{\parallel} + \bm{v}_{\perp}, \qquad \bm{v}_{\parallel} = (\bm{v}\cdot\hat{\bm{\tau}})\,\hat{\bm{\tau}}, \qquad \bm{v}_{\perp} = \bm{v} - (\bm{v}\cdot\hat{\bm{\tau}})\,\hat{\bm{\tau}} \label{eq:16-proj-def} $$

と一意に分解できる．行列で書けば $P_{\parallel} = \hat{\bm{\tau}}\hat{\bm{\tau}}^{\mathsf{T}}$，$P_{\perp} = \bm{1} - \hat{\bm{\tau}}\hat{\bm{\tau}}^{\mathsf{T}}$ である．これらが射影演算子であること，すなわち $P^2 = P$ は直接確かめられる:

$$ P_{\parallel}^2 = \hat{\bm{\tau}}\left( \hat{\bm{\tau}}^{\mathsf{T}}\hat{\bm{\tau}} \right)\hat{\bm{\tau}}^{\mathsf{T}} = \hat{\bm{\tau}} \cdot 1 \cdot \hat{\bm{\tau}}^{\mathsf{T}} = P_{\parallel}, \qquad P_{\perp}^2 = \bm{1} - 2P_{\parallel} + P_{\parallel}^2 = \bm{1} - P_{\parallel} = P_{\perp}. \label{eq:16-proj-idem} $$

また $P_{\parallel} P_{\perp} = P_{\parallel} - P_{\parallel}^2 = 0$ であり，2つの成分は互いに直交する部分空間に属する．$3M$ 次元空間が「経路に沿う1次元」と「経路に垂直な $3M-1$ 次元」に直和分解されるわけで，以下ではこの2つの空間でまったく別の力学を課すことになる．

16.5.2 素朴な弾性バンド法(PEB)と，その2つの病

MEPを数値的に求める最も自然な発想は，反応物と生成物の間に $P-1$ 個の中間構造(像, image)$\bm{R}_1, \dots, \bm{R}_{P-1}$ を並べ，隣り合う像をバネでつないだ「数珠」を作り，その全体のエネルギーを下げることである．両端 $\bm{R}_0 = \bm{x}_{\mathrm{R}}$，$\bm{R}_P = \bm{x}_{\mathrm{P}}$ は固定する．目的関数を

$$ S(\bm{R}_1, \dots, \bm{R}_{P-1}) = \sum_{i=0}^{P} E(\bm{R}_i) + \sum_{i=1}^{P} \frac{k}{2} \abs{\bm{R}_i - \bm{R}_{i-1}}^2 \label{eq:16-peb-obj} $$

と置き，$\partial S/\partial \bm{R}_i = \bm{0}$ を目指す．これを素朴な弾性バンド法(plain elastic band, PEB)と呼ぶ．バネは像が両端に落ち込んで散らばってしまうのを防ぐためのものである．像 $i$ に働く力は

\begin{align} \bm{F}_i^{\mathrm{PEB}} = -\frac{\partial S}{\partial \bm{R}_i} = \underbrace{-\nabla E(\bm{R}_i)}_{\text{真の力}\;\bm{F}_i^{\mathrm{true}}} + \underbrace{k\left( \bm{R}_{i+1} + \bm{R}_{i-1} - 2\bm{R}_i \right)}_{\text{バネ力}\;\bm{F}_i^{\mathrm{spring}}} \label{eq:16-peb-force} \end{align}

である(第2項は，$i$ を含む2つのバネ項 $\frac{k}{2}\abs{\bm{R}_{i+1}-\bm{R}_i}^2$ と $\frac{k}{2}\abs{\bm{R}_i - \bm{R}_{i-1}}^2$ をそれぞれ $\bm{R}_i$ で微分して符号を反転し，足し合わせたものである)．ところがこの素朴な方法は，実際にはMEPに収束しない．病は2つあり，しかも互いに逆方向である．

病(1):コーナーカット．バネ力の垂直成分が悪さをする．式\eqref{eq:16-peb-force}のバネ力は離散的な2階微分だから，像の間隔を $h$ として

$$ k\left( \bm{R}_{i+1} + \bm{R}_{i-1} - 2\bm{R}_i \right) \simeq k h^2 \frac{\dd^2 \bm{R}}{\dd \xi^2} = k h^2\, \frac{1}{\mathcal{R}}\, \hat{\bm{n}} \label{eq:16-peb-curv} $$

と書ける．ここで $\mathcal{R}$ は経路の曲率半径，$\hat{\bm{n}}$ は曲率中心を向く単位法線ベクトルである(曲線論の標準的な公式:弧長で2回微分すると法線方向を向き，大きさは曲率 $1/\mathcal{R}$ になる)．すなわちバネ力は経路を内側へ，まっすぐに引き伸ばそうとする張力として働く．MEPが曲がっている場所(反応の途中で運動の様式が切り替わる場所は必ず曲がっている)では，この張力が真の力の垂直成分と釣り合ってしまい，バンドは谷底より内側の高い斜面を通る．これをコーナーカットと呼ぶ．$k$ を大きくするほど深刻になる．

病(2):像のずり落ち．今度は真の力の平行成分が悪さをする．峠の近くでは経路方向にも大きな傾斜があるから，真の力の平行成分が像を谷底(両端)へ押し流す．バネはこれに抗するが，$k$ が小さいと像は両端に集まり，肝心の峠付近には像がほとんど残らない．峠の高さの分解能が失われるのである．$k$ を小さくするほど深刻になる．

ここに手詰まりがある．$k$ を大きくすれば病(1)が，小さくすれば病(2)が悪化し，両者を同時に治す $k$ は存在しない．しかし2つの病を見比べると，原因が力の「余計な成分」にあることが分かる:病(1)はバネ力の垂直成分，病(2)は真の力の平行成分である．ならば，それらを捨ててしまえばよい．

図16.6 素朴な弾性バンド法(PEB)の2つの病．曲線は谷(反応の峠を通る経路)の等高線，赤い点列がバネで結ばれた像である．(a) バネ定数が大きいと，バネ力の垂直成分がバンドを内側へ引き，谷底より高い経路(コーナーカット)に収束する．(b) バネ定数が小さいと，真の力の平行成分が像を両端の谷へ押し流し，峠付近の分解能が失われる．1つのバネ定数で両方を治すことはできない．

16.5.3 NEB法:力を「つつく」

NEB法のアイデアは驚くほど単純である．真の力からは垂直成分だけを，バネ力からは平行成分だけを取り出して足す:

$$ \bm{F}_i^{\mathrm{NEB}} = \bm{F}_i^{\mathrm{true}\,\perp} + \bm{F}_i^{\mathrm{spring}\,\parallel} = -\left[ \nabla E(\bm{R}_i) - \left( \nabla E(\bm{R}_i)\cdot\hat{\bm{\tau}}_i \right)\hat{\bm{\tau}}_i \right] + k\left( \abs{\bm{R}_{i+1} - \bm{R}_i} - \abs{\bm{R}_i - \bm{R}_{i-1}} \right) \hat{\bm{\tau}}_i . \label{eq:16-neb-force} $$

この力がゼロになるまで像を動かす．式\eqref{eq:16-neb-force}の各部分の役割を確認しよう．

第1項(真の力の垂直成分):像を谷底へ落とす．この項がゼロになった状態が，まさにMEPの定義\eqref{eq:16-mep-def}である．平行成分を捨てたので，像が経路に沿ってずり落ちることはない(病(2)の治療)．
第2項(バネ力の平行成分):像の間隔を均等にする．バネの寄与を最初から「経路方向の1次元問題」に限ったので，経路の形そのものには一切影響しない(病(1)の治療)．しかも，単に $k(\bm{R}_{i+1}+\bm{R}_{i-1}-2\bm{R}_i)$ を射影するのではなく，距離の差 $\abs{\bm{R}_{i+1}-\bm{R}_i} - \abs{\bm{R}_i - \bm{R}_{i-1}}$ を使っている点が実装上重要である．この形なら「前のバネが後ろのバネより伸びていれば前へ，逆なら後ろへ」という間隔の均等化だけを行い，経路の長さを縮めようとする余計な張力を持たない．

「つつく(nudge)」という名前は，余計な成分を落として必要な方向にだけそっと押す，この操作に由来する．ここで一つ，理論的な注意を述べておく．式\eqref{eq:16-neb-force}の力は，射影を含むためにもはや何らかの関数の勾配ではない(非保存力である)．したがって「エネルギーが下がったかどうか」で歩幅を決める直線探索(16.4.4項)は使えず，最適化は力だけを見る方式——減衰付きVerlet(速度が力と逆向きになったら速度をゼロにする「クエンチ」)，FIRE法，あるいは力に対するL-BFGS——で行う．NEB計算が構造最適化より収束しにくいと言われる理由の一つがこれである．

図16.7 NEB法における力の射影．像 $i$ において，真の力 $-\nabla E$(灰色細線)は接線 $\hat{\bm{\tau}}_i$ に平行な成分と垂直な成分に分解される．NEB法は真の力からは垂直成分(赤)だけを，バネ力からは平行成分(緑)だけを採用し，残りを捨てる．前者が像を谷底へ落とし，後者が像の間隔を均等に保つ．

16.5.4 接線ベクトルの取り方——上り坂側を使う

式\eqref{eq:16-neb-force}を実装するには，離散的な像の列から接線 $\hat{\bm{\tau}}_i$ を決めなければならない．最も素朴な定義は中心差分

$$ \hat{\bm{\tau}}_i^{\mathrm{cd}} = \frac{\bm{R}_{i+1} - \bm{R}_{i-1}}{\abs{\bm{R}_{i+1} - \bm{R}_{i-1}}} \label{eq:16-neb-tangent-cd} $$

である．ところがこの定義は，経路に垂直な方向の力が強い(谷が急峻な)場合にバンドに折れ曲がり(kink)を生じさせ，収束しないことが知られている．理由は次のように理解できる．中心差分の接線は $\bm{R}_i$ 自身の情報を含まないため，像 $i$ が経路から少しずれたとき，そのずれを戻す復元力が接線の推定に反映されない．エネルギーが急変する側と緩やかな側の情報が同じ重みで混ざることで，射影の基準がずれ，隣り合う像の間で系統的な食い違い(ジグザグ)が増幅されてしまうのである．

HenkelmanとJónssonが2000年に提案した改良版は，エネルギーが高い側の隣接像を向く差分だけを使う(上流差分, upwind差分):

$$ \bm{\tau}_i = \begin{cases} \bm{\tau}_i^{+} \equiv \bm{R}_{i+1} - \bm{R}_i, & E_{i+1} > E_i > E_{i-1} \;(\text{上り坂}), \\[4pt] \bm{\tau}_i^{-} \equiv \bm{R}_i - \bm{R}_{i-1}, & E_{i+1} < E_i < E_{i-1} \;(\text{下り坂}), \end{cases} \label{eq:16-neb-tangent-up} $$

としたうえで，$\bm{R}_i$ が局所的な極大・極小になっている場合($E_{i+1} > E_i < E_{i-1}$ または $E_{i+1} < E_i > E_{i-1}$)だけは，両側の差分をエネルギー差で重み付けして滑らかにつなぐ:

$$ \bm{\tau}_i = \begin{cases} \bm{\tau}_i^{+}\,\Delta E_i^{\max} + \bm{\tau}_i^{-}\,\Delta E_i^{\min}, & E_{i+1} > E_{i-1},\\[4pt] \bm{\tau}_i^{+}\,\Delta E_i^{\min} + \bm{\tau}_i^{-}\,\Delta E_i^{\max}, & E_{i+1} < E_{i-1}, \end{cases} \qquad \begin{aligned} \Delta E_i^{\max} &= \max\left( \abs{E_{i+1}-E_i},\, \abs{E_{i-1}-E_i} \right),\\ \Delta E_i^{\min} &= \min\left( \abs{E_{i+1}-E_i},\, \abs{E_{i-1}-E_i} \right). \end{aligned} \label{eq:16-neb-tangent-ext} $$

最後に $\hat{\bm{\tau}}_i = \bm{\tau}_i/\abs{\bm{\tau}_i}$ と規格化する．要点は「接線は必ずエネルギーが高い側を向く1本の差分から作る」ことであり，これによって折れ曲がりが消え，少ない像でも安定に収束するようになった．エネルギーが単調に変化する区間では場合分けの上2式のどちらかが必ず当てはまり，極値の像でのみ重み付き平均が使われる，という構造である．

16.5.5 Climbing Image NEB:鞍点を精密に射抜く

NEB法で得られるのは，像の列が描く折れ線としてのMEPである．しかし峠の高さ $E_a$ を知りたいのに，像がちょうど鞍点の上に乗る保証はない．像の間隔が $h$ なら，鞍点付近の放物線的な頂上を折れ線で近似する誤差は $O(h^2)$ であり，像を増やさない限り改善しない．

この問題をエレガントに解決するのがClimbing Image NEB(CI-NEB)である．数回の通常NEB反復のあと，最もエネルギーの高い像 $i_{\max}$ を選び，その像にだけ次の力を課す:

$$ \bm{F}_{i_{\max}}^{\mathrm{CI}} = -\nabla E(\bm{R}_{i_{\max}}) + 2\left( \nabla E(\bm{R}_{i_{\max}}) \cdot \hat{\bm{\tau}}_{i_{\max}} \right) \hat{\bm{\tau}}_{i_{\max}} . \label{eq:16-cineb} $$

すなわちバネ力を完全に切り，真の力の接線成分だけを符号反転する．なぜこれで鞍点に到達するのかを確かめよう．$\nabla E = (\nabla E)_{\parallel} + (\nabla E)_{\perp}$ と分解すると，$2(\nabla E \cdot \hat{\bm{\tau}})\hat{\bm{\tau}} = 2(\nabla E)_{\parallel}$ だから

\begin{align} \bm{F}_{i_{\max}}^{\mathrm{CI}} = -\left[ (\nabla E)_{\parallel} + (\nabla E)_{\perp} \right] + 2 (\nabla E)_{\parallel} = \underbrace{+ (\nabla E)_{\parallel}}_{\text{接線方向は}\,E\,\text{を上る}} \;\underbrace{- (\nabla E)_{\perp}}_{\text{垂直方向は}\,E\,\text{を下る}} \label{eq:16-cineb-decomp} \end{align}

となる．この像は，接線方向にはエネルギーを登り，垂直方向には下る．これはまさに1次の鞍点の特徴づけそのものである．そして力がゼロになる条件は

$$ \bm{F}_{i_{\max}}^{\mathrm{CI}} = \bm{0} \iff (\nabla E)_{\parallel} = \bm{0} \;\text{かつ}\; (\nabla E)_{\perp} = \bm{0} \iff \nabla E(\bm{R}_{i_{\max}}) = \bm{0} \label{eq:16-cineb-zero} $$

である(2つの成分は直交する部分空間に属するから，和がゼロなら各々がゼロ)．つまりCI像が止まる点は真の停留点であり，しかもそこに至る力学から1次の鞍点である．像の個数を増やさずに，遷移状態の座標とエネルギーを最適化と同じ精度($10^{-4}$ Hartree/Bohr など)で決められる．この像はバネから解放されているので間隔の均等性は失われるが，残りの像は通常のNEB力で動き続けるため経路全体は保たれる．

16.5.6 活性化エネルギーから反応速度へ

CI-NEBが与えた $E_a$ は，そのままでは「峠の高さ」という静的な量である．これを実験と比較できる速度に翻訳するのが遷移状態理論(transition state theory, TST)である．要点だけを述べる．反応物の谷の中で系が熱平衡にあり，峠を越えた系は二度と戻らないと仮定すると，単位時間あたりの反応確率(速度定数)は

$$ k(T) = \nu\, \exp\!\left( -\frac{E_a}{k_{\mathrm{B}} T} \right) \label{eq:16-tst} $$

と書ける．指数因子は「峠に到達できるだけのエネルギーを熱ゆらぎから受け取る確率」(Boltzmann因子)，前因子 $\nu$ は「峠に向かって試行する頻度」である．調和近似のもとでは，$\nu$ は反応物の極小における $3M$ 個の基準振動数と，鞍点における $3M-1$ 個の実振動数(虚振動数モードを除く)の比

$$ \nu = \frac{\prod_{j=1}^{3M} \nu_j^{\mathrm{min}}}{\prod_{j=1}^{3M-1} \nu_j^{\ddagger}} \label{eq:16-vineyard} $$

で与えられる(Vineyardの公式)．分子・分母の次元が1つずれているので $\nu$ は振動数の次元を持ち，典型的には $10^{12}$〜$10^{13}$ s$^{-1}$ ——原子の振動周期の逆数程度である．

例:障壁の高さが時間スケールを支配する

$\nu = 10^{13}$ s$^{-1}$，$T = 300$ K($k_{\mathrm{B}}T = 25.9$ meV)として，活性化エネルギーごとの平均待ち時間 $\tau = 1/k$ を計算してみる．

$E_a$	$E_a/k_{\mathrm{B}}T$	$\exp(-E_a/k_{\mathrm{B}}T)$	$\tau = 1/k$
0.2 eV	7.7	$4.4\times10^{-4}$	0.23 ns
0.5 eV	19.3	$4.0\times10^{-9}$	25 μs
1.0 eV	38.7	$1.6\times10^{-17}$	1.7 時間
1.5 eV	58.0	$6.3\times10^{-26}$	$5\times10^{4}$ 年

障壁が $0.5$ eV 増えるごとに時間スケールが $10^{8}$ 倍以上に伸びる．ここから2つの重要な教訓が得られる．第一に，第一原理分子動力学(16.6節)が直接追える時間は長くて数百 ps だから，$0.3$ eV 程度を超える障壁の反応は「ただ走らせるだけ」ではまず起こらない(稀事象問題;16.9節)．第二に，交換相関汎関数に由来する $E_a$ の誤差 $0.1$〜$0.2$ eV は，速度定数の $10^2$〜$10^3$ 倍の誤差に化ける．活性化エネルギーの計算値を報告するときは，この指数関数的な増幅を意識しなければならない．なお実際には，峠と谷での零点振動エネルギーの差($0.1$ eV に達することもある．特に水素が関わる反応)や，トンネル効果による低温での加速も補正として効いてくる．

物理的意味:NEB法が計算しているもの・していないもの

NEB法が与えるのは，0 K のポテンシャルエネルギー面上での最小エネルギー経路と鞍点である．有限温度で系が実際にたどる経路は，エントロピーの寄与によってMEPからずれうる(狭い峠より，少し高くても広い峠のほうが通りやすい)．この「自由エネルギー的な障壁」を扱うには，16.9節のメタダイナミクスのような自由エネルギー計算が必要になる．また，NEB法は反応の始状態と終状態を人間が与えることを前提とする．「どんな生成物ができるか」自体を予測する問題ではない点にも注意したい．

16.6 第一原理分子動力学

ここまでは，原子核を「エネルギーが下がる方向へ動かす」という数値最適化の対象として扱ってきた．これからは原子核を時間発展する古典粒子として扱う．すなわちNewtonの運動方程式

$$ M_I\, \frac{\dd^2 \bm{R}_I}{\dd t^2} = \bm{F}_I(\{\bm{R}_J\}) = -\frac{\partial E(\{\bm{R}_J\})}{\partial \bm{R}_I} \label{eq:16-newton-eom} $$

を数値積分する．ここで $M_I$ は原子核の質量である(原子単位系では電子質量が1なので，水素なら $M_{\mathrm{H}} = 1836$，炭素なら $M_{\mathrm{C}} = 21874$ という大きな数になる．この「$10^3$〜$10^5$ の質量比」こそBorn–Oppenheimer近似の根拠であり，同時に原子核を古典粒子として扱える理由でもある)．ポテンシャル $E$ を経験的な関数形で与えるのが古典分子動力学，各時刻に電子状態を解いて求めるのが第一原理分子動力学である．後者の決定的な利点は，化学結合の生成と切断を，あらかじめ何も仮定せずに追えることにある．

16.6.1 Born–Oppenheimer分子動力学の枠組み

最も素直な実装は，各時間ステップで電子系をきちんと基底状態まで収束させ，そこで得た力で核を動かす方法である．これをBorn–Oppenheimer分子動力学(BOMD)と呼ぶ．手順は次のとおりである．

時刻 $t$ の配置 $\{\bm{R}_I(t)\}$ でKohn–Sham方程式をSCF収束させ，$E$ と $n(\rr)$ を得る(第14章)．
Hellmann–Feynman力とPulay力(16.2節)から $\bm{F}_I(t)$ を計算する．
Verlet法(次項)で $\{\bm{R}_I(t+\Delta t)\}$ と $\{\bm{V}_I(t+\Delta t)\}$ を求める．
$t \to t + \Delta t$ として1に戻る．

計算コストの内訳は圧倒的に1が支配的である．したがってBOMDの実用性は，「1ステップあたり何回のSCF反復で済むか」に懸かっている．前ステップの密度を初期推定に使えば反復数は大きく減るが，あとで述べるようにこの操作は時間反転対称性を壊してエネルギーのドリフトを生むため，注意深い扱いを要する．

図16.8 Born–Oppenheimer分子動力学の計算ループ．1ステップごとに電子系を完全にSCF収束させ，そこで得た力で原子核を古典的に時間発展させる．計算量の大半は①が占める．全エネルギー(核の運動エネルギー+電子系の全エネルギー)が保存しているかどうかが，計算の健全性を測る最も鋭い指標である．

16.6.2 Verlet法の導出

運動方程式\eqref{eq:16-newton-eom}を数値積分する方法は無数にあるが，分子動力学ではほぼ例外なくVerlet法とその変種が使われる．まず導出しよう．

導出:Verlet法と局所誤差 $O(\Delta t^4)$

Step 1:前後に展開する．座標 $\bm{R}(t)$ を時刻 $t$ のまわりでTaylor展開する．加速度を $\bm{a} = \ddot{\bm{R}}$，その時間微分(躍度)を $\bm{b} = \dddot{\bm{R}}$，さらに4階微分を $\bm{c}$ と書くと

\begin{align} \bm{R}(t+\Delta t) &= \bm{R}(t) + \bm{V}(t)\,\Delta t + \frac{1}{2}\bm{a}(t)\,\Delta t^2 + \frac{1}{6}\bm{b}(t)\,\Delta t^3 + \frac{1}{24}\bm{c}(t)\,\Delta t^4 + O(\Delta t^5), \label{eq:16-verlet-fwd}\\[2pt] \bm{R}(t-\Delta t) &= \bm{R}(t) - \bm{V}(t)\,\Delta t + \frac{1}{2}\bm{a}(t)\,\Delta t^2 - \frac{1}{6}\bm{b}(t)\,\Delta t^3 + \frac{1}{24}\bm{c}(t)\,\Delta t^4 + O(\Delta t^5). \label{eq:16-verlet-bwd} \end{align}

2行目は1行目で $\Delta t \to -\Delta t$ と置いただけである．奇数次の項だけ符号が変わることに注意する．

Step 2:和を取る．2式を足すと，奇数次の項——特に $\Delta t^3$ の項——が完全に消える:

\begin{align} \bm{R}(t+\Delta t) + \bm{R}(t-\Delta t) = 2\bm{R}(t) + \bm{a}(t)\,\Delta t^2 + \frac{1}{12}\bm{c}(t)\,\Delta t^4 + O(\Delta t^6). \label{eq:16-verlet-sum} \end{align}

$\bm{V}$ と $\bm{b}$ の項が相殺し，$\frac{1}{24}\bm{c}\Delta t^4$ が2つ足されて $\frac{1}{12}\bm{c}\Delta t^4$ になった．$\bm{R}(t+\Delta t)$ について解けば

$$ \bm{R}(t+\Delta t) = 2\bm{R}(t) - \bm{R}(t-\Delta t) + \bm{a}(t)\,\Delta t^2 + O(\Delta t^4), \qquad \bm{a}_I(t) = \frac{\bm{F}_I(t)}{M_I} \label{eq:16-verlet-pos} $$

を得る．これがVerlet法である．1ステップあたりの誤差(局所打切り誤差)が $O(\Delta t^4)$ であることが重要で，これは「加速度しか使っていない2次の公式」としては1次分だけ得をしている．得をした理由は明快で，前後対称に展開して足したから3次項が消えたのである．

Step 3:速度を取り出す．式\eqref{eq:16-verlet-fwd}から\eqref{eq:16-verlet-bwd}を引くと，今度は偶数次が消えて

\begin{align} \bm{R}(t+\Delta t) - \bm{R}(t-\Delta t) = 2\bm{V}(t)\,\Delta t + \frac{1}{3}\bm{b}(t)\,\Delta t^3 + O(\Delta t^5) \quad\Longrightarrow\quad \bm{V}(t) = \frac{\bm{R}(t+\Delta t) - \bm{R}(t-\Delta t)}{2\Delta t} + O(\Delta t^2). \label{eq:16-verlet-vel} \end{align}

速度は座標より1段精度が落ち，しかも時刻 $t$ の速度を知るのに $t+\Delta t$ の座標が要る(運動エネルギーをその場で評価できない)．この不便さを解消するのが次の速度Verlet法である．

Verlet法は座標の漸化式であって，速度を陽に持たない．運動エネルギー(温度)を毎ステップ評価したい分子動力学では，数学的に等価でありながら速度を同時に更新する速度Verlet法のほうが便利である:

\begin{align} \bm{R}(t+\Delta t) &= \bm{R}(t) + \bm{V}(t)\,\Delta t + \frac{1}{2}\bm{a}(t)\,\Delta t^2, \label{eq:16-vv-pos}\\[2pt] \bm{V}(t+\Delta t) &= \bm{V}(t) + \frac{1}{2}\left[ \bm{a}(t) + \bm{a}(t+\Delta t) \right] \Delta t. \label{eq:16-vv-vel} \end{align}

速度の更新に新旧2つの加速度の平均(台形則)を使うのが要点である．この2式が式\eqref{eq:16-verlet-pos}と厳密に等価であることを確かめておこう．

導出:速度Verlet法とVerlet法の等価性

ステップ番号 $n$ で $\bm{R}_n = \bm{R}(n\Delta t)$ などと略記する．式\eqref{eq:16-vv-pos}を $n$ と $n-1$ について書くと

\begin{align} \bm{R}_{n+1} - \bm{R}_n &= \bm{V}_n \Delta t + \frac{1}{2}\bm{a}_n \Delta t^2, \label{eq:16-vv-eqA}\\[2pt] \bm{R}_n - \bm{R}_{n-1} &= \bm{V}_{n-1} \Delta t + \frac{1}{2}\bm{a}_{n-1} \Delta t^2. \label{eq:16-vv-eqB} \end{align}

式\eqref{eq:16-vv-eqA}から\eqref{eq:16-vv-eqB}を引くと

\begin{align} \bm{R}_{n+1} - 2\bm{R}_n + \bm{R}_{n-1} = \left( \bm{V}_n - \bm{V}_{n-1} \right)\Delta t + \frac{1}{2}\left( \bm{a}_n - \bm{a}_{n-1} \right)\Delta t^2 . \label{eq:16-vv-diff} \end{align}

ここで速度更新式\eqref{eq:16-vv-vel}を1ステップ前に適用すると $\bm{V}_n - \bm{V}_{n-1} = \frac{1}{2}(\bm{a}_{n-1} + \bm{a}_n)\Delta t$ である．これを代入すると

\begin{align} \bm{R}_{n+1} - 2\bm{R}_n + \bm{R}_{n-1} = \frac{1}{2}\left( \bm{a}_{n-1} + \bm{a}_n \right)\Delta t^2 + \frac{1}{2}\left( \bm{a}_n - \bm{a}_{n-1} \right)\Delta t^2 = \bm{a}_n \Delta t^2 . \label{eq:16-vv-equiv} \end{align}

$\bm{a}_{n-1}$ が完全に相殺し，式\eqref{eq:16-verlet-pos}そのものが再現された．すなわち速度Verlet法は，Verlet法の座標列を厳密に再現しつつ，同じ時刻の速度を副産物として与える．

実装上は，式\eqref{eq:16-vv-vel}に $\bm{a}(t+\Delta t)$ が現れることに対応して，次の3段階(キック・ドリフト・キック)に分けて書く:

\begin{align} \text{(i) キック:}\quad & \bm{V}(t+\tfrac{\Delta t}{2}) = \bm{V}(t) + \tfrac{\Delta t}{2}\,\bm{a}(t), \nonumber\\ \text{(ii) ドリフト:}\quad & \bm{R}(t+\Delta t) = \bm{R}(t) + \Delta t\, \bm{V}(t+\tfrac{\Delta t}{2}), \nonumber\\ \text{(iii) キック:}\quad & \bm{V}(t+\Delta t) = \bm{V}(t+\tfrac{\Delta t}{2}) + \tfrac{\Delta t}{2}\,\bm{a}(t+\Delta t). \label{eq:16-kdk} \end{align}

(ii) のあとで新しい配置の力を1回だけ計算し，(iii) に使う．1ステップあたりの力の計算(=SCF計算)はちょうど1回で済む．これが，精度の高いRunge–Kutta法(4次なら1ステップに4回の力評価)ではなくVerlet法が使われる第一の理由である．

16.6.3 なぜVerlet法は長時間安定なのか

局所誤差 $O(\Delta t^4)$ という精度だけなら，他にもっと良い積分法がある．それでも分子動力学がVerlet法を選ぶのは，長時間走らせたときにエネルギーが系統的にドリフトしないという，精度とは別次元の性質を持つからである．その根拠は2つある．

(1) 時間反転対称性．式\eqref{eq:16-verlet-pos}は $\bm{R}(t+\Delta t)$ と $\bm{R}(t-\Delta t)$ を入れ替えても同じ式である．したがって，ある時点で全速度を反転させて計算を続ければ，数値誤差の範囲で元の軌道を正確に逆にたどって戻る．Newton方程式そのものが持つ時間反転対称性を，離散化しても壊していないのである．この対称性があると，誤差は「行きと帰りで打ち消し合う」振動的な性質を持ち，一方向に蓄積しにくい．

(2) シンプレクティック性．より強い性質として，速度Verlet写像は正準変換(シンプレクティック写像)になっている．

数学ノート:シンプレクティック写像

1自由度の相空間 $(q, p)$ を考える(多自由度も同じ)．写像 $(q,p) \mapsto (Q,P)$ のJacobi行列を $M = \partial(Q,P)/\partial(q,p)$ とし，

$$ J = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \qquad \text{に対して} \qquad M^{\mathsf{T}} J M = J \label{eq:16-symplectic-def} $$

が成り立つとき，その写像をシンプレクティックであるという．両辺の行列式を取ると $(\det M)^2 \det J = \det J$，すなわち $\det M = \pm 1$(連続変形で恒等写像につながるので $+1$)——つまりシンプレクティック写像は相空間の体積を保つ(Liouvilleの定理)．ハミルトン方程式による厳密な時間発展はシンプレクティックであり，これを離散化しても保つ積分法をシンプレクティック積分法と呼ぶ．

その重要性は後退誤差解析にある．シンプレクティック積分法で得られる数値解は，元のハミルトニアン $H$ の近似解であると同時に，$H$ をわずかに変形した「影のハミルトニアン」

$$ \tilde{H} = H + \Delta t^2 H_2 + \Delta t^4 H_4 + \cdots \label{eq:16-shadow} $$

の厳密解(指数関数的に良い近似)になっている．$\tilde{H}$ は厳密に保存されるから，$H = \tilde{H} - O(\Delta t^2)$ の誤差は $O(\Delta t^2)$ の範囲で永遠に有界に振動するだけで，時間に比例して増えることがない．非シンプレクティックな高次公式(Runge–Kutta法など)は1ステップの精度は高くても，この保証がないためエネルギーが単調にドリフトする．

導出:速度Verlet法のシンプレクティック性と安定条件

Step 1:3つの写像に分解する．キック・ドリフト・キックの各段階\eqref{eq:16-kdk}を $(q, p)$($p = Mv$)の写像として書くと

$$ \text{キック}: \; (q, p) \mapsto (q,\; p + \tfrac{\Delta t}{2} F(q)), \qquad \text{ドリフト}: \; (q, p) \mapsto (q + \tfrac{\Delta t}{M} p,\; p) \label{eq:16-kdk-maps} $$

である．どちらも「一方の変数だけを，他方の関数だけ平行移動する」ずれ変換である．

Step 2:各写像がシンプレクティックであることを確かめる．キック写像のJacobi行列は，$c \equiv \frac{\Delta t}{2}\frac{\dd F}{\dd q}$ として $M_{\mathrm{K}} = \begin{pmatrix} 1 & 0 \\ c & 1\end{pmatrix}$ である．定義\eqref{eq:16-symplectic-def}の左辺を計算しよう．まず

$$ J M_{\mathrm{K}} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}\begin{pmatrix} 1 & 0 \\ c & 1\end{pmatrix} = \begin{pmatrix} c & 1 \\ -1 & 0 \end{pmatrix}, \qquad M_{\mathrm{K}}^{\mathsf{T}} \left( J M_{\mathrm{K}} \right) = \begin{pmatrix} 1 & c \\ 0 & 1\end{pmatrix}\begin{pmatrix} c & 1 \\ -1 & 0 \end{pmatrix} = \begin{pmatrix} c - c & 1 \\ -1 & 0 \end{pmatrix} = J . \label{eq:16-kick-symp} $$

確かにシンプレクティックである．ドリフト写像も同様に $M_{\mathrm{D}} = \begin{pmatrix} 1 & \Delta t/M \\ 0 & 1\end{pmatrix}$ で，同じ計算で $M_{\mathrm{D}}^{\mathsf{T}} J M_{\mathrm{D}} = J$ が示される．

Step 3:合成もシンプレクティックである．2つの写像を続けて行うとJacobi行列は積 $M = M_2 M_1$ になるから

$$ M^{\mathsf{T}} J M = M_1^{\mathsf{T}} \left( M_2^{\mathsf{T}} J M_2 \right) M_1 = M_1^{\mathsf{T}} J M_1 = J . \label{eq:16-symp-compose} $$

したがってキック・ドリフト・キックの合成である速度Verlet法はシンプレクティックである．

Step 4:調和振動子で安定条件を求める．$F = -M\omega^2 q$ の場合を実際に解く．$\theta \equiv \omega \Delta t$ とし，変数を $(q,\; p/(M\omega))$ に取ると，式\eqref{eq:16-vv-pos}, \eqref{eq:16-vv-vel}は

$$ \begin{pmatrix} q_{n+1} \\ p_{n+1}/(M\omega) \end{pmatrix} = \underbrace{\begin{pmatrix} 1 - \dfrac{\theta^2}{2} & \theta \\[6pt] -\theta\left( 1 - \dfrac{\theta^2}{4} \right) & 1 - \dfrac{\theta^2}{2} \end{pmatrix}}_{\displaystyle \equiv A(\theta)} \begin{pmatrix} q_{n} \\ p_{n}/(M\omega) \end{pmatrix} \label{eq:16-vv-harmonic} $$

という線形写像になる(式\eqref{eq:16-vv-pos}に $a_n = -\omega^2 q_n$ を代入すると1行目，式\eqref{eq:16-vv-vel}に $a_{n+1} = -\omega^2 q_{n+1}$ と1行目の結果を代入して整理すると2行目が出る)．行列式は

$$ \det A = \left( 1 - \frac{\theta^2}{2}\right)^2 + \theta^2 \left( 1 - \frac{\theta^2}{4}\right) = 1 - \theta^2 + \frac{\theta^4}{4} + \theta^2 - \frac{\theta^4}{4} = 1 \label{eq:16-vv-det} $$

で，確かに相空間の面積が保たれている．$2\times2$ 行列の固有値は $\lambda^2 - (\mathrm{Tr}A)\lambda + \det A = 0$ の根だから，$\det A = 1$ のもとで $\abs{\lambda} = 1$(振動解=安定)となる条件は判別式が負，すなわち $\abs{\mathrm{Tr}A} \le 2$ である．$\mathrm{Tr}A = 2 - \theta^2$ だから

$$ \abs{2 - \theta^2} \le 2 \qquad \Longleftrightarrow \qquad \theta = \omega \Delta t \le 2 . \label{eq:16-vv-stability} $$

$\omega\Delta t > 2$ では固有値が実数になり $\abs{\lambda} > 1$，すなわち振幅が指数関数的に発散する．時間刻みは系の最高振動数で決まる——これが次項の実務的な指針の数学的根拠である．

16.6.4 時間刻みの選び方

安定限界は $\omega_{\max}\Delta t < 2$ だが，限界ぎりぎりではエネルギー誤差が大きい(影のハミルトニアンの $\Delta t^2$ 補正が効きすぎる)．実用的な目安は，最高振動モードの周期 $T_{\min} = 2\pi/\omega_{\max}$ の $1/20$ 程度である．このとき $\omega\Delta t = 2\pi/20 = 0.31$ で安定限界から十分離れており，相対的なエネルギー誤差も $1\%$ 程度に収まる．

例:時間刻みの実際の見積もり

(a) 水素を含む系．C–H伸縮振動の波数は $\tilde{\nu} \simeq 3000$ cm$^{-1}$ である．振動数は $\nu = c\tilde{\nu} = (3.0\times10^{10}\ \mathrm{cm/s})\times(3000\ \mathrm{cm^{-1}}) = 9.0\times10^{13}$ Hz，周期は $T = 1/\nu = 11$ fs．したがって

$$\Delta t \simeq \frac{11\ \mathrm{fs}}{20} \simeq 0.5\ \mathrm{fs} \;(= 21\ \text{原子単位})$$

となる．第一原理MDで $\Delta t = 0.5$ fs が「水素を含む系の標準」とされるのはこの見積もりによる．

(b) 水素を含まない系．重い元素だけなら最高振動数は $1000$〜$1500$ cm$^{-1}$ 程度に下がり，$T \simeq 22$〜$33$ fs，$\Delta t = 1$〜$1.5$ fs が使える．

(c) 重水素置換．$\omega \propto 1/\sqrt{\mu}$($\mu$ は換算質量)だから，H を D に置き換えると最高振動数は $1/\sqrt{2}$ 倍になり，$\Delta t$ を $\sqrt{2}$ 倍にできる．同位体効果そのものを調べたいのでなければ，これは有効な計算の節約手段である．

(d) 到達できる時間．$\Delta t = 0.5$ fs，1ステップのSCF計算が10秒の系なら，1日で $8.6\times10^3$ ステップ $= 4.3$ ps しか進めない．第一原理MDの現実的な射程が $10$ ps〜$1$ ns 程度である，というのはこの計算から来ている．16.5節の表と見比べれば，$0.3$ eV を超える障壁の反応が「自然には起こらない」ことが分かるだろう．

16.6.5 エネルギー保存というモニタ

断熱的(熱浴なし)なBOMDでは，保存量

$$ E_{\mathrm{cons}} = \sum_I \frac{1}{2} M_I \abs{\bm{V}_I}^2 + E(\{\bm{R}_I\}) \label{eq:16-econs} $$

が時間によらず一定でなければならない．これは，計算の健全性を判定する最も鋭く，最も安価な指標である．$E_{\mathrm{cons}}$ が単調にドリフトするなら，原因は必ず次のいずれかである．

時間刻みが大きすぎる:$\Delta t$ を半分にしてドリフトが $1/4$ になれば($O(\Delta t^2)$)，これが原因である．
SCFの収束が緩い:未収束の密度から作った力は，真のエネルギー面の勾配ではない．すなわち非保存力であり，1ステップごとに正味の仕事をして系を暖める(あるいは冷ます)．構造最適化なら十分な収束条件($10^{-6}$ 程度)でも，MDでは $10^{-8}$〜$10^{-10}$ を要求することがあるのはこのためである．
力とエネルギーが整合していない:Pulay項の実装漏れ(16.2節)，実空間グリッドのeggbox誤差など．この場合，$\Delta t$ を小さくしてもドリフトは消えない．原因を切り分ける良い診断法である．
密度の外挿が時間反転対称性を壊している:SCFの初期推定に過去数ステップの密度を外挿して使うと反復数は減るが，「過去だけを見る」操作は時間反転対称性を破る．その結果，微小だが系統的なエネルギー散逸が生じる．これを避けるために，過去の情報を時間反転対称な形で使う密度外挿法(時間反転対称BOMD)が開発されている．

物理的意味:初期条件の与え方と等分配

目標温度 $T$ のシミュレーションを始めるとき，初速度をどう与えるべきか．素朴に「温度 $T$ に対応するMaxwell–Boltzmann分布」から引くと，走らせた直後に温度は$T/2$ 付近まで落ちる．理由はエネルギー等分配則である．調和的な系では，平衡状態で全エネルギーは運動エネルギーとポテンシャルエネルギーに等分される．極小構造(ポテンシャルエネルギーが最小)から出発すれば，最初に与えた運動エネルギーの半分がポテンシャルに移るからである．したがって，平衡化を急ぐには初期温度を $2T$ に取る．あわせて，全運動量(および孤立系なら全角運動量)をゼロにしておくこと——重心の並進運動は「温度」に数えたくない自由度であり，放っておくと系統的に温度を過大評価する．温度の定義に使う自由度数を $3N$ ではなく $3N-3$ とするのも同じ理由である．

16.7 温度制御:Nosé–Hoover法

前節のBOMDが生成するのは，エネルギー・体積・粒子数が一定のミクロカノニカル集団($NVE$)である．しかし実験は普通，温度が制御された環境で行われる．温度一定のカノニカル集団($NVT$)を生成するには，系を「熱浴」に接続しなければならない．ここで問われるのは，熱浴を決定論的な運動方程式としてどう書くかである．本節では，Noséが1984年に導入しHooverが1985年に整理した拡張系の方法を，拡張Lagrangianから運動方程式まで完全に導出する．

数学ノート:カノニカル分布とエルゴード仮説

温度 $T$ の熱浴に接した系が配位 $(\bm{R}, \bm{p})$ を取る確率密度は

$$ \rho(\bm{R},\bm{p}) = \frac{1}{Z}\, e^{-\beta \mathcal{H}(\bm{R},\bm{p})}, \qquad \mathcal{H} = \sum_I \frac{\abs{\bm{p}_I}^2}{2M_I} + U(\bm{R}), \qquad \beta = \frac{1}{k_{\mathrm{B}}T} \label{eq:16-canonical} $$

である(カノニカル分布)．運動量部分はGauss分布に分解でき，各運動量成分は $\braket{p_\alpha^2/2M} = k_{\mathrm{B}}T/2$ を満たす(エネルギー等分配則)．したがって自由度数を $N_f$ として，瞬間温度を

$$ T_{\mathrm{inst}}(t) \equiv \frac{1}{N_f k_{\mathrm{B}}} \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} \label{eq:16-Tinst} $$

と定義すれば，その長時間平均が熱力学的な温度に一致するはずである(分母の $N_f k_{\mathrm{B}}$ は，$\sum_I \abs{\bm{p}_I}^2/M_I = 2\times$運動エネルギー $= N_f k_{\mathrm{B}}T$ となるように選んである)．

ここで「分布を再現する」と「時間平均が集団平均に等しい」を結ぶのがエルゴード仮説である．分子動力学は本質的に時間平均を計算する道具なので，熱浴付きの運動方程式には (i) 正しい分布を不変分布として持つこと と (ii) その分布を実際に埋め尽くすこと(エルゴード性) の2つが要求される．この2つは別物であり，Nosé–Hoover法では (i) は証明できるが (ii) は自明ではない——後述する非エルゴード性の問題はここに起因する．

もう一つ，以下の証明で使う道具を用意しておく．関数 $f(s)$ が区間内でただ一つの零点 $s_0$ を持ち $f'(s_0) \neq 0$ のとき，デルタ関数の合成則

$$ \delta\left( f(s) \right) = \frac{\delta(s - s_0)}{\abs{f'(s_0)}} \label{eq:16-delta-comp} $$

が成り立つ($s_0$ の近傍で $f(s) \simeq f'(s_0)(s-s_0)$ と線形化し，デルタ関数のスケーリング則 $\delta(ax) = \delta(x)/\abs{a}$ を使えばよい)．

16.7.1 Noséの拡張系という発想

熱浴を表現する素朴な方法は「速度をときどき目標温度に合わせて一律にスケールする」ことだが，これは運動方程式ではないので時間反転対称性もなく，生成される分布もカノニカル分布ではない．Noséのアイデアは根本的に違う．熱浴を，系に付け加えたたった1つの新しい力学変数 $s$ として扱い，拡張された系全体をHamilton力学で走らせるのである．そして，拡張系のミクロカノニカル集団を物理変数へ射影すると，ちょうどカノニカル分布が現れるように $s$ の役割を設計する．

$s$ に持たせる役割は時間の刻みを伸縮させることである．仮想時間 $t'$ と実時間 $t$ を

$$ \dd t = \frac{\dd t'}{s} \label{eq:16-nose-time} $$

で結ぶ．$s > 1$ なら実時間はゆっくり進み(実速度は小さく=冷たく)，$s < 1$ なら速く進む(熱く)．この $s$ に運動エネルギー $\frac{Q}{2}(\dd s/\dd t')^2$($Q$ は熱浴の「質量」，時間の2乗×エネルギーの次元)と，対数型のポテンシャル $g k_{\mathrm{B}}T_d \ln s$ を与えたLagrangianが

$$ \mathcal{L}_{\mathrm{N}} = \sum_I \frac{M_I}{2}\, s^2 \left( \frac{\dd \bm{R}_I}{\dd t'} \right)^{\!2} - U(\{\bm{R}_I\}) + \frac{Q}{2}\left( \frac{\dd s}{\dd t'} \right)^{\!2} - g\, k_{\mathrm{B}} T_d \ln s \label{eq:16-nose-lag} $$

である．$T_d$ は目標温度，$g$ は後で決める定数(16.4節の勾配ベクトル $\bm{g}$ とは無関係な，ただのスカラーである)．第1項に $s^2$ が掛かっているのは，実速度が $\dd\bm{R}/\dd t = s\, \dd\bm{R}/\dd t'$ だから，これが実際の運動エネルギー $\frac{M}{2}(\dd\bm{R}/\dd t)^2$ を表すためである．対数ポテンシャル $\ln s$ という一見奇妙な選択の理由は，次項の分布の計算で明らかになる．

16.7.2 運動方程式の導出

導出:仮想時間のEuler–Lagrange方程式から実時間のHoover方程式へ

Step 1:仮想時間でEuler–Lagrange方程式を書く．以下 $' = \dd/\dd t'$ と略記する．$\bm{R}_I$ については

$$ \frac{\partial \mathcal{L}_{\mathrm{N}}}{\partial \bm{R}_I'} = M_I s^2 \bm{R}_I', \qquad \frac{\partial \mathcal{L}_{\mathrm{N}}}{\partial \bm{R}_I} = -\frac{\partial U}{\partial \bm{R}_I} = \bm{F}_I \quad\Longrightarrow\quad \frac{\dd}{\dd t'}\left( M_I s^2 \bm{R}_I' \right) = \bm{F}_I . \label{eq:16-nose-elR} $$

$s$ については，$\partial \mathcal{L}_{\mathrm{N}}/\partial s' = Q s'$ と

$$ \frac{\partial \mathcal{L}_{\mathrm{N}}}{\partial s} = \sum_I M_I s \abs{\bm{R}_I'}^2 - \frac{g k_{\mathrm{B}}T_d}{s} \quad\Longrightarrow\quad \frac{\dd}{\dd t'}\left( Q s' \right) = \sum_I M_I s \abs{\bm{R}_I'}^2 - \frac{g k_{\mathrm{B}}T_d}{s} \label{eq:16-nose-els} $$

が得られる(第1項は $s^2$ を $s$ で微分して $2s$，係数 $M_I/2$ と合わせて $M_I s$;第2項は $\dd(\ln s)/\dd s = 1/s$)．

Step 2:実時間の変数を定義する．式\eqref{eq:16-nose-time}より $\dd/\dd t' = (1/s)\,\dd/\dd t$ である．実速度と実運動量，そして摩擦係数を

$$ \bm{V}_I = \frac{\dd \bm{R}_I}{\dd t} = s\, \bm{R}_I', \qquad \bm{p}_I = M_I \bm{V}_I = M_I s\, \bm{R}_I', \qquad \zeta \equiv \frac{\dd \ln s}{\dd t} = \frac{1}{s}\frac{\dd s}{\dd t} \label{eq:16-nose-realvars} $$

と定義する．ここで重要な関係を1つ確認しておく:

$$ s' = \frac{\dd s}{\dd t'} = \frac{1}{s}\frac{\dd s}{\dd t} = \zeta . \label{eq:16-nose-sprime} $$

つまり$s$ の仮想時間微分が，そのまま摩擦係数 $\zeta$ である．

Step 3:座標の方程式を実時間に書き直す．式\eqref{eq:16-nose-realvars}より $M_I s^2 \bm{R}_I' = s \cdot (M_I s \bm{R}_I') = s\,\bm{p}_I$ である．これを式\eqref{eq:16-nose-elR}に代入し，$\dd/\dd t' = (1/s)\dd/\dd t$ を使うと

\begin{align} \frac{1}{s}\frac{\dd}{\dd t}\left( s\, \bm{p}_I \right) = \bm{F}_I \quad\Longrightarrow\quad \frac{1}{s}\left[ s\, \frac{\dd \bm{p}_I}{\dd t} + \bm{p}_I \frac{\dd s}{\dd t} \right] = \bm{F}_I \quad\Longrightarrow\quad \frac{\dd \bm{p}_I}{\dd t} + \bm{p}_I \underbrace{\frac{1}{s}\frac{\dd s}{\dd t}}_{= \zeta} = \bm{F}_I . \label{eq:16-nose-Rreal} \end{align}

2つ目の変形で積の微分法則を使い，3つ目で両辺を $s$ で割った．整理して

$$ \frac{\dd \bm{p}_I}{\dd t} = \bm{F}_I - \zeta\, \bm{p}_I . \label{eq:16-nh-p} $$

Step 4:熱浴の方程式を実時間に書き直す．式\eqref{eq:16-nose-els}の左辺は，式\eqref{eq:16-nose-sprime}を使って

$$ \frac{\dd}{\dd t'}\left( Q \zeta \right) = \frac{1}{s}\, Q\, \frac{\dd \zeta}{\dd t} \label{eq:16-nose-lhs} $$

である．右辺は $\bm{R}_I' = \bm{p}_I/(M_I s)$(式\eqref{eq:16-nose-realvars})より

$$ \sum_I M_I s\, \frac{\abs{\bm{p}_I}^2}{M_I^2 s^2} - \frac{g k_{\mathrm{B}}T_d}{s} = \frac{1}{s}\left[ \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} - g k_{\mathrm{B}}T_d \right] \label{eq:16-nose-rhs} $$

となる．両辺に共通の因子 $1/s$ が出たので，これを払うと

$$ \frac{\dd \zeta}{\dd t} = \frac{1}{Q}\left[ \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} - g k_{\mathrm{B}}T_d \right] . \label{eq:16-nh-zeta} $$

$s$ そのものが方程式から消え，$\zeta$ だけの閉じた形になったことに注意してほしい．

得られた式\eqref{eq:16-nh-p}, \eqref{eq:16-nh-zeta}がNosé–Hoover方程式である．まとめて書き，瞬間温度\eqref{eq:16-Tinst}を使って表すと

$$ \frac{\dd \bm{R}_I}{\dd t} = \frac{\bm{p}_I}{M_I}, \qquad \frac{\dd \bm{p}_I}{\dd t} = \bm{F}_I - \zeta\, \bm{p}_I, \qquad \frac{\dd \zeta}{\dd t} = \frac{N_f k_{\mathrm{B}}}{Q}\left[ T_{\mathrm{inst}}(t) - T_d \right] \quad (g = N_f) \label{eq:16-nh-eom} $$

物理的な読み方は明快である．$\zeta$ は符号を変えられる摩擦係数であり，$-\zeta\bm{p}_I$ という項で運動量を減衰(または増幅)させる．式\eqref{eq:16-nh-eom}の第3式によれば，系が目標より熱ければ($T_{\mathrm{inst}} > T_d$)$\zeta$ は増大してブレーキがかかり，冷たければ $\zeta$ は減少して(負になり)加速される．しかも $\zeta$ 自身が慣性($Q$)を持つため，応答は瞬時ではなく振動的である——これが「熱浴も1つの力学自由度である」ということの意味である．

図16.9 Nosé–Hoover熱浴の働き．(a) 瞬間温度を測って摩擦係数 $\zeta$ を駆動し，それが運動方程式に戻る負のフィードバック構造．$\zeta$ 自身が慣性(熱浴質量 $Q$)を持つため応答は振動的になる．(b) 目標温度からずれた初期条件から始めたときの温度の時間変化．単純な速度スケーリングと違い，行き過ぎと戻りを繰り返しながら緩和する．この振動の周期が $\sqrt{Q}$ に比例する．

16.7.3 保存量とカノニカル分布の証明

Nosé–Hoover方程式\eqref{eq:16-nh-eom}はHamilton方程式の形をしていない(摩擦項があるので相空間の体積を保たない)．それでも保存量は存在する．これがモニタとして極めて有用である．

導出:Nosé–Hoover系の保存量

次の量が時間によらず一定であることを示す:

$$ \mathcal{H}_{\mathrm{NH}} = \sum_I \frac{\abs{\bm{p}_I}^2}{2 M_I} + U(\{\bm{R}_I\}) + \frac{Q}{2}\zeta^2 + g\, k_{\mathrm{B}}T_d \ln s, \qquad \frac{\dd \ln s}{\dd t} = \zeta . \label{eq:16-nh-cons} $$

時間微分を素直に取る．$\dd\bm{R}_I/\dd t = \bm{p}_I/M_I$，$\partial U/\partial \bm{R}_I = -\bm{F}_I$ を使うと

\begin{align} \frac{\dd \mathcal{H}_{\mathrm{NH}}}{\dd t} &= \sum_I \frac{\bm{p}_I}{M_I}\cdot \frac{\dd \bm{p}_I}{\dd t} + \sum_I \frac{\partial U}{\partial \bm{R}_I}\cdot \frac{\dd \bm{R}_I}{\dd t} + Q \zeta \frac{\dd \zeta}{\dd t} + g k_{\mathrm{B}}T_d\, \frac{\dd \ln s}{\dd t} \nonumber\\[2pt] &= \sum_I \frac{\bm{p}_I}{M_I}\cdot\left( \bm{F}_I - \zeta \bm{p}_I \right) - \sum_I \bm{F}_I \cdot \frac{\bm{p}_I}{M_I} + Q\zeta \cdot \frac{1}{Q}\left[ \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} - g k_{\mathrm{B}}T_d \right] + g k_{\mathrm{B}}T_d\, \zeta . \label{eq:16-nh-cons-deriv} \end{align}

2行目では運動方程式\eqref{eq:16-nh-p}, \eqref{eq:16-nh-zeta}をすべて代入した．これを項ごとに整理する．第1項の $\bm{F}_I$ の寄与と第2項が打ち消し合う．残るのは

$$ -\zeta \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} + \zeta \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} - \zeta\, g k_{\mathrm{B}}T_d + \zeta\, g k_{\mathrm{B}}T_d = 0 \label{eq:16-nh-cons-zero} $$

である(第1の和と第2の和，$g k_{\mathrm{B}}T_d$ の2項がそれぞれ相殺する)．よって $\mathcal{H}_{\mathrm{NH}}$ は保存する．

$\mathcal{H}_{\mathrm{NH}}$ はエネルギーではない(熱浴が系に出し入れした熱を $\frac{Q}{2}\zeta^2 + g k_{\mathrm{B}}T_d\ln s$ が帳簿として記録している)が，数値積分の健全性を測る指標としては16.6節の $E_{\mathrm{cons}}$ とまったく同じ役割を果たす．$NVT$ 計算でもこの量のドリフトを必ず監視すべきである．

次に，この力学が本当にカノニカル分布を生むことを示す．ここで対数ポテンシャルの意味と $g$ の値が同時に決まる．

導出:拡張系のミクロカノニカル集団がカノニカル分布を与えること

Step 1:拡張系のHamiltonianを作る．Lagrangian\eqref{eq:16-nose-lag}から正準運動量を求めると，$\bm{p}_I' = \partial\mathcal{L}_{\mathrm{N}}/\partial\bm{R}_I' = M_I s^2 \bm{R}_I'$，$p_s = Q s'$ である．Legendre変換 $\mathcal{H}_{\mathrm{N}} = \sum \bm{p}_I'\cdot\bm{R}_I' + p_s s' - \mathcal{L}_{\mathrm{N}}$ を実行すると

$$ \mathcal{H}_{\mathrm{N}} = \sum_I \frac{\abs{\bm{p}_I'}^2}{2 M_I s^2} + U(\bm{R}) + \frac{p_s^2}{2Q} + g k_{\mathrm{B}}T_d \ln s \label{eq:16-nose-ham} $$

となる．$\bm{p}_I' = s\,\bm{p}_I$ という関係(式\eqref{eq:16-nose-realvars}と比較せよ)から，第1項は実運動量では $\sum \abs{\bm{p}_I}^2/2M_I$ である．

Step 2:拡張系のミクロカノニカル分配関数を書く．拡張系はHamilton力学に従うから，そのエネルギー $\mathcal{H}_{\mathrm{N}} = E$ 一定面上の一様分布(ミクロカノニカル分布)を仮定してよい:

$$ Z_{\mathrm{ext}} = \int \dd p_s \int \dd s \int \dd^{3N}p' \int \dd^{3N}R\; \delta\!\left( \mathcal{H}_{\mathrm{N}} - E \right). \label{eq:16-nose-Z} $$

Step 3:実運動量に変数変換する．$\bm{p}' = s\,\bm{p}$ とすると，$3N$ 個の成分それぞれに $s$ が掛かるからJacobianは $s^{3N}$，すなわち $\dd^{3N}p' = s^{3N}\, \dd^{3N}p$ である．物理系のHamiltonianを $\mathcal{H}(\bm{p},\bm{R}) = \sum\abs{\bm{p}_I}^2/2M_I + U$ と書くと

$$ Z_{\mathrm{ext}} = \int \dd p_s \dd^{3N}p\, \dd^{3N}R \int \dd s\; s^{3N}\, \delta\!\left( \mathcal{H} + \frac{p_s^2}{2Q} + g k_{\mathrm{B}}T_d \ln s - E \right). \label{eq:16-nose-Z2} $$

Step 4:$s$ 積分をデルタ関数で実行する．括弧の中を $f(s)$ と見ると，$f$ は $s$ について単調増加($\ln s$ の項による)で，零点は

$$ s_0 = \exp\!\left[ \frac{E - \mathcal{H} - p_s^2/2Q}{g k_{\mathrm{B}}T_d} \right], \qquad f'(s) = \frac{g k_{\mathrm{B}}T_d}{s} \label{eq:16-nose-s0} $$

である．合成則\eqref{eq:16-delta-comp}を使うと $\delta(f(s)) = s_0\,\delta(s-s_0)/(g k_{\mathrm{B}}T_d)$ だから，$s$ 積分は

$$ \int \dd s\; s^{3N} \delta(f(s)) = \frac{s_0^{3N+1}}{g k_{\mathrm{B}}T_d} = \frac{1}{g k_{\mathrm{B}}T_d} \exp\!\left[ \frac{(3N+1)\left( E - \mathcal{H} - p_s^2/2Q \right)}{g k_{\mathrm{B}}T_d} \right] \label{eq:16-nose-sint} $$

となる($s_0^{3N}\cdot s_0 = s_0^{3N+1}$ に注意)．

Step 5:$g$ を選ぶ．指数の中で $\mathcal{H}$ に掛かる係数が $-1/(k_{\mathrm{B}}T_d)$ になればカノニカル分布である．すなわち $(3N+1)/g = 1$，つまり

$$ g = 3N + 1 \qquad \Longrightarrow \qquad Z_{\mathrm{ext}} \propto \int \dd^{3N}p\, \dd^{3N}R\; e^{-\mathcal{H}/k_{\mathrm{B}}T_d} \times \int \dd p_s\, e^{-p_s^2/(2Q k_{\mathrm{B}}T_d)} . \label{eq:16-nose-canonical} $$

熱浴変数 $p_s$ の積分は $\mathcal{H}$ と分離した定数因子になる．したがって拡張系のミクロカノニカル集団を物理変数に射影すると，厳密にカノニカル分布が得られる．$\ln s$ というポテンシャルを選んだ理由がここで分かる:デルタ関数の零点が $\exp[\cdots]$ という指数関数の形になり，Boltzmann因子をちょうど作り出すからである．

Step 6:実時間に直すと $g = 3N$．上の計算は仮想時間 $t'$ での一様サンプリング(=拡張系の相空間の一様分布)に対応する．ところがHoover方程式\eqref{eq:16-nh-eom}は実時間 $t$ で書かれており，$\dd t = \dd t'/s$ だから，実時間での時間平均は仮想時間平均に対して重み $1/s$ が掛かる:

$$ \braket{A}_{t} = \frac{\int \dd t\, A}{\int \dd t} = \frac{\braket{A/s}_{t'}}{\braket{1/s}_{t'}} . \label{eq:16-nose-timeavg} $$

この $1/s$ の因子は，Step 4の $s_0^{3N+1}$ を $s_0^{3N}$ に変える．同じ議論を繰り返せば $(3N)/g = 1$，すなわち $g = 3N$ が正しい選択となる．より一般には，拘束(重心の固定など)を考慮した実際の自由度数 $N_f$ を使って $g = N_f$ とする．式\eqref{eq:16-nh-eom}をこの形で書いたのはそのためである．

16.7.4 熱浴質量の選び方と非エルゴード性

熱浴質量 $Q$．$Q$ は分布そのものには影響しない(上の証明に $Q$ は本質的に現れなかった)が，収束の速さを決める．線形化した解析から，熱浴の振動周期を $\tau_{\mathrm{th}}$ とすると

$$ Q = N_f\, k_{\mathrm{B}} T_d\, \tau_{\mathrm{th}}^2 \label{eq:16-nh-Q} $$

と選ぶのが標準である．$\tau_{\mathrm{th}}$ は系の代表的な振動周期と同程度(数十 fs〜数百 fs)に取る．$Q$ が大きすぎると熱浴の応答が遅く，実質的にミクロカノニカルなまま平衡化しない．$Q$ が小さすぎると熱浴が高周波で振動し，系の振動と共鳴して人工的な振る舞いを生む．

非エルゴード性という落とし穴．数学ノートで注意したとおり，正しい不変分布を持つことと，それを埋め尽くすことは別問題である．Nosé–Hoover法の有名な病理が，1個の調和振動子で現れる．この系(物理自由度2，熱浴自由度2の合計4次元相空間)にNosé–Hoover熱浴を付けても，軌道の大部分は規則的なトーラス上に閉じ込められ，カノニカル分布を再現しない．原因は自由度が少なすぎて力学が十分にカオス的でないことである．実際の分子・固体は自由度が多く非調和なので通常は問題ないが，剛性の高い分子や低温の固体，あるいは少数自由度の模型系では現実に起こりうる．

この問題の標準的な解決がNosé–Hooverチェーンである．「熱浴を熱浴で温度制御する」という再帰的な発想で，$M$ 段のチェーンでは

\begin{align} \frac{\dd \bm{p}_I}{\dd t} &= \bm{F}_I - \zeta_1 \bm{p}_I, \nonumber\\ \frac{\dd \zeta_1}{\dd t} &= \frac{1}{Q_1}\left[ \sum_I \frac{\abs{\bm{p}_I}^2}{M_I} - N_f k_{\mathrm{B}}T_d \right] - \zeta_1 \zeta_2, \nonumber\\ \frac{\dd \zeta_j}{\dd t} &= \frac{1}{Q_j}\left[ Q_{j-1}\zeta_{j-1}^2 - k_{\mathrm{B}}T_d \right] - \zeta_j \zeta_{j+1} \qquad (j = 2, \dots, M-1), \nonumber\\ \frac{\dd \zeta_M}{\dd t} &= \frac{1}{Q_M}\left[ Q_{M-1}\zeta_{M-1}^2 - k_{\mathrm{B}}T_d \right] \label{eq:16-nh-chain} \end{align}

とする．$j \ge 2$ の式では，1つ前の熱浴変数の「運動エネルギー」$Q_{j-1}\zeta_{j-1}^2$ が目標値 $k_{\mathrm{B}}T_d$(自由度1個分)からずれた分が駆動力になっている．連鎖の各段が余分な自由度を注入することで軌道がカオス的になり，エルゴード性が回復する．$M = 3$〜$5$ で十分であり，追加コストは無視できる(電子状態計算に比べれば $\zeta_j$ の積分は無料である)．現代の第一原理MDコードでは，$NVT$ 計算の既定値がNosé–Hooverチェーンになっていることが多い．

物理的意味:熱浴の選択と「壊してよいもの」

温度制御法にはほかにも，速度を一律に $\sqrt{T_d/T_{\mathrm{inst}}}$ 倍する単純スケーリング，緩和時間で緩やかに補正するBerendsen法，ランダム力と摩擦を加えるLangevin法などがある．判断の基準は「何を壊してよいか」である．単純スケーリングとBerendsen法は正しいカノニカル分布を生成しない(ゆらぎが過小になる)ため，平衡化の初期段階に限って使うべきである．Langevin法はエルゴード性が保証され実装も容易だが，ランダム力が原子の動力学そのものを乱すため，拡散係数や振動スペクトルなど動的な量の計算には向かない．Nosé–Hoover(チェーン)法は決定論的・時間反転対称で，分布も動力学も比較的よく保つため，動的な性質を調べる第一原理MDでは第一選択となる．

16.8 Car–Parrinello法

BOMDの計算コストは，毎ステップの完全なSCF収束が支配していた．1985年にCarとParrinelloが提案した方法は，この「SCFを解く」という手続きそのものを取り除いてしまう．発想は前節のNoséと同じ拡張系である:電子の波動関数に架空の質量を与え，原子核と同時に時間発展させる．うまく設計すれば，波動関数は瞬間ごとの基底状態のまわりを小さく速く振動するだけとなり，原子核から見れば実質的にBorn–Oppenheimer面上を動いているのと変わらなくなる．

16.8.1 拡張Lagrangianと運動方程式

Car–Parrinello法のLagrangianは次のとおりである．

$$ \mathcal{L}_{\mathrm{CP}} = \underbrace{\sum_i \mu \int \abs{\dot{\psi}_i(\rr)}^2 \dd^3 r}_{\text{電子の仮想運動エネルギー}} + \underbrace{\sum_I \frac{M_I}{2} \abs{\dot{\bm{R}}_I}^2}_{\text{核の運動エネルギー}} - \underbrace{E\left[ \{\psi_i\}, \{\bm{R}_I\} \right]}_{\text{Kohn–Shamエネルギー汎関数}} + \underbrace{\sum_{ij} \Lambda_{ij} \left( \int \psi_i^*\psi_j\, \dd^3 r - \delta_{ij} \right)}_{\text{規格直交拘束}} \label{eq:16-cp-lag} $$

$\mu$ は電子波動関数に与えた仮想質量(次元はエネルギー×時間$^2$，原子単位で $100$〜$1000$ 程度の数)，$\Lambda_{ij}$ は規格直交条件 $\int\psi_i^*\psi_j\dd^3 r = \delta_{ij}$ を課すLagrangeの未定乗数(エルミート行列)である．拘束が不可欠であることに注意したい．これがないと，変分エネルギーを下げるためにすべての軌道が最低エネルギー状態へ崩落してしまう(Pauli原理が失われる)．

数学ノート:場に対するEuler–Lagrange方程式

有限個の一般化座標 $q_a$ に対するEuler–Lagrange方程式は

$$ \frac{\dd}{\dd t}\left( \frac{\partial \mathcal{L}}{\partial \dot{q}_a} \right) = \frac{\partial \mathcal{L}}{\partial q_a} \label{eq:16-el-discrete} $$

である(作用 $\int\mathcal{L}\,\dd t$ の停留条件)．力学変数が連続無限個——すなわち各点の場の値 $\psi(\rr)$ ——である場合は，偏微分を汎関数微分に置き換えればよい:

$$ \frac{\dd}{\dd t}\left( \frac{\delta \mathcal{L}}{\delta \dot{\psi}_i^*(\rr)} \right) = \frac{\delta \mathcal{L}}{\delta \psi_i^*(\rr)} . \label{eq:16-el-field} $$

複素場では $\psi_i$ と $\psi_i^*$ を独立変数として扱う(実部と虚部を独立に動かすことと同値;第4章・16.2節の数学ノートと同じ約束である)．必要な汎関数微分は3つで，いずれも直接計算できる:

$$ \frac{\delta}{\delta \dot{\psi}_i^*(\rr)} \sum_j \mu \int \abs{\dot{\psi}_j}^2 \dd^3 r' = \mu\, \dot{\psi}_i(\rr), \qquad \frac{\delta E}{\delta \psi_i^*(\rr)} = \hat{H}_{\mathrm{KS}}\, \psi_i(\rr), \qquad \frac{\delta}{\delta \psi_i^*(\rr)} \sum_{jk}\Lambda_{jk}\!\int\!\psi_j^*\psi_k \dd^3 r' = \sum_j \Lambda_{ij}\psi_j(\rr). \label{eq:16-cp-funcderiv} $$

2番目は第10章で導いたKohn–Sham方程式の変分的な出発点そのものである(占有数 $f_i$ を含む場合は $f_i \hat{H}_{\mathrm{KS}}\psi_i$ となる)．

導出:Car–Parrinelloの運動方程式

Step 1:波動関数の方程式．式\eqref{eq:16-cp-funcderiv}を式\eqref{eq:16-el-field}に代入する．左辺は $\dd(\mu\dot\psi_i)/\dd t = \mu\ddot\psi_i$，右辺はLagrangian\eqref{eq:16-cp-lag}の第3項(符号が負)と第4項の寄与だから

$$ \mu\, \ddot{\psi}_i(\rr, t) = -\hat{H}_{\mathrm{KS}}\, \psi_i(\rr,t) + \sum_j \Lambda_{ij}\, \psi_j(\rr,t) . \label{eq:16-cp-eom-psi} $$

Step 2:原子核の方程式．$\bm{R}_I$ については通常のEuler–Lagrange方程式\eqref{eq:16-el-discrete}から

$$ M_I \ddot{\bm{R}}_I = -\frac{\partial E\left[ \{\psi_i\},\{\bm{R}_I\} \right]}{\partial \bm{R}_I} + \sum_{ij}\Lambda_{ij} \frac{\partial}{\partial \bm{R}_I}\!\int \psi_i^*\psi_j\, \dd^3 r . \label{eq:16-cp-eom-R} $$

最後の項は，基底関数が原子とともに動く場合にだけ現れる(重なり行列の座標微分=16.2節のPulay項と同じ起源である)．平面波基底では基底が原子位置を知らないので厳密にゼロであり，CP法が平面波基底とともに発展した理由の一つがこれである．

Step 3:未定乗数を決める．$\Lambda_{ij}$ は，拘束 $\int\psi_i^*\psi_j\dd^3 r = \delta_{ij}$ が時間発展の各瞬間で成り立つように決まる．この条件を2回微分した式に\eqref{eq:16-cp-eom-psi}を代入すると $\Lambda$ についての行列方程式が得られ，数値的にはSHAKE/RATTLEと呼ばれる反復法で毎ステップ解く．

Step 4:静止極限を確認する．$\ddot\psi_i = 0$(電子が止まっている)とすると式\eqref{eq:16-cp-eom-psi}は $\hat{H}_{\mathrm{KS}}\psi_i = \sum_j \Lambda_{ij}\psi_j$，すなわち $\Lambda$ を対角化する基底では $\hat{H}_{\mathrm{KS}}\psi_i = \varepsilon_i \psi_i$ ——Kohn–Sham方程式そのものである．すなわちCP動力学の「静止点」は基底状態であり，$-\hat{H}_{\mathrm{KS}}\psi_i + \sum_j\Lambda_{ij}\psi_j$ は「波動関数に働く力」と読める．この力をゼロにする代わりに，質量 $\mu$ を与えて振動させておくのがCP法の核心である．

したがって全体の保存量は

$$ E_{\mathrm{cons}}^{\mathrm{CP}} = \underbrace{\sum_i \mu \int \abs{\dot{\psi}_i}^2 \dd^3 r}_{\displaystyle E_{\mathrm{kin}}^{e}} + \sum_I \frac{M_I}{2}\abs{\dot{\bm{R}}_I}^2 + E\left[ \{\psi_i\},\{\bm{R}_I\}\right] \label{eq:16-cp-cons} $$

である．ここで$E_{\mathrm{kin}}^{e}$ は物理的な意味を持たない架空の量であることに注意したい．これが小さく，かつ時間とともに増大しないこと——それが「電子が基底状態に張り付いている」ことの数値的な判定基準になる．

16.8.2 断熱条件:電子はなぜ核に追随するのか

CP法が正しく機能するには，電子の仮想振動が核の運動より十分速い必要がある．速い自由度と遅い自由度が振動数的に分離していれば，両者の間のエネルギー交換は指数関数的に抑制される(断熱定理)からである．この条件を定量化しよう．

導出:電子の仮想振動数 $\omega \propto \sqrt{\Delta\varepsilon/\mu}$

Step 1:基底状態からのずれをパラメータ化する．Hamiltonianを凍結し(核が止まっている極限)，占有軌道 $\psi_i^0$ が空軌道 $\psi_a^0$ の方向へ少しだけ混ざった試行関数を考える:

$$ \psi_i = \frac{\psi_i^0 + c\,\psi_a^0}{\sqrt{1 + c^2}}, \qquad c \in \mathbb{R},\; \abs{c}\ll 1 . \label{eq:16-cp-mix} $$

分母は規格化のためである($\psi_i^0, \psi_a^0$ は規格直交)．この形なら拘束は自動的に満たされる．

Step 2:エネルギーを $c$ の2次まで展開する．$\hat{H}_{\mathrm{KS}}\psi_i^0 = \varepsilon_i\psi_i^0$，$\hat{H}_{\mathrm{KS}}\psi_a^0 = \varepsilon_a \psi_a^0$ を使うと

\begin{align} \braket{\psi_i | \hat{H}_{\mathrm{KS}} | \psi_i} = \frac{\varepsilon_i + c^2 \varepsilon_a}{1 + c^2} = \left( \varepsilon_i + c^2\varepsilon_a \right)\left( 1 - c^2 + O(c^4) \right) = \varepsilon_i + c^2\left( \varepsilon_a - \varepsilon_i \right) + O(c^4). \label{eq:16-cp-energy-c} \end{align}

1つ目の等号では交差項 $\braket{\psi_i^0|\hat{H}|\psi_a^0} = \varepsilon_a\braket{\psi_i^0|\psi_a^0} = 0$(直交性)が消えることを使い，2つ目では $1/(1+c^2)$ を幾何級数展開した．$\varepsilon_a > \varepsilon_i$ だから，これは $c = 0$ を底とする放物線である．

Step 3:仮想運動エネルギーを $\dot c$ で表す．式\eqref{eq:16-cp-mix}を時間微分すると，最低次で $\dot\psi_i = \dot{c}\,\psi_a^0 + O(c\dot c)$ だから

$$ \mu \int \abs{\dot\psi_i}^2 \dd^3 r = \mu\, \dot{c}^2 \int \abs{\psi_a^0}^2 \dd^3 r = \mu\, \dot{c}^2 . \label{eq:16-cp-kin-c} $$

Step 4:1自由度のLagrangianにまとめて解く．定数項を落とすと

$$ \mathcal{L}(c, \dot c) = \mu \dot{c}^2 - \left( \varepsilon_a - \varepsilon_i \right) c^2 . \label{eq:16-cp-lag-c} $$

Euler–Lagrange方程式\eqref{eq:16-el-discrete}を適用すると $\dd(2\mu\dot c)/\dd t = -2(\varepsilon_a - \varepsilon_i)c$，すなわち

$$ \ddot{c} = -\frac{\varepsilon_a - \varepsilon_i}{\mu}\, c \qquad \Longrightarrow \qquad \omega_{ia} = \sqrt{\frac{\varepsilon_a - \varepsilon_i}{\mu}} . \label{eq:16-cp-omega} $$

これは単振動である．電子の仮想振動数は，軌道エネルギー差の平方根に比例し，仮想質量の平方根に反比例する．(仮想運動エネルギーを $\frac{\mu}{2}\int\abs{\dot\psi}^2$ と定義する流儀では $\omega_{ia} = \sqrt{2(\varepsilon_a-\varepsilon_i)/\mu}$ となる．$\mu$ の定義が因子2だけ違うだけで物理は同じである．)

式\eqref{eq:16-cp-omega}から，電子の仮想振動数のスペクトルは次の2つの端で挟まれる:

$$ \omega_{\min}^{e} = \sqrt{\frac{E_{\mathrm{gap}}}{\mu}} \quad (\text{HOMO–LUMOギャップ}), \qquad \omega_{\max}^{e} \simeq \sqrt{\frac{E_{\mathrm{cut}}}{\mu}} \quad (\text{基底のカットオフ}). \label{eq:16-cp-spectrum} $$

ここから，CP法の設計を縛る2つの条件が同時に出る．

$$ \text{(断熱条件)}\quad \omega_{\min}^{e} = \sqrt{\frac{E_{\mathrm{gap}}}{\mu}} \;\gg\; \omega_{\max}^{\mathrm{ion}}, \qquad \text{(安定条件)}\quad \Delta t \;\lesssim\; \frac{2}{\omega_{\max}^{e}} = 2\sqrt{\frac{\mu}{E_{\mathrm{cut}}}} . \label{eq:16-cp-conditions} $$

前者は「$\mu$ を小さくせよ」，後者は「$\mu$ を小さくすると時間刻みも小さくせよ($\Delta t \propto \sqrt{\mu}$)」と言っている．この綱引きが $\mu$ の選択を決める．

図16.10 Car–Parrinello法の成立条件．(a) 電子の仮想振動数の帯($\sqrt{E_{\mathrm{gap}}/\mu}$ から $\sqrt{E_{\mathrm{cut}}/\mu}$ まで)が核の振動数の帯より十分高い位置になければならない．仮想質量 $\mu$ を下げると電子の帯は右(高振動数側)へ動いて安全になるが，同時に時間刻みを $\sqrt{\mu}$ に比例して小さくしなければならない．(b) 断熱条件が満たされていれば，CP軌道はBorn–Oppenheimer面のすぐ上を小さく速く振動し，核から見た平均的な運動はBOMDと変わらない．

例:典型的なパラメータの見積もり

半導体的な系($E_{\mathrm{gap}} = 4$ eV $= 0.147$ Hartree)，平面波カットオフ $E_{\mathrm{cut}} = 40$ Hartree(80 Ry)，水素を含む系($\omega_{\max}^{\mathrm{ion}} = 3000$ cm$^{-1} = 0.0137$ 原子単位)を考える．$\mu = 400$ 原子単位とすると，式\eqref{eq:16-cp-spectrum}より

$$ \omega_{\min}^{e} = \sqrt{\frac{0.147}{400}} = 0.019\ \text{a.u.}, \qquad \omega_{\max}^{e} = \sqrt{\frac{40}{400}} = 0.32\ \text{a.u.} $$

である．断熱ギャップの比は $\omega_{\min}^{e}/\omega_{\max}^{\mathrm{ion}} = 1.4$ ——決して大きくない．時間刻みは安定条件から $\Delta t \lesssim 2/0.32 = 6.3$ 原子単位 $= 0.15$ fs であり，実際のCP計算で $\Delta t = 5$〜$7$ 原子単位が使われるのと一致する．BOMDの $0.5$ fs に比べて3〜4倍細かいが，SCFを解かない分だけ1ステップは大幅に安い，というのがCP法の損益計算である．

この見積もりは，断熱条件がぎりぎりで成立していることも教えてくれる．分離を稼ぐ実用的な処方が2つある．(i) $\mu$ を下げる($\mu = 200$ なら $\omega_{\min}^{e}$ は $\sqrt{2}$ 倍になる．ただし $\Delta t$ も $1/\sqrt{2}$ 倍)．(ii) 水素を重水素で置き換えて $\omega_{\max}^{\mathrm{ion}}$ を $1/\sqrt{2}$ 倍にする．CP計算でしばしば重水素質量が使われるのはこのためである．なお，電子の仮想的な慣性は原子核の有効質量をわずかに増やす(質量繰り込み)ため，振動数の定量的な議論では $\mu \to 0$ の外挿が必要になることもある．

16.8.3 金属で破綻する理由と，BOMDとの使い分け

金属では定義により $E_{\mathrm{gap}} = 0$ である．すると式\eqref{eq:16-cp-spectrum}より $\omega_{\min}^{e} = 0$ ——電子の仮想振動数の帯が核の帯と重なってしまう．共鳴によって核から電子へエネルギーが流れ続け，$E_{\mathrm{kin}}^{e}$ が単調に増大する．物理的には「電子が加熱され，Born–Oppenheimer面から離れていく」ことを意味し，核の軌道も系統的に誤ったものになる．同じ問題は，狭ギャップ半導体，high-spin遷移金属錯体，化学反応の途中でギャップが閉じる場面(結合の切断は往々にしてこれを伴う)でも起こる．

対処法は2つある．第一は，電子系にもNosé–Hoover熱浴を付けて $E_{\mathrm{kin}}^{e}$ を小さな目標値に固定する2熱浴法である(核用の熱浴で温度 $T_d$ を，電子用の熱浴で $E_{\mathrm{kin}}^{e}$ を制御する)．これは有効だが，パラメータが増え，電子熱浴が動力学に与える影響の評価が難しい．第二は素直にBOMDを使うことである．両者の比較を整理しておこう．

表16.1 Born–Oppenheimer MD と Car–Parrinello MD の比較
	BOMD	CPMD
電子状態	毎ステップSCF収束	拡張Lagrangianで同時に時間発展
1ステップの費用	SCF反復 $10$〜$30$ 回	力の評価 $1$ 回相当
時間刻み	核の振動で決まる($0.5$〜$2$ fs)	電子の仮想振動で決まる($0.1$〜$0.2$ fs)
エネルギー面	厳密にBO面(SCF収束の範囲で)	BO面のすぐ上を振動
金属・小ギャップ系	問題なし(占有数のスメアリングで対応)	断熱分離が破綻(2熱浴などの工夫が必要)
局在基底との相性	良い	拘束と基底の座標依存性が絡み複雑
主な誤差源	SCF収束不足によるドリフト	断熱性の破れ，質量繰り込み

歴史的にはCP法が第一原理MDを実用化した立役者だが，現在では両者は使い分けられている．SCFの初期推定を時間反転対称に外挿する技術(16.6.5項)が発達し，BOMDの1ステップあたりのSCF反復数が大きく減ったため，大きな時間刻みが使えるBOMDの総コストがCPMDに匹敵するようになったからである．絶縁体・分子液体の長時間シミュレーションではCPMDが今も強く，金属や結合の生成消滅を伴う反応ではBOMDが選ばれる，というのが現在の大まかな棲み分けである．

16.9 発展的手法の概観:メタダイナミクスとQM/MM

ここまでの道具立てには，まだ2つの大きな限界が残っている．第一に，第一原理MDが直接追える時間は $10$ ps〜$1$ ns 程度であり，興味のある化学反応の多くはそれよりはるかに遅い(稀事象問題;16.5節の表を参照)．第二に，第一原理計算で扱える原子数はせいぜい $10^2$〜$10^3$ 個であり，酵素・電池電解液・触媒界面のような「反応中心は小さいが環境が巨大」な系にはそのままでは届かない．本節では，この2つの壁をそれぞれ突破する代表的な手法を，原理が分かる程度に紹介する．

16.9.1 メタダイナミクス:自由エネルギー面を埋めて描く

まず，何を計算したいのかを明確にしよう．反応の進行度を表す少数個の変数 $\bm{s} = (s_1(\bm{R}), \dots, s_d(\bm{R}))$——結合長，配位数，二面角，あるいはそれらの組合せ——を選ぶ．これを集団変数(collective variable, CV)と呼ぶ．CVの値が $\bm{s}$ である確率は，カノニカル分布\eqref{eq:16-canonical}を残りの自由度について積分して

$$ P(\bm{s}) = \frac{1}{Z}\int \dd \bm{R}\; \delta\!\left( \bm{s} - \bm{s}(\bm{R}) \right)\, e^{-\beta U(\bm{R})} \label{eq:16-cv-prob} $$

で与えられる．これに対応する自由エネルギー面を

$$ F(\bm{s}) = -k_{\mathrm{B}}T \ln P(\bm{s}) \quad (+\;\text{定数}) \label{eq:16-fes-def} $$

と定義する．$F(\bm{s})$ こそが求めたいものである．NEB法が与える $E_a$ が「0 Kのエネルギー障壁」であったのに対し，$F(\bm{s})$ はエントロピーを含む有限温度の障壁であり，実験の反応速度と直接比較すべき量である．ところが式\eqref{eq:16-cv-prob}を素朴なMDで推定しようとすると，障壁の向こう側は $e^{-\beta \Delta F}$ という指数的に小さい確率でしか訪れないため，統計が全く溜まらない．

LaioとParrinelloが2002年に提案したメタダイナミクスの発想は，訪れた場所に「砂」を積んで埋めてしまうことである．時刻 $t$ までに一定間隔 $\tau$ で訪れたCV値 $\bm{s}(t')$ に，幅 $\delta s$，高さ $W$ のGauss関数を堆積させ，それをバイアスポテンシャルとして系に加える:

$$ V_{\mathrm{bias}}(\bm{s}, t) = \sum_{k\tau\, <\, t} W \exp\!\left( -\frac{\abs{\bm{s} - \bm{s}(\bm{R}(k\tau))}^2}{2\,\delta s^2} \right) \label{eq:16-meta-bias} $$

系が感じるポテンシャルは $U(\bm{R}) + V_{\mathrm{bias}}(\bm{s}(\bm{R}), t)$ になる(原子に働く力には連鎖律で $-\partial V_{\mathrm{bias}}/\partial \bm{s}\cdot\partial\bm{s}/\partial\bm{R}$ が加わる)．谷にとどまるほどその谷が埋まっていくので，系は自分が来た道を「登れなく」なり，やがて障壁を越える．

驚くべきなのは次の点である．十分に長い時間ののち，堆積したバイアスは自由エネルギー面の符号を反転したものに収束する:

$$ V_{\mathrm{bias}}(\bm{s}, t \to \infty) \simeq -F(\bm{s}) + \text{定数}. \label{eq:16-meta-converge} $$

すなわち探索を加速する道具が，そのまま自由エネルギー面の測定器になっている．この関係の根拠は，次の自己無撞着性の議論で理解できる．バイアスが掛かった系でのCVの分布は，有効自由エネルギー $F(\bm{s}) + V_{\mathrm{bias}}(\bm{s})$ に対するBoltzmann分布

$$ P_{V}(\bm{s}) \propto \exp\!\left[ -\beta \left( F(\bm{s}) + V_{\mathrm{bias}}(\bm{s}) \right) \right] \label{eq:16-meta-biased} $$

である．もし $V_{\mathrm{bias}} = -F + C$ が実現したとすると，指数の中は定数になり $P_V(\bm{s})$ は一様になる．一様ということは，以後のGauss関数の堆積はCV空間のどこにも同じ割合で積まれる，すなわち $V_{\mathrm{bias}}$ に一様な定数を足すだけである．それは $V_{\mathrm{bias}} = -F + C'$ という形を保つ．つまり$V_{\mathrm{bias}} = -F + \text{定数}$ はこの堆積過程の不動点である．逆に，まだ埋まっていない谷があればそこに系が長く滞在して優先的に埋められるから，力学は不動点へ向かう．（厳密な収束の証明と誤差評価は2006年以降に整備された．）

図16.11 メタダイナミクスの原理．青い実線が真の自由エネルギー面 $F(s)$．系が滞在した位置に順次Gauss関数(黄色)が堆積し，谷が埋まっていく．十分に埋まると有効自由エネルギー $F(s) + V_{\mathrm{bias}}(s)$ は平らになり(緑破線)，系はCV空間を自由に拡散できるようになる．このとき $V_{\mathrm{bias}}(s) \simeq -F(s) + $定数であり，堆積したバイアスそのものが自由エネルギー面の測定値になる．

実用上の注意を3つ挙げる．

集団変数の選択がすべてを決める．反応の本質的な自由度がCVに含まれていないと，その「隠れた変数」の緩和が遅いために，埋め立てが見かけ上進んでも本当の自由エネルギー面は再現されない．CVの選び方に一般的な処方はなく，化学的な洞察が要求される．
収束の判定が難しい．式\eqref{eq:16-meta-converge}は $t\to\infty$ の主張であり，有限時間では $V_{\mathrm{bias}}$ が $-F$ のまわりを揺らぎ続ける．改良版のwell-tempered metadynamicsは，堆積するGauss関数の高さを既に積まれたバイアスに応じて $W \to W\,e^{-\beta V_{\mathrm{bias}}/(\Delta T/T)}$ のように減衰させることでこの揺らぎを抑え，$V_{\mathrm{bias}} \to -\frac{\Delta T}{T + \Delta T} F(\bm{s})$ という制御された形に収束させる($\Delta T$ は探索したいエネルギー範囲を決めるパラメータ)．
計算量は依然として大きい．第一原理MDと組み合わせる場合，必要なステップ数は数十〜数百 ps 相当になる．実務では，まず古典力場のMDで系を平衡化し，反応中心だけを量子力学的に扱い(次項)，そのうえでメタダイナミクスを走らせる，という三段構えが取られる．

16.9.2 QM/MM法:反応中心を量子力学で，環境を古典力学で

酵素の活性部位，溶液中の反応，固体中の欠陥——いずれも「化学結合の組み替えが起こるのは数十原子の領域だけで，残りの数万原子は結合の組み替えを伴わずに静電場と立体障害を提供している」という構造を持つ．ならば，系を2つに分割して，必要なところにだけ量子力学を使えばよい．これがQM/MM法である．

全エネルギーの分割には2つの流儀がある．加算的な定式化では

$$ E_{\mathrm{tot}} = E_{\mathrm{QM}}\left[ \text{QM領域};\, \{q_m\} \right] + E_{\mathrm{MM}}\left[ \text{MM領域} \right] + E_{\mathrm{QM/MM}}^{\mathrm{vdW + bonded}} \label{eq:16-qmmm-additive} $$

と書く．第1項は「MM領域の点電荷 $\{q_m\}$ が作る外場の中で解いたQM領域の全エネルギー」，第3項は境界をまたぐvan der Waals相互作用と結合項(古典力場で扱う)である．一方減算的(ONIOM型)な定式化では

$$ E_{\mathrm{tot}} = E_{\mathrm{MM}}\left[ \text{系全体} \right] + E_{\mathrm{QM}}\left[ \text{QM領域} \right] - E_{\mathrm{MM}}\left[ \text{QM領域} \right] \label{eq:16-qmmm-subtractive} $$

とする．「全体を古典で計算し，QM領域の古典的記述を量子的記述で置き換える」という読み方ができ，実装が容易である．

静電埋め込み．QM/MM法の質を決めるのは，環境がQM領域を分極させる効果をどう入れるかである．最も標準的な静電埋め込みでは，MM領域の点電荷をKohn–Shamハミルトニアンの外部ポテンシャルに直接加える:

$$ v_{\mathrm{ext}}(\rr) \;\longrightarrow\; v_{\mathrm{ext}}(\rr) - \sum_{m \in \mathrm{MM}} \frac{q_m}{\abs{\rr - \bm{R}_m}}, \qquad E_{nn} \;\longrightarrow\; E_{nn} + \sum_{I \in \mathrm{QM}} \sum_{m \in \mathrm{MM}} \frac{Z_I q_m}{\abs{\bm{R}_I - \bm{R}_m}} . \label{eq:16-qmmm-embed} $$

(電子の電荷が $-1$ であることから，電子と点電荷 $q_m$ の相互作用エネルギーは $-\int n(\rr) q_m/\abs{\rr-\bm{R}_m}\dd^3 r$ となり，1体ポテンシャルとしては上の第1式の形になる．)これにより，SCFのたびにQM領域の電子密度が環境の電場に応答して分極する．より粗い「機械的埋め込み」(QM計算を真空中で行い，静電相互作用は点電荷同士で古典的に扱う)では，この分極が全く入らない．酵素の活性部位のように強い電場がかかる環境では，両者の反応障壁の差は $10$ kcal/mol にも達しうる．

リンク原子法．境界が共有結合を横切るとき(タンパク質の主鎖と側鎖の間など)，QM領域の端の原子には切断された結合手が残ってしまう．最も単純で広く使われる処方がリンク原子法である:切断した結合の方向に沿った適切な距離のところに水素原子を1個置き，QM計算のときだけこの「リンク原子」を含めて閉殻構造を作る．守るべき経験則がいくつかある．

切断してよいのは非極性の単結合(典型的にはC–C)だけである．二重結合・芳香環・極性結合(C–N, C–Oなど)を切ると，電子構造の記述が大きく歪む．
リンク原子のすぐ近くのMM点電荷は，そのままでは非物理的に強くリンク原子を分極させる．境界近傍の電荷を削除する・隣接原子へ振り分ける・広がりを持たせる(smearing)といった補正が必要である．
境界は反応中心から十分離す．目安として，結合の組み替えが起こる原子から2〜3結合以上離す．

図16.12 QM/MM法の構成．反応中心だけをDFTで，周囲の巨大な環境を古典力場で扱う．境界が共有結合を横切る場合は，QM計算の際にリンク原子(水素)で結合手を埋める．MM領域の点電荷をKohn–Shamハミルトニアンの外部ポテンシャルに加える静電埋め込みにより，環境の電場によるQM領域の分極が取り込まれる．

物理的意味:なぜ酵素反応にこの組合せが要るのか

酵素は同じ化学反応を，水溶液中に比べて $10^{10}$ 倍以上速く進める．その理由は，活性部位が作り出す特異的な静電場が遷移状態を選択的に安定化するからだ，というのが現在の理解である．この描像を第一原理から検証するには，(i) 結合の切断・生成を扱えること(→ DFT)，(ii) タンパク質全体の静電場と立体構造を含めること(→ MM)，(iii) 有限温度でのゆらぎと自由エネルギーを扱えること(→ MD)，(iv) 秒スケールの反応を扱えること(→ メタダイナミクス)がすべて必要になる．本章で見てきた道具は，こうして一つの計算の中で組み合わされる．1998年と2013年のノーベル化学賞が，それぞれ密度汎関数法・量子化学の計算手法と，生体分子のマルチスケール模型に与えられたのは，この流れの到達点を示している．

16.10 まとめ・演習・参考文献

16.10.1 まとめ

ポテンシャルエネルギー面(PES):Born–Oppenheimer近似のもとで，全エネルギーは核座標 $3M$ 個の関数 $E(\{\bm{R}_I\})$ である．その勾配の符号を反転したものが力，2階微分がHessianであり，Hessianの固有値の符号で極小点(安定構造)と1次の鞍点(遷移状態)が分類される．本章の全手法は「この面の上をどう歩くか」の変奏である．
Hellmann–Feynmanの定理:$\Psi$ が厳密固有状態で規格化がパラメータによらないなら，$\dd E/\dd\lambda = \braket{\Psi|\partial\hat{H}/\partial\lambda|\Psi}$．波動関数の変化による寄与は変分原理により厳密に消える．$\lambda = \bm{R}_I$ と取ると，力は電子密度と原子核の古典静電相互作用だけで書ける(静電定理，式\eqref{eq:16-hf-force})．
Pulay力:基底関数が原子とともに動く場合，変分空間そのものが $\bm{R}_I$ に依存するため補正項\eqref{eq:16-pulay-force}が必要になる．平面波基底ではこの項は厳密にゼロであり，原子中心基底では必須である．
応力テンソル:一様ひずみ $\bm{\varepsilon}$ に対するエネルギー微分 $\sigma_{\alpha\beta} = \Omega^{-1}\partial E/\partial\varepsilon_{\alpha\beta}$．波動関数のスケーリング変換によって定義域の変化を吸収して計算する．トレースは圧力を与え，平衡構造では $2T = -E_{\mathrm{C}}$ というビリアル定理に帰着する．
構造最適化:最急降下法は条件数 $\kappa$ に比例するステップ数を要し実用にならない($((\kappa-1)/(\kappa+1))^n$ の収束)．Newton法は2次収束するがHessianが高価で，不定値だと鞍点へ落ちる．BFGSはセカント条件+対称性+ランク2更新から導かれ，曲率条件 $\bm{y}^{\mathsf{T}}\bm{s}>0$ のもとで正定値性を保つ．この曲率条件はWolfe条件付きの直線探索が自動的に保証する．信頼半径の拘束はレベルシフト $(H+\mu\bm{1})\bm{\Delta} = -\bm{g}$ を与える．収束判定は最大力 $10^{-3}$〜$10^{-4}$ Hartree/Bohr が標準だが，その意味は系のばね定数に依存する．
格子振動とフォノン(16.4.5項):平衡点まわりの2次展開で1次項が消えることから $E = E_0 + \frac12\sum\Phi u u$ が得られ，力の定数行列 $\Phi$ はHessianそのものである．剛体並進で力が変わらないことから総和則 $\sum_j \Phi_{i\alpha,j\beta} = 0$ が従う．周期系では平面波型の変位を代入することで運動方程式が動的行列\eqref{eq:16-ph-dyn-def}の固有値問題\eqref{eq:16-ph-eig}に分解し，単位胞あたり $n$ 原子なら $3n$ 本の分枝(うち3本が音響枝，$3n-3$ 本が光学枝)が現れる．ある $\bm{q}$ で $\omega^2 \lt 0$ ならその固有ベクトル方向にエネルギーが下がり，構造は動的に不安定である．$\Phi$ の計算法は有限変位法(スーパーセルが相互作用の到達距離を決める)とDFPT(線形応答としてスーパーセル不要で任意の $\bm{q}$)の2つ．調和自由エネルギーから零点振動エネルギー・比熱・(準調和近似で)熱膨張が得られる．
NEB法:バネで結んだ像の列に対し，真の力の垂直成分+バネ力の平行成分だけを使う(式\eqref{eq:16-neb-force})．これが素朴な弾性バンド法の2つの病(コーナーカットと像のずり落ち)を同時に治す．接線はエネルギーの高い側の差分から作る．CI-NEBは最高点の像のバネを切り接線方向の力を反転させることで，鞍点そのものを最適化と同じ精度で決定する．得られた $E_a$ は遷移状態理論 $k = \nu e^{-E_a/k_{\mathrm{B}}T}$ を通じて反応速度に翻訳される．
第一原理分子動力学:BOMDは毎ステップSCFを収束させてから核を動かす．Verlet法は前後対称なTaylor展開の和から得られ，局所誤差 $O(\Delta t^4)$，時間反転対称かつシンプレクティックである．後者により影のハミルトニアンが厳密に保存され，エネルギー誤差は有界に振動するだけでドリフトしない．安定条件は $\omega_{\max}\Delta t < 2$，実用的な時間刻みは最高振動周期の $1/20$(水素を含む系で $0.5$ fs)．エネルギー保存のモニタが計算の健全性の最良の指標である．
Nosé–Hoover法:熱浴を力学変数 $s$(質量 $Q$，対数ポテンシャル $g k_{\mathrm{B}}T_d\ln s$)として拡張系を作り，仮想時間から実時間へ変数変換すると $\dot{\bm{p}} = \bm{F} - \zeta\bm{p}$，$\dot\zeta = (\sum p^2/M - g k_{\mathrm{B}}T_d)/Q$ が得られる．保存量\eqref{eq:16-nh-cons}が存在し，拡張系のミクロカノニカル分布を射影するとカノニカル分布になる($g = N_f$，仮想時間サンプリングなら $N_f+1$)．少数自由度系での非エルゴード性はNosé–Hooverチェーンで解決される．
Car–Parrinello法:波動関数に仮想質量 $\mu$ を与えて核と同時に時間発展させる．電子の仮想振動数は $\omega \propto \sqrt{\Delta\varepsilon/\mu}$ であり，断熱条件 $\sqrt{E_{\mathrm{gap}}/\mu} \gg \omega_{\max}^{\mathrm{ion}}$ と安定条件 $\Delta t \lesssim 2\sqrt{\mu/E_{\mathrm{cut}}}$ が $\mu$ の選択を挟み撃ちにする．金属($E_{\mathrm{gap}}=0$)では断熱分離が成立せず破綻する．
発展的手法:メタダイナミクスは集団変数空間にGauss関数を堆積させて谷を埋め，$V_{\mathrm{bias}} \simeq -F(\bm{s})$ という不動点に到達することで，稀事象の加速と自由エネルギー面の再構成を同時に達成する．QM/MM法は反応中心をDFT，環境を古典力場で扱い，リンク原子で切断結合を埋め，静電埋め込みで環境による分極を取り込む．

本章で電子状態理論は，原子核の運動という「動く世界」に接続された．ただしここまで扱ってきたのはすべて平衡状態あるいはその近傍の話である．次章では，系の両端に異なる化学ポテンシャルを持つ電子浴を接続し，電流が定常的に流れる非平衡の問題——非平衡グリーン関数法による量子輸送——へ進む．そこでは，本章で構造を決めたナノ構造が，実際にどのような伝導特性を示すかが主題になる．

16.10.2 演習問題

演習16.1 最急降下法の条件数依存性とBFGS更新

(a) 2次形式\eqref{eq:16-sd-model}に対し，一般の出発点 $\bm{x}_0 = (u, v)^{\mathsf{T}}$ から厳密直線探索を行うときの歩幅 $\alpha_0$ を式\eqref{eq:16-sd-alpha}から求めよ．また，続く2ステップの変位が互いに直交する($\bm{g}_1^{\mathsf{T}}\bm{g}_0 = 0$)ことを示し，図16.3のジグザグがなぜ生じるかを説明せよ．

(b) $\kappa = 100$ の系で誤差を $10^{-3}$ 倍に減らすのに必要な最急降下法のステップ数を式\eqref{eq:16-sd-steps}で見積もれ．1ステップが10秒のSCF計算だとすると何時間かかるか．

(d) 1変数の2次関数 $E(x) = \frac{1}{2}h x^2 + bx$ にBFGS更新式\eqref{eq:16-bfgs}を1回だけ適用すると，$B_1$ が厳密なHessian $h$ に一致することを示せ．

ヒント: (a) 厳密直線探索の停留条件は「進んだ先での勾配が探索方向に直交すること」である．(c) 連鎖律より $\partial^2 E/\partial x'_i \partial x'_j = (H^{-1/2} H H^{-1/2})_{ij}$．(d) 1変数では $y = h\,s$ が厳密に成り立つので，式\eqref{eq:16-bfgs}の3項を直接計算すればよい．

演習16.2 Verlet法の性質

(a) 位置Verlet法\eqref{eq:16-verlet-pos}が時間反転($\Delta t \to -\Delta t$)に対して不変であることを示せ．またこの性質から，数値誤差を無視すれば軌道を逆にたどって初期条件に戻れることを説明せよ．

(b) 調和振動子に対する速度Verlet法の写像行列 $A(\theta)$(式\eqref{eq:16-vv-harmonic})を用いて，数値解の振動数 $\omega_{\mathrm{num}}$ が真の $\omega$ からどれだけずれるかを $\theta = \omega\Delta t$ の小ささについて2次まで求めよ．$\theta = 0.31$(周期の $1/20$)のとき，相対誤差は何%か．

ヒント: (b) $\det A = 1$ なので安定領域では固有値は $\lambda = e^{\pm i\theta_{\mathrm{num}}}$ と書け，$\mathrm{Tr}A = 2\cos\theta_{\mathrm{num}}$ である．$\cos\theta_{\mathrm{num}} = 1 - \theta^2/2$ を $\theta_{\mathrm{num}}$ について展開せよ．(c) 局所誤差は $O(\Delta t^4)$，影のハミルトニアンの補正は $O(\Delta t^2)$ である．

演習16.3 Nosé–Hoover熱浴

(a) Nosé–Hoover方程式\eqref{eq:16-nh-eom}において，$\zeta$ が長時間平均でゼロになることを，保存量\eqref{eq:16-nh-cons}が有界にとどまるという要請から論ぜよ．またそれが「長時間平均で $T_{\mathrm{inst}} = T_d$」を意味することを示せ．

(b) 自由度数 $g$ を誤って $N_f$ ではなく $3N$(重心の並進を除いていない)と取ったとき，$N = 100$ の系で平均温度は何%ずれるか見積もれ．

(d) 熱浴を1個だけ付けた1次元調和振動子がなぜエルゴード的でないのかを，拡張系の相空間の次元と保存量の個数の観点から論ぜよ．

ヒント: (a) 保存量の中の $g k_{\mathrm{B}}T_d \ln s$ の項に注目し，$\dd\ln s/\dd t = \zeta$ を時間積分せよ．$\ln s$ が線形に増大すれば保存量は有界でいられない．(b) $\braket{T_{\mathrm{inst}}}$ は $g/N_f$ 倍になる．(c) $k_{\mathrm{B}}T_d$ を Hartree に，$\tau_{\mathrm{th}}$ を原子単位時間($1$ fs $= 41.34$ a.u.)に直してから掛ければよい．(d) 4次元の相空間に2つの保存量があると，軌道は2次元トーラス上に閉じ込められうる．

演習16.4 CI-NEBと遷移状態理論

(a) CI力\eqref{eq:16-cineb}がゼロになる点では $\nabla E = \bm{0}$ であることを示した(式\eqref{eq:16-cineb-zero})．逆に，なぜその停留点が極小点や2次以上の鞍点ではなく1次の鞍点だと期待できるのか，力学的な安定性の観点から説明せよ．

(b) 活性化エネルギー $E_a = 0.8$ eV，前因子 $\nu = 10^{13}$ s$^{-1}$ の反応について，$300$ K と $500$ K での速度定数を計算し，その比を求めよ．

(d) NEB法の力\eqref{eq:16-neb-force}が何らかのスカラー関数の勾配になっていないことを，$\hat{\bm{\tau}}_i$ が $\{\bm{R}_j\}$ に依存することに注意して論ぜよ．またそのことが最適化アルゴリズムの選択にどう影響するか述べよ．

ヒント: (a) CI力の下では，接線方向は「エネルギーを上る」力学，垂直方向は「下る」力学になっている．安定な停留点になるためには，それぞれの方向でどんな曲率が必要か．(b) $k_{\mathrm{B}}T = 25.9$ meV($300$ K)，$43.1$ meV($500$ K)．(c) $e^{0.15/0.0259}$ を計算せよ．

16.10.3 参考文献

力とPulay補正

H. Hellmann, Einführung in die Quantenchemie (Deuticke, Leipzig, 1937).
R. P. Feynman, "Forces in Molecules," Phys. Rev. 56, 340 (1939). ——静電定理を明示した原論文．
P. Pulay, "Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules I," Mol. Phys. 17, 197 (1969). ——基底関数の座標依存性に由来する補正項の原論文．

応力テンソル

O. H. Nielsen and R. M. Martin, "Quantum-mechanical theory of stress and force," Phys. Rev. B 32, 3780 (1985);同 Phys. Rev. Lett. 50, 697 (1983).
R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004), Chap. 19. ——力と応力の実装を含む標準的な教科書．

構造最適化

J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. (Springer, 2006). ——BFGS・信頼領域法・Wolfe条件の標準的な数学的文献．
C. G. Broyden, J. Inst. Math. Appl. 6, 76 (1970);R. Fletcher, Comput. J. 13, 317 (1970);D. Goldfarb, Math. Comput. 24, 23 (1970);D. F. Shanno, Math. Comput. 24, 647 (1970). ——BFGS更新式の4つの独立な原論文．
H. B. Schlegel, "Estimating the Hessian for gradient-type geometry optimizations," Theor. Chim. Acta 66, 333 (1984). ——化学的知識に基づく初期Hessianの構成．

格子振動とフォノン

M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, 1954). ——格子動力学の古典的な標準書．動的行列・総和則・極性結晶の長距離項の原典的な扱い．
S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, "Phonons and related crystal properties from density-functional perturbation theory," Rev. Mod. Phys. 73, 515 (2001). ——DFPTの標準的な総説．16.4.5項(8)(b)の内容の全体像．
X. Gonze and C. Lee, "Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory," Phys. Rev. B 55, 10355 (1997). ——実装の詳細とLO–TO分裂の扱い．
K. Parlinski, Z. Q. Li, and Y. Kawazoe, "First-principles determination of the soft mode in cubic ZrO$_2$," Phys. Rev. Lett. 78, 4063 (1997). ——有限変位法(直接法)によるソフトモードの決定．
A. Togo and I. Tanaka, "First principles phonon calculations in materials science," Scr. Mater. 108, 1 (2015). ——有限変位法の実務的な解説と実装．

反応経路探索

H. Jónsson, G. Mills, and K. W. Jacobsen, "Nudged elastic band method for finding minimum energy paths of transitions," in Classical and Quantum Dynamics in Condensed Phase Simulations (World Scientific, 1998), p. 385.
G. Henkelman and H. Jónsson, "Improved tangent estimate in the nudged elastic band method," J. Chem. Phys. 113, 9978 (2000). ——上流差分による接線の改良．
G. Henkelman, B. P. Uberuaga, and H. Jónsson, "A climbing image nudged elastic band method for finding saddle points and minimum energy paths," J. Chem. Phys. 113, 9901 (2000). ——CI-NEBの原論文．
G. H. Vineyard, "Frequency factors and isotope effects in solid state rate processes," J. Phys. Chem. Solids 3, 121 (1957). ——調和遷移状態理論の前因子．

分子動力学と温度制御

L. Verlet, "Computer 'experiments' on classical fluids," Phys. Rev. 159, 98 (1967).
W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, J. Chem. Phys. 76, 637 (1982). ——速度Verlet法．
S. Nosé, "A unified formulation of the constant temperature molecular dynamics methods," J. Chem. Phys. 81, 511 (1984);Mol. Phys. 52, 255 (1984).
W. G. Hoover, "Canonical dynamics: Equilibrium phase-space distributions," Phys. Rev. A 31, 1695 (1985). ——実時間形式への書き換えと非エルゴード性の指摘．
G. J. Martyna, M. L. Klein, and M. Tuckerman, "Nosé–Hoover chains: The canonical ensemble via continuous dynamics," J. Chem. Phys. 97, 2635 (1992).
M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford University Press, 2010). ——拡張系の統計力学とシンプレクティック積分の詳細．
A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, "Time-reversible Born–Oppenheimer molecular dynamics," Phys. Rev. Lett. 97, 123001 (2006). ——密度外挿によるエネルギードリフトの解決．

Car–Parrinello法とマルチスケール

R. Car and M. Parrinello, "Unified approach for molecular dynamics and density-functional theory," Phys. Rev. Lett. 55, 2471 (1985).
G. Pastore, E. Smargiassi, and F. Buda, "Theory of ab initio molecular-dynamics calculations," Phys. Rev. A 44, 6334 (1991). ——断熱条件と仮想振動数の解析．
D. Marx and J. Hutter, Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods (Cambridge University Press, 2009). ——BOMDとCPMDを統一的に扱う標準的な教科書．
A. Laio and M. Parrinello, "Escaping free-energy minima," Proc. Natl. Acad. Sci. USA 99, 12562 (2002). ——メタダイナミクスの原論文．
A. Barducci, G. Bussi, and M. Parrinello, "Well-tempered metadynamics," Phys. Rev. Lett. 100, 020603 (2008).
A. Warshel and M. Levitt, "Theoretical studies of enzymic reactions," J. Mol. Biol. 103, 227 (1976). ——QM/MM法の出発点．
H. M. Senn and W. Thiel, "QM/MM methods for biomolecular systems," Angew. Chem. Int. Ed. 48, 1198 (2009). ——リンク原子・埋め込み法の実務的な総説．