密度汎関数法入門(1)-理論

はじめに
量子多体系の計算の難しさ
概要
Hohenberg-Konhの定理
Kohn-Sham equation

はじめに

先日の五月祭で密度汎関数法(Density Functional Theory, DFT)に関する展示をしたのですが、我ながらなかなかの力作に仕上がったので、ここで解説記事を書くことにしました。

密度汎関数法は、相互作用する量子多体系の問題を、その基底状態の粒子密度を再現するような相互作用のない問題に書き換えることで解く方法で、化学や固体物理における電子状態計算などによく用いられる手法です。

このシリーズ記事では、

最初にDFTの基礎理論と、数値計算の手法の簡略な紹介
現実の系への応用例として、原子の電子軌道の計算の説明

について書きます。特に、この記事では最初の基礎理論の紹介を行います。

なお、以下では電子の状態を考えることを仮定して話をするので、他の量子多体系を考えたい場合はCoulomb相互作用を仮定している部分を一般の相互作用に書き換えて議論を追ってみてください。

また、この記事を書く上で

Kohn, W. (1999). Nobel Lecture: Electronic structure of matter—wave functions and density functionals. Reviews of Modern Physics, 71(5), 1253.
常田貴夫著「密度汎関数法の基礎」

などを参考にしました。

続きは

密度汎関数法入門(2)-計算方法 - Airy's blog

密度汎関数法入門(3)-計算例(実装) - Airy's blog

密度汎関数法入門(4)-計算例(結果) - Airy's blog

です。

量子多体系の計算の難しさ

DFTについて紹介する前に、なぜDFTが必要なのかということを説明します。

例えば、古典系においてNewtonの運動方程式に従って運動するN粒子の問題を考えると、この場合は各時間におけるN個の粒子の位置と運動量を保持しておいて、それをNewtonの運動方程式に従って時間発展させれば良いことがわかります。

それに対し、量子力学でN粒子系の状態を決めるにはN体波動関数 $\Psi(r_1, r_2, \cdots, \r_N)$ の関数形を決めてやる必要があります。数値的に関数形を保持するには、空間内のすべての位置に対応する関数の値を保持しておくことはできないので、空間内の有限の領域を離散化して、各格子点での関数の値を覚えておくことになります。ここで、 $x,y,z$ 軸をそれぞれ[n]分割したとすると、一つの位置がとりうる格子点の個数は $n^3$ であり、N粒子波動関数の引数の組み合わせは $n^{3N}$ 通りになることがわかります。これは、N体波動関数の自由度がNの指数関数で発散することを示しており、したがって波動関数を直接扱う方程式であるSchroedinger方程式は数値計算には向かない方程式であるということになります。このように、基礎方程式を直接数値的に計算できず、何らかの工夫をしなければならないということが古典と量子の場合の大きな違いになります。

いかに説明するように、DFTはこの問題を回避して問題を解ける形にする工夫の一つです。

概要

DFTの基本的な理論は、以下の二つの定理からなります。

定理1(Hohenberg-Kohnの定理)

量子多体系の全ての情報は、基底状態の粒子密度に含まれる。

定理2(Kohn-Sham equation)

Coulomb相互作用で相互作用する電子の基底状態の電子密度は、相互作用のない1粒子Schroedinger方程式

$\left(-\frac{1}{2}\nabla^2 + \int dr' \frac{n(r')}{|r-r'|} + V_{ext}(r) + V_{xc}(r) \right) \psi_j(r) = \epsilon_j \psi_j(r)$

を解いて得られる電子密度に等しい。

以下では、それぞれの定理について証明を含めてより詳しく説明します。

Hohenberg-Konhの定理

Hohenberg-Kohnの定理の正確なステートメントは以下の通りです。

Hohenberg-Kohnの定理

あるポテンシャルの下での電子ガスの基底状態の電子密度がわかっているとき、基底状態の電子密度がそれと等しくなるようなポテンシャルは定数のシフトを除いて一位に定まる。

証明

今、考えている基底状態の電子密度を $n(r)$ とし、基底状態の電子密度が $n(r)$ となるような二つの異なるポテンシャル $V_1(r), V_2(r)$ があったとする。

ここで、Hamiltonianのうち、外場ポテンシャル以外の項(運動エネルギーと電子館相互作用)を $H_0$ とおき、

$H_1 = H_0 + V_1$

$H_2 = H_0 + V_2$

とする。 $H_1, H_2$ の基底状態をそれぞれ $\Psi_1, \Psi_2$ とおくと、

$( \Psi_1, H_1\Psi_1) \lt (\Psi_2, H_1\Psi_2)$

より

$(\Psi_1, H_0\Psi_1) + \int dr n(r)V_1(r) \lt (\Psi_2, H_0\Psi_2) + \int dr n(r)V_1(r)$

同様に

$( \Psi_2, H_2\Psi_2) \lt (\Psi_1, H_2\Psi_1)$

より

$(\Psi_2, H_0\Psi_2) + \int dr n(r)V_2(r) \lt (\Psi_1, H_0\Psi_1) + \int dr n(r)V_2(r)$

辺々足して

$(\Psi_2, H_0\Psi_2) + (\Psi_1, H_0\Psi_1) \lt (\Psi_2, H_0\Psi_2) + (\Psi_1, H_0\Psi_1)$

となるが、これは矛盾である。したがって、題意が示された。

ただし、証明の途中で基底状態が縮退していないことを仮定した。

Kohn-Sham equation

Coulomb相互作用で相互作用する電子の基底状態の電子密度は、以下のKohn-Sham equation

$\left(-\frac{1}{2}\nabla^2 + \int dr' \frac{n(r')}{|r-r'|} + V_{ext}(r) + V_{xc}(r) \right) \psi_j(r) = \epsilon_j \psi_j(r)$

を解いて得られる電子密度に等しい。ただし、未定義の量 $V_{xc}(r)$ については以下の証明の過程で定義する。

証明

最初に

$T_s[\tilde{n}] = \min_{\Psi \to \tilde{n}} [ (\Psi| T |\Psi)]$

とする。ここで、 $T$ は運動量演算子、 $\min_{\Psi \to \tilde{n}}$ は、電子密度が[\tilde{n}]となるような状態の上でminを取るという意味である。これは、先ほどのHohenberg-Kohnの定理の議論をnon-interatcing electron gasについて行うと $T_s$ が電子密度 $\tilde{n}$ の汎関数として定義可能であることがわかる。

次に、

$F[\tilde{n}] = \min_{\Psi \to \tilde{n}} [ (\Psi| H |\Psi)] = T_s[\tilde{n}] + \frac{1}{2}\int \frac{\tilde{n}(r)\tilde{n}(r')}{|r-r'|} dr dr' + E_{xc}[\tilde{n}]$

として、 $F[\tilde{n}]$ とexchange correlation functional $E_{xc}[\tilde{n}]$ を定義する。ここで、 $T_s、F$ が外場ポテンシャルの形に依存しないuniversalな汎関数であるということに注意してください。

このとき、Hamiltonian Hの基底状態エネルギーは

$E_{gnd} = min_\tilde{n} (F[\tilde{n}] + \int \tilde{n}(r)V_{ext}(r) dr ) = min_\tilde{n} (T_s[\tilde{n}] + \frac{1}{2}\int \frac{\tilde{n}(r)\tilde{n}(r')}{|r-r'|} dr dr' + E_{xc}[\tilde{n}]+ \int \tilde{n}(r)V_{ext}(r) dr)$

である。基底状態の電子密度を $n(r)$ とすると、変分法より

$\left (\frac{\delta T_s[\tilde{n}] }{\delta \tilde{n}(r)} + \int \frac{\tilde{n}(r')}{|r-r'|} dr' +\frac{\delta E_{xc}[\tilde{n}] }{\delta \tilde{n}(r)} + V_{ext}(r) \right)_{\tilde{n} = n}=0$

を得る。

ここで一度、相互作用のない電子ガスを考えると、

$\frac{\delta T_s[\tilde{n}] }{\delta \tilde{n}(r)} + V_{ext}(r) = 0$

を満たすような電子密度を再現するには、外場の中で運動する一粒子のSchroedinger方程式

$(-\frac{1}{2}\nabla^2 + V_{ext}(r))\psi_j(r) = \epsilon_j \psi_j(r)$

を解けば良いことがわかっている。

したがって、相互作用のある場合では、基底状態の電子密度を正しく再現するには、実効的な一粒子の方程式

$(-\frac{1}{2}\nabla^2 + \int \frac{n(r')}{|r-r'|}dr' + V_{ext}(r) + V_{xc}(r))\psi_j(r) = \epsilon_j \psi_j(r)$

を解けば良いことがわかる。ただし、exchange correlation potential $V_{xc}(r)$ は

$V_{xc}(r) = \frac{\delta E_{xc}[\tilde{n}]}{\delta \tilde{n}(r)}$

として定義される。

ここで得られた一粒子の方程式をKohn-Sham equationとよびます。

ここで、相互作用するN粒子の問題が相互作用しないN粒子の問題に書き換えられました。相互作用しないN粒子の問題は、実質的にはN個の一粒子問題と考えることができるので、問題を解ける形に落とし込むことができました。

また、DFTが優れている点は、上で示した二つの定理が厳密であるという点です。DFTではN体波動関数を扱わないので系の情報の一部を捨てていることになりますが、Hartree-Fock法などとは異なり、系の状態空間を制限するというタイプの近似が使われていません。したがって、exchange correlation functional $E_{xc}$ の形を良くすれば、原理的にはいくらでも厳密解に近い結果が得られることになります。