Kohn-Sham 方程式

Kohn-Sham方程式と基底状態密度

このようにHohenberg-Kohnの定理のみで基底状態を議論することは不可能に近 い。そのため、KohnとShamは多体の方程式そのままを解こうとするのではなく、

相互作用しない1電子問題に帰着させ、計算しやすくするアイデアを与えた[29]。

この手法を用いて、外部ポテンシャルに対応する正しい基底状態密度を与える計 算を行うことができる。また、応用の幅が広がるためこの手法は多くの研究者に よって用いられている。これはKohn-Sham方程式と呼ばれ、以下の２つの仮定の もとに作られている。

1. 厳密な基底状態の密度は相互作用のない基底状態の密度から表すことがで きる。

2. このハミルトニアンは運動エネルギーと有効ポテンシャルから成り立つ。

1番目の仮定については厳密な証明がなされていないが、多くの計算結果が実験結 果と整合した結果を与えるため、その正当性があると仮定されている。

N電子のKohn-Sham方程式は以下の式で与えられる。

HKSψi = ϵiψi (2.17)

n₀(r) =

∑N i

|ψi(r)|² (2.18)

Kohn-Sham方程式の固有関数の絶対値2乗の和が、正しい基底状態の密度n₀(r)に なると仮定している。このときのハミルトニアンHKS は以下で表される。

H_KS = −1

2∇²+V(r) (2.19)

ポテンシャルエネルギーVは今の段階でその具体的形式は与えられていない。ま た、運動エネルギーは

T_s= −1 2

∑N i=1

⟨ψi|∇²|ψi⟩ (2.20)

で与えられる。ここで、古典的な電子間クーロン相互作用を表すHartreeエネルギー E_H =

∫ n(r)n(r^′)

|r−r^′| drdr^′ (2.21) を導入し、Hohenberg-Kohnで与えられたエネルギー汎関数の形式を以下のように 書き換える。

E_KS[n] = T_s[n]+

∫

V_ext(r)n(r)dr+E_H[n]+E_xc[n]

= T_s[n]+E_ext[n]+E_H[n]+E_xc[n] (2.22) E_xcは交換相関汎関数と呼ばれる項で、Eext[n]は外部ポテンシャルによるエネル ギー項ある。この交換相関汎関数は、式(2.13)との対応を考えると、

E_xc[n]= T [n]−T_s[n]+E_int[n]−E_H[n] (2.23) で与えられる。これは、「多体相互作用した」真の運動エネルギーTと電子間相互 作用エネルギーE_intから「相互作用のない」運動エネルギーT_sとHartreeエネル ギーEHの差であり、多体相互作用から相互作用していないエネルギーのずれを表 している。電子の交換エネルギーと相関エネルギーがその部分に対応するため、こ の項は交換相関エネルギーと呼ばれる。

以上の式とKohn-Sham 方程式から得られる固有関数の規格化条件から、Kohn-Sham方程式(2.19)のVが決定される。Kohn-Sham方程式の固有関数は変分で停 留するために、

δ δψ^∗_i(r)



E_KS −∑

j

ϵj{⟨ψj|ψj⟩}





= δTs

δψ^∗_i +[δEext

δn + δEH

δn + δExc

δn ] δn

δψ^∗_i −ϵiψi =0 (2.24)

2.2. 第一原理バンド計算 15

を満たす。これらの汎関数微分は、

δTs

δψ^∗_i = −1

2∇²ψi (2.25)

δE_ext

δn = V_ext (2.26)

δE_H

δn = V_H =

∫ n(r)

|r−r^′|dr^′ (2.27) δE_xc

δn = Vxc (2.28)

δn

δψ^∗_i = ψi (2.29)

となることから、以上よりKohn-Sham方程式

H_KSψi = ϵiψi (2.30)

H_KS = −1

2∇²+V (2.31)

V = V_H+V_xc+V_ext (2.32)

を得る。Kohn-Sham方程式は1電子方程式のため、ほとんどの問題を現実的な時 間内で解くことができる。ポテンシャル項はV_xc以外既知となるため、Kohn-Sham ハミルトニアンにあるV_xcの形式がわかると方程式を解くことが出来る。ただし、

この交換相関項は多体相互作用から生まれる項に該当する。そのため、その関数 表記を特定するということはF_HK[n]の関数形を定めるのと同様の難しさがある。

厳密なKohn-Sham方程式を解くことは現在でも重要な問題となっているが、近似

的な関数形はすでに与えられ、それが今日第一原理バンド計算で用いられている。

交換相関汎関数

交換相関汎関数の形は古典的な長距離力を表すHartree項と分離できているため、

極めて局所的な汎関数として近似できる。これを式で表すと、

Exc =

∫

drn(r)ϵxc([n],r) (2.33) となる。ϵxc([n],r)は位置rにある1電子あたりのエネルギーである。交換相関汎 関数の形式を定めるためには、ϵxcを決めることになる。この一番簡単な形式は、

E_xc^LDA=

∫

drn(r)ϵxc^hom(n (r)) (2.34)

Ψ

_i

({ r }) n

₀

( r )

V

ext

H

Hohenberg-Kohn

ψ

_i

( r ) n

₀

( r )

V H

_KS

Kohn-Sham

図2.2: Kohn-Sham方程式とHohenberg-Kohnの定理との対応、Hohenberg-Kohnの

定理(左部分)から、同じ密度n(r)を持つような有効ハミルトニアンのKohn-Sham

方程式(右部分)に問題を置き換えることによって、Hoenberg-Kohnで計算できな

かったF_HKの部分を解消する。

となり、ϵxcがある位置の密度のみの関数で表される。これを局所密度近似(Local Density Approximation : LDA)という。LDAの交換相関エネルギー密度ϵxc^hom(n (r)) は一様電子ガスでは正しい形式となる。一様電子ガスの交換相関エネルギーは量 子モンテカルロ法を用いることで計算できる。この一様電子ガスで正しい表現が 実際の物質で現実に即した結果を与えることが、様々な物質の計算結果から知ら れている。しかし、なぜこのような簡単な形式でよく近似できるのか不思議に思 われる。局所密度近似が現実に即した結果を与えるのか説明を行う。

交換相関汎関数の式(2.34)を書きなおす。密度を一定に保ちながら、電荷をゼ ロから実際の値(Hartree単位系では1)まで仮想的に変化させることを考える。こ のときのエネルギーの変化は、運動エネルギーは密度が変化しないので一定にな るため、交換相関汎関数は電子相互作用のエネルギーからHartreeエネルギーの差 分になり、

Exc[n] =

∫ ₁

0

dλ⟨Ψλ|dVint

dλ |Ψλ⟩ −EH[n]

=

∫

drn(r)

∫

dr^′¯n_xc(r,r^′)

|r−r^′| (2.35)

と表せ、2行目の変換は形式的に書き換えている。ここで¯n_xc(r,r^′)は結合定数につ いて平均された正孔密度で、

¯n_xc(r,r^′)=

∫ 1 0

dλn^λ_xc(r,r^′) (2.36)

2.2. 第一原理バンド計算 17

で与えられる。n^λ_xcを交換相関正孔密度という。λの積分を ¯nにすることで押し付 けた形になっている。この式(2.35)で重要な部分は平均化された正孔密度の積分 が全積分になるため、角度平均が重要なパラメータになることである。この交換 相関正孔密度は、局所密度近似で与えられた場合と正確な場合では、「角度平均」

がよく一致する。例として、Ne原子の交換正孔[30]とSi結晶の交換相関正孔の 結果[31]がある。厳密解（または量子モンテカルロ法）とLDAにおける結果と比 較すると、NeとSiどちらの結果も角度平均で見るとLDAと厳密解は近い結果を 与える。このため、局所密度近似は一様な電子ガスだけではなく、現実の結晶構 造に対してもよい結果を与える。

以上から局所密度近似が現実の結果をよく近似することがわかったが、さらに 精度をあげるために、今まで位置交換相関エネルギー密度が密度n(r)の関数でし かなかったところを密度の傾きも考慮した計算を入れる方法がある。これを一般 化勾配近似(General Gradient Approximation : GGA)という。GGAはより現実の結 果に即するため、現在の多くの計算ではこの手法が使われている。局所密度近似 は密度の位置のみの関数であるため、表現方法はひと通りしかなかったが、密度 の傾きを考慮に入れると、表現方法は多くなる。そのため、様々な関数形式が提 案されている。GGAの形式は

E^GGA_xc =

∫

drn(r)ϵxc(n,|∇n|, ...) (2.37) となるが、本論文では、Perdew, Burke, Ernzerhofによる形式(PBE)[32]を用いる。

これらを交換エネルギーE^PBE_x ^−GGAと相関エネルギーE_c^PBE−GGAに分解して、

E_x^PBE−GGA =

∫

drn(r)ϵ_x^hom(n)F_x(s) (2.38) E_c^PBE−GGA =

∫

drn(r)(ϵc^hom(n)+H) (2.39) となる。ϵ_x(c)^homは電子ガスの交換(相関)エネルギーを表す。また、Fxは無次元量で あり、s= _2k^∇n_F,k_F = (⁹₄^π)^1/3r⁻¹_s 、rsは電子間の平均距離を表す。このF_xは

F_x(s)=1+κ− κ

1+µs²/κ (2.40)

と表される。κ = 0.804はLieb-Oxfordの制約を満たすように決められている。ま

た、µ=0.21951とすることで、局所近似の線形応答の表式になるようになってい

る。相関項の部分Hは H = e²

a₀γϕ³ (

1+ β

γt² 1+At² 1+At²+A²t⁴

)

(2.41) となる。ここでeは電気素量、a0はボーア半径、β=0.066725、γ= (1−ln2)/π²は 無次元定数、ϕ= ((1+ζ)^2/3+(1−ζ)^2/3)/2であり、ζ= (n^↑−n^↓)/nはスピン分極、t は無次元化された勾配t =|∇n|/(2ϕkTFn)、kTFはトーマスフェルミ波数。Aは

A= β γ



exp



−ϵc^hom

γϕ^{3 e}_a²₀



−1





−1

(2.42) から求まる。

Kohn-Sham方程式と軌道エネルギー

これまで、シュレディンガー方程式に似たKohn-Sham方程式(2.30)について説 明してきたが、この式は多体問題を簡略化させるために出てきた式であるため、こ の固有値には数学的以上の何の意味もない。しかし、Kohn-Sham方程式のエネル ギー固有値は、第一原理バンド計算において準粒子エネルギーと同じ扱いをする。

数学的な価値でしかないKohn-Sham方程式のエネルギー固有値に対して準粒子エ ネルギーと扱うことができる理由について説明する。そのために、軌道iの非整数 の占有数n_iを導入することで全エネルギーの占有数に対する微分を計算する。非 整数にしたので、電子密度関数と運動エネルギーを再定義し、

n(r) = ∑

i

ni|ψi(r)|² (2.43)

T_s[n] = ∑

i

n_i⟨ψi|∇²|ψi⟩=

∑N i=1

n_iϵi−

∫

V(r)n(r)dr (2.44) とする。Ts[n]の微分は、

∂Ts[n]

∂n_i = ϵi−V|ψi|²= ϵi− ∂

∂n_i (

E_H+E_xc+

∫

V_ext(r)n(r)dr )

(2.45) となる。全エネルギーは

E_KS = T_s[n]+

∫

V_ext(r)n(r)dr+E_H[n]+E_xc[n] (2.46)

2.2. 第一原理バンド計算 19

で与えられるので、これをn_iで微分をとると、

∂EKS

∂n_i = ϵi (2.47)

を得る。これはJanakの定理[33]と呼ばれるもので、Hartree-Fock近似のKoopmans の定理[34]に相当する。Koopmansの定理はN個の多体系の電子からi軌道の準 粒子を1つ抜き出した時のエネルギーの変化がHartree-Fock方程式の固有エネル ギーϵiに等しいというものである。Koopmansの定理は差分となるが、この形式は 微分になっている。Koopmansの定理は、Hartree-Fockの固有値を準粒子のエネル ギーとして意味付けられるため重要な要素であったので、このJanakの定理も同様 に固有値を準粒子のエネルギーと同様に扱えるのではないかと想像される。その ことを実際に確かめる。ここで絶対零度における化学ポテンシャルの定義を思い だすと、

µ= dE

dN (2.48)

となる。金属を考えると、部分的に満たされているバンドは最大エネルギーの部 分のみなので、式(2.47),(2.48)から

µ= ϵF (2.49)

が成り立つ。フェルミエネルギーは真空準位との差であることから、Kohn-Sham 方程式の最大固有値は実際に観測される値になる。フェルミエネルギーの等高線 を波数空間で形成したものがフェルミ面だが、フェルミ面については何も与えて はくれない。しかし、実際のフェルミ面とよく似た形になることが予想される。そ れは、フェルミ面内の体積は実際のものと同等であることと、フェルミ面の形状 は主に静電的なポテンシャルが影響していることにある。Kohn-Sham方程式から 得られるバンド分散は実際に角度分解光電子分光などで観測される結果に近いも のになっているものと期待される。i番目にある軌道を取り出すときのエネルギー を考えると、テーラー展開することで

∆i = EKS(ni =0)−EKS(ni =1) (2.50)

= −∂E

∂n_i + 1 2

∂²E

∂n²_i +· · · (2.51)

= −ϵi+ 1 2

∂²E

∂n²_i +· · · (2.52)

ドキュメント内 電気通信大学大学院電気通信学研究科 (ページ 33-49)