電子数の制御 －エネルギーギャップへの応用

CMP Technical Report No. 14
「電子数の制御」 −エネルギーギャップへの応用− 白井光雲
大阪大学・産業科学研究所
２００３年２月６日
Department of Computational Nanomaterials Design
ISIR, Osaka University
目次
1
はじめに
1
2
E(N + 1/N − 1) の計算
2.1 電子数 Nel 可変 . . .
2.2 Si のギャップ計算 . .
2.3 Etot (N ) からの評価 .
2.4 KS 固有値からの評価
2
2
2
4
6
3
4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
8
非整数の場合への拡張
10
まとめ
0
1
はじめに
これまでの pwm は電子数 Nel は与えられた原子に対して中性になるよう価電子
数が決められていた。今度はそれをわざと変えてみる。これはドープする問題や、
ギャップの問題では常に必要とされるテクニックである。
例えば有名な密度汎関数理論（DFT）における絶縁体のギャップ Eg の問題と取り
上げてみる。通常非金属に対しては、エネルギーギャップはバンド図を書きその価
電子帯上端と伝導帯下端との差で求められる（KS 固有値ギャップと呼ばれる）。Si
の例で示すと、バンドは図 1 のようになり、これよりギャップは Γ から X 点近傍ま
での垂直遷移分 0.68 eV と見積もられる。もちろんこれは実験値の 1.1 eV とはかけ
離れている。
KS 方程式における固有値 ² にはもともと軌道エネルギーといったような物理的意
味が付加されていない。そこでエネルギーギャップのように明確な観測量を計算す
るには全エネルギーに立ち上って考える。
そうすると、エネルギーギャップ Eg は、電子数 N の関数としての基底状態全エ
ネルギー E(N ) を使って
Eg = I − A
= E(N + 1) + E(N − 1) − 2E(N )
(1)
と表される。ここに
A = −(E(N + 1) − E(N ))
(2)
I = −(E(N ) − E(N − 1))
(3)
I はイオン化エネルギー、A は電子親和力である。
一方で、DFT における固有値 ² は Janak により
Z 1
E(N ) − E(N − 1) =
²N (N − 1 + n)dn
(4)
0
という意味付けがなされている。ここに ²N (M ) は M 個の電子系の N 番目の固有
値である。さらに最高占有 KS 固有状態に関しては、有限温度 DFT の枠で、その
T → 0 極限で
²N (N − 1 + n) = ²N (N )
(0 < n ≤ 1)
(5)
とできるので、
²N (N ) = E(N ) − E(N − 1)
(6)
が成り立つ。式（6）が成り立つのは厳密な密度汎関数の解に対してであり、した
がって通常は近似的にしか成り立たない。
1
さらに有用な式として Slater による「遷移状態」の概念がある。それによると近似
1
E(N + 1) − E(N ) ≈ ²LUMO (N + )
2
(7)
が成り立つ。
成立にはいろいろ条件が要求されるが、ともかく近似的には最高占有状態に関し
てはイオン化エネルギーという物理的意味が付加される。したがって I 、A（の負の
符号をとったもの）はそれぞれ ²V 、²C と対応付けられることになる。これらの関係
式よりエネルギーギャップが求まる。このことを試して見る。
用いたプログラムバージョンは
VersionNo
= ’0.62’
以降のものである。
2
2.1
E(N + 1/N − 1) の計算
電子数 Nel 可変
電子数 Nel を中性の時から変えるにはどうしたら良いか？例えば一つ増やすには
どうするのか？それは単に中性の状態の上に更に上のバンドの占有数を１とするこ
とで計算できる。
ただしそのままでは電荷の中性条件が破られてしまう。運動量表現での全エネル
ギーの表式ではそのうち G = 0 の成分が電気的中性条件に関係している。すなわち
電気的中性条件が満たされ、G = 0 の成分は有限値の値を得ている。電子を一個加
えたとき、中性条件を保つため、一様分布した正の電荷を付け加える。これにより
電気的中性条件が保たれるので、結果的には G = 0 の成分には何の変更もない。単
に N 個から N + 1 個への電気的中性条件と再解釈されるだけである。
しかしながら、このような取り扱いには微妙な問題が潜んでいる。著者自身はま
だこの問題を良く理解していない。興味ある読者は文献 [3] などを参照されたい。
2.2
Si のギャップ計算
可変電子数 Nel の計算例として式（1）を用いて Si のギャップを計算してみる。そ
こで問題となるのは実際の計算上での E(N + 1) の評価方法である。本来は固体に対
2
Si PseudoPotential
Energy \Ry\
0.7
4 13
444
4
4
3
4
3
61
1
4
4
0.6
3
2
3 45
2
16
13
1
1
2
1
1
1
1
2
1
4111
1
1
1
1
1
1
4
1
44 1 1
44 1 1
444 1 1
4442
4
4
4
1
1
1
0.5
0.4
722
1 2
4
2
1
2
2
41
0.3
4
2
2
1
2
1
22
222
222223
11 11
1 11
Σ
KS X Z W Q L
4
0.2
Γ
533
11
2
1
1 2
1
2
1
1 2
1
33
3
3
37
3 5
4
3 1
1
EF
5
4
5
4
5
5
1
4
1
Λ
4
1112
Γ
∆
5
5
55
3
X
図 1: Si のバンド図。ギャップ付近を拡大している。
しては N は非常に大きな数で、それに対してたった一個の電子を増やしたり減らし
たりしたときの全電子エネルギー変化を調べなくてはならない。ということは Si 結
晶に対してであれば、基本単位格子で電子数を変化させたのでは全く不十分で、本
来は巨視的サイズでの格子で電子数を変化させるべきである。勿論これは望むべく
して現実には出来ないことである。したがって単位格子のサイズを変えながら、つ
まり元の Nel を変えながら式（1）で評価されるエネルギー差がどのように収束する
かを調べることになる。
ここでは、計算する基本格子として Si2 、Si8 、Si16 、Si64 を用いる。用いた結晶の
データを表 1 にまとめる。違ったサイズの結晶格子で結晶の対称性が違う。そのた
め k 点のサンプリングが違ってくるが、本来は k 点の取り方は違ったスーパーセル
でも等価の取り方をすべきであるが、今回の計算上では特にそれを意識していない
が、だいたいはコンシステントになっているはずである。
違ったサイズの結晶格子のエネルギーを比較するとき、平面波カットオフエネル
ギー Ecut を揃える必要がある。そうすることで原子あたりの平面波数は同じになる。
それらの条件を表 2 に示す。Osaka2k では制御入力パラメータは直接的には Ecut で
はなく、k 空間中のカットオフ半径 AMax なので、それにより Ecut の値には多少の
ばらつきはあるがほとんど問題ではない。
3
表 1: 計算に用いた Si のスーパーセル。格子定数 a0 の単位はÅ。
cell Nel
SG
a0
notes
7
Si2
8 Oh F d3m 5.4307 primitive
Si8
32 P 1
C1
5.4307
Si16
64 Td2 F 43m 10.8614
Si64 256 P 1
C1
10.8614
表 2: Si のスーパーセル計算の計算条件 。既約ゾーン中の k 点サンプリング数 Nkp 、
k 空間中のカットオフ半径 AMax、対応する平面波カットオフエネルギー Ecut （Ry
単位）、平面波数 Npw
cell
Nkp
AMax
Ecut
Npw
Npw /Nat
notes
Si2
2
3.10 10.8067
169
84.50
Si8
4
5.37 10.8093
654
81.75
Si16
2
6.20 10.8067 1314
82.13
Si64
1
10.74 10.8093 5185
81.02
2.3
Etot (N ) からの評価
実際に電子数 Nel を変化させるには、入力ファイル*.para で、
OPTION BEGIN
nel_add=
1
OPTION END
とする。こうすることで単位格子中の電子数は中性のときのものより nel add だけ
増加（減少）する。
収束の安定性 電子数を１個だけ増加（減少）させたとき、全電子数は奇数個となるので収束の
安定性が多少気になる。
k 点が複数の場合で一番電子数が多い Si16 の場合を見てみる。中性の場合は
iter
====
1
2
3
Eel
(Ry/cell)
===============
9.5020763960
7.8643393756
7.8469341575
deE
(Ry/cell)
============
-4.7996E+01
-1.6377E+00
-1.7405E-02
Xsi
(Ry^2/cell)
============
3.6411E-02
2.6137E-04
1.6311E-05
4
nst/bk
aglmax
========
5/ 0
5/ 0
5/ 17
=============
0.620789166
0.235631620
0.025393999
4
5
6
7
8
7.8463003601
7.8462401006
7.8462253908
7.8462244210
7.8462249223
-6.3380E-04
-6.0260E-05
-1.4710E-05
-9.6981E-07
5.0126E-07
1.1798E-06
2.5122E-07
5.8234E-08
7.1566E-09
2.1610E-09
5/ 72
5/ 90
5/101
5/109
5/101
0.003731780
0.001997683
0.000666872
0.000321314
0.000162000
nst/bk
aglmax
========
5/ 0
5/ 1
5/ 25
5/ 61
5/ 78
5/109
5/118
5/115
5/ 92
=============
0.604762496
0.268947995
0.054527724
0.010716083
0.003834572
0.001290196
0.000363940
0.000289749
0.000092865
となるが、電子数を一つ増やした場合
iter
====
1
2
3
4
5
6
7
8
9
Eel
(Ry/cell)
===============
10.1262996980
8.4283979662
8.4050073862
8.4038122039
8.4036558186
8.4036320772
8.4036315020
8.4036273269
8.4036301836
deE
(Ry/cell)
============
-4.7465E+01
-1.6979E+00
-2.3391E-02
-1.1952E-03
-1.5639E-04
-2.3741E-05
-5.7514E-07
-4.1751E-06
2.8567E-06
Xsi
(Ry^2/cell)
============
3.6597E-02
4.1819E-04
2.1022E-05
3.9332E-06
3.6424E-07
8.5243E-08
1.6488E-08
2.2911E-09
8.5004E-10
とやはり共役勾配過程で失敗する回数（bk）が突然多くなるが、しかし全体として
エネルギーの収束度からいえば ∼ 10−6 Ry と実際の精度上ほとんど問題とならない。
Si64 の場合は電子数が多い分だけ収束は遅くなるが、k 点が一つだけなので安定
性はそれほど問題ではない。実際、この場合は中性のときの
iter
====
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Eel
(Ry/cell)
===============
39.4083409177
32.1600726705
31.9133824319
31.8767385030
31.8738947204
31.8733096013
31.8731268874
31.8730542716
31.8730191932
31.8730051178
31.8729997078
31.8729972185
31.8729960103
31.8729955442
31.8729953585
deE
(Ry/cell)
============
-1.8952E+02
-7.2483E+00
-2.4669E-01
-3.6644E-02
-2.8438E-03
-5.8512E-04
-1.8271E-04
-7.2616E-05
-3.5078E-05
-1.4075E-05
-5.4100E-06
-2.4892E-06
-1.2082E-06
-4.6614E-07
-1.8567E-07
Xsi
(Ry^2/cell)
============
3.6985E-02
6.6840E-04
1.3689E-04
1.5963E-05
3.6173E-06
1.1903E-06
3.4591E-07
1.9099E-07
9.8166E-08
3.2006E-08
1.1721E-08
7.3993E-09
3.2588E-09
9.7859E-10
4.7128E-10
nst/bk
aglmax
========
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
5/ 0
=============
0.620190720
0.407579895
0.582486021
0.142005929
0.010937489
0.004255341
0.002290733
0.001549264
0.001007941
0.000558122
0.000345341
0.000278348
0.000231856
0.000118026
0.000073081
nst/bk
aglmax
========
5/ 0
5/ 0
5/ 1
5/ 1
5/ 0
5/ 1
5/ 1
5/ 1
=============
0.617882140
0.382059141
0.287082363
0.037494410
0.015890865
0.007659149
0.006868536
0.006724500
に比べ、電子数を一だけ増加したときの結果
iter
====
1
2
3
4
5
6
7
8
Eel
(Ry/cell)
===============
39.7764154606
32.6375926661
32.4133675175
32.4003089035
32.3986107976
32.3980759982
32.3978752400
32.3977927478
deE
(Ry/cell)
============
-1.8939E+02
-7.1388E+00
-2.2423E-01
-1.3059E-02
-1.6981E-03
-5.3480E-04
-2.0076E-04
-8.2492E-05
Xsi
(Ry^2/cell)
============
3.6735E-02
7.0302E-04
6.5882E-05
8.5378E-06
3.2766E-06
1.3508E-06
4.2587E-07
1.9017E-07
5
9
10
11
12
13
14
15
32.3977506175
32.3977259772
32.3977130894
32.3977056343
32.3977002230
32.3976963765
32.3976935863
-4.2130E-05
-2.4640E-05
-1.2888E-05
-7.4551E-06
-5.4114E-06
-3.8464E-06
-2.7903E-06
1.6198E-07
8.0368E-08
3.1138E-08
2.4975E-08
1.7627E-08
1.4768E-08
7.5920E-09
5/
5/
5/
5/
5/
5/
5/
0
0
1
0
1
0
1
0.009852706
0.004512182
0.005479048
0.005225518
0.004465208
0.003853468
0.003258505
を比べると収束の安定性は多少劣っているが、特に悪くなっていることは無い。
結果 電子数を増加（減少）させたときの全エネルギーの変化からエネルギーギャップ
の結果を見てみる。
表 3: 全エネルギーの電子数依存とそれより求まったエネルギーギャップ。エネル
ギーは Ry 単位
cell
N
E(N )
E(N + 1)
E(N − 1)
Eg
Si2
8
-15.805116
-15.230271
-16.119406
0.2606
Si8
32
-63.220490
-62.640899
-63.623797
0.1764
Si16
64
-126.528012
-125.970607
-126.932190
0.1532
Si64
256
-505.623950
-505.099252
-506.088238
0.0604
全電子数 Nel が小さいとき、当然電子数変化が全エネルギーに与える影響は相対
的に大きくなる。理想的には、全電子数 Nel が大きくなり全エネルギーの電子数依
存性が無くなった極限で、なおかつ残るエネルギー差、式（1）がエネルギーギャッ
プを与える。
図 2 ではこの結果を全電子数の逆数 1/Nel の関数としてプロットしている。1/Nel →
0 の漸近値として 0.6 eV くらいであろうか。この値はよくいわれる LDA ギャップが
実験値 1.1 eV の約半分ということに一致している。 2.4
KS 固有値からの評価
一般に KS 固有値には物理的意味はないが、最高占有準位のそれはイオン化エネ
ルギーという意味が与えれれている。
そこで最高占有 KS 固有値からもギャップを評価して（表 4）、それを全エネルギー
の差より求めたもの（表 3）と比較してみる。式 2、3 で求まる、電子親和力および
イオン化ポテンシャルと、KS 固有値でのそれとはおおまかには一致していることが
わかる。
6
[eV]
4
E gap
3
2
1
0
0.00
0.05
0.10
0.15
1/Nel
図 2: 1/Nel に対するエネルギーギャップ
表 4: 全エネルギーの電子数依存から求まるイオン化エネルギーと最高占有 KS エネ
ルギー固有値の比較。最後の列は KS 固有値ギャップ。エネルギーは Ry 単位
cell
N
A
²c
I
²v
²c − ²v
Si2
8
0.5748
0.5861
0.3142
0.3983
0.1878
Si8
32
0.5795
0.5751
0.4033
0.3983
0.1768
Si16
64
0.5574
0.5628
0.4041
0.4356
0.1272
Si64
256
0.5246
0.5240
0.4642
0.4633
0.0607
ところで、ここでの KS 固有値はごく限られた特殊 k 点だけで評価している。Si
ではそのバンドギャップは Γ 点から X 点近傍の点への遷移に対応している。それゆ
え特殊 k 点サンプリングによる KS 固有値ギャップで求まる値はサンプリングが荒く
なればなるほどバンド図で求まるものより大きくなる。ギャップをバンド図でなく
全エネルギー値の差から求めようとするならば、k 点を細かくする必要がある。バ
ンド図を書かせれば、Si2 の場合でも KS 固有値ギャップとしてだいたい 0.07 Ry く
らいの値を得ている。
7
3
非整数の場合への拡張
これまでは変化させる電子数は整数に限られていたが、ここでその制限を取っ払
い任意実数を扱えるようにする。電子数 Nel を非整数だけ変化させるには、入力ファ
イル*.para で、
OPTION BEGIN
rnel_add=
-0.75
OPTION END
のようにする。こうすることで単位格子中の電子数は中性のときのものより rnel add
だけ増加（減少）する。そうなっているか、出力ファイル pwm *.etot などで
Electron Parameters
Nel0prim
=
255
real =
255.250000
**** -0.750000 electrons have been added
NEDIM
=
128
nband
=
128
などのようになっていることを確かめるべきである。
Si2 まずギャップの精度を出すためにはほとんど役立たないが、原理上の興味で Si2 を
試みる。今度は k サンプリングを 10 点と増やしている。
Si2 では Nel = 8 である。それに対して非整数個 n の電子を付け足す（−1 ≤ n ≤ 1）。
それぞれの n に対して KS 固有値 ²N (n) が求まる。また全エネルギー差からも
²N ≈
E(N + n) − E(N )
n
(8)
でその準位を評価してみる。なおこの式は序論で出てきた式に比べて、証明が与え
られていないので根拠に乏しいが、直感的にはもっともらしいので使ってみる。図
3 に非整数 n の場合の、KS 準位による価電子帯上端（VBM）、伝導帯下端のエネル
ギー（CBM）、およびそれを式（8）で評価したものを比較している。
なお、今回は k サンプリングを 10 点と増やしているので、式（1）から見積もれ
るギャップの値は前節のもの（0.26 Ry）より明らかに改善されている（0.21 Ry）。
まず図 3 で気づくことはデータがばらついていることである。滑らかな変化とは
なっていない。収束に問題があるのだろうか？
8
0.7
CB
VB
Etot
KS
Energy
0.6
0.5
0.4
0.3
0.0
1.0
n
図 3: n に対する固有値および全エネルギー差から求まった価電子帯上端（VBM）、
伝導帯下端のエネルギー（CBM）の準位
-505.0
0.55
ev
ec
-505.2
0.53
Ev
Ec
KS level
Etot
-505.4
-505.6
0.51
0.49
-505.8
0.47
-506.0
-506.2
-1.0
-0.5
0.0
0.5
0.45
1.0
-1
0
1
n
n
図 4: n に対する全エネルギー変化および KS 固有値による価電子帯上端（VBM）、
伝導帯下端のエネルギー（CBM）の準位。ev、ec は KS 固有値であり、Ev、Ec は
全エネルギー差から求めたもの。
9
Si64 つぎに Si64 の場合で調べてみる。この場合は k サンプリングは Γ の 1 点だけで
ある。
図 4 に n に対する全エネルギー変化および KS 固有値の結果が示される。これか
ら全エネルギーが n の滑らかな関数となっている。
KS 準位と全エネルギー差から求めた価電子帯上端（VBM）、伝導帯下端の（CBM）
のエネルギーは良く一致していることがわかる。すなわち近似として式（6）が良く
成り立っていることがわかる。
4
まとめ
電子数を可変とすることはうまくいっている。収束速度も悪くない。
その応用例としてのギャップの求め方を示した。ある程度大きなスーパーセルを
用いて、電子数１コ分の変化の影響が少なくなれば全エネルギー差から求まるイオ
ン化エネルギー、電子親和力と KS 固有値でのそれはほぼ一致している。
参考文献
[1] J. F. Janak, Phys. Rev. B18 7165 (1978).
[2] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Jr., Phys. Rev. Lett. 49
1691 (1982).
[3] G. Makov and M. C. Payne, Phys. Rev. B51 4014 (1995).
10