第6章Gauss基底による近似解と厳密解の比較

第3章で変分原理を，第5章で基底関数（Slater型軌道・Gauss型軌道・STO-3G）を学んだ．道具はそろった．本章ではその道具で，水素原子の基底状態エネルギーを，厳密解を知らないふりをして，たった1個のGauss関数から手計算で求める．水素原子は厳密解が分かっている数少ない系であり，「近似がどれくらい当たるのか」を答え合わせできる，またとない練習台になる．結果は $E_{\min}=\dfrac{8}{3\pi}\epsilon_{1s}=-11.55\ \mathrm{eV}$，厳密値 $-13.606\ \mathrm{eV}$ より15%も高い．この15%はどこから来たのか——原点でのカスプと遠方の減衰という2つの欠陥を図で確かめる．そして「分かりきった面倒な単純作業は機械にやらせるべき」という方針でSTO-3GをPySCFに任せ，誤差が6.69%まで縮むことを見る．最後に，計算機が黙って選んでいる RHF・UHF・ROHF の違いを，Li原子とLi$^+$イオンの実計算で確かめる．

この章で学ぶこと

球対称な試行関数に対する動径ラプラシアン $\nabla^2=\dfrac{1}{r^2}\dfrac{\dd}{\dd r}\!\left(r^2\dfrac{\dd}{\dd r}\right)$ の使い方と，$\nabla^2\phi_T=(-6a+4a^2r^2)\phi_T$ の導出
Gauss積分公式 $\displaystyle\int_0^\infty x^{2n}e^{-cx^2}\dd x$ を微分でつくる方法と，運動エネルギー・Coulomb項・規格化積分の完全な手計算
$E(\alpha)=A\alpha-B\sqrt{\alpha}$ の最小化から $\alpha_{\mathrm{opt}}=\dfrac{8}{9\pi}=0.2829$，$E_{\min}=\dfrac{8}{3\pi}\epsilon_{1s}=-11.55\ \mathrm{eV}$ を得ること
核カスプ条件 $\left.\dfrac{1}{\psi}\dfrac{\dd\psi}{\dd r}\right|_{r\to0}=-\dfrac{Z}{a_0}$ の導出と，Gauss関数が原点で必ず平らになるという構造的欠陥
遠方の減衰 $e^{-\alpha r^2}$ と $e^{-r/a_0}$ の違い，変分最適な $\alpha$ と重なり最大の $\alpha$ がずれる理由
PySCFによるSTO-3G計算の実行と出力の読み方，STO-1G／STO-3G／6-31G／厳密解の定量比較
閉殻と開殻，$\bm{\alpha}$スピンと$\bm{\beta}$スピン，RHF・UHF・ROHFの違いと使い分け
Li原子・Li$^+$イオンの計算例，電荷を変えると準位が下がる理由，Koopmansの定理の芽

前提：第3章（変分原理と試行関数），第5章（Slater型軌道・Gauss型軌道・STO-3G・重なり積分），第2章（水素原子の厳密解と $\epsilon_{1s}$）．Gauss積分・極座標のラプラシアン・体積要素は付録A，PySCFの導入は付録Bを参照．

6.1 厳密解を知らないつもりで水素原子を解く

6.1.1 なぜ答えを知っている問題を解き直すのか

第2章で，水素原子のSchrödinger方程式は厳密に解けること，そのエネルギー準位が

$$ \begin{equation} E_n=\frac{\epsilon_{1s}}{n^2},\qquad \epsilon_{1s}=-\frac{m\ee^4}{8\eps^2h^2}=-13.606\ \mathrm{eV} \label{eq:6-exact-level} \end{equation} $$

となることを学んだ．基底状態（$n=1$）の波動関数も

$$ \begin{equation} \psi_{1s}(r)=\frac{1}{\sqrt{\pi a_0^3}}\,e^{-r/a_0},\qquad a_0=\frac{\eps h^2}{\pi m\ee^2}=0.529\times10^{-10}\ \mathrm{m} \label{eq:6-exact-1s} \end{equation} $$

という簡単な形をしている．答えが分かっている問題を，わざわざ近似で解き直すのはなぜか．

なぜ？：答え合わせのできる練習台が要る

実際に計算したい系——ヘリウム原子，メタン分子，SrTiO$_3$結晶——は，どれも厳密には解けない．解けないから近似する．ところが近似計算には，「その答えが本当はどれくらい正しいのか」を教えてくれる仕組みがない．変分原理（第3章）が保証してくれるのは「真の基底状態エネルギー以上である」という不等式だけで，どれだけ上にいるのかは分からない．

そこで，厳密解が分かっている水素原子を使って，同じ近似法を走らせてみる．誤差が何%出るのか，どんな形で誤るのかを一度身体で覚えておけば，答えの分からない系に対しても「この基底関数ならこの程度」という感覚が持てる．ベンチマーク（benchmark）と呼ばれるこの発想は，計算科学の基本作法である．

6.1.2 問題設定：ハミルトニアンと試行関数

水素原子の電子1個に対するハミルトニアンは，原子核（陽子）を原点に固定して

$$ \begin{equation} \Ham=-\frac{\hbar^2}{2m}\nabla^2-\frac{\ee^2}{4\pi\eps}\frac{1}{r} \label{eq:6-hamiltonian} \end{equation} $$

である．第1項が電子の運動エネルギー，第2項が原子核との Coulomb 引力エネルギーである．原子核を止めてよいという断りは Born–Oppenheimer 近似と呼ばれ，第8章できちんと正当化する．

試行関数として，第5章で導入したGauss型軌道（Gaussian-type orbital, GTO）を1個だけ使う．

定義：本章で使う試行関数（1個のGauss関数）

$$ \begin{equation} \phi_T(r)=N\exp\!\left[-\alpha\left(\frac{r}{a_0}\right)^{\!2}\right] =N e^{-ar^2}, \qquad a\equiv\frac{\alpha}{a_0^2} \label{eq:6-trial} \end{equation} $$

ここで $\alpha$ は無次元の変分パラメータ，$a=\alpha/a_0^2$ は長さ$^{-2}$の次元をもつ量，$N$ は規格化定数である．$\phi_T$ は $\theta,\varphi$ に依存しない，すなわち球対称であることに注意してほしい．1s軌道を狙っているのだから，これでよい．

やることは，変分原理（第3章）そのままである．エネルギー期待値

$$ \begin{equation} \braket{E}(\alpha)=\frac{\bra{\phi_T}\Ham\ket{\phi_T}}{\braket{\phi_T|\phi_T}} \label{eq:6-energy-def} \end{equation} $$

を $\alpha$ の関数として求め，これを最小にする $\alpha$ を探す．変分原理により，どんな $\alpha$ を選んでも $\braket{E}(\alpha)\ge \epsilon_{1s}$ であるから，$\braket{E}$ を小さくすればするほど真の基底状態に近づく．

なぜ？：分母で割るのを忘れてはいけない

式\eqref{eq:6-energy-def}の分母 $\braket{\phi_T|\phi_T}$ は「規格化されていない試行関数を使うための保険」である．$\phi_T$ をあらかじめ規格化しておけば分母は 1 になるが，そのためには $N$ を $\alpha$ の関数として先に決めておかねばならず，$\alpha$ 微分のときに $N(\alpha)$ の微分まで面倒を見る羽目になる．分母を残したまま計算すると，最後に $N^2$ が分子・分母で綺麗に約分されて消える．規格化定数を決めずに済むのが，比の形にしておく実利である．

図6.1 (a) 水素原子の座標設定．試行関数が球対称なので，角度積分は $4\pi$ を出すだけで済む．(b) Coulombポテンシャルの井戸と厳密な1s準位，そして本章で使うGauss型試行関数の概形．Gauss関数は原点で「丸い山」になっており，厳密解のような尖りがない．この差が本章の主題である．

6.1.3 角度積分を先に済ませる

ブラケットを積分に書き戻す．3次元の体積要素は極座標で $\dd^3r=r^2\sin\theta\,\dd r\,\dd\theta\,\dd\varphi$ である（付録A）．分子は

$$ \begin{align} \bra{\phi_T}\Ham\ket{\phi_T} &=\int_0^{2\pi}\!\!\int_0^{\pi}\!\!\int_0^{\infty}\phi_T(r)\,\Ham\,\phi_T(r)\,r^2\sin\theta\,\dd r\,\dd\theta\,\dd\varphi \nonumber\\ &=\int_0^{2\pi}\dd\varphi\int_0^{\pi}\sin\theta\,\dd\theta\int_0^{\infty}\phi_T(r)\,\Ham\,\phi_T(r)\,r^2\,\dd r \nonumber\\ &=2\pi\cdot 2\cdot\int_0^{\infty}\phi_T(r)\,\Ham\,\phi_T(r)\,r^2\,\dd r =4\pi\int_0^{\infty}\phi_T(r)\,\Ham\,\phi_T(r)\,r^2\,\dd r \label{eq:6-angular} \end{align} $$

となる．$\phi_T$ が実関数なので複素共役は書かなくてよい．$\int_0^{2\pi}\dd\varphi=2\pi$，$\int_0^{\pi}\sin\theta\,\dd\theta=[-\cos\theta]_0^{\pi}=1-(-1)=2$ を使った．以後，この $4\pi$ は分子にも分母にも共通に現れるので，実は最後に約分される．それでも律儀に書いておく．

ハミルトニアン\eqref{eq:6-hamiltonian}を代入すると，計算すべきものは

$$ 4\pi N^2\int_0^{\infty}e^{-ar^2} \underbrace{\left(-\frac{\hbar^2}{2m}\nabla^2\right)}_{\text{第1項：運動エネルギー}} \!\!\!e^{-ar^2}\,r^2\dd r \;+\; 4\pi N^2\int_0^{\infty}e^{-ar^2} \underbrace{\left(-\frac{\ee^2}{4\pi\eps}\frac{1}{r}\right)}_{\text{第2項：Coulomb引力}} \!\!\!e^{-ar^2}\,r^2\dd r $$

の2つと，分母の規格化積分の合計3つである．ここからは面倒だが，一歩ずつ確実にやれば必ず終わる．順に片づけよう．

例題6.1　規格化定数 $N$

試行関数\eqref{eq:6-trial}を規格化する定数 $N$ を求めよ．ただし $\displaystyle\int_0^\infty r^2e^{-cr^2}\dd r=\frac{1}{4}\sqrt{\frac{\pi}{c^3}}$ を用いてよい．

解答　規格化条件は $\displaystyle\int\abs{\phi_T}^2\dd^3r=1$ である．角度積分を先に済ませると

$$ 1=4\pi N^2\int_0^\infty e^{-2ar^2}r^2\,\dd r =4\pi N^2\cdot\frac{1}{4}\sqrt{\frac{\pi}{(2a)^3}} =\pi N^2\frac{\sqrt{\pi}}{(2a)^{3/2}} . $$

したがって $N^2=(2a)^{3/2}/\pi^{3/2}=(2a/\pi)^{3/2}$，すなわち

$$ N=\left(\frac{2a}{\pi}\right)^{3/4} =\left(\frac{2\alpha}{\pi a_0^2}\right)^{3/4}. $$

これは第5章でGauss型軌道の規格化定数として求めたものと同じである．ただし本章の本計算では，6.1.2で述べたとおり $N$ を決めずに比の形で押し通す．$N$ が必要になるのは，最後に波動関数の形を図示するときだけである．

6.2 ラプラシアンの動径部分と運動エネルギー項

6.2.1 球対称関数に働くラプラシアン

直交座標でのラプラシアンは $\nabla^2=\dfrac{\partial^2}{\partial x^2}+\dfrac{\partial^2}{\partial y^2}+\dfrac{\partial^2}{\partial z^2}$ である．極座標に変換すると（付録A）

$$ \nabla^2=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial}{\partial r}\right) +\frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right) +\frac{1}{r^2\sin^2\theta}\frac{\partial^2}{\partial\varphi^2} $$

となるが，$\phi_T$ は $\theta,\varphi$ に依存しないので後ろの2項は消える．残るのは動径部分だけである．

数学ノート：動径ラプラシアンの2つの書き方

球対称な関数 $f(r)$ に対して

$$ \begin{equation} \nabla^2 f=\frac{1}{r^2}\frac{\dd}{\dd r}\left(r^2\frac{\dd f}{\dd r}\right) =\frac{\dd^2f}{\dd r^2}+\frac{2}{r}\frac{\dd f}{\dd r} \label{eq:6-laplacian-radial} \end{equation} $$

が成り立つ．等号の確認は積の微分をするだけである：

$$ \frac{1}{r^2}\frac{\dd}{\dd r}\left(r^2\frac{\dd f}{\dd r}\right) =\frac{1}{r^2}\left(2r\frac{\dd f}{\dd r}+r^2\frac{\dd^2 f}{\dd r^2}\right) =\frac{2}{r}\frac{\dd f}{\dd r}+\frac{\dd^2 f}{\dd r^2}. $$

右辺の形（$\dd^2/\dd r^2+ (2/r)\dd/\dd r$）のほうが具体的な関数に作用させるときは使いやすい．以後こちらを使う．

6.2.2 $\nabla^2\phi_T$ を求める

$\phi_T=Ne^{-ar^2}$ に\eqref{eq:6-laplacian-radial}を作用させる．まず1階微分は，合成関数の微分により

$$ \frac{\dd\phi_T}{\dd r}=N\cdot e^{-ar^2}\cdot(-2ar)=-2ar\,N e^{-ar^2} $$

である．次に2階微分は，$-2ar$ と $e^{-ar^2}$ の積の微分だから

$$ \frac{\dd^2\phi_T}{\dd r^2} =N\left[(-2a)e^{-ar^2}+(-2ar)\cdot(-2ar)e^{-ar^2}\right] =\left(-2a+4a^2r^2\right)Ne^{-ar^2}. $$

また

$$ \frac{2}{r}\frac{\dd\phi_T}{\dd r}=\frac{2}{r}\cdot(-2ar)Ne^{-ar^2}=-4a\,Ne^{-ar^2} $$

である（$r$ が約分されるので，$r\to0$ でも発散しない）．両者を足すと

$$ \begin{equation} \nabla^2\phi_T=\left(-2a+4a^2r^2-4a\right)Ne^{-ar^2} =\left(-6a+4a^2r^2\right)\phi_T . \label{eq:6-lap-phi} \end{equation} $$

6.2.3 必要なGauss積分をつくる

これから $\int_0^\infty r^2e^{-2ar^2}\dd r$ と $\int_0^\infty r^4e^{-2ar^2}\dd r$ が要る．第5章で使った基本形

$$ I_0(c)\equiv\int_0^\infty e^{-cx^2}\dd x=\frac{1}{2}\sqrt{\frac{\pi}{c}} $$

から，パラメータ $c$ で微分するという一手で全部つくれる．

数学ノート：Gauss積分を微分で増やす

$I_0(c)=\dfrac{1}{2}\sqrt{\pi}\,c^{-1/2}$ の両辺を $c$ で微分する．左辺では積分記号の下で微分してよく（被積分関数は $c\gt0$ で十分速く減衰する），

$$ \frac{\dd I_0}{\dd c}=\int_0^\infty\frac{\partial}{\partial c}e^{-cx^2}\dd x =-\int_0^\infty x^2e^{-cx^2}\dd x, \qquad \frac{\dd}{\dd c}\left(\frac{\sqrt{\pi}}{2}c^{-1/2}\right)=-\frac{\sqrt{\pi}}{4}c^{-3/2} $$

だから

$$ \begin{equation} \int_0^\infty x^2e^{-cx^2}\dd x=\frac{\sqrt{\pi}}{4}c^{-3/2}=\frac{1}{4}\sqrt{\frac{\pi}{c^3}} . \label{eq:6-gauss-2} \end{equation} $$

もう一度微分すれば

$$ \begin{equation} \int_0^\infty x^4e^{-cx^2}\dd x=\frac{3\sqrt{\pi}}{8}c^{-5/2}=\frac{3}{8}\sqrt{\frac{\pi}{c^5}} . \label{eq:6-gauss-4} \end{equation} $$

一般には $\displaystyle\int_0^\infty x^{2n}e^{-cx^2}\dd x=\frac{(2n-1)!!}{2(2c)^{n}}\sqrt{\frac{\pi}{c}}$ である（$(2n-1)!!=1\cdot3\cdot5\cdots(2n-1)$，$(-1)!!=1$）．

奇数乗のほうは置換 $t=x^2$ で初等的に出る．$\dd t=2x\,\dd x$ より

$$ \begin{equation} \int_0^\infty x^{2m+1}e^{-cx^2}\dd x=\frac{1}{2}\int_0^\infty t^{m}e^{-ct}\dd t=\frac{m!}{2c^{m+1}} . \label{eq:6-gauss-odd} \end{equation} $$

特に $m=0$ のとき $\displaystyle\int_0^\infty xe^{-cx^2}\dd x=\frac{1}{2c}$ である．講義スライドで「各自，試してご覧」と書かれていたのは，この\eqref{eq:6-gauss-odd}のことである．

本章で使うのは $c=2a$ の場合である．書き下しておこう．

$$ \begin{align} \int_0^\infty r^2e^{-2ar^2}\dd r&=\frac{1}{4}\sqrt{\frac{\pi}{(2a)^3}}=\frac{\sqrt{\pi}}{8\sqrt{2}\,a^{3/2}}, \label{eq:6-int-r2}\\[4pt] \int_0^\infty r^4e^{-2ar^2}\dd r&=\frac{3}{8}\sqrt{\frac{\pi}{(2a)^5}}=\frac{3\sqrt{\pi}}{32\sqrt{2}\,a^{5/2}}, \label{eq:6-int-r4}\\[4pt] \int_0^\infty r\,e^{-2ar^2}\dd r&=\frac{1}{4a}. \label{eq:6-int-r1} \end{align} $$

ここで $(2a)^{3/2}=2\sqrt{2}\,a^{3/2}$，$(2a)^{5/2}=4\sqrt{2}\,a^{5/2}$ を使った．$4\cdot2\sqrt2=8\sqrt2$，$8\cdot4\sqrt2=32\sqrt2$ である．

6.2.4 運動エネルギーの期待値（規格化前）

いよいよ第1項を計算する．\eqref{eq:6-lap-phi}を代入して

$$ \begin{align} \bra{\phi_T}\frac{\hat{p}^2}{2m}\ket{\phi_T} &=4\pi N^2\int_0^\infty e^{-ar^2}\left(-\frac{\hbar^2}{2m}\nabla^2\right)e^{-ar^2}\,r^2\dd r \nonumber\\ &=4\pi N^2\int_0^\infty e^{-ar^2}\left(-\frac{\hbar^2}{2m}\right)\left(-6a+4a^2r^2\right)e^{-ar^2}\,r^2\dd r \nonumber\\ &=4\pi N^2\left(-\frac{\hbar^2}{2m}\right)\int_0^\infty\left(-6a+4a^2r^2\right)e^{-2ar^2}\,r^2\dd r \nonumber\\ &=\frac{4\pi\hbar^2N^2}{m}\int_0^\infty\left(3a-2a^2r^2\right)e^{-2ar^2}\,r^2\dd r . \label{eq:6-kin-1} \end{align} $$

最後の行では $-\dfrac{\hbar^2}{2m}\times(-6a+4a^2r^2)=\dfrac{\hbar^2}{m}(3a-2a^2r^2)$ と整理した（符号を2回反転させているので，ここは特に慎重に）．

積分を\eqref{eq:6-int-r2}\eqref{eq:6-int-r4}で置き換える．

$$ \begin{align} \int_0^\infty\left(3a-2a^2r^2\right)e^{-2ar^2}r^2\dd r &=3a\int_0^\infty r^2e^{-2ar^2}\dd r-2a^2\int_0^\infty r^4e^{-2ar^2}\dd r \nonumber\\ &=3a\cdot\frac{\sqrt{\pi}}{8\sqrt2\,a^{3/2}}-2a^2\cdot\frac{3\sqrt{\pi}}{32\sqrt2\,a^{5/2}} \nonumber\\ &=\frac{3\sqrt{\pi}}{8\sqrt2\,a^{1/2}}-\frac{6\sqrt{\pi}}{32\sqrt2\,a^{1/2}} \nonumber\\ &=\frac{\sqrt{\pi}}{\sqrt{a}}\left(\frac{3}{8\sqrt2}-\frac{3}{16\sqrt2}\right) =\frac{\sqrt{\pi}}{\sqrt{a}}\cdot\frac{3}{16\sqrt2} =\frac{3\sqrt2}{32}\sqrt{\frac{\pi}{a}} . \label{eq:6-kin-2} \end{align} $$

括弧の中は $\dfrac{3}{8\sqrt2}-\dfrac{3}{16\sqrt2}=\dfrac{6-3}{16\sqrt2}=\dfrac{3}{16\sqrt2}$ であり，分母を有理化して $\dfrac{3}{16\sqrt2}=\dfrac{3\sqrt2}{32}$ とした．したがって

$$ \begin{equation} \bra{\phi_T}\frac{\hat{p}^2}{2m}\ket{\phi_T} =\frac{4\pi\hbar^2N^2}{m}\cdot\frac{3\sqrt2}{32}\sqrt{\frac{\pi}{a}} =\frac{3\sqrt2\,\pi\hbar^2N^2}{8m}\sqrt{\frac{\pi}{a}} . \label{eq:6-kin-result} \end{equation} $$

係数は $4\pi\cdot\dfrac{3\sqrt2}{32}=\dfrac{12\sqrt2\,\pi}{32}=\dfrac{3\sqrt2\,\pi}{8}$ である．運動エネルギーは正（$a\gt0$）であり，$a$ を大きくする（波動関数を狭くする）と $\sqrt{1/a}$ が小さくなる……ように見えるが，これは規格化していない量であることに注意してほしい．分母で割った後の振る舞いは6.5節で見る．

例題6.2　積分公式の使い方の確認

$\displaystyle\int_0^\infty r^4e^{-2ar^2}\dd r$ を，\eqref{eq:6-gauss-4}を使わずに，\eqref{eq:6-int-r2}を $a$ で微分することによって求めよ．

解答　\eqref{eq:6-int-r2}の両辺を $a$ で微分する．左辺は

$$ \frac{\dd}{\dd a}\int_0^\infty r^2e^{-2ar^2}\dd r=\int_0^\infty r^2(-2r^2)e^{-2ar^2}\dd r=-2\int_0^\infty r^4e^{-2ar^2}\dd r . $$

右辺は $\dfrac{\sqrt\pi}{8\sqrt2}\dfrac{\dd}{\dd a}a^{-3/2}=\dfrac{\sqrt\pi}{8\sqrt2}\left(-\dfrac{3}{2}\right)a^{-5/2}=-\dfrac{3\sqrt\pi}{16\sqrt2}a^{-5/2}$．したがって

$$ \int_0^\infty r^4e^{-2ar^2}\dd r=\frac{1}{2}\cdot\frac{3\sqrt\pi}{16\sqrt2}a^{-5/2}=\frac{3\sqrt{\pi}}{32\sqrt2\,a^{5/2}} , $$

となり\eqref{eq:6-int-r4}と一致する．ひとつのGauss積分を覚えておけば，あとは微分で無限に作れる——これがGauss型基底が計算に強い理由の一端である．

6.3 Coulomb引力の項

第2項は原子核との Coulomb 引力である．角度積分を済ませた形から出発する．

$$ \begin{align} \bra{\phi_T}\left(-\frac{\ee^2}{4\pi\eps}\frac{1}{r}\right)\ket{\phi_T} &=4\pi N^2\int_0^\infty e^{-ar^2}\left(-\frac{\ee^2}{4\pi\eps}\frac{1}{r}\right)e^{-ar^2}\,r^2\dd r \nonumber\\ &=4\pi N^2\left(-\frac{\ee^2}{4\pi\eps}\right)\int_0^\infty e^{-2ar^2}\,\frac{r^2}{r}\,\dd r \nonumber\\ &=-\frac{N^2\ee^2}{\eps}\int_0^\infty r\,e^{-2ar^2}\dd r . \label{eq:6-coul-1} \end{align} $$

ここで $4\pi\times\dfrac{1}{4\pi\eps}=\dfrac{1}{\eps}$ となって $4\pi$ が綺麗に消えたことに注目してほしい．SI単位系の $1/(4\pi\eps)$ という一見わずらわしい係数は，球対称な積分の $4\pi$ とちょうど打ち消し合うために置かれているのである．

残った積分は\eqref{eq:6-int-r1}そのもの，あるいは直接

$$ \int_0^\infty r\,e^{-2ar^2}\dd r =\left[-\frac{1}{4a}e^{-2ar^2}\right]_0^\infty =0-\left(-\frac{1}{4a}\right)=\frac{1}{4a} $$

である（$\dfrac{\dd}{\dd r}e^{-2ar^2}=-4ar\,e^{-2ar^2}$ を使えば，原始関数が $-\dfrac{1}{4a}e^{-2ar^2}$ であることは一目で分かる）．したがって

$$ \begin{equation} \bra{\phi_T}\left(-\frac{\ee^2}{4\pi\eps}\frac{1}{r}\right)\ket{\phi_T} =-\frac{N^2\ee^2}{4\eps\,a} . \label{eq:6-coul-result} \end{equation} $$

6.4 規格化積分（分母）

最後に式\eqref{eq:6-energy-def}の分母を計算する．角度積分は分子と全く同じで $4\pi$ が出る．

$$ \begin{align} \braket{\phi_T|\phi_T} &=\int_0^{2\pi}\!\!\int_0^\pi\!\!\int_0^\infty \phi_T(r)\phi_T(r)\,r^2\sin\theta\,\dd r\,\dd\theta\,\dd\varphi =4\pi\int_0^\infty N^2e^{-2ar^2}r^2\,\dd r \nonumber\\ &=4\pi N^2\cdot\frac{\sqrt\pi}{8\sqrt2\,a^{3/2}} =\frac{\pi^{3/2}N^2}{2\sqrt2\,a^{3/2}} =\frac{\sqrt2}{4}\left(\frac{\pi}{a}\right)^{3/2}N^2 . \label{eq:6-norm-result} \end{align} $$

ここで $\dfrac{4\pi}{8\sqrt2}=\dfrac{\pi}{2\sqrt2}$ を使い，さらに $\dfrac{1}{2\sqrt2}=\dfrac{\sqrt2}{4}$ と有理化した．$\pi\cdot\sqrt\pi=\pi^{3/2}$ である．

6.5 $E(\alpha)=A\alpha-B\sqrt{\alpha}$ の最小化

6.5.1 期待値を組み立てる

\eqref{eq:6-kin-result}\eqref{eq:6-coul-result}\eqref{eq:6-norm-result}を式\eqref{eq:6-energy-def}に入れる．

$$ \braket{E}= \dfrac{\dfrac{3\sqrt2\,\pi\hbar^2}{8m}\sqrt{\dfrac{\pi}{a}}\,N^2-\dfrac{\ee^2}{4\eps a}N^2} {\dfrac{\sqrt2}{4}\left(\dfrac{\pi}{a}\right)^{3/2}N^2} $$

まず $N^2$ が分子・分母で約分される．予告どおり，規格化定数は決めなくてよかったのである．次に2つの項を別々に処理する．

運動エネルギー項：

$$ \begin{align} \frac{\dfrac{3\sqrt2\,\pi\hbar^2}{8m}\left(\dfrac{\pi}{a}\right)^{1/2}}{\dfrac{\sqrt2}{4}\left(\dfrac{\pi}{a}\right)^{3/2}} &=\frac{3\sqrt2\,\pi\hbar^2}{8m}\cdot\frac{4}{\sqrt2}\cdot\left(\frac{\pi}{a}\right)^{1/2-3/2} \nonumber\\ &=\frac{3\pi\hbar^2}{2m}\cdot\left(\frac{\pi}{a}\right)^{-1} =\frac{3\pi\hbar^2}{2m}\cdot\frac{a}{\pi} =\frac{3\hbar^2}{2m}\,a . \label{eq:6-T-final} \end{align} $$

Coulomb項：

$$ \begin{align} \frac{-\dfrac{\ee^2}{4\eps a}}{\dfrac{\sqrt2}{4}\left(\dfrac{\pi}{a}\right)^{3/2}} &=-\frac{\ee^2}{4\eps a}\cdot\frac{4}{\sqrt2}\cdot\left(\frac{a}{\pi}\right)^{3/2} =-\frac{\ee^2}{\sqrt2\,\eps}\cdot\frac{a^{3/2}}{a\,\pi^{3/2}} =-\frac{\ee^2}{\sqrt2\,\pi^{3/2}\eps}\sqrt{a} . \label{eq:6-V-final} \end{align} $$

まとめると，次元付きの $a$ を使った表式は

$$ \begin{equation} \braket{E}(a)=\frac{3\hbar^2}{2m}\,a-\frac{\ee^2}{\sqrt2\,\pi^{3/2}\eps}\sqrt{a} \label{eq:6-energy-a} \end{equation} $$

である．ここに $a=\alpha/a_0^2$ を代入すれば，無次元パラメータ $\alpha$ の関数として

$$ \begin{equation} \braket{E}(\alpha)=A\,\alpha-B\sqrt{\alpha}, \qquad A=\frac{3\hbar^2}{2ma_0^2}, \qquad B=\frac{1}{a_0}\cdot\frac{\ee^2}{\eps}\cdot\frac{1}{\sqrt2\,\pi^{3/2}} \label{eq:6-energy-alpha} \end{equation} $$

を得る．$\sqrt{a}=\sqrt{\alpha}/a_0$ であることに注意（$a_0$ が1乗で $B$ の分母に入る）．

6.5.2 $A$ と $B$ を $\epsilon_{1s}$ で表す

$A$ と $B$ には物理定数が山ほど入っていて見通しが悪い．ここで一手間かけると，驚くほど綺麗になる．使うのはBohr半径と1s準位の定義（第2章）である．

$$ a_0=\frac{4\pi\eps\hbar^2}{m\ee^2}, \qquad \epsilon_{1s}=-\frac{\ee^2}{8\pi\eps a_0} \quad\Longleftrightarrow\quad \abs{\epsilon_{1s}}=\frac{\ee^2}{8\pi\eps a_0}=13.606\ \mathrm{eV} $$

導出：$A=3\abs{\epsilon_{1s}}$，$B=\dfrac{8}{\sqrt{2\pi}}\abs{\epsilon_{1s}}$

(i) $A$ について．Bohr半径の定義を $\hbar^2/m$ について解くと

$$ a_0=\frac{4\pi\eps\hbar^2}{m\ee^2} \quad\Longrightarrow\quad \frac{\hbar^2}{m}=\frac{a_0\ee^2}{4\pi\eps} \quad\Longrightarrow\quad \frac{\hbar^2}{ma_0^2}=\frac{\ee^2}{4\pi\eps a_0}=2\abs{\epsilon_{1s}} . $$

最後の等号は $\abs{\epsilon_{1s}}=\ee^2/(8\pi\eps a_0)$ の2倍が $\ee^2/(4\pi\eps a_0)$ であることによる．したがって

$$ A=\frac{3}{2}\cdot\frac{\hbar^2}{ma_0^2}=\frac{3}{2}\cdot2\abs{\epsilon_{1s}}=3\abs{\epsilon_{1s}} . $$

(ii) $B$ について．$\dfrac{\ee^2}{\eps a_0}=8\pi\abs{\epsilon_{1s}}$ だから

$$ B=\frac{8\pi\abs{\epsilon_{1s}}}{\sqrt2\,\pi^{3/2}} =\frac{8\abs{\epsilon_{1s}}}{\sqrt2\,\sqrt\pi} =\frac{8}{\sqrt{2\pi}}\abs{\epsilon_{1s}} =4\sqrt{\frac{2}{\pi}}\,\abs{\epsilon_{1s}} . $$

（$\dfrac{8}{\sqrt2\sqrt\pi}=\dfrac{8}{\sqrt{2\pi}}$，また $\dfrac{8}{\sqrt{2\pi}}=\dfrac{8}{\sqrt2\sqrt\pi}=\dfrac{8\sqrt2}{2\sqrt\pi}=4\sqrt{\dfrac{2}{\pi}}$．）

数値を入れておこう．$\abs{\epsilon_{1s}}=13.606\ \mathrm{eV}$ として

$$ \begin{equation} A=3\times13.606=40.82\ \mathrm{eV}, \qquad B=\frac{8}{\sqrt{2\pi}}\times13.606=3.1915\times13.606=43.42\ \mathrm{eV}. \label{eq:6-AB-num} \end{equation} $$

（$\sqrt{2\pi}=2.50663$，$8/2.50663=3.19154$．）以後は

$$ \begin{equation} \braket{E}(\alpha)=3\abs{\epsilon_{1s}}\,\alpha-\frac{8}{\sqrt{2\pi}}\abs{\epsilon_{1s}}\sqrt{\alpha} \label{eq:6-energy-eps} \end{equation} $$

という，$\abs{\epsilon_{1s}}$ を単位にした形で扱う．物理定数がすべて $\abs{\epsilon_{1s}}$ ひとつに押し込められた．

6.5.3 最小値を求める

$\braket{E}(\alpha)=A\alpha-B\sqrt\alpha$ を $\alpha$ で微分する．$\dfrac{\dd}{\dd\alpha}\sqrt\alpha=\dfrac{1}{2\sqrt\alpha}$ だから

$$ \begin{equation} \frac{\dd\braket{E}}{\dd\alpha}=A-\frac{B}{2}\frac{1}{\sqrt\alpha} . \label{eq:6-dEda} \end{equation} $$

これがゼロになる条件は $A=\dfrac{B}{2\sqrt\alpha}$，すなわち $\sqrt\alpha=\dfrac{B}{2A}$ である．両辺を2乗して

$$ \begin{equation} \alpha_{\mathrm{opt}}=\frac{B^2}{4A^2} . \label{eq:6-alpha-opt} \end{equation} $$

この点が最小値であることは，2階微分 $\dfrac{\dd^2\braket{E}}{\dd\alpha^2}=\dfrac{B}{4}\alpha^{-3/2}\gt0$（$B\gt0$，$\alpha\gt0$）から分かる．$\alpha\to0$ でも $\alpha\to\infty$ でも $\braket{E}$ は上がるので，これが唯一の最小点である．

最小値そのものは，$\alpha_{\mathrm{opt}}$ を代入して

$$ \begin{align} E_{\min}=\braket{E}(\alpha_{\mathrm{opt}}) &=A\cdot\frac{B^2}{4A^2}-B\cdot\sqrt{\frac{B^2}{4A^2}} =\frac{B^2}{4A}-B\cdot\frac{B}{2A} \nonumber\\ &=\frac{B^2}{4A}-\frac{B^2}{2A} =\frac{B^2-2B^2}{4A} =-\frac{B^2}{4A} . \label{eq:6-emin} \end{align} $$

6.5.4 数値を入れて検算する

ここが本章のヤマである．丁寧に代入しよう．$A=3\abs{\epsilon_{1s}}$，$B=\dfrac{8}{\sqrt{2\pi}}\abs{\epsilon_{1s}}$ だから，まず

$$ B^2=\frac{64}{2\pi}\abs{\epsilon_{1s}}^2=\frac{32}{\pi}\abs{\epsilon_{1s}}^2, \qquad 4A^2=4\cdot9\abs{\epsilon_{1s}}^2=36\abs{\epsilon_{1s}}^2 . $$

したがって

$$ \begin{equation} \alpha_{\mathrm{opt}}=\frac{B^2}{4A^2} =\frac{(32/\pi)\abs{\epsilon_{1s}}^2}{36\abs{\epsilon_{1s}}^2} =\frac{32}{36\pi}=\frac{8}{9\pi}=0.28294 . \label{eq:6-alpha-num} \end{equation} $$

$\abs{\epsilon_{1s}}$ が完全に消えて，純粋な数が残った．$9\pi=28.2743$，$8/28.2743=0.282942$ である．次元付きに戻せば $a_{\mathrm{opt}}=\dfrac{8}{9\pi a_0^2}$ となる．

エネルギーの最小値は

$$ \begin{equation} E_{\min}=-\frac{B^2}{4A} =-\frac{(32/\pi)\abs{\epsilon_{1s}}^2}{12\abs{\epsilon_{1s}}} =-\frac{32}{12\pi}\abs{\epsilon_{1s}} =-\frac{8}{3\pi}\abs{\epsilon_{1s}} =\frac{8}{3\pi}\epsilon_{1s} . \label{eq:6-emin-eps} \end{equation} $$

最後の等号は $\epsilon_{1s}=-\abs{\epsilon_{1s}}$（負の量）だから，マイナス符号を $\epsilon_{1s}$ に吸わせたものである．数値は

$$ \begin{equation} \frac{8}{3\pi}=\frac{8}{9.42478}=0.848826, \qquad E_{\min}=0.848826\times(-13.606)=-11.549\ \mathrm{eV}. \label{eq:6-emin-num} \end{equation} $$

定理：1個のGauss関数で解いた水素原子

試行関数 $\phi_T=N\exp[-\alpha(r/a_0)^2]$ を用いた変分計算の結果は

$$ \alpha_{\mathrm{opt}}=\frac{8}{9\pi}=0.28294\ \ (\text{無次元}), \qquad E_{\min}=\frac{8}{3\pi}\epsilon_{1s}=0.8488\,\epsilon_{1s}=-11.549\ \mathrm{eV}. $$

厳密解 $\epsilon_{1s}=-13.606\ \mathrm{eV}$ に対する誤差は

$$ \frac{\abs{\epsilon_{1s}}-\abs{E_{\min}}}{\abs{\epsilon_{1s}}}=1-\frac{8}{3\pi}=1-0.8488=0.1512 \quad\Longrightarrow\quad \textbf{15.1\%} $$

である．変分原理の予告どおり $E_{\min}\gt\epsilon_{1s}$（真の基底状態より上）になっていることを確認しておこう．

数学ノート：原子単位系で書くと2行で終わる

本書はSI単位系で通すが，原子単位系（Hartree単位系．$\hbar=m=\ee=4\pi\eps=1$，長さの単位が $a_0$，エネルギーの単位が $1\ \mathrm{Hartree}=27.2114\ \mathrm{eV}$）を使うと，\eqref{eq:6-energy-a}は

$$ E(a)=\frac{3}{2}a-2\sqrt{\frac{2a}{\pi}} $$

と書ける．$\dfrac{\dd E}{\dd a}=\dfrac32-\sqrt{\dfrac{2}{\pi}}\dfrac{1}{\sqrt a}=0$ より $\sqrt a=\dfrac23\sqrt{\dfrac2\pi}$，よって $a=\dfrac49\cdot\dfrac2\pi=\dfrac{8}{9\pi}\ [a_0^{-2}]$．エネルギーは

$$ E=\frac32\cdot\frac{8}{9\pi}-2\sqrt{\frac2\pi}\cdot\frac23\sqrt{\frac2\pi} =\frac{4}{3\pi}-\frac{8}{3\pi}=-\frac{4}{3\pi}=-0.42441\ \mathrm{Hartree}. $$

$-0.42441\times27.2114=-11.549\ \mathrm{eV}$．厳密解は $\epsilon_{1s}=-0.5\ \mathrm{Hartree}=-13.606\ \mathrm{eV}$ であり，比は $0.42441/0.5=0.84883=8/(3\pi)$ で一致する．SI単位系の計算がどこかで間違っていないかを確かめるときには，この原子単位系の短い計算が心強い検算になる．

エネルギー期待値 ⟨E⟩(α)=Aα−B√α の α 依存性のグラフ．曲線は α→0 で 0 から急に下がり，α=8/(9π)≈0.283 で最小値 −11.549 eV をとり，その後は増加して α=1.2 で約 +1.4 eV になる．水平の赤破線は厳密解 ε₁ₛ=−13.6 eV で，曲線はどこでもこれより上にある． — 図6.2 エネルギー期待値 $\braket{E}(\alpha)=A\alpha-B\sqrt\alpha$ の $\alpha$ 依存性．$\alpha\to0$（波動関数が無限に広がる）では $-B\sqrt\alpha$ が効いて $0$ に近づき，$\alpha$ が大きい（細くなる）と運動エネルギー $A\alpha$ が勝って正になる．両者が釣り合う $\alpha_{\mathrm{opt}}=8/(9\pi)$ で最小値 $-11.549\ \mathrm{eV}$ をとる．曲線は**どこでも破線（厳密解）より上**にあり，変分原理が守られている．最小点付近が平らであることにも注目してほしい（6.6.3項）．

物理的意味：$A\alpha$ と $-B\sqrt\alpha$ のせめぎ合い

$\alpha$ は波動関数の「細さ」を表す．$\alpha$ が大きいほど電子は核の近くに閉じ込められる．

運動エネルギー $A\alpha$（正，$\alpha$ の1乗）：閉じ込めると不確定性関係 $\Delta p\gtrsim\hbar/\Delta x$ により運動量のばらつきが増え，運動エネルギーが上がる．狭い箱に閉じ込められた粒子ほど激しく動く．
Coulomb引力 $-B\sqrt\alpha$（負，$\alpha$ の1/2乗）：核に近づくほど深い井戸の底に入れるので，エネルギーは下がる．

$\alpha$ のべきが違う（1乗と1/2乗）ので，両者は必ずどこかで釣り合う．もし両方が $\alpha$ の1乗なら，係数の大小で「無限に潰れる」か「無限に広がる」かのどちらかになり，原子は存在できない．原子が有限の大きさを持つのは，この指数のずれのおかげである．第2章で $a_0$ を導いたときと同じ物理が，別の顔で現れている．

例題6.3　運動エネルギーとポテンシャルエネルギーの内訳

最適な $\alpha_{\mathrm{opt}}=8/(9\pi)$ における $\braket{T}$ と $\braket{V}$ をそれぞれ求め，$2\braket{T}=-\braket{V}$（ビリアル定理）が成り立つことを確かめよ．また厳密解の値と比べよ．

解答　\eqref{eq:6-energy-eps}の第1項と第2項がそのまま $\braket{T}$ と $\braket{V}$ である．

$$ \braket{T}=A\alpha_{\mathrm{opt}}=3\abs{\epsilon_{1s}}\cdot\frac{8}{9\pi}=\frac{8}{3\pi}\abs{\epsilon_{1s}} =0.8488\times13.606=11.549\ \mathrm{eV}, $$ $$ \braket{V}=-B\sqrt{\alpha_{\mathrm{opt}}}=-B\cdot\frac{B}{2A}=-\frac{B^2}{2A} =-\frac{16}{3\pi}\abs{\epsilon_{1s}}=-23.098\ \mathrm{eV}. $$

（$B^2/(2A)=(32/\pi)\abs{\epsilon_{1s}}^2/(6\abs{\epsilon_{1s}})=(16/3\pi)\abs{\epsilon_{1s}}$．）確かに $2\braket{T}=23.098=-\braket{V}$ であり，$\braket{T}+\braket{V}=11.549-23.098=-11.549=E_{\min}$ である．

厳密解では $\braket{T}=13.606\ \mathrm{eV}$，$\braket{V}=-27.211\ \mathrm{eV}$（原子単位で $0.5$ と $-1.0$ Hartree）である．Gauss試行関数は運動エネルギーもCoulombエネルギーも，どちらも絶対値を15%ずつ過小評価している．ちなみに平均半径は $\braket{r}=\sqrt{2/(\pi a)}=\tfrac32a_0$ となり，厳密解と完全に一致する．ひとつの量が合っているからといって波動関数が正しいとは限らない．なおビリアル定理が成り立つのは偶然ではなく，$\alpha$ が波動関数の「スケール」を変えるパラメータであるとき，変分の停留条件がそのままビリアル定理になるからである．

6.6 なぜ15%も外れるのか — カスプ条件と遠方の減衰

変分原理が約束するのは $E_{\min}\ge\epsilon_{1s}$ だけである．実際に得られた $-11.549\ \mathrm{eV}$ は厳密値より $2.06\ \mathrm{eV}$ も高い．化学結合1本が数 eV であることを思えば致命的な誤差である．原因は試行関数の形そのものにあり，Gauss関数には指数関数型の厳密解に対して2つの構造的な欠陥がある．

6.6.1 核カスプ条件：原点で尖っていなければならない

厳密な1s波動関数\eqref{eq:6-exact-1s}を原点近傍で展開すると

$$ \psi_{1s}(r)=\frac{1}{\sqrt{\pi a_0^3}}\left(1-\frac{r}{a_0}+\frac{1}{2}\frac{r^2}{a_0^2}-\cdots\right) $$

となり，$r$ の1次の項が残る．すなわち $\dfrac{\dd\psi_{1s}}{\dd r}\Big|_{r=0}=-\dfrac{1}{a_0}\psi_{1s}(0)\ne0$ である．$r$ は正の量だから，原点を通って反対側へ抜ける直線に沿って見ると波動関数はV字を逆さにしたような尖りを持つ．これを核カスプ（nuclear cusp）と呼ぶ．

導出：核カスプ条件はSchrödinger方程式が要求する

カスプは経験則ではなく，方程式から必然的に出てくる．原点近傍で球対称な解を

$$ \psi(r)=\psi(0)\left(1+br+O(r^2)\right) $$

と置く．これをSchrödinger方程式 $\left[-\dfrac{\hbar^2}{2m}\nabla^2-\dfrac{Z\ee^2}{4\pi\eps r}\right]\psi=E\psi$ に入れる．動径ラプラシアン\eqref{eq:6-laplacian-radial}を使うと

$$ \nabla^2\psi=\frac{\dd^2\psi}{\dd r^2}+\frac{2}{r}\frac{\dd\psi}{\dd r} =\psi(0)\left[O(1)+\frac{2}{r}\left(b+O(r)\right)\right] =\psi(0)\frac{2b}{r}+O(1) . $$

左辺の $r\to0$ での発散項を集めると

$$ -\frac{\hbar^2}{2m}\cdot\frac{2b\,\psi(0)}{r}-\frac{Z\ee^2}{4\pi\eps r}\psi(0) =-\frac{\psi(0)}{r}\left[\frac{\hbar^2 b}{m}+\frac{Z\ee^2}{4\pi\eps}\right] . $$

右辺 $E\psi$ は $r\to0$ で有限だから，$1/r$ の係数はゼロでなければならない．よって

$$ \frac{\hbar^2b}{m}=-\frac{Z\ee^2}{4\pi\eps} \quad\Longrightarrow\quad b=-\frac{mZ\ee^2}{4\pi\eps\hbar^2}=-\frac{Z}{a_0} $$

（最後に $a_0=4\pi\eps\hbar^2/(m\ee^2)$ を使った）．すなわち

$$ \begin{equation} \left.\frac{1}{\psi}\frac{\dd\psi}{\dd r}\right|_{r\to0}=-\frac{Z}{a_0} \qquad\text{(核カスプ条件，Katoの定理)} \label{eq:6-cusp} \end{equation} $$

である．水素原子（$Z=1$）なら $-1/a_0$．この条件は加藤敏夫が1957年に一般の多電子系について厳密に証明した（文献[4]）．

なぜ？：尖っていないと何が困るのか

カスプ条件は「$-1/r$ の発散を運動エネルギー項の $2b/r$ で打ち消す」ことを要求している．Gauss関数は $\dfrac{\dd}{\dd r}e^{-ar^2}=-2are^{-ar^2}$ より $r=0$ で導関数が必ずゼロになり，$b=0$ である．つまりGauss関数は $1/r$ の発散を打ち消す手段を持たない．

その結果どうなるか．核のすぐそばは Coulomb ポテンシャルが最も深い一等地なのに，Gauss関数はそこで「平ら」になってしまい，電子密度を十分に積み上げられない．エネルギーを稼ぎ損ねるのである．15%の誤差の主犯はここにある．

原点近傍（r/a₀=0〜0.8，縦軸は 0.2 から拡大）の波動関数の比較．厳密解（赤）は r=0 で 0.564 から傾き −1/a₀ の折れ点（カスプ）を持ち，STO-3G（紫，0.455）と STO-1G（青，0.276）は r=0 で接線が水平になる． — 図6.3 原点近傍における波動関数の比較（縦軸は 0.20 から拡大）．厳密解（赤）は $r=0$ で傾き $-1/a_0$ の折れ点（カスプ）を持つ．Gauss関数を使う STO-1G（青）と STO-3G（紫）は導関数が必ずゼロになり，接線が水平である．Gauss関数を何本足しても，この平坦性は原理的に解消しない．

6.6.2 遠方の減衰：Gauss関数は落ちるのが速すぎる

第2の欠陥は $r$ が大きい側にある．束縛状態の波動関数は，ポテンシャルが無視できる遠方では

$$ -\frac{\hbar^2}{2m}\frac{\dd^2u}{\dd r^2}\simeq E\,u\quad(u=r\psi,\ E\lt0) \quad\Longrightarrow\quad u\propto e^{-\kappa r},\quad \kappa=\frac{\sqrt{2m\abs{E}}}{\hbar} $$

という単純指数関数で減衰する．水素の1s では $\kappa=1/a_0$ である．ところがGauss関数の減衰は $e^{-ar^2}$，すなわち指数の中身が $r$ の2乗である．$r$ が大きくなると

$$ \begin{equation} \frac{e^{-ar^2}}{e^{-r/a_0}}=\exp\!\left[-ar^2+\frac{r}{a_0}\right]\xrightarrow[r\to\infty]{}0 \label{eq:6-tail} \end{equation} $$

となり，Gauss関数のほうが圧倒的に速くゼロになる．第5章で「Gauss関数は裾が短い」と述べたことの定量的な意味がこれである．

波動関数の遠方での減衰の片対数プロット（r/a₀=0〜8）．厳密解（赤）は直線，STO-1G（青）は放物線で r≳2a₀ から急速に下へ曲がり，r=8a₀ で 4×10⁻⁹ になる．STO-3G（紫）は r≲6a₀ まで厳密解にほぼ重なる． — 図6.4 波動関数の遠方での減衰（片対数プロット）．縦軸が対数なので $e^{-r/a_0}$ の厳密解（赤）は直線，$e^{-\alpha r^2}$ の STO-1G（青）は**放物線**になり，$r\gtrsim2a_0$ から急速に下へ曲がる．$r=8a_0$ では厳密解より5桁も小さい．STO-3G（紫）は最も広がった成分（$\alpha_1=0.1098$）のおかげで $r\lesssim6a_0$ まで厳密解にほぼ重なる．分子間相互作用や負イオンなど裾が効く問題で基底関数を増やす必要があるのはこのためである．

6.6.3 変分最適な $\alpha$ と重なり最大の $\alpha$ はなぜ違うのか

第5章では，Slater型軌道とGauss関数の重なり積分を最大にするという基準でSTO-1Gの指数を決め，$\alpha_S=0.27095$ を得た．本章ではエネルギーを最小にする基準で $\alpha_E=8/(9\pi)=0.28294$ を得た．同じ「最良の1個のGauss関数」を求めているのに，答えが 4.4% ずれている．

なぜ？：「最良」は基準によって変わる

2つの基準は，$r$ 空間のどこを重視するかが違う．

重なり最大：$S=\int\psi_{1s}\phi_T\,\dd^3r=4\pi\int_0^\infty\psi_{1s}\phi_T r^2\dd r$．体積要素の $r^2$ の重みがあるので，$r\sim1$〜$2a_0$ の「電子が最も多くいる領域」で形が合っているかが効く．
エネルギー最小：Coulomb項 $\int\phi_T^2\,r^{-1}\,r^2\dd r$ の重みは $r^1$，運動エネルギー項は波動関数の曲率に敏感である．どちらもより内側を重視する．

エネルギー基準のほうが原点近傍を重く見るので，少しでもカスプ側の損失を取り返そうとして，より収縮した（$\alpha$ の大きい）関数を選ぶ．$\alpha_E=0.28294\gt\alpha_S=0.27095$ という大小関係はこれで説明がつく．

では，どちらを使うと結果はどれだけ違うのか．実際に数値を入れてみる．

表6.1 STO-1G の指数 $\alpha$ を2つの基準で決めたときの比較．エネルギーは\eqref{eq:6-energy-eps}，重なりは厳密な1s軌道との $S=\braket{\psi_{1s}|\phi_T}$（両者とも規格化）．
決め方	$\alpha$	$\braket{E}$ (Hartree)	$\braket{E}$ (eV)	重なり $S$
エネルギー最小（本章）	0.282942	$-0.424413$	$-11.5489$	0.97818
重なり最大（第5章）	0.270950	$-0.424218$	$-11.5436$	0.97840
差	4.4%	0.000195	0.0053	0.00022

$\alpha$ が 4.4% も違うのに，エネルギー差はわずか $0.005\ \mathrm{eV}$（誤差全体 $2.06\ \mathrm{eV}$ の 0.26%）である．図6.2で見たとおり最小点の近くで曲線が平らだからである．停留条件 $\dd\braket{E}/\dd\alpha=0$ が成り立つ点のまわりでは，ずれ $\delta\alpha$ に対してエネルギーの変化は

$$ \braket{E}(\alpha_{\mathrm{opt}}+\delta\alpha)-E_{\min} =\frac{1}{2}\left.\frac{\dd^2\braket{E}}{\dd\alpha^2}\right|_{\mathrm{opt}}(\delta\alpha)^2+\cdots $$

と2次でしか効かない．「試行関数が多少ずれていてもエネルギーはそこそこ当たる」という変分法の強みそのものである．裏を返せば，エネルギーが合っていても波動関数が正しいとは限らない．

6.6.4 それでもGauss関数を使い続ける理由

ここまで読むと「Slater型関数を使えばよいではないか」と思うだろう．実際，$\phi=Ne^{-\beta r}$ で同じ変分計算をすると $\beta=1/a_0$ で厳密解に到達する（演習6.1）．それが厳密解そのものなのだから当然である．

それでも量子化学計算がGauss基底を採用するのは，第5章で学んだGauss積の定理（Gaussian product theorem）があるからである．異なる中心 $\bm{R}_A,\bm{R}_B$ に置かれた2つのGauss関数の積は

$$ e^{-a\abs{\bm{r}-\bm{R}_A}^2}e^{-b\abs{\bm{r}-\bm{R}_B}^2} =K\,e^{-(a+b)\abs{\bm{r}-\bm{R}_P}^2}, \qquad \bm{R}_P=\frac{a\bm{R}_A+b\bm{R}_B}{a+b} $$

と，再び1個のGauss関数になる．Slater型関数ではこうはいかず，2中心・3中心・4中心の積分を解析的に処理できない．分子計算では4中心2電子積分が何百万個も必要になるから，この差は決定的である．

水素 1s 軌道の厳密解（赤）と STO-1G（青，α=8/(9π)）・STO-3G（紫）の比較（r/a₀=0〜4）．r≳1a₀ では STO-3G は厳密解とほとんど区別がつかず，差は原点近傍だけで目立つ． — 図6.5 水素1s軌道の厳密解と，STO-1G（Gauss 1本，$\alpha=8/(9\pi)$）・STO-3G（Gauss 3本の縮約）の比較．$r\gtrsim1a_0$ では STO-3G は厳密解とほとんど区別がつかず，差が目立つのは原点近傍だけである．そこがエネルギーを稼ぐ一等地であることが誤差の本質である．重なり積分 $S$ は STO-1G で 0.978，STO-3G で 0.99984．

6.7 STO-3Gを機械にやらせる

6.7.1 STO-3Gの試行関数

6.5節までで手計算した基底は，Gauss関数を1本しか使っていない．これはSTO-1Gと本質的に同じものである．第5章で導いたSTO-3Gなら，水素（$\zeta=1$）に対して

$$ \begin{align} \phi_T^{\text{STO-3G}}(r)=\ &0.444635\left(\frac{2\times0.109818}{\pi}\right)^{3/4}e^{-0.109818\,(r/a_0)^2} \nonumber\\ &+0.535328\left(\frac{2\times0.405771}{\pi}\right)^{3/4}e^{-0.405771\,(r/a_0)^2} \nonumber\\ &+0.154329\left(\frac{2\times2.22766}{\pi}\right)^{3/4}e^{-2.22766\,(r/a_0)^2} \label{eq:6-sto3g} \end{align} $$

という3本のGauss関数の固定した線形結合（縮約, contraction）を使う．3本の役割分担がはっきりしていることに注目してほしい．

表6.2 STO-3G（$\zeta=1$，水素）を構成する3本のGauss関数の役割．$\alpha$ は $a_0^{-2}$ 単位，$c_i=d_i(2\alpha_i/\pi)^{3/4}$ は規格化まで込みの係数．
#	指数 $\alpha_i$	縮約係数 $d_i$	$c_i$	広がりの目安 $1/\sqrt{\alpha_i}$	役割
1	0.109818	0.444635	0.06045	$3.02\,a_0$	広がった裾を作る
2	0.405771	0.535328	0.19397	$1.57\,a_0$	本体（電子の大半）
3	2.22766	0.154329	0.20056	$0.67\,a_0$	原点近くの尖りを演出

3本目の $\alpha_3=2.23$ は $\alpha_1=0.110$ の20倍も鋭く，原点付近に「山」を積み増して図6.3の紫の曲線を $0.455$ まで持ち上げている．逆に緩やかな $\alpha_1$ は，図6.4で $r\lesssim6a_0$ まで裾を保つ働きをする．

6.7.2 面倒な単純作業は機械にやらせる

では，このSTO-3Gで6.2〜6.5節と同じ変分計算を手でやろう——と考えると，たちまち行き詰まる．3本×3本＝9通りのGauss関数の積について，運動エネルギー・Coulomb引力・重なりの3種類，合計27個の積分が要る．分子になれば中心の数も増える．

物理的意味：人間がやるべき仕事，機械がやるべき仕事

望月先生は講義で次のように述べている．「そういう分かりきった面倒な単純作業は，機械にやらせるべきです．人間は，もっと creative なところに時間を使いましょう．」

ここで重要なのは，「面倒だから機械に投げる」のではなく，「分かりきっているから機械に投げる」という点である．6.2〜6.5節を自分の手で最後まで追い切った読者は，計算機の中で何が起きているかを知っている．その上で機械に任せるのと，何も分からずに任せるのとでは，出てきた数字の読み方がまるで違う．エラーが出たとき，答えが怪しいとき，頼れるのは手計算で身につけた感覚だけである．

6.7.3 PySCFで水素原子を計算する

ここからは PySCF（Python-based Simulations of Chemistry Framework，文献[6]）を使う．導入は付録Bにまとめてあるので，環境が整っていない読者は先にそちらを済ませてほしい．最低限，次の3つが確認できていればよい．

$ which python3         # どこに python があるのか
/usr/local/bin/python3
$ echo $SHELL           # Shell の種類は何か
/bin/zsh
$ pip show pyscf        # Python module PySCF があるのか
Name: pyscf
Version: 2.12.0
Summary: PySCF: Python-based Simulations of Chemistry Framework

コード6.1 計算環境の確認（付録B参照）

講義配布の pyscf_tutorial1.zip を展開すると calc_pyscf.py と各種の .xyz ファイルが入っている．水素原子の計算は次の1行である．

$ cd ~/Downloads/pyscf_tutorial1
$ python3 calc_pyscf.py --xyz XYZ_H.xyz --spin 1

コード6.2 水素原子をSTO-3Gで計算する．--xyz は読み込む構造ファイル，--spin は系全体の $2S$（不対電子の数）を指定する．水素原子は電子1個なので $2S=1$．

構造ファイル .xyz の中身はきわめて単純である．

1
Hydrogen atom
H       0.00000000      0.00000000      0.00000000

コード6.3 XYZ_H.xyz．1行目は原子数，2行目はメモ（何を書いてもよい），3行目以降が元素記号と $x,y,z$ 座標（単位はÅ）．分子構造は VESTA などの無料ソフトで可視化できる（付録B）．

実行すると次のような出力が得られる．

=== PySCF Calculation ===
XYZ file     : XYZ_H.xyz
Basis        : sto-3g
Charge       : 0
Spin (2S)    : 1
Method       : UHF
Total Energy : -0.4665818496 Hartree
             : -12.696339 eV
<S^2>        : 0.750000
Multiplicity : 2.0

MO energies (Hartree / eV):
  Alpha MO energies:
  MO  1:  -0.466582     -12.6963
  Beta  MO energies:
  MO  1:  -0.466582     -12.6963

コード6.4 水素原子・STO-3G の計算結果

数学ノート：出力の読み方

Total Energy：系の全エネルギー．$1\ \mathrm{Hartree}=27.2114\ \mathrm{eV}$ で換算されている（$-0.4665818\times27.2114=-12.6963$）．電子1個の系なので，軌道エネルギーと全エネルギーが一致する．
Method: UHF：使われた計算手法．6.8節で説明する．
<S^2>：全スピン角運動量の2乗の期待値．理論値は $S(S+1)=\frac12\cdot\frac32=0.75$．ぴったり 0.750000 なので，スピンの状態が正しく2重項（doublet）になっている（第7章）．
Multiplicity：スピン多重度 $2S+1=2$．
MO energies：分子軌道（molecular orbital）のエネルギー．Alpha は上向きスピン，Beta は下向きスピンの軌道である．STO-3G では水素の基底関数が1個だけなので，軌道も1本しか出てこない．

6.7.4 STO-1G と STO-3G の比較

結果が出そろった．手計算のSTO-1Gと，機械にやらせたSTO-3Gを並べてみよう．

表6.3 水素原子の基底状態エネルギー．誤差は厳密解 $\epsilon_{1s}=-0.5\ \mathrm{Hartree}=-13.606\ \mathrm{eV}$ に対する相対値 $1-\abs{E}/\abs{\epsilon_{1s}}$．6-31G以下はPySCFによる計算．
基底関数	Gauss関数の本数	$E$ (Hartree)	$E$ (eV)	誤差
STO-1G（手計算，$\alpha=8/9\pi$）	1	$-0.424413$	$-11.549$	15.1%
STO-3G	3（1本に縮約）	$-0.466582$	$-12.696$	6.69%
6-31G	4（2本に縮約）	$-0.498233$	$-13.558$	0.35%
6-311G	6（3本に縮約）	$-0.499810$	$-13.601$	0.04%
厳密解	—	$-0.500000$	$-13.606$	0

Gauss関数を1本から3本に増やしただけで誤差は 15.1% から 6.69% へ半減し，6-31G では 0.35% まで落ちる．基底関数を大きくすればエネルギーは単調に下がる——これは変分原理の直接の帰結である．基底の集合を広げれば試行関数の自由度が増え，最小値はより真値に近い側へ動くしかない．

例題6.4　誤差の見積もり

表6.3から，STO-1G から STO-3G への改善で「取り戻せたエネルギー」を eV 単位で求めよ．また，STO-3G の残差 $0.910\ \mathrm{eV}$ が化学的にどの程度の大きさか論ぜよ．

解答　STO-1G の誤差は $13.606-11.549=2.057\ \mathrm{eV}$，STO-3G の誤差は $13.606-12.696=0.910\ \mathrm{eV}$．取り戻した分は $2.057-0.910=1.147\ \mathrm{eV}$ で，もとの誤差の 55.8% である．

残差 $0.910\ \mathrm{eV}$ は $87.8\ \mathrm{kJ/mol}$ に相当する（$1\ \mathrm{eV}=96.485\ \mathrm{kJ/mol}$）．C–C単結合の結合エネルギーが約 $348\ \mathrm{kJ/mol}$，水素結合が $10$〜$40\ \mathrm{kJ/mol}$ であることを思えば，これは水素結合の数倍という無視できない大きさである．ただし救いがある．化学で問題になるのは絶対エネルギーではなくエネルギー差（反応熱，結合エネルギー，構造間の安定性の差）であり，似た系どうしで誤差が相殺されることが多い．この誤差相殺（error cancellation）に頼るのが実践的な量子化学計算の作法である．第8章・第10章で具体的に見る．

6.8 RHF・UHF・ROHFの違い

コード6.4の出力には Method : UHF と書かれていた．計算機は黙って「解き方」をひとつ選んでいたのである．理論的背景（Slater行列式・交換エネルギー）は第7章と第9章に譲り，ここでは準位図で見て分かる違いに絞る．

6.8.1 閉殻と開殻，$\bm{\alpha}$スピンと$\bm{\beta}$スピン

電子のスピンの $z$ 成分は $+\hbar/2$ か $-\hbar/2$ の2値しかとれない．前者を$\bm{\alpha}$スピン（up spin），後者を$\bm{\beta}$スピン（down spin）と呼ぶ．Pauliの排他原理により，ひとつの空間軌道には $\bm{\alpha}$ と $\bm{\beta}$ が1個ずつ，合計2個までしか入れない．

定義：閉殻と開殻

閉殻（closed shell）：すべての占有軌道が電子2個で埋まっている状態．$\bm{\alpha}$ 電子の数と $\bm{\beta}$ 電子の数が等しい（$N_\alpha=N_\beta$）．例：He, Ne, N$_2$, CH$_4$, H$_2$O．全スピン $S=0$，多重度 $2S+1=1$（一重項）．

開殻（open shell）：不対電子を持つ状態．$N_\alpha\ne N_\beta$．例：H, Li, O$_2$, ラジカル，多くの遷移金属化合物．$2S=N_\alpha-N_\beta$ が不対電子の数であり，これがPySCFの --spin オプションで指定する値である．

6.8.2 3つの解き方

定義：RHF・UHF・ROHF

RHF（Restricted Hartree–Fock，制限Hartree–Fock法）：$\bm{\alpha}$ 電子と $\bm{\beta}$ 電子が同一の空間軌道を使うと制限する．軌道エネルギーも $\bm{\alpha}$・$\bm{\beta}$ で同じになる．閉殻系（He, N$_2$, CH$_4$ など）に使う．
UHF（Unrestricted Hartree–Fock，非制限Hartree–Fock法）：$\bm{\alpha}$ と $\bm{\beta}$ に別々の空間軌道を許す．軌道エネルギーも別々に出る．開殻系（H, Li, O$_2$ など）に使う．閉殻系にも使える．PySCFのチュートリアルスクリプトでは既定値．
ROHF（Restricted Open-shell Hartree–Fock，開殻制限Hartree–Fock法）：開殻でありながら，$\bm{\alpha}$ と $\bm{\beta}$ の空間軌道（したがってエネルギー準位）を等しくする拘束を課す．二重占有される軌道は共通，不対電子だけが片方のスピンに入る．

図6.6 RHF・UHF・ROHF のエネルギー準位の違い．RHFは閉殻系（6電子）を扱い，$\bm{\alpha}$ と $\bm{\beta}$ は同じ空間軌道・同じエネルギーを共有する．UHFは開殻系（5電子）で両者を独立に最適化するので準位がずれる．ROHFは同じ開殻系を扱いながら，準位が一致するよう拘束を課す．UHFのほうが自由度が大きいので，必ず $E_{\mathrm{UHF}}\le E_{\mathrm{ROHF}}$ となる．

なぜ？：UHFで $\bm{\alpha}$ の準位が下がるのはなぜか

不対電子が $\bm{\alpha}$ 側にあるとき，$\bm{\alpha}$ 電子どうしは交換相互作用によって余分に安定化される．交換項は同じスピンの電子の間にしか働かない（第9章で完全に導出する）．$\bm{\alpha}$ 電子は「仲間が多い」ぶん得をし，$\bm{\beta}$ 電子は損をする．だからUHFでは $\bm{\alpha}$ の占有準位が $\bm{\beta}$ より低くなる．これはスピン分極（spin polarization）と呼ばれ，磁性を持つ物質を記述するには絶対に必要な自由度である．望月先生の研究（文献[12]）のような遷移金属酸化物の計算では，スピン分極を許すかどうかで金属か絶縁体かの答えが変わることさえある．

ROHFはこの自由度をわざと捨てる．捨てる代わりに，波動関数がスピン角運動量の固有関数であることが保証される（UHFでは一般に保証されない．これをスピン汚染（spin contamination）と呼び，<S^2> の値が $S(S+1)$ からずれることで検知できる）．

6.8.3 Li 原子を計算してみる

電子3個のLi原子（電子配置 $1s^2 2s^1$）は，不対電子が1個ある典型的な開殻系である．

$ python3 calc_pyscf.py --xyz XYZ_Li.xyz --spin 1

コード6.5 Li原子をSTO-3G・UHFで計算する．「これだけでできるってすごくないっすか？」（講義より）

Method       : UHF
Total Energy : -7.3155259813 Hartree
             : -199.065601 eV
<S^2>        : 0.750000
Multiplicity : 2.0

MO energies (Hartree / eV):
  Alpha MO energies:            # α は up spin のこと
  MO  1:  -2.369171    -64.4684
  MO  2:  -0.180124     -4.9014
  MO 3-5:  0.130126      3.5409   # 2p，3重縮退
  Beta  MO energies:            # β は down spin のこと
  MO  1:  -2.337858    -63.6164
  MO  2:   0.102253      2.7824
  MO 3-5:  0.190916      5.1951

コード6.6 Li原子・STO-3G・UHFの結果．MO 1 が 1s，MO 2 が 2s，MO 3–5 が3重に縮退した 2p である（STO-3Gの最小基底なので5本しかない）．

読み取れることを並べてみよう．

$\bm{\alpha}$ の 1s（$-64.4684\ \mathrm{eV}$）は $\bm{\beta}$ の 1s（$-63.6164\ \mathrm{eV}$）より $0.85\ \mathrm{eV}$ 低い．不対電子が $\bm{\alpha}$ 側にあるためのスピン分極である．
$\bm{\alpha}$ の 2s は $-4.9014\ \mathrm{eV}$ で占有．$\bm{\beta}$ 側の対応する準位は $+2.7824\ \mathrm{eV}$ と大きく上にあり，空である．
2p は3重に縮退している（$3.5409\ \mathrm{eV}$ が3つ）．原子は球対称なので $p_x,p_y,p_z$ が同じエネルギーになる．当然の結果だが，計算がそれを再現していることは，計算が正しく走っている良い証拠である．

次に，同じ系をROHFで解いてみる．

$ python3 calc_pyscf.py --xyz XYZ_Li.xyz --spin 1 --method rohf

コード6.7 Li原子をROHFで計算する

Method       : ROHF
Total Energy : -7.3155259813 Hartree
             : -199.065601 eV
<S^2>        : 0.750000
Multiplicity : 2.0

MO energies (Hartree / eV):
  MO  1:  -2.353491    -64.0417
  MO  2:  -0.038960     -1.0601
  MO 3-5:  0.160521      4.3680

コード6.8 Li原子・STO-3G・ROHFの結果．$\bm{\alpha}$・$\bm{\beta}$ の区別がなくなり，準位が1組にまとまっている．

注意：全エネルギーは同じでも軌道エネルギーは違う

UHF と ROHF の全エネルギーは $-7.3155259813$ Hartree で完全に一致している．これは STO-3G という最小基底では，Li のスピン分極を表現する余地が実質的にないためである（実際 <S^2> がぴったり 0.750000 で，スピン汚染がゼロ）．しかし軌道エネルギーは違う（UHFの 1s は $-64.47/-63.62$，ROHFは $-64.04$）．軌道エネルギーは計算手法に依存する補助的な量であり，観測量ではないことをここで肝に銘じてほしい．基底関数を 6-31G に拡張すると，UHF は $-7.4312358$ Hartree，ROHF は $-7.4312350$ Hartree となり，わずかに UHF が低くなる（自由度が多いぶん必ずこうなる）．

6.8.4 電荷を変える — Li$^+$ イオン

PySCFでは電荷も自由に指定できる．電子を1個取り去った Li$^+$（電子2個，閉殻）を計算してみよう．閉殻なので --spin 0 とすると，スクリプトは自動的にRHFを選ぶ．

$ python3 calc_pyscf.py --xyz XYZ_Li.xyz --spin 0 --charge 1

コード6.9 Li$^+$ イオンの計算

Charge       : 1
Spin (2S)    : 0
Method       : RHF
Total Energy : -7.1354476290 Hartree
             : -194.165420 eV
<S^2>        : 0.000000
Multiplicity : 1.0

MO energies (Hartree / eV):
  MO  1:  -2.736765    -74.4712
  MO  2:  -0.180035     -4.8990
  MO 3-5: -0.094515     -2.5719

コード6.10 Li$^+$・STO-3G・RHFの結果

Li 原子（UHF・ROHF）と Li⁺ イオン（RHF）の軌道エネルギー準位図．縦軸は −8〜+6 eV と −76〜−62 eV の 2 区間に分けてある．UHF では α（赤）と β（青）の準位がずれ，ROHF では一致する．Li⁺ では 1s が −74.5 eV まで下がる． — 図6.7 Li原子（UHF・ROHF）とLi$^+$イオン（RHF）の軌道エネルギー．縦軸は $-8$〜$+6\ \mathrm{eV}$ と $-76$〜$-62\ \mathrm{eV}$ の2区間に分けてある（途中を省略）．UHFでは $\bm{\alpha}$ と $\bm{\beta}$ の準位がずれ，ROHFでは一致する．Li$^+$ では遮蔽が弱まって全準位が下がり，1s は $-64.04\to-74.47\ \mathrm{eV}$ と $10\ \mathrm{eV}$ 以上沈む．

物理的意味：なぜ電荷を＋にすると準位が下がるのか

Li 原子では，2s 電子は 1s の2個の電子に核電荷 $Z=3$ を遮蔽されている．電子を1個取り去ると，残った電子から見た「実効的な核電荷」が増える．より強く引きつけられるので，軌道エネルギーはすべて下がる．図6.7で 1s が $10\ \mathrm{eV}$ 沈み，2s も $-1.06\to-4.90\ \mathrm{eV}$ と下がっているのはこの効果である．

しかし全エネルギーは上がることに注意してほしい．Li が $-199.066\ \mathrm{eV}$，Li$^+$ が $-194.165\ \mathrm{eV}$ で，$4.900\ \mathrm{eV}$ だけ高い．電子を1個引き剥がすにはエネルギーが要るのだから当然である．「準位が下がる」と「全エネルギーが上がる」は矛盾しない．軌道エネルギーの総和は全エネルギーではない——第9章で詳しく扱う重要な事実である．

例題6.5　Li のイオン化エネルギーとKoopmansの定理

コード6.6とコード6.10から Li の第一イオン化エネルギー $I$ を求めよ．また，Li の $\bm{\alpha}$ 2s 軌道のエネルギー $\varepsilon_{2s}=-4.9014\ \mathrm{eV}$ と比較せよ．実験値は $5.392\ \mathrm{eV}$ である．

解答　イオン化エネルギーは「Li$^+$ と Li の全エネルギーの差」として定義される．

$$ I=E(\mathrm{Li}^+)-E(\mathrm{Li}) =(-7.1354476)-(-7.3155260)=0.1800784\ \mathrm{Hartree}=4.900\ \mathrm{eV}. $$

一方，占有最高軌道（HOMO）のエネルギーは $\varepsilon_{2s}=-0.180124\ \mathrm{Hartree}=-4.9014\ \mathrm{eV}$ であり，

$$ \begin{equation} I\simeq-\varepsilon_{\mathrm{HOMO}} \label{eq:6-koopmans} \end{equation} $$

が $0.0013\ \mathrm{eV}$ の精度で成立している．この関係をKoopmansの定理（Koopmans' theorem，文献[11]）という．「電子を1個抜いても残りの軌道は変わらない」と仮定すれば，抜いた軌道のエネルギー分だけ全エネルギーが上がる，という素朴な論理である．実際には残りの電子が緩和するので厳密には成り立たないが，この例のように最小基底では緩和の余地が小さく，ほぼぴったり合う．

実験値 $5.392\ \mathrm{eV}$ と比べると STO-3G の $4.900\ \mathrm{eV}$ は 9% 低い．6-31G に拡張すると $I=5.327\ \mathrm{eV}$ となり，差は 1.2% まで縮む．基底関数を良くすれば実験に近づくという6.7節と同じ傾向がここでも見られる．

6.8.5 基底関数を 6-31G に変えるとどうなるか

スクリプトのオプションで基底関数を切り替えられる（付録B）．Li について STO-3G と 6-31G を比べてみよう．

表6.4 Li原子・Li$^+$イオンの計算結果の基底関数依存性（PySCF）．基底関数の数 $M$ は空間軌道の本数．実験値はNIST．
系・手法	基底	$M$	全エネルギー (Hartree)	$\varepsilon_{1s}^{\alpha}$ (eV)	$\varepsilon_{2s}^{\alpha}$ (eV)
Li, UHF	STO-3G	5	$-7.3155260$	$-64.468$	$-4.901$
Li, ROHF	STO-3G	5	$-7.3155260$	$-64.042$	$-1.060$
Li, UHF	6-31G	9	$-7.4312358$	$-67.416$	$-5.327$
Li, ROHF	6-31G	9	$-7.4312350$	$-67.192$	$-2.063$
Li$^+$, RHF	STO-3G	5	$-7.1354476$	$-74.471$	—
Li$^+$, RHF	6-31G	9	$-7.2354800$	$-75.990$	—

読み取れる傾向は3つある．

全エネルギーは必ず下がる．Li で $-7.3155\to-7.4312$ Hartree（$3.15\ \mathrm{eV}$ の低下）．変分原理の帰結である．
UHF は必ず ROHF 以下．6-31G では $-7.4312358$（UHF）$\lt-7.4312350$（ROHF）．差はごく小さいが符号は理論どおりである．
イオン化エネルギーは実験に近づく．$I=E(\mathrm{Li}^+)-E(\mathrm{Li})$ は STO-3G で $4.900\ \mathrm{eV}$，6-31G で $5.327\ \mathrm{eV}$．実験値 $5.392\ \mathrm{eV}$ に対する誤差は 9.1% → 1.2% と改善する．

なぜ？：全エネルギーは3 eVも動くのに，イオン化エネルギーは0.4 eVしか動かない

STO-3G から 6-31G への拡張で，Li の全エネルギーは $3.15\ \mathrm{eV}$，Li$^+$ の全エネルギーは $2.72\ \mathrm{eV}$ 下がる．ところが差であるイオン化エネルギーの変化は $0.43\ \mathrm{eV}$ にすぎない．基底関数の不足による誤差の大半が，両者で同じ向きに効いて相殺しているのである．これが例題6.4で触れた誤差相殺である．

だからこそ量子化学計算では，比較する2つの系に必ず同じ基底関数・同じ手法を使わなければならない．片方をSTO-3G，もう片方を6-31Gで計算して差を取れば，その差は物理ではなく基底関数の違いを見ているだけになる．当たり前のことのようだが，実際の研究でも起きる事故である．

6.9 まとめと演習

6.9.1 この章のまとめ

本章では，厳密解を知らないふりをして水素原子を1個のGauss関数で解き，その答えを厳密解と突き合わせた．

設定：ハミルトニアン\eqref{eq:6-hamiltonian}と球対称なGauss試行関数\eqref{eq:6-trial}に対し，期待値\eqref{eq:6-energy-def}を比の形のまま扱えば規格化定数 $N$ は自動的に約分される．
運動エネルギー：動径ラプラシアン\eqref{eq:6-laplacian-radial}から $\nabla^2\phi_T=(-6a+4a^2r^2)\phi_T$（式\eqref{eq:6-lap-phi}）．$6=2\times3$ の $3$ は空間の次元．Gauss積分は基本形を $c$ で微分すればいくらでも作れる．
3つの積分：運動エネルギー\eqref{eq:6-kin-result}，Coulomb引力\eqref{eq:6-coul-result}，規格化\eqref{eq:6-norm-result}．SI単位系の $1/(4\pi\eps)$ は角度積分の $4\pi$ と打ち消し合う．
最小化：$\braket{E}(\alpha)=A\alpha-B\sqrt\alpha$（式\eqref{eq:6-energy-alpha}）．$A=3\abs{\epsilon_{1s}}$，$B=\dfrac{8}{\sqrt{2\pi}}\abs{\epsilon_{1s}}$ と書けば物理定数はすべて $\abs{\epsilon_{1s}}$ に押し込められ，$\dd\braket{E}/\dd\alpha=0$ から

$$ \alpha_{\mathrm{opt}}=\frac{B^2}{4A^2}=\frac{8}{9\pi}=0.28294, \qquad E_{\min}=-\frac{B^2}{4A}=\frac{8}{3\pi}\epsilon_{1s}=0.8488\,\epsilon_{1s}=-11.549\ \mathrm{eV}. $$

誤差15.1%の正体：(i) 核カスプ条件\eqref{eq:6-cusp}をGauss関数は原理的に満たせない（原点で導関数が必ずゼロ）．(ii) 遠方で $e^{-\alpha r^2}$ は速く減衰しすぎる．孤立原子では (i) が支配的だが，分子・固体では (ii) も効く．
基準が違えば最適も違う：変分最適 $\alpha_E=0.28294$ と重なり最大 $\alpha_S=0.27095$ は 4.4% 違うが，エネルギー差は $0.005\ \mathrm{eV}$ にすぎない．停留点のまわりでエネルギーは2次でしか変わらないからである．エネルギーが合っていても波動関数が正しいとは限らない．
機械にやらせる：分かりきった単純作業は計算機に任せ，人間は creative な部分に時間を使う．ただし中身を知った上で任せること．PySCFによるSTO-3Gの結果は $-12.696\ \mathrm{eV}$（誤差6.69%），6-31Gなら $-13.558\ \mathrm{eV}$（誤差0.35%）．
RHF・UHF・ROHF：RHFは閉殻用，UHFは開殻用（スピン分極を許す），ROHFは開殻だが準位を同一に拘束．$E_{\mathrm{UHF}}\le E_{\mathrm{ROHF}}$．Li$^+$ では遮蔽が弱まって全準位が下がるが全エネルギーは上がり，Koopmansの定理 $I\simeq-\varepsilon_{\mathrm{HOMO}}$ が良く成り立つ．

物理的意味：本章で身につけるべき「感覚」

近似計算の結果を見たときに問うべきことは，「何%ずれているか」だけではない．どこがずれているかである．本章のGauss関数は，原点近傍というごく狭い領域の形を間違えただけで，エネルギーを15%外した．逆に平均半径 $\braket{r}$ は完璧に当てた．近似の良し悪しは，求めたい物理量がどの領域の波動関数に敏感かによって決まる．この視点は，密度汎関数法の汎関数を選ぶとき，擬ポテンシャルを選ぶとき，基底関数を選ぶとき——計算科学のあらゆる場面で効いてくる．

次章への橋渡し：本章までは電子1個の系だけを扱った．第7章では電子を2個以上に増やす．そこでPauliの排他原理を波動関数のレベルで正しく書く道具——Slater行列式——が必要になり，$\bm{\alpha}$ と $\bm{\beta}$ が本章のような「別々の箱」ではなく反対称性を通じて絡み合っていることが分かる．本章で顔を出したスピン分極・$\braket{\hat{S}^2}$・交換相互作用は，すべて第7章と第9章で正面から扱う．

6.9.2 演習問題

演習6.1　Slater型試行関数なら厳密解に到達する

試行関数 $\phi_T(r)=Ne^{-\beta r}$（$\beta\gt0$ は長さ$^{-1}$の次元をもつ変分パラメータ）を使って，6.2〜6.5節と同じ変分計算を実行し，

$$ \braket{E}(\beta)=\frac{\hbar^2\beta^2}{2m}-\frac{\ee^2\beta}{4\pi\eps} $$

となることを示せ．さらに $\dd\braket{E}/\dd\beta=0$ から $\beta_{\mathrm{opt}}=1/a_0$，$E_{\min}=\epsilon_{1s}$ を導け（厳密解に到達する）．

ヒント：$\nabla^2\phi_T=\left(\beta^2-\dfrac{2\beta}{r}\right)\phi_T$ をまず示す（$\dfrac{\dd\phi_T}{\dd r}=-\beta\phi_T$，$\dfrac{\dd^2\phi_T}{\dd r^2}=\beta^2\phi_T$，これに $\dfrac{2}{r}\dfrac{\dd\phi_T}{\dd r}=-\dfrac{2\beta}{r}\phi_T$ を足す）．積分は $\displaystyle\int_0^\infty r^ne^{-cr}\dd r=\dfrac{n!}{c^{n+1}}$ を $c=2\beta$，$n=0,1,2$ について使えばよい．分子・分母それぞれで $\int_0^\infty r^2e^{-2\beta r}\dd r=\dfrac{2}{(2\beta)^3}=\dfrac{1}{4\beta^3}$ などを丁寧に．最後に $\hbar^2/(ma_0^2)=2\abs{\epsilon_{1s}}$ と $\ee^2/(4\pi\eps a_0)=2\abs{\epsilon_{1s}}$ を使うと $E_{\min}=\abs{\epsilon_{1s}}-2\abs{\epsilon_{1s}}=-\abs{\epsilon_{1s}}$ が出る．

演習6.2　水素様イオンへの拡張

核電荷 $Ze_0$ をもつ水素様イオン（He$^+$, Li$^{2+}$, …）に対して，同じGauss試行関数\eqref{eq:6-trial}を使ったときの $\alpha_{\mathrm{opt}}$ と $E_{\min}$ を求めよ．厳密解は $E=Z^2\epsilon_{1s}$ である．誤差率は $Z$ に依存するか．

ヒント：ハミルトニアンの Coulomb 項が $-\dfrac{Z\ee^2}{4\pi\eps r}$ になるので，$B$ だけが $Z$ 倍になる（$A$ は変わらない）．\eqref{eq:6-alpha-opt}と\eqref{eq:6-emin}より $\alpha_{\mathrm{opt}}=\dfrac{(ZB)^2}{4A^2}=Z^2\cdot\dfrac{8}{9\pi}$，$E_{\min}=-\dfrac{(ZB)^2}{4A}=Z^2\cdot\dfrac{8}{3\pi}\epsilon_{1s}$．厳密解も $Z^2$ に比例するので，誤差率は $Z$ によらず 15.1% のままである．$Z$ が大きいほど波動関数は収縮する（$\alpha\propto Z^2$）が，形の悪さは相似形のまま残る，と理解するとよい．

演習6.3　最小点の平坦さを定量化する

(a) $\braket{E}(\alpha)=A\alpha-B\sqrt\alpha$ の2階微分を求め，$\alpha=\alpha_{\mathrm{opt}}$ での値を $A$ で表せ．
(b) その結果を使い，$\alpha$ を最適値から $\delta$（相対値，たとえば $\delta=0.044$）だけずらしたときのエネルギー上昇を評価する式を作り，$\alpha_S=0.27095$（$\delta=-0.0424$）の場合の上昇量を求めて表6.1の $0.0053\ \mathrm{eV}$ と比べよ．

ヒント：(a) $\dfrac{\dd^2\braket{E}}{\dd\alpha^2}=\dfrac{B}{4}\alpha^{-3/2}$ に $\sqrt{\alpha_{\mathrm{opt}}}=\dfrac{B}{2A}$ を入れると $\dfrac{2A^3}{B^2}$．(b) Taylor展開の2次項 $\dfrac{A^3}{B^2}\alpha_{\mathrm{opt}}^2\delta^2$ に $\alpha_{\mathrm{opt}}=\dfrac{B^2}{4A^2}$ を入れると $\dfrac{\abs{E_{\min}}}{4}\delta^2$ と驚くほど簡単になる．$\delta=-0.0424$ なら $\dfrac{11.549}{4}\times0.0424^2=0.0052\ \mathrm{eV}$ で表6.1と一致する．

演習6.4　STO-3G の原点での値

表6.2の値を使って，STO-3Gの試行関数\eqref{eq:6-sto3g}の $r=0$ における値 $\phi_T^{\text{STO-3G}}(0)$ を計算し，図6.3の $0.455\,a_0^{-3/2}$ を再現せよ．また，これが厳密解の値 $\psi_{1s}(0)=1/\sqrt{\pi}=0.5642\,a_0^{-3/2}$ の何%かを求めよ．さらに，$\dfrac{\dd\phi_T^{\text{STO-3G}}}{\dd r}\Big|_{r=0}=0$ であることを確かめよ．

ヒント：$r=0$ では指数関数がすべて 1 になるので，$\phi_T(0)=\sum_i d_i(2\alpha_i/\pi)^{3/4}=\sum_i c_i$ である．表6.2の $c_i$ を足すだけでよい（$0.06045+0.19397+0.20056$）．導関数については，各成分が $-2\alpha_i r\,c_ie^{-\alpha_ir^2}$ で $r=0$ を代入するとすべてゼロ．何本足しても $r$ の1次の項は現れない——これがカスプを作れない理由である．

演習6.5　誤差相殺を自分の手で確かめる

Li と Li$^+$ を STO-3G と 6-31G の両方で計算し，(a) それぞれの全エネルギー，(b) イオン化エネルギー $I=E(\mathrm{Li}^+)-E(\mathrm{Li})$ を求めよ．全エネルギーの基底関数依存性と，$I$ の基底関数依存性の大きさを比較し，誤差相殺がどれだけ効いているかを%で述べよ．実験値は $I_{\exp}=5.392\ \mathrm{eV}$ である．

ヒント：表6.4に答えの一部が載っている．全エネルギーは Li で $3.15\ \mathrm{eV}$，Li$^+$ で $2.72\ \mathrm{eV}$ 下がるのに対し，$I$ の変化は $0.43\ \mathrm{eV}$ にすぎない．基底関数の改善による低下量のうち $2.72/3.15=86\%$ が両者で共通に効いて相殺している．絶対エネルギーの誤差が数 eV あっても，エネルギー差の誤差は0.1 eVのオーダーに収まりうる——量子化学計算が実用になる理由の核心をここで体感してほしい．

6.9.3 参考文献

原田義也『量子化学（上巻）』裳華房（1992）
A. Szabo, N. S. Ostlund『新しい量子化学 — 電子構造の理論入門（上）』大野公男・阪井健男・望月祐志訳，東京大学出版会（1987）（原著：Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover, 1996）— 第3章にSTO-$n$Gの図（本章の図6.5に対応）とRHF/UHFの詳しい定式化がある．
W. J. Hehre, R. F. Stewart, J. A. Pople, "Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of Slater-Type Atomic Orbitals", J. Chem. Phys. 51, 2657 (1969). — STO-$n$G の原論文．本章で使った縮約係数・指数の出典．
T. Kato, "On the eigenfunctions of many-particle systems in quantum mechanics", Commun. Pure Appl. Math. 10, 151 (1957). — 核カスプ条件\eqref{eq:6-cusp}の厳密な証明．
S. F. Boys, "Electronic wave functions. I. A general method of calculation for the stationary states of any molecular system", Proc. R. Soc. London A 200, 542 (1950). — Gauss型基底の提唱．
Q. Sun et al., "PySCF: the Python-based simulations of chemistry framework", WIREs Comput. Mol. Sci. 8, e1340 (2018). — 本章で使った計算プログラム．
R. Ditchfield, W. J. Hehre, J. A. Pople, "Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules", J. Chem. Phys. 54, 724 (1971). — 6-31G 基底の原論文．
J. A. Pople, R. K. Nesbet, "Self-Consistent Orbitals for Radicals", J. Chem. Phys. 22, 571 (1954). — UHF の原論文．
C. C. J. Roothaan, "Self-Consistent Field Theory for Open Shells of Electronic Systems", Rev. Mod. Phys. 32, 179 (1960). — ROHF の定式化．
T. Koopmans, "Über die Zuordnung von Wellenfunktionen und Eigenwerten zu den einzelnen Elektronen eines Atoms", Physica 1, 104 (1934). — Koopmansの定理\eqref{eq:6-koopmans}．
P. Atkins, J. de Paula, Atkins' Physical Chemistry, 11th ed., Oxford University Press (2018).
Y. Mochizuki et al., "Theoretical Exploration of the Chemical Bonding in Layered Perovskites by COHP/COBI Analysis", J. Phys. Chem. C 125, 7959 (2021). — 局在基底への射影で固体の化学結合を解析する例．「基底関数の選択」は固体計算の解析にも直結する．