密度汎関数法
Hohenberg-Kohn の定理の記事で紹介したように、外部から加えられたポテンシャル中で相互作用しながら運動する電子系については、基底状態の密度 \(n_g(\b{r})\) が与えられれば、系のすべての性質が \(n_g(\b{r})\) から計算できる。そこで原理的には、基底エネルギーを \(n_g(\b{r})\) から計算する汎関数 \(E[n]\) に対して変分原理を適用すれば、基底状態の電子密度を変分的に求られる。これを利用して現実の物質の電子状態を計算するのが、密度汎関数法 (DFT) と呼ばれている計算手法だ。しかし DFT は \(E[n]\) の正しい形が全くわからないという問題を抱えている。
そこで \(E[n]\) の近似式が必要となる。現在一番使われているのが、Kohn-Sham 型の汎関数である。今回はこの Kohn-Sham 理論について解説する。
Kohn-Sham の理論
DFT のエネルギー汎関数 \(E[n]\) を構成しようとするときには、物理的直感に基づいて
\begin{align}
E[n] = T[n] + E_{xc}[n] + \int U(\b{r})n(\b{r}) d\b{r} + \frac{e^2}{4\pi\epsilon_0}\int \frac{n(\b{r})n'(\b{r})}{|\b{r}-\b{r}'|} d\b{r}d\b{r}'
\end{align}
とその構成要素に分割するのが常套手段だ。上式における各項は、
- \(T[n]\): 運動エネルギー汎関数
- \(E_{\mathrm{xc}}[n]\): 交換・相関エネルギー汎関数
- \(\int U(\b{r})n(\b{r}) d\b{r}\): 外部ポテンシャルによるポテンシャルエネルギー
- \(\frac{e^2}{4\pi\epsilon_0}\int \frac{n(\b{r})n'(\b{r})}{|\b{r}-\b{r}'|} d\b{r}d\b{r}'\): 電子間のクーロン力による (古典的な) ポテンシャルエネルギー
である。
前回は Thomas-Fermiの理論として、運動エネルギー汎関数 \(T[n]\) を一様電子の運動エネルギーによって近似する方法について説明した。しかし実用上は、Thomas-Fermiの運動エネルギー汎関数はそこまで良い計算結果を与えないとされている。
運動エネルギーは電子のエネルギーの大部分を占めるため、ここの精度を落とすのがまずいらしい。
そこで Kohn と Sham は、DFT を定式化する際に捨てられた
軌道という概念を復活させることを思いついた。Kohn-Sham の理論では、電子密度 \(n(\b{r})\) が、\(N\)電子の波動関数 \(\psi(\b{r}_1,\b{r}_2,\cdots,\b{r}_N)\)、
\begin{align}
\psi(\b{r}_1,\b{r}_2,\cdots,\b{r}_N)
= \phi_1(\b{r}_1) \phi_2(\b{r}_2) \cdots \phi_N(\b{r}_N) \tag{1}
\end{align}
によって導かれるものであることを
仮定する。ただし \(\phi_i(\b{r})\) は1電子の波動関数としての資格を持てるように、
\begin{align}
\int \phi_i(\b{r})\phi_j(\b{r}) d\b{r} = \delta_{ij}
\end{align}
を満たすように選ぶ。これが Kohn-Sham 理論の一番の肝である。
この波動関数は \(N\) 個の電子が全く独立に \(\phi_i(\b{r})\) という軌道で運動している状況を表しているので、その運動エネルギーは
\begin{align}
T_{\mathrm{KS}}[n_{\b{\phi}}] = -\sum_{i=1}^N \frac{\hbar^2}{2m} \int \phi_i^*(\b{r})\nabla^2 \phi_i(\b{r}) \b{r} \tag{1}
\end{align}
である。Kohn-Sham 理論では電子密度 \(n\) が 1 電子波動関数の組 \(\b{\phi}:=\{\phi_1, \phi_2, \cdots, \phi_N\}\) によって決まることを明示的に表すために、\(n_{\b{\phi}}\) と書いた。これは波動関数が (1) 式で表されているという仮定のもとで厳密な運動エネルギーを与えるので、一様な電子ガスから近似的に導いた Thomas-Fermi の汎関数よりも、「良い」運動エネルギー汎関数になっているといえるだろう。
\(\int U(\b{r})n(\b{r}) d\b{r}\) や \(\frac{e^2}{4\pi\epsilon_0}\int \frac{n(\b{r})n'(\b{r})}{|\b{r}-\b{r}'|} d\b{r}d\b{r}'\) で表されるポテンシャルエネルギーを計算するには、電子密度 \(n(\b{r})\) も必要となるので、ここでその表式を書いておこう。
これも運動エネルギーと同じように、電子が独立に運動しているという描像から、
\[
n_{\b{\phi}}(\b{r}) = \sum_{i=1}^N \left|\phi_i(\b{r})\right|^2 \tag{3}
\]
が得られる。(密度が本当にこれでいいのか不安に感じる人は、
電子密度の記事へどうぞ。)
Kohn-Sham 理論のエネルギー汎関数
ここまでで、エネルギー汎関数のうち、Kohn-Sham理論において運動エネルギーとポテンシャルエネルギーを計算する方法がわかったと思う。
あとは交換相関エネルギー \(E_{\mathrm{xc}}[n]\) という項を計算する方法を示せば、Kohn-Sham 理論で最小化するべきエネルギー汎関数が計算できることになる。
しかしこのページではこの汎関数 \(E_{\mathrm{xc}}[n]\) の具体的な形を導出することは避け、大まかな理論を追うことに注力したい。
というのも、\(E_{\mathrm{xc}}[n]\) は人が物理的直感に基づいて近似的な形を与えるしかない部分で、数え切れないほどのバリエーションが提案されているからだ。
そのどれもが、導出しようとすればそれなりの計算を要するようなので、Kohn-Sham 理論を紹介するというこの記事の趣旨からはずれてしまう。
僕がそれらの導出をまだしっかりと勉強しきれていないのもある。
ともかく、交換相関エネルギー汎関数 \(E_{\mathrm{xc}}[n]\) は (なにかとても複雑な議論によって) 与えられているものとしてしまおう。そうすると、
Kohn-Sham 理論とは、
\begin{align}
E[n_{\b{\phi}}] = T_{KS}[n_{\b{\phi}}] + E_{xc}[n_{\b{\phi}}] + \int U(\b{r})n_{\b{\phi}}(\b{r}) d\b{r} + \frac{e^2}{4\pi\epsilon_0}\int \frac{n_{\b{\phi}}(\b{r})n_{\b{\phi}}'(\b{r})}{|\b{r}-\b{r}'|} d\b{r}d\b{r}' \tag{4}
\end{align}
というエネルギー汎関数を、(2) 式の運動エネルギー \(\sum_{i=1}^N \frac{\hbar^2}{2m} \int \phi_i^*(\b{r}) \nabla^2 \phi_i(\b{r}) d\b{r}\)、(3) 式の電子密度 \(n_{\b{\phi}}(\b{r}) = \sum_{i=1}^N \left|\phi_i(\b{r})\right|^2 \) によって計算し、それを最小化するような関数の組 \(\{\phi_1, \phi_2, \cdots, \phi_N\}\) を見つけ出すことによって、基底状態のエネルギーを近似しようという試みである。
といえる。そしてこの最小化問題を効率的に解くための方程式が、
Kohn-Sham方程式
である。Kohn-Sham方程式は最後のセクションで導出するが、その前に、Kohn-Sham理論について注意するべき点を書こうと思う。
少しの注意
Kohn-Sham 理論では、電子密度が \(\psi(\b{r}_1,\b{r}_2,\cdots,\b{r}_N)
= \phi_1(\b{r}_1) \phi_2(\b{r}_2) \cdots \phi_N(\b{r}_N) \) という波動関数から導かれる形、
\[
n_{\b{\phi}}(\b{r}) = \sum_{i=1}^N \left|\phi_i(\b{r})\right|^2 \tag{3}
\]
であることを仮定する。
一見すると、背後に \(\psi(\b{r}_1,\b{r}_2,\cdots,\b{r}_N)
= \phi_1(\b{r}_1) \phi_2(\b{r}_2) \cdots \phi_N(\b{r}_N) \) という波動関数が実際に実現していて、その波動関数の電子密度として、\(n_{\b{\phi}}(\b{r})\) が得られているのだ、と理解するかもしれない。
しかしKohn-Sham理論が仮定することは、言葉どおり「\(\psi(\b{r}_1,\b{r}_2,\cdots,\b{r}_N)
= \phi_1(\b{r}_1) \phi_2(\b{r}_2) \cdots \phi_N(\b{r}_N) \) から"導かれる"電子密度」で「基底状態の電子密度」が記述できるということである。
実際の基底状態の波動関数 \(\Psi_g(\b{r}_1,\cdots,\b{r}_N)\) が、\(\phi_1(\b{r}_1) \phi_2(\b{r}_2) \cdots \phi_N(\b{r}_N) \) になるとは仮定していない。あくまで基底状態の電子密度が、\(n_{\b{\phi}}(\b{r})\) によって書き表せる、ということだけを仮定している。
そしてこの仮定に対する理論的裏付けは (僕が調べた限りでは) まだなされていないようだ。
そんなこんなで理論的背景には少しの不安はあるものの、実際に Kohn-Sham 理論を使って計算してみると、様々な物理系で実験値をよく再現するらしい。Kohn-Sham 理論に基づいた DFT 計算は、そのようなある種の経験則を元に広く使われているみたいだ。
もちろん Kohn-Sham の理論的裏付けを頑張って研究している人もたくさんいるはず。
Kohn-Sham 方程式の導出
目標は、関数 \(\phi_1,\cdots,\phi_N\) によって値が変わるエネルギー汎関数
\begin{align}
E[n_{\b{\phi}}] = T_{KS}[n_{\b{\phi}}] + E_{xc}[n_{\b{\phi}}] + \int U(\b{r})n_{\b{\phi}}(\b{r}) d\b{r} + \frac{e^2}{4\pi\epsilon_0}\int \frac{n_{\b{\phi}}(\b{r})n_{\b{\phi}}(\b{r}')}{|\b{r}-\b{r}'|} d\b{r}d\b{r}' \tag{4}
\end{align}
を最小化することだ。ただし \(\phi\) は何を使ってもいいわけではなく、
\begin{align}
\int \phi_i(\b{r})\phi_j(\b{r}) d\b{r} = \delta_{ij}
\end{align}
となっていてほしい。
\(E\) をもう少し書き下してみると、
\begin{align}
E[n_{\b{\phi}}] = \sum_{i=1}^N \int \phi_i^*(\b{r})\left(-\frac{\hbar^2\nabla^2}{2m} + U(\b{r})\right)\phi_i(\b{r}) d\b{r} + \sum_{i,j} \frac{e^2}{4\pi\epsilon_0}\int \frac{\phi_i^*(\b{r})\phi_i(\b{r})\phi_j^*(\b{r}')\phi_j(\b{r}')}{|\b{r}-\b{r}'|} d\b{r}d\b{r}' + E_{xc}[n_{\b{\phi}}]
\end{align}
となる。\(\int \phi_i(\b{r})\phi_j(\b{r}) d\b{r} = \delta_{ij}\) という拘束条件下におけるこの形の関数の最小化は
Hartree-Fock方程式の導出でも出てきた。(というか全く同じ形だ。) したがって
Hartree-Fock方程式の導出の記事 と全く同じように、関数 \(\phi_i^*\) による「微分」を考えてやることによって、以下が得られる。
Kohn-Sham 方程式:
\begin{align}
\left[-\frac{\hbar^2\nabla^2}{2m} + U(\b{r}) + \sum_{j}\int \frac{|\phi_{j}(\b{r}')|^2}{|\b{r}-\b{r}'|} d\b{r}' + \frac{\delta E_{xc}}{\delta n}\right] \phi_i(\b{r}) = \epsilon_i \phi_i(\b{r})
\end{align}
ただし、\(E_{xc}[n_{\b{\phi}}]\) の \(\phi_i^*\) による微分が、
\[
\frac{\delta E_{xc}[n_{\b{\phi}}]}{\delta \phi_i^*} = \frac{\delta E_{xc}[n_{\b{\phi}}]}{\delta n} \frac{\delta n}{\delta \phi_i^*} = \frac{\delta E_{xc}[n_{\b{\phi}}]}{\delta n} \phi_i
\]
と書けることを使った。
もちろん、\(E_{xc}[n_{\b{\phi}}]\) が運動エネルギーと同じように \(\{\phi_i\}\) を明示的に使うものである場合には、上の表式は少し異なるものになる。
最後に一つくらいは \(E_{xc}[n_{\b{\phi}}]\) の形を示して、具体的な Kohn-Sham 方程式を書いておこうと思う。局所密度近似 (local density approximation, LDA) と呼ばれる種類の交換相関汎関数は
\begin{align}
E_{xc, \mathrm{LDA}}[n_{\b{\phi}}] = \int \epsilon_{xc}(n_{\b{\phi}}(\b{r})) n_{\b{\phi}}(\b{r}) d\b{r}
\end{align}
の形を取る。ここで \(\epsilon_{xc}\) は密度 \(n\) の一様電子ガスが持つ交換エネルギー密度と相関エネルギー密度を足した量である。
\(\epsilon_{xc}\)の具体的な関数形をここで書くことは避ける。というのも、実は一様電子ガスですら相関エネルギーの解析解は得られていないのだ。そのため数値計算に基づいた近似式が色々提案されていて、それら近似式について個別の説明をすると話が逸れてしまう。通常は、\(\epsilon_{xc}(n)\) は \(n\) に関する初等関数 (多項式, log など) で表されているので、Kohn-Sham を理解するにあたっては、背後にはそういうものがあるのだなあ、と思っておけば良い。
このとき
\begin{align}
\frac{\delta E_{xc, \mathrm{LDA}}[n_{\b{\phi}}]}{\delta n}
= \epsilon_{xc}(n_{\b{\phi}}(\b{r})) + n_{\b{\phi}}(\b{r})\left.\frac{\partial \epsilon_{xc}}{\partial n}\right|_{n=n_{\phi}(\b{r})}
\end{align}
すると Kohn-Sham 方程式は
LDA-Kohn-Sham 方程式:
\begin{align}
\left[-\frac{\hbar^2\nabla^2}{2m} + U(\b{r}) + \sum_{j}\int \frac{|\phi_{j}(\b{r}')|^2}{|\b{r}-\b{r}'|} d\b{r}' + \epsilon_{xc}(n_{\b{\phi}}(\b{r})) + n_{\b{\phi}}(\b{r})\left.\frac{\partial \epsilon_{xc}}{\partial n}\right|_{n=n_{\phi}(\b{r})}\right] \phi_i(\b{r}) = \epsilon_i \phi_i(\b{r})
\end{align}
となる。
\(i=1,2,\cdots,N\) (\(N\)は電子数) だったので、Kohn-Sham 方程式を解くときには \(N\) 本の連立微分・積分方程式を解くことになる。
実際に解く際には、もちろん通常の微分方程式の数値解法を使うこともできるが、普通は \(\phi_i\) を平面波など適当な基底関数の和に展開して、行列の固有値問題に落としてしまうのが常套手段である。
その手続は Hartree-Fock の数値解法である
Roothaan 方程式 とほぼ同じなので、ここでは紹介しない。
将来気が向いたら書くかもしれないけど。