User Manual III

第一原理擬ポテンシャル法による
電子状態計算プログラム
「Osaka2002 nano」
III 「擬ポテンシャルの拡張編」
第 1.6 版
大阪大学・産業科学研究所
白井光雲
22 May 2006
History
Rev. 1.8 28 Jun 2013
1
目次
はじめに
第1章
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
擬ポテンシャルの理論
ノルム保存擬ポテンシャル
Kleinman-Bylander 近似 .
局所スピン密度近似 . . .
相対論的な一般式 . . . . .
スカラー相対論 . . . . . .
スピン軌道相互作用 . . .
内殻補正 . . . . . . . . . .
将来展望 . . . . . . . . . .
4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
第 2 章 原子擬ポテンシャル、atom の使い方：スピンなし
2.1 擬ポテンシャル生成 . . . . . . . . . . . . . . . .
2.1.1 出力ポテンシャルデータファイル . . . . .
2.2 参照原子の価電子と電子配置の取り方 . . . . . . .
2.2.1 いくつかの注意すべき例 . . . . . . . . . .
2.3 局所成分の取り方 . . . . . . . . . . . . . . . . . .
第3章
3.1
3.2
3.3
3.4
原子擬ポテンシャル、atom の使い方：拡張
スピンの擬ポテンシャル生成 . . . . . . . . .
内殻補正 . . . . . . . . . . . . . . . . . . . .
相対論的効果 . . . . . . . . . . . . . . . . .
GGA 補正 . . . . . . . . . . . . . . . . . . .
3.4.1 移植の原理・方針 . . . . . . . . . . .
3.4.2 ポテンシャルの作成 atm . . . . . .
3.4.3 GGA の使い方 . . . . . . . . . . . .
第 4 章 ポテンシャルデータの再構築
4.1 ポテンシャルディレクトリ . . . . .
4.2 GGA ポテンシャルの作成 atm . .
4.2.1 GGA ポテンシャルの作り方
4.2.2 inip . . . . . . . . . . . . .
4.2.3 pwm . . . . . . . . . . . . .
4.3 非線形内殻補正 . . . . . . . . . . .
2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
7
10
12
14
14
14
16
.
.
.
.
.
18
18
21
22
25
25
.
.
.
.
.
.
.
29
29
31
31
33
33
34
35
.
.
.
.
.
.
37
37
38
39
40
40
40
第5章
5.1
5.2
5.3
5.4
第6章
6.1
6.2
6.3
Osaka2k でのスピンの取り入れ方
スピンを入れた計算の流れ . . . . .
電子占有数のオプション . . . . . .
実例 . . . . . . . . . . . . . . . . .
5.3.1 O2 結晶 . . . . . . . . . . .
5.3.2 Si 結晶中の Fe 不純物 . . . .
関連した技法 . . . . . . . . . . . .
5.4.1 局在モーメントの半径 . . .
5.4.2 原子位置の最適化 . . . . . .
5.4.3 部分 DOS . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
42
42
43
44
44
46
48
48
48
48
内殻補正の計算
49
atom で内殻補正のための入力 . . . . . . . . . . . . . . . . . . . . . . . . . 49
pwm での実行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
非線形内殻補正のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
第 7 章 相対論の計算
7.1 相対論的な一般式 . . . . . . . . . . .
7.2 移植上いくつかの注意点 . . . . . . .
7.2.1 スピノールの表現 . . . . . . .
7.2.2 Kleinman-Bylander 型の使用
7.2.3 電子の交換相関エネルギー . .
7.2.4 波動関数の初期値 . . . . . . .
7.3 SCF 計算 . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
52
52
53
53
54
54
54
55
付 録 A 擬ポテンシャルテスト
58
付 録 B その他
60
参考文献
61
3
はじめに
擬ポテンシャル法も、スピン、内殻補正をはじめとしてポテンシャルに関して拡張され
てきている。その上、Osaka2k にも GGA も取り入れたりしたため、ポテンシャルに種類
が途端に膨れ上がってきた。そこでいくつもの擬ポテンシャルの型ができ、それらの間に
混乱がないように首尾一貫した取り扱うが求められる。
4
第1章
擬ポテンシャルの理論
擬ポテンシャルについては上巻でも説明されているが、ここではもう少し掘り下げて
みる。
1.1
ノルム保存擬ポテンシャル
固体の多くの物理的、化学的性質は広がった価電子状態により決まるという事実に基づ
いて、「元の急激に変化するポテンシャルに比べて緩やかに変化するポテンシャルを用い
て、価電子状態を再現することができないだろうか？」という問い掛けがなされた。これ
が擬ポテンシャル法という考えが生まれた原点であり、
「如何にうまく価電子状態をシミュ
レーとするか」が擬ポテンシャル法の永遠のテーマである。具体的にそのような擬ポテン
シャルに求められる条件は、まず
1. 価電子状態のエネルギーレベルが再現されること。
2. 平面波での計算が容易であるように波動関数が滑らかであること。つまり節なしの
波動関数になること。
となる。このような要請に従いいくつかの擬ポテンシャルモデルが提案されたが [1, 2]、
後になりそれが一般に次のように内殻状態への射影演算子として表されることがわかった
[3]。
∑
V̂ps =
Fc (r)|cihc|
(1.1)
c
ここに |ci は内殻状態である。Fc (r) は任意関数であり、それゆえ擬ポテンシャルの作り方
には事実上無限にあることになる。式（1.1）で与えられる擬ポテンシャルは、定義の示す
通り非局所関数であるばかりか、一般にはエネルギー依存性を持つ。それらに付随する計
算の煩雑さは様々な近似で解消されてきた。それらはいわゆる「モデル擬ポテンシャル」
と呼ばれるものである。それらに関するレビューは文献 [4] に詳しい。
これまでのところは、
「価電子状態」といっても、例えば同じ C 原子であっても、それが
ダイヤモンド構造を取るときと、グラファイト構造を取るときでは異なっており、それぞ
れについて擬ポテンシャルを作らねばならなかった。このことは大きな欠点で、擬ポテン
シャルを使う立場からすると、できれば、「例えば C であればプロトタイプの原子状態の
C で擬ポテンシャルを作ったならば、それがダイヤモンド結晶でもグラファイト結晶でも
通用するものであって欲しい」と願うのは当然である。このような性質を持つポテンシャ
ルは「transferable」であるといわれる。このような観点から、エネルギーの一致だけでな
く、価電子密度の一致にも注意が払われるようになった [5]。
5
ポテンシャルを transferable なものとし、擬ポテンシャルを「非経験的第一原理計算」た
らしめたのは、Hamann、Schlüter、Chiang による「ノルム保存」という考えである [6]。
それは
3. 価電子状態の擬波動関数は、あるコア半径 rc より外では真の波動関数と一致する。
4. 半径 rc 内での擬電荷密度と真の電荷密度の積分は同じであること。すなわち波動関
数のノルムは保存されること。
という条件である。参照原子に対して角運動量子数 l についての動径方向のシュレーディ
ンガー方程式は、
d2 ul
− 2 + Vl (r)ul (r) = εl ul (r)
(1.2)
dr
である。ここに Vl は遠心力項 l(l + 1)/r2 を含む。原子に対して各価電子軌道 wl からこの
ような擬波動関数を作り、それから、逆にシュレーディンガー方程式をポテンシャルにつ
ps
ps
いて解いたものが擬ポテンシャル Vl となる。実際にはこの得られたままの Vl は価電子
の電子間相互作用の成分が含まれていて、それは問題とする原子置かれている環境で変化
するものである。それでその寄与を差し引いた、
Vlps,ion (r) = Vlps (r) − VH (ρv (r)) − Vxc (ρv (r))
(1.3)
が裸のイオンポテンシャルとして様々な応用に用いられることになる。ここに VH はハー
トリー項で、Vxc は交換相関ポテンシャルである。こうして全ての原子に対してシステマ
テックに擬ポテンシャルが生成され [7]、以後、擬ポテンシャルは、ポテンシャルの作成
と、それを用いた KS 方程式を解くことに分離された。以降特に断りの無いかぎり、Vl あ
ps
ps,ion
るいは Vl はこの裸の Vl
を意味するものとする。
このような条件を満たすような擬波動関数の取り方は漸近的性質さえ満たしていれば
任意であるから、このようなノルム保存の擬ポテンシャルの作り方は無限にある。そこで
こうしてノルム保存以外に、波動関数を滑らかにする様々な付加条件を課して収束性を向
上させる「拡張ノルム保存」型が様々に考案された [8]。「OSAKA2000」で使われている
Troullier-Martins 型もこの延長上のもので、rc における波動関数の滑らかさ、解析的性質
の他、r = 0 における曲率が 0 の条件などを課して作られたものである [9]。
式（1.3）で得られた擬ポテンシャルは、軌道 l ごとに違うポテンシャルであるので、軌
道 l への射影演算子の形で表される。
V̂ ps =
∑
|liVl hl|
(1.4)
l
この形の擬ポテンシャルはエネルギー依存性がないので、この点で扱い易い。そのかわり
軌道 l 依存性を持つことになるのであるが、その結果式（1.4）が示す通り、非局所ポテン
6
シャルの形となる。1 非局所性を明示的に書くと
V ps (r, r0 ) = hr|V̂ ps |r0 i
∑
=
Ylm (r̂)[Vl (r)δ(r − r0 )]Ylm (r̂0 )
(1.5)
lm
となる。非局所的性格を持つ点では少し厄介であるが、球面調和関数に分解されるため、
平面波展開では部分波の方法が使え数値計算上は比較的取り扱いやすい。平面波展開の具
体的表式は
∫ ∞
4π(2l + 1)
Vlnl (q, q0 ) =
Vlnl (r)jl (qr)jl (q 0 r)r2 dr
(1.6)
Pl (cos γ)
Ω
0
である。ここに γ は q と q0 とのなす角である。
数値計算上は、式（1.4）は、球面調和関数の完備性を利用し、
∑
V̂ ps = V̂ loc +
|li∆Vlnl hl|
(1.7)
l
のように、局所成分 V loc (r) と、非局所成分 ∆Vlnl (r) に分けられる。ここに非局所成分
∆Vlnl (r) は
∆Vlnl (r) = Vlps (r) − V loc (r)
(1.8)
で与えられる。局所成分 V loc (r) の取り方は任意であるが、通常は原点から遠くの場所で
漸近的に ∼ Zv e2 /r となるように選ぶ。こうすることで、非局所成分 ∆Vlnl (r) は有限の半
径内だけで零でない有限値をとるので、非局所成分に関する実空間積分の計算が楽になる。
擬ポテンシャルに課せられる条件（4）は、半径 rc より外の点で、擬波動関数と真の波
動関数のエネルギーに関する対数微分が等しくなることを保証する。これにより擬ポテン
シャルはある参照エネルギーにおいて電子の散乱特性を正しく記述するばかりでなく、そ
の参照エネルギーからずれたところでもエネルギーのずれの 1 次のオーダーまでは正し
く記述されることが保証される。こうして異なった化学環境の中でも、この裸の擬ポテン
シャルを使い信頼できる結果を得ることができるのである。2
1.2
Kleinman-Bylander 近似
式（1.5）で定義された擬ポテンシャルは半局所的（V̂ SL ）と呼ばれることがある。位置
の角度成分に関しては非局所的であるが、動径方向成分に関しては局所的であるからであ
る。このため、その平面波展開の行列要素は式（1.6）のようにきちんと解析解で与えら
れるものの、展開平面波の数を Npw とすると、計算しなければならない行列要素の数は
Npw (Npw +1)/2 で、
（1.6）の積分回数は Npw の自乗に比例して増加する。また、Npw ×Npw
1
ポテンシャルが非局所的である場合、密度汎関数理論の元の形 Hohenberg-Kohn の定理に従うと、エネ
ルギー汎関数は密度 ρ(r) の汎関数でなく一密度行列 ρ(r, r0 ) の汎関数で与えられるものとなる [10, 11]。
2
擬波動関数と真の波動関数のエネルギーに関する対数微分が等しい条件は transferability の良いテスト
であることには確かだが、それを絶対的条件として過大評価することは危険である。[9] 参照。
7
のサイズの非局所ポテンシャル行列を波動関数に作用させるには Npw 2 回の浮動小数点演
算が必要となる。それゆえ大きな Npw の計算ではこの非局所ポテンシャルの計算がボト
ルネックとなってくる。
それに対して、Kleinman、Bylander は動径方向成分も非局所的に扱う完全非局所型を
提唱した [13]。それによると、
¯
¯
∑ ¯u0 ∆Vl ih∆Vl u0 ¯
lm
lm
KB
loc
V̂
= V̂ +
(1.9)
hu0lm |∆Vl | u0lm i
lm
ここに u0lm は擬ポテンシャルを作るときの原子の擬波動関数で、∆Vl は式（1.8）と同じ
ようにして V loc を切り取ったものである。これを Kleinman-Bylander 型の完全分離型非
局所ポテンシャルと呼ぶ。
V̂ KB |u0lm i = V̂ SL |u0lm i
(1.10)
と擬ポテンシャルを作った原子の場合に対しては、この V̂ KB は実効的に元の擬ポテンシャ
ルと同じとなる。
この完全非局所型ポテンシャルは、座標が右と左で完全に分離しているので、行列要素
として必要なものは Npw 個の h∆Vl u0lm |qi のセットで良く、積分計算の負担が激減できる。
加えて、非局所ポテンシャルのハミルトニアン行列の波動関数への作用も Npw のオーダー
の浮動小数点演算で済む。それゆえ計算効率から言うと、この Kleinman-Bylander 型完全
分離ポテンシャルは可能な限り使いたいところである。この「OSAKA2000」でも基本的
には全ての原子に対して Kleinman-Bylander 型を使うことになっている。
この Kleinman-Bylander 型完全分離ポテンシャル（略して KB ポテンシャル）の主な問
題は、式（1.9）の等号は正にこのポテンシャルを作った原子にしか成り立たず、他の化学
環境では、近似的にしか成り立たないということである。
V̂ KB |ui ' V̂ SL |ui
(1.11)
こうして元の擬ポテンシャルに比べ transferability が犠牲となる。さらに精度が少し落ち
るくらいならまだしも、時々真の束縛状態の下に偽の状態が現れることがあることが次第
に明らかになった。
この問題に対して、Kleinman-Bylander は「賢い V loc の選び方をすることで transferability を保つことができる」としている。その V loc の選び方としてどうすればよいか提
案者の Kleinman-Bylander の言葉に耳を傾けよう。彼らは具体的な条件として次の、
1. ∆Vl はできる限り短距離力でかつ小さくすること。
2. もし ∆Vl の内どれかが小さくできなければ、正の ∆Vl の値はできるだけ小さくする
こと。
¯
¯
3. ¯hu0lm |∆Vl |u0lm i¯ /hu0lm | |∆Vl | |u0lm i はできる限り 1 に近づけること。
という三つ挙げている [14]。3
3
この条件は文献 [17] によれば「経験的」とされているが、筆者は十分本質的なことを言い当てていると
思う。
8
KB ポテンシャルが参照状態とは異なる化学的環境ではもはや近似に過ぎないことは定
義から明白であるが、
「偽の状態」
（ゴースト準位）がでることは非常に困ったことである。
「偽の状態」とはもう少し具体的に述べると、真の価電子状態は擬ポテンシャル近似では節
の無い波動関数のはずで、従ってそれ以下のエネルギーのところには状態はないはずだ。
ところが節を持つ状態がそれよりも低いエネルギーで現れるのである。常識的には同じ境
界条件の元では「固有値エネルギーは節の数とともに増大する」（Wronskian の定理）は
ずなので、これは何かがおかしい。
この病理学的振る舞いの原因はその後の研究で明らかになった [15, 16, 17]。まず Wronskian の定理は局所ポテンシャルに対して成立するのであり（実際擬ポテンシャルを作る
ときの式（1.2）ではポテンシャルは局所的である）、KB ポテンシャルのような非局所ポ
テンシャルに対しては一般には成立しない [16]。また KB ポテンシャルを式（1.9）のよう
に射影演算子の形で見ると、この演算子ポテンシャルがある波動関数に作用するとき、た
またまそれを消してしまう（プロジェクトアウトする）ことが起こりうる。このときはこ
の波動関数は非局所ポテンシャルを感じることなく、局所成分のみの影響を受ける [15]。
局所成分がもし負に大きいならば、これは深い準位を作ることになる。これがゴーストの
正体である。
このゴーストがでるかどうかの条件は Gonze らによって詳しく分析された [16, 17]。そ
れによると、式（1.2）における波動関数 |∆Vl ul i を規格化したものを |uKB
l i として、また、
ElKB =
hu0l ∆Vl |∆Vl u0l i
(∆Vlrms )2
=
hu0l |∆Vl | u0l i
hu0l |∆Vl | u0l i
(1.12)
KB は ∆V̂ KB に対する
で定義されるパラメータ ElKB を KB エネルギーとする。|uKB
l i、El
l
固有関数、及びその固有値となる。こうして式（1.9）は
¯
¯
KB KB ¯
∆V̂lKB = ¯uKB
l iEl hul
(1.13)
の形に簡約される。この形にすると、KB エネルギー ElKB は局所ポテンシャルに対する非
局所ポテンシャルの影響の大きさを量るパラメータであることがわかる。しかし KB エネ
ルギーの定義（1.12）によると、ElKB が大きくなるのは ∆Vl が大きいとき（本質的に大
きい）ばかりでなく、たまたま分母が小さくなったときにもそうなる（人為的に大きい）。
ElKB が大きくなったのは、本質的に ∆Vl が大きいことによるものか、あるいは人為的な
原因によるものかを判断するには、∆Vlrms と ElKB の比を取ればよい。
ClKB =
∆Vlrms
0
= huKB
l |ul i
ElKB
(1.14)
この式よりわかることは、KB エネルギー ElKB が人為的に大きくなっているときは ClKB
0
（KB コサインと呼ぶ）は非常に小さくなっていて、従って |uKB
l i は |ul i とほぼ直交して
いる格好になり、非局所ポテンシャルの影響はほとんど受けないことになる。ElKB が正の
ときは局所ポテンシャルが負に大きいときであるから、このときゴーストが現れることに
なる。こうしてゴーストが現れる条件は
1. ElKB > 0 のとき かつ ClKB の絶対値が非常に小さいこと
9
と要約できる。この条件は Kleinman-Bylander の挙げた条件の３に相当する。では ElKB < 0
の場合は問題はないのかというと、文献 [16, 17] によれば、この場合にもゴーストが出る
条件が書かれているが、ClKB などの兼ね合いでどういう条件となるのか筆者には良くわ
からない。
ではこのゴーストが現れたときどうすればそれを解消できるか。簡単な解決法は擬ポテ
ンシャルを作るときのカットオフ半径を変えることである [16]。ElKB はカットオフ半径に
強く依存するため、適切なカットオフ半径を作ることで大抵は解決できる。それと局所成
分の取り方である。経験の示すところによると、多くの半導体では問題が起きるのは KB
エネルギーが人為的に大きくなる場合のみでその時のみカットオフ半径を変えればよい。
一方、遷移金属の場合は KB エネルギーが大きいのは本質的にそうであり、その場合は
ElKB が負にならないように局所ポテンシャルの取り方を調整すべきである。
1.3
局所スピン密度近似
多電子の交換相関相互作用の密度汎関数はわかっていない。ほとんどの場合は、局所密
度近似を使う。それに関しては第一巻で簡単に触れたが、そこではスピンの効果は入って
いなかった。ここではスピンも含めた局所スピン密度汎関数について述べる。4
スピンが入った場合でも、第一巻で述べた密度汎関数理論は、その密度 ρ(r) のところ
を、スピン密度 ρ↑ (r)、ρ↓ (r) に置き換えれば基本的には成り立つ [18]。つまり、基底状態
のエネルギーは、スピン密度に関するエネルギー汎関数の最小値で与えられる。そしてそ
の汎関数を最小にするスピン密度が基底状態のスピン密度を与える。5
スピンの違う状態は異なるポテンシャルを感じることになるが、それが良くわかってい
ないスピン密度汎関数の局所近似に頼っているので、その妥当性に関する議論が昔から果
てしなく延々と続くのである。こうした議論はもちろん本書の枠外のことであり、ここで
はもっぱら局所スピン密度近似に頼ることになる。
非相対論的な場合のスピンの効果は、元のイオン核のポテンシャルはスピン依存がない
ので、もっぱら交換相関ポテンシャルを通じて生じる。
交換相関相互作用の局所密度近似を使った密度汎関数の具体的表式を求める必要がある。
相関相互作用の部分については解析解はないので数値解に頼らざるを得ない。それを求め
る際の基礎データとして一様密度電子ガスのエネルギーをモンテカルロ法などで計算した
ものを使い、密度の関数としてパラメータフィッティングすることで得られる。ここでは、
Ceperley と Alder のシュミレーション結果 [19] を使い、Perdew と Zunger がパラメータ
化したものを使う [20]。
4
ここで言うスピンの効果は後に扱う相対論的効果の場合とは違う。ここのスピンの効果は多電子系におけ
る波動関数の対称性に基づく実効的相互作用である。
5
しかしスピンの場合はスピン密度とそれを与えるポテンシャルの間には 1 対 1 対応はない [18]。その場
合でも 7 ページの脚注 1 で述べられた Levy の拘束付き探索によりやはりスピン密度 ρ↑ (r)、ρ↓ (r) の汎関数
として最小化する方法は有効である。
10
それによると、電子密度 ρ に対して（あるいはスピン密度 ρ↑ 、ρ↓ に対して）、
ρ = (4πrs3 /3)−1
(1.15)
ζ = (ρ↑ − ρ↓ )/ρ.
(1.16)
で定義される平均電子半径 rs およびスピン分極率 ζ を変数として表される。以下、上付
き文字 i は、i = P で分極なし（すなわち ζ = 0）、i = F で 100 %の分極（ζ = 1）の場合
を示す。
交換エネルギー ²x は単に
0.916
rs
0.916
²Fx = −21/3
rs
²Px = −
(1.17)
(1.18)
と表される。
相関エネルギー ²c の方は複雑で、低密度では、
√
²ic = γ i /(1 + β1i rs + β2i rs )
(1.19)
で与えられ、対応する相関ポテンシャルは
µic = ²ic
√
(1 + 76 β1i rs + 43 β2i rs )
√
(1 + β1i rs + β2i rs )
(1.20)
となる。
一方、高密度側では、
²ic = Ai ln rs + Bi + Ci rs ln rs + Di rs ,
(1.21)
で与えられ、対応する相関ポテンシャルは
1
2
µic = Ai ln rs + (Bi − Ai ) + Ci rs ln rs ,
3
3
(1.22)
となる。
スピン分極率が中間領域（0 < ζ < 1）では、Barth と Hedin による内挿方式が良く用
いられる [18]。それによると
²c (rs , ζ) = ²Pc (rs ) + f (ζ)[²Fc (rs ) − ²Fc (rs )],
µc (rs , ζ) =
µPc (rs )
+
f (ζ)[µFc (rs )
−
(1.23)
µPc (rs )]
+[²Fc (rs ) − ²Pc (rs )][sgn(σ) − ζ]
df
,
dζ
(1.24)
ここに f (ζ) は
f (ζ) =
(1 + ζ)4/3 + (1 − ζ)4/3 − 2
24/3 − 2
で与えられる。
11
(1.25)
Table 1.1: Ceperley と Alder の相関エネルギーのシミュレーション結果のパラメータ化し
た値
Ceperley-Alder
Parameter
P
F
γ
-0.1423
-0.0843
β1
1.0529
1.3981
β2
0.3334
0.2611
A
0.0311
0.01555
B
-0.048
-0.0269
C
0.0020
0.0007
D
-0.0116
-0.0048
参考のため、表（1.1）にパラメータ化したときの値を記録しておく。
こうして参照原子で、スピンアップ、ダウンに対してそれぞれ実効的ポテンシャルが求
まる。ここでスピンが無い場合の式（1.3）に対応して、裸のイオンポテンシャルは
σ
Vlσps,ion (r) = Vlσps (r) − VH (r) − Vxc
(r)
(1.26)
σ
で求まる。式（1.26）で与えられる裸のイオンポテンシャルは、スピン依存の部分が Vxc
で完全には消去できず、一般にはスピン依存を持つ。このスピン依存を持ったポテンシャ
ルをそのまま使うこともできるが [21]、これは応用上好ましくないように思える。なぜな
ら、非相対論近似では、スピン依存は全て電子同士の交換相関項から生じるものですので、
スピン依存をあからさまに持ったイオンポテンシャルは transferability の点で信頼がおけ
ないと考えられるからだ。そこでそのスピン依存を持ったイオンポテンシャルをスピン占
有数で平均化させた
Vlps,ion (r) =
ρ↑ (r) ps
ρ↓ (r) ps
Vl↑ (r) +
V (r)
ρ(r)
ρ(r) l↓
(1.27)
が一般には使われる。しかし遷移金属の磁性などのスピンの微妙な効果はこれでどこまで
再現できるかは疑わしい。それで、あからさまなスピン依存性は持たせず、そのかわりス
ピン分極の依存性を持った項を付加した擬ポテンシャルでそれらの問題を大幅に改善され
ることが示された [22]。
1.4
相対論的な一般式
ノルム保存型の擬ポテンシャル法は相対論的領域へも拡張できる [23, 24]。擬ポテンシャ
ルの理論の完全性を期すため、以下相対論への拡張について述べる。
ディラック方程式は 4 成分波動関数ベクトルに対して書き下される。その 4 成分波動関
数ベクトルは上半分と下半分に分けることができ、エネルギーが正の解は、上半分は大成
12
分 G、下半分は小成分 F と呼ばれる。球対称の場合、全ハミルトニアンと可換な量子数
は、全角運動量 J = l + s、その z 成分 M 、およびヘリシティを示す演算子 K
K = β(Σ · l + 1)
の固有値 κ である（β = v/c）。κ と j の関係は

−(l + 1) = −(j + 1 ),
2
κ=
l = (j + 1 ),
2
(1.28)
for j = l +
for j = l −
1
2
1
2
(1.29)
となる。このとき、2 成分ベクトルである大成分、小成分は、角座標およびスピン成分に
関しては、
)
)
(
(
l + m + 1 1/2 m
l − m 1/2 m+1
J
σ↓
ΦM =
Yl σ ↑ +
Yl
2l + 1
2l + 1
1
1
J = l + ,M = m +
(1.30a)
2
2
(
)1/2
)1/2
(
l−m+1
l+m
0
ΦJM 0 =
Ylm−1 σ↑ −
Ylm σ↓
2l + 1
2l + 1
1
1
(1.30b)
J0 = l − , M0 = m −
2
2
とスピノル関数で表され、動径方向に関してはスカラー関数 G(r)、F (r) となる。コア領
域外での価電子を扱っている限り、波動関数の小成分 F は大成分 G に比べて α（微細構
造定数 = 1/137）の程度の大きさで、従って密度に関しては小成分は α2 のオーダーで無
視できる。それゆえ波動関数は大成分 G に関してのみのシュレーディンガー方程式を解く
ことで十分である。
( 2
)
d G κ(κ + 1)
−
−
G + [V (r) − ²]G = 0
(1.31)
dr2
r2
この方程式を解くことで得られる波動関数より、第 1.1 節の手続きを踏むことで擬波動関
数、そして擬ポテンシャルが得られる。
∑ [
]
V̂ ps =
|li V̄l (r) + VlSO (r)L · S hl|
(1.32)
l
ここに、V̄l (r) はスピンアップとダウンについて平均したイオンのポテンシャル
V̄l (r) =
]
1 [
lVl−1/2 (r) + (l + 1)Vl+1/2 (r)
2l + 1
(1.33)
で、VlSO (r) はスピン軌道相互作用のポテンシャル
VlSO (r) =
]
2 [
Vl+1/2 (r) − Vl−1/2 (r)
2l + 1
(1.34)
である。
ここで重要なことは、全電子計算の際は、相対論的補正は (Zα)2 の形で現れる。従って
重い原子では上で述べた波動関数の小成分を無視することはもはやできない。それに対し、
コア領域外の価電子のみを対象とする擬ポテンシャルに場合は、相対論的補正は α2 での
み効いてくるということである。
13
1.5
スカラー相対論
ポテンシャルの式（1.32）で、第二項のスピン軌道相互作用が無視できるときは、スピ
ンアップとダウン成分は混ざらず、それぞれ別々に解くことができ、このときの方程式は
スカラー相対論的方程式と呼ばれる。l は良い量子数のように取り扱え、方程式は l ごと
に立てられる。
1.6
スピン軌道相互作用
ポテンシャルの式（1.32）で、第二項のスピン軌道相互作用が無視できないときは、前
節のスカラー方程式と対応してベクトル方程式と呼ばれれ、スピンアップとダウン成分は
混ざり（同時に 1 だけ異なる l どうしも混ざる）、連立方程式を解かなければならない。
解き方として、スピン軌道相互作用ポテンシャルを摂動論として扱うものと、初めから
それを含めた形で解くものがある。
1.7
内殻補正
原子状態での電荷分布を内殻、価電子軌道に分解してみると Si の例では図 1.1 のように
なっている。
一方で、内殻電子が広がっているといわれているアルカリ金属原子では図 1.2 のように
なっている。これらの比較から、Si では内殻、価電子それぞれが空間的にはっきり分離さ
れていて、Na の場合はそれらの電荷密度分布の重なり合いが大きくなっているといえる
だろうか？図で見ると内殻、価電子の電荷密度を示す赤、青線の重なるところはだいたい
同じあたりにあるように見える。それよりも見た目で Si と比較して Na の特徴は、そもそ
も内殻、価電子の電荷密度のピーク値の比が全然違う。
密度汎関数理論では電子間相互作用によるポテンシャル Vee (r) は電子密度 ρ(r) だけの
関数で与えられる。擬ポテンシャル法では、便宜上全電子密度 ρ(r) のうち価電子 ρv (r) に
ついて理論を書き直したものである。
ハートリー項は電荷密度に関して線形であるので、ハートリーポテンシャルに関しては
この書き直しは簡単となる。
VH (ρ(r)) = VH (ρc (r) + ρv (r)) = VH (ρc (r)) + VH (ρv (r))
(1.35)
擬ポテンシャル構築時に、遮へいされたポテンシャル Vscr (r) からこのハートリーポテ
ps
ンシャルを除いて裸のイオン擬ポテンシャル Vion としている。
ps
Vion
(r) = Vscr (r) − VH (ρv (r))
(1.36)
ps
したがってこの Vion には内殻の VH (ρc ) の効果が入っているものと考えることが出来る。
交換相関ポテンシャルに関しても、単純に式（1.36）の処理を施し、裸のイオンポテン
シャルとしている。しかしながら交換相関項に関しては密度とは線形の関係にないので、
本来は正しくない。つまり
VXC (ρc (r) + ρv (r)) 6= VXC (ρc (r)) + VXC (ρv (r))
14
(1.37)
Si
20
chds_c
rho
chds_v
10
0
0
1
2
3
4
5
r
図 1.1: Charge densities of the true core and valence electrons of silicon. The charge
density in this figure are actually those given by 4πr2 ρ(r).
Na
1.0
15
chds_c
chds_v
0.8
10
rho
rho
0.6
0.4
5
0.2
0.0
0
0
1
2
3
4
0
5
1
2
3
4
5
r
r
図 1.2: Charge densities of the true core and valence electrons of Na. The charge
density in this figure are actually those given by 4πr2 ρ(r). The right-side figure shows
an expansion of the left-side one.
15
なのである。それを線形と仮定したことにより誤差が生じうる。内殻と価電子の電荷密度
に重なりが小さいときはこの誤差は無視できるが、その重なりが大きいときに重要となる。
ここに非線形の効果が現れる。
この誤差を解消する方法は、擬ポテンシャル構築時に内殻電荷による交換相関項を含め
てその遮へいポテンシャルから引いておく。すなわち
ps
(r) = Vscr (r) − VH (ρv (r)) − VXC (ρc (r) + ρv (r))
Vion
(1.38)
こうして得られた擬ポテンシャルを使うときは、内殻電子の電荷密度を常に含めた交換
相関ポテンシャルを当てはめる。
部分補正
原理上はこれで非線形交換相関補正の効果を正しく取り入れることが出来る。しかしな
がら内殻電子の電荷密度分布は図 1.1 や 1.2 にみられるように非常に局在した形のもので
ある。それゆえ平面波展開の方法においては、この ρc (r) を忠実に再現する点で大きな桎
梏となる。幸いこの内殻補正にとってはこの ρc (r) を忠実に再現することは必要ない。式
（1.37）でその非線形効果が著しいのは、ρc (r) と ρv (r) とが同じ程度の大きさの領域でだ
けである。ρc (r) が非常に大きくなる原子核のごく近傍ではその補正は重要でない。それ
ゆえそのような原子核のごく近傍の ρc (r) はそれを緩やかな関数に置き換えてもそれほど
問題ない。
そこで半径 rcc を境に以下のような近似を用いる。[25]

A sin(Br)/r if r < r
cc
(1.39)
ρcpartial (r) =
ρc (r)
if r > rcc
f ull
A、B は解析的に接続という条件より定める。
この部分補正のやり方として、近似式（1.39）よりもさらに改良されたものが提案され
ている。[26]
1.8
将来展望
… 未来のことは誰にもわかりませんが。
擬ポテンシャルで、遷移金属などを扱おうとすると、次第に内殻状態の微妙な影響を取
り入れざるを得ず、擬ポテンシャルを生成する技巧がますます複雑になってきたことを見
てきた。「OSAKA2000」の擬ポテンシャルは基本的にノルム保存型を扱っているが、こ
れとは異なったアプローチとしてウルトラソフト擬ポテンシャルというものが提案された
[27]。これは元々、擬ポテンシャルの収束性を向上させるためノルム保存という拘束を破
棄する方法だが、それによる transferability の低下を回避するため、殻内の電子密度によ
りポテンシャル補正するものである。これは一見、第 1.7 節の内殻補正 NLCC を思わせる
が、補正としては大きく違う。NLCC が内殻電荷密度を参照原子に対して一度求めたら固
体でもそれを固定して用いるのに対し、ここでの殻内の電子密度による補正では、殻内電
16
子密度と擬ポテンシャルは一体のものとして、それぞれの環境下でセルフコンシステント
に計算される。つまり、ノルム保存型では、ポテンシャルを生成する部分と KS 方程式を
解く部分は分離できたが、ウルトラソフト型では再びそれらは同時に解くものとみなされ
る。こうしてこの点では全電子計算と同じ格好で、ポテンシャルと電荷密度分布がセルフ
コンシステントに解かれるので、スピンの入った微妙な効果を扱うには向いているといえ
る。このウルトラソフト型を用いて、平面波法で扱うのが一番遅れたランタノイドやアク
チノイド化合物にも成功裏に応用されている [28]。
このように擬ポテンシャルの未来は再びポテンシャルと KS 方程式の解は一体として扱
われるようになり、ノルム保存型が滅びるのか、それともノルム保存型の中で、スピンの
効果も適切に取り扱えるようになれるのか筆者には予測できない。
この他、密度汎関数理論の枠内で、自己相互作用補正、非局所密度汎関数への拡張など
一言述べるべきテーマが多々あるが、この辺で。
17
第2章
2.1
原子擬ポテンシャル、atom の使い方：
スピンなし
擬ポテンシャル生成
ポテンシャル生成のプログラム atom の使い方については、上巻で簡単に触れましたが、
ここではさらに詳しく述べます。
まず、インプットデータファイル atom.dat のフォーマットは表のようになります。
Table 2.1: atom.dat のフォーマット
itype
ititle
ikerk
nameat icorr + ispp
znuc
zsh
ncore
nval
n
l
...
rcs
rcp
rsh
rmax
aa
zo↓
zo↑
evd
rcd
rcf
cf ac
bb
rcf ac
それぞれの行の意味を書くと、
一行目 計算種（itype）とタイトル（ititle）。計算種としては itype=ae: 全電子計算、
pg: 擬ポテンシャル生成、pe: 交換項に関する内殻補正をした擬ポテンシャル生
成、ph: ハートリー項まで含めた内殻補正をした擬ポテンシャル生成、pt: 擬ポテ
ンシャルテスト、pm: 擬ポテンシャルのテストと価電子の修正。
二行目 擬ポテンシャルの種類（ikerk）。どんな種類があるかは第一巻参照。
三行目 原子名（nameat ）と計算オプション（icorr + ispp）。icorr は電子相関汎関数
の種類。ispp に s の指定があればスピンを含む、r であれば相対論的計算。
四行目 核の電荷数（znuc）、内殻の電荷数（zsh）、内殻半径（rsh）、最大半径（rmax）、
動径方向の刻みのパラメータ（aa、bb ）。
五行目 内殻軌道の数（ncore ）と価電子軌道の数（nval ）。
六行目以下 n と l で指定される価電子軌道ごとの電子占有数（zo）。オプションとして
軌道エネルギー（evd）。
18
最終行 擬ポテンシャル生成のカットオフ半径（rc ）。およびオプションとして、内殻補
正パラメータ（cf ac と rcf ac）。
原子軌道はスピンを考慮しなければ、本来は主量子数 n と角運動量子数 l と共に、方位量
子数 m も与える必要があると疑問に思うかもしれません。開殻原子に対しては本来その通
りですが、このプログラム（他のものもほとんどそう）では、ポテンシャルは電子の寄与
も含めて完全な球対称を仮定しております。ですので、三つある p 軌道、px 、py 、pz に電
子が一個しかない場合でも、電子は例えば px だけを占めていると考えるのではなく、px 、
py 、pz それぞれを平等に 1/3 づつ占めると仮定しております。でなければ問題は本質的に
３次元問題となり解くことが急速に難しくなるでしょう。従って m は考慮しません。
全電子計算は、局所密度近似による動径方向シュレーディンガー方程式を解くことに
より、各固有値、軌道を得ます。プログラム atom のこの部分は、適当な境界値の元での
微分方程式の固有値を求めることに相当し、そのプログラムアルゴリズムの相当部分は、
Koonin の名著 [?] の第２章、３章の部分に相当します。プログラムの解読を試みられる方
は是非一読を勧めます。
四行目のパラメータは動径方向のグリッドを計算するためのものです。それは等間隔で
はなく、指数関数的なグリッドを使います。すなわち、
ri = a(exp{b(i − 1)} − 1)
(2.1)
で、この式の中の a や b はそれらの入力力計算されます。入力値を 0.0 にすれば、プログ
ラムが原子番号の大きさによってそれらのデフォルト値を与えます。通常はそれで十分で
す。
具体例として、C 原子を取り上げます。
pg
Carbon
tm2
n=C c=ca
0.0
0.0
0.0
1
2
2
0
2.00
0.00
2
1
2.00
0.00
1.25
1.25
0.00
0.00
0.0
0.0
0.0
と入力します。以下のようにポテンシャルが生成されます。
ここで出力ファイル atom.out をチェックしてみましょう。タイトルの後、
correlation = ca
nonspin-polarized
nuclear charge
number of core orbitals
number of valence orbitals
electronic charge
ionic charge
=
=
=
=
=
6.000000
1
2
6.000000
0.000000
と核の電荷、内殻軌道の数などの情報が打ち出されます。価電子の数は通常は中性原子の
ものを使いますが、イオン状態のものも計算できます。その時はその差がイオン電荷とし
て現れます。その後、
19
i
1
2
3
n
1
2
2
l
0
0
1
s
0.0
0.0
0.0
j
0.0
0.0
0.0
occ
2.0000
2.0000
2.0000
と、n、l で指定される軌道の占有数が打ち出されます。今の場合スピンなしですので、s、
j の量子数は関係ありません。動径方向のグリッドパラメータ
radial grid parameters
r(1) = .0 , r(2) =0.52E-05 , ... , r(1007) =119.49
a = 6.00 b = 80.00
が表示された後、プログラムは全電子計算に入ります。
全電子計算ではまず、各軌道に関する情報が示されます。例えば
n = 2
l
a
r
r
r
= 0 s = 0.0
extr
0.352
extr
0.139
zero
0.379
90/99 %
2.542
-0.815
1.216
3.889
は 2s 軌道に関するもので、波動関数のピークが二つでその波動関数の値が a extr で、その
位置が r extr で示されます。また波動関数が 0 の点が r zero で示されます。 r 90/99 %は
波動関数 |u(r)|2 の原点からの積分値がそれぞれ 90 および 99 %になる位置を示しており
ます。
その後、
C output data for orbitals
---------------------------nl
1s
2s
2p
s
0.0
0.0
0.0
occ
2.0000
2.0000
2.0000
eigenvalue
-19.89570500
-1.00194941
-0.39859847
kinetic energy
31.42793836
3.29299206
2.46684613
と各軌道のエネルギーが打ち出され、その全エネルギーの成分、
total energies
-------------sum of eigenvalues
=
-42.59250576
kinetic energy from ek
=
74.37555309
el-ion interaction energy =
-175.02012967
el-el interaction energy =
35.24955812
vxc
correction
=
-12.44704330
virial correction
=
0.47288553
exchange + corr energy
=
-9.45350385
kinetic energy from ev
=
74.37555098
potential energy
=
-149.22407541
--------------------------------------------total energy
=
-74.84852443
そしてビリアル定理の評価が続き、全電子計算が終了します。
次に擬ポテンシャル生成の過程に入ります。
20
pot energy
-67.23756234
-11.01409737
-9.25840513
C
pseudopotential generation using the Improved Troullier and Martins method
のタイトルの後、擬ポテンシャルに関するデータが出力されます。
nl
1s
2p
s
0.0
0.0
eigenvalue
-1.001949
-0.398598
rc
1.246587
1.246587
cdrc
0.360770
0.337930
delta
-0.488011
1.009600
と各準位、カットオフ半径がでます。このカットオフ半径は入力したときの値と微妙に違
うことがありますが、それは入力値に対してプログラムが若干その半径を修正しているの
です。cdrc はカットオフ半径内のその軌道の全電荷密度。delta は擬波動関数の規格化に
関するパラメータです。
その後、全電子計算の場合と同じく、各軌道の波動関数の最大値、その位置、零点など
その軌道に関する情報が打ち出されます。
最後に、得られた擬ポテンシャルが適切なものであるか、改めてその原子について解い
てみてテストします。
test of eigenvalues
の見出しの後、
rcut =
eval =
s
1.23
-1.00193
p
1.22
-0.39858
とカットオフ半径および固有値が打ち出されますので、それらが前に計算したものと一致
するかチェックします。
各擬ポテンシャルは、原点での値、および最小値とその位置
l
s
p
vps(0)
-1.890
-22.400
vpsmin
-8.558
-22.400
at r
0.923
0.000
によりだいたいの特性がわかります。
もう一度各軌道のエネルギー、占有数が
C output data for orbitals
---------------------------nl
1s
2p
s
0.0
0.0
occ
2.0000
2.0000
eigenvalue
-1.00193167
-0.39858306
kinetic energy
0.83899878
2.37880778
pot energy
-5.38595922
-6.17121542
のように計算されます。これより各軌道エネルギーは正しく再現されていますので得られ
た擬ポテンシャルは適切であることがわかります。最後に全エネルギーの内訳を打ち出し
て締めくくっています。
2.1.1
出力ポテンシャルデータファイル
生成された擬ポテンシャルのデータは前バージョンではファイル fort.10 に出力され、そ
れをリネームして使うと述べました。
「OSAKA2000」ver. 2 ではポテンシャルをスピンや、
相対論的にも拡張すべく、ファイルおよびフォーマットを変えました。擬ポテンシャルファ
21
イルは pseudo.dat#に出力されるものを使います。これはテキストファイルでなく、バイ
ナリーファイルです。また擬波動関数もバイナリーファイルとして別のファイル（fort.13）
に落としました。
pseudo.dat#どのようなデータが含まれているか説明しておきます。
最初に
nameat, icorr, irel, nicore,
atom-lda
2-JUN- 1
Improved Troullier - Martins potential
2s( 2.00) rc= 1.252p( 2.00) rc= 1.25
npotd, npotu, nr-1, a, b, zion
が記録されます。原子名の後、電子相関汎関数のタイプ icorr、スピン条件 irel、内殻補
正タイプ nicore が続きます。
icorr はページ 49 で述べられた通りです。irel はスピンがある場合 isp、相対論的な場
合 rel、それ以外では nrl となります。
nicore は非線形内殻補正のタイプで、四文字の変数です。最初の一文字は f か p でそれ
ぞれ完全、擬内殻補正を意味し、後の三文字は cec が交換相関項のみ、che はハートリー
項も含めた補正を意味します。非線形内殻補正がなければここには nc が入ります。
二行目はそのまま書かれているようなタイトル、三行目が、価電子のカットオフ半径お
よび占有数のデータ、四行目はダウン、アップのポテンシャルの数（npotd,npotu）、動径
方向グリッドの数 nr、式（2.1）のグリッドパラメータとなります。最後がイオンの電荷
zion です。
次に、動径方向グリッドの値 ri が i= 2 から nr まで打ち出されます。それからダウン、
アップのポテンシャルが npotd,npotu の数だけ打ちだされます。今の例では、ポテンシャ
ルはスピンが入っていませんので、どれか一つでよく、便宜上スピンダウンの部分に書か
れます。第??章で説明されるスピンがある場合は、スピンダウンとアップの両方にポテン
シャルが書かれますが、しかしそこで説明されるように最終的なイオンポテンシャルはス
ピンに関して平均化したものを出力しますので、両方とも同じ値です。相対論の場合は、
スピンダウンの部分に l に関する平均化ポテンシャルが、スピンアップにはスピン軌道相
互作用ポテンシャルが書き込まれます。
最後に、非線形内殻補正がある場合には内殻電荷密度 ρc に続き（内殻補正が無い場合
は全て 0）、価電子密度 ρv が打ち出されます。
2.2
参照原子の価電子と電子配置の取り方
ここまでは第一巻で述べられたことの繰り返しですが、これから s、p 価電子だけの典
型的な半導体以外の場合について述べる必要があります。s、p 価電子を持つもの以外の原
子を考える場合、まず行き当たる問題はどこまで価電子とみなすかということです。それ
に対する答えとして、最初に考えるのは原子の基底状態で開殻にある軌道でしょうが、そ
れでは不十分な場合がたくさんあります。
また式（1.7）の局所ポテンシャルをどう選ぶかも問題です。擬ポテンシャルの l に関す
る和が無限であれば、原理上どうとっても良いはずですが、実際には分解した l 成分数は
もちろん有限ですので、局所ポテンシャルの取り方により結果は多少なりとも異なります。
22
ことに Kleinman-Bylander 型ですと、それは元々近似式ですので、局所ポテンシャルの取
り方は重要ですので、次節で述べます。
どこまで価電子として含めるか、またそのうちどれを局所ポテンシャルとするかは、最
終的には、選択した l のセットよりも大きな l 成分を加えても結果が変わらないようなも
の、という基準で計算実験を通じて決まります。これでは「非経験的計算」とは言えない
という批判があると思いますが、（確かに今でもその取り方には議論はあるものの）だい
たいはこれくらいの取り方をすればよいという基準があります。
まずは、全原子に対してシステマテックな擬ポテンシャル作成のレシピを与えてある
Bachelet、Hamann、Schlüter の文献 [7] が参照になると思いますので、価電子の取り方は
だいたいそれに準じております。しかし軽元素に対して、その文献では d 電子も入れてお
りますが、実際には d 電子を入れなくとも計算結果にはほとんど影響しませんので、この
プログラムではほとんどの場合無視しております。こうして決められたデータが原子番号
順に psmap.dat に与えられております。但し局所ポテンシャルの選び方に関しては、この
文献 [7] では Kleinman-Bylander 型のことは考慮しておりませんので、次節で述べられる
ことを参考にして Kleinman-Bylander 型でも有効なように設定しております。
psmap.dat の読み方は、例えば原子番号 29 の銅のデータを見ますと、
29
2
0
11
0
3
1
Cu
とでております。
第一行目 は、原子番号 Z 、中性原子の価電子数 Nv 、擬ポテンシャルを作るべき軌道角運
動量の種類数 lmax 、原子記号となります。Cu の場合、基底状態の開殻構成は 3d10 4s1
ですので、価電子数は 11 です。擬ポテンシャルは 3d、4s、4p 軌道に対して作りま
す。つまり lmax = 3 です。
第二行目 が lmax 個ある l を順番に並べております。
第三行目 に局所ポテンシャルに選ばれた l = lloc が書かれております。
このようにして擬ポテンシャルの種類を指定したら、次に問題となるのは擬ポテンシャル
を生成するときの参照原子の電子配置です。これも Bachelet、Hamann、Schlüter の文献
[7] を参照にしましたが、本来このようにして作られた擬ポテンシャルが真に transferable
なものであれば、それは「作られたときの電子配置に依存しないはず」です。擬ポテン
シャルのオリジナルの論文を書いた Hamann、Schlüter、Chiang の論文 [6] でも、擬ポテ
ンシャルを作るときの電子配置を変えてみても、軌道エネルギーは全電子計算との相対誤
差が 0.1 %以下と素晴らしい一致を見せています。それゆえ原子の基底状態とそれほど大
きく異ならない限り、電子配置の問題は少ないと言えます。但しスピンが入った場合は少
し問題で、それに関しては第??章を参照して下さい。
例えば Cu 原子については（ファイル Cu.dat 参照）、価電子を s、p、d と取り、
pg
Copper
tm2
n=Cu c=ca
23
0.0
5
3
3
2
4
0
4
1
2.08
0.0
10.00
0.75
0.25
2.30
0.0
0.0
0.0
0.0
0.00
0.00
0.00
2.08
0.00
のようなデータを与えました。これにより d、s、p の擬ポテンシャルが作成されます。自分で
価電子の取り方を変えるときは、この atom.dat と psmap.dat がコンシステントとなるよう
に注意して下さい。擬ポテンシャルの並び方はこのデータですと（ファイル fort.10 の中で）、
d、s、p の順となります。この順番を上で述べた psmap.dat の第二行めの 2
0
1 の
並びと一致させます。
さらに、価電子として空の状態も指定できます。つまり、ncore 、nval 以下の行を
5
3
4
4
2.08
3
2
0
1
10.00
1.00
0.00
2.30
0.00
0.00
0.00
2.08
0.00
とすることもできます。このとき 4p は空の準位で、基底状態の電荷密度には寄与しませ
んが、原子に対するシュレーディンガー方程式を解くことでその波動関数を得ることがで
きます。但しそれが束縛状態になるとは限りません。擬波動関数が束縛状態でないと、擬
ポテンシャルを作ることに失敗しますので、このときは電子配置などの入力パラメータを
変えます。また、一度に全ての軌道に対する擬ポテンシャルを作る必要はありません。ま
ず、s と d の擬ポテンシャルを原子の基底状態の電子配位で求めて、それから p に対して
は別の電子配位で求めても構いません。
最終的には擬ポテンシャルのファイル*.pot（通常は fort.10 をそのままリネームしたも
の）は以下のようなフォーマットを持つようにします。
nr
}
r
u
pseudowavefunction s
···
marker wsp
}
r
u
pseudowavefunction p
···
marker wpp
}
r
v
pseudopotential s
···
marker vns
}
r
v
pseudopotential p
···
marker vnp
ヘッダーに動径方向の分割数 nr が書かれた後、波動関数、ポテンシャルの順に出力され
24
ます。
最後に各々の l 成分に対して、擬ポテンシャルのカットオフ半径データを与えなければ
なりません。第一巻ではカットオフ半径の与え方は最外殻ピーク位置よりも 1.2 1.6 くら
いを取ればよいと述べました。それより内側でも構いませんが、最外殻零点よりは大きく
なければなりません。実際には 3d 電子などは非常に核の近くに局在化しているため、カッ
トオフ半径はそれよりもかなり大きめに取らないと、Troullier-Martins 型でポテンシャル
に課せられている多くの条件を満たすことができなくなります。このあたりはいろいろ経
験が必要なところだと思いますが、筆者自身も 100 種に余る原子に対していろいろ試行し
ているわけではありません（そんなことは不可能です）。それで具体的値はユーザーに委
ねたいところですが、しかしともかくも何からの初期値くらいは準備しないと使えないと
思い、100 種に余る原子の atom.dat を用意しました。それらはディレクトリ atomscript
にあります。
もう一度強調しておきますが、これらのデータは「怠け者」のために用意したものであっ
て、いかなる意味においても値の的確さなどを保証するものではありません。あくまでも
ユーザー自身が調節、チェックして下さい。
2.2.1
いくつかの注意すべき例
水素原子に対しては、擬ポテンシャルでなく、裸のイオンポテンシャルを用いています。
2.3
局所成分の取り方
ここでは Mo を例に取り、局所ポテンシャルの取り方の影響を調べます。psmap.dat を
開けてみると、Mo は
42
2
0
6
0
3
1
Mo
となっております。すなわち擬ポテンシャルは、３つの d、s、p 軌道成分に対して作られ
ます。それらの動径依存性は図 2.1 に示されます。
これから d 成分は深く、これを局所成分とすると、s と p の ∆Vl は正に大きくなります。
これはページ 8 の Kleinman-Bylander の条件 2 に反することとなり、避けるべきです。
一方、s 成分を局所成分と選ぶと、∆Vd は負に大きく、また ∆Vp は正になったり負に
なったりすることがわかります。後者の場合 hu0l |∆Vl | u0l i がほとんど 0 になり人為的に大
きな KB エネルギーを与える可能性があります。psmap.dat での局所ポテンシャルの指定
は lloc = 0 を選択しております。
これにより作られたポテンシャルを Kleinman-Bylander 型のポテンシャルを使わず、元
の半局所的ポテンシャルのまま固体の Mo に適用して得たバンド図を図 2.2 に示します。
これから局所ポテンシャルとして lloc = 0 と取ろうが lloc = 2 と取ろうが結果はほとんど
変わりません。局所成分をどちらに取ろうが、電子は最終的には局所成分、非局所成分両
25
[Ry]
0
-5
-10
V
d
s
p
-15
-20
-25
0
1
2
3
4
5
r
図 2.1: Mo の擬ポテンシャル
方の影響を受けることの反映です。
一方、Kleinman-Bylander 型ポテンシャルを用いると、図 2.3 に示されますように劇的
な変化が生じます。lloc = 0 としたときは前の半局所ポテンシャルの場合と同じ結果を与
えますが、lloc = 2 と選んだ場合、−2 Ry くらいにほとんど分散の無い新たな準位ができ
ていることがわかります。この場合にも電子にかかるポテンシャルは局所と非局所の合わ
さったもののはずですが、第 1.2 節で述べられたように Kleinman-Bylander 型の場合は非
局所ポテンシャルが射影演算子で与えられますので、その演算子の張る部分空間と直交す
る波動関数に対しては実質上何の影響も与えなくなります。このとき電子は深い局所ポテ
ンシャルの影響のみを受け深いゴースト準位を作ります。
このゴースト準位を作る波動関数がどうなっているかを見てみましょう。図 2.4 に固体
Mo における Mo 原子の周りのゴースト準位の波動関数が示されております。波動関数は
第 1.2 節で述べられたような節を持つことが見て取れます。
このときの KB エネルギーや KB コサインの値をチェックしてみましょう。それは inip
の出力ファイル inip *.out に示されます。
==================== POTENTIAL TEST =======================
の見出しの後にリストされます。lloc = 0 の場合は、
---------- VKBNL ----------1 th atom =
mo
l
orb
sign
<u|DV|u>
2
d
-1
-4.1724
1
p
1
0.0212
EKB
-9.0109
0.9803
26
<u|DV^2|u>
37.5967
0.0207
COS_KB
-0.6805
0.1469
(a)
1
1 4
1 1
4
5
1
1
Energy \Ry\
4
1.5
1.0
1
7 4
1
2 5
2 1
3 2 2
1 1 1 1
6 4
4
4
4 4
1
5 3 3 3
1 3
1
1
1
3
1
3
1
1
3
3
1
3
3
1
3
9
1
3
7
1
3
4
5
1
Λ
N D P
4
1
4
5
5
4
3
5
3
1
1
1
0.5
0.0
(b)
Mo PseudoPotential
2.0
1
3
4 1
1
6 3
3
3
7
7
2 2 2 2
3
4
1
4
3
3
1 1 1
1
4 1
6
1
3
4
3 4
5
2
3
1 1 1
9 3
1
1 1
2.0
1.5
4 3
3
1 1 1
3 9
4
1
4 7
1 1
2 2
2
1
1
1
∆
Γ
G
H
N
Σ
EF
1.0
1 1
Mo PseudoPotential
1 3
8 3 1
7 4
1
1
4
1
1
5
2 2
3 2 2
1
1 1 1 1
6 4
4
4
4 4
3
5 3 3 1
3 1
1
1 1
1
1
3
3
1
3
3
1
3
3
3
1
0.0
1
3
9
1
3
4
5
7
3
1
1
3
3
4
3
1 4
3 3
1 1
1
4
1
1
2
2 2
1 1
1
1
Λ
N D P
5
1
1
4
5
4
4
5
3
1
0.5
Γ
1
8
1
4
1
6 3 3
3 2
7
2 2 2 3 7
1
4
3
1
1
1
4
6
1
3
4
3 4
5
3 1
1 1
3 1
9 1
5
1 4
4 1 1
∆
Γ
H
G
N
1
Σ
3
1 9
EF
41 7
2
1 1
Γ
図 2.2: Mo のバンド図。（a）は lloc = 0 の場合で、（b）は lloc = 2 の場合。
Mo PseudoPotential
(a)
Energy \Ry\
2.0
4
1
4 1 1
7 4
1.5
4
3 2 21 21 21
1 1
1
1
1
1
5
3
1
3
6 4
1.0
4
4
4 4
5 3 3 3 13
1
1
1 1
1
3
1
4
3
1
31
3
9
1
3
3
1
3
7
45
1
0.5
5
4
53
3
1
1
1
1
4
5
31
67 3 3
2 2 23 2 7
3
4
41
3
1 1 31
41 1
4 3 6 1
3 4
5 2
3
9 31 1 1 1
1 1
1.5
1.0
43
1 31 3
1 9
4
1
1 124 7
2 2
1
1
N D P
Λ
Γ
EF
H
G
N
Σ
1
3
13
1
3
31
65 34 43 43 43 4
1 1
1 1 1
31
3
1
31
31
3
1
3
3113
9
67
1
53145
1
1
1
1
1
1
1
541
53
1
41
51
5
4
1
51
5
3
1
31
1
1
4 2 2 2
2
3 3 3 3
41 41 1 31
431 31 43 41
4 1
4
4
3 3
9 31 31 1 1
7
6
6
7
8
3
1
864
3
4
31 43
1 3 3 9
413 41 141 1
3 43 4314 67
21
65 21 21 21
1 1 1
1 1 1
0.0
-0.5
-1.0
1
1 1
∆
7 4 3
8 3 1 1 1
3 2 21 21 21 5
1 1
864 431 431 314 431 4
0.5
-1.5
-2.0
0.0
Mo PseudoPotential
(b)
2.0
4
1
Γ
-2.5
1 1 1 1 1 1
N D P
1
1
1
Λ
Γ
1
∆
1
H
1 1 1 1 1 1
1 1 1 1
G
N
Σ
Γ
図 2.3: Mo のバンド図。Kleinman-Bylander 型を用いた場合。（a）は lloc = 0 の場合で、
（b）は lloc = 2 の場合。
となっております。各行は l とその名前の後、hu0l |∆Vl | u0l i、ElKB の符号、ElKB 、∆Vlrms 、
及び式（1.12）の ClKB をリストしています。これから d 軌道成分は負に大きい ElKB を持
ちますが、ClKB の大きさは小さくなく問題とならないことがわかります。
一方、lloc = 2 の場合は、
l
0
1
orb
s
p
sign
1
1
<u|DV|u>
0.3105
0.1020
EKB
5.4343
3.1387
<u|DV^2|u>
1.6874
0.3203
COS_KB
0.2390
0.1803
となり、s、p 軌道成分とも正に大きい ElKB を持つだけでなく、ClKB が小さくなっており、
危険であることを示しております。
このように、ElKB と ClKB を見比べることで大体の見当がつきます。大きい小さいの基
準は、ElKB に関しては問題とする価電子バンドの幅に比較して、ClKB に関しては文献 [17]
によれば 0.15 くらいを目安と考えれば良いでしょう。この大きさはしかしながら扱うポ
テンシャルの型で多少違います。文献 [17] では擬ポテンシャルとして Bachelet-HamannSchlüter のものを採用しております。「OSAKA2000」では Troullier-Martins 型を採用し
27
ϕ
r
図 2.4: Mo 結晶における、lloc = 2 と選んで Kleinman-Bylander 型を用いたときのゴース
ト準位の波動関数。Mo 原子を中心に波動関数が深く局在している様子を示したもの。
ておりますので、この条件の値は多少違うと見るべきです。
スピンなしの擬ポテンシャル作成のまとめ
価電子の与え方、及び局所成分の取り方は psmap.dat に与えられてい
る。個々の原子の擬ポテンシャルのカットオフ半径はデーターベース
atomscript の中にある該当する原子の atom.dat を持ってくる。作られ
た擬ポテンシャルが Kleinman-Bylander 型として用いても安全であるか
どうかは inip プログラムの中で簡便的にテストする。
Kleinman-Bylander 型の安全性の本当のテストには、固体のバンドを書いてみるなり、孤
立原子（大きなセル中に一個の原子を置く）での固有値を求めてみることが必要です。
28
第3章
原子擬ポテンシャル、atom の使い方：
拡張
スピンの擬ポテンシャル生成
3.1
スピンがある場合の擬ポテンシャルは第 2 章と同じように C 原子を具体例として取り上
げながら説明する。
pg
Carbon
tm2
n=C c=cas
0.0
0.0
0.0
1
2
2
0
1.00
1.00
2
1
2.00
0.00
1.25
1.25
0.00
0.00
0.0
0.0
0.0
と c=のところにオプション s を付け加えることで、スピン効果を含めることを宣言する。
価電子配置をどうするかが問題だが、ここではまず原子の基底状態 2s{ 12 ↓, 21 ↑}2p{ 22 ↓
, 0 ↑} を選んでいる。
出力ファイル atom.out をチェックする。タイトルの後、
correlation = ca
spin-polarized
とスピン分極の場合となっていることを確認する。
電子配置が
i
1
2
3
4
5
6
n
1
1
2
2
2
2
l
0
0
0
0
1
1
s
-0.5
0.5
-0.5
0.5
-0.5
0.5
j
0.0
0.0
0.0
0.0
0.0
0.0
occ
1.0000
1.0000
1.0000
1.0000
2.0000
0.0000
と、意図したものとなっているかどうかもチェックする。今度はスピンありなので、s が
ダウン、アップの順で埋めらる。
各軌道は n、l、s の組で指定される。以下全電子計算に入り各軌道の波動関数形状のデー
タが出力された後、軌道エネルギーが打ち出される。
nl
1s
1s
2s
2s
2p
2p
s
-0.5
0.5
-0.5
0.5
-0.5
0.5
occ
1.0000
1.0000
1.0000
1.0000
2.0000
0.0000
eigenvalue
-19.88158501
-19.81461085
-1.05891709
-0.87736265
-0.45168782
-0.28505567
kinetic energy
31.35723528
31.50593322
3.37159475
3.10801273
2.56725092
2.21056515
29
pot energy
-67.16015808
-67.32282483
-11.13670561
-10.71746050
-9.46881369
-8.70710876
今の場合、スピンが分極しているのは最外殻の 2p 軌道だけだが、それがスピン分極して
いるため、交換相関ポテンシャルのため、自身の軌道エネルギーばかりでなく、閉殻であ
る内殻電子軌道 1s でさえもわずかながらスピンアップとダウンではエネルギーに差が生
じている。1 ちなみに全エネルギーは
total energy
=
-74.93147822
とスピン分極なしの場合の値 −74.84852443（ページ 20 参照）に比べて下がっていること
がわかる。
次に擬ポテンシャル生成の過程に入る。
C
pseudopotential generation using the Improved Troullier and Martins method
のタイトルの後、擬ポテンシャルに関するデータが出力される。
nl
1s
1s
2p
2p
s
-0.5
0.5
-0.5
0.5
eigenvalue
-1.058917
-0.877363
-0.451688
-0.285056
rc
1.246587
1.246587
1.246587
1.246587
cdrc
0.367373
0.344094
0.352488
0.301124
delta
-0.477891
-0.516962
1.031753
0.943877
とスピンに関しては同じカットオフ半径が採用される。
最後に、得られた擬ポテンシャルが適切なものであるか、改めてその原子について解い
てみてテストする。このときのポテンシャルは、スピンがアップとダウンのもので占有数
での平均化を行ったもの式（1.27）を使っている。
今の場合は、
l
s
s
p
p
vps(0)
-1.901
-1.901
-22.474
-22.474
vpsmin
-8.557
-8.557
-22.474
-22.474
at r
0.923
0.923
0.000
0.000
となって、当然スピンアップとダウンでは同じとなっている。
電子配置は上の例では、原子基底状態 2s{ 21 ↓, 12 ↑}2p{ 22 ↓, 0 ↑} を考えてきましたが、非
分極状態 2s{ 12 ↓, 12 ↑}2p{ 22 ↓, 22 ↑} を考えることもできる。このときは対称性から、第 2.1
節のスピンを考慮しない場合の結果が再現されるはずだ。
さらに完全分極状態 2s{ 22 ↓, 0 ↑}2p{ 32 ↓, 0 ↑} を考えることもできる。どれをとってもス
ピン平均化ポテンシャルはほとんど変わらないはずだ。しかし完全分極状態の場合を実行
してみると（全エネルギーは −74.63751327 とスピン分極なしの場合の値に比べて少し上
がる）、擬波動関数を計算する過程で、2p↑ 状態が、
WARNING:
This orbital is not bound!
と束縛状態でないという警告を出す。全電子計算ではこの状態は負のエネルギーつまり束
縛状態となっているので、ここで生成された擬ポテンシャルの transferability は悪くなっ
1
このように最外殻のスピン状態が内殻のスピン分極に跳ね返る効果は一般に存在することになるが、特
に遷移金属ではこの効果が大きくなるため、それらの原子に対しては凍結内殻（frozen core）近似が適切で
なくなる。従ってこれらの原子は擬ポテンシャルで扱うことが一番難しいものとなっている。
30
たと考えられる。
いずれにせよ、
「OSAKA2000」では、スピンを含む場合、式（1.27）で与えられるスピ
ンに関する平均化したものをイオンポテンシャルとして用いている。従って価電子のスピ
ン分極は、純粋に電子の交換相関ポテンシャルからのみ生じることになる。
3.2
内殻補正
非線形内殻補正を行うためには、まずポテンシャル生成プログラム atom で、入力パラ
メータ itype を pe か ph と選ぶ。
内殻補正の仕方にはいくつかの選択がある。
完全内殻補正 得られた内殻電子密度をそのまま使う。cf ac < 0 か価電子が無い場合に
相当。
擬内殻補正 ある半径 rc 以内の内殻電子密度 ρc を、適当な滑らかな関数、Ar sin(Br)、
に置き換える。半径 rc の外では真の内殻電子密度を使う。このコア半径 rc は rcf ac
（正の数字）で与えられる。rcf ac < 0 とすると、内殻電子密度と価電子密度の比が
cf ac になる位置となる。
3.3
相対論的効果
相対論的効果を取り入れる場合の擬ポテンシャルは第 2 章と同じように C 原子を具体例
として取り上げながら説明する。
pg
Carbon
tm2
n=C c=car
0.0
0.0
0.0
1
2
2
0
1.00
1.00
2
1
2.00
0.00
1.25
1.25
0.00
0.00
0.0
0.0
0.0
と c=のところにオプション r を付け加えることで、相対論的効果を含めることを宣言する。
価電子配置をどうするかが問題でが、ここでは第??章と同様、とりあえず原子の基底状
態 2s{ 12 ↓, 12 ↑}2p{ 22 ↓, 0 ↑} を選んでみる。このような選び方をしても、atom は同じ n, l
に属するものはスピンダウン、アップの占有数の合計、zo = zo↓ + zo↑ に対して、ヘリシ
l+1
l
ティが負（κ < 0）のとき 2l+1
zo、正（κ > 0）で 2l+1
zo を割り振る。
出力ファイル atom.out をチェックしてみる。タイトルの後、
r e l a t i v i s t i c ! !
correlation = ca
spin-polarized
と相対論的計算になっていることを確認する。
電子配置は今度は
31
i
1
2
3
4
n
1
2
2
2
l
0
0
1
1
s
0.5
0.5
-0.5
0.5
j
0.5
0.5
0.5
1.5
occ
2.0000
2.0000
0.6667
1.3333
となる。各軌道は n、l、s、j の組で指定される。本来 l、s は良い量子数ではないが、簡
便上それらを使って出力する。正しくは (n, j, sign(κ)) の組で示されるべきでが。それだ
と、上の 4 つの軌道は、(1 12 ↑)、(2 21 ↑)、(2 21 ↓)、(2 32 ↑) となる。
以下全電子計算に入り各軌道の波動関数形状のデータが出力された後、軌道エネルギー
が打ち出される。
nl
1s
2s
2p
2p
s
0.5
0.5
-0.5
0.5
occ
2.0000
2.0000
0.6667
1.3333
eigenvalue
-19.89222766
-1.00237937
-0.39887080
-0.39821874
kinetic energy
0.00000000
0.00000000
0.00000000
0.00000000
pot energy
-67.27361940
-11.02195981
-9.26468624
-9.25605977
こうして (2 12 ↑) と (2 32 ↓) は同じ l 軌道から作られるが、j が違うためエネルギーが分裂す
ることがわかる。しかしその分裂は非常に小さく C 原子では相対論的効果は無視できるこ
とがわかる。
ちなみに全エネルギーは
total energy
=
-74.86537517
とスピン分極なしの場合の値 −74.84852443（ページ 20 参照）に比べてわずかに下がって
いることがわかる。
次に擬ポテンシャル生成の過程に入る。
C
pseudopotential generation using the Improved Troullier and Martins method
のタイトルの後、擬ポテンシャルに関するデータが出力される。
nl
1s
2p
2p
s
0.5
-0.5
0.5
eigenvalue
-1.002379
-0.398871
-0.398219
rc
1.246587
1.246587
1.246587
cdrc
0.361061
0.338271
0.337766
delta
-0.486426
1.011140
1.009403
とスピンに関しては同じカットオフ半径が採用される。
最後に、得られた擬ポテンシャルが適切なものであるか、改めてその原子について解い
てみてテストする。
今の場合は、
l
s
p
p
vps(0)
-1.905
-22.418
-22.402
vpsmin
-8.558
-22.418
-22.402
at r
0.923
0.000
0.000
32
となっている。第??章のスピンの場合と違って、相対論的な場合は、スピンのアップ、ダウ
ンでイオンの裸のポテンシャルは違う（今の例ではその違いは極くわずかに過ぎないが）。
atom はその後、イオンポテンシャルとして、各 l に対して、式（7.6）で与えられる平
均化ポテンシャルと、式（7.7）で与えられるスピン軌道相互作用のポテンシャルをこの順
に排出する。l = 0 の場合だけ平均化ポテンシャルのみとなる。スカラー相対論ではこの
平均化ポテンシャルだけでよいことになる。
3.4
GGA 補正
今回 GGA を Osaka2k に取り入れたので、その使い方およびテストを報告する。GGA を
取り入れた擬ポテンシャル作成の実際的なコードとしては、Martins らが従来の TroullierMartins 型のポテンシャル作成プログラム atom を発展させた atm を使っている。バージョ
ンは 5.693 のパッケージを使っている。
3.4.1
移植の原理・方針
移植したのは GGA としては最も広く用いられている PBE96 のスキムである。[?]
またこの GGA を平面波基底に移植する際、良く用いられている White-Bird のスキム
[?] を採用している。このスキムによると計算ステップは
1. 密度の変換: n(R) → n(G) (1×FFT)
∇n(G) = Gn(G) の計算
2. 密度勾配の実空間への戻し: Gn(G) → ∇n(R)
(3×FFT)
3. ３次元ベクトルの計算: F (R) = ∇n(R) ∂f /∂|∇n(R)|
4. F の逆格子空間への変換: F (R) → F (G)
(3×FFT)
5. F の勾配の絶対値計算
∇ · F (G) = G · F (G) の計算
6. |∇ · F (G)| の実空間への戻し: G · F (G) → ∇F (R)
(1×FFT)
となり、合計 FFT は８回必要となる。
この GGA を平面波基底に移植する際、良く用いられている White-Bird のスキム [?] を
採用したが、それが今一つ移植に問題がありうまく動いていなかった。時々、エネルギー
が跳ぶのである。この移植版では GGA の密度勾配に依存する部分の交換相関ポテンシャ
(2)
ル Vex は
∑ ∂F (r0 ) ∂∇n(r0 )
(2)
Vex
(r) =
(3.1)
∂∇n(r0 ) ∂n(r)
0
r
33
と書かれている。この式（3.1）を評価する際、右辺の積の第一項は同じ位置での微分同士
であるので問題ないが、第二項は違う位置同士の変分を取るのでややこしくなる。そこで
それを逆講師空間で評価してまた実空間に戻すという方法がとられている。すなわち
n(r) =
∑
n(G) exp(iG · r)
(3.2)
G
を介して偏微分を取ると解して
∑
[
]
∂∇n(r0 )
=i
G exp iG · (r0 − r)
∂n(r)
(3.3)
G
となる。したがって式（3.1）右辺の積の第一項 K(r) のフーリエ変換を K(G) として、
(2)
Vex
(r) = i
∑
G · K(−G) exp(−iG · r)
(3.4)
G
(2)
つまり Vex は G · K(−G) の逆フーリエ変換で与えられると解釈できる。
一方、今回試みた Balbas らの方法は White-Bird と違い、密度の空間微分を取るとき、
フーリエ変換した逆格子空間での微分表式から求めるのではなく、直接実空間で差分より
求めるものである。演算としてはより直接的である。すなわち
(2)
Vex
(r) = ∇r ·
∂∇fex
∂n(r0 )
(3.5)
とするものである。ここで注意したいのは、式（3.5）右辺の積の第一項と第二項とでは実
空間座標が違うことである。したがってこれもやはりローカルなポテンシャルと違い、空
間の一点 r でのポテンシャルを求めるのに他の全ての点 r0 での fex の微分を知らなければ
ならない点は WB のものと本質的に同じである。
ということで、BS、WB ともに本質的に同じ計算負荷であるが、しかしながら特に大き
なサイズの計算ではちょっとした移植の違いが大きな差を生むことがある。そのため、両者
をいつも比較したほうが良い。切り替えは、ソースファイル elecV GGA.f90 の MODULE
elecData PBE のはじめの方にある、
! type
CHARACTER(LEN=2) :: implem_type = ’WB’!’WB’ = White & Bird
!’BS’ = Baldas & Soler
を WB あるいは BS と切り替えることでできる。
3.4.2
ポテンシャルの作成 atm
atm は Martins のサイトで入手できる。解凍すると atom-5.693 というディレクトリが
現れる。最近はもう少しバージョンアップされているだろうが多分ほとんどマイナーな変
更にすぎない。
34
コンパイルは古い atom に比べて格段に整理されている。そのディレクトリの直下に
メインプログラム atom.f のみがあるが、それ以外の全てのサブルーチンは souce/ディ
レクトリに収められている。これらのサブルーチンは全てが、コンパイルしてライブラ
リ libatom.a として統合されることを想定している。これらのソースプログラムの中で
Osaka2k に適応するため唯一ソースコード pseud2.f を書き換える必要があった。これは
Osaka2k のホームページにあるものと置き換える。
このライブラリ libatom.a の作り方は全てそのディレクトリにある Makefile で行われ
るのでこれをターゲットマシンに則して書き換えれば良い。
make により libatom.a が作られることになるが、それを親ディレクトリに移す。良く
分からないがマシンによってはそのままではライブラリがアップデートされないマシンも
ある。その時は
ranlib libatom.a
としてやる。
最後にメインディレクトリにおいてメインプログラム atm とコンパイル、リンクして終
わりである。
3.4.3
GGA の使い方
ポテンシャルの準備
古い atom と同じように、原子名を使って、pseudo01.dat を nm.pot に fort.13 を
nm.pwf とリネームしデーターディレクトリ ppot/の元でさらにサブディレクトリ gga/
を作りそこに格納する。
例えば Si の例で示すと、インプットファイル atom.dat に
pe
Silicon
PRB 52 13281 (1995)
Ecut ~ 16Ry
tm2
n=Si c=pb
0.0
0.0
0.0
0.0
0.0
3
2
3
0
2.00
0.00
3
1
2.00
0.00
2.13
2.57
1.00
0.00
0.00
1.00
l=0,1,2 as local
0.0
-2.00
としてやる。そうすることで atom の中では、まず全電子計算で
Si output data for orbitals
---------------------------nl
s
occ
eigenvalue
1s
2s
2p
3s
3p
0.0
0.0
0.0
0.0
0.0
2.0000
2.0000
6.0000
2.0000
2.0000
-130.91497657
-10.20480819
-7.02579802
-0.79146107
-0.30063401
(-130.36911238)
( -10.14892693)
( -7.02876267)
( -0.79662741)
( -0.30705179)
となる。これはカッコ内の gga を使わなかったときの値と比較してある。擬ポテンシャル
データは
35
test of eigenvalues
rcut =
eval =
l
s
p
s
2.22
-0.79147
vps(0)
1.430
-2.423
p
2.46
-0.30063
vpsmin
-5.173
-5.192
at r
1.508
1.266
Si output data for orbitals
---------------------------nl
s
occ
eigenvalue
1s
2p
0.0
0.0
2.0000
2.0000
-0.79146671
-0.30063280
(-0.79661280)
(-0.30704485)
を与える。やはり lda の時の値をカッコ内に示す。
inip
inip では特に gga になって変わるものはない。
pwm
pwm.para で gga の指定はオプションにより行う。
OPTION BEGIN
gga ON
OPTION END
とする。なおこれにより共役勾配過程で、オプション
use_numeric_deriv ON
が自動的に付加される。
36
第4章
ポテンシャルデータの再構築
これまでに、スピン、内殻補正をはじめとしてポテンシャルに関して拡張してきた。そ
の上、GGA も取り入れたりしたため、ポテンシャルに種類が途端に膨れ上がって、ポテ
ンシャルディレクトリに少し整理が必要となってきた。そこで
バージョンが
VersionNo
= ’1.51’
VersionDate = ’30 Sep 2005’
からは、従来のポテンシャルディレクトリに修正を加えた。
4.1
ポテンシャルディレクトリ
従来はポテンシャルの種類としては、inip.para で
potential type (spin, NLCC, relativistic)
0
0
0
と、スピン、内殻補正、相対論補正の３つを区別してきた。それに応じたポテンシャルデー
タは spn、nlc、rel と三通り用意してきたが、しかしそれらの組み合わせも出てくるし
（例えばスピンと内殻補正）、さらに GGA も取り入れるとその組み合わせの数は途端に膨
れ上がる。また、取り扱う原子が複数の時、どれが優先されるかの指針も明確にしないと
混乱に陥る。
そこで、ポテンシャルディレクトリに図 4.1 のような階層構造を導入した。それぞれの
ディレクトリには表 4.1 のような役割を持たせている。
dir
nom
nlc
spn
scc
rel
Table 4.1: Potential directories.
potential type
normal potential
nonlinear core correction
spin polarization
spin polarization with nonlinear core correction
relativistic effects
図 4.1 を参照にしながら、その基本指針を述べる。
37
ppot/
nom
nlc
spn
scc
rel
gga/
nom
nlc
spn
scc
rel
図 4.1: Construction of the directory of potential data
1. LDA or GGA
LDA か GGA かは排他的であるので、まずこれでどのポテンシャルディレクト
リに属するかの大枠が決る。ある原子を GGA で扱いたいなら、他の全ての原
子も GGA で扱わねばならない。H だけは例外である。
2. 相対論的扱い
相対論的扱いをするときは、サブディレクトリ/rel にそのポテンシャルデー
ターを格納し、それを参照する。ある原子を相対論的で扱いたいなら、他の全
ての原子も相対論的で扱わねばならない。
3. スピン分極
スピン分極を考慮するときは、サブディレクトリ/spn にそのポテンシャルデー
ターを格納し、それを参照する。ある原子のスピン分極を扱いたいなら、他の
全ての原子もスピン分極するものとして扱わねばならない。
4. 非線形内殼補正
非線形内殼補正を考慮するときは、サブディレクトリ/nlc にそのポテンシャル
データーを格納し、それを参照する。非線形内殼補正に関しては特定の原子に
関してのみ扱うことができ、他の全ての原子は非線形内殼補正なしで取り扱え
る。スピンと非線形内殼補正を同時に考慮するときは、サブディレクトリ/scc
にそのポテンシャルデーターを格納し、それを参照する。
4.2
GGA ポテンシャルの作成 atm
GGA を扱うには、昔の擬ポテンシャル作成プログラム atom ではだめで、最近の atm を
使わなければならない。Osaka2k では ver. 5.693 のものを使っている。atm は Martins の
38
サイトで入手できる。1 解凍すると atom-5.693 というディレクトリが現れる。最近はも
う少しバージョンアップされているだろうが多分ほとんどマイナーな変更にすぎない。
コンパイルは古い atom に比べて格段に整理されている。そのディレクトリの直下に
メインプログラム atom.f のみがあるが、それ以外の全てのサブルーチンは souce/ディ
レクトリに収められている。これらのサブルーチンは全てが、コンパイルしてライブラ
リ libatom.a として統合されることを想定している。これらのソースプログラムの中で
Osaka2k に適応するため唯一ソースコード pseud2.f を書き換える必要があった。これは
Osaka2k のホームページにあるものと置き換える。
このライブラリ libatom.a の作り方は全てそのディレクトリにある Makefile で行われ
るのでこれをターゲットマシンに則して書き換えれば良い。
make により libatom.a が作られることになるが、それを親ディレクトリに移す。良く
分からないがマシンによってはそのままではライブラリがアップデートされないマシンも
ある。その時は
ranlib libatom.a
としてやる。
最後にメインディレクトリにおいてメインプログラム atm とコンパイル、リンクして終
わりである。
4.2.1
GGA ポテンシャルの作り方
GGA ポテンシャルを作るには atm へ、その入力ファイルで電子相関ポテンシャルの型
として pe としてやれば良い。
Si の例では、インプットファイル atom.dat に
pe
Silicon
PRB 52 13281 (1995)
Ecut ~ 16Ry
tm2
n=Si c=pb
0.0
0.0
0.0
0.0
0.0
3
2
3
0
2.00
0.00
3
1
2.00
0.00
2.13
2.57
1.00
0.00
0.00
1.00
l=0,1,2 as local
0.0
-2.00
としてやる。そうすることで atom の中では、まず全電子計算で
Si output data for orbitals
---------------------------nl
s
occ
eigenvalue
1s
2s
2p
3s
3p
0.0
0.0
0.0
0.0
0.0
2.0000
2.0000
6.0000
2.0000
2.0000
-130.91497657
-10.20480819
-7.02579802
-0.79146107
-0.30063401
(-130.36911238)
( -10.14892693)
( -7.02876267)
( -0.79662741)
( -0.30705179)
となる。これはカッコ内の gga を使わなかったときの値と比較してある。擬ポテンシャル
データは
1
http://bohr.inesc.pt/ jlm/pseudo.html
39
test of eigenvalues
rcut =
eval =
l
s
p
s
2.22
-0.79147
vps(0)
1.430
-2.423
p
2.46
-0.30063
vpsmin
-5.173
-5.192
at r
1.508
1.266
Si output data for orbitals
---------------------------nl
s
occ
eigenvalue
1s
2p
0.0
0.0
2.0000
2.0000
-0.79146671
-0.30063280
(-0.79661280)
(-0.30704485)
を与える。やはり lda の時の値をカッコ内に示す。
4.2.2
inip
inip.para で gga の指定はオプションにより行う。
OPTION BEGIN
gga ON
OPTION END
とする。この情報が以降の pwm などのプログラムに伝達される。
4.2.3
pwm
gga が使われていることは、inip *.inp から伝わってくる。なおこれにより共役勾配過
程で、オプション
use_numeric_deriv ON
が自動的に付加される。
4.3
非線形内殻補正
非線形内殻補正の指定は、その準備プログラム inip の入力で
potential type (spin, NLCC, relativistic)
0
1
0
のようにする。NLCC が 1 のとき交換相関項までに内殻補正を施し、2 のときはハートリー
項まで含めて内殻補正を行う。
はじめは、非線形内殻補正を行うとき、構成原子のどれでもひとつそれが入ると、他の
原子も全て NLCC のポテンシャルを使わないといけない構成であった。しかし、それで
は不都合が出てきたので、その制約を取り払った。指示した原子だけが NLCC の扱いを
受け、それ以外は内殻補正をしないポテンシャルを使う。
そのためには、inip.para でこれまでのように内殻補正を指定した後、オプションで
40
atomPnlCC=
na
atomPnlCC=
li
と、内殻補正したい原子の種類だけ指示する。したがって、いくら NLCC を 1 とセットし
ても、個々の原子に対する指示がない限り内殻補正はなされない。
41
第5章
Osaka2k でのスピンの取り入れ方
現執筆時点での Osaka2k、プログラムバージョンは
VersionNo
= ’1.66’
VersionDate = ’03 Apr 2006’
であるが、この時点でのでのスピンの扱いは以下のようにしている。
5.1
スピンを入れた計算の流れ
原子状態の擬ポテンシャルは Kleinman-Bylander (KB) 型を前提として、常磁性状態の
ものを参照状態として用いる。そうしないと、スピンアップとダウン状態で異なる波動関
数が出来、KB 型ポテンシャルにそのままでは適用しにくい。以下 KB 型ポテンシャルの
使用を前提として Osaka2k にスピンを取り入れる手続きを述べる。
1. ポテンシャルの準備
常磁性状態のものを参照状態として atom を動かす。例えば鉄の例だと次のようにな
るだろう。
pg
Iron
Guess
Ecut ~ 55 Ry
tm2
n=Fe c=cas
0.0
0.0
0.0
0.0
5
3
4
0
1.00
1.00
4
1
0.00
0.00
3
2
3.00
3.00
2.20
2.40
2.20
l=0 as local
0.0
0.0
できたポテンシャル、波動関数データをディレクトリー ppot/spn にセーブする。
2. inip
入力ファイル inip.para でポテンシャルの選択を
potential type (spin, NLCC, relativistic)
1
0
0
とする。それだけである。
3. pwm
スピンを使っているという指定は inip の中で行われており、pwm では特にスピンに
関する指定はない。ただし、注意すべきことは占有数に関する事項である。初期の
占有数に関する選択がいくつかある。これは多分例で示すのが一番である。
42
Te
T_elec_start0fermi
T_elec_end0fermi
iend_Tel_const
ikeep_Tel_const
maxIter0wfn
iter
図 5.1: Temperature control for Te
5.2
電子占有数のオプション
電子占有数を可変にするにはいろいろなやり方を試みたが、多分そのうちフェルミブロー
ドニングがベストであろう。それでこの技法について述べる。フェルミ分布はもちろん
f (ε) =
1
e(ε−µ)/kTe
+1
(5.1)
で、これにより各準位の電子占有数を決めてやる。これは計算技巧として導入されている
ものであるので、現実の物質での電子温度とは関係ない。
その分布を特徴づける電子温度 Te は与えるのであるが、典型的には図 5.1 のように次第
に変化させる。それぞれのパラメータのデフォルト値は表 5.1 に示されている。変えたい
ときはオプション入力することになる。pwm.para のオプション入力は典型的には
OPTION BEGIN
nband_extra=
4
fermi_broadening ON
T_elec_start0fermi=
0.1
T_elec_end0fermi=
0.001
iend_Tel_const=
25
OPTION END
のようになるだろう。
この意味は、まずフェルミブロードニングが意味あるようにするため計算するバンドの
数を増やす。そうしないと、計算するバンドの数は必要最低限のものしか取らない。1 次に
フェルミブロードニングを実行することを宣言する。その後、フェルミブロードニング温
度を 0.1 から 0.001 Ry まで変化させる。maxIter0wfn で指定される SCF 計算の繰り返し
の中で、一回目に電子温度は T elec start0fermi にセットされ iend Tel const 回目ま
で T elec end0fermi に下げられ、それ以降は maxIter0wfn 回目までその値に保たれる。
デフォルトでは iend Tel const は maxIter0wfn の 0.75 の整数に固定される。また前
の計算の波動関数を読み込むときは T elec start0fermi のデフォルト値は 0.01 である。
1
現在のバージョンでは、フェルミブロードニングを宣言すると同時にデフォルトで nband extra の値は
２にセットされる。
43
Table 5.1: Options of Fermi broadening. Temperature is given in Ry.
Option
Default
T elec start0fermi
0.1
iread=1
0.01
T elec end0fermi
0.001
ikeep Tel const
1
iend Tel const
0.75*maxIter0wfn
5.3
実例
早速、実例で示す。
5.3.1
O2 結晶
これは、このシリーズの中ですでに何回か取り上げられてる。結晶構造は
TITLE alpha_O2
DIMENSION 3
CELL
5.37500
3.42500
5.06600
90.00000
SYMMETRY NUMBER 12 QUALIFIER UNIQUE_b,CELL_1
ATOMS
NAME
O1
X
-0.10086
Y
0.00000
Z
POT
0.15960 o
47.90000
CHARGE
6.0000
90.00000
TEMP
0.5000
OCCUP
1.0000
SCAT
である。
inip.para で
Cutoff k radius (AMAX) given by lattice index without 2Pi
13.1
way to give sampling points (0:given manually, 1:calc)
1
number of k-sampling points
2
potential type (spin, NLCC, relativistic)
1
0
0
とパラメータを定めてやる。その結果、k 点サンプリングと平面波展開は
==================== K-Space Setup ==============================
k sampling point set
Nkpts =
3
No NM
index in p
in c
A/gmin i/o star
WTK
1
YP
-1 1 -1/ 4
-2 0 -1/ 4
0.35582
1
2
0.25000
2
YP
-1 1 -3/ 4
-2 0 -3/ 4
0.55738
1
2
0.25000
3
GN
-1 -1 -1/ 4
0 -2 -1/ 4
0.60300
1
4
0.50000
sum of wtk =
1.00000
==================== PW_Expansion ==============================
Cutoff in the reciprocal space
am
=
13.10000 (rel. units)
kcut =
11.58770 (ab^-1) with 2Pi
Ecut =
134.27490 (Ry)
UNIT of K
0.88456 (a.u.)
Planewave expansion
44
10
0
Tel
Stot
Eel
Te
10
10
S_tot
10
0
10
20
0
30
10
20
ITER
ITER
30
図 5.2: Temperature control for Te
with NHDIM = 6138
No Name
KB/IC in p
Nstr
1 YP
-1
1 -1/ 4
2
2 YP
-1
1 -3/ 4
2
3 GN
-1 -1 -1/ 4
4
Sum over WTK 1.000000
Inv
2
2
3
Eq
1
1
1
WTK
0.250000
0.250000
0.500000
NPW
6131
6128
6138
のようになる。
pwm.para で maxIter0wfn を 30 にセットしている。この結晶は価電子を 12 個持ってい
るので、デフォルトでは６つのバンドだけが計算される。そこでオプション入力で電子占
有数に関して
OPTION BEGIN
allow_posE ON
fermi_broadening ON
OPTION END
とフェルミブロードニングを宣言してやる。そうするとデフォルトで、電子占有数に関し
ては以下のような設定を取る。
--- Occupancy parameters --Initialize the occupation numbers
Nel0prim
=
12
real
=
12.000000
NEDIM
=
8
nband(d/u) =
8
8
extra =
2 bands have been added
Occupancies will be varied during SCF calculation
mixing parameter of Occ =
1.0000
Fermi broadening options
T_elec_start0fermi =
0.100
until ikeep_Tel_const =
T_elec_end0fermi
=
0.001
at
iend_Tel_const =
current T_elec0fermi =
0.10
1
22
その結果は、
Etot
==============
-63.67564427
Eel
=============
-31.15096064
delta E
==========
-8.305E-10
resid
==========
6.412E+00
TOTAL ENERGY
--ELECTRONIC ENERGY
Sum of eigenvalues
-9.59508430
-Hartree term
-25.72105314
Exchange correlation
-13.37964603
-Sum of Vxce
17.54482284
-----------------------------------------Electronic Energy
-31.15096064
Electrostatic Energy
-32.52468363
Present Total Energy
-63.67564427
45
(RY./CELL)
iter
===
31
となり、スピン統計に関しては
Spin statistics
Total m =
-2.00000
Total S =
-1.00000
ave m^2 =
0.00009
max S =
0.41808
local moments
kind/atom
m
1 o
-1.55E+00
( Bohr)
= Nu - Nd
min S =
radius
1.4000
-37.78726
v_frac
9.83E-02
となっている。
SCF 計算の過程はいつものように pwm *.etot で見れるが、スピン特有の情報はまた
ファイル pwm *.spn に出力されている。Total m とは
∆N = N↓ − N↑
(5.2)
のことであり、Total S とはその 1/2 である。局在モーメントに関しては、同じ kind に
属するものの中で和を取っている。
SCF 過程での電子エネルギーや、全スピンの変化は図 5.2 に示される。電子温度 Tel が
0.01 あたりから有限モーメントが生じ始めていることが分かる。
5.3.2
Si 結晶中の Fe 不純物
Si 中の格子位置の Fe 不純物の計算をする。Si64 のスーパーセルを用いる。Fe はそのスー
パーセルのちょうど中心
57
0.5000000
0.5000000
0.5000000
2
fe
に入れる。荷電状態として Fe+ の状態を計算する。中性だとスピンが出ないためそうした。
k 点サンプリングは Γ 点のみで、平面波展開は
==================== K-Space Setup ==============================
k sampling point set
Nkpts =
1
No NM
index in p
in c
A/gmin i/o star
WTK
1
GM
0 0 0/ 2
0 0 0/ 2
0.00000
1
1
1.00000
sum of wtk =
1.00000
==================== PW_Expansion ==============================
Cutoff in the reciprocal space
am
=
23.10000 (rel. units)
kcut =
7.07143 (ab^-1) with 2Pi
Ecut =
50.00510 (Ry)
UNIT of K
0.30612 (a.u.)
Planewave expansion
with NHDIM = 51627
No Name
KB/IC in p
Nstr Inv
Eq
WTK
NPW
1 GM
0
0
0/ 2
1
1
1
1.000000
51627
としている。
電子占有数に関しては、
OPTION BEGIN
allow_posE ON
nel_add=
-1
nband_extra=
4
46
22
10
1
10
0
20
10
18
0
10
20
10
30
0
10
20
iter
30
iter
図 5.3: Temperature control for Te in Si63 Fe
fermi_broadening ON
T_elec_start0fermi=
0.1
T_elec_end0fermi=
0.0005
iend_Tel_const=
25
spintargetatom=
57
OPTION END
としている。
--- Occupancy parameters --Initialize the occupation numbers
Nel0prim
=
259
real
=
259.000000
NEDIM
=
132
nband(d/u) =
132
132
extra =
2 bands have been added
Occupancies will be varied during SCF calculation
mixing parameter of Occ =
1.0000
Fermi broadening options
T_elec_start0fermi =
0.100
until ikeep_Tel_const =
T_elec_end0fermi
=
0.001
at
iend_Tel_const =
current T_elec0fermi =
0.10
1
25
SCF 計算の経過は図 5.3 に示される。電子エネルギー Eel および、全スピンの変化が、
電子温度とともに示されている。
Etot
==============
-537.54540704
Eel
=============
18.96009012
delta E
==========
-1.349E-03
resid
==========
1.752E-05
iter
===
31
=== Thermodynamic potentials (Ry) ===
p V
:
-S T
:
T_elec0fermi
Entropy
Internal energy
Enthalpy
Helmholtz
0.0000E+00
0.0000E+00
5.000000E-04
0.000000E+00
Spin statistics
Total m =
Total S =
ave m^2 =
max S =
local moments
kind/atom
1 si
-1.00000
-0.50000
0.00132
49.80555
m
1.83E-02
( Bohr)
:
:
:
= Nu - Nd
min S = -1968.61356
radius
2.5494
v_frac
8.01E-03
47
-537.545407
-537.545407
-537.545407
.....
16 si
57 fe
1.11E-01
-2.65E+00
2.5494
2.0480
7.90E-03
4.16E-03
と全モーメントが１、Fe の局在モーメントが 2.6 となっている。
関連した技法
5.4
5.4.1
局在モーメントの半径
局在モーメントを評価するのに、原子半径 r なるものを決めてやる。r の範囲内で ∆ρ
を積分する。問題はその半径をいくつに取るかである。それは私は分からない。いろいろ
やってみるしかないが、そこでいくつかの選択を残した。
1. auto ソースコード ctrlVar pwm.f90 の中の
CHARACTER(LEN=4) :: choiceMT = ’cutp’
変数で設定している。cutp で擬ポテンシャルを作るときのカットオフ半径から（あ
る係数を掛けて）計算する。
2. given fixd で原子で決まった値を用いる。その値はソースコード atomUtils.f90 の
中のルーチン RadiusOfSpinMom で与えられる。
5.4.2
原子位置の最適化
5.4.3
部分 DOS
謝辞
スピンを考慮した擬ポテンシャルの使い方につてい、物質材料研・佐々木泰三氏にご教
示願った。
48
第6章
6.1
内殻補正の計算
atom で内殻補正のための入力
擬ポテンシャル作成のための入力ファイル atom.dat は、のようになっている。内殻補
Table 6.1: atom.dat のフォーマット
itype
ititle
ikerk
nameat icorr + ispp
znuc
zsh
ncore
nval
n
l
...
rcs
rcp
rsh
rmax
aa
zo↓
zo↑
evd
rcd
rcf
cf ac
bb
rcf ac
正を行うときはこのうち、最初と最後の行が関係してくる。
一行目 計算種（itype）として内殻補正をするものは
itype=
• pe: 交換項に関する内殻補正をした擬ポテンシャル生成
• ph: ハートリー項まで含めた内殻補正をした擬ポテンシャル生成
. . .
最終行 擬ポテンシャル生成のカットオフ半径（rc ）データの後、内殻補正パラメータ
（cfac と rcfac）。
まず内殻補正は交換相関項のみで行うか、それともハートリー項まで含めて行うかを決
めなければならない。それが第一行目で指定される。
その次に、内殻電荷密度を真の値で記述するか、それとも部分的な値で記述するかを決
める（全補正か部分補正）。それが最終行で指定される。そこでリストされている内殻補
正パラメータ、cf ac と rcf ac の値により、全補正か部分補正かが決まる。それらの意味を
図 6.1 を参照しながら説明する。
全補正 cf ac < 0 内殻電荷分布は真のものをそのまま使う。
49
rc
rv
cfac
rcc
図 6.1: Charge densities of the model core and valence electrons of silicon.
部分補正 cf ac > 0
このとき内殻補正を行うカットオフ半径 rcc の指定方法は二つある。
• 直接その値を rcfac で指定することで与える。当然 rcfac の値は正である。こ
のとき cfac の値は無視される。
• rcfac を負に与えると、その値は無視され、与えられた正の cfac に対し、
ρc (rcc ) = cfac × ρv (rcc )
(6.1)
で与えられるカットオフ半径 rcc を計算してきてその値を用いる。
これで内殻補正したときの擬ポテンシャルを作成できる。このポテンシャルデータをディ
レクトリ ppot/nlc に格納する。
6.2
pwm での実行
メインプログラム pwm で内殻補正をするためには、その準備プログラム inip の入力で
potential type (spin, NLCC, relativistic)
0
1
0
のようにする。NLCC が 1 のとき交換相関項までに内殻補正を施し、2 のときはハートリー
項まで含めて内殻補正を行う。
pwm はその出力として、実空間での内殻電荷密度分布をファイル pwm *.cchg に出力す
る。そのフォーマットは通常の価電子電荷密度ファイル pwm *.rho と同じであるので、容
易に図示化できる。
50
added Jun 2005
はじめは、非線形内殻補正を行うとき、構成原子のどれでもひとつそれが入る
と、他の原子も全て NLCC のポテンシャルを使わないといけない構成であっ
た。しかし、それでは不都合が出てきたので、その制約を取り払った。指示し
た原子だけが NLCC の扱いを受け、それ以外は内殻補正をしないポテンシャ
ルを使う。
そのためには、inip.para でこれまでのように内殻補正を指定した後、オプションで
atomPnlCC=
na
atomPnlCC=
li
と、内殻補正したい原子の種類だけ指示する。
6.3
非線形内殻補正のまとめ
非線形内殻補正をうまく移植できた。
ハートリー項も含めた補正は擬ポテンシャル法にとって本来不必要なものである。それ
にも関わらずそれを取り込むと数値精度上の問題が発生するので、これを行うべきでない。
51
第7章
7.1
相対論の計算
相対論的な一般式
ノルム保存型の擬ポテンシャル法は相対論的領域へも拡張できる [23, 24]。擬ポテンシャ
ルの理論の完全性を期すため、以下相対論への拡張について述べる。
球対称ポテンシャルをもつ孤立原子に対する相対論的ディラック方程式の解は、全角運
動量 J = l + s、その z 成分 M 、およびヘリシティを示す演算子 K
K = β(Σ · l + 1)
の固有値 κ である。κ と j の関係は

−(l + 1) = −(j + 1 ),
2
κ=
l = (j + 1 ),
2
(7.1)
for j = l +
for j = l −
1
2
1
2
となる。角座標およびスピン成分に関する固有値 |JM i は、
(
)
)
(
l + m + 1 1/2 m
l − m 1/2 m+1
J
ΦM =
Yl σ ↑ +
Yl
σ↓
2l + 1
2l + 1
1
1
J = l + ,M = m +
2
2
)1/2
)
(
(
l−m+1
l + m 1/2 m
0
ΦJM 0 =
Ylm−1 σ↑ −
Yl σ↓
2l + 1
2l + 1
1
1
J0 = l − , M0 = m −
2
2
(7.2)
(7.3a)
(7.3b)
とスピノル関数で表され、動径方向に関してはスカラー関数 G(r) で与えられる。対応す
るポテンシャルは VJκ (r) である。以下、Jκ > 0 を J+、Jκ < 0 を J− と記する。
このポテンシャルを結晶中に導入する。元のポテンシャルの表式から、擬ポテンシャル
の演算子は
∑
V̂ ps =
|JM iVJ (r)hJM |
(7.4)
JM
となる。
実際にはこの式（7.4）よりも、|lmi で展開した表式の方が使い易いのでそのような表
式に直して用いる。一つが |li で展開したもので、
∑ [
]
V̂ ps =
|li V̄l (r) + VlSO (r)L · S hl|
(7.5)
l
52
となる。m には無関係である。ここに、V̄l (r) はスピンアップとダウンについて平均した
イオンのポテンシャル
V̄l (r) =
]
1 [
lVl−1/2 (r) + (l + 1)Vl+1/2 (r)
2l + 1
(7.6)
で、VlSO (r) はスピン軌道相互作用のポテンシャル
VlSO (r) =
]
2 [
Vl+1/2 (r) − Vl−1/2 (r)
2l + 1
(7.7)
である。これが文献 [23, 24] で使われている表式である。
もう一つ実用上有用なのは |lmi で展開したものである。スピンに関して対角成分は、
ps
V̂ave
=
l
∑
∑ ∑
l
|lmsiV̄l (r)hlms|
(7.8a)
m=−l s
非対角成分は
ps
V̂diff
=
l−1
1∑ ∑
[(l + m + 1)(l − m)]1/2 {|lm + 1βiVlSO (r)hlmα|
2
l
m=−l
+|lmαiVlSO (r)hlm + 1β|}
(7.8b)
で与えられる。Kleinman-Bylander の表式を使うときはこちらの方（7.8b）が有用である。
式（7.8）において、m の和の範囲が違うことに注意を払うべきである。
7.2
7.2.1
移植上いくつかの注意点
スピノールの表現
スピノール関数 χ(r) の表現
(
χ(r) =
ψ↓ (r)
ψ↑ (r)
)
(7.10)
それぞれの成分 ψσ (r) は平面波で展開される。N 個の平面波で展開されるとすると、χk (r)
は大きなベクトル
χk (r) = {ψk↓ (G1 ), ψk↓ (G2 ), . . . , ψk↓ (GN ); ψk↑ (G1 ), ψk↑ (G2 ), . . . , ψk↑ (GN )} (7.11)
で表現できる。つまりベクトル {χk (r)} は、前半の {ψk↓ (G)} と後半の {ψk↑ (G)} の二つ
で構成されていることになる。このことは格段特別のことではないように見えるが、しか
し対称操作などを施すとき、その順番などに注意を払わなければならないことを意味する。
53
7.2.2
Kleinman-Bylander 型の使用
上の理論は半局所擬ポテンシャルの表式である。実際のところは、擬ポテンシャルに
Kleinman-Bylander 型を用いている。上の表式が理解されれば、Kleinman-Bylander 型へ
の拡張は一筋縄でゆく。つまり、まず、
|VJ ϕJ i =
1
(ϕJ , VJ ϕJ )1/2
|VJ (r)ϕJ (r))
(7.12)
を作っておく。式（7.8）の演算子を働かせるところに、式（7.12）の |Vj+ li と |Vj− li を
ps
別々に作用させることで計算する。たとえば、V̂ave の場合であれば
|lmsiV̄l (r)hlms| →
]
1 [
(l + 1)(Yml )∗ Yml |VJ+ ϕJ+ ihVJ+ ϕJ+ | + l(Yml )∗ Yml |VJ− ϕJ− ihVJ− ϕJ− | (7.13)
2l + 1
のように置き換える。|VJ ϕJ i の部分は実際には平面波展開で表しておく。
7.2.3
電子の交換相関エネルギー
スピン偏極の効果はなしとしている。すなわち交換相関エネルギー Exc = Exc (ρ(r)) の
評価においては、
1
ρ̄(r) = (ρ↓ (r) + ρ↑ (r))
(7.14)
2
の平均密度を用いている。
7.2.4
波動関数の初期値
SCF 計算において波動関数の初期値は特別に意味はなく、どのような値から出発しても
最終状態は同じになることが求められる。にもかかわらず、SCF 計算において波動関数の
初期値は特別の形を準備した。スピノール関数は式（7.11）のように表現されている。こ
の表現に従って述べると、エネルギーの低い順で {χn (r)} を並べると、
χ2n+1 (r) = {ψn (r); 0}
χ2n (r) = { 0; ψn (r)}
···
(7.15)
のように設定した。言うまでもなく、SO 項が小さいときはスピンアップとダウンが混ざ
らない状態となることから、そのような状態を出発点にしようということである。
縮退の問題 もっとも SO 項を無視した場合でさえ、SCF 計算の中ではスピンアップとダ
ウンが混ざらない状態を保つことは困難である。そのような場合、式（7.15）の波動関数
を取ったとき、非対角項は 0 であるが、しかし対角項が同じであるので、対角化してしま
うと、結局スピンアップとダウンが混ざってしまう。縮退しているものはその任意線形結
合が解として許されるからである。
54
7.3
SCF 計算
例として、Sn にたいして SCF 計算をしてみる。相対論的計算は inip.para で
potential type (spin, NLCC, relativistic)
0
0
1
とするだけでよい。
計算条件はもっとも簡単なものにしておく。
==================== K-Space Setup ==============================
k sampling point set
Nkpts =
1
No NM
index in p
in c
A/gmin i/o star
WTK
1
GM
0 0 0/ 2
0 0 0/ 2
0.00000
1
1
1.00000
sum of wtk =
1.00000
==================== PW_Expansion ==============================
Cutoff in the reciprocal space
am
=
3.10000 (rel. units)
kcut =
2.75029 (ab^-1) with 2Pi
Ecut =
7.56410 (Ry)
UNIT of K
0.88719 (a.u.)
Planewave expansion
with NHDIM =
169
No Name
KB/IC in p
Nstr Inv
Eq
WTK
NPW
1 GM
0
0
0/ 2
1
1
1
1.000000
169
Sum over WTK 1.000000
つまり Γ 点だけの計算である。
初期波動関数はバンドごとにアップ、ダウン成分を
Spinors check
band
---1
2
3
4
5
6
7
8
---sum
down comp
---------1.000000
0.000000
1.000000
0.000000
1.000000
0.000000
1.000000
0.000000
---------4.000000
up comp
---------0.000000
1.000000
0.000000
1.000000
0.000000
1.000000
0.000000
1.000000
---------4.000000
sum
---------1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
---------8.000000
としている。
最初にスピンー軌道相互作用を取り入れない場合の計算を示す。これはオプション入力
が要され、デフォルトでは後で見るスピンー軌道相互作用を取り入れた計算が為される。
iter
====
1
2
3
4
5
6
7
8
9
10
11
12
Eel
(Ry/cell)
===============
0.9819582092
0.8419681574
0.8412449160
0.8401237017
0.8387485150
0.8375286405
0.8367551752
0.8361983172
0.8357845588
0.8355121913
0.8353458506
0.8352337365
deE
(Ry/cell)
============
-7.4861E+00
-1.3999E-01
-7.2324E-04
-1.1212E-03
-1.3752E-03
-1.2199E-03
-7.7347E-04
-5.5686E-04
-4.1376E-04
-2.7237E-04
-1.6634E-04
-1.1211E-04
Xsi
(Ry^2/cell)
============
2.4040E-02
5.5389E-05
1.1878E-05
2.9851E-05
4.1311E-05
2.3070E-05
2.5243E-05
1.4918E-05
9.0398E-06
8.6872E-06
4.4417E-06
4.4835E-06
55
nst/bk
aglmax
========
5/ 0
5/ 0
5/ 2
5/ 5
5/ 5
5/ 3
5/ 3
5/ 3
5/ 4
5/ 3
5/ 4
5/ 5
=============
0.666898244
0.169684019
0.030770060
0.058134918
0.082621556
0.070609754
0.029334778
0.027810265
0.033560215
0.022005066
0.016335470
0.011508869
13
14
15
16
17
18
19
20
0.8351610619
0.8351004259
0.8350628892
0.8350282313
0.8350084430
0.8349913082
0.8349792517
0.8349688763
-7.2675E-05
-6.0636E-05
-3.7537E-05
-3.4658E-05
-1.9788E-05
-1.7135E-05
-1.2056E-05
-1.0375E-05
1.4632E-06
2.4633E-06
7.6206E-07
1.1051E-06
4.3879E-07
7.4112E-07
2.7097E-07
4.8856E-07
5/
5/
5/
5/
5/
5/
5/
5/
4
3
5
2
3
3
6
4
0.007905886
0.010808862
0.006438418
0.009322284
0.004396429
0.004689096
0.003364602
0.003359226
==================== Final report of electronic convergence ===============
Etot
==============
-13.21763628
Eel
=============
0.83496888
delta E
==========
-1.038E-05
resid
==========
4.886E-07
TOTAL ENERGY
--ELECTRONIC ENERGY
Sum of eigenvalues
1.00056391
-Hartree term
-1.47763430
Exchange correlation
-4.32737767
-Sum of Vxce
5.63941694
-----------------------------------------Electronic Energy
0.83496888
Electrostatic Energy
-14.05260515
Present Total Energy
-13.21763628
iter
===
21
(RY./CELL)
と収束は遅い。
波動関数の最終値は
Spinors check
band
down comp
------------kpt= 1
npw=
169
1
0.811842
2
0.188158
3
0.352536
4
0.647463
5
0.903605
6
0.096395
7
0.994144
8
0.005857
------------sum
4.000000
up comp
----------
sum
----------
0.188158
0.811842
0.647464
0.352537
0.096395
0.903605
0.005856
0.994143
---------4.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
---------8.000000
となっている。また KS 準位は
===== KS levels =====
-0.435907
-0.435907
0.312019
0.312020
0.312021
0.312021
0.312129
と下の s 軌道は二重に、上の p 軌道は六重に縮退していることがわかる。
これにスピンー軌道相互作用を取り入れた場合の計算をしてみる。
iter
====
1
2
3
4
5
6
7
8
9
10
11
Eel
(Ry/cell)
===============
0.9740442653
0.8242424078
0.8156940385
0.8046616954
0.7999029669
0.7989699562
0.7987664897
0.7987118423
0.7986920063
0.7986863112
0.7986842216
deE
(Ry/cell)
============
-7.4940E+00
-1.4980E-01
-8.5484E-03
-1.1032E-02
-4.7587E-03
-9.3301E-04
-2.0347E-04
-5.4647E-05
-1.9836E-05
-5.6950E-06
-2.0896E-06
Xsi
(Ry^2/cell)
============
2.4381E-02
3.5656E-04
1.7372E-04
3.1925E-04
2.7730E-05
6.1192E-06
2.5239E-06
5.1525E-07
3.2390E-07
9.7063E-08
3.9853E-08
56
nst/bk
aglmax
========
5/ 0
5/ 0
5/ 3
5/ 1
5/ 1
5/ 3
5/ 4
5/ 7
5/ 8
5/ 9
5/ 6
=============
0.671785927
0.175507448
0.153832383
0.417714126
0.045297567
0.033879033
0.013956457
0.006414234
0.006117360
0.002505639
0.001821033
0.312168
12
0.7986832066
-1.0150E-06
2.3673E-08
13
0.7986832195
1.2818E-08
1.5507E-08
CG process is stopped because increase in Eel 1.2818E-08
5/
5/
8
5
0.001392111
0.000640278
==================== Final report of electronic convergence ===============
Etot
==============
-13.25392193
Eel
=============
0.79868322
delta E
==========
1.282E-08
resid
==========
1.551E-08
TOTAL ENERGY
--ELECTRONIC ENERGY
Sum of eigenvalues
0.97973474
-Hartree term
-1.48285503
Exchange correlation
-4.29494025
-Sum of Vxce
5.59674376
-----------------------------------------Electronic Energy
0.79868322
Electrostatic Energy
-14.05260515
Present Total Energy
-13.25392193
Spinors check
band
down comp
------------kpt= 1
npw=
169
1
0.187825
2
0.812175
3
0.500019
4
0.499993
5
0.257136
6
0.742864
7
0.464676
8
0.535299
------------sum
3.999987
iter
===
13
(RY./CELL)
up comp
----------
sum
----------
0.812175
0.187825
0.499981
0.500007
0.742864
0.257136
0.535324
0.464701
---------4.000013
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
---------8.000000
KS 準位は
===== KS levels =====
-0.434053
-0.434053
0.290174
0.290186
0.312027
と、今度は p 軌道は二重と、四重縮退に分かれる。
57
0.312027
0.321713
0.321714
付 録A
擬ポテンシャルテスト
第 2.2 節で述べられたように、理想的な、transferable な擬ポテンシャルは生成された
ときと原子の電子配置によらないはずです。擬ポテンシャル生成プログラム atom で作ら
れたポテンシャルでは、異なる電子配置を使ってテストすることができます。
例えば、スピンなしの C 原子の場合ですと、励起状態 (2s1 )(2p3 ) は全電子計算をさせ
ると、
C output data for orbitals
---------------------------nl
s
occ
eigenvalue
kinetic energy
1s
0.0
2.0000
-19.95681799
31.47302470
2s
0.0
1.0000
-1.03415562
3.31872486
2p
0.0
3.0000
-0.42823707
2.50323341
--------------------------------------------total energy
=
-74.24387371
pot energy
-67.28576197
-11.05623425
-9.33465366
となります。
一方、擬ポテンシャルを通常のように基底状態 (2s2 )(2p2 ) で作り、その擬ポテンシャル
により励起状態 (2s1 )(2p3 ) を計算させて見ます。
そのためには入力ファイル atom.dat を
pg
Carbon
tm2
n=C c=ca
0.0
0.0
0.0
1
2
2
0
2.00
0.00
2
1
2.00
0.00
1.25
1.25
0.00
0.00
pt
Carbon
tm2
n=C c=ca
0.0
0.0
0.0
0
2
2
0
1.00
0.00
2
1
3.00
0.00
0.0
0.0
0.0
0.0
0.0
0.0
と２回分の入力を入れます。２回目の計算は純粋に擬ポテンシャルでの計算ですので、内
殻軌道数 ncore は 0 とします。またカットオフ半径は不要です。
すると atom は初めに基底状態配置を (2s2 )(2p2 ) を用いて擬ポテンシャルを生成し、そ
のデータをファイル pseudo.dat01 に落とします。その計算過程はページ 20 に記述されて
いる通りです。
58
次にその擬ポテンシャルのデータを用いて、(2s1 )(2p3 ) 状態を計算します。やはり擬ポ
テンシャルのセクションの結果は
nl
s
occ
eigenvalue
kinetic energy
1s
0.0
1.0000
-1.03309977
0.84527275
2p
0.0
3.0000
-0.42854824
2.41690260
--------------------------------------------total energy
=
-10.08726669
excitation energy
=
pot energy
-5.39786323
-6.22634564
0.60398389
これを最初の (2s1 )(2p3 ) に関する全電子計算と比較すると、各軌道エネルギーはほぼ 10−3
くらいの相対誤差で一致していることがわかります。これは非常に良い transferability を
示していると見ることができます。
また擬ポテンシャルの扱いの中で、(2s1 )(2p3 ) 状態の基底状態 (2s2 )(2p2 ) に対する励起
エネルギーは 0.60398389 (Ry) となりますが、これは全電子計算の場合（同じくページ 20）
と比較して
-74.24387371-(-74.84852443)=0.60465072
となり、相対誤差 10−3 で非常に良い一致を示します。
スピン分極エネルギー このことを利用して、原子状態でのスピン分極エネルギーを求め
ることができます。まず、分極なしの状態 (2s{ 12 ↓, 12 ↑})(2p{ 12 ↓, 12 ↑}) で擬ポテンシャル
を作り、それを用いてスピン分極状態 (2s{ 21 ↓, 12 ↑})(2p{ 22 ↓, 0 ↑}) を計算させます。
入力ファイル atom.dat に
pg
Carbon
tm2
n=C c=cas
0.0
0.0
0.0
1
2
2
0
1.00
1.00
2
1
1.00
1.00
1.25
1.25
0.00
0.00
pt
Carbon
tm2
n=C c=cas
0.0
0.0
0.0
0
2
2
0
1.00
1.00
2
1
2.00
0.00
0.0
0.0
0.0
0.0
0.0
0.0
と記述します。
プログラムは分極なしの状態 (2s{ 21 ↓, 21 ↑})(2p{ 12 ↓, 12 ↑}) で擬ポテンシャルを作り、次
にスピン分極状態 (2s{ 21 ↓, 12 ↑})(2p{ 22 ↓, 0 ↑}) を計算します。二段階目の最後に
excitation energy
=
-0.09206819
と両者のエネルギー差を打ち出します。つまりスピン分極により、−0.0920 [Ry] のエネル
ギー利得があったことがわかります。
59
付 録B
その他
まだまだ続く。
60
References
[1] J. C. Phillips and L. Kleinman, Phys. Rev. 116, 287 (1959)
[2] M. H. Cohen and V. Heine, Phys. Rev. 122, 1821 (1961)
[3] B. J. Austin, V. Heine, and L. J. Sham, Phys. Rev. 127, 276 (1962)
[4] V. Heine, Solid Stete Physics 24, (1970) p. 1, and all the volume.
[5] W. C. Topp and J. J. Hopfield, Phys. Rev. B 7, 1295 (1973)
[6] D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett. 43, 1494 (1979)
[7] G. B. Bachelet, D. R. Hamann, and M. Schlüter, Phys. Rev. B 26, 4199 (1982)
[8] E. L. Shirley, D. C. Allan, R. M. Martin, and J. D. Joannopoulos, Phys. Rev. B
40, 3652 (1989)
[9] N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991)
[10] T. L. Gilbert, Phys. Rev. B 12, 2111 (1975)
[11] G. A. Henderson, Phys. Rev. A 23, 19 (1981)
[12] M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979)
[13] L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982)
[14] D. M. Bylander and L. Kleinman, Phys. Rev. B 41, 907 (1990)
[15] P. E. Blöchl, Phys. Rev. B 41, 5414 (1990)
[16] X. Gonze, P. Käckell, and M. Scheffler, Phys. Rev. B 41, 12264 (1990)
[17] X. Gonze, R. Stumpf, and M. Scheffler, Phys. Rev. B 44, 8503 (1991)
[18] U. von Barth and L. Hedin, J. Phys. C: Solid State Phys. 5, 1629 (1972)
[19] D. M. Ceperley and B. Alder, Phys. Rev. Lett. 45, 566 (1980)
[20] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981)
[21] A. Zunger, Phys. Rev. B 22, 649 (1980)
[22] S. C. Watson and E. A. Carter, Phys. Rev. B 58, R13309 (1998)
61
[23] L. Kleinman, Phys. Rev. B 21, 2630 (1980)
[24] G. B. Bachelet and M. Schlüter, Phys. Rev. B 25, 2103 (1982)
[25] S. G. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982)
[26] D. Porezag, M. R. Pederson, and A. Y. Liu, Phys. Rev. B 60, 14132 (1999)
[27] D. Vanderbilt, Phys. Rev. B 41, 7892 (1990)
[28] C. J. Pickard, B. Winkler, R. K. Lee, M. C. Payne, M. H. Lee, J. S. Lin, J. A. White,
V. Milman, and D. Vanderbilt, Phys. Rev. Lett. 85, 5122 (2000)
62