密度汎関数理論

近似の数理　第9章 §9.7　★★★ 厳密な存在定理と、具体形が分からない汎関数を閉じる C 型の近似の境目

設定

$N$ 個の電子からなる系のハミルトニアンを、原子単位で

$$H_v=T+W+V,\qquad T=-\frac12\sum_{i=1}^N\Delta_i,\qquad W=\sum_{i<j}\frac{1}{\lvert x_i-x_j\rvert },\qquad V=\sum_{i=1}^Nv(x_i)$$

と書く（スピンは記号から省く）。$T$ と $W$ はどの原子・分子でも共通で、系ごとに違うのは原子核がつくる外部ポテンシャル $v$ だけである。反対称な波動関数 $\psi$ の 1 粒子密度を

$$\rho_\psi(x)=N\int\lvert \psi(x,x_2,\dots,x_N)\rvert ^2\,\mathrm{d} x_2\cdots\mathrm{d} x_N$$

とすると、$\langle \psi,\,V\psi\rangle =\int v\rho_\psi\,\mathrm{d} x$ となる。外部ポテンシャルのエネルギーは密度だけで決まる。第5章の Hartree–Fock 近似は波動関数を Slater 行列式に制限する A 型の近似だった。DFT は、$3N$ 変数の波動関数の代わりに 3 変数の密度 $\rho$ を基本変数にとる。

Hohenberg–Kohn の定理

証明

$v_1-v_2$ が定数でないと仮定して矛盾を導く。基底エネルギーを $E_1,E_2$ とする。

第1段（$\psi_1$ と $\psi_2$ は異なる）。もし $\psi_1=\psi_2=\psi$（位相を除いて）なら、二つの固有値方程式を引いて$\bigl(\sum_i(v_1-v_2)(x_i)-(E_1-E_2)\bigr)\psi=0$。$\psi$ が正の測度の集合上で $0$ にならない（ほとんど至る所で $0$ でない）なら、$\sum_i g(x_i)=E_1-E_2$（$g=v_1-v_2$）がほとんど至る所で成り立つ。$x_2,\dots,x_N$ を固定して $x_1$ だけ動かすと $g(x_1)$ が定数となり、仮定に反する。

第2段（変分原理）。$\psi_2$ は $H_{v_1}$ の基底状態ではない（$H_{v_1}$ の基底状態は縮退がなく $\psi_1$ だけで、第1段より $\psi_2\ne\psi_1$）。したがって変分原理（第3章）により厳密な不等号で

$$E_1<\langle \psi_2,\,H_{v_1}\psi_2\rangle =\langle \psi_2,\,H_{v_2}\psi_2\rangle +\int(v_1-v_2)\rho\,\mathrm{d} x=E_2+\int(v_1-v_2)\rho\,\mathrm{d} x .$$

1 と 2 の役割を入れ替えると $E_2<E_1+\int(v_2-v_1)\rho\,\mathrm{d} x$。二つを足すと $E_1+E_2<E_1+E_2$ となり矛盾する。

定理から、$\rho$ が $v$ を（定数を除いて）決め、$v$ が $H_v$ を決め、$H_v$ が基底状態 $\psi$ を決める。したがって基底状態のすべての性質は $\rho$ の汎関数であり、特に

$$F_{\mathrm{HK}}[\rho]=\langle \psi_\rho,\,(T+W)\psi_\rho\rangle$$

が定義でき、$E_v=\min_\rho\bigl\{F_{\mathrm{HK}}[\rho]+\int v\rho\bigr\}$ という密度についての変分原理（Hohenberg–Kohn の第2定理）が成り立つ。問題は $F_{\mathrm{HK}}$ の定義域である。$F_{\mathrm{HK}}$ は、ある $v$ の基底状態の密度として現れる $\rho$（$v$ 表現可能な密度）の上でしか定義されない。どの $\rho$ が $v$ 表現可能かは分かっておらず、表現可能でない密度が存在することも示されている。変分原理の探索範囲が分からないのでは、変分法として使えない。

Levy–Lieb の制約付き探索と Lieb の凸解析

証明

波動関数全体についての $\inf$ を、「まず密度を固定して $\psi$ を動かし、次に密度を動かす」二段に分ける。$\langle \psi,\,H_v\psi\rangle =\langle \psi,\,(T+W)\psi\rangle +\int v\rho_\psi$ の第2項は密度を固定すれば定数なので

$$\inf_\psi\langle \psi,\,H_v\psi\rangle =\inf_\rho\Bigl[\inf_{\psi:\rho_\psi=\rho}\langle \psi,\,(T+W)\psi\rangle +\int v\rho\Bigr]=\inf_\rho\Bigl[F_{\mathrm{LL}}[\rho]+\int v\rho\Bigr].$$

$\rho$ が $v$ の基底状態 $\psi_\rho$ の密度なら、$\rho_\psi=\rho$ を満たす任意の $\psi$ について$\langle \psi,\,(T+W)\psi\rangle +\int v\rho=\langle \psi,\,H_v\psi\rangle \ge\langle \psi_\rho,\,H_v\psi_\rho\rangle =F_{\mathrm{HK}}[\rho]+\int v\rho$。したがって $F_{\mathrm{LL}}[\rho]\ge F_{\mathrm{HK}}[\rho]$、$\psi=\psi_\rho$ で等号がとられる。

この定式化は $v$ 表現可能性の問題を回避する。探索範囲は $N$ 表現可能な密度で、これは簡単な条件で特徴付けられている。Hohenberg–Kohn の定理（一意性）は、この変分原理を使うためにはもはや必要ない。

Lieb はさらにこの構造を凸解析の言葉で整理した。主張だけを述べる。基底エネルギー $E_v$ を $v$ の関数と見ると、$E_v$ は $v$ についてのアフィン関数 $\langle \psi,\,H_v\psi\rangle$ の下限なので凹関数である。そこで Legendre–Fenchel 変換

$$F_{\mathrm L}[\rho]=\sup_v\Bigl\{E_v-\int v\rho\,\mathrm{d} x\Bigr\}$$

を定義すると、$F_{\mathrm L}$ は凸で下半連続な汎関数になり、$E_v=\inf_\rho\{F_{\mathrm L}[\rho]+\int v\rho\}$ が成り立つ（$\rho\in L^1\cap L^3$、$v\in L^{3/2}+L^\infty$ という対になる関数空間の上で）。$F_{\mathrm L}$ は $F_{\mathrm{LL}}$ の凸包（混合状態にまで探索を広げた汎関数）に一致する。$E_v$ と $F_{\mathrm L}$ の関係は、熱力学で内部エネルギーと自由エネルギーが Legendre 変換で結ばれるのと同じ構造である。凸性があるので、汎関数の性質（下半連続性、最小値の存在）を一般論から議論できる。ただし $F_{\mathrm L}$ の微分可能性は自明でなく、Kohn–Sham 方程式の導出で使う「汎関数微分」の存在は数学的にはなお微妙な点として残っている。

Kohn–Sham 方程式

$F_{\mathrm{LL}}$ は存在するが、その具体形は分からない。$N$ 電子問題を解かずに $F_{\mathrm{LL}}[\rho]$ を計算する方法がないからである。Kohn と Sham は、分からない部分をできるだけ小さくする書き換えを考えた。相互作用のない $N$ 電子系の運動エネルギーを

$$T_s[\rho]=\inf\bigl\{\langle \Phi,\,T\Phi\rangle \ :\ \Phi\ \text{は Slater 行列式で}\ \rho_\Phi=\rho\bigr\}$$

と定義し、古典的な静電エネルギー（Hartree エネルギー）を $J[\rho]=\frac12\iint\rho(x)\rho(y)/\lvert x-y\rvert \,\mathrm{d} x\,\mathrm{d} y$ とする。そして

$$E_{\mathrm{xc}}[\rho]:=F_{\mathrm{LL}}[\rho]-T_s[\rho]-J[\rho]\tag{9.39}$$

を交換相関汎関数と定義する。これは定義なので、ここまでは何も近似していない。$E_{\mathrm{xc}}$ には、Pauli 原理による交換エネルギー、電子が互いを避ける相関エネルギー、そして相互作用による運動エネルギーの増分 $\langle \psi,\,T\psi\rangle -T_s$ が押し込まれている。

エネルギー $T_s[\rho]+J[\rho]+E_{\mathrm{xc}}[\rho]+\int v\rho$ を最小にする $\rho$ が、ある Slater 行列式 $\Phi=\det(\phi_1,\dots,\phi_N)$ の密度$\rho=\sum_i\lvert \phi_i\rvert ^2$ として実現され（非相互作用 $v$ 表現可能性の仮定）、$E_{\mathrm{xc}}$ が汎関数微分 $v_{\mathrm{xc}}=\delta E_{\mathrm{xc}}/\delta\rho$ をもつと仮定する。軌道 $\phi_i$ について正規直交性の制約のもとで変分をとると（Lagrange 乗数 $\varepsilon _i$）

$$\Bigl(-\frac12\Delta+v(x)+v_{\mathrm H}(x)+v_{\mathrm{xc}}(x)\Bigr)\phi_i=\varepsilon _i\phi_i,\qquad v_{\mathrm H}(x)=\int\frac{\rho(y)}{\lvert x-y\rvert }\mathrm{d} y,\qquad\rho=\sum_{i=1}^N\lvert \phi_i\rvert ^2\tag{9.40}$$

が得られる。これが Kohn–Sham 方程式で、相互作用のない電子が有効ポテンシャル $v_s=v+v_{\mathrm H}+v_{\mathrm{xc}}$ の中を運動する 1 体の問題の形をしている。$v_s$ が $\rho$ に依存するので、Hartree–Fock と同じく自己無撞着に解く。

構造をまとめると、Kohn–Sham の枠組みは $E_{\mathrm{xc}}$ が厳密なら厳密であり、近似はすべて $E_{\mathrm{xc}}$ の選び方に集約される。運動エネルギーの大部分（$T_s$）を軌道で厳密に扱うので、Thomas–Fermi 理論（第6章）のように運動エネルギーを密度の局所式で近似したときの大きな誤差（原子の殻構造が出ない、分子が結合しない）は起きない。未知の部分を、全エネルギーの中で比較的小さな $E_{\mathrm{xc}}$ に押し込んだことが、この書き換えの要点である。

交換相関汎関数の近似：C 型の閉じ方

$E_{\mathrm{xc}}$ の厳密な形は存在が保証されるだけで分からない。実際の計算では近似形を使う。

• 局所密度近似（LDA）：$E_{\mathrm{xc}}^{\mathrm{LDA}}[\rho]=\int\rho(x)\,e_{\mathrm{xc}}\bigl(\rho(x)\bigr)\,\mathrm{d} x$。$e_{\mathrm{xc}}(n)$ は密度 $n$ の一様電子気体の 1 電子あたりの交換相関エネルギーである。交換部分は解析的に求まり $e_{\mathrm x}(n)=-C_{\mathrm x}n^{1/3}$、$C_{\mathrm x}=\frac34(3/\pi)^{1/3}\approx0.7386$（検算スクリプト）。相関部分は Ceperley と Alder の量子モンテカルロ計算を補間した式が使われる。
• 一般化勾配近似（GGA）：密度の勾配も入れて $E_{\mathrm x}^{\mathrm{GGA}}=-C_{\mathrm x}\int\rho^{4/3}F_{\mathrm x}(s)$ とする。$s\propto\lvert \nabla\rho\rvert /\rho^{4/3}$ は無次元の勾配、$F_{\mathrm x}$ は増大因子。
• 混成汎関数：Hartree–Fock 型の厳密交換エネルギーを一定の割合で混ぜる。B3LYP の 3 つの混合パラメータは、分子の熱化学データに合うように決められた。

どれも小さいパラメータの極限として導かれたものではなく、C 型の閉じ方である。それでも、正当化の完全な欠如ではない部分が二つある。

第一に、LDA はある極限で厳密である。 一様電子気体ではもちろん厳密で、密度がゆっくり変化する極限（密度の変化の長さが局所的な Fermi 波長に比べて長い）では、LDA が厳密な汎関数の主要項を与えることがLewin、Lieb、Seiringer によって数学的に証明されている。LDA は、B 型の極限で正しい式を、その極限から遠い原子・分子（密度が原子核の近くで急変し、外側で指数的に減衰する）に外挿して使っている。外挿がなぜ予想以上にうまくいくかについては、LDA の交換相関の「穴」が和則（穴に含まれる電子がちょうど 1 個であること）を満たしている、という説明がなされている。

第二に、厳密条件が正当化の代わりの制約として使われる。 厳密な $E_{\mathrm{xc}}$ が満たすべき性質のうち証明できるものがいくつかあり、近似汎関数をそれに合わせて作る。代表例は次の二つである。

• 座標の一様スケーリング：交換エネルギーは $\rho_\gamma(x)=\gamma^3\rho(\gamma x)$ について $E_{\mathrm x}[\rho_\gamma]=\gamma E_{\mathrm x}[\rho]$ を満たす（例題 9.8）。
• Lieb–Oxford の下界：電子間反発の期待値と古典的な静電エネルギーの差について$\langle \psi,\,W\psi\rangle -J[\rho_\psi]\ge-C_{\mathrm{LO}}\int\rho_\psi^{4/3}$ がすべての $\psi$ で成り立ち、$C_{\mathrm{LO}}\le1.68$（Lieb と Oxford の評価）。$E_{\mathrm{xc}}$ に対する下界を与える。

Perdew–Burke–Ernzerhof の GGA（PBE）は、経験的なパラメータを使わずに厳密条件だけから増大因子 $F_{\mathrm x}(s)=1+\kappa-\kappa/(1+\mu s^2/\kappa)$ を決めている。上限 $1+\kappa$ の $\kappa=0.804$ は、完全にスピン偏極した場合に局所的に Lieb–Oxford の界を満たす最大の値として選ばれている。実際、$2^{1/3}(1+\kappa)C_{\mathrm x}=2^{1/3}\times1.804\times0.7386\approx1.679$ となる（検算スクリプト）。

DFT の実際の評価は、厳密に解ける極限（一様電子気体、1 電子系、強い圧縮の極限など）での一致、厳密条件の充足、そして量子化学の高精度計算や実験値からなるベンチマーク集合での誤差統計によって行われる。LDA から GGA、混成汎関数へと進む系列は、パラメータを加えるごとにベンチマークの誤差を減らしてきたが、それは「ある極限で誤差が $0$ に行く」という形の保証ではない。特定の系（強相関電子系、分散力、自己相互作用誤差が目立つ系）では、どの近似汎関数も系統的に誤ることが知られている。