第1章多電子問題と第一原理計算
物質は原子核と電子の集まりであり,その振る舞いは量子力学によって支配されている.したがって「原子番号と電子数」という最小限の情報だけを入力として,Schrödinger方程式を解けば,原理的にはあらゆる物質の性質——結晶構造,結合エネルギー,磁性,光学スペクトル——を実験に頼らず予言できるはずである.これが第一原理計算(first-principles calculation)の思想である.しかし多電子系のSchrödinger方程式は,そのままでは絶望的に解けない.本章では,まず解くべき問題(多電子ハミルトニアン)を厳密に書き下し,以後の全章で用いる単位系(Hartree原子単位)とBorn–Oppenheimer近似を丁寧に導入する.そのうえで「なぜ波動関数を直接扱う方法が破綻するのか」を定量的に見積もり,その困難を回避する諸手法——とりわけ本書の主題である密度汎関数理論(DFT)——の見取り図を与える.次章以降で個々の理論を精密に構築していくための,出発点となる章である.
- 第一原理計算とは何か:入力は原子番号と電子数だけ,という思想と,そこから予測できる物性の範囲
- 多電子系のハミルトニアンをSI単位系で書き下し,5つの項の物理的意味を説明できるようになる
- 水素原子のSchrödinger方程式の無次元化を全ステップ実行し,Bohr半径 $a_0$ とHartreeエネルギー $E_h$ を導出する.原子単位とSI単位の換算表を使えるようになる
- Born–Oppenheimer(断熱)近似の導出:どの項が,なぜ落とされるのかを式の上で特定できるようになる
- 多電子波動関数の基本性質(不可弁別性・反対称性・規格直交性)の起源を理解する
- 「次元の呪い」の定量的見積りと,第一原理計算手法(CI・DFT・QMC・多体グリーン関数法)の分類
1.1 物質科学と第一原理計算
20世紀初頭に完成した量子力学は,物質科学に対してひとつの壮大な約束を与えた.原子核と電子の間に働く力はCoulomb相互作用としてすべて分かっており,その運動を支配する方程式(Schrödinger方程式,相対論的にはDirac方程式)も分かっている.ならば,物質を構成する原子の種類(原子番号)と数さえ指定すれば,あとは方程式を解くだけで物質のあらゆる性質が予言できるはずである.Diracは1929年,量子力学の完成直後にこの状況を見抜き,化学の全体と物理学の大部分を支配する数学的法則はすでに完全に知られており,残る困難はただその方程式が複雑すぎて厳密に解けない点にある,と指摘した[1].つまり物質科学の理論的課題は,この時点で「法則の発見」から「方程式の解法」へと移ったのである.
実験で測定されたパラメータ(格子定数,有効質量,経験ポテンシャルなど)を一切使わず,量子力学の基本法則と物理定数のみから出発して物質の性質を計算する方法を,第一原理計算(first-principles calculation)あるいはab initio計算と呼ぶ.「第一原理」とは,それ以上遡れない基本原理——ここでは量子力学とCoulomb相互作用——だけを仮定するという意味である.入力は本質的に
- 系に含まれる原子核の種類 $\{Z_I\}$(原子番号)と数,
- 電子数 $N_e$(中性系なら $N_e=\sum_I Z_I$),
- 必要に応じて原子核の初期配置 $\{\bm{R}_I\}$(結晶なら格子と基底),
だけであり,調整可能な経験パラメータは存在しない.この「パラメータのなさ」こそが第一原理計算の予言能力の源泉である.同じ枠組み・同じ近似を,半導体にも金属にも分子にも表面にも適用でき,実験に先立って未知物質の性質を計算できる.
1.1.1 何が計算できるか
第一原理計算の最も基本的な出力は系の全エネルギー(total energy)$E$ である.一見地味な量だが,全エネルギーが原子配置の関数 $E(\{\bm{R}_I\})$ として計算できると,そこから驚くほど多くの物性が導かれる.
- 安定構造と格子定数:$E(\{\bm{R}_I\})$ を最小化する配置が平衡構造である.結晶の格子定数,分子の結合長・結合角が,実験情報なしに決まる(現代のDFT計算では格子定数は典型的に実験値の1%程度の精度で再現される).
- 結合・凝集のエネルギー論:分子の解離エネルギー,固体の凝集エネルギー,欠陥形成エネルギー,表面エネルギー,相安定性(多形間のエネルギー差).
- 力学的性質:エネルギーの2階微分から弾性定数,体積弾性率,フォノン(格子振動)分散,さらに熱力学関数.
- 原子に働く力と動力学:$\bm{F}_I=-\partial E/\partial\bm{R}_I$ を使えば構造最適化や第一原理分子動力学が実行でき,化学反応の経路や活性化エネルギーが追跡できる.
- 電子状態そのもの:バンド構造,状態密度,電子密度分布,仕事関数.磁性体ならスピン密度と磁気モーメント,磁気秩序(強磁性か反強磁性か).
- 応答とスペクトル:誘電率,光学吸収,X線分光,振動分光(赤外・Raman),NMRシフトなど,外場に対する応答として定式化される諸量.
物質科学の研究において,これらの計算は大きく3つの役割を果たす.第一に実験の解析・解釈:測定されたスペクトルや構造データの背後にある電子状態を同定する.第二に機構の理解に基づく物質設計:なぜこの物質は硬いのか,なぜ磁性を持つのかという機構を電子論のレベルで理解し,望みの機能を持つ物質の設計指針を得る.第三に物質探索:実験合成に先立って候補物質の安定性や物性を計算でスクリーニングする(機能から構造を逆に求める逆問題的アプローチ).いずれの場合も,方程式は共通であり,図1.1に示す一本の流れ——入力(原子番号)→ 多電子Schrödinger方程式 → 制御された近似 → 物性予測——を通る.
この流れの中で本書が主題とするのは「制御された近似」の段階,とくに密度汎関数理論(density functional theory; DFT)である.DFTは,多電子波動関数の代わりに電子密度 $\rho(\rr)$ という3変数の関数を基本変数に据えることで,多電子問題を実用可能な計算量に落とし込む.その考え方は§1.7で概観し,厳密な定式化は後の章で行う.まずは,解くべき方程式そのものを正確に書き下すことから始めよう.
1.2 多電子系のハミルトニアン
この節では,原子・分子・固体を記述する非相対論的な多体ハミルトニアンをSI単位系で書き下し,各項の物理的意味を確認する.この章の後半以降は原子単位系(§1.3)に移行するが,単位系の切り替えを明確にするため,出発点はあえてSI単位系で書く.
1.2.1 系の設定
考える系は,$N_e$ 個の電子と $N_n$ 個の原子核からなる.記号を次のように定める.
- 電子:位置 $\rr_i$($i=1,\dots,N_e$),質量 $m_e$,電荷 $-e$($e \gt 0$ は素電荷)
- 原子核:位置 $\bm{R}_I$($I=1,\dots,N_n$),質量 $M_I$,電荷 $+Z_I e$($Z_I$ は原子番号)
仮定するのは次の4点である.(i) 非相対論的量子力学(相対論効果は重元素で重要になるが,本書では必要な箇所で補正として言及する).(ii) 粒子間の相互作用は静的なCoulomb相互作用のみ(電磁場の量子化・遅延効果は無視).(iii) 原子核は点電荷(核の内部構造は考えない).(iv) 外部電磁場は当面かけない.これらの仮定の下で,系のすべての情報は波動関数
$$ \Psi(\rr_1,\dots,\rr_{N_e},\bm{R}_1,\dots,\bm{R}_{N_n},t) $$に含まれ,その時間発展は時間依存Schrödinger方程式
\begin{equation} i\hbar\frac{\partial}{\partial t}\Psi = \hat{H}\Psi \label{eq:tdse} \end{equation}で決まる.ここで $\hbar = h/2\pi$ は換算Planck定数である.
1.2.2 ハミルトニアンの5つの項
ハミルトニアン $\hat{H}$ は,運動エネルギーとCoulombポテンシャルエネルギーの和として,SI単位系で次のように書ける.
\begin{equation} \hat{H} = \underbrace{-\sum_{i=1}^{N_e}\frac{\hbar^2}{2m_e}\nabla_i^2}_{\hat{T}_e} \underbrace{-\sum_{I=1}^{N_n}\frac{\hbar^2}{2M_I}\nabla_I^2}_{\hat{T}_n} \underbrace{-\sum_{i=1}^{N_e}\sum_{I=1}^{N_n}\frac{Z_I e^2}{4\pi\varepsilon_0\abs{\rr_i-\bm{R}_I}}}_{\hat{V}_{en}} +\underbrace{\sum_{i \lt j}^{N_e}\frac{e^2}{4\pi\varepsilon_0\abs{\rr_i-\rr_j}}}_{\hat{V}_{ee}} +\underbrace{\sum_{I \lt J}^{N_n}\frac{Z_I Z_J e^2}{4\pi\varepsilon_0\abs{\bm{R}_I-\bm{R}_J}}}_{\hat{V}_{nn}} \label{eq:h-si} \end{equation}ここで $\nabla_i^2 = \partial^2/\partial x_i^2+\partial^2/\partial y_i^2+\partial^2/\partial z_i^2$ は $i$ 番目の電子の座標に関するLaplacian,$\varepsilon_0$ は真空の誘電率である.5つの項の意味を順に確認する.
- $\hat{T}_e$:電子の運動エネルギー.古典力学の $\bm{p}^2/2m_e$ に対応する演算子で,量子力学の処方 $\bm{p}\to-i\hbar\nabla$ により $-\hbar^2\nabla^2/2m_e$ となる.波動関数が空間的に急激に変化する(強く局在する)ほど大きくなる量であり,電子が原子核に「落ち込まない」で有限の広がりを保つのはこの項のためである.物質の安定性と化学結合を理解する鍵となる項で,第2章のビリアル定理で主役を演じる.
- $\hat{T}_n$:原子核の運動エネルギー.形は $\hat{T}_e$ と同じだが,質量 $M_I$ が電子質量の数千倍から数十万倍であるため,係数 $\hbar^2/2M_I$ は圧倒的に小さい.この事実がBorn–Oppenheimer近似(§1.4)の根拠になる.
- $\hat{V}_{en}$:電子–原子核間の引力.電荷 $-e$ と $+Z_Ie$ の間のCoulomb引力なので符号は負である.電子を物質内に束縛する項であり,これがなければ物質は存在しない.電子の側から見ると,原子核の配置 $\{\bm{R}_I\}$ が作る「外場」とみなせるので,後の章では外部ポテンシャル $v_{\mathrm{ext}}(\rr) = -\sum_I Z_I e^2/(4\pi\varepsilon_0\abs{\rr-\bm{R}_I})$ と呼ぶ.DFTの中心概念「外部ポテンシャルと電子密度の一対一対応」の「外部ポテンシャル」とはこの項のことである.
- $\hat{V}_{ee}$:電子間の反発.同符号電荷どうしなので正である.この項こそが多電子問題を「多体問題」たらしめる元凶である.もしこの項がなければ,ハミルトニアンは各電子ごとの項の単純な和になり,問題は1電子問題に分離できてしまう(変数分離).電子間反発があるために,ある電子の運動は他のすべての電子の位置に依存し,波動関数は $3N_e$ 個の変数が絡み合った関数になる.
- $\hat{V}_{nn}$:原子核間の反発.原子核の配置を固定して電子状態を解く文脈(§1.4)では単なる定数として加わるが,構造安定性や原子に働く力を論じるときには本質的な寄与をする.
和の記法について一つ注意しておく.電子間反発の $\sum_{i \lt j}$ は「すべての電子対について1回ずつ足す」という意味である.同じ和は,順序対にわたる和を使って
\begin{equation} \sum_{i \lt j}^{N_e}\frac{e^2}{4\pi\varepsilon_0 \abs{\rr_i-\rr_j}} = \frac{1}{2}\sum_{i=1}^{N_e}\sum_{\substack{j=1 \\ j\neq i}}^{N_e}\frac{e^2}{4\pi\varepsilon_0 \abs{\rr_i-\rr_j}} \label{eq:pair-sym} \end{equation}とも書ける.右辺の制限なし二重和では,同一の電子対 $\{i,j\}$ が $(i,j)$ と $(j,i)$ の2回数えられ,しかも被和項は $i$ と $j$ の入れ替えについて対称(分母の $\abs{\rr_i-\rr_j}$ は入れ替えで不変)なので,2で割れば各対を1回ずつ数えたことになる.文献によって両方の記法が使われるので,係数 $1/2$ の有無に注意されたい.$j=i$ の項(電子の自己相互作用)は物理的に存在しないので和から除かれている——この当たり前の事実が,後にDFTの近似汎関数の「自己相互作用誤差」問題として再登場することを予告しておく.
物理的意味:このハミルトニアンの普遍性
式 \eqref{eq:h-si} の驚くべき点は,その普遍性にある.水素分子も,DNAも,鉄の結晶も,半導体デバイスも,このハミルトニアンの $\{Z_I\}$,$N_e$,$N_n$ を変えるだけですべて記述される.物質ごとに異なる「力の法則」は存在せず,物質の個性はすべて,同一の方程式の解の多様性として現れる.金属と絶縁体の違い,磁性の有無,硬さや色の違いは,方程式に入れた仮定の違いではなく,解である波動関数の性質の違いなのである.だからこそ,この方程式を(近似的にでも)解く一般的な方法が確立すれば,それは特定の物質ではなく物質科学全体に対する道具になる.
1.2.3 定常状態:時間依存方程式から時間に依存しない方程式へ
式 \eqref{eq:tdse} は時間依存の方程式だが,$\hat{H}$ 自身は時間を含まない.この場合,解は「定常状態」の重ね合わせに分解でき,各定常状態は時間に依存しないSchrödinger方程式を満たす.以後の章の出発点となる式なので,変数分離の手続きを省略せずに実行しておこう.
導出:変数分離による定常状態の抽出
全座標(電子と核)をまとめて $\bm{X}=(\rr_1,\dots,\rr_{N_e},\bm{R}_1,\dots,\bm{R}_{N_n})$ と書く.解を空間部分と時間部分の積
$$ \Psi(\bm{X},t) = \psi(\bm{X})\,f(t) $$と仮定して(変数分離の仮定),式 \eqref{eq:tdse} に代入する.左辺の時間微分は $\psi(\bm{X})$ に作用せず,右辺の $\hat{H}$ は空間座標の微分と掛け算だけからなるので $f(t)$ に作用しない.よって
\begin{equation} i\hbar\,\psi(\bm{X})\frac{\dd f(t)}{\dd t} = f(t)\,\hat{H}\psi(\bm{X}) \label{eq:sep1} \end{equation}となる.$\psi(\bm{X})f(t)\neq 0$ の領域で両辺を $\psi(\bm{X})f(t)$ で割ると(ゼロ点は連続性で埋められる),
\begin{equation} i\hbar\,\frac{1}{f(t)}\frac{\dd f(t)}{\dd t} = \frac{\hat{H}\psi(\bm{X})}{\psi(\bm{X})} \label{eq:sep2} \end{equation}を得る.左辺は $t$ のみの関数,右辺は $\bm{X}$ のみの関数である.任意の $t$ と $\bm{X}$ で両者が等しくあり続けるためには,両辺が $t$ にも $\bm{X}$ にも依存しない定数でなければならない(左辺を $t$ で偏微分すればゼロ,右辺を $\bm{X}$ で偏微分してもゼロ,という論法で確かめられる).この分離定数を $E$ と置く.$E$ はエネルギーの次元を持つ.すると式 \eqref{eq:sep2} は2本の方程式に分かれる.時間部分は
\begin{equation} i\hbar\frac{\dd f}{\dd t} = E f \label{eq:f-ode} \end{equation}である.これは1階の線形常微分方程式で,変数分離形 $\dd f/f = -(iE/\hbar)\,\dd t$ に書き直して両辺を積分すると($\int \dd f/f = \ln f$ を使った),
$$ \ln f(t) = -\frac{iE}{\hbar}t + \text{const.} \quad\Longrightarrow\quad f(t) = f(0)\,e^{-iEt/\hbar} $$となる.定数 $f(0)$ は $\psi$ の規格化に吸収できるので $f(0)=1$ と取る.空間部分は
\begin{equation} \hat{H}\psi(\bm{X}) = E\psi(\bm{X}) \label{eq:tise} \end{equation}すなわち $\hat{H}$ の固有値問題である.以上より,時間を含まないハミルトニアンに対する変数分離解は
$$ \Psi(\bm{X},t) = \psi(\bm{X})\,e^{-iEt/\hbar} $$と定まった.この状態では確率密度 $\abs{\Psi}^2 = \abs{\psi}^2\abs{e^{-iEt/\hbar}}^2 = \abs{\psi}^2$ が時間によらないので定常状態と呼ぶ.一般の初期条件に対する解は,固有状態の完全性により定常状態の線形重ね合わせ $\Psi = \sum_n c_n \psi_n e^{-iE_n t/\hbar}$ で表される.
この導出で得られた式 \eqref{eq:tise} が,本書全体で解こうとする方程式——時間に依存しないSchrödinger方程式——である.とくに固有値が最小の状態 $\psi_0$ を基底状態(ground state),その固有値 $E_0$ を基底状態エネルギーと呼ぶ.絶対零度における物質の構造・結合・磁性は基底状態の性質であり,DFTが直接の対象とするのもまず基底状態である.
次節では,式 \eqref{eq:h-si} に現れる $\hbar$,$m_e$,$e$,$4\pi\varepsilon_0$ という定数の群れを一掃する.数式が見やすくなるだけでなく,電子の世界の「自然な物差し」が何かが明らかになる.
1.3 原子単位系(Hartree原子単位)
式 \eqref{eq:h-si} には $\hbar$,$m_e$,$e$,$4\pi\varepsilon_0$ という4つの定数が繰り返し現れる.SI単位系はメートル・キログラム・秒という人間サイズの物差しで作られているため,原子の世界を記述すると $\hbar = 1.054572\times10^{-34}\,\mathrm{J\cdot s}$ のような極端な数値を引きずり続けることになる.そこで発想を逆転し,電子にとって自然な物差しで単位系を組み直す.これがHartree原子単位系(Hartree atomic units; 以下たんに原子単位,a.u.)である.
定義:Hartree原子単位系
次の4つの定数の値を1と定める単位系をHartree原子単位系という.
$$ \hbar = 1,\qquad m_e = 1,\qquad e = 1,\qquad 4\pi\varepsilon_0 = 1 $$すなわち,作用の単位を $\hbar$,質量の単位を電子質量 $m_e$,電荷の単位を素電荷 $e$,Coulomb力の結合定数 $k_e=1/4\pi\varepsilon_0$ を1に取る.長さ・エネルギー・時間などの単位はこれら4つから組み立てられる導出単位であり,長さの単位はBohr半径 $a_0$,エネルギーの単位はHartreeエネルギー $E_h$ となる(本節で導出する).本書では以後,断りのない限りこの単位系を用いる.
「定数を1と置く」という操作はしばしば天下り的に導入されるが,実際には方程式の無次元化という系統的な手続きの帰結である.もっとも簡単な多電子系の「原型」である水素原子を例に,この手続きを1ステップずつ実行してみよう.すると Bohr半径とHartreeエネルギーが,人為的な定義としてではなく,方程式自身が持っている固有のスケールとして自然に浮かび上がる.
1.3.1 水素原子のSchrödinger方程式の無次元化
これから行うことは次の通りである:水素原子の時間に依存しないSchrödinger方程式(SI単位系)に対し,長さを未知のスケール $a$ で測り直す変数変換を施し,方程式を無次元の形に整理する.その際「係数がすべて1になる」ように $a$ を選ぶと,その $a$ がBohr半径に他ならないことを示す.
導出:無次元化によるBohr半径とHartreeエネルギーの出現
ステップ0(出発点).水素原子は,原点に固定した陽子(電荷 $+e$)のCoulomb場の中の1個の電子である(陽子の運動は §1.4 の精神で無視する.厳密には換算質量を使うべきだが,その補正は $m_e/M_p\sim 5\times10^{-4}$ の程度で,ここでの議論には影響しない).時間に依存しないSchrödinger方程式は,SI単位系で
\begin{equation} \left[-\frac{\hbar^2}{2m_e}\nabla^2 - \frac{e^2}{4\pi\varepsilon_0\, r}\right]\psi(\rr) = E\,\psi(\rr) \label{eq:hatom-si} \end{equation}である.$r=\abs{\rr}$ は原点(陽子)からの距離である.
ステップ1(長さのスケール変換).長さの次元を持つ未定の定数 $a$ を導入し,無次元の座標 $\tilde{\rr}$ を
$$ \rr = a\,\tilde{\rr} \qquad\Longleftrightarrow\qquad \tilde{\rr} = \rr/a $$で定義する.微分演算子がどう変換されるかを,$x$ 成分について連鎖則で確かめる.$\tilde{x} = x/a$ なので $\partial\tilde{x}/\partial x = 1/a$ であり,
\begin{align} \frac{\partial}{\partial x} &= \frac{\partial \tilde{x}}{\partial x}\frac{\partial}{\partial \tilde{x}} = \frac{1}{a}\frac{\partial}{\partial \tilde{x}}, \label{eq:chain1}\\[2pt] \frac{\partial^2}{\partial x^2} &= \frac{1}{a}\frac{\partial}{\partial \tilde{x}}\left(\frac{1}{a}\frac{\partial}{\partial \tilde{x}}\right) = \frac{1}{a^2}\frac{\partial^2}{\partial \tilde{x}^2} \label{eq:chain2} \end{align}となる(2行目では,定数 $1/a$ が微分の外に出せることを使った).$y,z$ 成分も同様なので,Laplacianは
\begin{equation} \nabla^2 = \frac{1}{a^2}\tilde{\nabla}^2,\qquad \tilde{\nabla}^2 \equiv \frac{\partial^2}{\partial\tilde{x}^2}+\frac{\partial^2}{\partial\tilde{y}^2}+\frac{\partial^2}{\partial\tilde{z}^2} \label{eq:lap-scale} \end{equation}と変換される.また距離は $r = \abs{a\tilde{\rr}} = a\tilde{r}$($a \gt 0$)である.これらを式 \eqref{eq:hatom-si} に代入すると
\begin{equation} \left[-\frac{\hbar^2}{2m_e a^2}\tilde{\nabla}^2 - \frac{e^2}{4\pi\varepsilon_0\,a}\,\frac{1}{\tilde{r}}\right]\psi = E\,\psi \label{eq:hatom-scaled1} \end{equation}となる.運動エネルギー項の係数 $\hbar^2/(m_e a^2)$ とポテンシャル項の係数 $e^2/(4\pi\varepsilon_0 a)$ は,どちらもエネルギーの次元を持つ(全体がエネルギー $E$ と等置されているのだから当然である).
ステップ2(エネルギーの無次元化).方程式全体を,運動エネルギー項の係数から作ったエネルギースケール
$$ \varepsilon_a \equiv \frac{\hbar^2}{m_e a^2} $$で割る.第1項の係数は $\dfrac{\hbar^2}{2m_e a^2}\Big/\dfrac{\hbar^2}{m_e a^2} = \dfrac{1}{2}$ となる.第2項の係数は
\begin{equation} \frac{e^2}{4\pi\varepsilon_0 a}\Big/\frac{\hbar^2}{m_e a^2} = \frac{e^2}{4\pi\varepsilon_0 a}\cdot\frac{m_e a^2}{\hbar^2} = \frac{m_e e^2\, a}{4\pi\varepsilon_0 \hbar^2} \label{eq:coeff-ratio} \end{equation}となる(割り算を逆数の掛け算に直し,$a^2/a = a$ とまとめた).これは無次元の数である.よって方程式は
\begin{equation} \left[-\frac{1}{2}\tilde{\nabla}^2 - \left(\frac{m_e e^2\, a}{4\pi\varepsilon_0\hbar^2}\right)\frac{1}{\tilde{r}}\right]\psi = \frac{E}{\varepsilon_a}\,\psi \label{eq:hatom-scaled2} \end{equation}という無次元の形になった.ここまで $a$ は任意である.
ステップ3(スケールの決定=Bohr半径).無次元係数がちょうど1になるように $a$ を選ぶ.すなわち
\begin{equation} \frac{m_e e^2\, a}{4\pi\varepsilon_0\hbar^2} = 1 \qquad\Longrightarrow\qquad a = \frac{4\pi\varepsilon_0\hbar^2}{m_e e^2} \equiv a_0 \label{eq:a0-def} \end{equation}この長さ $a_0$ がBohr半径である.数値を入れると $a_0 = 5.291772\times10^{-11}\,\mathrm{m} = 0.529177\,$Å となる(数値計算は下の例を参照).無次元化の要請だけから,方程式に内在する長さスケールが一意に決まったことに注意してほしい.前期量子論でBohrが軌道条件から導いた半径と同じものが,ここでは「Coulomb項と運動エネルギー項が同じ土俵に乗る長さ」として現れている.
ステップ4(エネルギースケール=Hartreeエネルギー).$a=a_0$ のときのエネルギースケール $\varepsilon_{a_0}$ を計算する.式 \eqref{eq:a0-def} の逆数 $\dfrac{1}{a_0} = \dfrac{m_e e^2}{4\pi\varepsilon_0\hbar^2}$ を代入して,
\begin{align} \varepsilon_{a_0} = \frac{\hbar^2}{m_e a_0^2} &= \frac{\hbar^2}{m_e}\left(\frac{1}{a_0}\right)^2 \notag\\ &= \frac{\hbar^2}{m_e}\left(\frac{m_e e^2}{4\pi\varepsilon_0\hbar^2}\right)^2 \notag\\ &= \frac{\hbar^2}{m_e}\cdot\frac{m_e^2 e^4}{(4\pi\varepsilon_0)^2\hbar^4} \notag\\ &= \frac{m_e e^4}{(4\pi\varepsilon_0)^2\hbar^2} \equiv E_h \label{eq:eh-def} \end{align}1行目から2行目では $1/a_0$ に式 \eqref{eq:a0-def} を代入し,2行目から3行目では括弧の2乗を展開し,3行目から4行目では $m_e^2/m_e = m_e$,$\hbar^4/\hbar^2$ を約分して $\hbar^2$ を分母に残した.この $E_h$ がHartreeエネルギーであり,数値は $E_h = 4.359745\times10^{-18}\,\mathrm{J} = 27.2114\,\mathrm{eV}$ である.
ステップ5(無次元化された方程式).$a=a_0$,$\tilde{E}\equiv E/E_h$ と書くと,式 \eqref{eq:hatom-scaled2} は
\begin{equation} \left[-\frac{1}{2}\tilde{\nabla}^2 - \frac{1}{\tilde{r}}\right]\psi = \tilde{E}\,\psi \label{eq:hatom-au} \end{equation}となる.物理定数が一切現れない,純粋に数学的な固有値問題である.これはまさに「$\hbar=m_e=e=4\pi\varepsilon_0=1$ と置いた」水素原子の方程式に他ならない.つまり原子単位系とは,長さを $a_0$ で,エネルギーを $E_h$ で測ることにした無次元化の略記法である.
以上の導出で分かったことをまとめる.(1) 「定数を1と置く」ことは,長さを $a_0$,エネルギーを $E_h$ を単位として測ることと厳密に同値である.(2) $a_0$ と $E_h$ は人為的な約束ではなく,Coulomb場中の電子という問題に内在する固有スケールである.実際,水素原子の基底状態の波動関数の広がりは $a_0$ 程度,束縛エネルギーは $E_h$ 程度になる.式 \eqref{eq:hatom-au} の基底状態固有値は $\tilde{E}_0 = -1/2$(第2章の演習で扱う)なので,水素原子の基底状態エネルギーは
$$ E_0 = -\frac{1}{2}\,E_h = -13.606\ \mathrm{eV} $$であり,よく知られたRydbergエネルギー $13.6\,$eV(水素のイオン化エネルギー)が再現される.$1\,\mathrm{Ry} = E_h/2$ という関係も,ここから読み取れる(文献にはRydberg原子単位を使うものもあるので注意).
例:$a_0$ と $E_h$ の数値をCODATA定数から計算する
基本定数の値(CODATA推奨値,有効数字7桁)は次の通りである.
$$ \hbar = 1.054572\times10^{-34}\ \mathrm{J\,s},\quad m_e = 9.109384\times10^{-31}\ \mathrm{kg},\quad e = 1.602177\times10^{-19}\ \mathrm{C},\quad k_e \equiv \frac{1}{4\pi\varepsilon_0} = 8.987552\times10^{9}\ \mathrm{kg\,m^3\,s^{-2}\,C^{-2}} $$Bohr半径.式 \eqref{eq:a0-def} を $k_e$ で書き直すと $a_0 = \hbar^2/(k_e m_e e^2)$.分子と分母を別々に計算する.
$$ \hbar^2 = (1.054572\times10^{-34})^2 = 1.112122\times10^{-68}\ \mathrm{J^2\,s^2} $$ $$ k_e m_e e^2 = (8.987552\times10^{9})\times(9.109384\times10^{-31})\times(1.602177\times10^{-19})^2 = 2.101607\times10^{-58}\ \mathrm{kg^2\,m^3\,s^{-2}} $$次元を確認すると $\dfrac{\mathrm{J^2\,s^2}}{\mathrm{kg^2\,m^3\,s^{-2}}} = \dfrac{\mathrm{kg^2\,m^4\,s^{-4}}\cdot\mathrm{s^2}}{\mathrm{kg^2\,m^3\,s^{-2}}} = \mathrm{m}$ となり,確かに長さである.数値は
$$ a_0 = \frac{1.112122\times10^{-68}}{2.101607\times10^{-58}} = 5.291772\times10^{-11}\ \mathrm{m} = 0.529177\ \text{Å} $$Hartreeエネルギー.$E_h = \hbar^2/(m_e a_0^2)$ に代入する.$a_0^2 = 2.800285\times10^{-21}\,\mathrm{m^2}$,$m_e a_0^2 = 9.109384\times10^{-31}\times2.800285\times10^{-21} = 2.550889\times10^{-51}\,\mathrm{kg\,m^2}$ なので,
$$ E_h = \frac{1.112122\times10^{-68}}{2.550889\times10^{-51}} = 4.359745\times10^{-18}\ \mathrm{J} $$これを素電荷で割れば電子ボルト単位になる:
$$ E_h = \frac{4.359745\times10^{-18}}{1.602177\times10^{-19}} = 27.2114\ \mathrm{eV} $$さらにAvogadro数 $N_A = 6.022141\times10^{23}\,\mathrm{mol^{-1}}$ を掛ければモル単位になる:
$$ E_h N_A = 4.359745\times10^{-18}\times 6.022141\times10^{23} = 2625.50\ \mathrm{kJ/mol} = 627.5\ \mathrm{kcal/mol} $$すなわち $1\ \mathrm{Hartree} = 27.2114\ \mathrm{eV} = 627.5\ \mathrm{kcal/mol}$ である.この3通りの表現は,物理(eV)と化学(kcal/mol)の文献を行き来する際に頻繁に必要になる.
数学ノート:単位の無次元化の一般論
上で行った操作は,任意の物理方程式に適用できる一般的な手続きである.手順を整理しておく.
- 方程式に現れる次元定数を列挙する.水素原子の場合は $\hbar$,$m_e$,$e^2/4\pi\varepsilon_0$ の3つ(最後の組み合わせは常にこの形でしか現れないので1つと数える).
- 関与する基本次元を数える.ここでは質量 M,長さ L,時間 T の3つ(電荷は $e^2/4\pi\varepsilon_0$ が $\mathrm{ML^3T^{-2}}$ の次元を持つため,独立には現れない).
- 各変数を特性スケールで割って無次元変数を作る.$\rr = a\tilde{\rr}$,$E = \varepsilon\tilde{E}$ のように,スケール($a$,$\varepsilon$ など)は未定のまま代入する.
- 1つの項の係数で方程式全体を割り,残った無次元係数を1と置いてスケールを決める.独立に選べるスケールの数は基本次元の数に等しい.
定数の数(3)と基本次元の数(3)が一致する水素原子では,すべての係数を1にでき,方程式には無次元パラメータが残らない.これは深い意味を持つ:非相対論的なCoulomb系には調整可能な無次元定数が存在しない.どの原子・分子・固体も(原子番号という整数を除けば)同じ一つの方程式に従うのであり,これが§1.2で述べたハミルトニアンの普遍性の次元解析的な表現である.一方,定数を1つ追加すると事情が変わる.光速 $c$ を加えると,無次元の組み合わせ
$$ \alpha = \frac{e^2}{4\pi\varepsilon_0\hbar c} \simeq \frac{1}{137.036} $$(微細構造定数)が作れてしまい,これはどんな単位系を選んでも消せない.$\alpha$ は「相対論効果の強さ」を表す真のパラメータであり,原子単位系でも $c = 1/\alpha \simeq 137.036$ という値を持ち続ける.相対論を無視できるかどうかの判定にこの数が現れることを,すぐ後(§1.3.3)で見る.
1.3.2 多電子ハミルトニアンの原子単位表示
水素原子で確立した処方を,多電子ハミルトニアン \eqref{eq:h-si} 全体に適用しよう.これから行うことは,すべての座標(電子も核も)を $a_0$ でスケールし,方程式を $E_h$ で割るという機械的な操作である.項の種類は運動エネルギー項とCoulomb項の2種類しかないので,それぞれ1つずつ変換を確認すれば十分である.
導出:多電子ハミルトニアンの無次元化
全座標を $\rr_i = a_0\tilde{\rr}_i$,$\bm{R}_I = a_0\tilde{\bm{R}}_I$ とスケールする.
運動エネルギー項.式 \eqref{eq:lap-scale} と同様に $\nabla_i^2 = a_0^{-2}\tilde{\nabla}_i^2$ なので,
\begin{align} -\frac{\hbar^2}{2m_e}\nabla_i^2 &= -\frac{\hbar^2}{2m_e a_0^2}\tilde{\nabla}_i^2 \notag\\ &= -\frac{1}{2}\left(\frac{\hbar^2}{m_e a_0^2}\right)\tilde{\nabla}_i^2 = -\frac{E_h}{2}\tilde{\nabla}_i^2 \label{eq:kin-scale} \end{align}となる.最後の等号で,式 \eqref{eq:eh-def} の途中で示した関係 $E_h = \hbar^2/(m_e a_0^2)$ を使った.核の運動エネルギー項も同様で,$-\dfrac{\hbar^2}{2M_I}\nabla_I^2 = -\dfrac{m_e}{M_I}\cdot\dfrac{\hbar^2}{2m_e a_0^2}\tilde{\nabla}_I^2 = -\dfrac{E_h}{2(M_I/m_e)}\tilde{\nabla}_I^2$ となる(分母分子に $m_e$ を掛けて $E_h$ をくくり出した).すなわち原子単位系では核の質量は「電子質量の何倍か」という無次元数 $M_I/m_e$ で表される.
Coulomb項.代表として電子間反発の1つの項を変換する.$\abs{\rr_i-\rr_j} = a_0\abs{\tilde{\rr}_i-\tilde{\rr}_j}$ なので,
\begin{equation} \frac{e^2}{4\pi\varepsilon_0\abs{\rr_i-\rr_j}} = \frac{e^2}{4\pi\varepsilon_0 a_0}\cdot\frac{1}{\abs{\tilde{\rr}_i-\tilde{\rr}_j}} \label{eq:coul-scale1} \end{equation}となる.ここで係数 $e^2/(4\pi\varepsilon_0 a_0)$ を計算する.$1/a_0 = m_e e^2/(4\pi\varepsilon_0\hbar^2)$(式 \eqref{eq:a0-def})を代入して,
\begin{equation} \frac{e^2}{4\pi\varepsilon_0 a_0} = \frac{e^2}{4\pi\varepsilon_0}\cdot\frac{m_e e^2}{4\pi\varepsilon_0\hbar^2} = \frac{m_e e^4}{(4\pi\varepsilon_0)^2\hbar^2} = E_h \label{eq:coul-scale2} \end{equation}となり(最後の等号は式 \eqref{eq:eh-def} の定義そのもの),Coulomb項の自然な大きさも $E_h$ であることが分かる.電子–核引力項,核間反発項も分母の距離が $a_0$ でスケールされるだけなので全く同様である.
まとめ.ハミルトニアン全体が $\hat{H} = E_h\times(\text{無次元の演算子})$ の形になったので,Schrödinger方程式 $\hat{H}\psi = E\psi$ の両辺を $E_h$ で割れば,無次元の方程式 $\tilde{H}\psi = \tilde{E}\psi$($\tilde{E}=E/E_h$)が得られる.
以後はチルダを省き,「長さは $a_0$ 単位,エネルギーは $E_h$ 単位」と約束して書く.得られた多電子ハミルトニアンの原子単位表示は次の通りである.これは本書全体の基礎方程式なので,強調枠に入れておく.
ここで $M_I$ は電子質量を単位とした核の質量(無次元;たとえば陽子なら $M_I = 1836.15$)である.式 \eqref{eq:h-si} と見比べてほしい.物理定数がすべて消え,構造——運動エネルギーと3種類のCoulomb相互作用——だけが残った.以後の章の式は基本的にこの形で書かれる.
1.3.3 導出単位:時間・速度,そして換算表
質量・電荷・作用・Coulomb定数を固定すると,他のすべての物理量の単位が自動的に決まる.主要なものを導出しておこう.
時間の原子単位.定常状態の位相因子 $e^{-iEt/\hbar}$ の指数は無次元でなければならないから,時間の自然な単位は $t_0 = \hbar/E_h$ である.数値は
\begin{equation} t_0 = \frac{\hbar}{E_h} = \frac{1.054572\times10^{-34}\ \mathrm{J\,s}}{4.359745\times10^{-18}\ \mathrm{J}} = 2.418884\times10^{-17}\ \mathrm{s} \simeq 24.2\ \text{アト秒} \label{eq:time-au} \end{equation}である.これは電子が原子内を「一周」する程度の時間スケールであり,電子ダイナミクスがアト秒($10^{-18}\,$s)の世界であることを示している.逆に言えば $1\,\mathrm{fs} = 41.34\,$a.u. であり,第一原理分子動力学の時間刻み(典型的には $0.5$–$1\,$fs)が原子単位でどの程度かの感覚はここから得られる.
速度の原子単位.長さを時間で割って $v_0 = a_0/t_0 = a_0 E_h/\hbar$ である.この積を整理すると見通しの良い形になる.$a_0 = 4\pi\varepsilon_0\hbar^2/(m_e e^2)$ と $E_h = m_e e^4/[(4\pi\varepsilon_0)^2\hbar^2]$ を掛けると,
\begin{align} a_0 E_h &= \frac{4\pi\varepsilon_0\hbar^2}{m_e e^2}\cdot\frac{m_e e^4}{(4\pi\varepsilon_0)^2\hbar^2} = \frac{e^2}{4\pi\varepsilon_0} \label{eq:a0Eh} \end{align}($\hbar^2$,$m_e$,$4\pi\varepsilon_0$ の1乗,$e^2$ が約分された).よって
\begin{equation} v_0 = \frac{a_0 E_h}{\hbar} = \frac{e^2}{4\pi\varepsilon_0\hbar} = \left(\frac{e^2}{4\pi\varepsilon_0\hbar c}\right) c = \alpha c = \frac{c}{137.036} = 2.187691\times10^{6}\ \mathrm{m/s} \label{eq:velocity-au} \end{equation}となる.速度の原子単位は光速の $1/137$,すなわち微細構造定数 $\alpha$ がここで顔を出す.水素原子の基底状態の電子の典型的速さがちょうど $1\,$a.u. $=\alpha c$ であることから,軽い原子の電子は十分非相対論的であり,式 \eqref{eq:h-au} の非相対論的取り扱いが正当化される.一方,原子番号 $Z$ の原子の最内殻($1s$)電子の速さは約 $Z\,$a.u. になるため,$Z$ が大きい元素(たとえば金 $Z=79$ では $v/c \simeq 79/137 \simeq 0.58$)では相対論効果が無視できない.重元素を扱う際にスカラー相対論補正やスピン–軌道相互作用が必要になる理由は,この簡単な見積りに尽きている.
主要な物理量について,原子単位とSI単位の対応を表1.1にまとめる.
| 物理量 | 原子単位での値 | 定義(SI定数による表式) | SI単位系での値 |
|---|---|---|---|
| 質量 | 1 | $m_e$(電子質量) | $9.109384\times10^{-31}\ \mathrm{kg}$ |
| 電荷 | 1 | $e$(素電荷) | $1.602177\times10^{-19}\ \mathrm{C}$ |
| 作用(角運動量) | 1 | $\hbar = h/2\pi$ | $1.054572\times10^{-34}\ \mathrm{J\,s}$ |
| Coulomb定数 | 1 | $k_e = 1/4\pi\varepsilon_0$ | $8.987552\times10^{9}\ \mathrm{kg\,m^3\,s^{-2}\,C^{-2}}$ |
| 長さ | 1 | $a_0 = 4\pi\varepsilon_0\hbar^2/(m_e e^2)$(Bohr半径) | $5.291772\times10^{-11}\ \mathrm{m} = 0.529177$ Å |
| エネルギー | 1 | $E_h = m_e e^4/(4\pi\varepsilon_0\hbar)^2$(Hartree) | $4.359745\times10^{-18}\ \mathrm{J} = 27.2114\ \mathrm{eV}$ |
| 時間 | 1 | $\hbar/E_h$ | $2.418884\times10^{-17}\ \mathrm{s}$ |
| 速度 | 1 | $v_0 = a_0E_h/\hbar = \alpha c$ | $2.187691\times10^{6}\ \mathrm{m\,s^{-1}}$ |
| 磁束密度 | 1 | $\hbar/(e a_0^2)$ | $2.350518\times10^{5}\ \mathrm{T}$ |
| 磁気双極子モーメント | 1 | $e\hbar/m_e = 2\mu_B$ | $1.854802\times10^{-23}\ \mathrm{J\,T^{-1}}$ |
磁束密度の単位 $2.35\times10^{5}\,$T(23.5万テスラ)は,実験室で作れる定常磁場(数十T)がいかに「弱い」摂動かを教えてくれる.また磁気モーメントの原子単位がBohr磁子 $\mu_B$ の2倍であることは,スピン磁気モーメントを扱う章で符号・係数を確認する際に役立つ.
例:エネルギースケールの換算に慣れる
実際の研究で頻出する換算をいくつか挙げる.$1\,E_h = 27.2114\,\mathrm{eV} = 627.5\,\mathrm{kcal/mol} = 2625.5\,\mathrm{kJ/mol}$ を基準にする.
- 化学結合:H$_2$分子の解離エネルギー(実験値)は $4.75\,\mathrm{eV} = 4.75/27.2114 = 0.1746\,E_h$.典型的な共有結合は $0.1$–$0.4\,E_h$ 程度である.
- 化学精度:反応エネルギーの目標精度としてよく引かれる「化学精度」$1\,\mathrm{kcal/mol}$ は $1/627.5 = 1.59\times10^{-3}\,E_h = 0.0434\,\mathrm{eV}$.全エネルギーが数十〜数千 $E_h$ の系でこの精度を出す必要があることが,第一原理計算の数値的な厳しさの源である.
- 室温:$k_BT$($T=300\,$K)は $0.02585\,\mathrm{eV} = 9.50\times10^{-4}\,E_h \simeq 1\,\mathrm{kcal/mol}$.熱揺らぎと化学精度が同程度であることは偶然ではなく,化学反応の速度論が室温の熱励起で決まっているためである.
- DFT計算の収束判定:自己無撞着計算のエネルギー収束条件は $10^{-6}$–$10^{-8}\,E_h$ 程度に設定されることが多い.$10^{-6}\,E_h \simeq 3\times10^{-5}\,\mathrm{eV}$ である.
物理的意味:原子単位系は「電子の世界の等身大の物差し」
原子単位系の効用は,式が短くなることだけではない.この単位系では,原子スケールの現象に関わる量はおおむね1程度の数になる.結合長は数 a.u.,価電子の束縛エネルギーは $0.1$–$1$ a.u.,価電子の速さは $\sim 1$ a.u..計算結果が $10^{6}$ や $10^{-9}$ といった値を示したら,それは何かが桁違いに大きい・小さいという物理的なシグナルであり,単位の換算ミスや数値誤差の検出にも役立つ.逆に,無次元化して消せなかった数——核の質量比 $M_I/m_e \sim 10^{3\text{–}5}$ と光速 $c = 137$——は,この世界の「本当の」大小関係を表しており,それぞれBorn–Oppenheimer近似(次節)と非相対論近似の妥当性を支配している.単位系を整えるという一見事務的な作業が,実は近似の構造を浮かび上がらせているのである.
1.4 Born–Oppenheimer(断熱)近似
多電子問題を解く最初の,そして最も基本的な近似がBorn–Oppenheimer近似(断熱近似)である[3].骨子は単純で,「核は電子よりはるかに重いので,電子から見れば核はほぼ止まっている.だからまず核を固定して電子だけの問題を解き,核の運動はその後で考えればよい」というものである.しかしこの直観を式にする過程には学ぶべき構造が多く,またどの項を落としたのかを明示しておかないと,この近似が破れる状況(光化学反応,電子–格子相互作用など)を正しく認識できない.本節では省略なしに導出する.単位系は前節で導入した原子単位を用いる.
1.4.1 質量比の見積り
出発点は質量比である.最も軽い核である陽子でも
\begin{equation} \frac{m_e}{M_p} = \frac{9.109384\times10^{-31}\ \mathrm{kg}}{1.672622\times10^{-27}\ \mathrm{kg}} = \frac{1}{1836.15} \simeq 5.45\times10^{-4} \label{eq:mass-ratio} \end{equation}であり,たとえば炭素原子核($^{12}$C)なら $M/m_e = 12\times1822.89 \simeq 2.19\times10^4$,遷移金属では $10^5$ のオーダーに達する($1\,\mathrm{u} = 1822.89\,m_e$).式 \eqref{eq:h-au} を見ると,核の運動エネルギー項だけが $1/M_I \lesssim 10^{-3}$ という小さな係数を持ち,他の項はすべて $O(1)$ である.この非対称性を利用する.
1.4.2 電子ハミルトニアンと断熱基底
これから行うことは次の3段階である.(1) ハミルトニアンを「核の運動エネルギー」と「それ以外」に分割し,後者(電子ハミルトニアン)の固有関数系を核配置ごとに用意する.(2) 全波動関数をその固有関数系で展開して全Schrödinger方程式に代入し,核の波動関数に対する厳密な連立方程式を導く.(3) 連立を引き起こしている項(非断熱結合項)が $1/M_I$ に比例することを確認して落とす.
まず分割する.電子座標をまとめて $\rr \equiv (\rr_1,\dots,\rr_{N_e})$,核座標をまとめて $\bm{R} \equiv (\bm{R}_1,\dots,\bm{R}_{N_n})$ と略記し,
\begin{equation} \hat{H} = \hat{T}_n + \hat{H}_e(\bm{R}), \qquad \hat{T}_n = -\sum_{I=1}^{N_n}\frac{1}{2M_I}\nabla_I^2 \label{eq:h-split} \end{equation} \begin{equation} \hat{H}_e(\bm{R}) = -\frac{1}{2}\sum_{i=1}^{N_e}\nabla_i^2 - \sum_{i=1}^{N_e}\sum_{I=1}^{N_n}\frac{Z_I}{\abs{\rr_i-\bm{R}_I}} + \sum_{i \lt j}^{N_e}\frac{1}{\abs{\rr_i-\rr_j}} + \sum_{I \lt J}^{N_n}\frac{Z_I Z_J}{\abs{\bm{R}_I-\bm{R}_J}} \label{eq:h-elec} \end{equation}と書く.$\hat{H}_e(\bm{R})$ を電子ハミルトニアンと呼ぶ.ここで重要なのは,$\hat{H}_e$ の中では $\bm{R}$ が微分されない,すなわち $\bm{R}$ は演算子ではなくパラメータとして入っているという点である(核間反発 $\hat{V}_{nn}$ は電子座標を含まない定数項だが,慣例に従い $\hat{H}_e$ に含めておく.こうすると後で出てくるポテンシャルエネルギー面が「全エネルギー」の面になり便利である).
核配置 $\bm{R}$ を1つ固定するごとに,電子ハミルトニアンの固有値問題
を考えることができる.これを電子Schrödinger方程式と呼び,その固有関数 $\Phi_k(\rr;\bm{R})$ を断熱電子状態,固有値 $E_k(\bm{R})$ を第 $k$ 断熱ポテンシャルエネルギー面(potential energy surface; PES)と呼ぶ.セミコロンの後の $\bm{R}$ は「パラメータとしての依存性」を表す記法である.各 $\bm{R}$ において $\{\Phi_k\}$ は電子座標の関数として完全正規直交系をなす:
\begin{equation} \braket{\Phi_l|\Phi_k}_{\rr} \equiv \int \Phi_l^*(\rr;\bm{R})\,\Phi_k(\rr;\bm{R})\,\dd\rr = \delta_{lk} \quad(\text{任意の }\bm{R}\text{ で}) \label{eq:phi-ortho} \end{equation}ここで $\int\dd\rr$ は全電子座標についての積分(スピンを含めるときはスピン和も込み;§1.5),添字 $\rr$ は「電子座標についてのみ積分し,$\bm{R}$ は残る」ことを明示している.
1.4.3 Born–Huang展開と厳密な連立方程式
全波動関数 $\Psi(\rr,\bm{R})$ を,各核配置での断熱電子状態で展開する:
\begin{equation} \Psi(\rr,\bm{R}) = \sum_k \chi_k(\bm{R})\,\Phi_k(\rr;\bm{R}) \label{eq:born-huang} \end{equation}これをBorn–Huang展開と呼ぶ[3,4].展開係数 $\chi_k(\bm{R})$ は核座標の関数であり,「電子状態 $k$ の上での核の波動関数」と解釈される.ここまでは近似ではない:固定した $\bm{R}$ ごとに $\{\Phi_k(\cdot\,;\bm{R})\}$ は電子座標の関数空間の完全系だから,任意の $\Psi$ はこの形に必ず書ける(その代償として係数 $\chi_k$ が $\bm{R}$ に依存する).
導出:核波動関数に対する厳密な連立方程式
全Schrödinger方程式 $(\hat{T}_n + \hat{H}_e)\Psi = E\Psi$ に展開 \eqref{eq:born-huang} を代入する.まず $\hat{H}_e$ の作用は簡単である.$\hat{H}_e$ は電子座標にしか作用しない($\chi_k(\bm{R})$ は素通りする)ので,固有値方程式 \eqref{eq:elec-se} を使って
\begin{equation} \hat{H}_e\sum_k \chi_k\Phi_k = \sum_k \chi_k\,\hat{H}_e\Phi_k = \sum_k \chi_k\,E_k(\bm{R})\,\Phi_k \label{eq:he-action} \end{equation}となる.次に $\hat{T}_n$ の作用を計算する.$\Phi_k(\rr;\bm{R})$ は $\bm{R}$ にパラメータ依存するので,$\nabla_I$ は $\chi_k$ と $\Phi_k$ の両方に作用することに注意する.積の微分則を2回使う.まず1階微分は
\begin{equation} \nabla_I(\chi_k\Phi_k) = (\nabla_I\chi_k)\Phi_k + \chi_k(\nabla_I\Phi_k) \label{eq:prod1} \end{equation}である.この両辺にもう一度 $\nabla_I\cdot$ を作用させ,右辺の各項に再び積の微分則を使うと
\begin{align} \nabla_I^2(\chi_k\Phi_k) &= \nabla_I\cdot\big[(\nabla_I\chi_k)\Phi_k\big] + \nabla_I\cdot\big[\chi_k(\nabla_I\Phi_k)\big] \notag\\ &= (\nabla_I^2\chi_k)\Phi_k + (\nabla_I\chi_k)\cdot(\nabla_I\Phi_k) + (\nabla_I\chi_k)\cdot(\nabla_I\Phi_k) + \chi_k(\nabla_I^2\Phi_k) \notag\\ &= (\nabla_I^2\chi_k)\Phi_k + 2(\nabla_I\chi_k)\cdot(\nabla_I\Phi_k) + \chi_k(\nabla_I^2\Phi_k) \label{eq:prod2} \end{align}となる(2行目で同じ交差項が2回現れるのでまとめて係数2にした).したがって
\begin{equation} \hat{T}_n(\chi_k\Phi_k) = -\sum_I\frac{1}{2M_I}\Big[(\nabla_I^2\chi_k)\Phi_k + 2(\nabla_I\chi_k)\cdot(\nabla_I\Phi_k) + \chi_k(\nabla_I^2\Phi_k)\Big] \label{eq:tn-action} \end{equation}である.第1項は「核の波動関数だけが微分される」項,第2・第3項は「電子波動関数の核座標依存性が微分される」項である.後者こそが,電子状態間の遷移(非断熱遷移)を引き起こす.
式 \eqref{eq:he-action} と \eqref{eq:tn-action} を全Schrödinger方程式に入れ,左から $\Phi_l^*(\rr;\bm{R})$ を掛けて電子座標で積分する(状態 $l$ への射影).直交性 \eqref{eq:phi-ortho} により,$\Phi_k$ がそのまま残っている項では $k=l$ だけが生き残る:
\begin{align} &\braket{\Phi_l|\textstyle\sum_k \chi_k E_k\Phi_k}_{\rr} = E_l(\bm{R})\,\chi_l(\bm{R}), \qquad \braket{\Phi_l|\textstyle\sum_k (\nabla_I^2\chi_k)\Phi_k}_{\rr} = \nabla_I^2\chi_l(\bm{R}) \notag \end{align}一方,$\nabla_I\Phi_k$,$\nabla_I^2\Phi_k$ を含む項では直交性が使えず,行列要素
\begin{equation} \bm{F}_{lk}^{I}(\bm{R}) \equiv \braket{\Phi_l|\nabla_I\Phi_k}_{\rr}, \qquad G_{lk}^{I}(\bm{R}) \equiv \braket{\Phi_l|\nabla_I^2\Phi_k}_{\rr} \label{eq:coupling-def} \end{equation}が残る.$\bm{F}_{lk}^I$ を1階非断熱結合ベクトル(derivative coupling),$G_{lk}^I$ を2階非断熱結合と呼ぶ.以上をまとめると,核の波動関数 $\chi_l$ に対する方程式は
\begin{equation} \Big[-\sum_I\frac{1}{2M_I}\nabla_I^2 + E_l(\bm{R})\Big]\chi_l(\bm{R}) \;-\;\sum_k\sum_I\frac{1}{2M_I}\Big[2\bm{F}_{lk}^{I}\cdot\nabla_I + G_{lk}^{I}\Big]\chi_k(\bm{R}) = E\,\chi_l(\bm{R}) \label{eq:coupled} \end{equation}となる.これはまだ厳密である.すべての電子状態 $k$ の核波動関数 $\chi_k$ が第2の和を通じて互いに結合しており,多電子問題の完全な複雑さがここに保存されている.
得られた式 \eqref{eq:coupled} の構造を読み取ろう.もし第2の和(非断熱結合項)がなければ,各電子状態 $l$ ごとに方程式が独立になり,「核は $E_l(\bm{R})$ というポテンシャルの中を運動する」という描像が厳密に成立する.結合項はすべて $1/2M_I$ という小さな係数を伴っている.そこで次の近似を置く.
定義:断熱近似とBorn–Oppenheimer近似
式 \eqref{eq:coupled} において,
- 断熱近似:異なる電子状態を結ぶ非対角項($k\neq l$ の $\bm{F}_{lk}^I$,$G_{lk}^I$)をすべて無視する.系は1枚のポテンシャルエネルギー面 $E_l(\bm{R})$ 上に留まり続ける.
- Born–Oppenheimer近似:さらに対角項($k=l$)も無視する.核の方程式は \begin{equation} \Big[-\sum_I\frac{1}{2M_I}\nabla_I^2 + E_l(\bm{R})\Big]\chi_l(\bm{R}) = E\,\chi_l(\bm{R}) \label{eq:bo-final} \end{equation} となる.
本書で単に「BO近似」というときは後者を指す.対角項だけを残す中間の近似(対角Born–Oppenheimer補正,DBOC)は高精度分子計算で用いられる.
1.4.4 落とした項はどれほど小さいか
非断熱結合項を落とす操作の妥当性を,2つの角度から確認する.
(i) 質量比による抑制.結合項は $\dfrac{1}{2M_I}\big[2\bm{F}^I_{lk}\cdot\nabla_I + G^I_{lk}\big]$ の形で,あらわに $1/M_I \sim 10^{-3}$–$10^{-5}$ の因子を持つ.一方,$\bm{F}^I_{lk}$ や $G^I_{lk}$ は電子波動関数が核座標に対してどれだけ敏感かを表す量で,電子構造が結合長程度($\sim$ 数 a.u.)のスケールで変化することから,典型的には $O(1)$(a.u.)である.したがって結合項のエネルギー寄与は電子エネルギー間隔($O(0.1$–$1)\,E_h$)に比べて $10^{-3}$ 以下に抑えられる.
(ii) エネルギーギャップによる抑制.より精密には,非対角結合 $\bm{F}^I_{lk}$ の大きさは電子状態間のエネルギー差に反比例することが示せる.これは重要な関係式なので導出しておく.
導出:非断熱結合とエネルギーギャップの関係(Hellmann–Feynman型公式)
電子固有値方程式 \eqref{eq:elec-se} の両辺を核座標 $\bm{R}_I$ で微分する(核座標はパラメータなので,これは通常の偏微分である).積の微分則を左辺・右辺それぞれに使うと
\begin{equation} (\nabla_I\hat{H}_e)\Phi_k + \hat{H}_e(\nabla_I\Phi_k) = (\nabla_I E_k)\Phi_k + E_k(\nabla_I\Phi_k) \label{eq:hf-step1} \end{equation}となる.ここで $\nabla_I\hat{H}_e$ は,$\hat{H}_e$ の中で $\bm{R}_I$ にあらわに依存する項(電子–核引力と核間反発)だけを微分した演算子である.左から $\Phi_l^*$($l\neq k$)を掛けて電子座標で積分すると,各項は
\begin{align} \braket{\Phi_l|(\nabla_I\hat{H}_e)|\Phi_k} + \braket{\Phi_l|\hat{H}_e|\nabla_I\Phi_k} &= (\nabla_I E_k)\underbrace{\braket{\Phi_l|\Phi_k}}_{=0} + E_k\braket{\Phi_l|\nabla_I\Phi_k} \label{eq:hf-step2} \end{align}となる(右辺第1項は直交性 \eqref{eq:phi-ortho} で消えた).左辺第2項では $\hat{H}_e$ のエルミート性を使い,$\hat{H}_e$ を左の $\Phi_l$ に作用させる:$\braket{\Phi_l|\hat{H}_e|\nabla_I\Phi_k} = \braket{\hat{H}_e\Phi_l|\nabla_I\Phi_k} = E_l\braket{\Phi_l|\nabla_I\Phi_k}$($E_l$ は実数).これを移項して整理すると
\begin{equation} \braket{\Phi_l|(\nabla_I\hat{H}_e)|\Phi_k} = (E_k - E_l)\braket{\Phi_l|\nabla_I\Phi_k} \;\;\Longrightarrow\;\; \bm{F}^I_{lk} = \frac{\braket{\Phi_l|(\nabla_I\hat{H}_e)|\Phi_k}}{E_k(\bm{R}) - E_l(\bm{R})} \quad(l\neq k) \label{eq:f-offdiag} \end{equation}が得られる.分子はポテンシャルの核座標微分(力の演算子)の行列要素で $O(1)$ の量,分母は電子状態間のエネルギーギャップである.
式 \eqref{eq:f-offdiag} から,BO近似の適用限界が定量的に読み取れる.電子状態間のギャップが大きい(基底状態が励起状態からよく離れている)ほど結合は小さく,BO近似は良い.逆に2つのPESが接近・交差する領域($E_k \simeq E_l$,いわゆる円錐交差や擬交差)では $\bm{F}^I_{lk}$ が発散的に増大し,BO近似は破綻する.光励起後の分子の無輻射緩和,視物質の異性化などの光化学はまさにこの破綻領域の物理である.また金属では占有・非占有状態のギャップがゼロなので,電子–格子相互作用(フォノンによる電気抵抗,超伝導)を論じるときには落とした項を摂動として拾い直す必要がある.本書の範囲(基底状態の構造・結合エネルギー)ではBO近似は極めて良い近似であり,以後これを前提とする.
なお対角項についても一言.1階の対角結合 $\bm{F}^I_{ll}$ は,電子波動関数が実関数に取れる場合(時間反転対称性がある場合)には恒等的にゼロである.これは規格化条件から従う:$\braket{\Phi_l|\Phi_l} = 1$ の両辺を $\bm{R}_I$ で微分すると
\begin{equation} 0 = \nabla_I\braket{\Phi_l|\Phi_l} = \braket{\nabla_I\Phi_l|\Phi_l} + \braket{\Phi_l|\nabla_I\Phi_l} = 2\braket{\Phi_l|\nabla_I\Phi_l} = 2\bm{F}^I_{ll} \label{eq:f-diag-zero} \end{equation}となる(実関数なら $\braket{\nabla_I\Phi_l|\Phi_l} = \braket{\Phi_l|\nabla_I\Phi_l}$ が成り立つことを使った).よって断熱近似で実際に落ちる対角項は2階の $G^I_{ll}$ のみで,これは $O(1/M_I)$ のPES補正(DBOC)を与えるにすぎない(章末演習1.3).
1.4.5 ポテンシャルエネルギー面という概念
BO近似の帰結として,「原子核は,電子が作るポテンシャルエネルギー面 $E_0(\bm{R})$ の上を運動する」という,化学と物質科学の根幹をなす描像が確立した.$E_0(\bm{R})$ は $3N_n$ 変数の関数であり,その地形がその物質の構造と動力学のすべてを支配する:
- 極小点 $\nabla_{\bm{R}}E_0 = 0$(かつHessianが正定値)が安定構造・準安定構造である.最深の極小が基底構造,他の極小は多形や異性体に対応する.
- 極小のまわりの曲率(Hessian行列 $\partial^2 E_0/\partial R_{I\alpha}\partial R_{J\beta}$)が振動数・フォノンを与える.
- 極小どうしを結ぶ峠(鞍点)が化学反応や構造相転移の遷移状態であり,峠の高さが活性化エネルギーである.
- 面の勾配 $\bm{F}_I = -\nabla_{I}E_0(\bm{R})$ が原子に働く力であり,構造最適化と第一原理分子動力学の駆動力である.
第一原理計算の実務は,その大部分が「電子Schrödinger方程式 \eqref{eq:elec-se} を(DFTで近似的に)解いて $E_0(\bm{R})$ とその勾配を評価する」ことに帰着する.図1.2に,最も簡単な場合として二原子分子のPES(この場合は核間距離 $R$ の1変数関数なのでポテンシャル曲線)を模式的に示す.
例:電子・振動・回転のエネルギースケール比
BO近似の内部構造は,分子のスペクトルに現れる3つのエネルギースケールの分離として観測できる.これを質量比だけから見積もってみよう.原子単位で考える.
電子励起:電子は $a_0 \sim 1$ の領域に閉じ込められているので,運動エネルギーは $\sim\hbar^2/(m_e a_0^2) = 1\,E_h$.よって電子準位の間隔は $\Delta E_{\mathrm{el}} \sim 1\,E_h$($\sim$数 eV).
分子振動:PESの極小まわりの曲率(ばね定数)$k$ は電子構造で決まるので,原子単位で $k \sim E_h/a_0^2 \sim 1$.核(質量 $M$)の調和振動数は $\omega = \sqrt{k/M}$ だから,
$$ \hbar\omega = \hbar\sqrt{\frac{k}{M}} \sim \sqrt{\frac{\hbar^2 E_h}{M a_0^2}} = \sqrt{\frac{m_e}{M}\cdot\underbrace{\frac{\hbar^2}{m_e a_0^2}}_{=E_h}\cdot E_h} = \left(\frac{m_e}{M}\right)^{1/2} E_h $$(根号の中で分母分子に $m_e$ を掛けて $E_h$ を2つ作った).$m_e/M \sim 10^{-4}$ とすれば $\hbar\omega \sim 10^{-2}\,E_h \sim 0.1\,\mathrm{eV}$.実測の分子振動(数百〜数千 cm$^{-1}$,すなわち $0.01$–$0.4\,$eV)と整合する.
分子回転:慣性モーメント $I \sim M a_0^2$ の回転準位間隔は $\hbar^2/I \sim (m_e/M)\,E_h \sim 10^{-4}\,E_h \sim$ meV.
すなわち $\Delta E_{\mathrm{el}} : \hbar\omega : \Delta E_{\mathrm{rot}} \sim 1 : (m_e/M)^{1/2} : (m_e/M)$ という階層が,質量比のみから導かれる.BornとOppenheimerの原論文[3]は,まさに $(m_e/M)^{1/4}$ を展開パラメータとする系統的摂動論としてこの階層を導いた.
物理的意味:「断熱」の名の由来
核がゆっくり動くとき,電子は各瞬間の核配置に対する固有状態(の1つ)に「即座に」調整され,状態のラベル $k$ を保ったまま追従する——これは量子力学の断熱定理(ハミルトニアンをゆっくり変化させると系は固有状態に留まり続ける)の状況そのものである.「断熱近似」という名称はここに由来する.熱の出入りがないという熱力学の断熱とは直接の関係はないので混同しないように.なお,電子が基底状態に留まり続けるという仮定は,絶対零度近くの構造・結合の議論では良いが,光励起や高温の非平衡ダイナミクスでは破れる.それらは「非断熱動力学」として現代でも活発な研究対象である.
以上で,多電子問題は「核配置をパラメータとする電子だけの問題 \eqref{eq:elec-se}」に整理された.以後の章では,断りのない限り電子Schrödinger方程式 \eqref{eq:elec-se} を対象とし,$\hat{H}_e$ を単に $\hat{H}$,$E_0(\bm{R})$ を単に $E$ と書く.次節では,この電子波動関数 $\Phi(\rr_1,\dots,\rr_{N_e})$ が満たすべき基本的な対称性を確認する.
1.5 多電子波動関数の基本性質
前節までで,解くべき方程式は電子Schrödinger方程式 \eqref{eq:elec-se} に絞り込まれた.しかし方程式を書いただけでは問題は確定しない.「どのような関数を解の候補として許すか」——すなわち境界条件と対称性——を指定して初めて固有値問題は定まる.多電子系の場合,この指定は数学的な便宜ではなく,電子が互いに区別できないという物理的事実から導かれる非常に強い制限になる.本節ではその制限を,前提を明示しながら一歩ずつ導く.ここで得られる反対称性は第5章(Slater行列式とHartree–Fock法)と第6章(第二量子化)の全体を支配する原理であり,詳しい帰結——交換積分,Fermiホール,Pauliの排他律の定量的な表現——はそちらで展開する.本節はその土台を据えることに徹する.
記号を整理しておく.以後,電子数を単に $N$ と書く($N \equiv N_e$).また核座標 $\bm{R}$ への依存は,§1.4の意味で常にパラメータとして背後にあるものとし,いちいち書かない.多電子波動関数は $\Psi$ と書く(§1.4の断熱電子状態 $\Phi_k(\rr;\bm{R})$ に,以下で導入するスピン変数を明示し,状態の添字 $k$ とパラメータ $\bm{R}$ を省いたものである).
1.5.1 スピンと合成変数
式 \eqref{eq:h-elec} の電子ハミルトニアンを見直すと,そこにはスピンを表す演算子が一切現れていない.非相対論的な近似では,Coulomb相互作用も運動エネルギーもスピンには目もくれないからである.それにもかかわらず,電子状態を論じるうえでスピンを省くことはできない.理由は2つある.第一に,電子はスピン $1/2$ を持つ粒子であり,1個の電子の状態を完全に指定するには空間座標だけでなくスピンの向きも要る.第二に,これから導く反対称性の要請はスピンを含めた変数について課されるため,スピンを落とすと要請そのものが書けない.そこで空間座標とスピン座標をまとめた変数を導入する.
定義:合成変数 $x=(\rr,\sigma)$ と積分記号
1個の電子の状態変数を,空間座標 $\rr$ とスピン変数 $\sigma$ の組
$$ x \equiv (\rr,\sigma),\qquad \sigma \in \{\uparrow,\downarrow\} $$と定義する.$\sigma=\uparrow$ はスピン $z$ 成分 $s_z=+1/2$,$\sigma=\downarrow$ は $s_z=-1/2$ を表す(原子単位では $\hbar=1$ なので $s_z$ は $\pm1/2$ という無次元の数である).$x$ についての積分記号は,空間積分とスピン和の両方を意味する略記とする:
\begin{equation} \int \dd x \;(\cdots) \;\equiv\; \sum_{\sigma=\uparrow,\downarrow}\int \dd^3 r\;(\cdots) \label{eq:spin-int} \end{equation}$N$ 電子系の波動関数は $\Psi(x_1,x_2,\dots,x_N)$ と書く.これは $3N$ 個の連続変数と $N$ 個の2値変数の関数である.
ハミルトニアンがスピンに作用しないという事実は,無害どころか有用である.演算子 $\hat{H}_e$ は全スピン演算子 $\hat{\bm{S}}=\sum_i\hat{\bm{s}}_i$ の各成分と可換なので,固有関数は $\hat{H}_e$,$\hat{S}^2$,$\hat{S}_z$ の同時固有関数に選べる.すなわち「エネルギー $E$,全スピン $S$,その $z$ 成分 $M_S$」という量子数で状態を分類できる.この分類は第5章で開殻系を扱うときに本質的になる.逆に,スピン–軌道相互作用のような相対論的補正を加えるとこの可換性は失われ,スピンと軌道は分離できなくなる(§1.3.3で見たように重元素で重要になる効果である).
1.5.2 確率解釈・規格化・固有関数の直交性
Bornの確率解釈を多電子系に拡張する.$\abs{\Psi(x_1,\dots,x_N)}^2\,\dd^3r_1\cdots\dd^3r_N$ は,「電子1が $\rr_1$ のまわりの微小体積 $\dd^3r_1$ にスピン $\sigma_1$ で見出され,かつ電子2が $\rr_2$ のまわりに $\sigma_2$ で見出され,$\dots$」という同時確率である.全空間・全スピン配置にわたって足し合わせれば確率の総和は1でなければならないから,規格化条件
\begin{equation} \int \dd x_1\int \dd x_2\cdots\int \dd x_N\;\abs{\Psi(x_1,\dots,x_N)}^2 = 1 \label{eq:psi-norm} \end{equation}が課される.式 \eqref{eq:spin-int} の約束により,これは $3N$ 重の空間積分と $2^N$ 通りのスピン配置の和を意味する.この条件を満たす関数の集合(2乗可積分な関数の空間)が,固有値問題 \eqref{eq:elec-se} の解の候補である.
次に,異なるエネルギー固有値に属する固有関数どうしが直交することを確かめる.これは量子力学の一般論として既習かもしれないが,多電子系でも同じ論法がそのまま通用すること,とくに部分積分の境界項が落ちる理由を確認しておく価値がある(この直交性は§1.4の式 \eqref{eq:phi-ortho} で既に使った).
導出:$\hat{H}_e$ のエルミート性と固有関数の規格直交性
ステップ1(エルミート性).演算子 $\hat{A}$ がエルミートであるとは,任意の(規格化可能な)関数 $f,g$ に対して
\begin{equation} \int \dd x_1\cdots \dd x_N\; f^*\,(\hat{A}g) = \int \dd x_1\cdots \dd x_N\; (\hat{A}f)^*\,g \label{eq:herm-def} \end{equation}が成り立つことをいう.以下これを $\braket{f|\hat{A}g} = \braket{\hat{A}f|g}$ と略記する.$\hat{H}_e$ は「実関数の掛け算」と「Laplacian」からできている.前者(Coulombポテンシャル)は明らかにエルミートである:掛け算は順序を変えても同じで,ポテンシャルが実数なら複素共役をとっても変わらないからである.問題は後者である.
$-\tfrac12\nabla_i^2$ について式 \eqref{eq:herm-def} の左辺と右辺の差を計算する.まずベクトル解析の恒等式
$$ \nabla_i\cdot\big(f^*\nabla_i g\big) = \nabla_i f^*\cdot\nabla_i g + f^*\nabla_i^2 g $$を使う(積の微分則をベクトル場の発散に適用しただけである).これを $f^*\nabla_i^2 g$ について解いて,$i$ 番目の電子の座標について積分すると
\begin{align} \int \dd^3 r_i\; f^*\nabla_i^2 g &= \int \dd^3 r_i\;\nabla_i\cdot\big(f^*\nabla_i g\big) - \int \dd^3 r_i\;\nabla_i f^*\cdot\nabla_i g \notag\\ &= \oint_{S_\infty} \big(f^*\nabla_i g\big)\cdot \dd\bm{S} \;-\; \int \dd^3 r_i\;\nabla_i f^*\cdot\nabla_i g \label{eq:herm-lap} \end{align}となる.2行目ではGaussの発散定理を用いて体積積分を無限遠の閉曲面 $S_\infty$ 上の面積分に書き換えた.ここで境界項が消える理由を明示しておく.束縛状態の波動関数は $r_i\to\infty$ で指数関数的に減衰する($f\sim e^{-\kappa r_i}$,$\kappa=\sqrt{2\abs{E}}\gt 0$).一方,半径 $R$ の球面の面積は $4\pi R^2$ という多項式的な増大にすぎない.したがって面積分は $R^2 e^{-2\kappa R}\to 0$ となり,$R\to\infty$ で消える.指数関数的減衰が多項式的増大に必ず勝つ,というのがここでの要点である(散乱状態のように減衰しない関数を扱う場合には,周期境界条件を課して表面そのものを無くすという別の処方が必要になる.第3章・第12章で使う).
したがって
$$ \int \dd^3 r_i\; f^*\nabla_i^2 g = -\int \dd^3 r_i\;\nabla_i f^*\cdot\nabla_i g $$が成り立つ.同じ計算を $f$ と $g$ の役割を入れ替えて行うと $\int (\nabla_i^2 f)^* g = -\int \nabla_i f^*\cdot\nabla_i g$ となり(複素共役をとっても実の微分演算は変わらない),両者の右辺は完全に同一である.よって
$$ \int \dd^3 r_i\; f^*\nabla_i^2 g = \int \dd^3 r_i\;(\nabla_i^2 f)^*\,g $$すなわちLaplacianはエルミートであり,結局 $\hat{H}_e$ 全体がエルミートである.
ステップ2(固有値が実であること).$\hat{H}_e\Psi_n = E_n\Psi_n$,$\braket{\Psi_n|\Psi_n}=1$ とする.エルミート性を $f=g=\Psi_n$ に適用すると
$$ \braket{\Psi_n|\hat{H}_e\Psi_n} = E_n\braket{\Psi_n|\Psi_n} = E_n, \qquad \braket{\hat{H}_e\Psi_n|\Psi_n} = E_n^*\braket{\Psi_n|\Psi_n} = E_n^* $$となる(2つ目の式では,左側の関数に複素共役がかかるので $E_n^*$ が出る).両者は等しいので $E_n = E_n^*$,すなわち $E_n$ は実数である.エネルギーが実数であるという当たり前の事実が,エルミート性から保証されている.
ステップ3(直交性).$E_m \neq E_n$ の2つの固有関数を考える.エルミート性を $f=\Psi_m$,$g=\Psi_n$ に適用すると
\begin{align} E_n\braket{\Psi_m|\Psi_n} &= \braket{\Psi_m|\hat{H}_e\Psi_n} \notag\\ &= \braket{\hat{H}_e\Psi_m|\Psi_n} = E_m^*\braket{\Psi_m|\Psi_n} = E_m\braket{\Psi_m|\Psi_n} \label{eq:eig-real} \end{align}となる(最後の等号でステップ2の結果 $E_m^*=E_m$ を使った).移項して
$$ (E_n - E_m)\braket{\Psi_m|\Psi_n} = 0 $$を得る.仮定より $E_n-E_m\neq0$ だから $\braket{\Psi_m|\Psi_n}=0$ でなければならない.縮退している場合($E_m=E_n$,$m\neq n$)にはこの論法は使えないが,同じ固有値に属する固有関数の張る部分空間の中でGram–Schmidt直交化を行えばよい.以上より,固有関数系はつねに
\begin{equation} \braket{\Psi_m|\Psi_n} = \int \dd x_1\cdots\dd x_N\;\Psi_m^*\,\Psi_n = \delta_{mn} \label{eq:orthonormal} \end{equation}を満たすように選べる.これを規格直交性という.
1.5.3 不可弁別性から交換位相へ
ここからが本節の核心である.古典力学では,粒子に番号を貼り,その軌道を追跡することで個々の粒子を区別できる.量子力学ではそれができない.電子の波束は重なり合い,いったん重なった2つの電子のうち「どちらがどちらか」を原理的に決める方法は存在しない.すべての電子は完全に同一の性質(質量・電荷・スピン)を持ち,互いに交換しても物理的状況は変わらない——これを不可弁別性(indistinguishability)という.この物理的要請を式にする.
まず,ハミルトニアン自身が電子の番号の付け替えについて対称であることを確認しておく.2電子系(一般の $N$ でも同じである)で書き下すと,式 \eqref{eq:h-elec} は
\begin{equation} \hat{H}_e = -\frac{1}{2}\nabla_1^2 - \frac{1}{2}\nabla_2^2 + v_{\mathrm{ext}}(\rr_1) + v_{\mathrm{ext}}(\rr_2) + \frac{1}{\abs{\rr_1-\rr_2}} + \text{(定数)} \label{eq:h-symmetric} \end{equation}である($v_{\mathrm{ext}}(\rr) = -\sum_I Z_I/\abs{\rr-\bm{R}_I}$ は§1.2で導入した外部ポテンシャル,定数は核間反発).この式で番号 $1$ と $2$ をそっくり入れ替えてみると,第1項と第2項が入れ替わり,第3項と第4項が入れ替わり,第5項は $\abs{\rr_2-\rr_1}=\abs{\rr_1-\rr_2}$ なのでそのまま残る.全体としてはまったく同じ演算子に戻る.つまりハミルトニアンは電子の番号付けを知らない.番号は我々が計算の便宜のために勝手に貼ったラベルにすぎないのである.
そこで交換操作を演算子として定義する.交換演算子 $\hat{P}_{ij}$ を,波動関数の第 $i$ 引数と第 $j$ 引数を入れ替える操作
$$ \hat{P}_{ij}\,\Psi(\dots,x_i,\dots,x_j,\dots) \equiv \Psi(\dots,x_j,\dots,x_i,\dots) $$で定義する.定義から直ちに $\hat{P}_{ij}^2=\hat{1}$(2回入れ替えれば元に戻る)であり,上で確かめたハミルトニアンの対称性は $[\hat{P}_{ij},\hat{H}_e]=0$ と表現される.
不可弁別性の要請は,観測量である確率密度が交換で変わらないこと,すなわち
\begin{equation} \abs{\Psi(x_1,x_2,\dots,x_N)}^2 = \abs{\Psi(x_2,x_1,\dots,x_N)}^2 \label{eq:indist} \end{equation}と書かれる(表記を簡単にするため,以下は最初の2つの引数の交換で議論する.他の対でも同じである).もしこの等式が破れていれば,確率分布を精密に測定することで「1番の電子」と「2番の電子」に物理的な区別を与えられてしまい,電子が同種粒子であるという前提と矛盾する.ここから位相の制限が出る.
導出:交換位相は $\pm1$ に限られる
ステップ1.式 \eqref{eq:indist} は2つの複素数の絶対値が等しいことを述べている.絶対値が等しい2つの複素数は,絶対値1の位相因子だけ異なる.したがってある実数 $\alpha$ が存在して
\begin{equation} \Psi(x_1,x_2,\dots) = e^{i\alpha}\,\Psi(x_2,x_1,\dots) \label{eq:phase} \end{equation}と書ける.ここで $\alpha$ は座標に依存しない定数と仮定する.理由は,交換という操作が「どの場所で行うか」に依存しない粒子の種類そのものの性質だからである(この仮定をゆるめた場合については下の注意を見よ).
ステップ2.式 \eqref{eq:phase} は変数名の付け方によらない関係式だから,$x_1$ と $x_2$ の名前を入れ替えた形
$$ \Psi(x_2,x_1,\dots) = e^{i\alpha}\,\Psi(x_1,x_2,\dots) $$もまた成り立つ.これを式 \eqref{eq:phase} の右辺に代入すると,
\begin{align} \Psi(x_1,x_2,\dots) &= e^{i\alpha}\,\Psi(x_2,x_1,\dots) \notag\\ &= e^{i\alpha}\Big(e^{i\alpha}\,\Psi(x_1,x_2,\dots)\Big) \notag\\ &= e^{2i\alpha}\,\Psi(x_1,x_2,\dots) \label{eq:phase-twice} \end{align}となる.1行目から2行目へは上の関係式を代入しただけ,2行目から3行目へは指数法則 $e^{i\alpha}e^{i\alpha}=e^{2i\alpha}$ を使っただけである.
ステップ3.式 \eqref{eq:phase-twice} は「$\Psi$ を2回交換すると元の $\Psi$ に戻る」という自明な事実を式にしたものである.恒等的にゼロでない $\Psi$ について両辺を比べると
\begin{equation} e^{2i\alpha} = 1 \qquad\Longrightarrow\qquad e^{i\alpha} = \pm 1 \label{eq:pm-one} \end{equation}が従う.最後の含意は,$z\equiv e^{i\alpha}$ と置けば $z^2=1$,すなわち $(z-1)(z+1)=0$ から $z=\pm1$ が出ることによる.連続的に選べるはずだった位相が,2つの値のどちらかに離散化されてしまった点が重要である.
$e^{i\alpha}=+1$ を選ぶ粒子(交換について対称)をボソン,$e^{i\alpha}=-1$ を選ぶ粒子(交換について反対称)をフェルミオンと呼ぶ.どちらであるかは粒子の種類ごとに決まっており,実験事実として,また相対論的な場の量子論におけるスピン統計定理[5]の帰結として,半整数スピンを持つ粒子はフェルミオン,整数スピンを持つ粒子はボソンである.電子はスピン $1/2$ なのでフェルミオンであり,したがって電子系の波動関数は次を満たさなければならない.
すなわち任意の2電子の交換について符号が反転する.これが多電子波動関数に課される基本的な制限であり,式 \eqref{eq:elec-se} を解くときには,この対称性を持つ関数の中だけで解を探さなければならない.一般の並べ替え(置換)$P$ に対しては,$P$ を互換の積に分解したときの互換の個数の偶奇に応じて符号が決まり,
\begin{equation} \Psi(x_{P(1)},x_{P(2)},\dots,x_{P(N)}) = \mathrm{sgn}(P)\,\Psi(x_1,x_2,\dots,x_N) \label{eq:antisym-N} \end{equation}と書ける.この $\mathrm{sgn}(P)$ という記号の意味を,次の数学ノートで自己完結的に説明しておく.第5章のSlater行列式は,まさにこの符号つきの和として定義されるからである.
数学ノート:置換・互換・置換の符号
置換.$\{1,2,\dots,N\}$ からそれ自身への1対1写像 $P$ を置換という.$P$ は「$1$ を $P(1)$ に,$2$ を $P(2)$ に,$\dots$ 送る」という並べ替えの指示である.置換の全体は写像の合成を積として群をなし($N$ 次対称群 $S_N$),その要素数は $N!$ 個である($1$ の行き先が $N$ 通り,$2$ の行き先が残り $N-1$ 通り,$\dots$ と数えればよい).$N=10$ でも $10!=3\,628\,800$ 通りある.
互換.2つの要素だけを入れ替え,他は動かさない置換を互換(transposition)という.任意の置換は互換の積として書ける.手続きは単純で,「$1$ の行き先が $1$ でなければ,$1$ と $P(1)$ を入れ替える互換を掛けて $1$ を正しい位置に固定する.次に $2$ について同じことをする.$\dots$」という操作を繰り返せば,高々 $N-1$ 回の互換で恒等置換にたどり着く.
符号.分解の仕方は一意ではない(たとえば同じ互換を2回余分に掛けても結果は変わらない).しかし互換の個数の偶奇は分解の仕方によらず一定である.これは次のように分かる.多項式 $\displaystyle \Delta(t_1,\dots,t_N)=\prod_{i \lt j}(t_i-t_j)$ を考え,置換 $P$ に対して変数を $t_i \to t_{P(i)}$ と置き換える操作を考えると,因子 $(t_i-t_j)$ は別の因子 $\pm(t_k-t_l)$ に移るだけなので,$\Delta$ は $\pm\Delta$ に移る.1回の互換でこの符号がちょうど $-1$ 倍になることを直接確かめられるので,$P$ を $m$ 個の互換に分解すれば符号は $(-1)^m$ である.ところが $\Delta \to \pm\Delta$ という結果は $P$ だけで決まり分解には依らないから,$(-1)^m$ も分解に依らない.この値を $\mathrm{sgn}(P)\equiv(-1)^m$ と書き,置換の符号という.偶数個の互換に分解できる置換を偶置換($\mathrm{sgn}=+1$),奇数個なら奇置換($\mathrm{sgn}=-1$)と呼ぶ.
性質.定義から直ちに $\mathrm{sgn}(PQ)=\mathrm{sgn}(P)\,\mathrm{sgn}(Q)$,$\mathrm{sgn}(P^{-1})=\mathrm{sgn}(P)$ が従う($P$ と $Q$ の分解を並べればよい).行列式の定義 $\det A=\sum_P \mathrm{sgn}(P)\prod_i A_{i,P(i)}$ にこの符号が現れることを思い出せば,式 \eqref{eq:antisym-N} を満たす関数を作る自然な道具が行列式であることが予感できるだろう.これがSlater行列式(第5章)である.
数学ノート:位相を定数と仮定してよいか
式 \eqref{eq:phase} で位相 $\alpha$ を定数と仮定した点は,初学者が引っかかりやすい箇所である.仮に位相が座標に依存して $\alpha(x_1,x_2)$ であったとしても,ステップ2・3の議論から $\alpha(x_1,x_2)+\alpha(x_2,x_1)=2\pi n$($n$ は整数)が要求される.さらに,交換操作が状態の重ね合わせ(量子力学の線形性)と両立するためには,$\alpha$ はどの状態に対しても共通の定数でなければならない.加えて,3次元空間では2粒子を入れ替える経路が連続変形で「ほどける」ため,位相は経路によらず $\pm1$ に限られることが位相幾何学的に示される.2次元空間ではこの「ほどけ」が起こらず,$e^{i\alpha}$ が $\pm1$ 以外の値を取る粒子(エニオン)が原理的に許される.分数量子Hall効果の準粒子がその実例だが,本書が扱う3次元のバルク電子系ではつねに $-1$ である.
1.5.4 反対称性の直接的な帰結
式 \eqref{eq:antisym2} は単なる符号の約束ではない.そこから物質の構造を決定づける結論が直ちに出る.
Pauliの排他律.式 \eqref{eq:antisym2} で $x_i = x_j = x$,すなわち2つの電子が空間座標もスピンもまったく同じ状態にある場合を考える.左辺と右辺は同じ関数値を指すから
\begin{equation} \Psi(\dots,x,\dots,x,\dots) = -\,\Psi(\dots,x,\dots,x,\dots) \quad\Longrightarrow\quad 2\Psi(\dots,x,\dots,x,\dots)=0 \quad\Longrightarrow\quad \Psi(\dots,x,\dots,x,\dots)=0 \label{eq:pauli-zero} \end{equation}が従う(右辺を左辺に移項して $2\Psi=0$,両辺を2で割った).すなわち2個の電子が同一の合成変数 $x$ を持つ確率はゼロである.これがPauliの排他律の波動関数レベルでの表現である.よく「同じ量子状態に2個の電子は入れない」と述べられるが,その原型はこの2行の計算に尽きている.第5章では,この帰結が「同一のスピン軌道を2重に占有できない」という軌道描像の言明に翻訳されることを見る.
交換ホール.式 \eqref{eq:pauli-zero} は $x_i=x_j$ という1点だけの話ではない.波動関数は連続関数だから,$x_i \to x_j$ の近傍でも $\Psi$ は小さくなければならない.合成変数が一致するとは「同じ位置」かつ「同じスピン」ということなので,結論は同じスピンを持つ2電子は互いに近づけないとなる.各電子のまわりには,同スピン電子の存在確率が抑えられた領域——交換ホール(exchange hole, Fermiホール)——が付き従うのである.これは電荷どうしのCoulomb反発とはまったく起源の異なる「実効的な反発」であり,純粋に波動関数の対称性から生じる.交換ホールは電子間反発のエネルギーを下げる方向に働き,その利得が交換エネルギーである.定量的な扱いは第5章(Hartree–Fock交換項)と第7章(一様電子ガスの交換エネルギー)で行う.なお,スピンが反平行($\sigma_i\neq\sigma_j$)なら $x_i\neq x_j$ なので式 \eqref{eq:pauli-zero} は何も要求しない.反平行スピンの2電子が同じ場所にいることは,反対称性の観点からは許される(実際にはCoulomb反発によって避け合う.この効果が相関であり,第7章・第11章の主題となる).
物理的意味:反対称性が物質の体積を作る
反対称性という一見形式的な要請が,我々の身のまわりの物質の最も基本的な性質——物質が有限の体積を持ち,押しても簡単には潰れないこと——を支えている.もし電子がボソンであったなら,すべての電子が最低エネルギーの1つの軌道に落ち込み,原子は殻構造を失い,周期表も化学結合も存在しない.反対称性が電子を次々と上のエネルギー準位へ追いやることで,原子には殻ができ,元素に個性が生まれ,化学が可能になる.金属中の自由電子が室温の熱エネルギーよりはるかに大きい Fermiエネルギー(数 eV)を持つのも同じ理由であり,白色矮星や中性子星を重力崩壊から支える縮退圧もまた同じ起源である.式 \eqref{eq:antisym2} のマイナス符号1つが,これらすべての背後にある.
1.5.5 電子密度の定義
本節の最後に,DFTの主役となる量——電子密度——を多電子波動関数から定義しておく.ここでも不可弁別性が本質的な役割を果たす.
まず「$i$ 番目の電子を合成変数 $x$ の状態に見出す確率密度」を,他の電子の変数をすべて積分して消すことで定義する:
$$ P_i(x) \equiv \int \dd x_1\cdots\dd x_{i-1}\,\dd x_{i+1}\cdots \dd x_N\; \abs{\Psi(x_1,\dots,x_{i-1},x,x_{i+1},\dots,x_N)}^2 $$しかし電子には番号がないのだから,$P_i$ そのものは観測できる量ではない.観測できるのは「番号を問わず,そこに電子が見出される確率密度」であり,それは各 $i$ についての和
$$ n(x) = \sum_{i=1}^{N} P_i(x) $$である.ここで反対称性が効く.式 \eqref{eq:antisym2} より $\abs{\Psi}^2$ は任意の引数の入れ替えについて不変である(符号 $-1$ が2乗で消える).したがって $P_i(x)$ の被積分関数において第1引数と第 $i$ 引数を入れ替え,さらに積分変数の名前を付け替えると,$P_i(x)$ は $P_1(x)$ と完全に一致する.$N$ 個の項がすべて等しいので和は単に $N$ 倍になり,
$$ n(x) = N\,P_1(x) = N\int \dd x_2\cdots \dd x_N\;\abs{\Psi(x,x_2,\dots,x_N)}^2 $$を得る.最後に,スピンを区別せず空間分布だけを見るためにスピン和を取れば,電子密度が定義される.
定義:電子密度
\begin{equation} \rho(\rr) \;\equiv\; N \sum_{\sigma}\int \dd x_2\cdots \dd x_N\; \abs{\Psi\big((\rr,\sigma),x_2,\dots,x_N\big)}^2 \label{eq:rho-def} \end{equation}$\rho(\rr)\,\dd^3r$ は,微小体積 $\dd^3 r$ の中に見出される電子の個数の期待値である.スピン和を取らずに $\sigma$ を固定したものをスピン密度 $\rho_\sigma(\rr)$ と呼び,$\rho=\rho_\uparrow+\rho_\downarrow$,$m=\rho_\uparrow-\rho_\downarrow$(磁化密度)と分解する.
この定義が電子数を正しく数えていることを確かめよう.式 \eqref{eq:rho-def} を全空間で積分すると
\begin{align} \int \dd^3 r\;\rho(\rr) &= N\sum_\sigma \int \dd^3 r\int \dd x_2\cdots\dd x_N\;\abs{\Psi\big((\rr,\sigma),x_2,\dots,x_N\big)}^2 \notag\\ &= N\int \dd x_1 \int \dd x_2\cdots\dd x_N\;\abs{\Psi(x_1,x_2,\dots,x_N)}^2 \notag\\ &= N\cdot 1 = N \label{eq:rho-norm} \end{align}となる.1行目から2行目へは,空間積分とスピン和をまとめて式 \eqref{eq:spin-int} の $\int\dd x_1$ に書き直しただけである.2行目から3行目へは規格化条件 \eqref{eq:psi-norm} を使った.したがって $\int\rho\,\dd^3r=N$,すなわち電子密度の積分は電子数に等しい.当然の結果だが,この関係式は後にDFTで粒子数を保つ拘束条件として,またLagrange未定乗数(化学ポテンシャル)を導入する足場として繰り返し使われる(第4章,第9章,第10章).
ここで立ち止まって,$\Psi$ と $\rho$ の情報量の差を意識しておきたい.$\Psi$ は $3N$ 個の連続変数の複素関数であるのに対し,$\rho$ は電子が何個あろうと3変数の実関数である.$\rho$ は $\Psi$ から作られる派生量にすぎず,明らかに情報を失っている.にもかかわらず,基底状態に限れば $\rho$ が系のすべてを決めてしまう——これがHohenberg–Kohnの定理(第9章)であり,DFTが成立する理由である.その意味を実感するために,次節では $\Psi$ を直接扱う方針がどれほど絶望的かを数値で確かめる.
1.6 なぜ直接解けないのか:次元の呪い
Born–Oppenheimer近似によって,解くべき問題は電子だけの固有値問題 \eqref{eq:elec-se} に還元された.核の自由度が消え,方程式は $3N$ 個の変数の偏微分方程式になった.ここまでくれば,あとは計算機に任せて数値的に解けばよさそうに思える.ところが,そうはいかない.本節では「解けない」という言葉の意味を,精神論ではなく具体的な数値で確かめる.結論を先に述べておくと,困難は方程式が数学的に難解であることではなく,解を書き留めるために必要な情報量そのものが,物理的に利用可能な資源を天文学的に超えてしまうことにある.
1.6.1 波動関数は $3N$ 次元空間の関数である
最初に,初学者が最も誤解しやすい点を確認しておく.多電子波動関数 $\Psi(x_1,x_2,\dots,x_N)$ は,「3次元空間に置かれた $N$ 個の関数」ではない.それは $3N$ 次元の配位空間(configuration space)の上で定義された,たった1個の関数である.$N=1$ ならこれは見慣れた3次元の軌道だが,$N=2$ で早くも6次元,$N=10$ で30次元の関数になる.
もし電子間反発 $\hat{V}_{ee}$ が存在しなければ,事情はまったく違っていた.$\hat{V}_{ee}=0$ のときハミルトニアンは1電子ハミルトニアンの単純な和 $\hat{H}=\sum_i \hat{h}(x_i)$ になり,変数分離ができて解は積の形
$$ \Psi(x_1,\dots,x_N) = \varphi_1(x_1)\varphi_2(x_2)\cdots\varphi_N(x_N) $$に書ける(反対称化はまだ考えないでおく).このとき必要な情報は「3次元関数 $N$ 個」であり,$N$ に対して線形にしか増えない.ところが $\hat{V}_{ee}$ があると,ある電子がどこにいるかが他のすべての電子の位置に依存するようになり,この積の形が壊れる.$\Psi$ は $3N$ 個の変数が分離不可能に絡み合った関数となり,必要な情報量は $N$ に対して指数関数的に増える.これが多体問題の本質的な困難であり,$\hat{V}_{ee}$ を「元凶」と呼んだ(§1.2.2)理由である.
1.6.2 実空間グリッドによる見積り
情報量を具体的に数えるために,最も素朴な数値表現を考える.1個の電子の座標を,一辺 $L$ の立方体の中の等間隔グリッドで表すことにし,各デカルト成分を $p$ 点で離散化する.すると1電子あたり $p^3$ 個のグリッド点があり,$N$ 電子の配位空間では
\begin{equation} (\text{格子点の総数}) = \underbrace{p^3\times p^3\times\cdots\times p^3}_{N\ \text{個}} = p^{3N} \label{eq:grid-count} \end{equation}個の点に,それぞれ $\Psi$ の値を1つずつ記録しなければならない.倍精度実数1個は8バイトなので,必要な記憶容量は
\begin{equation} (\text{必要バイト数}) = 8\times p^{3N}\ \mathrm{B} \label{eq:grid-bytes} \end{equation}である.指数の肩に $N$ が乗っている点に注意してほしい.$p^{3N} = e^{3N\ln p}$ と書けば,これが $N$ の指数関数であることがはっきりする.電子を1個増やすたびに必要量は $p^3$ 倍になるのである.
例:$p=10$ のグリッドで水分子の波動関数を保存する
まず,この見積りがどれほど甘いかを確認する.$p=10$ とし,原子や小分子が収まる一辺 $L=10\,a_0$($\simeq 5.3$ Å)の箱を取ると,グリッド間隔は $h = L/p = 1\,a_0 = 0.53$ Å である.電子密度は核の近傍で $0.1\,a_0$ 程度のスケールで激しく変化するから,この刻みは実用にはまったく足りない.それでもなお,以下の数値になる.
水分子 H$_2$O の電子数は $N = 8+1+1 = 10$ である.式 \eqref{eq:grid-count} より格子点数は
$$ p^{3N} = 10^{3\times 10} = 10^{30} $$個.式 \eqref{eq:grid-bytes} より必要な記憶容量は
$$ 8\times 10^{30}\ \mathrm{B} = \frac{8\times10^{30}}{10^{12}}\ \mathrm{TB} = 8\times 10^{18}\ \mathrm{TB} \simeq 10^{19}\ \mathrm{TB} $$である($1\,\mathrm{TB}=10^{12}$ Bを使って割り算しただけである).この数がどれほど非常識かを比較で示す.
- 現在の最高性能クラスのスーパーコンピュータの主記憶は $10$ PB $=10^{16}$ B 程度である.必要量はその $8\times10^{30}/10^{16} = 8\times10^{14}$ 倍,すなわちスーパーコンピュータ約800兆台分である.
- 人類が保有するデジタルデータの総量は $10^{23}$ B のオーダー(数百ゼタバイト)と見積もられている.必要量はその $8\times10^{30}/10^{23} \simeq 10^{8}$ 倍,1億倍である.
- 仮に記憶できたとしても,ハミルトニアンを1回作用させるだけで格子点数と同じオーダー,すなわち $10^{30}$ 回程度の演算が要る.$10^{18}$ 演算毎秒(エクサスケール)の計算機を使っても $10^{30}/10^{18} = 10^{12}$ 秒,1年 $=3.16\times10^{7}$ 秒で割って約3万年かかる.しかも固有値問題を解くにはこれを何度も繰り返す必要がある.
さらに,実用的な精度を得るには $h\simeq 0.1\,a_0$,すなわち $p\simeq 100$ が必要である.すると格子点数は $100^{3N} = 10^{6N}$ となり,$N=10$ では $10^{60}$ 個——上の数字をもう一度30桁上回る.スピン配置の数 $2^N = 1024$ 通りを掛ける効果など,この爆発の前ではまったく問題にならない.
電子数を変えたときの様子を表1.2にまとめる.$p=10$ という極端に粗いグリッドでも,実用的に保存できるのは電子数4〜5個までであることが読み取れる.ヘリウム原子(2電子)やリチウム原子(3電子)ならこの素朴な方法でも解けるし,実際そのような直接数値解は高精度参照値として作られている.しかし炭素原子(6電子)ですでに現実的でなくなり,水分子(10電子)は絶望的である.
| 電子数 $N$ | 対応する系の例 | 格子点数 $10^{3N}$ | 必要記憶容量 | 実現可能性 |
|---|---|---|---|---|
| 1 | H 原子 | $10^{3}$ | 8 kB | 自明 |
| 2 | He 原子,H$_2$ 分子 | $10^{6}$ | 8 MB | ノートPCで可能 |
| 3 | Li 原子 | $10^{9}$ | 8 GB | PCで可能 |
| 4 | Be 原子 | $10^{12}$ | 8 TB | 大型計算機なら可能 |
| 5 | B 原子 | $10^{15}$ | 8 PB | 最大級のスパコンの限界 |
| 6 | C 原子 | $10^{18}$ | 8 EB | 不可能 |
| 10 | H$_2$O 分子 | $10^{30}$ | $8\times10^{18}$ TB | 絶望的 |
| 20 | Ca 原子 | $10^{60}$ | — | — |
| 100 | 小さな有機分子 | $10^{300}$ | — | — |
最後の行を味わってほしい.観測可能な宇宙に存在する原子の総数は $10^{80}$ 個程度と見積もられている.原子1個を1バイトの記憶素子に使えたとしても,100電子系の波動関数を書き留めるには220桁足りない.これは技術の進歩で埋まる差ではない.Moore則に従って計算機の性能が2年で2倍になったとしても,扱える電子数は $p^3=10^3$ 倍ごとに1個増えるだけで,電子を1個増やすのに約20年かかる計算になる.
1.6.3 基底関数展開でも事情は変わらない:組合せ論的爆発
「グリッドが素朴すぎるのだ,賢い基底関数を使えばよい」と考えるのは自然である.実際,量子化学の標準的なやり方はそうする.しかし結論から言えば,指数関数的増大は基底の取り方を変えても消えない.それは表現方法の問題ではなく,問題そのものの性質だからである.
典型的な手続きはこうである.まず1電子の状態を張る基底として $M$ 個のスピン軌道 $\{\varphi_1,\dots,\varphi_M\}$ を用意する($M \gt N$).次に,そこから $N$ 個を選んで反対称な $N$ 電子関数(Slater行列式,第5章)を作る.1つの選び方 $I=\{i_1 \lt i_2 \lt \dots \lt i_N\}$ が1つの行列式 $\Phi_I$ を与える.反対称関数の完全系はこれらの行列式で張られるので,厳密な波動関数は
\begin{equation} \Psi = \sum_{I} C_I\,\Phi_I \label{eq:ci-expand} \end{equation}と展開できる.これを配置間相互作用(configuration interaction; CI)展開といい,すべての $I$ を含めたものを完全CI(full CI)と呼ぶ.完全CIは,与えられた基底の範囲内では厳密解である.
問題は展開の項数である.$M$ 個から $N$ 個を選ぶ組合せの数だから,
\begin{equation} D = \binom{M}{N} = \frac{M!}{N!\,(M-N)!} \label{eq:ci-dim} \end{equation}個の行列式が必要になる.$D$ はCI展開の次元であり,同時にこの基底で表したハミルトニアン行列 $H_{IJ}=\braket{\Phi_I|\hat{H}_e|\Phi_J}$ の次元でもある.具体的な数値を見よう.
例:完全CIの次元
- 水分子,cc-pVDZ基底.この基底は H$_2$O に対して24個の空間軌道を与える.電子は10個で,$\uparrow$ スピン5個,$\downarrow$ スピン5個に分かれる.スピンごとに独立に選ぶので,行列式の数は $$ D = \binom{24}{5}\times\binom{24}{5} = 42\,504^2 = 1.81\times10^{9} $$ すなわち約18億個である.倍精度の係数 $C_I$ を保存するだけで $8\times1.81\times10^9 = 14.5$ GB になる.これは実際に計算されたことのある規模だが,基底関数としては最小限のものであり,この基底での完全CIエネルギーはまだ厳密解から数 kcal/mol ずれている.
- もう少し大きな系.$M=100$ 個のスピン軌道に $N=10$ 電子なら $\displaystyle\binom{100}{10}=1.73\times10^{13}$,$M=200$,$N=20$ なら $\displaystyle\binom{200}{20}=1.61\times10^{27}$ である.
- ベンゼン C$_6$H$_6$.電子数42,中程度の基底で102個の空間軌道を取ると,$\displaystyle\binom{102}{21}^2 \simeq 1.1\times10^{43}$ 個の行列式になる.
さらに,ハミルトニアン行列を密行列として保存するには $D^2$ 個の要素が要る.$D=10^{9}$ でも $D^2=10^{18}$ 個であり,これも不可能である.実務ではDavidson法のような反復対角化を使い,行列を保存せずに「行列とベクトルの積」だけを計算することで $O(D)$ の記憶量に抑える.しかしその $D$ 自身が指数関数的に増える以上,壁は動かない.
数学ノート:二項係数の増大とStirlingの公式
$\displaystyle\binom{M}{N}$ が $N$ に対して指数関数的に増えることを,きちんと示しておこう.最も「けち」な基底,すなわち電子数のちょうど2倍しかスピン軌道を用意しない場合($M=2N$)を考える.これは基底を減らせるだけ減らした状況である.
階乗の漸近形としてStirlingの公式
\begin{equation} n! \simeq \sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}\qquad(n\gg1) \label{eq:stirling} \end{equation}を使う(証明は $\ln n! = \sum_{k=1}^n \ln k \simeq \int_1^n \ln x\,\dd x = n\ln n - n + 1$ という積分近似が骨格で,係数 $\sqrt{2\pi n}$ はより精密な評価から出る).これを $\displaystyle\binom{2N}{N} = \frac{(2N)!}{(N!)^2}$ に代入する.分子と分母を別々に書くと
\begin{align} \binom{2N}{N} &\simeq \frac{\sqrt{2\pi (2N)}\,\big(2N/e\big)^{2N}}{\Big[\sqrt{2\pi N}\,\big(N/e\big)^{N}\Big]^{2}} \notag\\[2pt] &= \frac{\sqrt{4\pi N}\;2^{2N}N^{2N}e^{-2N}}{2\pi N\;N^{2N}e^{-2N}} \notag\\[2pt] &= \frac{2\sqrt{\pi N}\;4^{N}}{2\pi N} = \frac{4^{N}}{\sqrt{\pi N}} \label{eq:central-binom} \end{align}1行目から2行目へは,分子で $(2N)^{2N}=2^{2N}N^{2N}$ と分解し,分母では2乗を展開して $(\sqrt{2\pi N})^2 = 2\pi N$,$\big[(N/e)^N\big]^2 = N^{2N}e^{-2N}$ とした.2行目から3行目へは $N^{2N}e^{-2N}$ を約分し,$\sqrt{4\pi N}=2\sqrt{\pi N}$,$2^{2N}=4^N$ と書き直し,最後に $2\sqrt{\pi N}/(2\pi N) = 1/\sqrt{\pi N}$ とまとめた.
結果は明快である.多項式因子 $1/\sqrt{\pi N}$ を無視すれば $\displaystyle\binom{2N}{N}\sim 4^{N}$,すなわち電子を1個増やすごとに次元が4倍になる.基底を大きくすればこの底はさらに大きくなる.$N=20$ での近似値は $4^{20}/\sqrt{20\pi}=1.39\times10^{11}$ で,厳密値 $\displaystyle\binom{40}{20}=1.38\times10^{11}$ とよく一致する.
1.6.4 指数の壁
Kohnはノーベル賞受賞講演[6]でこの状況を指数の壁(exponential wall)と呼び,多体波動関数は電子数が $10^{3}$ を超えるような系では「正当な科学的概念ではない」とまで述べた.原理的には存在する量でも,書き留めることも検証することもできないなら,科学の道具として機能しないという趣旨である.図1.3に,この壁の高さと,後で見るDFTの計算量($O(N^3)$)との対比を示す.
もう一つ,事態を厳しくしている要素がある.要求される相対精度である.ベンゼン分子の全電子エネルギーは約 $-232\,E_h$ である.一方,化学的に意味のある議論をするには「化学精度」$1\ \mathrm{kcal/mol}=1.59\times10^{-3}\,E_h$(§1.3.3の例)が必要である.両者の比を取ると
$$ \frac{1.59\times10^{-3}}{232} = 6.9\times10^{-6} $$となり,全エネルギーを相対精度 $10^{-6}$ 程度で求めなければならない.しかも問題を難しくしているのは,この微小なエネルギー差が「たまたま小さい」のではなく,大きな正の量(運動エネルギーと電子間反発)と大きな負の量(電子–核引力)の差として現れることである.巨大な数どうしの引き算で微小な差を正確に出す——これが第一原理計算の数値的な難しさの正体である.
物理的意味:「解けない」とはどういうことか
ここで言う「解けない」は,方程式の解が存在しないという意味でも,数学的に難解だという意味でもない.式 \eqref{eq:elec-se} の解は確かに存在し,原理的には一意に定まっている.問題は,その解を書き留めることも,他人に伝えることも,正しいかどうか検証することもできないという点にある.$10^{300}$ 個の数値からなる対象は,たとえ神が教えてくれたとしても人類には受け取る場所がない.したがって多体問題の研究とは,「厳密解を求める努力」ではなく,「知りたい物理量を,扱える大きさの情報から取り出す方法を設計する努力」である.この視点の転換こそが,次節で分類する諸手法の共通の出発点であり,とりわけDFTの発想の核心である.$\Psi$ が $3N$ 変数であるのに対し,我々が実際に知りたい量——全エネルギー,電子密度,力,状態密度——はいずれもはるかに小さな情報量しか持たない.ならば,$\Psi$ を経由せずに直接それらを求める道はないか.この問いへの一つの答えがDFTである.
指数の壁を回避する戦略は,大きく3つに分類できる.(i) 基本変数を取り替える:$\Psi$ の代わりに電子密度 $\rho(\rr)$(3変数)や1体Green関数 $G(\rr,\rr';\omega)$(7変数)を基本変数に据える.DFTと多体グリーン関数法がこれにあたる.(ii) 確率的にサンプリングする:$3N$ 次元空間を全部覆うのを諦め,乱数で重要な領域だけを拾う.量子モンテカルロ法である.(iii) 波動関数の形を賢く制限する:$\Psi$ の関数形にパラメータの少ない構造を仮定し,その範囲で最良の解を探す.Hartree–Fock法,結合クラスター法,テンソルネットワーク法がこれである.次節でこれらを整理し,本書がDFTを選ぶ理由を述べる.
1.7 第一原理計算手法の分類
指数の壁を前にして,20世紀後半の物理学と化学はいくつもの迂回路を開拓してきた.本節ではそれらを俯瞰し,それぞれが「何を犠牲にして何を得ているか」を明確にする.分類の軸は単純で,何を基本変数に選ぶかである.この選択が計算量と精度のすべてを決める.ここで得られる見取り図は,本書のどの章がどこに位置するかを示す地図にもなる.
1.7.1 波動関数理論
基本変数を波動関数 $\Psi$ のままにし,その関数形に構造を仮定して変分的に最良の解を探す立場である.歴史的にはこれが最初に発展し,量子化学の主流をなしてきた.
Hartree–Fock(HF)法.$\Psi$ を1個のSlater行列式に限定する.すなわち「各電子は,他の電子が作る平均的な場の中を独立に動く」という平均場近似である.$N$ 電子の問題が $N$ 個の1電子方程式に帰着するので,計算量は基底関数の数 $M$ に対して $O(M^4)$(2電子積分の数),実効的には $O(N^3)$–$O(N^4)$ にとどまる.反対称性は厳密に満たされるので交換効果は正確に取り込まれるが,それ以外の電子どうしの避け合い——相関——が完全に欠落する.この欠落分
$$ E_{\mathrm{corr}} \equiv E_{\text{厳密}} - E_{\mathrm{HF}} $$を相関エネルギーと呼ぶ.相関エネルギーは全エネルギーの1%程度と小さいが,結合エネルギーや反応エネルギーと同程度の大きさ(分子で数 eV)を持つため,無視することはできない.HF法の詳細は第5章で全面的に扱う.
post-HF法.HFの行列式を出発点に,励起配置を系統的に加えて相関を取り戻す.取り込む励起の範囲によって階層があり,2次摂動論(MP2,$O(N^5)$),結合クラスター法(CCSD,$O(N^6)$;CCSD(T),$O(N^7)$),そして完全CI($O(e^N)$)と続く.CCSD(T)は小分子で化学精度を達成でき,「黄金律」と呼ばれる基準手法である[10].この階層の長所は系統的改善が可能であること,すなわち計算資源を注げば厳密解に近づけることが保証されている点である.短所は,そのコストが $N$ の高いべき,最終的には指数関数で増えることであり,扱えるのはせいぜい数十原子までである.
1.7.2 密度汎関数理論
基本変数を波動関数から電子密度 $\rho(\rr)$ に取り替える立場である.式 \eqref{eq:rho-def} で定義した $\rho$ は,電子が何個あろうと3変数の実関数であり,情報量は $N$ に依存しない.この一点だけで,指数の壁は原理的に消える.
ただしこれは「情報を捨てて楽になった」のではない.基底状態に限れば $\rho$ が系のすべての性質を決定するという定理(Hohenberg–Kohnの定理[7],第9章)が背後にあり,原理的には厳密な理論である.難しさは別の場所に移る:$\rho$ からエネルギーを与える汎関数 $E[\rho]$ の正確な形が誰にも分からないのである.実用上は近似汎関数を用いるため,DFTは「原理的には厳密,実際には近似に依存する」理論となる.計算量は $O(N^3)$(固有値問題の対角化が律速),工夫すれば $O(N)$ まで下げられ,数百から数万原子の系が扱える.精度は近似汎関数の質で決まり,標準的なLDA/GGAでは格子定数が1%程度,凝集エネルギーが $0.1$–$0.3$ eV程度の誤差になる.バンドギャップの系統的な過小評価など,既知の弱点もある(第11章).
1.7.3 量子モンテカルロ法
$3N$ 次元の積分を,格子で覆うのではなく乱数で標本抽出する.多次元積分の誤差が次元によらず標本数 $M_s$ だけで決まり $\propto 1/\sqrt{M_s}$ となる,というモンテカルロ法の一般的性質を利用する立場である.変分モンテカルロ(VMC)は試行波動関数の期待値を確率的に評価し,拡散モンテカルロ(DMC)は虚時間発展によって基底状態を射影する[11].計算量は $O(N^3)$–$O(N^4)$ 程度で,並列化効率が非常に高い.
最大の障害は符号問題である.フェルミオンの波動関数は反対称性のため正負両方の値を取り,確率分布として直接扱えない.実用上は波動関数の節面(ゼロになる面)をあらかじめ固定する固定節近似で回避するが,これは制御しにくい誤差を残す.それでも精度は高く,DFTの近似汎関数を検証するための参照データ(たとえば一様電子ガスの相関エネルギー,第7章)を供給する重要な役割を担っている.
1.7.4 多体グリーン関数法
基本変数を1体Green関数 $G(\rr,\rr';\omega)$ に取る立場である.$G$ は7変数の関数で,$\Psi$ よりはるかに小さく,$\rho$ よりは大きい.この中間的な選択の見返りとして,$\rho$ からは直接得にくい励起状態の情報——電子を1個付け加える/取り去るのに要するエネルギー,すなわち準粒子スペクトル——が自然に得られる.摂動展開の最低次を取ったものがGW近似[12]であり,遮蔽されたCoulomb相互作用 $W$ を使うことで金属でも収束する展開になっている(遮蔽の理論は第8章で扱う).光吸収スペクトルには電子–正孔相互作用を含むBethe–Salpeter方程式(BSE)を組み合わせる[13].計算量は $O(N^3)$–$O(N^4)$ だが前係数が大きく,通常はDFTの結果を出発点とした補正として使われる.バンドギャップの誤差は $0.1$–$0.3$ eV程度まで改善する.
1.7.5 比較と本書の選択
以上を表1.3にまとめ,図1.4に「扱える系の大きさ」と「精度」の平面上での位置関係を示す.
| 手法 | 基本変数 | 計算量 | 典型的な精度 | 得意分野 | 主な限界 |
|---|---|---|---|---|---|
| Hartree–Fock | 1個のSlater行列式 | $O(N^3)$–$O(N^4)$ | 結合エネルギー誤差 $\sim1$ eV(相関を欠く) | 定性的描像,post-HFの出発点 | 相関エネルギーがゼロ,金属で破綻 |
| MP2 / CCSD(T) | 励起配置を加えた波動関数 | $O(N^5)$ / $O(N^7)$ | $\sim1$ kcal/mol(化学精度) | 小分子の熱化学・反応 | 数十原子が上限,周期系に不向き |
| 完全CI | 全Slater行列式の重ね合わせ | $O(e^{N})$ | 基底の範囲で厳密 | ベンチマーク | 電子10個・軌道20個程度が限界 |
| 密度汎関数理論 | 電子密度 $\rho(\rr)$(3変数) | $O(N^3)$(線形法で $O(N)$) | 格子定数 $\sim1\%$,凝集エネルギー $0.1$–$0.3$ eV | 固体・表面・分子動力学,数百〜数万原子 | 汎関数の近似に依存,ギャップの過小評価,強相関系 |
| 量子モンテカルロ | 確率標本としての $\abs{\Psi}^2$ | $O(N^3)$–$O(N^4)$,高並列 | 凝集エネルギー $\sim0.05$ eV | 参照データ,強相関系 | 符号問題・固定節近似,力の計算が困難 |
| 多体グリーン関数法(GW/BSE) | 1体Green関数 $G$(7変数) | $O(N^3)$–$O(N^4)$(大きな前係数) | バンドギャップ $0.1$–$0.3$ eV | 準粒子スペクトル,光学応答 | 出発点(通常DFT)に依存,全エネルギーは苦手 |
図1.4から読み取れる傾向は明快である.手法は左上(小さい系・高精度)から右下(大きい系・低精度)へ並んでおり,精度と系のサイズはトレードオフの関係にある.本書がDFTを軸に据えるのは,次の3つの理由による.
- 計算量が多項式的である.$O(N^3)$,工夫すれば $O(N)$ であり,指数の壁を原理的に回避している.数百から数万原子という,物質科学で実際に問題になる規模を扱える唯一の第一原理手法である.
- 原理は厳密である.後で見るように,DFTは近似から出発する理論ではない.まず厳密な定理(第9章)を立て,そのうえで未知の部分(交換相関エネルギー)だけを近似する.したがって「どこが近似でどこが厳密か」が明確であり,改善の方向も定まっている.
- 物質科学の共通言語である.バンド構造,状態密度,力と応力,構造最適化,分子動力学,分光スペクトルといった実務上の主要な出力が,すべてこの枠組みの中で計算できる(第IV部〜第VI部).
もちろん万能ではない.強相関系(遷移金属酸化物,$f$ 電子系),van der Waals力が支配する系,励起状態などでは標準的な近似汎関数が破綻し,他の手法との組み合わせが必要になる.その限界も含めて理解することが本書の目標である.
1.7.6 DFTの2本柱(予告)
以後の章の地図として,DFTの骨格をここで一度だけ提示しておく.証明も導出も行わない.式の形と用語に慣れ,「どの章で何が正当化されるのか」を把握するのが目的である.DFTは次の2つの柱で立っている.
第1の柱:全エネルギーは電子密度の汎関数である.Hohenberg–Kohnの定理[7](第9章)は,基底状態の電子密度 $\rho(\rr)$ が外部ポテンシャル $v_{\mathrm{ext}}$ を(定数分を除いて)一意に決めることを主張する.$v_{\mathrm{ext}}$ が決まればハミルトニアン \eqref{eq:h-elec} が決まり,その基底状態も全物性も決まる.したがって全エネルギーは $\rho$ だけの関数,すなわち汎関数(関数を入力として数を返す写像.第4章で数学的に定義する)として書ける.
第2の柱:多体効果を一箇所に押し込め,一体問題に見せかける.Kohn–Shamの構成[8](第10章)は,「相互作用のない $N$ 電子系で,その密度が真の系の密度に一致するもの」を導入する.厄介な多体効果はすべて交換相関エネルギー $E_{xc}[\rho]$ という1つの項に押し込め,残りは厳密に計算できる形にする.この結果,解くべき方程式は1電子のSchrödinger型方程式になる.
この2本柱を式にすると,次の2つになる.まず全エネルギー汎関数である.
各項の意味は次の通りである.$T_s[\rho]$ は相互作用のない参照系の運動エネルギーで,真の運動エネルギー $T$ とは異なる(その差は $E_{xc}$ に吸収される).第2項は外部ポテンシャル(原子核が作るCoulomb場)によるエネルギーで,$\rho$ の1次の汎関数である.第3項はHartreeエネルギー(古典的な静電エネルギー)
\begin{equation} E_{\mathrm{H}}[\rho] = \frac{1}{2}\int \dd^3 r\int \dd^3 r'\;\frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}} \label{eq:ehartree} \end{equation}である(係数 $1/2$ は,式 \eqref{eq:pair-sym} と同じく電子対を2度数えないためのものである).そして第4項 $E_{xc}[\rho]$ が交換相関エネルギー,すなわちそれ以外のすべてである.定義上は「厳密なエネルギーから他の3項を引いた残り」であり,その正確な形は知られていない.DFTの実用上の全問題はこの1項に集約されている.近似の階層(LDA,GGA,ハイブリッド)は第11章の主題である.なお本章では外部ポテンシャルを $v_{\mathrm{ext}}$,交換相関エネルギーを $E_{xc}$ と書くが,これらを厳密に導出する第10章・第11章では,記号を簡潔にするためそれぞれ $v(\rr)$,$E_{xc}$ と書く.指すものは同じである.
次に,$E[\rho]$ を粒子数一定の条件のもとで最小化することから導かれるのがKohn–Sham方程式である.
密度は得られた軌道から
\begin{equation} \rho(\rr) = \sum_{i}^{\text{占有}} \abs{\phi_i(\rr)}^2 \label{eq:rho-ks} \end{equation}と再構成される.式 \eqref{eq:ks-eq} の右辺の有効ポテンシャルが $\rho$ に依存し,その $\rho$ は式 \eqref{eq:rho-ks} で軌道から作られるのだから,これは自分自身を含む方程式である.したがって適当な $\rho$ から出発して,$v_{\mathrm{eff}}$ を作り,方程式を解き,新しい $\rho$ を作り……という反復で解く.この手続きを自己無撞着(self-consistent field; SCF)計算といい,第10章と第14章で詳しく扱う.$\delta E_{xc}/\delta\rho$ という記号は汎関数微分であり,第4章で定義する新しい数学的道具である.
物理的意味:何が起きたのか
式 \eqref{eq:h-elec} と式 \eqref{eq:ks-eq} を見比べてほしい.前者は $3N$ 変数の連立した多体方程式,後者は3変数の1電子方程式である.電子間相互作用 $1/\abs{\rr_i-\rr_j}$ という「2つの座標を結ぶ項」が消え,代わりに1つの座標だけの関数 $v_{\mathrm{eff}}(\rr)$ が現れている.指数の壁は,こうして形式的には取り払われた.もちろんタダではない.代償は $E_{xc}[\rho]$ という未知の汎関数に押し込められており,そこには多体問題の困難がそっくり残っている.DFTが偉大なのは,困難を消したからではなく,困難を一箇所に隔離し,しかもその一箇所が小さい(全エネルギーの数%)ことを示した点にある.小さい項なら,粗い近似でも全体の精度はそれほど損なわれない.これがDFTが実用理論として成功した本質的な理由である.
ここで提示した式 \eqref{eq:e-rho},\eqref{eq:ks-eq},\eqref{eq:rho-ks} は,本書の残りの部分の目的地である.第2章から第8章までは,これらの式がどこから来るのかを理解するための準備——化学結合の描像,電子ガスの性質,密度で書いたエネルギーの原型,交換と相関の意味——を積み上げる道のりであり,第9章と第10章でこれらの式を厳密に導出する.以下の対応を頭に置いて読み進めるとよい.
- 「エネルギーを密度で書く」という発想の原型と汎関数微分の数学 → 第4章(Thomas–Fermi模型)
- $T_s$ にあたる相互作用のない系の運動エネルギー → 第3章(自由電子ガス)
- $E_{xc}$ の「交換」の部分の正体 → 第5章(Hartree–Fock法),第7章(一様電子ガス)
- $E_{xc}$ の「相関」の部分 → 第7章,第11章
- 式 \eqref{eq:e-rho} の厳密な根拠 → 第9章(Hohenberg–Kohnの定理)
- 式 \eqref{eq:ks-eq} の導出と $\varepsilon_i$ の意味 → 第10章(Kohn–Sham方程式)
1.8 まとめと本書の道筋
1.8.1 本章のまとめ
- 第一原理計算とは,経験パラメータを一切使わず,原子番号 $\{Z_I\}$・電子数 $N$・原子配置 $\{\bm{R}_I\}$ だけを入力として量子力学の基本法則から物性を予言する方法である.出力の中心は全エネルギー $E(\{\bm{R}_I\})$ であり,そこから安定構造・力・弾性定数・フォノン・反応障壁などが導かれる.
- 多電子系のハミルトニアンは,電子と核の運動エネルギー,電子–核引力,電子間反発,核間反発の5項からなる(式 \eqref{eq:h-si}).この形は水素分子から遷移金属結晶まで普遍であり,物質の個性はすべて解の多様性として現れる.
- Hartree原子単位系($\hbar=m_e=e=4\pi\varepsilon_0=1$)は,水素原子のSchrödinger方程式を無次元化する手続きから自然に導かれる.長さの単位はBohr半径 $a_0 = 4\pi\varepsilon_0\hbar^2/(m_ee^2) = 0.529177$ Å,エネルギーの単位はHartree $E_h = m_ee^4/(4\pi\varepsilon_0\hbar)^2 = 27.2114$ eV $=627.5$ kcal/mol である(表1.1).無次元化して消せない数——核質量比 $M_I/m_e\sim10^{3\text{–}5}$ と $c=1/\alpha=137$——が,それぞれBorn–Oppenheimer近似と非相対論近似の妥当性を支配する.
- Born–Oppenheimer近似は,Born–Huang展開 \eqref{eq:born-huang} という厳密な出発点から,$1/M_I$ に比例する非断熱結合項を落とすことで得られる(式 \eqref{eq:coupled} → \eqref{eq:bo-final}).非対角結合は電子状態間のエネルギーギャップに反比例するので(式 \eqref{eq:f-offdiag}),ギャップが閉じる状況——円錐交差,光化学,金属の電子–格子相互作用——では破綻する.近似の結果,核が運動する舞台としてのポテンシャルエネルギー面 $E_0(\bm{R})$ が確立する.
- 多電子波動関数は,電子の不可弁別性から交換位相が $\pm1$ に限られ(式 \eqref{eq:pm-one}),電子はフェルミオンなので反対称でなければならない(式 \eqref{eq:antisym2}).ここからPauliの排他律 \eqref{eq:pauli-zero} と交換ホールが直ちに従う.電子密度は $\rho(\rr)=N\sum_\sigma\int\dd x_2\cdots\dd x_N\abs{\Psi}^2$ で定義され,$\int\rho\,\dd^3r=N$ を満たす.
- 波動関数は $3N$ 次元配位空間の関数であり,必要な情報量は電子数に対して指数関数的に増える.1次元あたり10点という粗いグリッドでも水分子(10電子)で $10^{30}$ 個の数値が必要で,これは人類の全ストレージの1億倍にあたる.基底関数展開に切り替えても完全CIの次元は $\binom{M}{N}\sim 4^N$ と指数関数的に増える.これが指数の壁である.
- 壁を回避する戦略は,基本変数を取り替える(DFT・多体グリーン関数法),確率的にサンプリングする(量子モンテカルロ),波動関数の形を制限する(HF・結合クラスター)の3つに分類される(表1.3,図1.4).本書は $O(N^3)$ の計算量で数百〜数万原子を扱えるDFTを軸に据える.
- DFTは2本の柱で立つ.(1) 全エネルギーは電子密度の汎関数である(式 \eqref{eq:e-rho}).(2) 多体効果を交換相関エネルギー $E_{xc}[\rho]$ に集約し,一体問題(Kohn–Sham方程式 \eqref{eq:ks-eq})に還元する.この2式が本書の目的地である.
1.8.2 本書の道筋
本章で提示した式 \eqref{eq:e-rho} と \eqref{eq:ks-eq} に到達し,さらにそれを使いこなすまでの道筋を図1.5に示す.本書は6つの部からなり,理論の構築(第I部〜第III部),固体への適用(第IV部),数値実装(第V部),応用(第VI部)という順に進む.
各部の狙いを短く述べておく.
- 第I部(第1〜3章)多電子問題の基礎.本章で立てた問題を,2つの側面から掘り下げる.第2章では座標スケーリングによってビリアル定理 $2T+U=\sum_n\bm{R}_n\cdot\bm{F}_n$ を導き,化学結合が生じるときに運動エネルギーとポテンシャルエネルギーがどう変化しなければならないかという一般則を得る.第3章では相互作用のない電子の集団(自由電子ガス)を厳密に解き,Fermi球・状態密度・運動エネルギー密度といった,以後あらゆる場面で使う基本量を手に入れる.
- 第II部(第4〜8章)密度汎関数理論への道.「エネルギーを密度で書く」という発想が生まれ,鍛えられていく過程をたどる.第4章のThomas–Fermi模型はその最初の実例であり,同時に汎関数微分という数学的道具を導入する場でもある.第5章のHartree–Fock法は反対称性の帰結を完全に展開して交換の正体を暴き,第6章の第二量子化はそれを扱いやすい言語に翻訳する.第7章のジェリウムモデルで交換・相関エネルギーを密度の関数として具体的に求め,第8章では電子系が外場をどう遮蔽するかを線形応答理論で調べる.
- 第III部(第9〜11章)密度汎関数理論の核心.第9章でHohenberg–Kohnの定理を証明し,第10章でKohn–Sham方程式を変分的に導出する.第11章では未知の $E_{xc}[\rho]$ をどう近似するか——LDAがなぜうまくいくのか,GGAやハイブリッド汎関数は何を改善するのか——を論じる.
- 第IV部(第12〜13章)固体の電子状態.周期系に特有の構造(Blochの定理,Brillouinゾーン,バンド構造)と,局所的な原子配置が状態密度の形をどう決めるか(Green関数,モーメント,Friedelモデル)を扱う.
- 第V部(第14〜15章)数値実装.Kohn–Sham方程式を実際に計算機で解く方法.基底関数の選択,一般化固有値問題,自己無撞着ループの収束,そして内殻電子を消去して計算量を劇的に減らす擬ポテンシャル法を扱う.
- 第VI部(第16〜18章)応用.Hellmann–Feynman力と応力から構造最適化・第一原理分子動力学へ,非平衡Green関数による量子輸送へ,そして内殻励起・光学応答といった分光スペクトルの計算へと進む.
次章への橋渡し.本章では方程式を書き下し,それが直接には解けないことを確認した.しかし「解けない」ことと「何も言えない」ことは別である.実は,式 \eqref{eq:h-au} のハミルトニアンが「運動エネルギー項は座標の $-2$ 乗,Coulomb項は $-1$ 乗のスケーリングをする」という構造を持つ,というだけの情報から,厳密で強力な定理が1つ導ける.それがビリアル定理であり,そこから「なぜ原子が集まると結合するのか」という問いに対する一般的な答え——凝集の際には運動エネルギーは必ず増加し,ポテンシャルエネルギーは必ず減少する——が得られる.波動関数を1度も解かずに得られるこの結論を,次章で導く.
1.8.3 演習問題
演習1.1:原子単位の導出単位に慣れる
表1.1に載っていない導出単位を,自分で組み立てて数値を求めよ.
- 振動分光でよく使われる波数の単位 $1\ \mathrm{cm^{-1}}$ は何 $E_h$ か.逆に $1\,E_h$ は何 $\mathrm{cm^{-1}}$ か.
- 電場の原子単位 $E_h/(e\,a_0)$ をSI単位($\mathrm{V/m}$)で表せ.得られた値を,レーザーで達成される電場や誘電破壊の閾値($\sim10^{8}\ \mathrm{V/m}$)と比較して意味を述べよ.
- 圧力の原子単位 $E_h/a_0^3$ をPaおよびGPaで表せ.地球中心の圧力(約 $360$ GPa)は原子単位でいくらか.
ヒント:(1) 光子のエネルギーは $E=hc\tilde{\nu}$($\tilde{\nu}$ は波数).$h=2\pi\hbar$ に注意.(2) 電位の単位は $E_h/e = 27.2114$ V であるから,これを長さの単位 $a_0$ で割ればよい.(3) 圧力はエネルギー密度と同じ次元を持つ.$a_0^3 = 1.4818\times10^{-31}\ \mathrm{m^3}$ を使え.いずれの場合も,まず次元が合っていることを確認してから数値を入れること.
演習1.2:指数の壁を自分で見積もる
- 実用的な精度が出る $p=100$ のグリッドを使う場合,倍精度の波動関数を $10^{16}$ バイト(10 PB)以内に収められる電子数 $N$ の最大値を求めよ.
- 最も基底を切り詰めた場合($M=2N$ 個のスピン軌道)の完全CI次元 $\binom{2N}{N}$ が $10^{12}$ を超えるのは,$N$ がいくつからか.式 \eqref{eq:central-binom} の近似式を使い,常用対数を取って解け.
- 計算機の性能と記憶容量が2年で2倍になる(Moore則)と仮定する.$p=10$ のグリッド表現で扱える電子数が1個増えるのに何年かかるか.
ヒント:(1) 必要バイト数は $8\times10^{6N}$.両辺の常用対数を取れば1次不等式になる.(2) $\log_{10}\binom{2N}{N} \simeq N\log_{10}4 - \tfrac12\log_{10}(\pi N)$.$\log_{10}4=0.602$ を使い,まず第2項を無視した粗い解を求めてから代入して補正するとよい.(3) 電子1個の追加は必要量を $10^3$ 倍にする.$2^{t/2}=10^{3}$ を $t$ について解けばよい.
演習1.3:対角Born–Oppenheimer補正(DBOC)
§1.4.4で,断熱近似において実際に落とされるのは2階の対角項 $G^I_{ll}=\braket{\Phi_l|\nabla_I^2\Phi_l}$ だけであることを見た.この項の性質を調べる.電子波動関数 $\Phi_l$ は実関数に取れるものとする.
- 式 \eqref{eq:f-diag-zero} の $\braket{\Phi_l|\nabla_I\Phi_l}=0$ の両辺をさらに $\nabla_I$ で微分することにより, $$G^I_{ll} = -\braket{\nabla_I\Phi_l|\nabla_I\Phi_l}$$ を示せ.
- これを式 \eqref{eq:coupled} に代入し,対角項だけを残した核の方程式が,ポテンシャルエネルギー面を $$E_l(\bm{R}) \;\longrightarrow\; E_l(\bm{R}) + \sum_I \frac{1}{2M_I}\braket{\nabla_I\Phi_l|\nabla_I\Phi_l}$$ と置き換えたものになることを示せ.またこの補正項が常に非負であることを述べよ.
- この補正の大きさが $O(m_e/M_I)$ であることを,$\braket{\nabla_I\Phi_l|\nabla_I\Phi_l}$ が原子単位で $O(1)$ であることを根拠に議論せよ.H$_2$ とD$_2$(重水素分子)ではどちらで補正が大きいか.この違いは実験的にどのような形で現れうるか.
ヒント:(1) 積の微分則 $\nabla_I\braket{\Phi_l|\nabla_I\Phi_l} = \braket{\nabla_I\Phi_l|\nabla_I\Phi_l} + \braket{\Phi_l|\nabla_I^2\Phi_l}$ を使う.(2) 式 \eqref{eq:coupled} の結合項には全体に $-$ 符号が付いていることに注意.$\braket{f|f}\ge0$ は内積の性質である.(3) $M_{\mathrm{D}}\simeq2M_{\mathrm{H}}$ である.同位体置換で電子状態は変わらないが核質量だけが変わる,という点が鍵になる.
演習1.4:2電子の反対称波動関数と電子密度
規格直交な2つのスピン軌道 $\varphi_a(x)$,$\varphi_b(x)$($\braket{\varphi_a|\varphi_b}=\delta_{ab}$)から作った2電子波動関数
$$ \Psi(x_1,x_2) = \frac{1}{\sqrt{2}}\Big[\varphi_a(x_1)\varphi_b(x_2) - \varphi_b(x_1)\varphi_a(x_2)\Big] $$について,以下を確かめよ.
- $\Psi$ は式 \eqref{eq:antisym2} の反対称性を満たし,かつ規格化条件 \eqref{eq:psi-norm} を満たすこと(係数 $1/\sqrt2$ がその役割を果たしていること).
- 式 \eqref{eq:rho-def} の定義から,この状態の電子密度が $\rho(\rr) = \sum_\sigma\big[\abs{\varphi_a(\rr,\sigma)}^2 + \abs{\varphi_b(\rr,\sigma)}^2\big]$ となること.また $\int\rho\,\dd^3r = 2$ となること.
- $\varphi_a = \varphi_b$ と置くと $\Psi \equiv 0$ となること.これはPauliの排他律 \eqref{eq:pauli-zero} の具体的な現れである.
- (発展)2電子を同じ位置 $\rr$ に置き,スピンを両方 $\uparrow$ に取ると $\Psi=0$ になるが,一方を $\uparrow$,他方を $\downarrow$ に取ると一般には $\Psi\neq0$ であることを,$\varphi_a=\phi_a(\rr)\eta_\uparrow$,$\varphi_b=\phi_b(\rr)\eta_\downarrow$ の場合について具体的に確かめよ.交換ホールが同スピン電子の間にだけ生じる理由を説明せよ.
ヒント:(1) $\abs{\Psi}^2$ を展開すると4項出る.うち2項は $\abs{\varphi_a(x_1)}^2\abs{\varphi_b(x_2)}^2$ 型で積分すると1,残り2項は交差項で,直交性 $\braket{\varphi_a|\varphi_b}=0$ により消える.(2) 定義式で $x_2$ について積分する際に同じ直交性を使う.因子 $N=2$ を忘れないこと.(4) スピン関数の直交性 $\sum_\sigma\eta^*_\uparrow(\sigma)\eta_\downarrow(\sigma)=0$ を使う.
参考文献
- P. A. M. Dirac, "Quantum mechanics of many-electron systems", Proc. R. Soc. Lond. A 123, 714 (1929).
- E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, "CODATA recommended values of the fundamental physical constants: 2018", Rev. Mod. Phys. 93, 025010 (2021). 本章の物理定数の数値はこれに基づく.
- M. Born and J. R. Oppenheimer, "Zur Quantentheorie der Molekeln", Ann. Phys. (Leipzig) 84, 457 (1927).
- M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford University Press, 1954), Appendix VIII.
- W. Pauli, "The connection between spin and statistics", Phys. Rev. 58, 716 (1940).
- W. Kohn, "Nobel Lecture: Electronic structure of matter—wave functions and density functionals", Rev. Mod. Phys. 71, 1253 (1999). 「指数の壁」の議論はこの講演による.
- P. Hohenberg and W. Kohn, "Inhomogeneous electron gas", Phys. Rev. 136, B864 (1964).
- W. Kohn and L. J. Sham, "Self-consistent equations including exchange and correlation effects", Phys. Rev. 140, A1133 (1965).
- A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (McGraw-Hill, 1989).
- T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory (Wiley, 2000).
- W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, "Quantum Monte Carlo simulations of solids", Rev. Mod. Phys. 73, 33 (2001).
- L. Hedin, "New method for calculating the one-particle Green's function with application to the electron-gas problem", Phys. Rev. 139, A796 (1965).
- G. Onida, L. Reining, and A. Rubio, "Electronic excitations: density-functional versus many-body Green's-function approaches", Rev. Mod. Phys. 74, 601 (2002).
- R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004).