第8章Born–Oppenheimer近似と構造最適化
第7章までに，変分原理と摂動論を武器にして，水素原子やヘリウム原子の「エネルギー」を計算してきた．しかし，そこで求めた数値はいったい何のエネルギーだったのだろうか．電子だけのエネルギーなのか，原子核も含めた系全体のエネルギーなのか．ここをあいまいにしたままでは，これから扱う分子や固体の計算——結合長を求める，格子定数を求める，どちらの結晶構造が安定かを判定する——を正しく解釈できない．本章ではまず「全エネルギー」の中身を五つの項に分解して定義し，そのうち原子核の運動エネルギーだけを捨てるという大胆な近似，すなわちBorn–Oppenheimer近似を導入する．この近似は現代の量子化学計算・密度汎関数理論計算のほぼすべてが立脚する大前提であり，これによって初めて「ポテンシャルエネルギー曲面」という概念が生まれ，「原子に働く力」が定義でき，「構造最適化」というアルゴリズムが可能になる．後半では，実際に酸素分子 $\mathrm{O}_2$ の結合長を $0.9\ \text{Å}$ から $2.15\ \text{Å}$ まで変えながら全エネルギーを計算し，教科書に出てくるあの原子間ポテンシャル曲線を自分の手で描く．さらに Hellmann–Feynman の定理を完全に証明し，力の計算から構造緩和までを一本の線でつなぐ．
- 分子系の全ハミルトニアンが五つの項($\hat{T}_{\mathrm{Nu}},\ \hat{T}_{\mathrm{el}},\ V_{\mathrm{Nu\text{-}el}},\ V_{\mathrm{Nu\text{-}Nu}},\ V_{\mathrm{el\text{-}el}}$)からなること
- 陽子・中性子は電子の約1836倍重いという事実から，$\Psi(\bm{r},\bm{R})\simeq\Psi_{\mathrm{el}}(\bm{r},\bm{R})\Psi_{\mathrm{Nu}}(\bm{R})$ という積の近似が導かれること
- 核座標 $\bm{R}$ が電子波動関数に変数ではなくパラメータとして入ること，そのため完全な変数分離はできないこと
- 全エネルギー $E\simeq\min_{\bm{R}}E_{\mathrm{el}}(\bm{R})$ が絶対零度($0\ \mathrm{K}$)の内部エネルギーであり，有限温度ではエントロピーが要ること
- 断熱近似=「イオンは電子基底状態のポテンシャルエネルギー曲面上のみを動く」という言い換えと，2段階の最小化
- $\mathrm{O}_2$ の全エネルギー曲線を実際に計算する手順(xyzファイル・
sed・forループ・grep・awk・matplotlib) - 三重項酸素がなぜ基底状態か，なぜ液体酸素が磁石に引かれるか
- Hellmann–Feynmanの定理 $\bm{F}_I=-\bra{\Psi_{\mathrm{el}}}\nabla_I(\Ham_{\mathrm{el}}+V_{\mathrm{Nu\text{-}Nu}})\ket{\Psi_{\mathrm{el}}}$ の完全な証明とPulay力への注意
- Hessian行列・力の定数・振動数，そして最急降下法と準Newton法(BFGS)による構造緩和
8.1 「全エネルギー」とは何のエネルギーか
第2章から第7章まで，われわれは「エネルギー」という言葉をあまり吟味せずに使ってきた．水素原子の $1s$ 準位 $\epsilon_{1s}=-13.606\ \mathrm{eV}$，ヘリウム原子の基底状態エネルギー，Gauss基底で求めた近似値——これらは何のエネルギーだったのか．電子の内部エネルギーだろうか．それとも原子核まで含めた系全体の自由エネルギーだろうか．この問いにきちんと答えることが，本章の出発点である．
8.1.1 原子の場合:三つの項
原子番号 $Z$ の原子核ひとつのまわりに $n$ 個の電子がある系を考えよう．核は原点に固定してあるとする．この系のエネルギーに寄与する項は三つである．
第一に，電子の運動エネルギー(kinetic energy of electrons)．1個の電子の運動エネルギー演算子は $\hat{\bm{p}}^2/2m=-\hbar^2\nabla^2/2m$ だから，$n$ 個ぶんの和をとって
$$ \hat{T}_{\mathrm{el}}=\sum_{i=1}^{n}\left(-\frac{\hbar^2}{2m}\nabla_i^2\right). $$第二に，電子が原子核から受けるポテンシャルエネルギー．電荷 $+Z\ee$ の核と電荷 $-\ee$ の電子の間のCoulombエネルギーは負(引力)で，
$$ V_{\mathrm{Nu\text{-}el}}=-\sum_{i=1}^{n}\frac{Z\ee^{2}}{4\pi\eps\, r_i},\qquad r_i=\abs{\bm{r}_i}. $$第三に，電子同士の反発．これは正(斥力)で，
$$ V_{\mathrm{el\text{-}el}}=\frac{1}{2}\sum_{i=1}^{n}\sum_{j\neq i}\frac{\ee^{2}}{4\pi\eps\abs{\bm{r}_i-\bm{r}_j}} . $$三つを足したものが原子の全ハミルトニアンである:
$$ \begin{equation} \Ham_{\text{原子}} =\sum_{i=1}^{n}\left(-\frac{\hbar^2}{2m}\nabla_i^2\right) -\sum_{i=1}^{n}\frac{Z\ee^{2}}{4\pi\eps\, r_i} +\frac{1}{2}\sum_{i=1}^{n}\sum_{j\neq i}\frac{\ee^{2}}{4\pi\eps\abs{\bm{r}_i-\bm{r}_j}} . \label{eq:8-atom-hamiltonian} \end{equation} $$数学ノート:二重和の $\tfrac12$ と $j\neq i$
電子間反発の項でしばしば $\sum_j\sum_i$ とだけ書かれるが，これは記法の省略である．正しくは同じ電子どうしのペア($i=j$)を除き($\abs{\bm{r}_i-\bm{r}_i}=0$ で発散してしまう)，さらに $(i,j)$ と $(j,i)$ を二重に数えないように $\tfrac12$ を掛ける．$\tfrac12\sum_i\sum_{j\neq i}$ と $\sum_{i\lt j}$ は同じものである．$n=3$ なら和の項数は $\tfrac12\times3\times2=3$ 個，すなわち $(1,2),(1,3),(2,3)$ の3ペアで，確かに一致する．本書ではこの流儀に統一する．
8.1.2 分子の場合:核と核の反発が加わる
式 \eqref{eq:8-atom-hamiltonian} は原子核が1個の場合であった．原子核が2個以上になると，核どうしのCoulomb反発が新たに現れる．核 $I$ の電荷は $+Z_I\ee$，位置は $\bm{R}_I$ だから，
$$ V_{\mathrm{Nu\text{-}Nu}}=\frac{1}{2}\sum_{I=1}^{N}\sum_{J\neq I}\frac{Z_IZ_J\ee^{2}}{4\pi\eps\abs{\bm{R}_I-\bm{R}_J}} . $$この項は電子座標をまったく含まない．核の位置 $\bm{R}=\{\bm{R}_1,\dots,\bm{R}_N\}$ さえ決めてしまえば，それはただの数(定数)である．あとで見るように，この「電子から見れば定数」という性質が本章の議論全体の鍵になる．
さらに，核が動くならその運動エネルギーも要るはずである．核 $I$ の質量を $M_I$，運動量演算子を $\hat{\bm{P}}_I=-i\hbar\nabla_I$ とすれば
$$ \hat{T}_{\mathrm{Nu}}=\sum_{I=1}^{N}\frac{\hat{\bm{P}}_I^{2}}{2M_I}=\sum_{I=1}^{N}\left(-\frac{\hbar^{2}}{2M_I}\nabla_I^{2}\right). $$以上をすべて集めると，$N$ 個の原子核と $n$ 個の電子からなる系の厳密なハミルトニアンが得られる．これが本章の出発点であり，以後何度も参照する式である．
五つの項の意味を図8.1(a)に描いた．ここで注意してほしいのは，この式そのものには何の近似も入っていないことである．相対論効果とスピン–軌道相互作用を無視した非相対論的な扱いという限定はあるが，その枠内では式 \eqref{eq:8-full-hamiltonian} は厳密である．にもかかわらず，この式をそのまま解いた計算というものを，われわれは第7章までに一度も行っていない．核の運動エネルギー $\hat{T}_{\mathrm{Nu}}$ が，いつのまにか消えていたのである．
ん？ なぜ原子核の運動エネルギーは考えなくてよいのだろうか．
例題8.1 水の分子でハミルトニアンの項数を数える
水分子 $\mathrm{H_2O}$ について，式 \eqref{eq:8-full-hamiltonian} の五つの項がそれぞれ何個の和からなるか数えよ．また，波動関数 $\Psi$ は何個の実変数の関数か(スピンは除く)．
解答:原子核は O が1個と H が2個で $N=3$．電子数は $8+1+1=10$ で $n=10$．
- $\hat{T}_{\mathrm{Nu}}$:核ごとに1項なので $N=3$ 項．
- $\hat{T}_{\mathrm{el}}$:電子ごとに1項なので $n=10$ 項．
- $V_{\mathrm{Nu\text{-}el}}$:核と電子のすべての組で $N\times n=3\times10=30$ 項．
- $V_{\mathrm{Nu\text{-}Nu}}$:核のペアの数で $\dbinom{3}{2}=3$ 項．
- $V_{\mathrm{el\text{-}el}}$:電子のペアの数で $\dbinom{10}{2}=\dfrac{10\times9}{2}=45$ 項．
合計 $3+10+30+3+45=91$ 項である．変数の数は，核が $3\times3=9$ 個，電子が $3\times10=30$ 個で，合わせて $39$ 個．たった3原子の分子でも39次元空間の偏微分方程式であり，まともに解けるはずがない．だからこそ，変数を減らす近似が必要になる．
8.2 Born–Oppenheimer近似
8.2.1 核は重い，という一言から
原子核をつくる陽子と中性子は，電子に比べて桁違いに重い．正確には
$$ \begin{equation} \frac{m_{\mathrm{p}}}{m}=\frac{1.6726\times10^{-27}\ \mathrm{kg}}{9.1094\times10^{-31}\ \mathrm{kg}}=1836.15\simeq1840 \label{eq:8-mass-ratio} \end{equation} $$である．式 \eqref{eq:8-mass-ratio} の比は，質量数 $A$ の原子核ならさらに $A$ 倍になる．たとえば酸素($A=16$)の核は電子の約 $2.9\times10^{4}$ 倍も重い．
なぜ？：重いものは「止まっている」とみなしてよいのか
質量が大きいというだけでは「動かない」ことにはならない．重要なのは速さである．核と電子はCoulomb力でつながっているから，やりとりする運動量 $P$ は同程度と考えてよい．すると速さは $v=P/M$ だから，質量に反比例して
$$\frac{v_{\text{核}}}{v_{\text{電子}}}\sim\frac{m}{M}\sim\frac{1}{1836}$$となる．電子が原子のまわりを1周する間に，核は $1/1836$ 周ぶんしか動かない．つまり電子の立場から見れば，核は事実上その場に釘づけになっている．逆に核の立場から見ると，電子は速すぎて個々の位置を追えず，なめらかに広がった「電子雲」としてしか感じられない．この二つの見え方を別々の問題に分けてしまおう，というのがBorn–Oppenheimer近似の発想である．
例題8.2 時間スケールを実際に見積もる
(1) 水素原子の $1s$ 電子の速さは $v=\alpha c=2.19\times10^{6}\ \mathrm{m/s}$($\alpha$ は微細構造定数)である．Bohr半径 $a_0=0.529\times10^{-10}\ \mathrm{m}$ の円周を1周するのに要する時間を求めよ．
(2) 酸素分子の振動数は実測で $\tilde{\nu}=1580\ \mathrm{cm^{-1}}$ である．振動の周期を求め，(1)と比べよ．
解答:(1) 周回時間は
$$T_{\mathrm{el}}=\frac{2\pi a_0}{v}=\frac{2\pi\times0.529\times10^{-10}}{2.19\times10^{6}}=1.52\times10^{-16}\ \mathrm{s}=0.152\ \mathrm{fs}.$$(2) 波数 $\tilde\nu$ の振動数は $\nu=c\tilde\nu=2.998\times10^{10}\ \mathrm{cm/s}\times1580\ \mathrm{cm^{-1}}=4.74\times10^{13}\ \mathrm{Hz}$ だから，周期は
$$T_{\mathrm{vib}}=\frac{1}{\nu}=2.11\times10^{-14}\ \mathrm{s}=21.1\ \mathrm{fs}.$$比は $T_{\mathrm{vib}}/T_{\mathrm{el}}\simeq139$．核が1回伸び縮みするあいだに，電子は約140周もしている．電子は核の動きに瞬時に追随できると考えてよい．
8.2.2 波動関数を積の形に分ける
いま述べた描像を数式にしよう．核の座標をまとめて $\bm{R}=\{\bm{R}_1,\bm{R}_2,\dots,\bm{R}_N\}$，電子の座標をまとめて $\bm{r}=\{\bm{r}_1,\bm{r}_2,\dots,\bm{r}_n\}$ と書く．系全体の波動関数 $\Psi(\bm{r},\bm{R})$ を，電子の波動関数と核の波動関数の積で近似する．
定義:Born–Oppenheimer近似
核の波動関数を $\Psi_{\mathrm{Nu}}(\bm{R})$，電子の波動関数を $\Psi_{\mathrm{el}}(\bm{r},\bm{R})$ とするとき，系全体の波動関数を
$$ \begin{equation} \Psi(\bm{r},\bm{R})\ \simeq\ \Psi_{\mathrm{el}}(\bm{r},\bm{R})\,\Psi_{\mathrm{Nu}}(\bm{R}) \label{eq:8-bo-product} \end{equation} $$と積の形に書く近似をBorn–Oppenheimer近似(Born–Oppenheimer approximation，BO近似)という．M. Born と J. R. Oppenheimer が1927年に提出した(文献[4])．
この式でいちばん大切なのは，$\Psi_{\mathrm{el}}$ が $\bm{r}$ だけでなく $\bm{R}$ にも依存していることである．核をどこに置くかによって電子の感じるポテンシャルが変わるのだから，電子状態が核配置に依存するのは当たり前である．しかしその依存のしかたが，$\bm{r}$ の場合とはまったく違う．
注意:「変数」と「パラメータ」を混同しないこと
$\Psi_{\mathrm{el}}(\bm{r},\bm{R})$ において，$\bm{r}$ は変数だが $\bm{R}$ はパラメータである．これは記法上の区別ではなく，次の二つの意味で本質的である．
- 電子のSchrödinger方程式を解くとき，微分は $\bm{r}$ についてだけ行い，$\bm{R}$ は定数として扱う．
- 規格化は $\bm{r}$ についてだけ行う:$\displaystyle\int\abs{\Psi_{\mathrm{el}}(\bm{r},\bm{R})}^{2}\dd\bm{r}=1$ がすべての $\bm{R}$ について成り立つ．
$\bm{R}$ の値を一つ決めるごとに，電子の問題が一つ立ち上がって一つの答えが出る．$\bm{R}$ を少しずつ変えれば，答えの列 $E_{\mathrm{el}}(\bm{R})$ が得られる．この「$\bm{R}$ をずらしながら何度も解く」というのが，8.5節で $\mathrm{O}_2$ に対して実際にやることそのものである．
8.3 数式による定式化
ここからは式 \eqref{eq:8-bo-product} を式 \eqref{eq:8-full-hamiltonian} に代入して，何が起こるかを一段ずつ確かめる．計算は長いが，一つひとつは高校数学の積の微分にすぎない．
8.3.1 記号を短くする
毎回すべての和を書くのは煩雑なので，式 \eqref{eq:8-full-hamiltonian} の五つの項をそれぞれ一つの記号で表す．Schrödinger方程式 $\Ham\Psi=E\Psi$ は
$$ \begin{equation} \Bigl[\hat{T}_{\mathrm{Nu}}+\hat{T}_{\mathrm{el}}+V_{\mathrm{Nu\text{-}Nu}}(\bm{R})+V_{\mathrm{el\text{-}el}}(\bm{r})+V_{\mathrm{Nu\text{-}el}}(\bm{r},\bm{R})\Bigr]\Psi(\bm{r},\bm{R})=E\,\Psi(\bm{r},\bm{R}) \label{eq:8-shorthand} \end{equation} $$と書ける．括弧内の各項が何の変数に依存するかに注目してほしい．$V_{\mathrm{Nu\text{-}Nu}}$ は $\bm{R}$ だけ，$V_{\mathrm{el\text{-}el}}$ は $\bm{r}$ だけ，そして $V_{\mathrm{Nu\text{-}el}}$ だけが $\bm{r}$ と $\bm{R}$ の両方を含む．この最後の項が，これから見る「完全には変数分離できない」という事情の元凶である．
8.3.2 積を代入して，核の運動エネルギーを展開する
$\Psi=\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}}$ を代入する．ポテンシャル項は掛け算するだけなので何も起こらない．$\hat{T}_{\mathrm{el}}$ は $\bm{r}$ について微分するので $\Psi_{\mathrm{Nu}}(\bm{R})$ を素通りする:
$$ \hat{T}_{\mathrm{el}}\bigl[\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}}\bigr] =\Psi_{\mathrm{Nu}}(\bm{R})\,\hat{T}_{\mathrm{el}}\Psi_{\mathrm{el}}(\bm{r},\bm{R}). $$問題は $\hat{T}_{\mathrm{Nu}}$ である．$\nabla_I$ は核座標 $\bm{R}_I$ についての微分だから，$\Psi_{\mathrm{el}}$ も $\Psi_{\mathrm{Nu}}$ も両方とも微分されてしまう．積の微分法則を2回使うと
$$ \begin{equation} \nabla_I^{2}\bigl[\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}}\bigr] =\Psi_{\mathrm{el}}\,\nabla_I^{2}\Psi_{\mathrm{Nu}} +2\bigl(\nabla_I\Psi_{\mathrm{el}}\bigr)\!\cdot\!\bigl(\nabla_I\Psi_{\mathrm{Nu}}\bigr) +\Psi_{\mathrm{Nu}}\,\nabla_I^{2}\Psi_{\mathrm{el}} \label{eq:8-laplacian-product} \end{equation} $$となる．念のため1次元で確かめておこう．$f(R)g(R)$ に対して $(fg)'=f'g+fg'$，もう一度微分して $(fg)''=f''g+2f'g'+fg''$．式 \eqref{eq:8-laplacian-product} は $f=\Psi_{\mathrm{Nu}}$，$g=\Psi_{\mathrm{el}}$ としたものにほかならない．
したがって $\hat{T}_{\mathrm{Nu}}$ の作用は
$$ \hat{T}_{\mathrm{Nu}}\bigl[\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}}\bigr] =\Psi_{\mathrm{el}}\bigl[\hat{T}_{\mathrm{Nu}}\Psi_{\mathrm{Nu}}\bigr] \;\underbrace{-\sum_{I}\frac{\hbar^{2}}{M_I}\bigl(\nabla_I\Psi_{\mathrm{el}}\bigr)\!\cdot\!\bigl(\nabla_I\Psi_{\mathrm{Nu}}\bigr) -\sum_{I}\frac{\hbar^{2}}{2M_I}\Psi_{\mathrm{Nu}}\nabla_I^{2}\Psi_{\mathrm{el}}}_{\displaystyle \equiv\ \hat{\Lambda}\ \text{(非断熱項)}} $$と書ける．第2項と第3項をまとめて非断熱項(nonadiabatic coupling terms)あるいは微分結合項という．
なぜ？：非断熱項を捨ててよい理由
非断熱項 $\hat{\Lambda}$ には，どちらの項にも $1/M_I$ が掛かっている．核質量は電子質量の数千倍から数十万倍だから，これだけで係数が小さい．それに加えて，$\nabla_I\Psi_{\mathrm{el}}$ と $\nabla_I^{2}\Psi_{\mathrm{el}}$ は「核をわずかに動かしたとき電子波動関数がどれだけ変わるか」を表す量である．8.2節で見たように電子は核の動きに瞬時に追随して，なめらかに形を変えるだけだから，これらは大きな量にはならない．小さい係数 × 小さい量なので，$\hat{\Lambda}$ 全体を落としてよい，というのがBO近似の中身である．
逆に言えば，二つの電子状態のエネルギーが接近する場所(円錐交差，conical intersection)では $\Psi_{\mathrm{el}}$ が $\bm{R}$ のわずかな変化で激変するので，$\nabla_I\Psi_{\mathrm{el}}$ が発散的に大きくなり，BO近似は破綻する．光化学反応や無輻射失活はまさにその舞台である．
$\hat{\Lambda}$ を捨てると，式 \eqref{eq:8-shorthand} は
$$ \begin{equation} \bigl[\hat{T}_{\mathrm{Nu}}\Psi_{\mathrm{Nu}}(\bm{R})\bigr]\Psi_{\mathrm{el}}(\bm{r},\bm{R}) +\Bigl[\hat{T}_{\mathrm{el}}\Psi_{\mathrm{el}} +V_{\mathrm{Nu\text{-}el}}\Psi_{\mathrm{el}} +V_{\mathrm{Nu\text{-}Nu}}\Psi_{\mathrm{el}} +V_{\mathrm{el\text{-}el}}\Psi_{\mathrm{el}}\Bigr]\Psi_{\mathrm{Nu}}(\bm{R}) =E\,\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}} \label{eq:8-after-sub} \end{equation} $$となる．
8.3.3 両辺を $\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}}$ で割る
変数分離の常套手段どおり，式 \eqref{eq:8-after-sub} の両辺を $\Psi(\bm{r},\bm{R})=\Psi_{\mathrm{el}}\Psi_{\mathrm{Nu}}$ で割ってみる．第1項は $\Psi_{\mathrm{el}}$ が約分され，第2項は $\Psi_{\mathrm{Nu}}$ が約分されて
を得る．ここで第7章のヘリウム原子のときと同じ状況に出会う．この式は完全な変数分離形になっていない．第1項は $\bm{R}$ だけの関数だからよいが，第2項は $\bm{r}$ と $\bm{R}$ の両方を含んでいる．もし第2項が $\bm{r}$ だけの関数なら「$\bm{R}$ の関数 + $\bm{r}$ の関数 = 定数」となって，それぞれが定数に等しいと結論できたのだが，そうはならない．核由来のエネルギーは核の位置だけで決まるのに，電子由来のエネルギーは電子の位置だけでなく核の位置にも依存してしまう——これが分離を妨げている．
8.3.4 二つの方程式に分ける
それでもBO近似は使える．そのからくりは，$\bm{R}$ を固定した電子の固有値問題を先に解いてしまうことにある．$\bm{R}$ を止めたときの電子のSchrödinger方程式を
$$ \begin{equation} \Bigl[\hat{T}_{\mathrm{el}}+V_{\mathrm{Nu\text{-}el}}(\bm{r},\bm{R})+V_{\mathrm{el\text{-}el}}(\bm{r})+V_{\mathrm{Nu\text{-}Nu}}(\bm{R})\Bigr]\Psi_{\mathrm{el}}(\bm{r},\bm{R}) =E_{\mathrm{el}}(\bm{R})\,\Psi_{\mathrm{el}}(\bm{r},\bm{R}) \label{eq:8-el-eigen} \end{equation} $$と書く．左辺の演算子は，$\bm{R}$ を定数とみなせば $\bm{r}$ だけの演算子である．だからこれは(第7章までに扱ってきたのとまったく同じ形の)電子だけの固有値問題であり，その固有値 $E_{\mathrm{el}}(\bm{R})$ は$\bm{r}$ を含まない，$\bm{R}$ だけの関数になる．積分してしまえば電子座標は消えるからである．
これを式 \eqref{eq:8-divide} に戻すと第2項は $E_{\mathrm{el}}(\bm{R})$ そのものになり，両辺に $\Psi_{\mathrm{Nu}}$ を掛け直して
$$ \begin{equation} \Bigl[\hat{T}_{\mathrm{Nu}}+E_{\mathrm{el}}(\bm{R})\Bigr]\Psi_{\mathrm{Nu}}(\bm{R})=E\,\Psi_{\mathrm{Nu}}(\bm{R}) \label{eq:8-nuc-eigen} \end{equation} $$という核だけのSchrödinger方程式が得られる．注目すべきは，電子系のエネルギー $E_{\mathrm{el}}(\bm{R})$ が，核にとってのポテンシャルエネルギーの役割を果たしていることである．電子の詳細はすべて $E_{\mathrm{el}}(\bm{R})$ という一つの関数に押し込められた．分子振動論も，格子振動(フォノン)も，分子動力学法も，すべてこの式 \eqref{eq:8-nuc-eigen} から出発する．
8.3.5 全エネルギーとは何か
さて，量子化学計算やDFT計算が出力する「Total Energy(全エネルギー)」とは，式 \eqref{eq:8-nuc-eigen} の $E$ ではない．実際の計算では核の運動エネルギー $\hat{T}_{\mathrm{Nu}}$ をさらに $0$ とおいてしまう．すると $E=E_{\mathrm{el}}(\bm{R})$ となり，あとは核配置 $\bm{R}$ を動かしてこれを最小にすればよい．
定義:全エネルギー(total energy)
核の運動エネルギーを $0$ とみなしたうえで，電子系のエネルギーを変分原理によって最小化し，さらに核配置についても最小化して得られる値
$$ \begin{equation} E\ \simeq\ \min_{\bm{R}}\ E_{\mathrm{el}}(\bm{R}) \qquad\text{ただし}\quad E_{\mathrm{el}}(\bm{R})=\min_{\Psi_{\mathrm{el}}}\frac{\bra{\Psi_{\mathrm{el}}}\Ham_{\mathrm{el}}+V_{\mathrm{Nu\text{-}Nu}}\ket{\Psi_{\mathrm{el}}}}{\braket{\Psi_{\mathrm{el}}|\Psi_{\mathrm{el}}}} \label{eq:8-total-energy} \end{equation} $$を全エネルギー(total energy)という．内側の最小化は第3章の変分原理そのものである．
核の運動エネルギーをゼロとおくということは，核が熱運動もせず，零点振動すらしていないということである．つまり全エネルギーは絶対零度($0\ \mathrm{K}$)における系の内部エネルギーである．ここは学生がいちばん誤解しやすいところなので，強調しておきたい．
注意:全エネルギーは $0\ \mathrm{K}$ の量である
第一原理計算で「この構造のほうが $0.1\ \mathrm{eV}$ 安定だ」と言うとき，それは $0\ \mathrm{K}$・零点振動なしという条件下での比較である．有限温度で物質がどちらの相を選ぶかを論じるには，自由エネルギー
$$F=E+E_{\mathrm{ZPE}}-TS$$を比べなければならない．$E_{\mathrm{ZPE}}$ は零点振動エネルギー(8.7節のHessianから求まる)，$S$ には振動エントロピー・配置エントロピー・電子エントロピー・磁気エントロピーなどが含まれる．室温 $300\ \mathrm{K}$ では $k_{\mathrm{B}}T=0.0259\ \mathrm{eV}$ であり，$TS$ 項が $0.1\ \mathrm{eV}$ 程度の差をひっくり返すことは珍しくない．全エネルギーだけで相の安定性を断定してはいけない．
物理的意味:なぜBO近似はここまで成功しているのか
量子化学計算も密度汎関数理論計算も，ほぼ例外なくBO近似を大前提としている．「分子の形」「結合長」「格子定数」「フォノン分散」「反応経路」といった，われわれが物質を語るときに使うほとんどすべての概念は，$E_{\mathrm{el}}(\bm{R})$ という関数が存在することを前提にしている．BO近似がなければ，そもそも「分子の構造」という言葉すら定義できない——核の位置は電子と量子力学的にもつれ合ったままで，「原子がここにある」と言えなくなってしまうからである．BO近似は単なる計算の便法ではなく，化学の言葉そのものを可能にしている近似なのである．
8.4 断熱近似とポテンシャルエネルギー曲面
8.4.1 言い換えとしての断熱近似
前節の結論を，核の側から言い換えてみよう．核が従う方程式 \eqref{eq:8-nuc-eigen} を見ると，核は $E_{\mathrm{el}}(\bm{R})$ というポテンシャルの中を動く粒子にほかならない．しかもここで使った $E_{\mathrm{el}}(\bm{R})$ は，各 $\bm{R}$ における電子系の基底状態のエネルギーである．核がどれほどゆっくり動こうとも，電子はつねに基底状態に居続けて，励起状態に飛び移ったりはしない．
定義:断熱近似
断熱近似(adiabatic approximation)とは，イオン(原子核)は電子の基底状態のポテンシャルエネルギー曲面上のみを動くと考える近似のことである．BO近似は断熱近似とも呼ばれる．厳密には両者はイコールではない(BO近似では非断熱項 $\hat\Lambda$ を全部落とすが，断熱近似では対角成分 $\bra{\Psi_{\mathrm{el}}}\hat\Lambda\ket{\Psi_{\mathrm{el}}}$ を残す流儀もある)が，本質的には同じものと考えてよい．
定義:ポテンシャルエネルギー曲面
核配置 $\bm{R}$ に対して電子基底状態のエネルギー $E_{\mathrm{el}}(\bm{R})$ を対応させる関数を，ポテンシャルエネルギー曲面(potential energy surface, PES)という．$N$ 原子分子では，並進3自由度と回転3自由度(直線分子は2)を除いて $3N-6$ 次元(直線分子は $3N-5$ 次元)の曲面である．二原子分子なら $3\times2-5=1$ 次元，つまり結合長 $R$ 一つの関数になる．8.5節で描く曲線がまさにこれである．
8.4.2 2段階の最小化
全エネルギーを求める手続きは，式 \eqref{eq:8-total-energy} が示すとおり二重の最小化である．
- 内側(電子の問題):核を $\bm{R}$ に固定したまま，変分原理によって電子系のエネルギー $E_{\mathrm{el}}(\bm{R})$ を最小化する．この段階では核間反発 $V_{\mathrm{Nu\text{-}Nu}}(\bm{R})$ は定数なので，まだ下げようがない．
- 外側(核の問題):核を動かして，$E_{\mathrm{el}}(\bm{R})$ そのものを最小にする $\bm{R}$ を探す．これが構造緩和である．
この入れ子構造を図8.2に示した．上段が各 $\bm{R}$ における電子の問題，下段がその答えを並べて得られるPESと，その上を転がる核である．
PESの上には，最小点(安定構造)のほかにも，局所的な窪み(準安定構造)や，二つの窪みを結ぶ峠(鞍点，saddle point)がある．鞍点は化学反応の遷移状態に対応し，その高さが活性化エネルギーを与える．第13章・第14章で扱う固体の構造安定性の議論も，すべてこのPESの上での話である．担当教員の研究でも，$\mathrm{SrTiO_3}$ 表面の多数の候補構造についてPESの局所最小点を網羅的に探索している(文献[12])．
8.5 実践:O₂の結合長を系統的に変える
理屈はここまでにして，実際に手を動かそう．酸素分子 $\mathrm{O}_2$ の結合長を $0.9\ \text{Å}$ から $2.15\ \text{Å}$ まで $0.05\ \text{Å}$ 刻みで変え，それぞれについて全エネルギーを計算して，ポテンシャルエネルギー曲線を描く．使うのは PySCF(文献[8])と，Unix系OSの標準的なコマンドだけである．作業ディレクトリの準備は付録Bを参照してほしい．
8.5.1 構造ファイルを作る
分子の構造はxyzファイルという単純な書式で表す．1行目に原子数，2行目にコメント，3行目以降に「元素記号 $x$ $y$ $z$」を Å 単位で並べるだけである．cat コマンドでファイルの中身を表示してみよう．
$ cat XYZ_O2.xyz
2
Converted from SDF_O2.sdf (Angstrom)
O 0.00000000 0.00000000 0.00000000
O 1.20000000 0.00000000 0.00000000コード8.1 cat でxyzファイルを表示する．1個目の O が原点，2個目が $x$ 軸上 $1.2\ \text{Å}$ にある．
結合長を変えるには，1.20000000 という文字列を別の数値に置き換えればよい．文字列置換には sed(stream editor)を使う．書式は
sed -e 's/置換前/置換後/g' 入力ファイル名 > 保存ファイル名コード8.2 sed の基本形．s は substitute(置換)，g は global(行内のすべて)を意味する．> は出力をファイルに書き出すリダイレクトである．
$ sed -e 's/1.20000000/2.0000000/g' XYZ_O2.xyz > XYZ_O2_2p0.xyz
$ cat XYZ_O2_2p0.xyz
2
Converted from SDF_O2.sdf (Angstrom)
O 0.00000000 0.00000000 0.00000000
O 2.0000000 0.00000000 0.00000000コード8.3 結合長 $2.0\ \text{Å}$ の構造ファイルができた．ファイル名の 2p0 の p は小数点 . の代わりである(ファイル名に . を入れると拡張子と紛らわしいため)．
8.5.2 1点ずつ計算してみる
まずは2通りだけ計算する．付録Bで用意した calc_pyscf.py は，xyzファイル・スピン多重度・基底関数を指定して SCF 計算を実行するスクリプトである．
$ python3 calc_pyscf.py --xyz XYZ_O2.xyz --spin 2 --basis 6-31g
========================================
PySCF Calculation
========================================
XYZ file : XYZ_O2.xyz
Basis : 6-31g
Charge : 0
Spin (2S) : 2
Theory : scf
----------------------------------------
Method : UHF
Total Energy : -4069.348307 eV
<S^2> : 2.033075
Multiplicity : 3.0コード8.4 結合長 $1.2\ \text{Å}$ の $\mathrm{O}_2$(高スピン)の計算．全エネルギーは $-4069.348\ \mathrm{eV}$．
同じことを $2.0\ \text{Å}$ の構造に対して行うと $-4061.466\ \mathrm{eV}$ になる．$1.2\ \text{Å}$ のほうが $7.88\ \mathrm{eV}$ も低い．--spin 2 の意味(なぜ $\mathrm{O}_2$ は $2S=2$ なのか)は8.6節で説明する．
注意:全エネルギーの絶対値そのものに意味はない
$-4069\ \mathrm{eV}$ という値は，$\mathrm{O}_2$ の16個の電子をすべて核から無限遠まで引き離すのに必要なエネルギーに相当する巨大な量である．一方，化学結合の強さは数 $\mathrm{eV}$，構造の違いによるエネルギー差は $0.01\ \mathrm{eV}$ 程度のこともある．すなわち意味があるのはつねに差だけであり，しかも $10^{-5}$ 以下の相対精度で差をとることになる．だからこそ，比較する二つの計算では基底関数・計算手法・収束条件をそろえなければならない．
8.5.3 ループ処理で26通りの構造を一気に作る
26通りを手で打つのは苦行である．「同じ作業のくり返しは機械にやらせるべき」だ．シェルの for ループを使う．
$ for l in {0..25} ; do
> bond_length=`printf "%.6f\n" $(( 0.9 + l * 0.05 ))`
> rename=`echo "$bond_length" | tr '.' 'p'`
> sed -e "s/1.20000000/$bond_length/g" XYZ_O2.xyz > XYZ_O2_"$rename".xyz
> doneコード8.5 結合長 $0.9,0.95,\dots,2.15\ \text{Å}$ の26個の構造ファイルを生成する．
読み方を1行ずつ確認しよう．
{0..25}は $0,1,2,\dots,25$ という数の列に展開される．変数lがこれらを順にとる．- ドルマーク
$は変数の値を取り出す記号である．$(( ... ))は算術演算を行う書き方で，0.9 + l * 0.05を計算する．printf "%.6f\n"で小数点以下6桁に整形する． tr '.' 'p'は文字.をpに置き換えるコマンド．$1.9 \to$1p900000となり，これをファイル名に使う．- 最後の
sedで，もとの1.20000000を新しい結合長に置換して保存する．
注意:シェルによって小数の扱いが違う
コード8.5 の $(( 0.9 + l * 0.05 )) は，macOS の既定シェルである zsh では小数として計算される．しかし Linux でよく使われる bash の算術展開は整数しか扱えず，この行はエラーになる．bash を使う場合は
bond_length=`awk -v l=$l 'BEGIN{printf "%.6f", 0.9 + l*0.05}'`のように awk や bc に計算させればよい．自分のシェルは echo $SHELL で確認できる．
実行後に ls すると，XYZ_O2_0p900000.xyz から XYZ_O2_2p150000.xyz まで26個のファイルができている．中身を確認すると，確かに座標が書き換わっている．
$ cat XYZ_O2_1p900000.xyz
2
Converted from SDF_O2.sdf (Angstrom)
O 0.00000000 0.00000000 0.00000000
O 1.900000 0.00000000 0.00000000コード8.6 $1.20000$ だったところが $1.900000$ に変わっている．
8.5.4 ループ処理で26通りの計算を回す
同じ要領で，26個の構造ファイルを順に PySCF に流し込む．出力は output_0p900000 のようなファイルに保存する．
$ for l in {0..25} ; do
> bond_length=`printf "%.6f\n" $(( 0.9 + l * 0.05 ))`
> rename=`echo "$bond_length" | tr '.' 'p'`
> python3 calc_pyscf.py --xyz XYZ_O2_"$rename".xyz --spin 2 --basis 6-31g > output_"$rename"
> doneコード8.7 26通りの電子状態計算を一括実行する．1点あたり数秒なので，全体でも1分程度で終わる．
8.5.5 出力から必要な数値だけ抜き出す(grep と awk)
26個の出力ファイルから全エネルギーだけを拾い出したい．文字列を含む行を抜き出すコマンドが grep である．ファイル名の * は「任意の文字列」を意味するワイルドカードである．
$ grep "Total Energy" output_*
output_0p900000:Total Energy : -4060.401153 eV
output_0p950000:Total Energy : -4064.025729 eV
output_1p000000:Total Energy : -4066.428637 eV
(中略)
output_2p150000:Total Energy : -4060.057679 eVコード8.8 grep で必要な行だけを抜き出す．
数値だけがほしいので，さらに awk(オーク)につなぐ．awk は行をスペース区切りの「列」に分解し，$4 と書けば4列目を取り出す．上の行を数えると，1列目 Total，2列目 Energy，3列目 :，4列目が数値である．縦棒 | はパイプで，左のコマンドの出力を右のコマンドの入力につなぐ．
$ grep "Total Energy" output_* | awk '{print $4}'
-4060.401153
-4064.025729
-4066.428637
(中略)
-4060.057679コード8.9 awk で4列目の数値だけを取り出す．
プロットするには「1列目に結合長，2列目に全エネルギー」という2列のファイルがほしい．ここまでの道具を組み合わせればよい．
$ for l in {0..25} ; do
> bond_length=`printf "%.6f\n" $(( 0.9 + l * 0.05 ))`
> rename=`echo "$bond_length" | tr '.' 'p'`
> total_energy=`grep "Total Energy" output_$rename | awk '{print $4}'`
> echo $bond_length $total_energy
> done > energies_O2.txtコード8.10 2列のデータファイル energies_O2.txt をつくる．done の後ろの > でループ全体の出力をまとめてファイルに書き出している．
8.5.6 プロットする
Python に慣れているなら matplotlib，慣れていないなら Excel でよい．ここでは matplotlib のスクリプトを示す．
import numpy as np
import matplotlib.pyplot as plt
r, e = np.loadtxt("energies_O2.txt", unpack=True)
fig, ax = plt.subplots(figsize=(5.0, 3.6))
ax.plot(r, e, "o-", color="#14684a", markersize=4)
ax.set_xlabel("Bond length (Angstrom)")
ax.set_ylabel("Total energy (eV)")
ax.grid(alpha=0.3)
i = e.argmin()
ax.annotate("min: %.2f A" % r[i], xy=(r[i], e[i]),
xytext=(r[i] + 0.35, e[i] + 2.0),
arrowprops=dict(arrowstyle="->"))
fig.tight_layout()
fig.savefig("O2_curve.png", dpi=200)コード8.11 energies_O2.txt を読み込んでプロットする Python スクリプト．unpack=True で列ごとに配列に分けられる．
得られた数値を表8.1に，グラフを図8.3に示す．教科書によく出てくる，あの原子間ポテンシャルの曲線になっている．
| $R$ (Å) | 高スピン (eV) | 低スピン (eV) | $R$ (Å) | 高スピン (eV) | 低スピン (eV) |
|---|---|---|---|---|---|
| 0.90 | −4060.401153 | −4058.258513 | 1.55 | −4066.497268 | −4064.227915 |
| 0.95 | −4064.025729 | −4061.839113 | 1.60 | −4065.920407 | −4063.659827 |
| 1.00 | −4066.428637 | −4064.207275 | 1.65 | −4065.337272 | −4063.086134 |
| 1.05 | −4067.944972 | −4065.697019 | 1.70 | −4064.753701 | −4062.512530 |
| 1.10 | −4068.823058 | −4066.555558 | 1.75 | −4064.174595 | −4061.943789 |
| 1.15 | −4069.246134 | −4066.965069 | 1.80 | −4063.604065 | −4061.383916 |
| 1.20 | −4069.348307 | −4067.058685 | 1.85 | −4063.045509 | −4060.836211 |
| 1.25 | −4069.226505 | −4066.932459 | 1.90 | −4062.501654 | −4060.303309 |
| 1.30 | −4068.949681 | −4066.654577 | 1.95 | −4061.974597 | −4059.787227 |
| 1.35 | −4068.566101 | −4066.272645 | 2.00 | −4061.465863 | −4059.289414 |
| 1.40 | −4068.109170 | −4065.819518 | 2.05 | −4060.976468 | −4058.810817 |
| 1.45 | −4067.601987 | −4065.317846 | 2.10 | −4060.506997 | −4058.351961 |
| 1.50 | −4067.060767 | −4064.783486 | 2.15 | −4060.057679 | −4057.913023 |
物理的意味:曲線の形が語ること
曲線の左側(短距離側)が急峻に立ち上がるのは，核どうしのCoulomb反発と，Pauliの排他原理による内殻電子の反発(第7章)のためである．右側(長距離側)がゆるやかなのは，結合が伸びるにつれて結合性軌道の重なり(第5章)が指数関数的に失われていくためである．この非対称な井戸が，熱膨張という日常的な現象の起源になっている．温度が上がって振動振幅が大きくなると，井戸が非対称なぶん平均の核間距離が長いほうへずれるのである．
例題8.3 放物線内挿で平衡結合長・力の定数・振動数を求める
表8.1の高スピンの値のうち $R=1.15,\ 1.20,\ 1.25\ \text{Å}$ の3点を使い，(1) 平衡結合長 $R_{\mathrm{e}}$，(2) 力の定数 $k$，(3) 振動の波数 $\tilde\nu$ を求めよ．酸素原子の質量は $M=15.995\ u$，$1\,u=1.6605\times10^{-27}\ \mathrm{kg}$ とする．
解答:刻み幅を $h=0.05\ \text{Å}$，中央の点を $R_0=1.20\ \text{Å}$ とし，$E_-=E(R_0-h)$，$E_0=E(R_0)$，$E_+=E(R_0+h)$ と書く．3点を通る2次関数 $E(R)=a(R-R_0)^2+b(R-R_0)+E_0$ の係数は，$E_\pm=ah^2\pm bh+E_0$ から
$$ a=\frac{E_+-2E_0+E_-}{2h^{2}},\qquad b=\frac{E_+-E_-}{2h} $$と決まる．数値を入れると
$$ \begin{align} E_+-2E_0+E_-&=-4069.226505-2(-4069.348307)+(-4069.246134)=0.223975\ \mathrm{eV},\notag\\ E_+-E_-&=-4069.226505-(-4069.246134)=0.019629\ \mathrm{eV}.\notag \end{align} $$(1) 頂点の位置は $E'(R)=2a(R-R_0)+b=0$ より $R_{\mathrm{e}}=R_0-b/(2a)$．$b/(2a)=\dfrac{h}{2}\cdot\dfrac{E_+-E_-}{E_+-2E_0+E_-}$ だから
$$ \begin{equation} R_{\mathrm{e}}=R_0-\frac{h}{2}\cdot\frac{E_+-E_-}{E_+-2E_0+E_-} =1.20-\frac{0.05}{2}\times\frac{0.019629}{0.223975}=1.20-0.00219=1.1978\ \text{Å}. \label{eq:8-parabola-vertex} \end{equation} $$(2) 力の定数は $k=E''(R_{\mathrm{e}})=2a=(E_+-2E_0+E_-)/h^{2}$．単位をSIに直す．$h=0.05\ \text{Å}=5.0\times10^{-12}\ \mathrm{m}$，$1\ \mathrm{eV}=1.6022\times10^{-19}\ \mathrm{J}$ なので
$$ k=\frac{0.223975\times1.6022\times10^{-19}\ \mathrm{J}}{(5.0\times10^{-12}\ \mathrm{m})^{2}} =\frac{3.5885\times10^{-20}}{2.50\times10^{-23}}=1.435\times10^{3}\ \mathrm{N/m}. $$(3) 二原子分子の換算質量は $\mu=M_AM_B/(M_A+M_B)=M/2$ だから
$$ \mu=\frac{15.995\times1.6605\times10^{-27}}{2}=1.3280\times10^{-26}\ \mathrm{kg}. $$角振動数と波数は
$$ \omega=\sqrt{\frac{k}{\mu}}=\sqrt{\frac{1435.4}{1.3280\times10^{-26}}}=3.288\times10^{14}\ \mathrm{rad/s}, \qquad \tilde\nu=\frac{\omega}{2\pi c}=\frac{3.288\times10^{14}}{2\pi\times2.998\times10^{10}}=1.75\times10^{3}\ \mathrm{cm^{-1}} . $$実験値は $R_{\mathrm{e}}=1.2075\ \text{Å}$，$k=1177\ \mathrm{N/m}$，$\tilde\nu=1580\ \mathrm{cm^{-1}}$(文献[11])である．結合長は $0.8\%$ の誤差でよく合っているが，力の定数は $22\%$，振動数は $10\%$ 過大評価している．これは Hartree–Fock 法が電子相関を無視しているために結合が硬くなりすぎることによる(第9章・第10章)．
数学ノート:格子の粗さがもたらす誤差
例題8.3の放物線内挿の公式 \eqref{eq:8-parabola-vertex} は3階微分の効果を無視している．$h$ の2次まで展開すると，中心差分による1階微分には $\dfrac{h^{2}}{6}E'''$ の誤差が乗る．$\mathrm{O}_2$ の場合，表8.1から $E'''\simeq-6.6\times10^{2}\ \mathrm{eV/\text{Å}^{3}}$ と見積もられるので，この誤差は $0.28\ \mathrm{eV/\text{Å}}$ にもなる．実際，8.8節で行う構造緩和(解析的な微分を使う)では $R_{\mathrm{e}}=1.1948\ \text{Å}$ が得られ，放物線内挿の $1.1978\ \text{Å}$ とは $0.003\ \text{Å}$ ずれる．曲線を細かく刻めばこの差は縮まるが，それより解析的な力を使うほうが賢い．それが8.7節の主題である．
8.6 三重項酸素 — 高スピンと低スピン
8.5節でさりげなく --spin 2 と指定したが，これは何を意味しているのか．実は酸素分子は，化学を学ぶうえで極めて重要な例外的性質をもっている．
8.6.1 酸素原子と酸素分子の電子配置
酸素の原子番号は $8$，電子配置は $1s^{2}2s^{2}2p^{4}$ である．$2p$ 軌道は3本($2p_x,2p_y,2p_z$)あり，そこに4個の電子を入れる．Hundの規則(第7章)により，1本には対をつくって2個，残り2本には平行スピンで1個ずつ入る．したがって不対電子が2個あり，スピン量子数は $S=\tfrac12+\tfrac12=1$，多重度は $2S+1=3$ で三重項である(項記号 $^{3}\mathrm{P}$)．
では $\mathrm{O}_2$ 分子ではどうか．2個の酸素原子から16個の電子が集まる．原子軌道が重なって分子軌道(molecular orbital, MO)をつくる様子を表8.2にまとめた．詳しい対称性の議論は第11章・第12章で行うので，ここでは結果だけを使う．
| 分子軌道 | 性格 | 電子数 | 備考 |
|---|---|---|---|
| $\sigma_{1s}$ | 結合性 | 2 | 内殻．結合性と反結合性が打ち消し合う |
| $\sigma^{*}_{1s}$ | 反結合性 | 2 | |
| $\sigma_{2s}$ | 結合性 | 2 | |
| $\sigma^{*}_{2s}$ | 反結合性 | 2 | |
| $\sigma_{2p}$ | 結合性 | 2 | |
| $\pi_{2p}$(2重縮退) | 結合性 | 4 | |
| $\pi^{*}_{2p}$(2重縮退) | 反結合性 | 2 | 2本に1個ずつ，平行スピン |
| $\sigma^{*}_{2p}$ | 反結合性 | 0 | 空 |
ここが決定的である．最高被占軌道である $\pi^{*}_{2p}$ は2本が縮退しており，そこに入る電子は2個しかない．Hundの規則に従えば，2個の電子は別々の軌道に平行スピンで入る．したがって $\mathrm{O}_2$ の基底状態も
$$ \begin{equation} S=\tfrac{1}{2}+\tfrac{1}{2}=1,\qquad \text{多重度}=2S+1=3 \label{eq:8-multiplicity} \end{equation} $$の三重項である．式 \eqref{eq:8-multiplicity} の $2S$ を PySCF は指定させるので，--spin 2 となる．これが8.5節の指定の正体である．項記号は $^{3}\Sigma_{g}^{-}$ と書くが，この記号の読み方(なぜ $\Sigma$ で，なぜ $g$ で，なぜ $-$ なのか)は第12章で対称性と群論を学んでから説明する．
ついでに結合次数も数えておこう．価電子だけで数えると，結合性軌道の電子は $\sigma_{2s}(2)+\sigma_{2p}(2)+\pi_{2p}(4)=8$ 個，反結合性軌道の電子は $\sigma^{*}_{2s}(2)+\pi^{*}_{2p}(2)=4$ 個だから
$$ \text{結合次数}=\frac{8-4}{2}=2 . $$$\mathrm{O}_2$ は二重結合である．$1s$ 由来の $\sigma_{1s}$ と $\sigma^{*}_{1s}$ は互いに打ち消し合うので数に入れても答えは変わらない．
物理的意味:なぜ液体酸素は磁石に引かれるのか
高校化学で習うLewis構造式では，$\mathrm{O}_2$ は $\mathrm{O}\!=\!\mathrm{O}$ と書かれ，すべての電子が対をつくっている．この描像が正しければ酸素分子は反磁性で，磁石にはまったく反応しないはずである．ところが実際には，液体酸素を磁石の極の間に注ぐと，青い液体が磁極に吸い寄せられて留まる．これは酸素が常磁性(paramagnetic)であることの，目で見える証拠である．
その理由が表8.2の $\pi^{*}_{2p}$ にある．ここに2個の不対電子が平行スピンで入っているため，分子全体が正味のスピン磁気モーメント($S=1$)をもつ．Lewis構造式では説明できず，分子軌道法で初めて説明できた——この一点だけでも，分子軌道法が化学結合論に革命をもたらしたことがわかる．計算の --spin 2 という一行は，この物理を反映しているのである．
8.6.2 高スピンと低スピンを比べる
では，無理やり2個の電子を同じ $\pi^{*}_{2p}$ 軌道に押し込んで対にしたら(一重項，$S=0$，低スピン状態)，エネルギーはどうなるだろうか．8.5節とまったく同じ手順を --spin 0 で繰り返せばよい．結果を表8.1の第3列と図8.4に示す．
図8.4から読み取れることは三つある．
- 全体的にエネルギーが上昇した．平衡点での差は $-4067.059-(-4069.348)=2.290\ \mathrm{eV}$ である．三重項のほうが安定であり，$\mathrm{O}_2$ の基底状態は確かに三重項である．
- 曲線の形はほとんど変わらない．放物線内挿すると，低スピンの平衡結合長は $R_{\mathrm{e}}=1.196\ \text{Å}$ で，高スピンの $1.198\ \text{Å}$ とほぼ同じである．これは当然で，スピンの向きを変えても $\pi^{*}_{2p}$ に入る電子の総数は2個のままだから，結合次数は $2$ で変わらないのである．
- したがって，結合長の違いは目で見てもわからない．「結合長はどうなっているのだろう」という問いに答えるには，曲線を26点も計算するのではなく，力を使って一気に最小点に落とす方法がほしい．それが次節以降の主題である．
なぜ？：平行スピンのほうが安定になる理由(第9章の予告)
2個の電子を「別々の軌道に平行スピンで置く」ほうが「同じ軌道に反平行スピンで置く」より安定になる理由は，二つある．第一に，同じ軌道に押し込むと2個の電子が空間的に重なり，Coulomb反発 $V_{\mathrm{el\text{-}el}}$ が増える．第二に，そしてこちらがより本質的だが，平行スピンの電子は Pauli の排他原理によって互いに近づけず，そのぶんCoulomb反発が下がる．この後者の効果を交換エネルギー(exchange energy)という．Slater行列式(第7章)から交換エネルギーが自動的に出てくる仕組みは第9章で詳しく扱う．図8.4の $2.29\ \mathrm{eV}$ という差は，まさに交換エネルギーの大きさを目で見ているのである．
注意:「低スピン計算」は真の一重項状態ではない
ここで --spin 0 として得た解は，閉殻の単一Slater行列式(RHF解)であって，実験でいう $\mathrm{O}_2$ の第一励起状態 $a^{1}\Delta_{g}$ とは別物である．実験値では $a^{1}\Delta_{g}$ は基底状態より $0.98\ \mathrm{eV}$ 高いだけで，計算の $2.29\ \mathrm{eV}$ とはかなり違う．$^{1}\Delta_{g}$ を正しく記述するには二つのSlater行列式の重ね合わせが必要で，単一行列式では原理的に表せないからである．
また，UHF で --spin 0 を指定すると，初期値によっては $\braket{\hat{S}^{2}}\simeq1$ という「スピン対称性の破れた解」に落ちることがあり，これも純粋な一重項ではない．単一行列式近似の限界を意識せずに数値を鵜呑みにしてはいけない．
注意:Hartree–Fockは結合の解離をまったく再現できない
図8.3の曲線は右端($2.15\ \text{Å}$)でもまだ急な傾きをもっており，解離極限に達していない．真の解離極限は，酸素原子を1個だけ計算して2倍すれば得られる．UHF/6-31G で計算すると $E(\mathrm{O})=-2034.876\ \mathrm{eV}$，すなわち $2E(\mathrm{O})=-4069.752\ \mathrm{eV}$ となる．ところが分子の最安定値は $-4069.350\ \mathrm{eV}$ であるから，
$$D_{\mathrm{e}}=2E(\mathrm{O})-E(\mathrm{O_2})=-0.40\ \mathrm{eV}\ \lt\ 0 .$$なんと，この計算は$\mathrm{O}_2$ 分子は2個の酸素原子より不安定であると予言している．分極関数を加えた基底(6-31G* や cc-pVTZ)にすると $D_{\mathrm{e}}=+1.3\sim1.4\ \mathrm{eV}$ と正になるが，実験値 $5.21\ \mathrm{eV}$ には遠く及ばない．
これは Hartree–Fock 法の有名な弱点で，原因は電子相関の欠如にある．平衡付近の曲線の形(結合長・曲率)はそこそこ再現できても，深さ(結合エネルギー)はまったく駄目なのである．この失敗をどう克服するかが，第9章以降で交換相関エネルギーを論じる動機になる．「よく再現できた」と喜ぶときは，何が再現できて何ができていないかを必ず点検してほしい．
8.7 Hellmann–Feynmanの定理と原子間力
26通りの構造を計算して曲線を描くのは，二原子分子だからできたことである．$N$ 原子の分子では $3N-6$ 次元の曲面を格子状に走査しなければならず，$N=10$ でも24次元，各方向10点でも $10^{24}$ 回の計算になって完全に不可能である．必要なのは，1回の電子状態計算から「どちらに動けばエネルギーが下がるか」を直接知る方法である．それが原子間力の計算であり，理論的な支柱が Hellmann–Feynman の定理である．
8.7.1 電子系のエネルギーを核座標の関数とみる
変分原理に基づき，核の位置を $\bm{R}$ に固定したうえで電子系のエネルギーを最小化して，電子の波動関数 $\Psi_{\mathrm{el}}(\bm{r},\bm{R})$ が決まったとする．そのとき得られるエネルギーは
$$ \begin{equation} \epsilon_{\mathrm{el}}(\bm{R})=\bra{\Psi_{\mathrm{el}}}\Ham_{\mathrm{el}}+V_{\mathrm{Nu\text{-}Nu}}\ket{\Psi_{\mathrm{el}}}, \label{eq:8-eps-el} \end{equation} $$ここで電子ハミルトニアンと核間反発は
$$ \Ham_{\mathrm{el}}=\sum_{i=1}^{n}\frac{\hat{\bm{p}}_i^{2}}{2m} -\sum_{I}\sum_{i}\frac{Z_I\ee^{2}}{4\pi\eps\abs{\bm{R}_I-\bm{r}_i}} +\frac{1}{2}\sum_{i}\sum_{j\neq i}\frac{\ee^{2}}{4\pi\eps\abs{\bm{r}_i-\bm{r}_j}}, \qquad V_{\mathrm{Nu\text{-}Nu}}=\frac{1}{2}\sum_{I}\sum_{J\neq I}\frac{Z_IZ_J\ee^{2}}{4\pi\eps\abs{\bm{R}_I-\bm{R}_J}} $$である．$\Psi_{\mathrm{el}}$ は規格化されている($\braket{\Psi_{\mathrm{el}}|\Psi_{\mathrm{el}}}=1$)としてよい．$\epsilon_{\mathrm{el}}(\bm{R})$ は，8.5節で $\mathrm{O}_2$ について実際に計算した「全エネルギー」そのものであり，核座標 $\bm{R}$ の関数である．
式 \eqref{eq:8-eps-el} の $\epsilon_{\mathrm{el}}(\bm{R})$ こそが，核にとってのポテンシャルである．力学では，力はポテンシャルの勾配の符号を変えたもの，$\bm{F}=-\nabla U$ であった．核にとってのポテンシャルが $\epsilon_{\mathrm{el}}(\bm{R})$ なのだから，これを核の位置で微分すれば原子間力が得られそうである．
なぜ？：$\bm{F}=-\nabla U$ がしっくり来ない人へ
ばねを思い出そう．自然長からの伸びを $x$ とすると，ばねの弾性エネルギーは $U=\tfrac12kx^{2}$ である．このとき
$$F=-\frac{\dd U}{\dd x}=-\frac{\dd}{\dd x}\left(\frac{1}{2}kx^{2}\right)=-kx$$となり，確かにHookeの法則が出てくる．マイナス符号は「力はエネルギーが減る方向を向く」ことを表している．$x\gt0$(伸びている)なら力は縮む向き($F\lt0$)，というわけである．図8.3の曲線を坂道だと思えば，ボールは坂を下る向きに力を受ける——それだけのことである．
8.7.2 定理とその証明
問題は，$\epsilon_{\mathrm{el}}(\bm{R})$ の中で $\bm{R}$ に依存しているのがハミルトニアンだけでなく波動関数 $\Psi_{\mathrm{el}}$ もだ，という点である．$\bm{R}$ を動かせば電子分布も変わるのだから，素朴に微分すれば $\Psi_{\mathrm{el}}$ の微分の項が出てきてしまい，扱いが厄介になる．ところが，それらの項はきれいに消えてくれる．それが次の定理である．
定理:Hellmann–Feynmanの定理
ハミルトニアン $\Ham(\lambda)$ がパラメータ $\lambda$ に依存し，その規格化された固有状態を $\ket{\psi(\lambda)}$，固有値を $E(\lambda)$ とする:
$$\Ham(\lambda)\ket{\psi(\lambda)}=E(\lambda)\ket{\psi(\lambda)},\qquad \braket{\psi(\lambda)|\psi(\lambda)}=1 .$$このとき，固有値の $\lambda$ 微分はハミルトニアンの微分の期待値だけで与えられる:
$$ \begin{equation} \frac{\dd E(\lambda)}{\dd\lambda}=\bra{\psi(\lambda)}\frac{\partial\Ham(\lambda)}{\partial\lambda}\ket{\psi(\lambda)} . \label{eq:8-hf-theorem} \end{equation} $$波動関数の微分 $\partial\ket{\psi}/\partial\lambda$ は現れない．H. Hellmann(1937，文献[5])と R. P. Feynman(1939，文献[6])にちなむ．
導出:Hellmann–Feynmanの定理の証明
固有値は期待値として $E(\lambda)=\bra{\psi}\Ham\ket{\psi}$ と書ける($\braket{\psi|\psi}=1$ を使った)．これを $\lambda$ で微分する．積が三つ($\bra{\psi}$，$\Ham$，$\ket{\psi}$)あるので，積の微分法則から3項が出る:
$$ \begin{align} \frac{\dd E}{\dd\lambda} &=\left(\frac{\partial}{\partial\lambda}\bra{\psi}\right)\Ham\ket{\psi} +\bra{\psi}\frac{\partial\Ham}{\partial\lambda}\ket{\psi} +\bra{\psi}\Ham\left(\frac{\partial}{\partial\lambda}\ket{\psi}\right). \label{eq:8-hf-proof-1} \end{align} $$第1項と第3項を処理する．$\Ham$ は Hermite 演算子なので，右に作用させても左に作用させても固有値 $E$ を出す:
$$\Ham\ket{\psi}=E\ket{\psi},\qquad \bra{\psi}\Ham=E\bra{\psi}.$$これを使うと，第1項と第3項はどちらも $E$ を括り出せて
$$ \begin{align} \left(\frac{\partial}{\partial\lambda}\bra{\psi}\right)\Ham\ket{\psi} +\bra{\psi}\Ham\left(\frac{\partial}{\partial\lambda}\ket{\psi}\right) &=E\left[\left(\frac{\partial}{\partial\lambda}\bra{\psi}\right)\ket{\psi} +\bra{\psi}\left(\frac{\partial}{\partial\lambda}\ket{\psi}\right)\right]\notag\\ &=E\,\frac{\partial}{\partial\lambda}\braket{\psi|\psi} \label{eq:8-hf-proof-2} \end{align} $$となる．ここで再び積の微分法則を(今度は逆向きに)使ったことに注意してほしい．ところが規格化条件から $\braket{\psi|\psi}=1$，すなわち $\lambda$ によらない定数である．定数を微分すればゼロだから
$$ \begin{equation} \frac{\partial}{\partial\lambda}\braket{\psi|\psi}=\frac{\partial}{\partial\lambda}(1)=0 . \label{eq:8-hf-proof-3} \end{equation} $$よって式 \eqref{eq:8-hf-proof-2} は式 \eqref{eq:8-hf-proof-3} によりゼロとなる．すなわち式 \eqref{eq:8-hf-proof-1} の第1項と第3項は合わせてゼロになり，残るのは第2項だけである:
$$\frac{\dd E}{\dd\lambda}=\bra{\psi}\frac{\partial\Ham}{\partial\lambda}\ket{\psi}. \qquad\blacksquare$$証明の急所は，波動関数の微分が個別にはゼロでないのに，規格化条件のおかげで足すとゼロになることである．$\psi$ が実関数なら $\partial\braket{\psi|\psi}/\partial\lambda=2\braket{\psi|\partial_\lambda\psi}$ なので，$\braket{\psi|\partial_\lambda\psi}=0$，つまり波動関数の変化はつねに自分自身と直交する方向に起きる，という幾何学的な意味づけもできる．
導出:変分的に最適化した近似波動関数でも成り立つ
実際の計算では，$\Psi_{\mathrm{el}}$ は厳密な固有関数ではなく，有限個のパラメータ $\bm{c}=(c_1,c_2,\dots)$(たとえばLCAO係数)をもつ試行関数を変分的に最適化したものである．この場合でも定理は成り立つ．エネルギーを $E(\bm{c},\lambda)$ と書き，最適解を $\bm{c}^{*}(\lambda)$ とすると，全微分は連鎖律により
$$ \frac{\dd E}{\dd\lambda} =\frac{\partial E}{\partial\lambda}\bigg|_{\bm{c}} +\sum_{k}\underbrace{\frac{\partial E}{\partial c_k}}_{=\,0}\frac{\dd c_k}{\dd\lambda} =\frac{\partial E}{\partial\lambda}\bigg|_{\bm{c}} $$となる．変分の停留条件 $\partial E/\partial c_k=0$ がすべての $k$ について成り立っているからである．したがって，変分的に完全に最適化された波動関数であれば，係数の変化を追わなくてよい．Hartree–Fock法やKohn–Sham法のSCF計算が収束していれば，この条件は満たされている．
8.7.3 Hellmann–Feynman力
定理をわれわれの問題に適用しよう．パラメータ $\lambda$ を核 $I$ の座標 $\bm{R}_I$ の各成分とすればよい．力は $\bm{F}_I=-\nabla_I\epsilon_{\mathrm{el}}(\bm{R})$ だから
これをHellmann–Feynman力(Hellmann–Feynman force)という．まるで電子系のエネルギーをポテンシャルエネルギーとして扱っているように見えるが，それがまさにBO近似(断熱近似)の意味するところなのである．
式 \eqref{eq:8-hf-force} の中身を具体的に書き下してみよう．$\Ham_{\mathrm{el}}$ のうち $\bm{R}_I$ を含むのは核–電子引力の項だけである．$\nabla_I\dfrac{1}{\abs{\bm{R}_I-\bm{r}_i}}=-\dfrac{\bm{R}_I-\bm{r}_i}{\abs{\bm{R}_I-\bm{r}_i}^{3}}$ を使うと
$$ \nabla_I\Ham_{\mathrm{el}} =-\sum_{i}\frac{Z_I\ee^{2}}{4\pi\eps}\left(-\frac{\bm{R}_I-\bm{r}_i}{\abs{\bm{R}_I-\bm{r}_i}^{3}}\right) =+\frac{Z_I\ee^{2}}{4\pi\eps}\sum_{i}\frac{\bm{R}_I-\bm{r}_i}{\abs{\bm{R}_I-\bm{r}_i}^{3}} . $$この期待値をとるとき，$n$ 個の電子についての和 $\sum_i$ は，電子密度(electron density)
$$ \rho(\bm{r})=n\int\abs{\Psi_{\mathrm{el}}(\bm{r},\bm{r}_2,\dots,\bm{r}_n)}^{2}\,\dd^{3}r_2\cdots\dd^{3}r_n $$による積分に置き換わる．$\rho(\bm{r})\dd^{3}r$ は「点 $\bm{r}$ のまわりの微小体積に電子が見つかる個数の期待値」であり，$\int\rho\,\dd^{3}r=n$ である．また核間反発の微分は
$$ \nabla_I V_{\mathrm{Nu\text{-}Nu}} =-\sum_{J\neq I}\frac{Z_IZ_J\ee^{2}}{4\pi\eps}\frac{\bm{R}_I-\bm{R}_J}{\abs{\bm{R}_I-\bm{R}_J}^{3}} $$である．以上をまとめると
$$ \begin{equation} \bm{F}_I= \underbrace{\frac{Z_I\ee^{2}}{4\pi\eps}\int \rho(\bm{r})\,\frac{\bm{r}-\bm{R}_I}{\abs{\bm{r}-\bm{R}_I}^{3}}\,\dd^{3}r}_{\text{電子雲が核 }I\text{ を引く力}} \;+\; \underbrace{\frac{Z_I\ee^{2}}{4\pi\eps}\sum_{J\neq I}Z_J\,\frac{\bm{R}_I-\bm{R}_J}{\abs{\bm{R}_I-\bm{R}_J}^{3}}}_{\text{他の核が核 }I\text{ を押す力}} \label{eq:8-hf-electrostatic} \end{equation} $$を得る．
物理的意味:Feynmanの静電定理
式 \eqref{eq:8-hf-electrostatic} をよく見てほしい．右辺には $\hbar$ も $m$ も現れず，電子の運動エネルギーも電子間反発も出てこない．あるのは電荷分布がつくる古典的な静電気力だけである．つまり，
核に働く力は，電子密度 $\rho(\bm{r})$ と他の核とがつくる静電気力にほかならない．
量子力学が必要なのは，正しい $\rho(\bm{r})$ を求めるためだけであって，いったん $\rho$ が求まれば，あとは高校物理のCoulombの法則で力が計算できる．Feynman は学部生のときにこれを見出し，卒業論文にした(文献[6])．「化学結合とは，核と核の間に電子密度がたまり，その電子雲が両側の核を内向きに引っ張っている状態である」という直感的な結合の描像は，この定理から自然に出てくる．
例題8.4 水素原子の核に働く力はゼロである
水素原子の $1s$ 電子は球対称な密度分布 $\rho(\bm{r})=\abs{\phi_{1s}(r)}^{2}$ をもつ．式 \eqref{eq:8-hf-electrostatic} を使って，原点にある核に働く力がゼロであることを示せ．
解答:核は1個なので第2項(核間反発)はない．第1項は $\bm{R}_I=\bm{0}$ として
$$\bm{F}=\frac{\ee^{2}}{4\pi\eps}\int\rho(r)\,\frac{\bm{r}}{r^{3}}\,\dd^{3}r .$$被積分関数のうち $\rho(r)/r^{3}$ は $r=\abs{\bm{r}}$ だけの関数(スカラー)であり，$\bm{r}$ はベクトルである．$\bm{r}\to-\bm{r}$ という変換を考えると，$r$ は変わらないのでスカラー部分は不変，一方ベクトル $\bm{r}$ は符号を変える．積分領域(全空間)もこの変換で不変だから，積分の値は
$$\bm{F}=-\bm{F}\quad\Longrightarrow\quad \bm{F}=\bm{0}.$$球対称な電荷分布はその中心にある電荷に力を及ぼさない，というGaussの法則でおなじみの事実に対応する．孤立原子が自発的に動き出さないという当たり前の結論が，きちんと出てくる．
なお，この原子を別の原子に近づけると密度が球対称でなくなり(分極する)，核は結合の相手側へ引かれるようになる．これが結合が生じるときの力の起源である．
例題8.5 $\mathrm{O}_2$ に働く力の大きさを数値で実感する
表8.1の高スピンのデータから，$R=0.95\ \text{Å}$ における原子に働く力を中心差分で見積もり，SI単位に直せ．また $R=1.20\ \text{Å}$ ではどうなるか．
解答:中心差分は $\dfrac{\dd E}{\dd R}\simeq\dfrac{E(R+h)-E(R-h)}{2h}$，$h=0.05\ \text{Å}$．$R=0.95\ \text{Å}$ では
$$\frac{\dd E}{\dd R}\simeq\frac{-4066.428637-(-4060.401153)}{0.10}=\frac{-6.027484}{0.10}=-60.3\ \mathrm{eV/\text{Å}} .$$力は $F=-\dd E/\dd R=+60.3\ \mathrm{eV/\text{Å}}$，すなわち核間距離を伸ばす向きである．SI単位では
$$F=60.3\times\frac{1.6022\times10^{-19}\ \mathrm{J}}{1.0\times10^{-10}\ \mathrm{m}}=9.7\times10^{-8}\ \mathrm{N}\simeq0.10\ \mathrm{\mu N}.$$原子1個にかかる力としては桁違いに大きい．原子スケールでは $\mathrm{nN}$ ($10^{-9}\ \mathrm{N}$)が力の目安であり，$0.1\ \mathrm{\mu N}$ はその100倍である．それだけ $0.95\ \text{Å}$ という結合長が「押し込まれすぎ」ということである．
一方 $R=1.20\ \text{Å}$ では
で，$F=-0.196\ \mathrm{eV/\text{Å}}$ とほぼゼロ(わずかに縮む向き)．平衡点のすぐ近くにいることが力からも確認できる．実際に解析的な微分で計算すると $R=1.20\ \text{Å}$ での力は $-0.47\ \mathrm{eV/\text{Å}}$ であり，平衡は $R_{\mathrm{e}}=1.1948\ \text{Å}$ にある．
注意:Pulay力 — 局在基底では補正が要る
8.7.2節の変分版の証明には，隠れた前提がある．「変分パラメータ $\bm{c}$ について最適化してあれば $\partial E/\partial c_k=0$」というのは正しいが，変分空間そのものが $\bm{R}$ に依存していないことが必要である．ところが，量子化学計算で使うGauss型基底関数(第5章・第6章)は原子核に貼りついて一緒に動く．核を動かすと基底関数の中心も動き，変分空間が変形してしまう．
基底関数を $\chi_\mu(\bm{r}-\bm{R}_{I(\mu)})$ と書けば，エネルギーの全微分には
$$\frac{\dd E}{\dd\bm{R}_I}=\underbrace{\bra{\Psi}\nabla_I\Ham\ket{\Psi}}_{\text{Hellmann–Feynman項}}+\underbrace{2\,\mathrm{Re}\bra{\nabla_I\Psi}\Ham-E\ket{\Psi}}_{\text{Pulay項}}$$という余分な寄与が現れる．これをPulay力(Pulay force，基底関数の不完全性による力)といい，P. Pulay が1969年に指摘した(文献[7])．無視すると力が数十%も狂うことがあるので，実際のプログラムは必ずこの補正を含めて微分を評価している．
一方，固体計算でよく使う平面波基底は原子位置に依存しない(空間全体に一様に張られている)ので，原理的にPulay力は現れない．ただし，平面波の数を打ち切るカットオフエネルギーを固定したまま格子定数(セル体積)を変えると，基底の数が実効的に変わるため類似の補正(Pulay応力)が必要になる．基底の完全性が上がればどちらの補正も小さくなる，という点は共通である．
8.7.4 Hessian行列 — 力の次の階
力は1階微分だった．2階微分を集めた行列をHessian行列(Hessian matrix)といい，成分は
$$ \begin{equation} \Phi_{I\alpha,J\beta}=\frac{\partial^{2}\epsilon_{\mathrm{el}}}{\partial R_{I\alpha}\,\partial R_{J\beta}} =-\frac{\partial F_{I\alpha}}{\partial R_{J\beta}} \qquad(\alpha,\beta=x,y,z) \label{eq:8-hessian} \end{equation} $$である($3N\times3N$ の対称行列)．物理的には力の定数(force constant)であり，「核 $J$ を $\beta$ 方向にずらしたとき，核 $I$ が $\alpha$ 方向にどれだけ力を受けるか」を表す．平衡点のまわりで $\epsilon_{\mathrm{el}}$ をTaylor展開すると，1階の項は力がゼロなので消え，
$$ \epsilon_{\mathrm{el}}(\bm{R}_{\mathrm{e}}+\bm{u})=\epsilon_{\mathrm{el}}(\bm{R}_{\mathrm{e}}) +\frac{1}{2}\sum_{I\alpha}\sum_{J\beta}\Phi_{I\alpha,J\beta}\,u_{I\alpha}u_{J\beta}+O(u^{3}) $$となる．これは多次元のばねにほかならない．質量で重みづけした行列 $D_{I\alpha,J\beta}=\Phi_{I\alpha,J\beta}/\sqrt{M_IM_J}$ を対角化すれば，固有値が $\omega^{2}$，固有ベクトルが基準振動モードを与える．二原子分子では自由度が1本だけ残り，$\mu=M_AM_B/(M_A+M_B)$ を換算質量として
$$ \begin{equation} \omega=\sqrt{\frac{k}{\mu}},\qquad k=\frac{\dd^{2}\epsilon_{\mathrm{el}}}{\dd R^{2}}\bigg|_{R=R_{\mathrm{e}}} \label{eq:8-omega} \end{equation} $$となる．例題8.3で計算したのは，まさに式 \eqref{eq:8-omega} の $k$ と $\omega$ である．
式 \eqref{eq:8-hessian} で定義したHessianの用途は，振動スペクトルの予測だけではない．固有値の符号を見れば，その構造が本当に安定かどうかを判定できる．すべて正なら極小(安定構造)，1つだけ負なら鞍点(遷移状態)，負が複数あればさらに不安定な停留点である．負の固有値をもつモード(虚振動)は，構造がどちらに壊れたがっているかを教えてくれる．担当教員の研究でも，この虚振動を手がかりに新しい強誘電体を探索している(文献[13])．
8.8 構造緩和のアルゴリズム
原子間力が求まれば，そのベクトルに従って原子核を動かせばよい．そうすると全エネルギーがより減少する方向に構造が変化する．これが構造緩和(structural relaxation)あるいは構造最適化(geometry optimization)である．手順を図8.7のフローチャートにまとめた．
8.8.1 どちらへ，どれだけ動かすか
フローチャートの⑥「原子を少し動かす」をどう実装するかが，最適化アルゴリズムの中身である．もっとも素朴なのは最急降下法(steepest descent)で，力の向きにそのまま動かす:
$$ \begin{equation} \bm{R}^{(n+1)}=\bm{R}^{(n)}+\alpha\,\bm{F}^{(n)},\qquad \alpha\gt0 . \label{eq:8-steepest} \end{equation} $$式 \eqref{eq:8-steepest} は坂道をボールが転がり落ちるのと同じで，必ずエネルギーは下がる($\alpha$ が十分小さければ)．しかし収束は遅い．細長い谷に入ると，谷底に沿って進まず左右の斜面をジグザグに往復してしまうからである．
もっと賢いやり方は，2階微分の情報を使うことである．エネルギーを現在位置のまわりで2次まで展開し
$$ \epsilon_{\mathrm{el}}(\bm{R}^{(n)}+\Delta\bm{R})\simeq \epsilon_{\mathrm{el}}(\bm{R}^{(n)})-\bm{F}^{(n)}\!\cdot\!\Delta\bm{R} +\frac{1}{2}\Delta\bm{R}^{\mathsf{T}}\Phi\,\Delta\bm{R}, $$これを $\Delta\bm{R}$ について最小化する．$\partial/\partial\Delta\bm{R}$ をとってゼロとおけば $-\bm{F}^{(n)}+\Phi\Delta\bm{R}=0$，すなわち
$$ \begin{equation} \Delta\bm{R}=\Phi^{-1}\bm{F}^{(n)} . \label{eq:8-newton} \end{equation} $$式 \eqref{eq:8-newton} がNewton法である．もしエネルギーがぴったり2次関数なら，この1歩で最小点に到達する．谷が細長くても，Hessian $\Phi$ がその細長さを知っているので，正しく谷底へ向かってくれる．
ただし $\Phi$ を毎回計算するのは高くつく($3N\times3N$ 個の2階微分が要る)．そこで，これまでの歩みの履歴から $\Phi$ を推定しながら進む準Newton法(quasi-Newton method)が使われる．基本になるのは割線条件(secant condition)である．1次元の $f'(x)$ を思い出すと $f''\simeq\Delta f'/\Delta x$ だった．多次元でも同じで，2ステップ間の変位 $\bm{s}=\bm{R}^{(n+1)}-\bm{R}^{(n)}$ と勾配の変化 $\bm{y}=\nabla\epsilon^{(n+1)}-\nabla\epsilon^{(n)}$ に対して，近似Hessian $B$ は
$$ \begin{equation} B_{n+1}\,\bm{s}=\bm{y} \label{eq:8-secant} \end{equation} $$を満たすべきである．この割線条件 \eqref{eq:8-secant} を満たす更新式のうち，もっとも広く使われるのがBFGS法(Broyden–Fletcher–Goldfarb–Shanno法，文献[10])で，
$$ \begin{equation} B_{n+1}=B_{n}+\frac{\bm{y}\bm{y}^{\mathsf{T}}}{\bm{y}^{\mathsf{T}}\bm{s}}-\frac{B_{n}\bm{s}\,\bm{s}^{\mathsf{T}}B_{n}}{\bm{s}^{\mathsf{T}}B_{n}\bm{s}} \label{eq:8-bfgs} \end{equation} $$である．式 \eqref{eq:8-bfgs} が割線条件を満たすことは，右から $\bm{s}$ を掛けるだけで確かめられる:
$$ B_{n+1}\bm{s}=B_{n}\bm{s}+\bm{y}\,\frac{\bm{y}^{\mathsf{T}}\bm{s}}{\bm{y}^{\mathsf{T}}\bm{s}}-B_{n}\bm{s}\,\frac{\bm{s}^{\mathsf{T}}B_{n}\bm{s}}{\bm{s}^{\mathsf{T}}B_{n}\bm{s}} =B_{n}\bm{s}+\bm{y}-B_{n}\bm{s}=\bm{y}. \qquad\checkmark $$$\bm{y}\bm{y}^{\mathsf{T}}$ は列ベクトルと行ベクトルの積で $3N\times3N$ 行列になること，$\bm{y}^{\mathsf{T}}\bm{s}$ はスカラー(内積)になることに注意すれば，行列の積の計算だけである．
さらに実用上は，2次近似が信用できる範囲を制限する信頼半径(trust radius)が併用される．1歩の変位がこの半径を超えないように切り詰め，実際のエネルギー低下が予測どおりなら半径を広げ，外れたら狭める，という制御を行う．
8.8.2 O₂の構造緩和を実行する
それでは，わざと $0.9\ \text{Å}$ という短すぎる結合長から出発して，構造緩和が本当に平衡点を見つけてくれるか試してみよう．付録Bの calc_pyscf_relax.py は，PySCF と geomeTRIC(文献[9])を組み合わせて図8.7のループを回すスクリプトである．
$ cd ~/Downloads/pyscf_tutorial3/O2
$ python3 calc_pyscf_relax.py --xyz XYZ_O2_0p900000.xyz --spin 2 --basis 6-31gコード8.12 結合長 $0.9\ \text{Å}$ の構造から構造緩和を開始する．
実行すると，次のようなログが流れる．
Step 0 : Gradient = 1.701e+00/1.701e+00 (rms/max) Energy = -149.2169914759
Step 1 : Displace = 6.061e-02/6.061e-02 Trust = 1.000e-01 Grad = 6.080e-01 E = -149.4655345067 (-2.485e-01)
Step 2 : Displace = 3.370e-02/3.370e-02 Trust = 1.414e-01 Grad = 2.844e-01 E = -149.5208309326 (-5.530e-02)
Step 3 : Displace = 2.962e-02/2.962e-02 Trust = 2.000e-01 Grad = 1.001e-01 E = -149.5416474412 (-2.082e-02)
Step 4 : Displace = 1.608e-02/1.608e-02 Trust = 2.828e-01 Grad = 2.795e-02 E = -149.5454546188 (-3.807e-03)
Step 5 : Displace = 6.231e-03/6.231e-03 Trust = 3.000e-01 Grad = 4.194e-03 E = -149.5458288490 (-3.742e-04)
Step 6 : Displace = 1.100e-03/1.100e-03 Trust = 3.000e-01 Grad = 2.184e-04 E = -149.5458380000 (-9.151e-06)
Step 7 : Displace = 6.044e-05/6.044e-05 Trust = 3.000e-01 Grad = 8.388e-07 E = -149.5458380253 (-2.523e-08)
Converged!コード8.13 構造緩和のログ(要点のみ抜粋)．エネルギーはHartree単位，勾配はHartree/Bohr単位，変位はÅ単位である．
各ステップでの結合長も含めて表8.3に整理した．
| Step | 結合長 (Å) | 全エネルギー (Ha) | $\Delta E$ (Ha) | 勾配 (Ha/Bohr) |
|---|---|---|---|---|
| 0 | 0.900000 | −149.2169914759 | — | $1.701\times10^{0}$ |
| 1 | 1.021220 | −149.4655345067 | $-2.485\times10^{-1}$ | $6.080\times10^{-1}$ |
| 2 | 1.088629 | −149.5208309326 | $-5.530\times10^{-2}$ | $2.844\times10^{-1}$ |
| 3 | 1.147876 | −149.5416474412 | $-2.082\times10^{-2}$ | $1.001\times10^{-1}$ |
| 4 | 1.180044 | −149.5454546188 | $-3.807\times10^{-3}$ | $2.795\times10^{-2}$ |
| 5 | 1.192506 | −149.5458288490 | $-3.742\times10^{-4}$ | $4.194\times10^{-3}$ |
| 6 | 1.194706 | −149.5458380000 | $-9.151\times10^{-6}$ | $2.184\times10^{-4}$ |
| 7 | 1.194827 | −149.5458380253 | $-2.523\times10^{-8}$ | $8.388\times10^{-7}$ |
出力として，構造緩和された分子構造が XYZ_O2_0p900000-finish.xyz というファイルで得られる．
$ cat XYZ_O2_0p900000-finish.xyz
2
Optimized geometry, E = -149.5458380253 Ha
O -0.14741338 0.00000000 0.00000000
O 1.04741338 -0.00000000 -0.00000000コード8.14 緩和後の構造．結合長は $1.04741338-(-0.14741338)=1.19483\ \text{Å}$．座標の原点がずれているのは，最適化の途中で重心を固定しているためで，結合長そのものには影響しない．
結果をまとめよう．
- 結合長は $0.9\ \text{Å}$ から $1.195\ \text{Å}$ になった．図8.3の曲線の谷底にぴたりと落ち込んでいる．
- Step 0 から Step 7 まで，合計8回のSCF計算が実行された．原子を動かしたのは7回である．26通りの構造を総当たりした8.5節に比べ，計算回数は $1/3$ である．原子数が増えればこの差は天文学的に開く．
- 全エネルギーは $-149.5458380253\ \mathrm{Ha}$．$1\ \mathrm{Ha}=27.2114\ \mathrm{eV}$ を掛けて $-4069.3495\ \mathrm{eV}$ となり，格子上の最小値 $-4069.3483\ \mathrm{eV}$($R=1.20\ \text{Å}$)よりわずかに $1.3\ \mathrm{meV}$ だけ低い．谷底を格子で捉えきれていなかったぶんが，きちんと拾われている．
- 勾配は $1.7\to0.61\to0.28\to0.10\to0.028\to4.2\times10^{-3}\to2.2\times10^{-4}\to8.4\times10^{-7}$ と，最後は1ステップごとに桁が2つ以上落ちている．これはNewton法系のアルゴリズムに特徴的な超1次収束である．
例題8.6 最初の1歩はなぜ $0.12\ \text{Å}$ なのか
(1) 表8.3の Step 0 の勾配 $1.701\ \mathrm{Ha/Bohr}$ を $\mathrm{eV/\text{Å}}$ に換算せよ($1\ \mathrm{Ha}=27.2114\ \mathrm{eV}$，$1\ \mathrm{Bohr}=0.529177\ \text{Å}$)．
(2) 例題8.3で求めた力の定数 $k=89.6\ \mathrm{eV/\text{Å}^{2}}$ を使って Newton 法の1歩 $\Delta R=F/k$ を見積もり，実際の移動量 $1.021220-0.900000=0.1212\ \text{Å}$ と比べよ．
解答:(1) 単位換算係数は
$$1\ \mathrm{Ha/Bohr}=\frac{27.2114\ \mathrm{eV}}{0.529177\ \text{Å}}=51.42\ \mathrm{eV/\text{Å}}$$だから，$1.701\times51.42=87.5\ \mathrm{eV/\text{Å}}$．SI単位では $1.40\times10^{-7}\ \mathrm{N}$ である．
(2) Newton 法の1歩は
$$\Delta R=\frac{F}{k}=\frac{87.5\ \mathrm{eV/\text{Å}}}{89.6\ \mathrm{eV/\text{Å}^{2}}}=0.98\ \text{Å}.$$これは実際の移動量の8倍もある．$0.9\ \text{Å}$ から $0.98\ \text{Å}$ 動かせば $1.88\ \text{Å}$ になってしまい，最小点を大きく通り過ぎる．原因は，$k$ が平衡点での曲率であり，$0.9\ \text{Å}$ という壁のように急な領域では実際の曲率はずっと大きいからである．2次近似が通用しない遠方でNewton法を鵜呑みにすると暴走する．
そこで信頼半径が効く．ログの Step 1 は Trust = 1.000e-01，すなわち原子1個あたりの変位を $0.1\ \text{Å}$ に制限しており，実際の変位は $0.0606\ \text{Å}$ である．2個の原子が逆向きに動くので結合長は $2\times0.0606=0.1212\ \text{Å}$ 伸びた——これが観測された値と一致する．その後，予測が当たるにつれて Trust は $0.1\to0.141\to0.2\to0.283\to0.3$ と広げられ，収束が加速している．
注意:構造緩和は「いちばん近い窪み」に落ちるだけ
構造緩和が見つけるのは，初期構造から力に沿って下ったところにある局所最小点であって，系全体の最安定構造(大域最小点)である保証はまったくない．$\mathrm{O}_2$ のような1次元のPESなら問題は起きないが，原子数が増えると局所最小点の数は指数関数的に増える．実際の物質探索では，多数の初期構造から出発する，乱数で構造を生成する，進化的アルゴリズムを使う，といった構造探索の枠組みが別に必要になる(文献[12])．「構造最適化した」という言葉は，「大域最小を見つけた」という意味ではないことを，つねに意識してほしい．
8.9 まとめと演習
8.9.1 この章のまとめ
- 分子・固体の厳密なハミルトニアン \eqref{eq:8-full-hamiltonian} は五つの項からなる:核の運動エネルギー，電子の運動エネルギー，核–電子引力，核–核反発，電子–電子反発．
- 陽子は電子の約1836倍重いので，核は電子よりはるかにゆっくり動く．そこで波動関数を $\Psi\simeq\Psi_{\mathrm{el}}(\bm{r},\bm{R})\Psi_{\mathrm{Nu}}(\bm{R})$ と積に分けるのがBorn–Oppenheimer近似である．$\bm{R}$ は $\Psi_{\mathrm{el}}$ の変数ではなくパラメータとして入る．
- 積を代入して核のラプラシアンを展開すると非断熱項 $\hat\Lambda$ が現れるが，$1/M_I$ に比例して小さいので落とす．残った式は $\bm{R}$ を固定した電子の固有値問題 \eqref{eq:8-el-eigen} と，$E_{\mathrm{el}}(\bm{R})$ をポテンシャルとする核の方程式 \eqref{eq:8-nuc-eigen} に分かれる．
- 全エネルギーとは，核の運動エネルギーを $0$ とみなした $E\simeq\min_{\bm{R}}E_{\mathrm{el}}(\bm{R})$ のことで，$0\ \mathrm{K}$ の内部エネルギーである．有限温度では零点振動とエントロピーを別途考える必要がある．
- 断熱近似は「イオンは電子基底状態のポテンシャルエネルギー曲面(PES)上のみを動く」という言い換えであり，全エネルギー計算は電子→核の2段階最小化になる．
- $\mathrm{O}_2$ の結合長を $0.9$–$2.15\ \text{Å}$ で変えて計算すると，教科書どおりの原子間ポテンシャル曲線が得られる．最小は $R=1.20\ \text{Å}$(放物線内挿で $1.198\ \text{Å}$)で，実験値 $1.2075\ \text{Å}$ をよく再現する．ただし結合エネルギーはHartree–Fockでは全く再現できない．
- $\mathrm{O}_2$ の基底状態は $\pi^{*}_{2p}$ に不対電子2個をもつ三重項($S=1$，$^{3}\Sigma_{g}^{-}$，PySCFでは
spin = 2)である．一重項より $2.29\ \mathrm{eV}$ 安定で，この差の主因は交換エネルギー(第9章)である．液体酸素が磁石に引かれるのはこの不対電子のためである． - Hellmann–Feynmanの定理 \eqref{eq:8-hf-theorem} により，エネルギーのパラメータ微分はハミルトニアンの微分の期待値だけで書ける．規格化条件 $\braket{\psi|\psi}=1$ の微分がゼロになることが証明の急所である．核座標に適用すると原子間力 \eqref{eq:8-hf-force} が得られ，その中身は電子密度と他の核がつくる古典的な静電気力である．
- 原子中心のGauss基底では変分空間が $\bm{R}$ に依存するため，Pulay力の補正が必要である．平面波基底では原理的に不要．
- 2階微分の行列がHessian(力の定数)で，対角化すれば振動モードと振動数が得られ，固有値の符号で安定性が判定できる．
- 構造緩和は「SCF→力→収束判定→原子移動」のループである．最急降下法よりNewton法・準Newton法(BFGS)が速く，信頼半径で暴走を防ぐ．$\mathrm{O}_2$ では $0.9\ \text{Å}$ から8回のSCFで $1.195\ \text{Å}$ に収束した．
8.9.2 演習問題
演習8.1 質量比から時間スケールへ
(a) 陽子と電子の質量比を用いて，両者の運動量が同程度であるときの速度比を求めよ．(b) 運動エネルギーが同程度であるときの速度比を求め，(a) と比べよ．(c) 炭素核($A=12$)についても (b) を計算せよ．(d) これらの見積もりのうち，Born–Oppenheimer近似の正当化にとってより本質的なのはどちらか，理由とともに述べよ．
ヒント：(a) $v=P/M$ より $v_{\text{核}}/v_{\text{el}}=m/M$．(b) $\tfrac12Mv^{2}$ が等しいなら $v\propto1/\sqrt{M}$ なので比は $\sqrt{m/M}$．$\sqrt{1/1836}=1/42.8$．(c) $M=12\times1836\,m$．(d) 核と電子はCoulomb力で結ばれており，力積のやりとりを通じて運動量が同程度になる，という描像を考えよ．どちらの見積もりでも「核のほうが桁違いに遅い」という結論は変わらないことも確認せよ．
演習8.2 非断熱項を書き下す
(a) BO近似で捨てた非断熱項 $\hat\Lambda$ を，$\Psi_{\mathrm{el}}$，$\Psi_{\mathrm{Nu}}$，$\nabla_I$，$M_I$ を使って書き下せ．(b) $\hat\Lambda$ の各項が $1/M_I$ に比例することを確認し，それが小さい理由を説明せよ．(c) 二つの電子状態のエネルギーが接近する点(円錐交差)でBO近似が破綻する理由を，$\nabla_I\Psi_{\mathrm{el}}$ の振る舞いから説明せよ．
ヒント：(a) 式 \eqref{eq:8-laplacian-product} の第2項と第3項に $-\hbar^{2}/2M_I$ を掛けて和をとる．(b) 電子波動関数は核の位置に「なめらかに」依存する．(c) 摂動論(第4章)では，状態の変化量の分母にエネルギー差 $E_1-E_0$ が現れた．この分母がゼロに近づくとどうなるか．
演習8.3 粗い格子で放物線内挿するとどうなるか
表8.1の高スピンの値のうち $R=1.10,\ 1.20,\ 1.30\ \text{Å}$ の3点($h=0.10\ \text{Å}$)を使って，例題8.3と同じ手順で $R_{\mathrm{e}}$，$k$，$\tilde\nu$ を求めよ．刻み幅 $h=0.05\ \text{Å}$ の結果および構造緩和の結果 $1.1948\ \text{Å}$ と比較し，差の原因を論じよ．
ヒント：$E_+-2E_0+E_-=0.923875\ \mathrm{eV}$，$E_+-E_-=-0.126623\ \mathrm{eV}$ となるはずである．$k=(E_+-2E_0+E_-)/h^{2}$ で $h$ が2倍になると分母は4倍．答えは $R_{\mathrm{e}}=1.207\ \text{Å}$，$k=1.48\times10^{3}\ \mathrm{N/m}$，$\tilde\nu=1.77\times10^{3}\ \mathrm{cm^{-1}}$．刻みが粗いほど3次以上の非調和項を拾ってしまうことを，図8.6を見ながら確認せよ．
演習8.4 調和振動子でHellmann–Feynmanの定理を検算する
1次元調和振動子 $\Ham=\dfrac{\hat{p}^{2}}{2m}+\dfrac{1}{2}m\omega^{2}\hat{x}^{2}$ のエネルギー固有値は $E_n=\left(n+\tfrac12\right)\hbar\omega$ である．パラメータを $\lambda=\omega$ として定理 \eqref{eq:8-hf-theorem} を適用し，$\braket{n|\hat{x}^{2}|n}$ を求めよ．さらに，その結果からポテンシャルエネルギーの期待値が全エネルギーのちょうど半分になること(ビリアル定理)を確かめよ．
ヒント：$\partial\Ham/\partial\omega=m\omega\hat{x}^{2}$，$\dd E_n/\dd\omega=(n+\tfrac12)\hbar$．両者を等置して $\braket{\hat{x}^{2}}=\dfrac{(n+1/2)\hbar}{m\omega}$．ポテンシャルの期待値は $\tfrac12m\omega^{2}\braket{\hat{x}^{2}}$ である．波動関数の具体形(Hermite多項式)を一度も使わずに答えが出ることに注目せよ．
演習8.5 Pulay力はどんなときに現れるか
(a) 原子中心のGauss基底を用いた計算で，Hellmann–Feynman力だけでは正しい力が得られない理由を，8.7.2節の変分版の証明のどの前提が破れるかを指摘して説明せよ．(b) 平面波基底ではPulay力が現れない理由を述べよ．(c) 基底関数の数を無限に増やしていくと，Pulay力はどうなると期待されるか．
ヒント：(a) 「変分空間そのものが $\bm{R}$ に依存しない」という前提．基底関数 $\chi_\mu(\bm{r}-\bm{R}_{I})$ は核と一緒に動く．(b) 平面波 $e^{i\bm{G}\cdot\bm{r}}$ に核座標は入っていない．(c) 完全基底に近づくと，基底を動かしても張られる空間は変わらなくなる．ただしセル体積を変える場合はカットオフの扱いに注意．
演習8.6 酸素分子の磁性とスピン状態
(a) 表8.2をもとに $\mathrm{O}_2$ の $S$，多重度，結合次数を導け．(b) 液体酸素が磁石に引き寄せられる理由を，Lewis構造式では説明できず分子軌道法で説明できることを含めて述べよ．(c) 8.6節の低スピン計算($\braket{\hat{S}^{2}}=0$)が，実験でいう $a^{1}\Delta_{g}$ 状態(基底状態より $0.98\ \mathrm{eV}$ 高い)と一致しない理由を述べよ．(d) 図8.4で，高スピンと低スピンの平衡結合長がほとんど変わらないのはなぜか．
ヒント：(a) $\pi^{*}_{2p}$ は2重縮退で電子2個．Hundの規則．結合次数 $=$(結合性電子数 $-$ 反結合性電子数)$/2$．(c) $^{1}\Delta_{g}$ は2つのSlater行列式の重ね合わせでしか表せない(第7章・第9章)．(d) 結合次数を決めるのは軌道の占有数であって，スピンの向きではない．
8.9.3 参考文献
- 原田 義也『量子化学(上巻)』裳華房(2007).
- A. Szabo, N. S. Ostlund(大野公男・阪井健男・望月祐志 訳)『新しい量子化学 — 電子構造の理論入門(上)』東京大学出版会(1987)(原著:Modern Quantum Chemistry, Dover, 1996).
- P. Atkins, J. de Paula, Atkins' Physical Chemistry, 11th ed., Oxford University Press (2018).
- M. Born, J. R. Oppenheimer, "Zur Quantentheorie der Molekeln", Ann. Phys. 389, 457 (1927).
- H. Hellmann, Einführung in die Quantenchemie, Franz Deuticke, Leipzig und Wien (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. Theory", Mol. Phys. 17, 197 (1969).
- Q. Sun et al., "PySCF: the Python-based simulations of chemistry framework", WIREs Comput. Mol. Sci. 8, e1340 (2018).
- L.-P. Wang, C. Song, "Geometry optimization made simple with translation and rotation coordinates", J. Chem. Phys. 144, 214108 (2016)(geomeTRIC).
- J. Nocedal, S. J. Wright, Numerical Optimization, 2nd ed., Springer (2006).
- K. P. Huber, G. Herzberg, Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules, Van Nostrand Reinhold (1979)($\mathrm{O}_2$ の実験値:$R_{\mathrm{e}}=1.20752\ \text{Å}$，$\omega_{\mathrm{e}}=1580.19\ \mathrm{cm^{-1}}$，$D_0=5.115\ \mathrm{eV}$).
- Y. Mochizuki et al., "Theoretical Exploration of the Surface Structures of $\mathrm{SrTiO_3}$", Chem. Mater. 35, 2047 (2023).
- Y. Mochizuki et al., "$\mathrm{Li_2SrNb_2O_7}$: A New Ferroelectric", Chem. Mater. 33, 1257 (2021).
- C. Kittel(宇野良清 ほか 訳)『固体物理学入門(上)』第8版，丸善(2005).