密度汎関数理論入門 — 目次 第III部 密度汎関数理論の核心 / 第10章

第10章コーン・シャム方程式

前章では，ホーヘンベルグ・コーンの定理とLevyの制約付き探索によって，基底状態エネルギーが電子密度 $\rho(\rr)$ の汎関数の最小化で厳密に得られることを確立した．しかし理論を実用にするには，普遍汎関数 $F[\rho]=T[\rho]+V_{ee}[\rho]$ の中身，とりわけ運動エネルギー汎関数 $T[\rho]$ の具体形が必要である．第4章で見たトーマス・フェルミ(TF)模型はこれを局所近似で置き換えたが，殻構造も化学結合も再現できなかった．本章では，この障害を回避するKohnとShamの独創的な戦略——「同じ密度を持つ相互作用しない仮想系」を導入し，運動エネルギーの大部分を軌道を使って厳密に評価する——を学び，密度汎関数理論の中心方程式であるコーン・シャム(KS)方程式を，変分計算を一切省略せずに導出する．これは本書全体のクライマックスであり，現代の第一原理計算のほぼすべてがこの方程式の上に立っている．未知の部分はすべて交換相関エネルギー $E_{xc}[\rho]$ という一つの項に集約され，その近似(LDA・GGAなど)は次章の主題となる．

この章で学ぶこと
  • TF模型の失敗の主因が運動エネルギーの局所近似にあること，そして「非相互作用の仮想系」で $T_s[\rho]$ を厳密に評価するというKohn–Shamの戦略
  • エネルギー汎関数の分解 $E[\rho]=T_s+\int v\rho+E_{\mathrm{H}}+E_{xc}$ と，交換相関エネルギー $E_{xc}$ の定義(=「残り全部」)
  • 軌道の規格直交拘束付き変分によるKS方程式 $\bigl[-\tfrac{1}{2}\nabla^2+v_{\mathrm{eff}}(\rr)\bigr]\phi_i=\varepsilon_i\phi_i$ の完全な導出と，ラグランジュ乗数行列の対角化(正準形)
  • 密度に関するオイラー方程式からの別導出と，$v_{\mathrm{eff}}$ の同定，KS系の存在仮定($v_s$-表示可能性)
  • 全エネルギーの表式 $E=\sum_i\varepsilon_i-E_{\mathrm{H}}-\int v_{xc}\rho+E_{xc}$(二重数え補正)——固有値の和は全エネルギーではない
  • 自己無撞着場(SCF)ループの構造と，スピン密度汎関数法への拡張
  • KS固有値の意味:Janakの定理 $\partial E/\partial n_i=\varepsilon_i$ の完全証明，$\varepsilon_{\mathrm{HOMO}}=-I$，$\Delta$SCF，光電子分光との比較とバンドギャップ問題

本章でもハートリー原子単位系($\hbar=m_e=e=4\pi\varepsilon_0=1$，第1章参照)を用いる．エネルギーの単位はハートリー(Ha)で，$1\ \mathrm{Ha}=27.2114\ \mathrm{eV}$ である．電子数は $N$，外部ポテンシャル(原子核がつくるクーロンポテンシャル)は $v(\rr)$ と書く．

10.1 コーン・シャムの戦略

まず，なぜ「軌道」を持ち込む必要があるのかを，TF模型の失敗の分析からはっきりさせる．そのうえで，Kohn–Shamの仮想系のアイデアとエネルギー汎関数の分解を導入する．

10.1.1 トーマス・フェルミ模型はなぜ失敗したか

第4章で導入したTF模型は，運動エネルギー汎関数を一様電子ガスの結果(第3章のフェルミ球)で局所近似した:

$$ \begin{equation} T^{\mathrm{TF}}[\rho] = C_F\int \dd^3r\,\rho^{5/3}(\rr),\qquad C_F=\frac{3}{10}\left(3\pi^2\right)^{2/3}. \label{eq:10-tf-kinetic} \end{equation} $$

この近似の帰結は深刻だった．原子の殻構造(動径密度分布のピーク構造)が消え，分子の化学結合は一切安定に存在できず(Tellerの定理)，負イオンも記述できない．原因を定量的に見るために，アルゴン原子の運動エネルギーを方法別に比較してみる．

表10.1 アルゴン原子($N=18$)の運動エネルギーの比較(単位:Ha)．KS-LDAは本章で導入するコーン・シャム法に局所密度近似(第11章)を組み合わせたもの．
方法運動エネルギー (Ha)ハートリー・フォック(HF)との差
ハートリー・フォック526.82—
トーマス・フェルミ489.95約 $-37$ Ha(約 $-7\%$)
コーン・シャム(LDA)525.95約 $-0.9$ Ha(約 $-0.2\%$)

TF近似の誤差は約 37 Ha，率にして 7% である．「たった7%」と思ってはいけない．ビリアル定理(第2章)により，クーロン系の平衡状態では $E=-T$ が成り立つから，運動エネルギーは全エネルギーと同じ巨大なスケール(Arで約 527 Ha $\approx$ 14300 eV)を持つ．一方，化学結合のエネルギーはせいぜい 0.1〜0.3 Ha(数 eV)である．運動エネルギーの 7%(約 1000 eV!)の誤差は，化学結合のエネルギースケールの数千倍であり，これでは結合の有無すら論じられない．つまり，

物理的意味:誤差のスケール比較

全エネルギーの大部分を占める運動エネルギーを，化学精度(〜0.001 Ha)どころか 1 Ha の精度でも近似できないことが，TF系理論の致命傷である．逆に言えば，運動エネルギーさえ高精度に評価できれば，密度汎関数理論は実用理論になり得る．表10.1のKS-LDAの行(誤差 0.2%)は，本章で述べる方法がまさにそれを達成したことを示している．

それでは，なぜTF近似はこれほど運動エネルギーを誤るのか．運動エネルギーは波動関数の「曲がり具合」($|\nabla\phi|^2$)で決まる量であり，殻構造・節(ノード)・軌道の直交性といった，密度 $\rho(\rr)$ の局所値だけからは読み取れない情報に強く依存するからである．ここに「軌道を経由して運動エネルギーを評価する」という発想の必然性がある．

10.1.2 相互作用しない仮想系(コーン・シャム補助系)

W. KohnとL. J. Shamは1965年，次のような大胆な迂回路を提案した(文献[2])．

定義:コーン・シャム補助系(KS系)

相互作用する $N$ 電子系の基底状態密度 $\rho(\rr)$ が与えられたとき，互いに全く相互作用せず，共通の局所ポテンシャル $v_s(\rr)$ の中を運動する $N$ 個の電子からなる仮想系であって，その基底状態密度が同じ $\rho(\rr)$ に一致するものを，コーン・シャム補助系と呼ぶ．相互作用がないので，この仮想系の基底状態は1電子方程式

$$ \begin{equation} \left[-\frac{1}{2}\nabla^2+v_s(\rr)\right]\phi_i(\rr)=\varepsilon_i\phi_i(\rr) \label{eq:10-aux-schrodinger} \end{equation} $$

の解である規格直交軌道 $\{\phi_i\}$ を，パウリの排他律に従ってエネルギーの低い順に $N$ 個占有した単一スレーター行列式(第5章)で与えられ，密度と運動エネルギーは

$$ \begin{equation} \rho(\rr)=\sum_{i=1}^{N}\abs{\phi_i(\rr)}^2, \qquad T_s=\sum_{i=1}^{N}\int \dd^3r\,\phi_i^*(\rr)\left(-\frac{1}{2}\nabla^2\right)\phi_i(\rr) \label{eq:10-ts-def} \end{equation} $$

となる．$T_s$ を非相互作用運動エネルギーと呼ぶ．和は占有軌道についてとる．ここで $i$ はスピンも込めた量子数であり($N$ 本のスピン軌道)，スピン縮退した閉殻系では空間軌道 $N/2$ 本に占有数2を与える書き方も同値である(10.6節)．

要点は次の通りである．$T_s$ は式 \eqref{eq:10-ts-def} のとおり軌道の汎関数として書かれているが，軌道は(仮想系の1電子方程式を通じて)密度から決まるので，$T_s$ は密度の陰関数的な汎関数 $T_s[\rho]$ でもある．より正確には，第9章のLevyの制約付き探索と同じ流儀で

$$ \begin{equation} T_s[\rho]=\min_{\Phi\to\rho}\braket{\Phi|\hat{T}|\Phi} \label{eq:10-ts-levy} \end{equation} $$

と定義できる．ここで $\Phi$ は密度 $\rho$ を与える単一スレーター行列式全体を走る．つまり「$\rho$ を再現する行列式のうち運動エネルギーが最小のもの」の運動エネルギーが $T_s[\rho]$ である．$T_s[\rho]$ の $\rho$ による陽な表式は不要であることに注意してほしい．実際の計算では軌道表式 \eqref{eq:10-ts-def} を使って厳密に評価すればよい．これがTF近似 \eqref{eq:10-tf-kinetic} との決定的な違いである．

$T_s$ は真の運動エネルギー $T$ と厳密には一致しない($T_s \le T$．相互作用する真の波動関数は単一行列式より複雑だからである)．しかし後で見るように，その差 $T-T_s$ は相関エネルギーと同程度(原子で 1 Ha 以下)の小さな量であり，これを未知項に繰り込んでも致命傷にならない．「大きな量は厳密に，小さな量だけを近似する」——これがKS戦略の核心である．

現実の相互作用系 電子間クーロン反発 1/|r−r'| あり 外部ポテンシャル v(r)・多体波動関数 Ψ コーン・シャム補助系(仮想) 電子間相互作用なし(独立粒子) 有効ポテンシャル veff(r)・軌道 {φi} 写像 同一の基底状態密度 ρ(r) を共有
図10.1 コーン・シャム法の概念図．現実の相互作用系(左)と同じ基底状態密度 $\rho(\rr)$ を持つように，有効ポテンシャル $v_{\mathrm{eff}}(\rr)$ 中の独立粒子系(右)を構成する．運動エネルギーの大部分は右の系で厳密に評価される．

10.1.3 エネルギー汎関数の分解と交換相関エネルギーの定義

これから，第9章で導入した全エネルギー汎関数を，KS系を基準にして厳密に書き換える．近似はまだ一切入らない．

第9章の結果から，外部ポテンシャル $v(\rr)$ 中の $N$ 電子系の基底状態エネルギーは

$$ \begin{equation} E[\rho]=F[\rho]+\int \dd^3r\,v(\rr)\rho(\rr), \qquad F[\rho]=T[\rho]+V_{ee}[\rho] \label{eq:10-hk-energy} \end{equation} $$

の最小化で得られる．$T[\rho]$ は真の運動エネルギー，$V_{ee}[\rho]$ は真の電子間相互作用エネルギーである．ここで $F[\rho]$ に，KS系で厳密に計算できる二つの量を「足して引く」．一つは非相互作用運動エネルギー $T_s[\rho]$，もう一つは古典的クーロン反発(ハートリー)エネルギー

$$ \begin{equation} E_{\mathrm{H}}[\rho]=\frac{1}{2}\iint \dd^3r\,\dd^3r'\, \frac{\rho(\rr)\rho(\rr')}{\abs{\rr-\rr'}} \label{eq:10-hartree-def} \end{equation} $$

である．すなわち恒等変形

\begin{align} F[\rho]&=T[\rho]+V_{ee}[\rho]\notag\\ &=T_s[\rho]+E_{\mathrm{H}}[\rho] +\underbrace{\bigl(T[\rho]-T_s[\rho]\bigr)+\bigl(V_{ee}[\rho]-E_{\mathrm{H}}[\rho]\bigr)}_{\displaystyle \equiv\,E_{xc}[\rho]} \label{eq:10-f-decomp} \end{align}

を行う．1行目から2行目へは $T_s+E_{\mathrm{H}}$ を加えて同じものを引いただけであり，等号は厳密である．こうして定義される

$$ \begin{equation} E_{xc}[\rho]\equiv\bigl(T[\rho]-T_s[\rho]\bigr)+\bigl(V_{ee}[\rho]-E_{\mathrm{H}}[\rho]\bigr) \label{eq:10-exc-def} \end{equation} $$

を交換相関エネルギー(exchange-correlation energy)と呼ぶ．これを使うと全エネルギー汎関数は

$$ \begin{equation} E[\rho]=T_s[\rho]+\int \dd^3r\,v(\rr)\rho(\rr)+E_{\mathrm{H}}[\rho]+E_{xc}[\rho] \label{eq:10-ks-decomp} \end{equation} $$

と分解される．各項の意味を整理しよう．

物理的意味:KS分解は「誤差の隔離」である

式 \eqref{eq:10-ks-decomp} の右辺4項のうち，初めの3項は厳密に計算でき，しかも大きさの順でいえば圧倒的にこの3項が支配的である．たとえばアルゴン原子では全エネルギー約 $-527$ Ha に対し，交換エネルギーは約 $-30$ Ha，相関エネルギーは約 $-0.7$ Ha，$T-T_s$ も 1 Ha 以下にすぎない．TF模型が運動エネルギー(数百 Ha)そのものを近似したのに対し，KS法では近似の対象を全体の数%以下の $E_{xc}$ だけに隔離する．$E_{xc}$ の近似が多少粗くても(例えば10%の誤差でも)，全エネルギーへの影響は数 Ha 程度に抑えられ，しかもエネルギー差を取る際に大部分が相殺する．近似理論の設計として，これ以上ないほど賢い分業である．

注意:ハートリー項は自己相互作用を含む

$E_{\mathrm{H}}[\rho]$ は密度分布 $\rho$ が「自分自身」と相互作用するエネルギーであり，電子1個の系($N=1$，例:水素原子)でもゼロにならない．1個の電子が自分自身とクーロン反発することは物理的にあり得ないから，これは自己相互作用(self-interaction)という偽の寄与である．厳密な $E_{xc}$ は定義 \eqref{eq:10-exc-def} により($N=1$ では $V_{ee}=0,\ T=T_s$ なので)$E_{xc}=-E_{\mathrm{H}}$ となってこれを完全に打ち消す．しかし第11章で見るLDAやGGAなどの近似汎関数は打ち消しが不完全で，自己相互作用誤差と呼ばれる系統誤差が残る．バンドギャップの過小評価(10.7節)や局在軌道の記述の失敗の根源の一つである．

以上で本章の土台が整った．得られたものをまとめると:全エネルギーは \eqref{eq:10-ks-decomp} の形に厳密に書け，未知なのは $E_{xc}[\rho]$ だけである．次節では，この汎関数を軌道について変分し，KS方程式を導出する．

10.2 KS方程式の導出

この節が本章の，そして本書全体のクライマックスである．エネルギー汎関数 \eqref{eq:10-ks-decomp} を，軌道の規格直交条件を拘束としてラグランジュ未定乗数法で変分し，KS方程式を導く．使う数学は第4章で導入した汎関数微分と，これから数学ノートで整備する「軌道汎関数の連鎖律」だけである．1行も飛ばさずに進む．

10.2.1 変分問題の設定

式 \eqref{eq:10-ks-decomp} の $E[\rho]$ において，$T_s$ は軌道 $\{\phi_i\}$ の陽な汎関数，残りの3項は $\rho$ の汎関数であり，$\rho$ 自身は \eqref{eq:10-ts-def} により軌道から組み立てられる．したがって $E$ 全体を軌道の汎関数 $E[\{\phi_i,\phi_i^*\}]$ と見なし，軌道について最小化すればよい．ただし軌道は勝手な関数ではいけない．KS補助系はフェルミ粒子系であり，その基底状態は $N$ 本の規格直交軌道を占有した単一スレーター行列式だから，拘束条件

$$ \begin{equation} \int \dd^3r\,\phi_i^*(\rr)\phi_j(\rr)=\delta_{ij} \qquad (i,j=1,\dots,N) \label{eq:10-orthonorm} \end{equation} $$

を課す必要がある．対角成分($i=j$)は規格化条件で，これにより $\int\rho\,\dd^3r=N$(電子数保存)が自動的に保証される．非対角成分($i\neq j$)は直交条件で，パウリの排他律により「$N$ 個の電子が互いに異なる状態を占める」ことを表現している．直交性を課さないと，全電子が同じ最低軌道に落ち込む解(ボソン的な解)が最小値になってしまう．

拘束条件 \eqref{eq:10-orthonorm} は $i,j$ の組ごとに1本，合計 $N^2$ 本ある(複素共役の関係 $\braket{\varphi_j|\varphi_i}=\braket{\varphi_i|\varphi_j}^*$ により，独立な複素方程式は対角 $N$ 本と上三角 $N(N-1)/2$ 本の計 $N(N+1)/2$ 本にまで減る．対角成分はもともと実数なので，実数の条件として数え直せば $N+N(N-1)=N^2$ 本となり，第5章5.6.1項の数え方と一致する)．ラグランジュ未定乗数法では，拘束1本につき乗数を1個導入するから，乗数は行列 $\varepsilon_{ij}$ になる．そこでラグランジュ汎関数を

$$ \begin{equation} \Omega[\{\phi_i,\phi_i^*\}] = E[\{\phi_i,\phi_i^*\}] -\sum_{k=1}^{N}\sum_{l=1}^{N}\varepsilon_{lk} \left(\int \dd^3r\,\phi_k^*(\rr)\phi_l(\rr)-\delta_{kl}\right) \label{eq:10-lagrangian} \end{equation} $$

と構成する(乗数の添字を $\varepsilon_{lk}$ と書いたのは，後で行列表記がきれいになるようにするための便宜である)．この $\Omega$ を，拘束なしの独立変数としての $\phi_i^*$ について停留にする:

$$ \begin{equation} \frac{\delta \Omega}{\delta \phi_i^*(\rr)}=0 \qquad (i=1,\dots,N). \label{eq:10-stationary} \end{equation} $$

数学ノート:複素関数による変分——なぜ $\phi$ と $\phi^*$ を独立に扱えるか

軌道 $\phi_i(\rr)$ は複素関数である．複素数 $z=x+iy$ を変数に持つ実数値関数 $f(z,z^*)$ の停留条件を考えよう．本来の独立変数は実部 $x$ と虚部 $y$ の2つであり，停留条件は $\partial f/\partial x=0$ かつ $\partial f/\partial y=0$ である．ここで形式的な微分(Wirtinger微分)

$$ \begin{equation} \frac{\partial f}{\partial z}\equiv\frac{1}{2}\left(\frac{\partial f}{\partial x}-i\frac{\partial f}{\partial y}\right), \qquad \frac{\partial f}{\partial z^*}\equiv\frac{1}{2}\left(\frac{\partial f}{\partial x}+i\frac{\partial f}{\partial y}\right) \label{eq:10-wirtinger} \end{equation} $$

を定義すると，2条件 $\partial f/\partial x=\partial f/\partial y=0$ は2条件 $\partial f/\partial z=\partial f/\partial z^*=0$ と完全に同値である(上の2式を足し引きすれば互いに移り合う)．つまり「$z$ と $z^*$ をあたかも独立変数のように扱って別々に微分する」流儀は，実部・虚部で微分する正攻法の言い換えにすぎない．さらに $f$ が実数値なら $\partial f/\partial z=(\partial f/\partial z^*)^*$ が成り立つので，片方の条件 $\partial f/\partial z^*=0$ だけ課せば十分である(もう片方はその複素共役として自動的に成立する)．

汎関数の場合も全く同様で，実数値汎関数 $\Omega[\phi,\phi^*]$ に対して $\delta\Omega/\delta\phi_i^*=0$ を課せば，$\delta\Omega/\delta\phi_i=0$ は自動的に満たされる．以下ではこの流儀で $\phi_i^*$ についての変分だけを計算する．

数学ノート:軌道汎関数の連鎖律

汎関数微分そのものは第4章で導入した($\delta A=\int \dd^3r\,\frac{\delta A}{\delta \rho(\rr)}\delta\rho(\rr)$ を満たす $\frac{\delta A}{\delta\rho(\rr)}$ が $A$ の汎関数微分である)．ここで必要になるのは，「$\rho$ の汎関数 $A[\rho]$ を，$\rho$ を組み立てている軌道 $\phi_i^*$ で微分する」ための連鎖律である．

密度は $\rho(\rr')=\sum_j \phi_j^*(\rr')\phi_j(\rr')$ だから，$\phi_i^*$ だけを $\phi_i^*+\delta\phi_i^*$ と変化させたときの密度の一次変化は

$$ \begin{equation} \delta\rho(\rr')=\delta\phi_i^*(\rr')\,\phi_i(\rr') \label{eq:10-drho} \end{equation} $$

である($\phi_j$ 側は動かしていないので，和のうち $j=i$ の項の $\phi_i^*$ 因子だけが変化する)．これを $A$ の一次変化に代入すると

\begin{align} \delta A &=\int \dd^3r'\,\frac{\delta A}{\delta\rho(\rr')}\,\delta\rho(\rr') \notag\\ &=\int \dd^3r'\,\frac{\delta A}{\delta\rho(\rr')}\,\phi_i(\rr')\,\delta\phi_i^*(\rr'). \label{eq:10-chain-1} \end{align}

1行目は汎関数微分の定義，2行目は \eqref{eq:10-drho} の代入である．一方，$A$ を $\phi_i^*$ の汎関数と見たときの微分の定義は $\delta A=\int \dd^3r'\,\frac{\delta A}{\delta\phi_i^*(\rr')}\delta\phi_i^*(\rr')$ だから，両者を見比べて

$$ \begin{equation} \frac{\delta A}{\delta \phi_i^*(\rr)} =\frac{\delta A}{\delta \rho(\rr)}\,\phi_i(\rr) \label{eq:10-chain-rule} \end{equation} $$

を得る．これが軌道汎関数の連鎖律である．同じことを $\delta\rho(\rr')/\delta\phi_i^*(\rr)=\phi_i(\rr)\,\delta(\rr-\rr')$ と書き，デルタ関数を挟んだ積分

$$ \begin{equation} \frac{\delta A}{\delta \phi_i^*(\rr)} =\int \dd^3r'\,\frac{\delta A}{\delta \rho(\rr')}\,\frac{\delta\rho(\rr')}{\delta\phi_i^*(\rr)} =\int \dd^3r'\,\frac{\delta A}{\delta \rho(\rr')}\,\phi_i(\rr')\,\delta(\rr-\rr') \label{eq:10-chain-delta} \end{equation} $$

として理解してもよい(最後の等号でデルタ関数の積分を実行すれば \eqref{eq:10-chain-rule} に戻る)．

10.2.2 各項の汎関数微分

これから，$\Omega$ を構成する5つの項——$T_s$，$\int v\rho$，$E_{\mathrm{H}}$，$E_{xc}$，拘束項——の $\phi_i^*(\rr)$ による汎関数微分を1つずつ計算する．

導出:$\delta T_s/\delta\phi_i^*(\rr)$

$T_s$ の定義 \eqref{eq:10-ts-def} を積分変数 $\rr'$ で書くと

$$ \begin{equation} T_s=\sum_{j=1}^{N}\int \dd^3r'\,\phi_j^*(\rr')\left(-\frac{1}{2}\nabla'^2\right)\phi_j(\rr'). \label{eq:10-ts-rewrite} \end{equation} $$

$T_s$ は $\phi_j^*$ について線形なので，$\phi_i^*$ による汎関数微分は「$\phi_i^*(\rr')$ をデルタ関数 $\delta(\rr'-\rr)$ で置き換える」操作に等しい(第4章:線形汎関数 $A=\int f g$ の微分は $\delta A/\delta f(\rr)=g(\rr)$)．和のうち生き残るのは $j=i$ の項だけであり，

\begin{align} \frac{\delta T_s}{\delta \phi_i^*(\rr)} &=\int \dd^3r'\,\delta(\rr'-\rr)\left(-\frac{1}{2}\nabla'^2\right)\phi_i(\rr') \notag\\ &=-\frac{1}{2}\nabla^2\phi_i(\rr). \label{eq:10-dts} \end{align}

1行目はデルタ関数への置き換え，2行目はデルタ関数の性質 $\int \dd^3r'\,\delta(\rr'-\rr)g(\rr')=g(\rr)$ による積分の実行である．

導出:$\delta\bigl(\int v\rho\bigr)/\delta\phi_i^*(\rr)$

外部ポテンシャル項 $A_{\mathrm{ext}}=\int \dd^3r'\,v(\rr')\rho(\rr')$ は $\rho$ の線形汎関数なので，$\rho$ による微分は

$$ \begin{equation} \frac{\delta A_{\mathrm{ext}}}{\delta\rho(\rr)}=v(\rr). \label{eq:10-dext-rho} \end{equation} $$

連鎖律 \eqref{eq:10-chain-rule} を適用して

$$ \begin{equation} \frac{\delta A_{\mathrm{ext}}}{\delta \phi_i^*(\rr)} =\frac{\delta A_{\mathrm{ext}}}{\delta\rho(\rr)}\,\phi_i(\rr) =v(\rr)\,\phi_i(\rr). \label{eq:10-dext} \end{equation} $$

導出:$\delta E_{\mathrm{H}}/\delta\phi_i^*(\rr)$

まず $E_{\mathrm{H}}$ の $\rho$ による微分を求める．定義 \eqref{eq:10-hartree-def} で $\rho\to\rho+\delta\rho$ と置き，$\delta\rho$ の一次の項だけを拾うと($\delta\rho$ の二次の項は変分では捨てる)

\begin{align} \delta E_{\mathrm{H}} &=\frac{1}{2}\iint \dd^3r_1\,\dd^3r_2\, \frac{\delta\rho(\rr_1)\,\rho(\rr_2)}{\abs{\rr_1-\rr_2}} +\frac{1}{2}\iint \dd^3r_1\,\dd^3r_2\, \frac{\rho(\rr_1)\,\delta\rho(\rr_2)}{\abs{\rr_1-\rr_2}} \label{eq:10-dh-expand} \\ &=2\times\frac{1}{2}\iint \dd^3r_1\,\dd^3r_2\, \frac{\delta\rho(\rr_1)\,\rho(\rr_2)}{\abs{\rr_1-\rr_2}} \label{eq:10-dh-symm} \\ &=\int \dd^3r_1\,\delta\rho(\rr_1) \underbrace{\int \dd^3r_2\,\frac{\rho(\rr_2)}{\abs{\rr_1-\rr_2}}}_{\displaystyle \equiv\,v_{\mathrm{H}}(\rr_1)}. \label{eq:10-dh-var} \end{align}

\eqref{eq:10-dh-expand} は積の変分(ライプニッツ則)である．\eqref{eq:10-dh-symm} へは，第2項の積分変数の名前を入れ替えた($\rr_1\leftrightarrow\rr_2$．定積分の値は積分変数の名前によらないから許される):クーロン核 $1/\abs{\rr_1-\rr_2}$ は $\rr_1$ と $\rr_2$ の入れ替えに対して対称なので，名前を替えた第2項は第1項と完全に同じ形になり，因子2が $1/2$ を打ち消す．\eqref{eq:10-dh-var} では $\rr_2$ 積分を内側にまとめた．汎関数微分の定義と見比べて

$$ \begin{equation} \frac{\delta E_{\mathrm{H}}}{\delta\rho(\rr)} =v_{\mathrm{H}}(\rr) =\int \dd^3r'\,\frac{\rho(\rr')}{\abs{\rr-\rr'}}. \label{eq:10-vh-def} \end{equation} $$

$v_{\mathrm{H}}(\rr)$ をハートリーポテンシャルと呼ぶ．電荷分布 $\rho$ がつくる古典的静電ポテンシャルそのものである．連鎖律 \eqref{eq:10-chain-rule} により

$$ \begin{equation} \frac{\delta E_{\mathrm{H}}}{\delta \phi_i^*(\rr)} =v_{\mathrm{H}}(\rr)\,\phi_i(\rr). \label{eq:10-dh} \end{equation} $$

導出:$\delta E_{xc}/\delta\phi_i^*(\rr)$

$E_{xc}[\rho]$ の陽な形は知らないが，それでも困らない．$\rho$ の汎関数であることさえ認めれば，その汎関数微分を

$$ \begin{equation} v_{xc}(\rr)\equiv\frac{\delta E_{xc}[\rho]}{\delta \rho(\rr)} \label{eq:10-vxc-def} \end{equation} $$

と記号として定義し，連鎖律 \eqref{eq:10-chain-rule} を適用すればよい:

$$ \begin{equation} \frac{\delta E_{xc}}{\delta \phi_i^*(\rr)} =v_{xc}(\rr)\,\phi_i(\rr). \label{eq:10-dxc} \end{equation} $$

$v_{xc}(\rr)$ を交換相関ポテンシャルと呼ぶ．$E_{xc}$ の近似形(第11章)を決めれば $v_{xc}$ は具体的に計算できる関数になる．たとえばLDAのように $E_{xc}=\int \varepsilon_{xc}(\rho(\rr))\rho(\rr)\dd^3r$ という局所形なら，$v_{xc}(\rr)=\dfrac{\dd}{\dd\rho}\bigl[\rho\,\varepsilon_{xc}(\rho)\bigr]\Big|_{\rho=\rho(\rr)}$ という「ただの関数」になる．

導出:拘束項の汎関数微分

拘束項 $C\equiv\sum_{k,l}\varepsilon_{lk}\bigl(\int \dd^3r'\,\phi_k^*(\rr')\phi_l(\rr')-\delta_{kl}\bigr)$ は $\phi_k^*$ について線形である．$\phi_i^*$ による微分では和のうち $k=i$ の項だけが生き残り，デルタ関数への置き換えと積分の実行により

\begin{align} \frac{\delta C}{\delta\phi_i^*(\rr)} &=\sum_{l}\varepsilon_{li}\int \dd^3r'\,\delta(\rr'-\rr)\,\phi_l(\rr') \notag\\ &=\sum_{l=1}^{N}\varepsilon_{li}\,\phi_l(\rr). \label{eq:10-dconstraint} \end{align}

10.2.3 停留条件と非正準形のKS方程式

以上の部品 \eqref{eq:10-dts}, \eqref{eq:10-dext}, \eqref{eq:10-dh}, \eqref{eq:10-dxc}, \eqref{eq:10-dconstraint} を停留条件 \eqref{eq:10-stationary} に代入して組み立てる:

\begin{align} 0=\frac{\delta \Omega}{\delta\phi_i^*(\rr)} &=-\frac{1}{2}\nabla^2\phi_i(\rr) +v(\rr)\phi_i(\rr) +v_{\mathrm{H}}(\rr)\phi_i(\rr) +v_{xc}(\rr)\phi_i(\rr) -\sum_{l=1}^{N}\varepsilon_{li}\,\phi_l(\rr). \label{eq:10-stationary-full} \end{align}

3つのポテンシャルをまとめて有効ポテンシャル

$$ \begin{equation} v_{\mathrm{eff}}(\rr)\equiv v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr) \label{eq:10-veff-def} \end{equation} $$

を定義し，1電子ハミルトニアン(KSハミルトニアン)を $\hat{h}_{\mathrm{KS}}\equiv-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}(\rr)$ と書けば，\eqref{eq:10-stationary-full} は

$$ \begin{equation} \hat{h}_{\mathrm{KS}}\,\phi_i(\rr)=\sum_{l=1}^{N}\varepsilon_{li}\,\phi_l(\rr) \qquad (i=1,\dots,N) \label{eq:10-ks-noncanonical} \end{equation} $$

となる．右辺に乗数行列が現れており，まだ固有値方程式の形をしていない．これを非正準形のKS方程式と呼ぶ．ここから対角化によって見慣れた固有値方程式(正準形)に持ち込む．

10.2.4 乗数行列のエルミート性とユニタリ対角化(正準形)

まず乗数行列 $\varepsilon_{li}$ がエルミート行列であることを示す．式 \eqref{eq:10-ks-noncanonical} の両辺に $\phi_m^*(\rr)$ を掛けて全空間で積分すると

\begin{align} \int \dd^3r\,\phi_m^*\,\hat{h}_{\mathrm{KS}}\,\phi_i &=\sum_{l=1}^{N}\varepsilon_{li}\int \dd^3r\,\phi_m^*\,\phi_l \notag\\ &=\sum_{l=1}^{N}\varepsilon_{li}\,\delta_{ml} \notag\\ &=\varepsilon_{mi}. \label{eq:10-eps-matrix-element} \end{align}

1行目から2行目へは規格直交条件 \eqref{eq:10-orthonorm}(停留点では拘束が満たされている)を使い，2行目から3行目へはクロネッカーデルタで和をつぶした．つまり乗数行列は $\hat{h}_{\mathrm{KS}}$ の行列要素 $\varepsilon_{mi}=\braket{\phi_m|\hat{h}_{\mathrm{KS}}|\phi_i}$ に等しい．ここで $\hat{h}_{\mathrm{KS}}$ はエルミート演算子である．実際，$v_{\mathrm{eff}}(\rr)$ は実関数(実密度から作られる実ポテンシャルの和)だから掛け算演算子としてエルミートであり，$-\frac{1}{2}\nabla^2$ は2回の部分積分でエルミート性が示せる(境界項は，束縛状態の軌道が無限遠で十分速くゼロに近づくため消える．第2章参照)．したがって

$$ \begin{equation} \varepsilon_{mi} =\braket{\phi_m|\hat{h}_{\mathrm{KS}}|\phi_i} =\braket{\hat{h}_{\mathrm{KS}}\phi_m|\phi_i} =\braket{\phi_i|\hat{h}_{\mathrm{KS}}|\phi_m}^* =\varepsilon_{im}^*, \label{eq:10-eps-hermite} \end{equation} $$

すなわち $\varepsilon$ はエルミート行列である(2番目の等号がエルミート性，3番目は内積の共役対称性)．

次に，占有軌道のユニタリ変換に対する不変性を確認する．$N\times N$ ユニタリ行列 $U$($U^\dagger U=\bm{1}$)で占有軌道を混ぜた新しい軌道

$$ \begin{equation} \phi_k'(\rr)=\sum_{j=1}^{N}\phi_j(\rr)\,U_{jk} \label{eq:10-unitary-rot} \end{equation} $$

を考える．密度は

\begin{align} \rho'(\rr) &=\sum_{k=1}^{N}\phi_k'^*(\rr)\,\phi_k'(\rr) \notag\\ &=\sum_{k=1}^{N}\sum_{j=1}^{N}\sum_{j'=1}^{N}U_{jk}^*\,\phi_j^*(\rr)\,\phi_{j'}(\rr)\,U_{j'k} \notag\\ &=\sum_{j,j'=1}^{N}\phi_j^*(\rr)\,\phi_{j'}(\rr)\underbrace{\sum_{k=1}^{N}U_{j'k}\,U_{jk}^*}_{=(UU^\dagger)_{j'j}=\delta_{j'j}} \notag\\ &=\sum_{j=1}^{N}\abs{\phi_j(\rr)}^2=\rho(\rr). \label{eq:10-rho-invariant} \end{align}

1行目は定義，2行目は \eqref{eq:10-unitary-rot} の代入，3行目は $k$ についての和を先に実行できるよう並べ替え，$UU^\dagger=\bm{1}$ を使って4行目に至る．全く同じ計算(微分演算子 $-\frac{1}{2}\nabla^2$ を挟むだけ)で $T_s$ も不変であることが分かる．つまり密度・$T_s$・したがって全エネルギーは，占有軌道をユニタリ変換しても変わらない．物理的に意味があるのは個々の軌道ではなく，占有軌道が張る $N$ 次元部分空間である．

この自由度を使って乗数行列を対角化する．$\varepsilon$ はエルミートなので，ユニタリ行列 $U$ と実対角行列 $\varepsilon_{\mathrm{d}}=\mathrm{diag}(\varepsilon_1,\dots,\varepsilon_N)$ により

$$ \begin{equation} \varepsilon=U\,\varepsilon_{\mathrm{d}}\,U^\dagger, \qquad\text{すなわち}\qquad \varepsilon\,U=U\,\varepsilon_{\mathrm{d}} \label{eq:10-eps-diag} \end{equation} $$

と対角化できる(エルミート行列のスペクトル定理．固有値 $\varepsilon_k$ は実数)．行ベクトル表記 $\bm{\phi}=(\phi_1,\dots,\phi_N)$ を使うと，非正準形 \eqref{eq:10-ks-noncanonical} は $\hat{h}_{\mathrm{KS}}\,\bm{\phi}=\bm{\phi}\,\varepsilon$ とまとめて書ける．両辺に右から $U$ を掛けると

\begin{align} \hat{h}_{\mathrm{KS}}\,\bm{\phi}\,U &=\bm{\phi}\,\varepsilon\,U \notag\\ &=\bm{\phi}\,U\,\varepsilon_{\mathrm{d}}. \label{eq:10-diag-step} \end{align}

1行目は $\hat{h}_{\mathrm{KS}}$ が(関数にはたらく演算子であって軌道ラベルには作用しないので)数行列 $U$ と可換であることを使って左辺を書き換えたもの，2行目は \eqref{eq:10-eps-diag} の代入である．そこで $\bm{\phi}'=\bm{\phi}U$，すなわち \eqref{eq:10-unitary-rot} の新軌道を使えば，\eqref{eq:10-diag-step} は成分ごとに

$$ \begin{equation} \hat{h}_{\mathrm{KS}}\,\phi_k'(\rr)=\varepsilon_k\,\phi_k'(\rr) \label{eq:10-canonical-step} \end{equation} $$

となる．しかも \eqref{eq:10-rho-invariant} で見たとおり，この変換で $\rho$ は不変だから，$\rho$ から作られる $v_{\mathrm{H}}$, $v_{xc}$, したがって $\hat{h}_{\mathrm{KS}}$ 自身も変わらない．プライムを落として書き直せば，最終結果に到達する．

$$ \begin{equation} \left[-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}(\rr)\right]\phi_i(\rr)=\varepsilon_i\,\phi_i(\rr), \qquad v_{\mathrm{eff}}(\rr)=v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr), \qquad v_{xc}(\rr)=\frac{\delta E_{xc}[\rho]}{\delta\rho(\rr)} \label{eq:10-ks-equation} \end{equation} $$

これがコーン・シャム方程式(の正準形)である．得られたものを言葉でまとめよう:全エネルギー \eqref{eq:10-ks-decomp} を最小にする軌道は，有効ポテンシャル $v_{\mathrm{eff}}$ 中の1電子シュレーディンガー方程式の固有関数として(適切なユニタリ変換のもとで)選べる．$N$ 電子の多体問題が，形式上は「1電子問題を $N$ 本解く」問題に厳密に置き換わった．

物理的意味:KS方程式の構造

注意:どの $N$ 個の軌道を占有するか

変分の停留条件が与えるのは「占有軌道はKS方程式の固有関数から選べる」ことまでである．エネルギー最小の解では，固有値の低い順に $N$ 本を占有する(アウフバウ原理)．もし占有すべき最低の空軌道と最高の占有軌道の固有値が一致(縮退)する場合は，分数占有を許す拡張(10.7節のJanakの枠組み)が自然な処方を与える．

10.3 密度のオイラー方程式との整合(別導出)

前節では軌道について変分してKS方程式を得た．この節では視点を変え，密度そのものについて変分する第9章流のオイラー方程式から，同じ $v_{\mathrm{eff}}$ にたどり着くことを示す．二つの導出を突き合わせることで，KS法の論理構造——特に，どこに存在仮定が潜んでいるか——がはっきり見える．

10.3.1 相互作用系のオイラー方程式

第9章で見たとおり，基底状態密度は電子数拘束 $\int\dd^3r\,\rho(\rr)=N$ のもとで $E[\rho]=F[\rho]+\int v\rho$ を最小化する．拘束のラグランジュ乗数を $\mu$ として，汎関数

$$ \begin{equation} \Omega_\rho[\rho]=F[\rho]+\int \dd^3r\,v(\rr)\rho(\rr)-\mu\left(\int \dd^3r\,\rho(\rr)-N\right) \label{eq:10-omega-rho} \end{equation} $$

の停留条件 $\delta\Omega_\rho/\delta\rho(\rr)=0$ をとる．各項の微分は，$F$ の微分が定義により $\delta F/\delta\rho$，外部ポテンシャル項が線形汎関数なので $v(\rr)$，拘束項が同じく線形なので $\mu$ であり，

$$ \begin{equation} \frac{\delta F[\rho]}{\delta \rho(\rr)}+v(\rr)=\mu \label{eq:10-euler-int} \end{equation} $$

を得る．これが相互作用系のオイラー方程式であり，乗数 $\mu$ は化学ポテンシャル(電子を1個加えるのに要するエネルギーの連続版)の意味を持つ(第4章のTF理論の $\mu_{\mathrm{TF}}$ と同じ論理である)．ここで $F$ をKS分解 \eqref{eq:10-f-decomp} に従って $F=T_s+E_{\mathrm{H}}+E_{xc}$ と書き直すと，\eqref{eq:10-vh-def} と \eqref{eq:10-vxc-def} により

$$ \begin{equation} \frac{\delta T_s}{\delta \rho(\rr)}+v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr)=\mu. \label{eq:10-euler-int2} \end{equation} $$

10.3.2 KS系のオイラー方程式と $v_{\mathrm{eff}}$ の同定

今度は，局所ポテンシャル $v_s(\rr)$ の中の相互作用しない $N$ 電子系を考える．この系のエネルギー汎関数は運動エネルギーと外場の項しかない:

$$ \begin{equation} E_s[\rho]=T_s[\rho]+\int \dd^3r\,v_s(\rr)\rho(\rr). \label{eq:10-es-def} \end{equation} $$

同じ電子数拘束のもとで変分すると(乗数を $\mu_s$ とする)，オイラー方程式は

$$ \begin{equation} \frac{\delta T_s}{\delta \rho(\rr)}+v_s(\rr)=\mu_s. \label{eq:10-euler-ks} \end{equation} $$

さて，KS法の要求は「この非相互作用系の基底状態密度が，相互作用系の基底状態密度 $\rho(\rr)$ に一致すること」であった．同じ $\rho$ が \eqref{eq:10-euler-int2} と \eqref{eq:10-euler-ks} の両方を満たすなら，両式の $\delta T_s/\delta\rho$ は同じ密度で評価した同じ量だから，2式を辺々引き算して

\begin{align} 0&=\Bigl[v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr)-\mu\Bigr]-\Bigl[v_s(\rr)-\mu_s\Bigr] \notag\\ \Longrightarrow\quad v_s(\rr)&=v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr)+(\mu_s-\mu). \label{eq:10-veff-identify} \end{align}

定数 $\mu_s-\mu$ はポテンシャルの原点の付け替えにすぎず，固有値全体と化学ポテンシャルを同じだけずらすだけで密度や全エネルギーに影響しない．よって定数を除いて

$$ \begin{equation} v_s(\rr)=v_{\mathrm{eff}}(\rr)=v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr) \label{eq:10-veff-identify2} \end{equation} $$

が同定される．10.2節で軌道変分から得た有効ポテンシャル \eqref{eq:10-veff-def} と完全に一致した．すなわち，KS方程式を自己無撞着に解いて得た密度は，密度に関する変分条件 $\delta E/\delta\rho=\mu$(電子数一定のもとでの $\delta E=0$)を自動的に満たす．

物理的意味:この別導出から見えること

10.3.3 隠れた仮定:$v_s$-表示可能性(KS法はAnsatzである)

この導出は，KS法が立脚する仮定を最も鮮明に見せてくれる．それは「相互作用系の基底状態密度 $\rho(\rr)$ を基底状態密度として持つような非相互作用系(局所ポテンシャル $v_s$)が存在する」という仮定である．これを密度の非相互作用 $v$-表示可能性($v_s$-表示可能性)と呼ぶ．第9章で相互作用系の $v$-表示可能性を論じたのと同型の問題が，非相互作用系側でも生じているわけである．

任意の(物理的に妥当な)密度に対してこの仮定が成り立つかどうかは，一般には証明されていない．したがって厳密に言えば，KS法は定理ではなくAnsatz(仮設)である．ただし実用上この仮定が破綻して問題になることは稀であり，また高精度計算で得た密度から $v_{\mathrm{eff}}$ を逆算する「逆コーン・シャム法」の研究により，典型的な原子・分子・固体の密度については対応する $v_s$ が数値的に構成できることが確かめられている．本書では以後，この仮定を認めて進む．

10.4 全エネルギーの表式

KS方程式を解いて固有値 $\varepsilon_i$ と軌道が得られたとして，全エネルギー $E$ をどう計算するか．初学者が必ず一度は誤解する点——「固有値を足せば全エネルギーになるのではないか」——を，ここで完全に清算する．これから (i) 固有値の和が何に等しいかを証明し，(ii) それを使って全エネルギーの正しい表式を導く．

10.4.1 固有値の和は何か

導出:$\sum_i\varepsilon_i=T_s+\int v_{\mathrm{eff}}\,\rho$ の証明

KS方程式 \eqref{eq:10-ks-equation} の両辺に左から $\phi_i^*(\rr)$ を掛ける:

$$ \begin{equation} \phi_i^*(\rr)\left(-\frac{1}{2}\nabla^2\right)\phi_i(\rr) +v_{\mathrm{eff}}(\rr)\,\phi_i^*(\rr)\,\phi_i(\rr) =\varepsilon_i\,\phi_i^*(\rr)\,\phi_i(\rr). \label{eq:10-sum-step1} \end{equation} $$

全空間で積分すると，右辺は規格化条件 $\int\abs{\phi_i}^2\dd^3r=1$ により $\varepsilon_i$ となる:

$$ \begin{equation} \int \dd^3r\,\phi_i^*\left(-\frac{1}{2}\nabla^2\right)\phi_i +\int \dd^3r\,v_{\mathrm{eff}}(\rr)\abs{\phi_i(\rr)}^2 =\varepsilon_i. \label{eq:10-sum-step2} \end{equation} $$

占有軌道 $i=1,\dots,N$ について辺々足し合わせる:

\begin{align} \sum_{i=1}^{N}\varepsilon_i &=\sum_{i=1}^{N}\int \dd^3r\,\phi_i^*\left(-\frac{1}{2}\nabla^2\right)\phi_i +\int \dd^3r\,v_{\mathrm{eff}}(\rr)\sum_{i=1}^{N}\abs{\phi_i(\rr)}^2 \label{eq:10-sum-step3} \\ &=T_s+\int \dd^3r\,v_{\mathrm{eff}}(\rr)\,\rho(\rr). \label{eq:10-eigsum} \end{align}

\eqref{eq:10-sum-step3} では和と積分の順序を交換し(有限和なので常に許される)，\eqref{eq:10-eigsum} では第1項に $T_s$ の定義 \eqref{eq:10-ts-def}，第2項に密度の定義 \eqref{eq:10-ts-def} を使った．

固有値の和は「運動エネルギー+有効ポテンシャル中のポテンシャルエネルギー」である．ところが全エネルギー \eqref{eq:10-ks-decomp} に入っているのは外場 $v$ とのエネルギーと $E_{\mathrm{H}}+E_{xc}$ であって，$\int v_{\mathrm{eff}}\rho$ ではない．差を正しく処理しよう．

10.4.2 二重数え補正と全エネルギー

まず補題として，ハートリーポテンシャルと密度の積分が $E_{\mathrm{H}}$ のちょうど2倍になることを示す:

\begin{align} \int \dd^3r\,v_{\mathrm{H}}(\rr)\,\rho(\rr) &=\int \dd^3r\left[\int \dd^3r'\,\frac{\rho(\rr')}{\abs{\rr-\rr'}}\right]\rho(\rr) \label{eq:10-vhrho-1} \\ &=\iint \dd^3r\,\dd^3r'\,\frac{\rho(\rr)\,\rho(\rr')}{\abs{\rr-\rr'}} \label{eq:10-vhrho-2} \\ &=2E_{\mathrm{H}}. \label{eq:10-vhrho} \end{align}

\eqref{eq:10-vhrho-1} は $v_{\mathrm{H}}$ の定義 \eqref{eq:10-vh-def} の代入，\eqref{eq:10-vhrho-2} は積分順序をまとめただけ，\eqref{eq:10-vhrho} は $E_{\mathrm{H}}$ の定義 \eqref{eq:10-hartree-def}(係数 $1/2$ 付き)との比較である．

これで準備が整った．全エネルギー \eqref{eq:10-ks-decomp} から出発し，$T_s$ を固有値和で置き換える:

\begin{align} E&=T_s+\int \dd^3r\,v\rho+E_{\mathrm{H}}+E_{xc} \label{eq:10-etot-1} \\ &=\left[\sum_{i=1}^{N}\varepsilon_i-\int \dd^3r\,v_{\mathrm{eff}}\,\rho\right] +\int \dd^3r\,v\rho+E_{\mathrm{H}}+E_{xc} \label{eq:10-etot-2} \\ &=\sum_{i=1}^{N}\varepsilon_i -\int \dd^3r\,\bigl(v+v_{\mathrm{H}}+v_{xc}\bigr)\rho +\int \dd^3r\,v\rho+E_{\mathrm{H}}+E_{xc} \label{eq:10-etot-3} \\ &=\sum_{i=1}^{N}\varepsilon_i -\int \dd^3r\,v_{\mathrm{H}}\,\rho -\int \dd^3r\,v_{xc}\,\rho+E_{\mathrm{H}}+E_{xc} \label{eq:10-etot-4} \\ &=\sum_{i=1}^{N}\varepsilon_i -2E_{\mathrm{H}}+E_{\mathrm{H}} -\int \dd^3r\,v_{xc}\,\rho+E_{xc}. \label{eq:10-etot-5} \end{align}

\eqref{eq:10-etot-2} は \eqref{eq:10-eigsum} を $T_s$ について解いて代入，\eqref{eq:10-etot-3} は $v_{\mathrm{eff}}$ の定義 \eqref{eq:10-veff-def} の展開，\eqref{eq:10-etot-4} では $\int v\rho$ が相殺し，\eqref{eq:10-etot-5} で補題 \eqref{eq:10-vhrho} を使った．整理して:

$$ \begin{equation} E=\sum_{i=1}^{N}\varepsilon_i -E_{\mathrm{H}}[\rho] -\int \dd^3r\,v_{xc}(\rr)\,\rho(\rr) +E_{xc}[\rho] \label{eq:10-etot-final} \end{equation} $$

これがKS法の全エネルギーの標準的表式である(二重数え補正付きの表式と呼ばれる)．得られた教訓を強調しておく:

物理的意味:なぜ固有値の和は全エネルギーではないのか

各KS電子は，全電子(自分も含む)がつくる平均場 $v_{\mathrm{H}}$ の中でのエネルギー $\varepsilon_i$ を持つ．電子 $i$ と電子 $j$ の静電反発は，$\varepsilon_i$ の中に($j$ がつくる場として)1回，$\varepsilon_j$ の中に($i$ がつくる場として)もう1回，合計2回数えられてしまう．だから固有値を単純に足すと電子間反発が二重に入り，$-E_{\mathrm{H}}$ で1回分引き戻さなければならない．交換相関についても同様に，ポテンシャル経由で数えた $\int v_{xc}\rho$ を引いて，正しいエネルギー $E_{xc}$ を足し直す．$\sum_i\varepsilon_i$(バンドエネルギーと呼ばれる)はあくまで補助量である——ただし，エネルギー差や摂動応答には有用な情報を含む(第13章のフリーデル模型や第16章の力の計算で再会する)．

例:2つの極限での検算

(1) 相互作用のない極限．もし電子間相互作用が存在しなければ $E_{\mathrm{H}}=E_{xc}=0$，$v_{xc}=0$ であり，\eqref{eq:10-etot-final} は $E=\sum_i\varepsilon_i$ に帰着する．独立粒子系では固有値の和が全エネルギーそのものであり，表式は正しい極限を持つ．

(2) 水素原子($N=1$)．厳密汎関数では，10.1節の注意で述べたとおり $E_{xc}[\rho]=-E_{\mathrm{H}}[\rho]$ であり，さらに $v_{xc}(\rr)=-v_{\mathrm{H}}(\rr)$(汎関数微分をとればよい)．このとき $v_{\mathrm{eff}}=v+v_{\mathrm{H}}-v_{\mathrm{H}}=v=-1/r$ となり，KS方程式は水素原子のシュレーディンガー方程式そのものになって $\varepsilon_1=-0.5$ Ha を与える．全エネルギーは

$$ E=\varepsilon_1-E_{\mathrm{H}}-\int v_{xc}\rho+E_{xc} =\varepsilon_1-E_{\mathrm{H}}+2E_{\mathrm{H}}-E_{\mathrm{H}} =\varepsilon_1=-0.5\ \mathrm{Ha} $$

と正しい値になる(途中で $\int v_{xc}\rho=-\int v_{\mathrm{H}}\rho=-2E_{\mathrm{H}}$ を使った)．一方，LDAなどの近似汎関数では $E_{xc}\neq-E_{\mathrm{H}}$ なので，水素原子ですら全エネルギーに誤差($+0.02$ Ha 程度)が出る．自己相互作用誤差の最も簡単な実例である．

10.5 自己無撞着場(SCF)ループ

KS方程式 \eqref{eq:10-ks-equation} は普通の固有値問題に見えるが，ポテンシャル $v_{\mathrm{eff}}$ が解(の密度)に依存するという点で本質的に非線形である．この節では，この非線形問題を反復で解く標準手続き——自己無撞着場(self-consistent field, SCF)ループ——の構造を理解する．数値的な詳細(混合法の各種，収束加速)は第14章に譲り，ここでは論理の骨格を押さえる．

10.5.1 不動点問題としてのKS方程式

手続きを抽象的に書くと次のようになる．密度 $\rho_{\mathrm{in}}$ を1つ与えると，$v_{\mathrm{eff}}[\rho_{\mathrm{in}}]$ が決まり，KS方程式を解いて軌道が得られ，そこから新しい密度 $\rho_{\mathrm{out}}$ が組み上がる．この一連の操作を写像

$$ \begin{equation} \rho_{\mathrm{out}}=\mathcal{K}[\rho_{\mathrm{in}}] \label{eq:10-ks-map} \end{equation} $$

と書こう．求めたい解は，この写像の不動点

$$ \begin{equation} \rho^{\star}=\mathcal{K}[\rho^{\star}] \label{eq:10-fixed-point} \end{equation} $$

である．入力密度と出力密度が一致すること($\rho_{\mathrm{in}}=\rho_{\mathrm{out}}$)を自己無撞着(self-consistency)といい，このとき「軌道が作る密度」と「ポテンシャルを作った密度」が矛盾なく噛み合う．10.3節で示したとおり，自己無撞着に達した解こそが変分条件 $\delta E/\delta\rho=\mu$ を満たす基底状態密度である．逆に言えば，収束していない途中段階の密度・エネルギーは変分的な意味を持たないので，収束判定を甘くしてはならない．

10.5.2 ループの構成要素

実際のループは次の要素からなる(図10.2)．

  1. 初期密度の用意:通常は孤立原子の密度を重ね合わせた $\rho^{(0)}(\rr)=\sum_I \rho_{\mathrm{atom}}^{I}(\rr-\bm{R}_I)$ を使う．分子や固体の密度は原子密度の重ね合わせからそれほど大きくずれないため，良い初期値になる．
  2. 有効ポテンシャルの構成:$v_{\mathrm{H}}$ は定義 \eqref{eq:10-vh-def} の積分を直接実行してもよいが，実装ではポアソン方程式 $$ \begin{equation} \nabla^2 v_{\mathrm{H}}(\rr)=-4\pi\rho(\rr) \label{eq:10-poisson} \end{equation} $$ を解くことが多い(微分方程式のほうが高速に解ける場合が多い)．この式は，恒等式 $\nabla^2\frac{1}{\abs{\rr-\rr'}}=-4\pi\delta(\rr-\rr')$(クーロン核が3次元ラプラス方程式のグリーン関数であること．原点以外で $1/r$ が調和関数であることと，原点を囲む微小球でガウスの定理を使うことから従う)を \eqref{eq:10-vh-def} の積分の中で使えば直ちに得られる．$v_{xc}$ は採用した近似汎関数から局所的に計算する(第11章)．
  3. KS方程式の対角化:基底関数(平面波・局在基底など．第14・15章)でハミルトニアン行列を作り，固有値問題を解いて $\{\varepsilon_i,\phi_i\}$ を得る．計算コストの大半はここに集中する．
  4. 新密度の構成と収束判定:$\rho_{\mathrm{out}}(\rr)=\sum_i^{\mathrm{occ}}\abs{\phi_i(\rr)}^2$ を作り，例えば残差ノルム $$ \begin{equation} \Delta=\int \dd^3r\,\bigl[\rho_{\mathrm{out}}(\rr)-\rho_{\mathrm{in}}(\rr)\bigr]^2 \label{eq:10-residual} \end{equation} $$ が閾値(実用上 $10^{-6}$〜$10^{-10}$ 程度，単位系と実装による)を下回ったか，および全エネルギーの変化が十分小さいかで収束を判定する．
  5. 密度混合:収束していなければ次の入力密度を作る．最も単純なのは線形混合 $$ \begin{equation} \rho_{\mathrm{in}}^{(n+1)}=(1-\alpha)\,\rho_{\mathrm{in}}^{(n)}+\alpha\,\rho_{\mathrm{out}}^{(n)}, \qquad 0 \lt \alpha \le 1 \label{eq:10-mixing} \end{equation} $$ である．

なぜ $\rho_{\mathrm{out}}$ をそのまま次の入力にせず，わざわざ混ぜるのか．単純代入($\alpha=1$)は振動・発散することが多いからである．直観はこうだ:ある場所で入力密度が少なすぎると，そこの有効ポテンシャルは(電子の遮蔽が足りず)深くなりすぎ，出力密度は逆に多すぎになる．次の反復では逆向きに行きすぎる．系が大きいほど，また金属のように電荷が動きやすい系ほど，この「行きすぎ」の振り子は大きくなり(長波長の電荷振動，いわゆるチャージ・スロッシング)，反復は収束しない．混合率 $\alpha$ を小さくして一歩の修正を控えめにすれば振り子は減衰する——ただし収束は遅くなる．この安定性と速度のトレードオフを賢く解決する方法(Pulay混合・Broyden法・Kerkerプレコンディショナなど)は第14章で詳述する．

初期密度 ρ⁽⁰⁾(原子密度の重ね合わせ) 有効ポテンシャルの構成 ∇²vH = −4πρin, vxc = δExc/δρ KS方程式を対角化 [−½∇² + veff] φi = εi φi 新しい密度 ρout = Σᵢ |φi|² 収束判定 ∫(ρout−ρin)² d³r < 閾値 ? No 密度混合 (1−α)ρin+αρout 新しい ρin Yes 自己無撞着解(ρin = ρout = ρ*) 全エネルギー E = Σεi − EH − ∫vxc ρ + Exc 密度・力・各種物性の解析へ
図10.2 自己無撞着場(SCF)ループのフローチャート．入力密度から有効ポテンシャルを作り，KS方程式を解いて出力密度を得る．入力と出力が一致(自己無撞着)するまで，混合した密度で反復する．収束後の全エネルギーは二重数え補正付きの式 \eqref{eq:10-etot-final} で評価する．

例:SCF反復の典型的な振る舞い

小さな分子や単純な半導体では，線形混合($\alpha\sim0.3$)でも10〜30回程度の反復でエネルギーが $10^{-8}$ Ha レベルまで収束するのが典型である．一方，フェルミ準位近傍に状態が密集する磁性金属や，細長い金属スラブでは，単純混合は容易に発散し，Kerker型の前処理や準ニュートン法(第14章)が必須になる．「SCFが回らない」ときにパラメータを闇雲にいじるのではなく，どの物理(状態密度・長波長遮蔽・スピン自由度)が振動を起こしているかを考える習慣をつけたい．

10.6 スピン密度汎関数法

ここまでの定式化は，スピンについて暗黙に「上向きと下向きが同数で同じ空間分布を持つ」(スピン無偏極)と仮定してきた．しかし磁性体や奇数電子系ではこの仮定は成り立たない．この節では，KS法をスピン分極系に拡張したスピン密度汎関数法(spin-density functional theory)の形式を与える．この拡張はvon BarthとHedinによって整備された(文献[6])．

10.6.1 スピン密度を独立変数にとる

スピン $\sigma=\uparrow,\downarrow$ ごとの電子密度(スピン密度)

$$ \begin{equation} \rho^{\sigma}(\rr)=\sum_{i=1}^{N_\sigma}\abs{\phi_{i\sigma}(\rr)}^2, \qquad \rho(\rr)=\rho^{\uparrow}(\rr)+\rho^{\downarrow}(\rr), \qquad N_\uparrow+N_\downarrow=N \label{eq:10-spin-density} \end{equation} $$

を独立変数の組として，全エネルギー汎関数を

$$ \begin{equation} E[\rho^{\uparrow},\rho^{\downarrow}] =T_s[\rho^{\uparrow},\rho^{\downarrow}] +\int \dd^3r\,v(\rr)\,\rho(\rr) +E_{\mathrm{H}}[\rho] +E_{xc}[\rho^{\uparrow},\rho^{\downarrow}] \label{eq:10-lsda-energy} \end{equation} $$

と書く．ここで重要な非対称性に注意する．外部ポテンシャル項とハートリー項は全密度 $\rho$ だけに依存する(クーロン力はスピンを見ない)．一方，$T_s$ と $E_{xc}$ はスピンごとの内訳に依存する．$T_s$ は各スピンチャネルの軌道から

$$ \begin{equation} T_s=\sum_{\sigma=\uparrow,\downarrow}\sum_{i=1}^{N_\sigma} \int \dd^3r\,\phi_{i\sigma}^*(\rr)\left(-\frac{1}{2}\nabla^2\right)\phi_{i\sigma}(\rr) \label{eq:10-ts-spin} \end{equation} $$

と計算され，$E_{xc}$ がスピン内訳に依存するのは，交換がパウリの排他律により同種スピン間でだけ働くからである(第5章・第7章)．スピン分極の度合いは

$$ \begin{equation} \zeta(\rr)=\frac{\rho^{\uparrow}(\rr)-\rho^{\downarrow}(\rr)}{\rho(\rr)}, \qquad -1\le\zeta\le 1 \label{eq:10-zeta} \end{equation} $$

で測る($\zeta=0$:無偏極，$\zeta=\pm1$:完全偏極)．磁化密度は $m(\rr)=\rho^{\uparrow}(\rr)-\rho^{\downarrow}(\rr)$ であり，全磁気モーメントは $M=\int m\,\dd^3r=N_\uparrow-N_\downarrow$(ボーア磁子単位)となる．

10.6.2 スピン依存KS方程式

導出は10.2節の反復である．軌道 $\phi_{i\sigma}$ の規格直交拘束(直交性は同一スピン内でのみ課せばよい．異なるスピンの軌道はスピン関数の直交性で自動的に直交する)のもとで \eqref{eq:10-lsda-energy} を $\phi_{i\sigma}^*$ について変分すると，10.2節と全く同じ手順(各項の汎関数微分→乗数行列の対角化)により，スピンチャネルごとのKS方程式

$$ \begin{equation} \left[-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}^{\sigma}(\rr)\right]\phi_{i\sigma}(\rr) =\varepsilon_{i\sigma}\,\phi_{i\sigma}(\rr), \qquad v_{\mathrm{eff}}^{\sigma}(\rr)=v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}^{\sigma}(\rr), \qquad v_{xc}^{\sigma}(\rr)=\frac{\delta E_{xc}[\rho^{\uparrow},\rho^{\downarrow}]}{\delta \rho^{\sigma}(\rr)} \label{eq:10-spin-ks} \end{equation} $$

が得られる．新しく計算が必要になるのは $E_{xc}$ の微分だけである．チェックすべき点を確認しておこう:$v_{\mathrm{H}}$ は全密度から作られ両スピンに共通，$v_{xc}^{\sigma}$ はスピンごとに異なる．この $v_{xc}^{\uparrow}\neq v_{xc}^{\downarrow}$ という差こそが，交換相互作用による有効磁場(交換分裂)の起源である．占有は両スピンチャネルを通した共通のフェルミ準位で決める:全 $\{\varepsilon_{i\sigma}\}$ を低い順に並べ，下から $N$ 個を占有する．$N_\uparrow$ と $N_\downarrow$ は入力ではなく，SCFの結果として決まる(これが磁気モーメントの第一原理予測を可能にする)．

注意:無偏極の場合との整合と，いつスピン分極計算が必須か

閉殻系・非磁性系では自己無撞着解が $\rho^{\uparrow}=\rho^{\downarrow}=\rho/2$ を満たし，\eqref{eq:10-spin-ks} の2本の方程式は同一になって，スピン無偏極のKS方程式(空間軌道に占有数2)に帰着する．逆に次の場合はスピン分極形式が必須である:

また，$E_{xc}[\rho^{\uparrow},\rho^{\downarrow}]$ の具体的近似(局所スピン密度近似LSDA，スピン依存GGA)は第11章で構成する．磁場が空間的に回転する非共線磁性への拡張(2成分スピノル軌道と $2\times2$ 密度行列)もあるが，本書の範囲を超えるので言及にとどめる．

10.7 KS固有値の意味

10.2節の導出を振り返ると，KS固有値 $\varepsilon_i$ は規格直交拘束のラグランジュ乗数として方程式に入ってきた量である．乗数は変分問題を解くための道具であって，それ自体が観測量である保証はどこにもない．にもかかわらず，実際の第一原理計算では $\varepsilon_i$ を並べてバンド構造を描き，光電子分光の実験スペクトルと重ね合わせるのが日常である．この慣行はどこまで正当化されるのか．本節では (i) Janakの定理 $\partial E/\partial n_i=\varepsilon_i$，(ii) 厳密汎関数における $\varepsilon_{\mathrm{HOMO}}=-I$，(iii) $\Delta$SCF法とSlaterの遷移状態近似，(iv) 実験との比較とバンドギャップ問題，の順に，証明を一切省略せずに調べる．

10.7.1 原則:KS軌道と固有値は補助量である

まず原則を確認する．ここまでに証明したことは次の2点だけである．

KS法が保証するのはこの2つ，すなわち $\rho$ と $E$ だけである．軌道 $\phi_i$ は「密度を再現するための足場」であって現実の電子の波動関数ではないし，固有値 $\varepsilon_i$ も同様に，それ単独では現実の系の励起エネルギーやイオン化エネルギーを表さない．第5章のハートリー・フォック(HF)法にはKoopmansの定理——$\varepsilon_i^{\mathrm{HF}}$ が「軌道緩和を無視したときの $i$ 番目の電子のイオン化エネルギーの符号を変えたもの」に等しい——があったが，KS法にはそれと同型の定理は一般には存在しない．$\hat{h}_{\mathrm{KS}}$ はHFのフォック演算子と違い，「$N-1$ 電子が作る場」ではなく「$N$ 電子が作る場」を含んでいるからである．

注意:KS固有値について言ってはいけないこと

それでも $\varepsilon_i$ には明確な意味がある．ただしそれは「電子のエネルギー」ではなく，全エネルギーの占有数に関する微係数という形をとる．これがJanakの定理である．

10.7.2 分数占有数への拡張

Janakの定理を述べるには，「軌道 $i$ の占有数 $n_i$ を 0 と 1 の間で連続的に変えたときのエネルギー変化」を考えなければならない．そこでまず，KS形式を分数占有数へ拡張する．密度と非相互作用運動エネルギーを

$$ \begin{equation} \rho(\rr)=\sum_{i} n_i\abs{\phi_i(\rr)}^2, \qquad T_s=\sum_{i} n_i\int \dd^3r\,\phi_i^*(\rr)\left(-\frac{1}{2}\nabla^2\right)\phi_i(\rr), \qquad 0\le n_i\le 1,\quad \sum_i n_i=\nu \label{eq:10-frac-density} \end{equation} $$

と定義し直す．和は(空軌道も含めて)すべての軌道について走らせるが，$n_i=0$ の軌道は寄与しない．全電子数 $\nu=\sum_i n_i$ は整数でなくてもよいものとする．上限 $n_i\le1$ はパウリの排他律であり(スピンを含めた軌道1本に電子1個まで)，これを外すと全電子が最低軌道に落ちてしまう．全エネルギーは形の上では前と同じく

$$ \begin{equation} E\bigl[\{n_i\},\{\phi_i\}\bigr] =T_s+\int \dd^3r\,v(\rr)\rho(\rr)+E_{\mathrm{H}}[\rho]+E_{xc}[\rho] \label{eq:10-frac-energy} \end{equation} $$

と書かれる．$n_i$ は $T_s$ の中と，$\rho$ を通して残りの3項の中に，二重に入り込んでいることに注意しておく．

数学ノート:分数占有数は何を意味するか

「電子 0.5 個」は一見すると非物理的だが，次の3つの文脈で明確な意味を持つ．

(1) 開いた系(アンサンブル)．電子浴(リザーバー)と電子をやりとりできる系を考えると，粒子数は保存量ではなく，系の状態は粒子数の異なる状態の統計混合になる．粒子数 $M$ と $M+1$ の基底状態を重み $1-\omega$ と $\omega$ で混ぜた混合状態(密度演算子 $\hat{\Gamma}=(1-\omega)\ket{\Psi_M}\bra{\Psi_M}+\omega\ket{\Psi_{M+1}}\bra{\Psi_{M+1}}$)の平均粒子数は $M+\omega$ であり，平均密度は $\rho=(1-\omega)\rho_M+\omega\rho_{M+1}$ となる．分数電子数はこの平均値のことである．この定式化をアンサンブルDFTと呼ぶ(文献[5])．

(2) 有限温度．金属の実際の計算では，フェルミ準位近傍の軌道にフェルミ・ディラック分布 $n_i=[e^{(\varepsilon_i-\mu)/k_BT}+1]^{-1}$ を与える「スメアリング」を行う(第14章)．このとき $n_i$ は自動的に分数になる．

(3) 形式的な補間．いま必要なのは，$n_i$ を連続変数と見なしたときの $E$ の微係数だけである．式 \eqref{eq:10-frac-density}–\eqref{eq:10-frac-energy} はどんな $\{n_i\}\in[0,1]$ に対しても意味を持つ数式なので，微分は数学的に何の問題もなく定義できる．定理の証明にはこの立場だけで十分である．

分数占有の場合のKS方程式を確認しておこう．10.2節の連鎖律 \eqref{eq:10-chain-rule} は，いまや密度が $\rho(\rr')=\sum_j n_j\phi_j^*(\rr')\phi_j(\rr')$ なので，$\delta\rho(\rr')/\delta\phi_i^*(\rr)=n_i\phi_i(\rr)\delta(\rr-\rr')$ と修正され，

$$ \begin{equation} \frac{\delta A}{\delta \phi_i^*(\rr)}=n_i\,\frac{\delta A}{\delta\rho(\rr)}\,\phi_i(\rr) \label{eq:10-frac-chain} \end{equation} $$

となる．同様に $T_s$ の微分も \eqref{eq:10-dts} に $n_i$ が掛かる．したがって \eqref{eq:10-stationary-full} に対応する停留条件は全体に $n_i$ が掛かった形

$$ \begin{equation} \frac{\delta E}{\delta\phi_i^*(\rr)} =n_i\left[-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}(\rr)\right]\phi_i(\rr) =n_i\,\hat{h}_{\mathrm{KS}}\,\phi_i(\rr) \label{eq:10-frac-dE} \end{equation} $$

になり，拘束項と合わせて対角化すれば($n_i\neq0$ の軌道について $n_i$ で割って)

$$ \begin{equation} \hat{h}_{\mathrm{KS}}\phi_i=\varepsilon_i\phi_i \label{eq:10-frac-ks} \end{equation} $$

という，形の上では全く同じKS方程式が得られる．占有数は $v_{\mathrm{eff}}$ を作る密度 \eqref{eq:10-frac-density} を通してのみ方程式に効く．以下では，この自己無撞着解が任意の $\{n_i\}$ に対して得られているものとする．

10.7.3 Janakの定理

定理:Janakの定理(1978年，文献[4])

式 \eqref{eq:10-frac-density}–\eqref{eq:10-frac-energy} で定義された全エネルギー $E$ を，占有数 $\{n_j\}$ を与えたうえでKS方程式 \eqref{eq:10-frac-ks} を自己無撞着に解いて得た値と見なす．このとき，$n_i$ についての微係数は $i$ 番目のKS固有値に等しい:

$$ \begin{equation} \frac{\partial E}{\partial n_i}=\varepsilon_i. \label{eq:10-janak} \end{equation} $$

証明に入る前に，この主張の構造を確認しておく．$n_i$ を変えると，(a) 式 \eqref{eq:10-frac-energy} に陽に現れる係数 $n_i$ が変わるだけでなく，(b) 密度が変わって $v_{\mathrm{eff}}$ が変わり，その結果軌道 $\phi_j$ 自体もすべて変化する．一見すると (b) の寄与を全部追跡しなければならないように見える．しかし，軌道は $E$ を停留にするように決まっているため，(b) の寄与は丸ごと消える．この「停留点ではパラメータ微分から陰的な依存性が落ちる」という仕組みを，先に数学ノートで独立に整理しておく．

数学ノート:停留点でのパラメータ微分(包絡線定理)

2変数関数 $f(x,y)$ を考え，各 $x$ ごとに $y$ を「$f$ を停留にする値」$y^\star(x)$ に選ぶ．すなわち

$$ \begin{equation} \left.\frac{\partial f}{\partial y}\right|_{(x,\,y^\star(x))}=0 \label{eq:10-envelope-cond} \end{equation} $$

とし，$F(x)\equiv f\bigl(x,y^\star(x)\bigr)$ とおく．$F$ を $x$ で微分すると，合成関数の微分則により

\begin{align} \frac{\dd F}{\dd x} &=\left.\frac{\partial f}{\partial x}\right|_{y=y^\star} +\left.\frac{\partial f}{\partial y}\right|_{y=y^\star}\frac{\dd y^\star}{\dd x} \label{eq:10-envelope-1} \\ &=\left.\frac{\partial f}{\partial x}\right|_{y=y^\star}. \label{eq:10-envelope-2} \end{align}

\eqref{eq:10-envelope-1} は全微分の公式，\eqref{eq:10-envelope-2} は停留条件 \eqref{eq:10-envelope-cond} により第2項が消えたことによる．つまり停留点で評価する限り，$y$ が $x$ にどう依存するかを知らなくても，$x$ による偏微分だけ計算すれば全微分が得られる．$\dd y^\star/\dd x$ は掛かる係数がゼロなので，いくら複雑でも構わない．

これは包絡線定理と呼ばれる一般的な原理であり，第16章のヘルマン・ファインマンの定理(核座標で微分するとき波動関数の変化を追わなくてよい，という定理)も同じ仕組みである．以下の証明では $x\to n_i$，$y\to\{\phi_j,\phi_j^*\}$(無限次元だが論理は同じ)と読み替える．

導出:Janakの定理 $\partial E/\partial n_i=\varepsilon_i$

ステップ1:全微分を陽な項と陰な項に分ける．$E$ は $\{n_j\}$ と $\{\phi_j,\phi_j^*\}$ の関数(汎関数)であり，$\phi_j$ は $n_i$ に依存する．全微分は

\begin{align} \frac{\dd E}{\dd n_i} &=\left(\frac{\partial E}{\partial n_i}\right)_{\{\phi\}} +\sum_{k}\int \dd^3r\left[ \frac{\delta E}{\delta \phi_k^*(\rr)}\frac{\partial \phi_k^*(\rr)}{\partial n_i} +\frac{\delta E}{\delta \phi_k(\rr)}\frac{\partial \phi_k(\rr)}{\partial n_i} \right] \label{eq:10-janak-1} \end{align}

と書ける．第1項は軌道を固定したまま係数 $n_i$ だけを動かす寄与(陽な依存性)，第2項は軌道の変化を通じた寄与(陰な依存性)である．$\phi_k$ と $\phi_k^*$ の両方を書いたのは，10.2節の数学ノートのとおり実部・虚部の2自由度を数え落とさないためである．

ステップ2:陰な項が消えることを示す．停留条件 \eqref{eq:10-frac-dE} と，そのエルミート共役(実数値汎関数なので $\delta E/\delta\phi_k=\bigl(\delta E/\delta\phi_k^*\bigr)^*$)を代入する．正準形では $\hat{h}_{\mathrm{KS}}\phi_k=\varepsilon_k\phi_k$ だから $\delta E/\delta\phi_k^*=n_k\varepsilon_k\phi_k$，$\delta E/\delta\phi_k=n_k\varepsilon_k\phi_k^*$($\varepsilon_k$ は実数)である．よって \eqref{eq:10-janak-1} の第2項は

\begin{align} \sum_{k}\int \dd^3r\;n_k\varepsilon_k \left[\phi_k(\rr)\frac{\partial \phi_k^*(\rr)}{\partial n_i} +\phi_k^*(\rr)\frac{\partial \phi_k(\rr)}{\partial n_i}\right] &=\sum_{k}n_k\varepsilon_k\,\frac{\partial}{\partial n_i}\int \dd^3r\,\phi_k^*(\rr)\phi_k(\rr) \label{eq:10-janak-2} \\ &=\sum_{k}n_k\varepsilon_k\,\frac{\partial}{\partial n_i}(1) \;=\;0. \label{eq:10-janak-3} \end{align}

\eqref{eq:10-janak-2} では積の微分則 $\partial(\phi_k^*\phi_k)/\partial n_i=\phi_k(\partial\phi_k^*/\partial n_i)+\phi_k^*(\partial\phi_k/\partial n_i)$ を逆向きに使って括弧の中を1つの微分にまとめ，さらに微分と積分の順序を交換した．\eqref{eq:10-janak-3} では規格化条件 $\int\abs{\phi_k}^2\dd^3r=1$ が $n_i$ の値によらず恒等的に成り立つことを使った．定数を微分すればゼロである．これが包絡線定理の具体的な現れであり，軌道の変化を一切追跡せずに済む理由である．

(非正準形のまま計算しても結論は同じである．$\delta E/\delta\phi_k^*=\sum_l\Lambda_{lk}\phi_l$ と書くと，第2項は $\sum_{k,l}\Lambda_{lk}\,\partial\braket{\phi_k|\phi_l}/\partial n_i$ となり，規格直交条件 $\braket{\phi_k|\phi_l}=\delta_{kl}$ が恒等的に成り立つのでやはりゼロになる．)

ステップ3:陽な項を計算する．残ったのは，軌道を固定したまま $n_i$ で偏微分する項だけである．式 \eqref{eq:10-frac-energy} の4項を順に微分する．まず $T_s$ は $n_i$ について陽に線形だから

$$ \begin{equation} \left(\frac{\partial T_s}{\partial n_i}\right)_{\{\phi\}} =\int \dd^3r\,\phi_i^*(\rr)\left(-\frac{1}{2}\nabla^2\right)\phi_i(\rr). \label{eq:10-janak-ts} \end{equation} $$

残りの3項は $\rho$ を通してのみ $n_i$ に依存する．軌道を固定したときの密度の変化は，\eqref{eq:10-frac-density} から

$$ \begin{equation} \left(\frac{\partial \rho(\rr)}{\partial n_i}\right)_{\{\phi\}}=\abs{\phi_i(\rr)}^2 \label{eq:10-janak-drho} \end{equation} $$

である．したがって，$\rho$ の汎関数 $A[\rho]$ については連鎖律により

$$ \begin{equation} \left(\frac{\partial A}{\partial n_i}\right)_{\{\phi\}} =\int \dd^3r\,\frac{\delta A}{\delta\rho(\rr)}\left(\frac{\partial\rho(\rr)}{\partial n_i}\right)_{\{\phi\}} =\int \dd^3r\,\frac{\delta A}{\delta\rho(\rr)}\abs{\phi_i(\rr)}^2 \label{eq:10-janak-chain} \end{equation} $$

となる．これを3項それぞれに適用する．$\delta\bigl(\int v\rho\bigr)/\delta\rho=v$，$\delta E_{\mathrm{H}}/\delta\rho=v_{\mathrm{H}}$(式 \eqref{eq:10-vh-def})，$\delta E_{xc}/\delta\rho=v_{xc}$(式 \eqref{eq:10-vxc-def})だから

\begin{align} \left(\frac{\partial}{\partial n_i}\int \dd^3r\,v\rho\right)_{\{\phi\}}&=\int \dd^3r\,v(\rr)\abs{\phi_i(\rr)}^2, \label{eq:10-janak-vext}\\ \left(\frac{\partial E_{\mathrm{H}}}{\partial n_i}\right)_{\{\phi\}}&=\int \dd^3r\,v_{\mathrm{H}}(\rr)\abs{\phi_i(\rr)}^2, \label{eq:10-janak-vh}\\ \left(\frac{\partial E_{xc}}{\partial n_i}\right)_{\{\phi\}}&=\int \dd^3r\,v_{xc}(\rr)\abs{\phi_i(\rr)}^2. \label{eq:10-janak-vxc} \end{align}

ステップ4:足し合わせてKS方程式を使う．\eqref{eq:10-janak-ts} と \eqref{eq:10-janak-vext}–\eqref{eq:10-janak-vxc} を合計すると

\begin{align} \frac{\dd E}{\dd n_i} &=\int \dd^3r\,\phi_i^*\left(-\frac{1}{2}\nabla^2\right)\phi_i +\int \dd^3r\,\bigl[v+v_{\mathrm{H}}+v_{xc}\bigr]\,\phi_i^*\phi_i \label{eq:10-janak-4} \\ &=\int \dd^3r\,\phi_i^*(\rr) \left[-\frac{1}{2}\nabla^2+v_{\mathrm{eff}}(\rr)\right]\phi_i(\rr) \label{eq:10-janak-5} \\ &=\int \dd^3r\,\phi_i^*(\rr)\,\varepsilon_i\,\phi_i(\rr) \label{eq:10-janak-6} \\ &=\varepsilon_i\int \dd^3r\,\abs{\phi_i(\rr)}^2 \;=\;\varepsilon_i. \label{eq:10-janak-7} \end{align}

\eqref{eq:10-janak-4} では $\abs{\phi_i}^2=\phi_i^*\phi_i$ と書き直して3つの積分を1つにまとめ，\eqref{eq:10-janak-5} では有効ポテンシャルの定義 \eqref{eq:10-veff-def} を使った．\eqref{eq:10-janak-6} はKS方程式 \eqref{eq:10-frac-ks} の代入，\eqref{eq:10-janak-7} は $\varepsilon_i$ が定数なので積分の外に出し，規格化条件を使った．以上で $\partial E/\partial n_i=\varepsilon_i$ が示された．∎

物理的意味:固有値は「限界エネルギー」である

Janakの定理は，$\varepsilon_i$ を微分量として特徴づける．すなわち $\varepsilon_i$ は，軌道 $i$ に無限小の電子を追加(または除去)するときの，電子1個あたりのエネルギー変化である．経済学の用語を借りれば「限界費用」であり，「電子1個分のエネルギー」という有限量ではない．この区別が，$\sum_i\varepsilon_i\neq E$(10.4節)という事実と矛盾なく共存する理由である．実際，Janakの定理を全占有数について積分すると

$$ \begin{equation} E(N)-E(0)=\int_0^1\dd\lambda\,\frac{\dd E}{\dd\lambda} =\int_0^1\dd\lambda\sum_{i=1}^{N}\varepsilon_i(\lambda) \label{eq:10-janak-integral-all} \end{equation} $$

($n_i=\lambda$ を $0\to1$ まで一斉に上げる経路をとった)となり，$\varepsilon_i$ が占有数に依存するため単純な $\sum_i\varepsilon_i$ にはならない．固有値の和が全エネルギーでないのは，まさにこの $\lambda$ 依存性(=電子を入れると自分が作った場が変わること)のためである．

10.7.4 積分形:最高占有準位とイオン化エネルギー

Janakの定理は微分の形をしているので，占有数について積分すれば有限のエネルギー差が得られる．軌道 $i$ の占有数だけを $1$ から $0$ まで下げ，他はそのままにする経路を考えると，微積分学の基本定理により

\begin{align} E(n_i{=}1)-E(n_i{=}0) &=\int_0^1 \dd n\,\frac{\partial E}{\partial n_i}\bigg|_{n_i=n} \notag\\ &=\int_0^1 \dd n\,\varepsilon_i(n) \label{eq:10-janak-int} \end{align}

となる．2行目でJanakの定理 \eqref{eq:10-janak} を代入した．$\varepsilon_i(n)$ は「占有数を $n$ にして自己無撞着に解いたときの固有値」であり，$n$ に依存することに注意する．軌道 $i$ から電子を1個抜くのに要するエネルギー(イオン化エネルギー)を $I_i\equiv E(n_i{=}0)-E(n_i{=}1)$ と定義すれば

$$ \begin{equation} I_i=-\int_0^1 \dd n\,\varepsilon_i(n). \label{eq:10-Ii-integral} \end{equation} $$

ここで「$\varepsilon_i$ は $n$ によらない」と近似すれば $I_i\simeq-\varepsilon_i(1)$ となり，HF法のKoopmansの定理(第5章)と同じ形の関係が得られる．しかしこれはあくまで近似である．実際には電子を抜くと残った電子が感じる遮蔽が減るので $v_{\mathrm{eff}}$ は深くなり，$\varepsilon_i(n)$ は $n$ が減るにつれて下がる(より負になる)．この軌道緩和の効果を無視した分だけ $-\varepsilon_i(1)$ は $I_i$ を過小評価する．

ところが，最高占有準位(HOMO)については，厳密汎関数を使えば近似抜きで厳密な関係が成り立つ．鍵になるのは，厳密な $E$ が電子数 $\nu$ の関数として折れ線になるという事実である．

導出:全エネルギーの電子数依存性は折れ線になる(Perdew–Parr–Levy–Balduz)

電子浴と電子をやりとりできる系を考える(10.7.2節の数学ノート(1))．粒子数が確定していない状態は，粒子数 $M$($M=0,1,2,\dots$)の $N$ 電子基底状態 $\Psi_M$ を重み $w_M\ge0$ で混ぜた統計混合

$$ \begin{equation} \hat{\Gamma}=\sum_{M}w_M\ket{\Psi_M}\bra{\Psi_M}, \qquad \sum_M w_M=1,\qquad \sum_M w_M M=\nu \label{eq:10-ensemble} \end{equation} $$

で記述される(第1式は確率の規格化，第2式は平均粒子数が $\nu$ であるという条件)．このアンサンブルのエネルギー期待値は，$\ket{\Psi_M}$ が $M$ 電子ハミルトニアンの基底状態であることから

$$ \begin{equation} \mathcal{E}[\hat\Gamma]=\mathrm{Tr}\bigl[\hat{\Gamma}\hat{H}\bigr]=\sum_M w_M E(M) \label{eq:10-ensemble-energy} \end{equation} $$

である．基底状態エネルギー $E(\nu)$ は，拘束 \eqref{eq:10-ensemble} のもとでこの線形関数 $\sum_M w_M E(M)$ を $\{w_M\}$ について最小化した値として定義される．

ここからが要点である．「非負の重み $w_M$，和が1，平均が $\nu$」という条件を満たす点 $\bigl(\nu,\sum_M w_M E(M)\bigr)$ の全体は，平面上の点列 $\{(M,E(M))\}$ の凸結合の集合，すなわちそれらの点の凸包にほかならない．その最小値(下側の境界)は，点列の下側凸包である．実在の系では $E(M)$ は $M$ について凸，すなわち

$$ \begin{equation} E(M-1)+E(M+1)\ge 2E(M) \qquad\Longleftrightarrow\qquad I\ge A \label{eq:10-convexity} \end{equation} $$

が成り立つ($I=E(M{-}1)-E(M)$ はイオン化エネルギー，$A=E(M)-E(M{+}1)$ は電子親和力．$I\ge A$ は「電子を1個抜くほうが1個加えるより高くつく」という経験則であり，実測でも例外は知られていない)．凸な点列の下側凸包は，隣り合う点を結んだ折れ線そのものである．したがって $\nu=M+\omega$($0\le\omega\le1$)に対し，最小を与える重みは $w_M=1-\omega$, $w_{M+1}=\omega$，他はゼロであり，

$$ \begin{equation} E(M+\omega)=(1-\omega)E(M)+\omega E(M+1) \label{eq:10-pwl} \end{equation} $$

が得られる．すなわち $E(\nu)$ は整数点を頂点とする折れ線である．∎

折れ線であることから，微係数(化学ポテンシャル)$\mu(\nu)=\partial E/\partial\nu$ は各区間で定数となり，整数点で跳ぶ:

\begin{align} \mu(\nu)&=E(N)-E(N-1)=-I \qquad (N-1\lt\nu\lt N), \label{eq:10-mu-left} \\ \mu(\nu)&=E(N+1)-E(N)=-A \qquad (N\lt\nu\lt N+1). \label{eq:10-mu-right} \end{align}

一方，$\nu$ を $N-1$ から $N$ へ増やす過程では，加わる電子は最高占有軌道(HOMO)に入っていくから，Janakの定理により $\partial E/\partial\nu=\varepsilon_{\mathrm{HOMO}}(\nu)$ である．式 \eqref{eq:10-mu-left} の右辺は $\nu$ によらない定数だから，$\varepsilon_{\mathrm{HOMO}}(\nu)$ もこの区間全体で定数でなければならない．とくに $\nu\to N$ の極限，つまり通常の $N$ 電子系のHOMO固有値について

$$ \begin{equation} \varepsilon_{\mathrm{HOMO}}=-I \label{eq:10-homo-ip} \end{equation} $$

が厳密に成り立つ(文献[5], [10])．これがKS固有値に与えられた唯一の厳密な物理的意味である．同じ議論を $\nu\gt N$ 側で行えば，$N+1$ 電子系のHOMO(=$N$ 電子系で見れば最低空軌道に対応する準位)について $\varepsilon_{\mathrm{HOMO}}(N{+}1)=-A$ が言える．ただし後で見るように，これは $N$ 電子系の $\varepsilon_{\mathrm{LUMO}}$ とは一致しない．

E(ν) 電子数 ν N−1 N N+1 E(N−1) E(N) E(N+1) 傾き = −I ( = εHOMO ) 傾き = −A ν = N で微分が不連続 跳び = I − A = Eg 厳密汎関数:折れ線 LDA/GGA:下に凸の曲線 HF:上に凸の曲線
図10.3 全エネルギーの電子数依存性．厳密汎関数では $E(\nu)$ は整数点を頂点とする折れ線になり(式 \eqref{eq:10-pwl})，各区間の傾きが $-I$ および $-A$，整数点での傾きの跳びが基本ギャップ $E_g=I-A$ を与える．LDA/GGAなどの近似汎関数は $\rho$ のなめらかな関数なので $E(\nu)$ もなめらかな下に凸の曲線になり，整数点の折れ目を再現できない．逆にHF近似の $E(\nu)$ は各区間で上に凸に曲がる．この曲がる向きがバンドギャップ誤差の向きを決めることを10.7.6節で見る．

物理的意味:折れ線性の破れ=非局在化誤差

図10.3の折れ線は「電子は整数個でいるのが損でも得でもない」ことを意味する．近似汎関数の凸曲線は折れ線より下にあるので，分数電荷を持つほうがエネルギー的に得になってしまう．この誤差は非局在化誤差(delocalization error)と呼ばれ，次のような破綻を生む．

HF近似はこの逆で，$E(\nu)$ が各区間で折れ線より上に凸に外れる(局在化誤差と呼ばれる)．交換を厳密に扱う一方で遮蔽(相関)を欠くため，分数電荷を過剰に嫌うのである．下に凸・上に凸という曲がり方の違いが，バンドギャップの誤差の向きをそのまま決めることを，10.7.6節で図10.3に戻って示す．

例:$\varepsilon_{\mathrm{HOMO}}=-I$ は近似汎関数では大きく破れる

ヘリウム原子のイオン化エネルギーは $I=24.59$ eV である．式 \eqref{eq:10-homo-ip} によれば，厳密汎関数のKS計算は $\varepsilon_{\mathrm{HOMO}}=-24.59$ eV を与えなければならない．ところがLDA計算では $\varepsilon_{\mathrm{HOMO}}\simeq-15.5$ eV にしかならず，$I$ の 6 割程度しかない．原子・分子で「LDA/GGAのHOMO準位はイオン化エネルギーの半分程度」という傾向は広く知られている．

原因は $v_{xc}$ の遠方での振る舞いにある．厳密な交換相関ポテンシャルは，$N$ 電子系から電子1個が遠ざかったとき，残った $N-1$ 個の正イオンからのクーロン引力を再現しなければならないので $v_{xc}(\rr)\to-1/r$($r\to\infty$)と減衰する(文献[10])．これに対しLDAの $v_{xc}\propto-\rho^{1/3}$ は，密度が指数関数的に減衰する遠方で指数関数的にゼロになってしまう．ポテンシャルの尾が浅すぎるため，HOMOの束縛が弱くなりすぎるのである．一方，次項で述べる $\Delta$SCF法は全エネルギーの差を使うので，この欠陥の影響をほとんど受けない(LDAでも $I\simeq24.3$ eV 程度と良好な値を与える)．

10.7.5 $\Delta$SCF法とSlaterの遷移状態

前項で見たとおり，固有値をそのまま読んでイオン化エネルギーとするのは近似汎関数では危険である．では実際にはどうするか．答えは「固有値ではなく全エネルギーを使う」ことである．

(a) $\Delta$SCF法．イオン化エネルギーの定義そのものに戻り，$N$ 電子系と $N-1$ 電子系のSCF計算を別々に行って差をとる:

$$ \begin{equation} I^{\Delta\mathrm{SCF}}=E_{\mathrm{SCF}}(N-1)-E_{\mathrm{SCF}}(N). \label{eq:10-dscf} \end{equation} $$

これを$\Delta$SCF法と呼ぶ．この方法が優れている理由は3つある．

短所は，状態ごとに1回ずつSCF計算が必要なこと，そして内殻準位のように「特定の軌道から電子を抜いた励起状態」を扱うには占有数を人為的に固定する拘束付きSCFが要ることである．また周期系では，荷電したセルを扱うために背景電荷の補正が必要になる．内殻励起への応用は第18章で詳しく扱う．

(b) Slaterの遷移状態近似．$\Delta$SCF は2回の計算を要する．1回の計算で同等の精度を出そうというのがSlaterの遷移状態(transition state)近似である．発想は「積分 \eqref{eq:10-Ii-integral} を中点で評価する」という，きわめて単純なものである．

導出:Slaterの遷移状態と，その誤差の評価

式 \eqref{eq:10-Ii-integral} の被積分関数 $\varepsilon_i(n)$ を，区間の中点 $n=1/2$ のまわりでテイラー展開する:

$$ \begin{equation} \varepsilon_i(n)=\varepsilon_i(\tfrac12) +\left(n-\tfrac12\right)\varepsilon_i'(\tfrac12) +\frac{1}{2}\left(n-\tfrac12\right)^2\varepsilon_i''(\tfrac12) +\frac{1}{6}\left(n-\tfrac12\right)^3\varepsilon_i'''(\tfrac12) +\cdots \label{eq:10-ts-taylor} \end{equation} $$

($'$ は $n$ による微分)．これを $0$ から $1$ まで積分する．変数変換 $u=n-\frac12$ を行うと積分区間は $-\frac12\le u\le\frac12$ という原点対称な区間になり，各項の積分は

\begin{align} \int_{-1/2}^{1/2}\dd u&=1, \qquad \int_{-1/2}^{1/2}u\,\dd u=0, \label{eq:10-ts-mom1}\\ \int_{-1/2}^{1/2}u^2\,\dd u&=\left[\frac{u^3}{3}\right]_{-1/2}^{1/2}=\frac{2}{3}\cdot\frac{1}{8}=\frac{1}{12}, \qquad \int_{-1/2}^{1/2}u^3\,\dd u=0 \label{eq:10-ts-mom2} \end{align}

となる．奇数次のべきの積分が消えるのは，対称区間における奇関数の積分だからである(これが中点を選ぶことの御利益である)．したがって

\begin{align} \int_0^1\dd n\,\varepsilon_i(n) &=\varepsilon_i(\tfrac12)\cdot 1 +\varepsilon_i'(\tfrac12)\cdot 0 +\frac{1}{2}\varepsilon_i''(\tfrac12)\cdot\frac{1}{12} +\frac{1}{6}\varepsilon_i'''(\tfrac12)\cdot 0+\cdots \label{eq:10-ts-int1} \\ &=\varepsilon_i(\tfrac12)+\frac{1}{24}\varepsilon_i''(\tfrac12)+O\!\left(\varepsilon_i''''\right). \label{eq:10-ts-int2} \end{align}

式 \eqref{eq:10-Ii-integral} と合わせて

$$ \begin{equation} I_i=-\varepsilon_i(\tfrac12)-\frac{1}{24}\varepsilon_i''(\tfrac12)-\cdots \;\simeq\;-\varepsilon_i\!\left(n_i=\tfrac12\right). \label{eq:10-slater-ts} \end{equation} $$

すなわち，軌道 $i$ の占有数を $1/2$ に固定して自己無撞着に解き，その固有値の符号を変えたものがイオン化エネルギーである．この半占有状態を遷移状態と呼ぶ．

誤差の比較．他の評価法と並べてみよう．

とくに，$\varepsilon_i(n)$ が $n$ の1次関数であれば $\varepsilon_i''=0$ となり，遷移状態近似は厳密になる．∎

物理的意味:なぜ半占有なのか

$\varepsilon_i(n)$ が $n$ についてほぼ線形であることは，図10.3で見た近似汎関数の $E(\nu)$ がほぼ放物線であることと同じ事実の言い換えである．実際，$E(n_i)=an_i^2+bn_i+c$ と2次式で近似すればJanakの定理より $\varepsilon_i(n_i)=\partial E/\partial n_i=2an_i+b$ となり，$n_i$ の厳密な1次関数になる．このとき中点則は誤差ゼロである．

直観的にはこうである．電子を抜く途中では，系は「まだ電子がある状態」と「もう電子がない状態」の中間にいる．$-\varepsilon_i(1)$ は「抜く前」の場を，$-\varepsilon_i(0)$ は「抜いた後」の場を使った評価であり，それぞれ緩和を無視する側と過剰に見積もる側に外れる．半占有はそのちょうど中間の場を使うので，緩和の効果を一次のオーダーで正しく取り込む．$\Delta$SCFが2回の計算で緩和を厳密に扱うのに対し，遷移状態は1回の計算で緩和を一次まで取り込む方法だと言える．

例:内殻準位の束縛エネルギー

X線光電子分光(XPS)で測る内殻準位の束縛エネルギー(たとえば炭素の $1s$ で約 285 eV)は，KS固有値をそのまま読むと数十 eV も浅く出る．内殻に空孔ができると価電子が強く引き寄せられ，非常に大きな緩和が起こるからである．$\Delta$SCF(内殻軌道の占有数を 0 に固定した拘束付きSCF)や遷移状態(占有数 $1/2$)を使うと，この緩和が取り込まれ，絶対値で 0.5 eV 程度の精度に到達できる．異なる化学環境にある同種原子の束縛エネルギー差(化学シフト)であればさらに精度は上がる．詳しくは第18章で扱う．

10.7.6 実験との比較:光電子分光とバンドギャップ問題

理論的な整理が済んだところで，実際に $\varepsilon_i$ を実験と比べるとどうなるかを見る．結論を先に言えば，占有準位は数 eV の精度で使えるが，バンドギャップは系統的に，しかも大幅に過小評価される．この非対称な振る舞いには明確な理由があり，それが本節の締めくくりとなる汎関数微分の不連続である．

(a) 占有準位と光電子分光．紫外・X線光電子分光(UPS/XPS)や角度分解光電子分光(ARPES)は，試料に光子を当てて飛び出す電子の運動エネルギーを測り，束縛エネルギー $E_b=E(N-1,\text{終状態})-E(N)$ を決める実験である．これを $-\varepsilon_i$ と比べると，次のことが知られている．

固体で相対的にうまくいく理由は，Janakの定理の積分形 \eqref{eq:10-Ii-integral} から理解できる．マクロな固体で1個の電子を非局在なバンド状態から抜いても，密度の変化は $\delta\rho\sim1/N\to0$ であり，$v_{\mathrm{eff}}$ はほとんど変わらない．したがって $\varepsilon_i(n)$ の $n$ 依存性(軌道緩和)が消え，$I_i=-\int_0^1\varepsilon_i(n)\dd n\to-\varepsilon_i$ が良い近似になる．逆に，孤立原子や内殻準位のように電子の抜き差しで場が大きく変わる場合には，この論法は使えない．

注意:準粒子と $v_{xc}$ の違い

そもそも光電子分光が測るのは，$N$ 電子系から電子を1個取り去った $N-1$ 電子系の励起状態のエネルギーである．多体論の言葉では，この量は準粒子方程式

$$ \begin{equation} \left[-\frac{1}{2}\nabla^2+v(\rr)+v_{\mathrm{H}}(\rr)\right]\psi_i(\rr) +\int \dd^3r'\,\Sigma\bigl(\rr,\rr';E_i\bigr)\,\psi_i(\rr')=E_i\,\psi_i(\rr) \label{eq:10-qp} \end{equation} $$

で決まる．$\Sigma(\rr,\rr';E)$ は自己エネルギーと呼ばれ，$v_{xc}(\rr)$ と違って (i) 非局所($\rr$ と $\rr'$ の2点関数)であり，(ii) エネルギー依存性を持ち，(iii) 一般に複素数(虚部が準粒子の有限寿命を与える)である．KS方程式は $\Sigma$ を局所・実・エネルギー非依存の $v_{xc}(\rr)$ で置き換えたものと見なせるから，$\varepsilon_i$ が $E_i$ に一致する理由はもともとない．$\Sigma$ を遮蔽クーロン相互作用の1次で近似する $GW$ 近似は，この差を系統的に埋める標準的手法である．

(b) バンドギャップの過小評価．半導体・絶縁体では事情がもっと深刻である．KS固有値から読んだギャップ $E_g^{\mathrm{KS}}=\varepsilon_{\mathrm{LUMO}}-\varepsilon_{\mathrm{HOMO}}$(固体では伝導帯下端と価電子帯上端の差)は，実験の基本ギャップ $E_g=I-A$ を一貫して小さく見積もる．

表10.2 LDAで計算したバンドギャップ $E_g^{\mathrm{KS}}$ と実験値の比較(単位:eV)．代表的な文献値であり，計算条件により多少ばらつく．
物質LDA の $E_g^{\mathrm{KS}}$実験値 $E_g$誤差
Si0.51.17約 $-57\%$
Ge0.0(金属と誤予測)0.74約 $-100\%$
GaAs0.31.52約 $-80\%$
C(ダイヤモンド)4.15.48約 $-25\%$
MgO4.87.83約 $-39\%$
LiF8.914.20約 $-37\%$

「LDAは近似が粗いから」で片付けてはいけない．驚くべきことに，厳密な $E_{xc}$ を使ったとしても $E_g^{\mathrm{KS}}\neq E_g$ である．両者の差を与えるのが，次に導く汎関数微分の不連続 $\Delta_{xc}$ である(文献[7], [8])．

導出:$E_g=E_g^{\mathrm{KS}}+\Delta_{xc}$(汎関数微分の不連続)

道具はすべて 10.3 節と 10.7.4 節で用意済みである．相互作用系のオイラー方程式 \eqref{eq:10-euler-int2} を再掲する:

$$ \begin{equation} \frac{\delta T_s}{\delta\rho(\rr)}+v(\rr)+v_{\mathrm{H}}(\rr)+v_{xc}(\rr)=\mu. \label{eq:10-dd-euler} \end{equation} $$

これを，電子数 $\nu$ を整数 $N$ に下から近づけた極限($\nu=N-\eta$, $\eta\to0^+$)と，上から近づけた極限($\nu=N+\eta$)の2通りで評価し，辺々引き算する．以下，上付き $\mp$ でそれぞれの極限を表す．

ステップ1:密度は連続，化学ポテンシャルは不連続．式 \eqref{eq:10-pwl} のアンサンブル描像では，$\nu=N\pm\eta$ の密度は $\rho_{N\pm\eta}=(1\mp\eta)\rho_N\pm\eta\,\rho_{N\pm1}$ の形をしており，$\eta\to0$ で $\rho_N$ に連続的に近づく．したがって $\rho$ の関数である $v(\rr)$(そもそも $\rho$ によらない)と $v_{\mathrm{H}}(\rr)$(密度の線形汎関数)は，両極限で同じ値をとる．一方 $\mu$ は \eqref{eq:10-mu-left}, \eqref{eq:10-mu-right} により

$$ \begin{equation} \mu^{-}=-I,\qquad \mu^{+}=-A,\qquad \mu^{+}-\mu^{-}=I-A=E_g \label{eq:10-dd-mu} \end{equation} $$

と跳ぶ．

ステップ2:$\delta T_s/\delta\rho$ の跳びはKSギャップである．非相互作用系のオイラー方程式 \eqref{eq:10-euler-ks} を $\delta T_s/\delta\rho$ について解くと

$$ \begin{equation} \frac{\delta T_s}{\delta\rho(\rr)}\bigg|_{\rho_N}=\mu_s-v_s(\rr) \label{eq:10-dd-ts} \end{equation} $$

である．ここで $v_s(\rr)$ は密度 $\rho_N$ を再現する局所ポテンシャルであり，$\rho_N$ が両極限で共通なのだから $v_s$ も共通である(10.3.2節で見たHK定理による一意性)．残るのは非相互作用系の化学ポテンシャル $\mu_s$ であるが，これは「$v_s$ 中の独立粒子系に無限小の電子を出し入れするときのエネルギー」だから，下から近づければ最高占有準位の $\varepsilon_{\mathrm{HOMO}}$，上から近づければ最低空準位の $\varepsilon_{\mathrm{LUMO}}$ である(相互作用がないので，追加した電子は空いている最低の準位にそのまま入る)．よって

$$ \begin{equation} \frac{\delta T_s}{\delta\rho}\bigg|^{+}-\frac{\delta T_s}{\delta\rho}\bigg|^{-} =\mu_s^{+}-\mu_s^{-} =\varepsilon_{\mathrm{LUMO}}-\varepsilon_{\mathrm{HOMO}} \equiv E_g^{\mathrm{KS}}. \label{eq:10-dd-ts-jump} \end{equation} $$

ステップ3:引き算する．\eqref{eq:10-dd-euler} の上側極限から下側極限を引くと，$v$ と $v_{\mathrm{H}}$ は相殺し

\begin{align} \left[\frac{\delta T_s}{\delta\rho}\bigg|^{+}-\frac{\delta T_s}{\delta\rho}\bigg|^{-}\right] +\left[v_{xc}^{+}(\rr)-v_{xc}^{-}(\rr)\right] &=\mu^{+}-\mu^{-} \label{eq:10-dd-sub} \\ E_g^{\mathrm{KS}}+\underbrace{\left[v_{xc}^{+}(\rr)-v_{xc}^{-}(\rr)\right]}_{\displaystyle \equiv\,\Delta_{xc}} &=E_g \label{eq:10-dd-sub2} \end{align}

となる．\eqref{eq:10-dd-sub2} では \eqref{eq:10-dd-ts-jump} と \eqref{eq:10-dd-mu} を代入した．ここで重要な副産物がある:$E_g^{\mathrm{KS}}$ も $E_g$ も数(位置によらない)なのだから，差 $\Delta_{xc}$ も位置 $\rr$ によらない定数でなければならない．すなわち，電子数が整数を横切るとき，交換相関ポテンシャルは空間的に一様な定数だけ跳ぶ．整理して:

$$ \begin{equation} E_g=E_g^{\mathrm{KS}}+\Delta_{xc}, \qquad \Delta_{xc}=\lim_{\eta\to0^+}\left[ \frac{\delta E_{xc}}{\delta\rho(\rr)}\bigg|_{N+\eta} -\frac{\delta E_{xc}}{\delta\rho(\rr)}\bigg|_{N-\eta}\right]. \label{eq:10-dd-result} \end{equation} $$

∎

$$ \begin{equation} E_g=I-A=\bigl(\varepsilon_{\mathrm{LUMO}}-\varepsilon_{\mathrm{HOMO}}\bigr)+\Delta_{xc} =E_g^{\mathrm{KS}}+\Delta_{xc} \label{eq:10-gap-key} \end{equation} $$
エネルギー より深い占有準位 εHOMO(N) = −I εLUMO(N)(N電子系のKS空準位) −A = εHOMO(N+1) EgKS Δxc Eg = I − A ν が整数 N を横切ると，vxc(r) が空間的に一様な定数 Δxc だけ跳ぶ LDA/GGA は Δxc ≡ 0 → EgKS しか得られず，ギャップは系統的に過小評価
図10.4 基本ギャップ $E_g=I-A$ とコーン・シャムギャップ $E_g^{\mathrm{KS}}=\varepsilon_{\mathrm{LUMO}}-\varepsilon_{\mathrm{HOMO}}$ の関係．両者の差が汎関数微分の不連続 $\Delta_{xc}$ である(式 \eqref{eq:10-gap-key})．$-A$ は $N+1$ 電子系のHOMO準位であり，$N$ 電子系のLUMO準位より $\Delta_{xc}$ だけ高い．

物理的意味:なぜLDA・GGAでは $\Delta_{xc}=0$ になるのか

LDAやGGAの交換相関エネルギーは

$$ E_{xc}^{\mathrm{LDA/GGA}}[\rho]=\int \dd^3r\;f\bigl(\rho(\rr),\nabla\rho(\rr)\bigr) $$

という形をしており，$f$ は $\rho$ と $\nabla\rho$ のなめらかな関数である．したがってその汎関数微分

$$ v_{xc}(\rr)=\frac{\partial f}{\partial\rho}-\nabla\cdot\frac{\partial f}{\partial(\nabla\rho)} $$

も $\rho$ と $\nabla\rho$ のなめらかな関数であり，密度が連続に変化する限り $v_{xc}$ も連続に変化する．上の導出のステップ1で見たとおり，電子数が整数を横切るとき密度は連続である．ゆえにLDA/GGAでは $\Delta_{xc}$ は恒等的にゼロであり，これらの汎関数は $E_g^{\mathrm{KS}}$ しか与えられない．仮にLDAが厳密な密度を再現できたとしても，ギャップは $\Delta_{xc}$ だけ足りないままである．Siでは $\Delta_{xc}$ は 0.6〜0.9 eV と見積もられており，表10.2の食い違い(0.67 eV)のほぼ全体を説明してしまう．

これに加えて，自己相互作用誤差(10.1節)と非局在化誤差(10.7.4節)も占有準位を押し上げてギャップを狭める方向に働く．したがって「バンドギャップ問題」は，汎関数の精度不足というよりKS形式そのものの構造に由来する部分が大きい．処方箋としては，(i) 厳密交換を一部混ぜて汎関数を軌道依存(非局所)にするハイブリッド汎関数——非局所性ゆえに $\Delta_{xc}\neq0$ となる——，(ii) 自己エネルギーを直接計算する $GW$ 近似，(iii) 有限系なら $\Delta$SCF，がある．(i) は第11章，(ii)(iii) は第18章で扱う．

物理的意味:図10.3で読むギャップ誤差の向き——下に凸のLDA/GGAは過小，上に凸のHFは過大

$\Delta_{xc}$ の導出は厳密だが，いささか抽象的でもある．同じ結論を，図10.3の $E(\nu)$ 曲線の曲がり方から直観的に読み直そう．鍵になるのは，実験のバンド端は「弦の傾き」，計算のバンド端は「接線の傾き」という区別である．

実験のバンド端=弦の傾き．価電子帯上端(VBM)と伝導帯下端(CBM)の位置は，定義により整数電子数のエネルギー差で決まるから，図10.3で隣り合う整数点を結んだ弦(差分)の傾きである:

$$ -I=\frac{E(N)-E(N-1)}{N-(N-1)} \;=\;\text{VBMの位置}, \qquad -A=\frac{E(N+1)-E(N)}{(N+1)-N} \;=\;\text{CBMの位置}. $$

これは式 \eqref{eq:10-mu-left}, \eqref{eq:10-mu-right} の言い換えである．

計算のバンド端=接線の傾き．一方，KS計算やHF計算が固有値として出力するバンド端は，Janakの定理により $E(\nu)$ の接線(微分)の傾きである:

$$ \varepsilon_{\mathrm{VBM}}=\lim_{\nu\to N^-}\frac{\partial E}{\partial\nu}, \qquad \varepsilon_{\mathrm{CBM}}=\lim_{\nu\to N^+}\frac{\partial E}{\partial\nu}. $$

(10.7.3節の導出は $E_{xc}$ の具体形を一切使っていないので，この読み方は近似汎関数にもHFにもそのまま通用する．)厳密汎関数では $E(\nu)$ が各区間で直線だから弦と接線は一致し，バンド端もギャップも正しく出る．近似で $E(\nu)$ が曲がると，曲がった分だけ接線が弦からずれる——そしてずれる向きは曲がる向きで決まる．

物理的な由来も対照的である．LDA/GGAの下に凸は自己相互作用誤差の名残——電子が自分自身と反発するため，電荷を分数に薄めて広げるほうが得になる——であり，HFの上に凸は遮蔽の欠如——電子を出し入れするコストを，他の電子の応答(遮蔽)を含まない裸のクーロン相互作用で見積もってしまう——である．

この見方は処方箋の理解にも直結する．ハイブリッド汎関数(第11章)がギャップを改善するのは，下に凸のLDA/GGAと上に凸のHFを混ぜると曲率が部分的に相殺して $E(\nu)$ が折れ線に近づくからである．DFT+U(第11章)の補正項 $\frac{U}{2}\sum_i n_i(1-n_i)$ も，分数占有にペナルティを課して下に凸を直線へ引き戻す働きと読める．

注意:「バンドギャップ」という言葉の3つの意味

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

10.8.1 まとめ

次章への橋渡し．本章で得られた枠組みは，「厳密な形式 + 近似はすべて $E_{xc}[\rho]$ に集約」という明快な分業である．方程式の形も，全エネルギーの表式も，SCFの手続きも，$E_{xc}$ の中身が何であるかを一切問わずに確定した．裏を返せば，DFT計算の精度は $E_{xc}$ の近似の質だけで決まる．第11章では，いよいよこの最後の未知項に踏み込む:一様電子ガスの厳密な結果(第7章)を局所的に持ち込む局所密度近似(LDA)，それが期待以上に働く理由(交換相関ホールの総和則)，密度勾配を取り込む一般化勾配近似(GGA)，そして厳密交換を混ぜるハイブリッド汎関数までを構成する．

10.8.2 演習問題

演習10.1 スピン分極系のJanakの定理

10.6節のスピン密度汎関数法の枠組みで，占有数を $n_{i\sigma}$($\sigma=\uparrow,\downarrow$)に一般化し，$\rho^\sigma=\sum_i n_{i\sigma}\abs{\phi_{i\sigma}}^2$，$\rho=\rho^\uparrow+\rho^\downarrow$ とする．このとき

$$ \frac{\partial E}{\partial n_{i\sigma}}=\varepsilon_{i\sigma} $$

が成り立つことを，10.7.3節の導出をなぞって示せ．

ヒント:軌道を固定したときの偏微分では $\partial\rho^{\sigma'}/\partial n_{i\sigma}=\delta_{\sigma\sigma'}\abs{\phi_{i\sigma}}^2$，$\partial\rho/\partial n_{i\sigma}=\abs{\phi_{i\sigma}}^2$ である．$E_{\mathrm{H}}$ と $\int v\rho$ は全密度 $\rho$ にしか依存しないので $v_{\mathrm{H}}$, $v$ がそのまま現れ，$E_{xc}$ からは $v_{xc}^\sigma=\delta E_{xc}/\delta\rho^\sigma$ が現れる．合計すると $v_{\mathrm{eff}}^\sigma$ が組み上がる．陰な項が消える論法(規格化条件が $n_{i\sigma}$ によらず成り立つこと)はそのまま使える．

演習10.2 遷移状態近似が厳密になる場合

(1) 全エネルギーが占有数 $n_i$ の2次式 $E(n_i)=an_i^2+bn_i+c$ で表されるとき，Janakの定理から $\varepsilon_i(n_i)$ を求め，Slaterの遷移状態近似 $I_i=-\varepsilon_i(1/2)$ が厳密に成り立つことを示せ．

(2) 同じ2次式のもとで，$\Delta$SCF の値 $E(0)-E(1)$ と遷移状態の値が一致することを確かめよ．

(3) Koopmans型 $-\varepsilon_i(1)$ と両端平均 $-\frac12[\varepsilon_i(0)+\varepsilon_i(1)]$ の誤差を，$a$ を使って表せ．

ヒント:(1) $\varepsilon_i(n)=2an+b$ は $n$ の1次関数なので，$\int_0^1\varepsilon_i\dd n$ は中点値に等しい．(3) $\varepsilon_i''=0$ なので，本文の誤差評価($-\frac{1}{12}\varepsilon_i''$ と $+\frac{1}{24}\varepsilon_i''$)からは両端平均も中点則もともに厳密になるはずである——実際に計算して確認せよ．一方Koopmans型は $I_i$ を $a$ だけ過小評価する($a\gt0$ は $E(n_i)$ が下に凸であることに対応し，この $a$ が軌道緩和エネルギーにあたる)．

演習10.3 ポテンシャルの定数シフト

外部ポテンシャルを $v(\rr)\to v(\rr)+C$($C$ は定数)と変えたとき，次を示せ．

(1) 自己無撞着な密度 $\rho$ は変わらず，固有値だけが $\varepsilon_i\to\varepsilon_i+C$ とずれる．

(2) 全エネルギーの表式 \eqref{eq:10-etot-final} を使うと $E\to E+NC$ となり，これは $\int(v+C)\rho=\int v\rho+NC$ という直接の評価と一致する．

ヒント:(1) 定数を足しても固有関数は変わらないので密度は不変．密度が不変なら $v_{\mathrm{H}}$ も $v_{xc}$ も不変で，自己無撞着性が保たれる．(2) 式 \eqref{eq:10-etot-final} の4項のうち変化するのは $\sum_i\varepsilon_i$ だけであり，$N$ 個の固有値がそれぞれ $C$ ずれる．この演習は，10.3.2節で「定数 $\mu_s-\mu$ は捨ててよい」と述べたことの検算にもなっている．

演習10.4 $\mathrm{H}_2^+$ の解離と非局在化誤差

水素原子1個の全エネルギーを，電子数 $\nu\in[0,1]$ の関数として $E_{\mathrm{at}}(\nu)$ と書く($E_{\mathrm{at}}(0)=0$，$E_{\mathrm{at}}(1)=-0.5$ Ha)．$\mathrm{H}_2^+$ を無限に引き離した極限では，2つの原子に合計1個の電子を $\nu$ と $1-\nu$ に分配する自由度がある．

(1) 厳密汎関数では $E_{\mathrm{at}}(\nu)$ が $\nu$ の1次関数(図10.3の折れ線)であることを使い，全エネルギー $E_{\mathrm{at}}(\nu)+E_{\mathrm{at}}(1-\nu)$ が $\nu$ によらないことを示せ．すなわち $(1,0)$ と $(0.5,0.5)$ は縮退している．

(2) 近似汎関数では $E_{\mathrm{at}}(\nu)$ が下に凸のなめらかな曲線になる．このとき $\nu=1/2$ の配置が $\nu=1$ の配置より低いエネルギーを持つことを，凸性の定義から示せ．

(3) (2) の結果が，LDA/GGAによる $\mathrm{H}_2^+$ の解離曲線の破綻(結合距離を伸ばしてもエネルギーが正しい解離極限に収束しない)をどう説明するか，200字程度で述べよ．

ヒント:(2) 凸関数の定義 $E(\lambda x_1+(1-\lambda)x_2)\le\lambda E(x_1)+(1-\lambda)E(x_2)$ を $x_1=1$, $x_2=0$, $\lambda=1/2$ に適用すると $E_{\mathrm{at}}(1/2)\le\frac12E_{\mathrm{at}}(1)$．これを2倍する．(3) 1電子系なのに誤りが出る点(自己相互作用誤差，10.1節の注意)にも触れるとよい．

10.8.3 参考文献

  1. P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Phys. Rev. 136, B864 (1964). ——DFTの出発点(第9章)．
  2. W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140, A1133 (1965). ——本章の主題．KS方程式の原論文．
  3. M. Levy, Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem, Proc. Natl. Acad. Sci. USA 76, 6062 (1979). ——制約付き探索(式 \eqref{eq:10-ts-levy} の流儀)．
  4. J. F. Janak, Proof that $\partial E/\partial n_i=\varepsilon_i$ in density-functional theory, Phys. Rev. B 18, 7165 (1978). ——10.7.3節の定理の原論文．
  5. J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Jr., Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy, Phys. Rev. Lett. 49, 1691 (1982). ——アンサンブルDFT，$E(\nu)$ の折れ線性，$\varepsilon_{\mathrm{HOMO}}=-I$．
  6. U. von Barth and L. Hedin, A local exchange-correlation potential for the spin polarized case, J. Phys. C: Solid State Phys. 5, 1629 (1972). ——スピン密度汎関数法(10.6節)．
  7. J. P. Perdew and M. Levy, Physical Content of the Exact Kohn-Sham Orbital Energies: Band Gaps and Derivative Discontinuities, Phys. Rev. Lett. 51, 1884 (1983).
  8. L. J. Sham and M. Schlüter, Density-Functional Theory of the Energy Gap, Phys. Rev. Lett. 51, 1888 (1983). ——[7]と同時発表．式 \eqref{eq:10-gap-key} の起源．
  9. J. C. Slater, Quantum Theory of Molecules and Solids, Vol. 4: The Self-Consistent Field for Molecules and Solids (McGraw-Hill, 1974). ——遷移状態法．
  10. C.-O. Almbladh and U. von Barth, Exact results for the charge and spin densities, exchange-correlation potentials, and density-functional eigenvalues, Phys. Rev. B 31, 3231 (1985). ——$v_{xc}\to-1/r$ の漸近形と $\varepsilon_{\mathrm{HOMO}}=-I$．
  11. D. E. Eastman, F. J. Himpsel, and J. A. Knapp, Experimental Exchange-Split Energy-Band Dispersions for Fe, Co, and Ni, Phys. Rev. Lett. 44, 95 (1980). ——3d遷移金属のARPESと計算の比較．
  12. R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, 1989). ——化学寄りの標準教科書．Janakの定理・遷移状態の記述が詳しい．
  13. R. M. Dreizler and E. K. U. Gross, Density Functional Theory: An Approach to the Quantum Many-Body Problem (Springer, 1990). ——数学的に厳密な記述．$v$-表示可能性の議論．
  14. R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004). ——固体物理寄りの標準教科書．本章の全範囲をカバーする．