KKR-Green 関数法によるバンド計算
佐藤和則（１９９８年６月）
小倉昌子（２００５年２月改訂）

目次
1

はじめに

3

2

プログラムの入手、コンパイルと実行
2.1 プログラムの入手 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 コンパイル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 実行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5

3

KKR-Green 関数法
3.1 Single-Site の問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Schrödinger 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 漸近形 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Wronskian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.4 自由空間の Green 関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Traditional KKR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Muffin-tin potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Cell-centered expantion . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 構造定数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.4 KKR-matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Crystal Green’s Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 積分方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Dyson-equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 周期境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Density of States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.5 Contour integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 散乱問題の形式解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Lippmann-Schwinger equation . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 散乱振幅 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Born Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.4 Multiple Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
6
6
7
7
8
9
9
9
10
11
11
11
12
12
13
14
14
14
16
16
17

4

密度汎関数法4.1 Hohenberg-Kohn の定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 多電子問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 基礎関数としての密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18
18
18
19

1

4.2
4.3
4.4

4.1.3 変分原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Kohn-Sham 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
局所密度近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
逐次近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20
20
22
22

5

プログラムの構造

23

6

計算の実行
6.1 エネルギーの原点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 インプットファイル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 配列の大きさを指定するパラメータ . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 アウトプット . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25
25
25
29
29

7

計算例
7.1 Density of States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 Beryllium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.2 Iron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.3 Nickel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 スーパーセル法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 格子定数、格子緩和 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33
33
34
34
35
35
35

8

参考文献

36

図目次
1
2
3
4
5
6
7
8
9
10

Muffin-tin potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Contour integration of the Green’s Function . . . . . . . . . . . . . . . . . . . .
Born Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
T-matrix of Multiple Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . .
バンド計算の流れ図 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
エネルギーの原点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Be の状態密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Fe の状態密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Ni の状態密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Super Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9
15
17
18
23
25
33
34
35
36

表目次
1
2
3
4
5

原子単位 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
主なサブルーチンの概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
インプットファイル入力パラメータ（１） . . . . . . . . . . . . . . . . . . . . .
インプットファイル入力パラメータ（２） . . . . . . . . . . . . . . . . . . . . .
配列に関するパラメータ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

6
24
27
28
29

1

はじめに

物質の多彩な性質を理解するために、固体の電子状態について多くの理論的な研究がなされている。固体中での電子状態を、周期的に並んだポテンシャル中の一電子状態ととらえるバンド理論は多くの物性を説明することができたが、用いる有効ポテンシャルを何らかの実験結果とつき合わせて決めるなど、経験的な要素が少なからず含まれており、未知の物理量の予測ということになると不安な点が多い。また、電子間には静電相互作用が働くわけであるからセルフコンシス
テントな計算をしていないのは明らかに問題である。非経験的に第一原理から電子状態をセルフ
コンシステントに計算することは原子や少数の原子からなる分子について Hartree-Fock 法を用いて行なわれてきているが、これを固体にそのまま応用すると、交換相互作用の計算が煩雑過ぎることがわかる。そこで、J. C. Slater らは、交換相互作用の計算にパラメーター α を導入した
Xα 法を開発した。これにより固体についても第一原理からセルフコンシステントな計算ができるようになったが、結晶格子と原子番号のみを入力パラメーターとし、調節可能なパラメーターを含まない計算が可能になったのは、密度汎関数法が確立されてからである。Kohn-Sham の理論に局所密度近似を導入することで、望みの精度でセルフコンシステントな解を得ることができるようになった。この冊子では、赤井らによって開発された密度汎関数法に基づくグリーン関数
法 (KKR 法) を用いたバンド計算コードの使用法を解説する。バンド計算の手法としてグリーン関数法は非常に高速な方法であるということと、系のグリーン関数を求める方法であるので不純物問題や不規則合金の電子状態の計算にも自然に応用できるという特色がある。また、スーパー
セルで不純物をシミュレートしてやれば、格子緩和を考慮した電子状態を計算できるが、その計算の高速性から緩和を含めた不純物の電子状態の系統的な議論を進めることができる。

2

プログラムの入手、コンパイルと実行

Fortran で書かれているソースコードを手に入れ、コンパイルして実行ファイルを作るまでの手順を基本的な UNIX のコマンドを紹介しながら説明する。

2.1

プログラムの入手

KKR プログラムは WEB ページ：http://sham.phys.sci.osaka-u.ac.jp/∼kkr/ 上で公開されている。“Source codes and samples” と書かれたところをクリックし、パスワードを入力する。
（パ
スワードについてはサイトの管理者に問い合わせてください）プログラムの全体は、アーカイブ、
圧縮して cmd3.tar.gz となっているのでこれをダウンロードする。
次に、解凍とアーカイブしたファイルの取り出しを行う。ファイルの管理のことを考えてホームディレクトリの下に kkr という名前のディレクトリをつくってそこで作業を行なうことにする。
~> mkdir kkr~> cd kkr~/kkr> mv ../cmd3.tar.gz .
~/kkr> gzip -d cmd3.tar.gz~/kkr> tar xvf cmd3.tar.. はカレントディレクトリの真上のディレクトリを指す。mkdir は新しいディレクトリを作るコ
マンド、cd はカレントディレクトリを変更するコマンド、mv はファイルを移動するコマンドで3

ある。gzip -d はコマンド gzip で圧縮されたファイル (******.gz) を解凍するコマンドである。また、tar コマンドはアーカイブしたファイル (******.tar) からもとのファイルを取り出すコマンドである。解凍により cmd3.tar.gz は cmd3.tar に置き換えられる。tar コマンドにより cmd3.tar
内のファイルが取り出される。カレントディレクトリにあるファイルの一覧は、ls コマンドで得られる。ls と入力すると、~/kkr>ls
cmd3.tar
cmdtext.txt
cmdtext.txt.bak
cpa2002v006d

kkrnote.tex
memo.txt
readme.txt
sample_FeCoNi.tar.gz

sample_InMnAs.tar.gz
samples.tar.gzのようなファイルが、ディレクトリ kkr の下にできていることが確認できる。ディレクトリ cpa2002v006dの下にプログラムが置かれている。kkrnote.tex は KKR 法について memo.txt はプログラムの使用法について書かれたものなので参考にしていただきたい。readme.txt はプログラムの更新記録である。sample***.tar.gz には計算例が収められている。cmd3.tar はもう必要ないので消去しておく。

~/kkr> rm cmd3.tar
ディレクトリ cpa2002v006d に移動して ls と入力すると、~/kkr>cd cpa2002v006d~/kkr/cpa2002v006d>ls
data
gpd
f77_clock
gpdos.finintel_clock

makefile
readmesourceutilのようなファイルができていることが確認できる。makefile は、次の節で説明する make コマン
ドのためのファイルである。ディレクトリ source の下にプログラムのソースファイル******.f がおかれる。in, data のディレクトリはそれぞれ、インプットファイルを置く、計算して得られたポテンシャルを保存するために作った。
（ここで書かれていないその他のファイルは必ずしも必要なものではありません。これらの詳細について興味のある方は制作者に問い合わせてください）

2.2

コンパイル

KKR プログラムは、source ディレクトリ以下に置かれた 100 以上のサブルーチンからなっている。対象とする結晶や計算の規模に応じて一部分を変更してコンパイルしなおすことが多いので、コンパイルには make コマンドを使う。makefile にソースファイル名と実行ファイル名を書いておくと make コマンドはそれらのファイルをコンパイル、リンクして指定した名前の実行ファ
イルをつくる。二度目以降のコンパイル時には、make を実行すると変更したファイルのみコンパイルされ、実行ファイルがアップデートされる。
^/kkr/cpa2002v006d> make
これで、実行ファイル specx が得られる。プログラムを追加した場合は makefile にそのファイル
名を加えなければならない。

4

また、サブルーチン uclock、udate は計算機に依存する。これらはソースファイル中のコメントに従って計算機に適合するように変更する必要がある。いずれも時刻を取得するルーチンで
ある。

2.3

実行

結晶格子や格子点に置く原子の種類などを書いたインプットファイルを入力することでバンド計算が開始される。インプットファイルの内容と書き方については他の章で説明する。計算結果
は標準出力に書かれるのでファイルとして残しておくためには次のようにする。

~/kkr/cpa2002v006d> specx < in/input.in >output.out &不等号 hi で標準入力、標準出力をそれぞれ、in/input.in, output.out に切替えている。最後の&はコマンドをバックグラウンドで実行させることを意味する。計算には時間がかかる場合が多いのと、進行状況を確かめられるのでバックグラウンドで実行するのが便利である。フォアグラ
ウンドで実行した場合は計算が終るまで入力待ちのプロンプトは出てこないので他の作業はできない。フォアグラウンドのジョブをバックグラウンドにまわすには次のように Cntr+Z でジョブを中断し bg と入力すれば、~/kkr/cpa2002v006d> specx <in/input.in >output.out^Z
Suspended~/kkr/cpa2002v006d> bg[job. No] specx <in/input.in >output.out &~/kkr/cpa2002v006d>となって、入力待ちの状態になる。また、現在進行中のジョブを中止したい時には、kill コマン
ドをもちいる。

~/kkr/cpa2002v006d> kill %job No.
ジョブについての情報は、jobs, ps, top などのコマンドで得ることができる。

5

3

KKR-Green 関数法

3.1

Single-Site の問題

原点に置いた、ただ一つのポテンシャルによる散乱を考える。ただし、このポテンシャルは遠方で十分速く 0 に近づくと仮定されている。

3.1.1

Schrödinger 方程式

定常的な系の状態を記述する Schrödinger 方程式
"

#

h̄2 2
∇ + U (r) ψ(r) = Eψ(r)Hψ(r) = −2me(1)

を、散乱に対応する境界条件の下で解く。電子状態の計算に都合の良い単位は原子単位である。
この単位で、h̄2 = 1, e2 = 2, me = 1/2 である。

量
長さ
エネルギー

表 1: 原子単位単位h̄2= 1(bohr) = 0.529(Å)
me e2
me e4
= 1(Ry) = 13.6(eV)2h̄2

この単位で Schrödinger 方程式は次のようになる。

[−∇2 + V (r)]ψ = Eψ.

(2)

さらに、球対称なポテンシャルによる散乱を考えるならば、球座標を用いて次のように変数を分離することができる。

ψ(r) = Rl (r)Ylm (θ, φ),
·

(3)

¸

1 d
d
l(l + 1)
− 2 r2 ++ V (r) − E Rl = 0,r dr dr
r2

"

(4)

#

1
∂
l(l + 1)∂
1
∂2+
Ylm (θ, φ) = 0.
sin
θ
+
2
2
2
2
r sin θ ∂θ∂θ r sin θ ∂φr2

(5)

角度部分の解は球面調和関数である。また、Rl (r) = ul (r)/r とおくと、"#

l(l + 1)d2
+ V (r) − E ul (r) = 0
− 2+
dr
r2

(6)

と書くこともできる。
動径波動関数は次のように規格化する。
Z RM T0[rPl (r; E)]2 dr = 1
6(7)

3.1.2

漸近形

ポテンシャルの外では動径波動関数の厳密解は次のように書ける。
√√
Pl (r; E) = Al (E)jl ( Er) + Bl (E)nl ( Er)(8)jl , nl は球ベッセル、球ノイマン関数でそれぞれ、原点で正則な解と原点で非正則な解を表す。よっPて、ポテンシャルの外で全体の波動関数は、ψ(r) = L CL (E)Pl (r; E)YL (θ, φ) とかける。ここ
で、L = (l, m) である。散乱を表す境界条件は、遠方で波動関数が次の漸近形を持つことである。
√

ψ(r) → ei Ez +f (θ) i√Er
e
r

(9)

このことから、Al (E) と Bl (E) はパラメーター ηl (E) 用いて表すことができるとわかる。
√√
Pl (r; E) = eiηl [cos ηl jl ( Er) + sin ηl nl ( Er)]√lπ
→ eiηl sin( Er −+ ηl )(10)2
このように、ηl はポテンシャルによる遠方での波動関数の位相のずれを表す。また、次のように定義される第一種ハンケル関数を用いると、(1)hl= jl + inl√√(1)
Pl (r; E) = jl ( Er) + ieηl sin ηl hl ( Er)
√
√
√
(1)
= jl ( Er) − i Etl (E)hl ( Er)

(11)
(12)
(13)

のようにかけて、散乱問題は次のように定義される t マトリックスを求める問題に帰着する。

1tl (E) = − √ eiηl sin ηl
E
3.1.3

(14)

WronskianSchrödinger 方程式の帰結である Wronskian の関係式から解を構成することを考える。
"#

d2
l(l + 1)
+ V (r) − E rPl (r; E) = 0
− 2+
dr
r2
# Ã √
"
!
jl ( Er)l(l + 1)
d2
√
−E r
− 2+
= 0
dr
r2
nl ( Er)

(15)
(16)

が満たされていることから次の積分を部分積分を繰り返して評価できる。
Z r
0
Z ∞
r

√
√
rjl ( Er)V (r)rPl (r; E)dr = W [rjl ( Er), rPl (r, E)]r

(17)

√
√
1
rnl ( Er)V (r)rPl (r; E)dr = − √ eiηl cos ηl − W [rnl ( Er), rPl (r, E)]r
E

(18)

二つ目の等式中の定数は Pl (r; E) の漸近形を考慮して得られた。また、Wronskian の定義は次のようである。

W [a(r), b(r)] = a(r)db(r) da(r)
−
b(r)
dr
dr
7

(19)

√
√
一番目の式に nl ( E) を、二番目の式に jl ( E) をかけて足し算すると、√ Z ∞ √√√
jl ( Er< )nl ( Er> )V (r0 )Pl (r0 ; E)r02 dr0
Pl (r; E) = eiηl cos ηl jl ( Er) − E
0

Z ∞ √
√
√
√
iηl
→ e cos ηl jl ( Er) − Enl ( Er)jl ( Er0 )V (r0 )Pl (r0 ; E)r02 dr00√√= eiηl [cos ηl jl ( Er) − sin ηl nl ( Er)]√√= Cl (E)jl ( Er) − Sl (E)nl ( Er)√√(1)= [Cl (E) − iSl (E)] jl ( Er) + iSl (E)hl ( Er)

(20)

ここで、

Cl (E) = eiηl cos ηl√√=
EW [rPl (r; E), rnl ( Er)]r>RM T

(21)

iηl

Sl (E) = esin ηl√√= − EW [rPl (r; E), rjl ( Er)]r>RM T

(22)

である。
つまり、t マトリックスは Wronskian から次のように計算される。

1
Sl (E)
tl (E) = − √E Cl (E) − iSl (E)
3.1.4

(23)

自由空間の Green 関数

自由空間の Green 関数は次の微分方程式を満たす。

[∇2 + E]g(r, r 0 ) = δ(r − r 0 )(24)

この解は Fourier 変換によって求めることができる。

g(r, r 0 ) = g(r − r 0 )√0e(i E|r −r |)= −4π|r − r 0 |さらに部分波分解は次のように与えられることがわかっている。
√ X √
√
(1)
jl ( Er< )YL (r< )hl ( Er> )YL (r> )g(r, r 0 ) = −i E

(25)

(26)

L

ところで、非斉次方程式[∇2 + E]φL (r) = V (r)φL (r)(27)の一般解はこの方程式の特解と、対応する斉次方程式[∇2 + E]φ0L (r) = 0(28)の一般解の和で表される。Green 関数の定義から特解を書き下せるので、一般解は規格化定数を
ZL (E) と書いておくと、Z

φL (r) = ZL (E)φ0L (r) +g(r − r 0 )V (r0 )φL (r 0 )dr 0

とかける。

8

(29)

3.2

Traditional KKR

Korringa, Kohn, Rostoker によって導かれたバンド計算法 (KKR 法) をまとめる。
3.2.1

Muffin-tin potential

いままで、ただ一つのポテンシャルによる散乱を考えてきた。ここでは、固体やクラスターなどのように散乱中心が複数あり、その中で多重散乱をうけている電子を考える。このような問題
は、ポテンシャルの形状について次のマフィンティン近似を導入することで厳密に解くことができる。ポテンシャルを各原子球の中心を原点として、
Ãv(r) for r < RM T
0 for r > RM T

v(r) =!

(30)と近似する。RM T は重ならないようにとった原子球の半径である。
rm
m
Rm
O

n
Rn

rn

R

MT

v(r)

r

図 1: Muffin-tin potential

3.2.2

Cell-centered expantion

全系の波動関数 ψ(r) は次の Schrödinger 方程式をみたす。

[−∇2 + V (r)]ψ(r) = Eψ(r)

9

(31)

ψ(r) は自由空間の Green 関数を用いて次のようにかける。
Zg(r, r 0 )V (r 0 )ψ(r 0 )dv 0

ψ(r) =Z

0 )g(Rm + rm , Rn + rn

=

XZ

=

n

X
n

0 )ψ(R + r 0 )dv 0v(Rn + rnnnn0 )v(r 0 )ψ(r 0 )dv 0g(rm , Rn − Rm + rn
n
n
n

(32)

最後の変形では m 番目のセルに原点をもってきた。m 番目のセル内の波動関数を φm
L (r) とすると、
m[−∇2 + v(r)]φm
L (r) = EφL (r)

(33)

Z

m
0
φm
L (r) = ZL (E)φL (r) +00g(r, r 0 )v(r 0 )φm
L (r )dv(34)であった。そこで全系の波動関数を m 番目のセル内で、φm
L (r) で展開する。

ψ(r) =X

CLm φmL (rm )

(35)

L以上の式を代入して、
XZ

CLm [ZLm (E)φ0L (rm ) +0 )v(r 0 )φm (r 0 )dv 0 ]g(rm , rm
m
m L
m

L

=

XZ
n

3.2.3

0 )v(r 0 )g(rm , Rn − Rm + rnnX0 )dv 0CLn φnL (rn
n

(36)

L

構造定数

0 ) を部分波に分解してみる。
g(rm , Rn − Rm + rn

√ X √0 ) = −i Eg(rm , Rn − Rm + rn
jl ( Erm )YL (rm )
L

√

(1)
×hl (

0 + R − R |)Y (r 0 + R − R )
E|rn
n
m
n
m
L n
√ X √
= −i Ejl ( Erm )YL (rm )
L

×4π

X

√
0
00
(1)
il−l +l CLL0 L00 hl00 ( E|Rn − Rm |)L0 ,L00√
0 )×YL00 (Rn − Rm )jl0 ( Ern0 )YL0 (rn
X

mn 0
0
φ0L (rm )gLL
0 φL0 (rn )(37)√√ X l−l0 +l00
(1)i
CLL0 L00 hl00 ( E|Rn − Rm |)YL00 (Rn − Rm )= −i4π E

(38)

=

LL0

mn
gLL
0

L00
mn を構造定数と呼ぶ。これは結晶格子による量であり、格子点にあるポテ
ここにあらわれた gLL
0
ンシャルによらない。

10

3.2.4

KKR-matrix

Green 関数の部分波分解ができたので波動関数の展開係数が満たすべき方程式を導くことがで
きる。
X

Z

CLm [ZLm (E)φ0L (rm )0 )v(r 0 )φm (r 0 )dv 0 ]g(rm , rm
m
m
m L

+

L

=

X Z X
n6=m

Z

CLm ZLm (E)φ0L (rm ) =L0 )v(r 0 )φn (r 0 )dv 0 C n
φ0L0 (rn
n L0
n L n

(40)

L

X
X

X X

φ0L (rm )

Z

mn
gLL
0

n6=m L0

L

=

(39)

LL0
0 )v(r 0 )φm (r 0 )dv 0g(rm , rm
m
m L
m

+
X

X

0 )dv 0CLn φnL (rn
n

mn 0
0
0
φ0L (rm )gLL
0 φL0 (rn )v(rn )X X

φ0L (rm )

L

n6=m

mn n
n
gLL
0 YL0 (E)CL0

L0

結局次のような永年方程式を得る。
X

mn n
n
[δLL0 δmn ZLm (E) − gLL
0 YL0 (E)] CL0 = 0.

(41)

L0 n

ここで出てきた、ZLm (E), YLm0 (E) は前節の考察から、

ZLm (E) = Clm (E) − iSlm (E)(42)

YLm0 (E) = Slm (E)(43)

であることがわかる。つまり、自明でない解が存在するためには、
mn ndet [δLL0 δmn − gLL
0 tL0 (E)] = 0(44)が成り立たなくてはならない。

3.3
3.3.1

Crystal Green’s Function
積分方程式

前章では自由空間の Green 関数を使って、部分波の方法により多重散乱をうけた電子の波動関数を構成した。ここでは、系の Green 関数 G(r, r 0 ) の構成を試みる。
まず、G が満たす積分方程式をつくる。G は次の Schrödinger 方程式を満たす。

[∇2 + E − V (r)]G(r, r 0 ) = δ(r − r 0 )(45)

次のように考えれば、[∇2 + E]G(r, r 0 ) = V (r)G(r, r 0 ) + δ(r − r 0 )(46)

f (r, r 0 )=G は次のように書けることがわかる。
Z

G(r, r 0 ) =g(r, r 00 )f (r 00 , r 0 )dv 00
Z=g(r, r 0 ) +g(r, r 00 )V (r00 )G(r 00 , r 0 )dv 00
11

(47)

3.3.2Dyson-equation

先に見たように自由空間の Green 関数は次のように展開された。

g(r + Rm , r 0 + Rn ) = δmn g(r, r 0 ) +X

mn
0JL (r)gLL
0 JL0 (r )

(48)

LL0

√ X
g(r, r 0 ) = −i E
JL (r< )HL (r> )

(49)

L

ここで、

√
JL (r) = jl ( Er)YL (r)√(1)
HL (r) = hl ( Er)YL (r)

(50)
(51)

である。
同じように、G を部分波に展開してその展開係数を求めることを考える。G(r, r 0 ) も Schrödinger
方程式の解であるから、Pl (r)YL (r) で展開できる。g の展開との類推から、
0
G(r + Rm , r 0 + Rn ) = δmn Gms (r, r ) +X

0
JL (r)Gmn
LL0 JL0 (r )

(52)

LL0

√ X
0
JL (r< )HL (r> )Gms (r, r ) = −i E

(53)

L

√
JL (r) = Pl ( Er)YL (r)√HL (r) = Ql ( Er)YL (r)

(54)
(55)

mn
と展開できるとして、係数 Gmn
LL0 を決定する。直接、積分方程式に代入することで GLL0 にたいする代数方程式を得る。
mn
Gmn
LL0 = gLL0 +

X

ml l
ln
gLL
00 tL00 GL00 L0

(56)

L00 l

3.3.3

周期境界条件mn
結晶のように散乱中心が周期的にならんでいる場合 Gmn
LL0 , gLL0 の Fourier 変換は次のように与えることができる。

Z

Gmn
LL0

=

mn
gLL
0

=

1
GLL0 (k)e−ik(Rm −Rn ) dk
τ 1stBZ
Z
1
gLL0 (k)e−ik(Rm −Rn ) dk
τ 1stBZ(57)(58)

これを先に得た方程式に代入して、X

gLL00 (k)tL00 GL00 L0 (k)(59)[δL00 L − gL00 L (k)tL ] gL00 L0 (k)(60)

GLL0 (k) = gLL0 (k) +L00を得る。これを解くと、GLL0 (k) =X
L00が得られる。

12

3.3.4

Density of States

セルフコンシステントな計算を進めるため、また実験と比較する物理量を計算するために最も重要なのは電子の密度分布である。密度分布が系の Green 関数から直接得られることを示すために Green 関数の固有関数による展開を考える。Schrödinger 方程式 Hφ = Eφ の固有値 En に属する固有関数を φn と書く。

Hφn = En φn(61)[E − H]G(r, r 0 ) = δ(r − r 0 )(62)Green 関数の満たす方程式は、

G を次のように展開して、E − H をかけると、X

G(r, r 0 ) =Gn (r 0 )φn (r)

(63)

n

X

(E − H)G(r, r 0 ) =Gn (r 0 )(E − H)φn (r)nδ(r − r 0 )X=

Gn (r 0 )(E − En )φn (r)

n

よって、

Gn (r 0 ) = φ∗ (r 0 )/(E − En )

(64)

以上より、
X φ∗ (r 0 )φn (r)nG(r, r 0 ) =n

E − En(65)を得る。E = En で発散してしまうが、E をわずかに複素数にすることで意味を持たせることが
できる。
恒等式、

1
1
= P. − iπδ(x)x + i²
x

(66)

を使えば、

G(r, r 0 ) =X φ∗ (r 0 )φn (r)
n

n

= P.

E + i² − En
X φ∗ (r 0 )φn (r)nn

E − En− iπX

δ(E − En )φ∗n (r 0 )φn (r)

(67)

n

とかけて、密度分布と Green 関数のあいだに次の関係が成り立つことがわかる。

ρ(r, E) =X

δ(E − En )φ∗n (r)φn (r)n1= − ImG(r, r)
π
13

(68)

3.3.5

Contour integration

Green 関数のエネルギーに関する積分は、状態密度の複雑な構造のため大変難しい。実軸上で定義された Green 関数を複素エネルギー平面に接続し、実軸から離れて眺めると、状態密度の構造はならされて精度の良いエネルギー積分が実行できることを示す。

G(r, r; z) =X φ∗ (r)φn (r)
n

z−E

n

=

Z ∞ X

δ(E 0 − En )

−∞ n
Z
1 ∞

= −

π

−∞

φ∗n (r)φn (r) 0dEz − E0ImG(r, r 0 ; E 0 )/(z − E 0 )dE 0(69)

このように、解析接続された Green 関数は実軸上を除いて複素平面上いたるところで解析的である。

n(r) =X|φn (r)|2

n

=

Z EF
−∞

ρ(r, E)dE

Z

1 EF
ImG(r, r; E)dE
π −∞
Z
1
= −
ImG(r, r; z)dz
π A
Z
1
= −
ImG(r, r; z)dz
π C
= −

(70)

ここで、z = E + iΓ で Green 関数をみると、ImG(r, r; z) = ImG(r, r; E + iΓ)ZΓ
1 ∞ImG(E 0 ) 0
dE 0
=
π −∞
(E − E)2 − Γ2(71)つまり、実軸上にあったデルタ関数は幅 Γ だけひろがり、その結果状態密度の複雑な構造はなら
される。

3.4

散乱問題の形式解Schrödinger 方程式を特定の表示を使わず状態ベクトルとそれに作用する演算子とで書いておくと散乱問題の解を形式的に書き下すことができる。これにより散乱問題に対する見通しの良い描像がえられる。

3.4.1

Lippmann-Schwinger equation

全系のハミルトニアン H を

H = H0 + V
14

(72)

Contour Integration of the Green’s Function
Im E

Im E

Ef

Ef

C

A

Re E

B

Re E
C’

δ (E-En)Im G(z)2Γ
z=E+i ΓE
図 2: Contour integration of the Green’s Function

H0 = −∇2(73)

H|ψi = E|ψi(74)と書いておく。Schrödinger 方程式は、である。ところで、H0 の固有状態を |φi と書くと、Schrödinger 方程式の解が次のように形式的に書けることが直接の代入によって確かめられる。

|ψi = |φi +1
V |ψi
E − H0(75)

ここに出てきた特異な演算子に対しては、エネルギーをわずかに複素数にすることで意味を持たせることができる。

1
V |ψ ± i
E − H0±= |φi + G±
0 V |ψ i|ψ ± i = |φi +(76)

この方程式を Lippmann-Schwinger の方程式と呼ぶ。この方程式は特別な表示によらない。

15

3.4.2

散乱振幅

Lippmann-Schwinger 方程式を位置表示で書いてみる。
Z

hr|G0 |r0 ihr0 |V |ψidv 0hr|ψi = hr|φi +(77)

G0 の位置表示は Helmholtz 方程式の Green 関数であることがわかり、ポテンシャルが局所的であることから、√

hr|G0 |r 0 i
V̂01 ei E|r −r |= −4π |r − r 0 |(78)= δ(r − r 0 )V (r)(79)

であるので、

hr|ψi=

1
hr|φi −
4π

Z

√

0

ei E|r −r |V (r 0 )hr 0 |ψidv 0|r − r 0 |√1 ei Er→ hr|φi −hφ|V |ψi
4π r

(80)

となる。先に定義した散乱振幅 f (θ, φ) は、

f (θ, φ) = −1
hφ|V |ψi
4π(81)であることがわかり、これは t-matrix と同じものである。

3.4.3

Born Series

t-matrix の意味するところを直観的に知るために次のようなことを考える。
1
E − H + i²
1
E − H0 + i²

G =
G0 =

(82)
(83)

であるが、演算子に対する恒等式、

1
1
−
A B

=
=

1
1
(B − A)A
B
1
1(B − A)
B
A(84)

をもちいると、次のような式が成り立つことがわかる。

G0 = G − GV G0(85)

G = G0 + G 0 V G(86)

これを用いて、Lippmann-Schwinger 方程式は、|ψi = |φi + G0 V |ψi(87)= |φi + GV |φi

(88)

16

となる。よって散乱振幅は、

1 0
hφ |(V + V GV )|φi
4π1
= − hφ0 |T |φi4π= V + V GV

f (θ, φ) = −

T

(89)
(90)

G を逐次に展開してみると、
G = G0 + G0 V G(91)= G0 + G0 V G0 + G0 V G0 V G0 + ...

(92)= V + V GV

(93)

= V + V G0 V + V G0 V G0 V + ...

(94)= V + V G0 T

(95)であるから、

T

φ’

φ’

φ’

=

V

ψ

φ’

V
V G0

+

φ

+

φ

G
0
V
V
V G
0

+ ....

φ

図 3: Born Series
つまり、t-matrix はポテンシャルによる多段階の散乱過程を記述している。

3.4.4

Multiple Scattering

散乱中心がたくさんある時。

V

N
X

=

Vα

(96)

α=1

全体の t-matrix である T は、T

= V + V G0 V + ...

(97)

=

(98)

X

Vα +

α

=

X
α

X

Vα G0 Vβ + ...

α,β

Tα +

XX

T α G0 T β +α β6=αXX X
α β6=α γ6=β17

Tα G0 Tβ G0 Tγ + ...

(99)

のように、一つの粒子による散乱、２つの粒子による二重散乱、３個の粒子による散乱、というように、多重散乱の過程を表している。

Tα

T

=

+

Σα

Σ Σ
α β=α

Tβ
Tα

Tγ
Tα

+

Σ Σ
Σ
α β=α γ=βTβ+ ....

図 4: T-matrix of Multiple Scattering4

密度汎関数法3N の自由度を持つ多電子波動関数を扱う代わりに、３つの自由度しか持たない密度分布の汎関数として多電子系を記述することで問題は簡単になる。

4.1

Hohenberg-Kohn の定理

多電子系の基底状態の全エネルギーは電子密度の汎関数であり、電子密度について全エネル
ギーを最小化することで基底状態の全エネルギーが与えられることを示す。

4.1.1

多電子問題

外場 v(r) 中で互いに相互作用しながら運動している多電子系を考える。この系の HamiltonianHは、

H = T +V +U
T

=

X

(100)

−∇2i

(101)

v(ri )

(102)

i

V

=

X
i

18

U

=

X

2
|r − rj |
i≤j i

(103)

と書け、Schrödinger 方程式は、

HΦ = EΦ

(104)である。基底状態 Ψ は、次のような変分問題を解くことで与えられた。

HΨ = E0 Ψ(105)

E0 [Ψ] = MinhΨ|H|Ψi(106)hΨ|Ψi = N

(107)

さらに電子密度は、

ρ(r) = Ψ∗ Ψ

(108)で与えられるから、次のような mapping は明らかである。

{v} → {Ψ} → {ρ}(109)

つまり、多電子波動関数や電子密度、またエネルギー期待値は外部ポテンシャルの汎関数である。

4.1.2

Ψ = Ψ[v](110)ρ = ρ[v](111)

E0 = E0 [v](112)基礎関数としての密度

逆に、密度の汎関数としてポテンシャルを書くことができるか考える。ρ が異なる v, v 0 を与えたとする。v, v 0 が与える Schrödinger 方程式の解を Ψ, Ψ0 とする。それぞれの基底状態のエネル
ギー E, E 0 を調べてみる。

E 0 = hΨ0 |H 0 |Ψ0 i< hΨ|H 0 |Ψi= hΨ|H − v + v 0 |ΨiZ(v 0 − v)ρdv

= E+

(113)

E = hΨ|H|Ψi< hΨ0 |H|Ψ0 i= hΨ0 |H 0 − v 0 + v|Ψ0 i
Z= E0 +(v − v 0 )ρdv

(114)

上式の和は E 0 + E < E + E 0 を与え矛盾である。よって、ρ によって v は一意的に定まる。v から Ψ が決定されることから次のような mapping が証明された。

{ρ} → {v} → {Ψ}
19

(115)

つまり、外部ポテンシャルや多電子波動関数、また、さまざまな期待値は ρ の汎関数である。

v = v[ρ](116)Ψ = Ψ[ρ](117)hAi = hΨ[ρ]|Â|Ψ[ρ]i= A[ρ]
4.1.3

(118)

変分原理

ある外部ポテンシャルに対して、基底状態の全エネルギー汎関数 Ev は、

Ev = Ev [ρ]Z=ρvdv + F [ρ](119)

F [ρ] = T [ρ] + U [ρ](120)であり、正しい電子密度を与えれば Ev は基底状態の全エネルギーを与える。Ev が正しい電子密度で極小値をとることは次のことからわかる。あるポテンシャル v に対する基底状態 Ψ は、粒子数一定の条件で全エネルギーを最小化することから求められた。Ψ0 が異なったポテンシャルについての基底状態を表す時、同じポテンシャルに対して、
Z
0

0

0

0

0

Ev [Ψ ] = hΨ |v|Ψ i + hΨ |T + U |Ψ i =≥ Ev [Ψ]vρ0 dv + F [ρ0 ]Z= hΨ|v|Ψi + hΨ|T + U |Ψi =vρdv + F [ρ](121)が成り立つから、Ev [ρ] ≤ Ev [ρ0 ](122)である。ここで F [ρ] は、ポテンシャルに依らない普遍的な汎関数であるから、具体的に F [ρ] の表式が与えられれば多電子問題は密度に関する変分問題として簡単に解けるということになる。

4.2

Kohn-Sham 方程式

Hohenberg-Kohn が導いた変分問題は、有効ポテンシャル中を運動する相互作用のない多体系の Schrödinger 方程式をセルフコンシステントに解く問題にすりかえることができる。
全エネルギー汎関数をZ

Ev [ρ] = T [ρ] +vρdv + U [ρ]

Z

= T+

vρdv +Z= Ts +

12Z

ρ(r)ρ(r0 )drdr0 + EXC|r − r0 |vρdv + EH + EXC + (T − Ts )Z

= Ts +vρdv + EH + ẼXC
20

(123)

と変形する。ここで、EH , EXC , Ts はそれぞれ、Hartree エネルギー、交換-相関エネルギー、独立粒子の運動エネルギーである。電子間の静電相互作用 (Hartree エネルギー) 以外の多体効果は
全て交換-相関エネルギーにまとめてしまった。
Ts は、一体の波動関数 ψi から、
X

Ts = hΨ| −= −XZ∇2i |Ψii

ψi∗ ∇2i ψi dv(124)

ψi∗ ψi(125)iと計算できる。また密度は、
X

ρ =

i

つぎに、

Ev [ρ] = EB [ρ] + ED [ρ]

(126)

Z

EB [ρ] = Ts +Vef f ρdv

Z

Z

ED [ρ] = −Vef f ρdv +(127)vρdv + EH + ẼXC

(128)と、有効ポテンシャル Vef f を導入する。

δEv
=0
δρ

(129)

は、次の二つの連立方程式となる。

δEB
δρ
δED
δρ

= 0

(130)

= 0

(131)

一つめの式から有効ポテンシャル中を運動する一体の Schrödinger 方程式が得られ、二つめの式から有効ポテンシャルを密度から計算する式が得られる。まとめると、−∇2i ψi (r) + Vef f (r)ψi (r) = ²i ψi (r)Z

Vef f (r) = v(r) +X

ρ(r) =(132)
ρ(r0 )|r − r0 |dv 0 + VXC (r)(133)|ψi |2(134)

δ ẼXC [ρ]
δρ

(135)

i

VXC (r) =これらの方程式を Kohn-Sham 方程式と呼ぶ。

21

4.3

局所密度近似一様電子ガスについて、

ẼXC [ρ]/N= εXC (ρ)(136)

はわかっている。そこで、一様でない電子系の交換-相関エネルギーを、局所的には一様な電子ガスで近似できるとして次のように書く。
Z

ẼXC [ρ] =ρ(r)εXC (ρ(r)) dv(137)

よって、交換-相関ポテンシャルは、

VXC (r) = εXC (ρ(r)) + ρ(r)dεXC (ρ(r))
ρ(r)

と書ける。

4.4

逐次近似

Kohn-Sham の理論によればバンド計算は次のように進めることになる。
• Vef f を仮定。
• 一体 Schrödinger 方程式を解く。
P• fermi level まで、

2i |ψi | を実行して ρ(r) を求める。

• ρ(r) に対して Poisson 方程式を解き、Hartree ポテンシャルを求める。
• Vef f を LDA に依って求める。
• 新しい Vef f を古い Vef f と比べて Vef f を変更する。
以上を収束が得られるまで繰り返す。

22

(138)

START23

図 5: バンド計算の流れ図

Construct starting potenpotenv.f
tial.

gsdatp.f atmicv.f

Atomic calculation.

No

Starting potential
exists?

gtchst.f

Structural constants.

Configure lattice vectors
and the reciprocal lattice
vectors.
genlpt.f madlng.fchklat.f prmvec.f
bzmesh.fWhich System ?

specx.f readin.f

Yes

Check neutrality andneutrl.f
shift EF.

gensm.f chrdnc.f rrcnv.f

Energy integral of theGreen function and get
charge and spin densities.

spckkr.f kkrsed.f
cinvrx.f bzmsmb.f

Back scattering.

Radial wave functions
and the atomic phase
phasea.f radial.f
shifts.

Define complex energy
drvmesh.f
contour.

Self-consistent
iteration starts.

Fix the core state conficorcnf.f
gulation.

No

END

dsenum.f hypera.f

Display results.

Iteration loop ends.

Yes

erranc.f

Check convergence.

banden.f totalw.f

Total energy calculation.

potenv.f

Solve the Poissn equation
and add the exchangecorrelation potential toconstruct new potential.

5プログラムの構造

KKR コードは次のような流れ図にそって計算を実行している。主なザブルーチンについて、その内容を簡単にまとめておく。

表 2: 主なサブルーチンの概要
サブルーチン内容

specx.freadin.f
spmain.fchklat.fprmvec.fbzmesh.fgenlpt.fmadlng.f
gtchst.fgsdatp.f
atmicv.f
corcnf.f
drvmsh.fphasea.fradial.ffczero.f
cstate.f
spckkr.fkkrsed.fcinvrx.f
bzmsmb.fchrdnc.fneutrl.fpotenv.fbanden.f
totalw.ferranc.fdsenum.f
hypera.f

パラメーター、配列の定義、readin.f, spmain.f の呼びだし。
インプットファイルの読み込み。
KKR のメインプログラム。
与えられたマフィンティン半径が適当か調べる。
基本並進ベクトル、逆格子ベクトルを生成する。
第一ブリルアンゾーンの既約部分に k 点を定義する。
実空間、逆格子空間の格子点を生成する。
マーデルング定数の計算。
構造定数の計算。
ポテンシャルの初期値がなければ原子の計算から初期値を決める。
原子の電子状態の計算。
core の電子配置の設定。
Green 関数の積分路の定義。
t-matrix の計算。
動径 Schrödinger 方程式を解く。
Single Site Green’s function の計算に使う。
core 電子について動径 Schrödinger 方程式を解く。
Crystal Green’s Function の計算。
構造定数の計算。
(1 − gt)−1 の計算。
k 空間での積分。
電荷密度分布の計算。
チャージニュートラリティーを計算して、EF のシフトの大きさを決める。
poissn 方程式を解く。局所近似に基づいて交換相関ポテンシャルを求める。
エネルギー固有値の足しあげから、運動エネルギーを計算する。
全エネルギーの計算。
収束の判定し、新しい入力ポテンシャルを作る。
計算結果の表示。
超微細場の計算。

24

6

計算の実行

ここでは、実際に計算を始めるために必要なインプットファイルの書き方、配列の大きさに関する KKR プログラム内部のパラメーター、得られるアウトプットの見方を説明する。

6.1

エネルギーの原点

KKR 法ではマフィンティン近似が採用されており、球間には一定のポテンシャルを仮定している。そこをエネルギーの原点にとる。スピン分極がある時はスピンアップとスピンダウンでエ
ネルギーの原点が違う。概念図を図にしておく。

EZ

Ef

0

EW
EZ

band state

Energy contour (ewidth)

Re E

V(r)

valence statecore state

Nucleus

図 6: エネルギーの原点

6.2

インプットファイル

インプットファイルにはおもに結晶格子に関する情報が書かれる。鉄ニッケル合金の例を次にあげる。先頭が c で始まる行はコメント行である。区切りは空白またはカンマである。それぞれのフィールドの内容をコメント行に略号で示してある。また、各フィールドの内容を表にまとめておく (readin.f 参照)。値に数値を入力するところにはだいたいの目安を与えておいたが、それらは計算の進行状況を見て最適な値に設定すべきである。

25

c----------------------- input data ------------------------c
go
file
go
data/fenic ----------------------------------------------------------c
brvtypac/ab/a alpha beta gamma
bcc
5.3056 ,,,,,,c ----------------------------------------------------------c
edelt
ewidthreltypsdftypmagtyprecord
0.001
1.2
nrl
mjw
mag
init
c ----------------------------------------------------------c
outtypbzqltymaxitrpmix
update480
0.024
c ----------------------------------------------------------cntyp
1c ----------------------------------------------------------c
type
ncmp
rmt
field
l_max
anclr
conc
FeNi
2
0
0
2
26
60
28
40
c ----------------------------------------------------------cnatm
1c ----------------------------------------------------------c
atmicx
atmtyp000FeNic -----------------------------------------------------------

次に、必要と思われるものに解説を加える。適当でないパラメーターを与えた場合、それぞれのサブルーチンでエラーメッセージが出されるのでそれを参考にして正しく修正する。

• go dsp, dos,spc のモードでは、record, outtyp, maxitr はそれぞれ自動的に 2nd, quit, 1 に設定される。
• file カレントディレクトリからの相対パス名/ファイル名. 拡張子 のように指定できる。
• brvtyp, 格子定数 brvtyp から明白なものに関しては、c/a, alpha などの入力は省略できる。
• edelt このプログラムで計算される Green 関数は、複素エネルギー上面の正の実軸がブランチカットになっており Green 関数は定義されない。そこで複素平面上の虚数部分の正の側から実軸に近づいた極限として実エネルギーでの Green 関数が定義される。数値計算上は有限の虚数部分を与えており、その大きさが edelt である。
• ewidth Green 関数が計算されるエネルギー領域の広さである（図 6 参照）。ewidth の大きさは価電子帯を完全にカバーするものでなければならない。価電子帯の幅がわからない場
合には適当な ewidth を用いて荒い計算をしておき、状態密度を出力して価電子帯の幅を見積もるようにすればよい。また、Green 関数の計算に必要な Wronskian や動経波動関数はエネルギーが EW±EZ の範囲にとられたサンプリングポイントで計算され、チェビシェフ
多項式を用いて近似される。それを用いて Energy Contour 上で Green 関数が計算されるため、Energy Contour が EW±EZ の範囲に収まるようにしなくてはならない。EW,EZ は自動的に設定されるが、計算の過程でフェルミレベルは動いていくので、Energy Contour26表 3: インプットファイル入力パラメータ（１）
フィールド

値

意味、単位

go

go
ngo
dos
dsp
spc
ファイル名
fcc
bcc
hcp(hex)
sc
bct
st
fco
bco
bso
so
bsm
sm
trc
rhb(trg)fcta（格子定数）
c/a 比
b/a 比
α
β
γ
∼0.001—nrl
sra
mjw, vbh, vwn
mag, nmag
nmag-mag(rvrs)
kick
init
2nd
1st

（処理の指定）バンド計算を行なう。
何もしない。
状態密度の表示を行なう。
既にある計算の結果を表示する。
Bloch spectral function（dispersion relation）の表示。
指定したファイル名でポテンシャルデータを記録する。
（ブラベー格子の指定）face centered cubic
body centered cubic
hexagonal close packed
simple cubic
body centered tetragonal
simple tetragonal
face centered orthorhombic
body centered orthorhombic
base centered orthorhombic
simple orthorhombic
base centered monoclinic
simple monoclinic
triclinic
rhombohedral (trigonal)
face centered tetragonal
原子単位。

file
brvtyp

a
c/a
b/a
alpha
beta
gamma
edelt
ewidthreltypsdftypmagtyprecord単位は度。

フェルミエネルギーに付与される虚数部分。単位は Ry。
Energy Contour の幅。図 6 参照。単位は Ry。
相対論を考慮しない。
相対論を考慮する。
交換相関ポテンシャルの指定。
スピンを考慮して magnetic な計算を行なう。
非磁性の計算を行う。
磁化の符号を反転する。
人為的に対称性を壊し、磁性状態へ移行させる。
原子の計算から入力ポテンシャルを決める。
最新のポテンシャルデータから計算を続ける。
２番目に新しいデータから計算を続ける。

27

表 4: インプットファイル入力パラメータ（２）
フィールド

値

意味、単位

outtypupdatequit
0, 1, 2,...
t, s, m, l, h, u
≥10.01-0.03
≥1
タイプ名
≥1
≥0
field の大きさ
≥0
≥0
—
≥1
位置ベクトル
タイプ名

ポテンシャルデータを更新する。
計算結果は捨てて更新をしない。
大きいほどブリルアンゾーンのメッシュは細かい。
nfqlty.f で与えられる整数に対応する。
iteration の最大回数。
ポテンシャルの混合のパラメター。
単位格子中の独立なサイトの種類の数。
個々のサイトの種類の名前。
そのサイトを構成する原子の種類の数。
マフィンティン半径 (a 単位)。
外部磁場 (Ry)。
l=0 から l=l max まで考慮する。
原子番号。
そのサイト中でその原子が占める割合。
単位格子中の原子の個数。
原子の位置。
その位置にあるサイトの種類の指定。

bzqltymaxitr
pmix
ntyp
type
ncmp
rmt
field
l max
anclr
conc
natm
atmicx
atmtypが計算の途中で EW±EZ からはみ出す可能性がある。その場合計算は終了するが、go=2nd
でもう一度計算すれば新しい EW,EZ が設定される。

• record 異常終了などで計算結果が失われてしまわないように、前回のポテンシャルのデー
タはバックアップされ、今回得られた新しい計算結果と一緒に保存されている。2nd は最新のデータを、1st はその前のデータをイニシャルポテンシャルとして用いることを示す。
• bzqlty 与えた数字に応じて第一ブリルアンゾーンの既約部分に k 点 (逆格子空間でのサン
プリングポイント) がとられる。数字に対応していくつの k 点がとられるかは場合による。
• rmt 与えたマフィンティン半径では球どうしが重なってしまう場合、マフィンティン半径は与えた半径比で球が接するように設定し直される。0 を与えた時または省略した時は半径比として原子半径比がとられる。
• anclr 0 は空格子点をおくことを意味する。
• conc concentration はサイト内で規格化されるので総和はいくらになってもよい。
• atmicx a を単位とする。直交座標で指定する時は例えば 0.5, 0.5, 0.5 などと書き、基本ベクトルで指定する時は 0.5a 0.5b, 0.5c などと書く。
また、ひとつのファイルの中に複数のインプットファイルを続けて書いておくと計算は順番に進められる。

28

6.3

配列の大きさを指定するパラメータ

KKR プログラムで、主な配列はプログラム specx.f で宣言されるので、specx.f 中のパラメー
タを変更してやることで、必要な大きさの配列が確保できるようになっている。specx.f 中で指定できるパラメータについて表にまとめる。
表 5: 配列に関するパラメータ

6.4

パラメータ

意味

natmmxncnpmxmsizemx
mxlmx
nk1mx,nk3mx
msex
ngmx
nrpmx
ngpmx
msr
mse0
tol

原子の数の最大値。
原子の種類の数の最大値。
P
2
Green 関数行列の最大サイズ。 natm
i=1 (lmax,i + 1) 以上。
mxl-1 の角運動量までで展開する。mxl の最大値。
nk1mx+nk3mx が k 点の最大数。
エネルギーメッシュの点数の最大値。go=dos ではこの点数を使う。
チェビシェフ展開の次数。
実空間の格子点の最大数。
逆空間の格子点の最大数。
動径座標のメッシュの点数。
（deta 文）go=dos 以外でのエネルギーメッシュの点数。
（deta 文）ポテンシャルが tol まで収束したら計算を終了する。

アウトプット

さきほど例としてあげた鉄ニッケル合金のインプットファイルを使って計算した結果出力されるアウトプットを見てみる。
2-Feb-2005
meshr mse
400
35

ng

mxl
15

3

data read in
go=go
file=data/feni
brvtyp=bcc a= 5.30560 c/a=1.00000 b/a=1.00000
alpha= 0.0 beta= 0.0 gamma= 0.0
edelt= 0.0010 ewidth= 1.200 reltyp=nrl
sdftyp=mjw
magtyp=mag
record=init outtyp=update bzqlty=4
maxitr= 80 pmix=0.02400
ntyp= 1 natm= 1 ncmpx= 2
complex energy mesh1( -1.2000, 0.0000)
4( -1.1970, 0.0121)
7( -1.1710, 0.0368)
10( -1.0562, 0.0779)
13( -0.7577, 0.1158)
16( -0.3646, 0.1104)
19( -0.1156, 0.0708)
22( -0.0276, 0.0360)
25( -0.0057, 0.0165)
28( -0.0011, 0.0073)
31( -0.0002, 0.0032)
34( 0.0000, 0.0014)

2( -1.1998, 0.0030)
5( -1.1928, 0.0186)
8( -1.1478, 0.0489)
11( -0.9792, 0.0930)
14( -0.6240, 0.1199)
17( -0.2591, 0.0988)
20( -0.0735, 0.0575)
23( -0.0165, 0.0280)
26( -0.0033, 0.0126)
29( -0.0006, 0.0055)
32( -0.0001, 0.0024)
35( 0.0000, 0.0010)

3( -1.1990, 0.0070)
6( -1.1849, 0.0267)
9( -1.1111, 0.0629)
12( -0.8787, 0.1063)
15( -0.4891, 0.1179)
18( -0.1764, 0.0850)
21( -0.0455, 0.0458)
24( -0.0098, 0.0216)
27( -0.0019, 0.0096)
30( -0.0004, 0.0042)
33( -0.0001, 0.0018)

29

file to be accessed=data/feni
createdlattice constantbravais=bcc
a= 5.30560
c/a=1.0000
b/a=1.0000
alpha= 90.00
beta= 90.00
gamma= 90.00
primitive translation vectors
a=( -0.50000 0.50000 0.50000)
b=( 0.50000 -0.50000 0.50000)
c=( 0.50000 0.50000 -0.50000)
type of site
type=FeNi

rmt=0.43301 field= 0.000
lmxtyp= 2
component= 1 anclr= 26. conc= 0.6000
component= 2 anclr= 28. conc= 0.4000

atoms in the unit cell
position=
0.00000000

0.00000000

0.00000000

type=FeNi***msg in spmain...new ew, ez generated
ew=
0.09998 ez=
0.80300
preta= 0.35542 eta= 0.35542
isymop= 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
last= 243

np= 15nt= 141

nrpt= 169

nk=

29

itr= 1rms error = -1.872
itr= 2rms error = -2.390
itr= 3rms error = -2.879
itr= 4rms error = -3.357
itr= 5rms error = -3.890
itr= 6rms error = -3.879
itr= 7rms error = -4.878
itr= 8rms error = -4.850
itr= 9rms error = -5.574
itr= 10rms error = -6.235
interval= 10
cpu time=
0.00 secnl
cnf
energy----------------------------------1s
2.000-508.5203
2s
2.000-59.20742p
6.000-51.18063s
2.000-6.80273p
6.000-4.45633d
6.000-0.66964s
2.000-0.4930
itr= 1rms error = -1.724
itr= 2rms error = -2.239
itr= 3rms error = -2.657
itr= 4rms error = -2.967
itr= 5rms error = -3.645
itr= 6rms error = -3.695
itr= 7rms error = -4.268
itr= 8rms error = -5.061
itr= 9rms error = -5.356
itr= 10rms error = -5.470
itr= 11rms error = -6.321
interval= 11
cpu time=
0.01 secnl
cnf
energy----------------------------------1s
2.000-595.8142
2s
2.000-70.7050

30

nd=

1

2p
3s
3p
3d
4s

6.000
2.000
6.000
8.000
2.000

-61.8163-7.9863-5.2730-0.7806-0.5252

record 1 will be overlaied by input andrecord 2 will be replaced by new output.
core configuration for Z= 26
state 1s 2s 2p 3s 3p 3d 4s 4p 4d 5s 5p 4f 5d 6s 6p 5f 6d 7sup
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
down
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
core configuration for Z= 28
state 1s 2s 2p 3s 3p 3d 4s 4p 4d 5s 5p 4f 5d 6s 6p 5f 6d 7sup
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
down
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0***** self-consistent iteration starts *****
FeNi
itr= 1 neu= -1.4601 moment= 0.0905 te= -2718.3225808 err= 0.119
itr= 2 neu= -0.8991 moment= 0.1151 te= -2718.3221114 err= -0.277
itr= 3 neu= -0.2974 moment= 0.1168 te= -2718.3157206 err= -0.149
itr= 4 neu= 0.1787 moment= 0.1296 te= -2718.3197036 err= -0.627
itr= 5 neu= 0.5050 moment= 0.1933 te= -2718.3214574 err= -0.407
itr= 6 neu= 0.6526 moment= 0.2873 te= -2718.3261849 err= -0.483
itr= 7 neu= 0.5415 moment= 0.4203 te= -2718.3258145 err= -0.530
itr= 8 neu= 0.3129 moment= 0.5384 te= -2718.3258235 err= -1.056
itr= 9 neu= 0.1143 moment= 0.7051 te= -2718.3265769 err= -0.801
itr= 10 neu= -0.0203 moment= 0.8951 te= -2718.3295190 err= -0.873
itr= 11 neu= -0.1221 moment= 1.0540 te= -2718.3323905 err= -0.967
itr= 12 neu= -0.1746 moment= 1.2526 te= -2718.3352375 err= -0.844
itr= 13 neu= -0.1626 moment= 1.4605 te= -2718.3369449 err= -0.608
itr= 14 neu= -0.1454 moment= 1.6426 te= -2718.3395758 err= -0.754
itr= 15 neu= -0.1518 moment= 1.6880 te= -2718.3406928 err= -0.650
itr= 16 neu= -0.1663 moment= 1.7401 te= -2718.3410295 err= -0.607
itr= 17 neu= -0.1522 moment= 1.7988 te= -2718.3415185 err= -1.255
itr= 18 neu= -0.1303 moment= 1.7861 te= -2718.3414079 err= -1.326
itr= 19 neu= -0.0791 moment= 1.7878 te= -2718.3410749 err= -1.206
itr= 20 neu= -0.0246 moment= 1.7774 te= -2718.3410051 err= -1.553
itr= 21 neu= 0.0213 moment= 1.7734 te= -2718.3409834 err= -1.224
itr= 22 neu= 0.0494 moment= 1.7695 te= -2718.3410544 err= -1.282
itr= 23 neu= 0.0583 moment= 1.7667 te= -2718.3410975 err= -1.829
itr= 24 neu= 0.0518 moment= 1.7646 te= -2718.3410481 err= -1.769
itr= 25 neu= 0.0361 moment= 1.7609 te= -2718.3409914 err= -1.879
itr= 26 neu= 0.0172 moment= 1.7587 te= -2718.3409593 err= -2.014
itr= 27 neu= 0.0005 moment= 1.7597 te= -2718.3409526 err= -1.849
itr= 28 neu= -0.0111 moment= 1.7612 te= -2718.3409625 err= -2.108
itr= 29 neu= -0.0169 moment= 1.7625 te= -2718.3409716 err= -2.097
itr= 30 neu= -0.0176 moment= 1.7634 te= -2718.3409749 err= -2.114
itr= 31 neu= -0.0147 moment= 1.7639 te= -2718.3409743 err= -2.287
itr= 32 neu= -0.0102 moment= 1.7640 te= -2718.3409723 err= -2.403
itr= 33 neu= -0.0055 moment= 1.7640 te= -2718.3409707 err= -2.315
itr= 34 neu= -0.0017 moment= 1.7638 te= -2718.3409703 err= -2.499
itr= 35 neu= 0.0009 moment= 1.7637 te= -2718.3409703 err= -2.727
itr= 36 neu= 0.0024 moment= 1.7635 te= -2718.3409704 err= -2.832
itr= 37 neu= 0.0029 moment= 1.7633 te= -2718.3409703 err= -2.844
itr= 38 neu= 0.0026 moment= 1.7632 te= -2718.3409701 err= -2.749
itr= 39 neu= 0.0019 moment= 1.7631 te= -2718.3409699 err= -2.813
itr= 40 neu= 0.0010 moment= 1.7631 te= -2718.3409698 err= -2.999
itr= 41 neu= 0.0002 moment= 1.7631 te= -2718.3409698 err= -3.004
interval= 41
cpu time=
17.79 sec
sdftyp=mjw
reltyp=nrl
dmpc=0.024
FeNi
itr= 41 neu 0.0002 chr,spn 8.8000 1.7631 intc,ints 0.9632 -0.0234
rms err= -3.164 -3.274 -3.004 -3.177
ef=
0.7260630
0.7343370 def=
2.6962356 12.2016065

31

total energy=-2718.3409698*** type-FeNi
Fe (z= 26.0) ***core charge in the muffin-tin sphere =17.9787766
valence charge in the cell (spin up ) =0.19826(s)0.19326(p)4.36075(d)valence charge in the cell (spin down) =0.19758(s)0.20985(p)1.84235(d)total charge= 24.98082
valence charge (up/down)=
4.75227
2.24978
spin moment=
2.50249 orbital moment=
0.00000
core level (spin up )-507.1266873 Ry(1s)-5.5302084 Ry(3s)core level (spin down)-507.1176964 Ry(1s)-5.3171341 Ry(3s)-57.8685202 Ry(2s)-3.1878588 Ry(3p)-49.8296484 Ry(2p)-57.7567489 Ry(2s)-2.9808461 Ry(3p)-49.7425563 Ry(2p)hyperfine field of FeNi
-239.803 KG (core= -250.980 KG valence=
11.177 KG )core contribution-21.109 KG(1s)-548.149 KG(2s)318.277 KG(3s)charge density at the nucleus
11820.1880 (core=
11814.5387 valence=core contribution
10701.4478(1s)972.7303(2s)5.6493 )140.3607(3s)*** type-FeNi
Ni (z= 28.0) ***core charge in the muffin-tin sphere =17.9914848
valence charge in the cell (spin up ) =0.21450(s)0.20389(p)4.50139(d)valence charge in the cell (spin down) =0.22638(s)0.23660(p)3.74662(d)total charge= 27.12086
valence charge (up/down)=
4.91978
4.20959
spin moment=
0.71019 orbital moment=
0.00000
core level (spin up )-594.3771108 Ry(1s)-6.5907736 Ry(3s)core level (spin down)-594.3686783 Ry(1s)-6.5162732 Ry(3s)-69.2884479 Ry(2s)-3.8813868 Ry(3p)-60.3956628 Ry(2p)-69.2422357 Ry(2s)-3.8088998 Ry(3p)-60.3587004 Ry(2p)hyperfine field of FeNi
-166.611 KG (core=-77.215 KG valence=-89.396 KG )core contribution-9.105 KG(1s)-204.777 KG(2s)136.667 KG(3s)charge density at the nucleus
14827.9567 (core=14822.0421 valence=core contribution
13396.9377(1s)1243.1999(2s)cpu used

5.9145 )181.9046(3s)17.94 sec

日付、meshr, mse, ng, mxl の表示の後読み込んだインプットファイルの内容が確認される。次に用いる複素平面上のエネルギーメッシュが表示される。実軸 0 がフェルミレベルである。その後結晶格子や原子の情報が表示される。インプットファイルに書かれた設定が不適当である場合、修正した値が印字される。ew、ez、k 点の数 (nk)、独立な構造定数のブロック数 (nd) などが計算されて表示される。specx.f で設定した値を越えるとエラーメッセージを出して終了するのでspecx.f のパラメーターを設定し直してやり直す。
record=init のとき、または与えたポテンシャルファイルが存在しない時は原子の電子状態計
算から初期ポテンシャルが構成される。

32

次に計算で採用される core の電子配置が示される。core の電子についてはその原子核のまわりに十分良く局在していると考えてバンド計算には入れず、原子の問題として扱う。core 電子のレベルが Energy Contour にはいってしまった場合は、そのレベルは自動的に valence として扱われる。
iteration loop が始まると、ループの回数、チャージニュートラリティー、モーメント、全エネ
ルギー、エラーの値がループ毎に表示され計算が収束していく様子がわかる。エラーが設定値より小さくなるとセルフコンシステントな解が求められたとして計算を終了する。maxitr 以上ループしても収束が得られない場合は計算を打ち切る。
フェルミレベル、フェルミレベルの状態密度などの出力の後、それぞれの原子について各軌道の占有率、core の軌道のエネルギーレベルなど基本的な計算結果、そして得られた電子状態から計算した物理量が出力される。valence の電子のエネルギーバンドが完全に積分されていることに注意しなくてはならない。

7

計算例

7.1

Density of States

正しい計算が行なわれたかどうかの確認のために状態密度を表示してみることは重要である。
それだけでなく、状態密度とその各サイトへの射影、角運動量についての射影は固体の電子状態、化学結合を直接反映するので、その物性の直観的な理解に必要不可欠である。
状態密度の構造を見るためには Green 関数を複素平面上で実軸にそって眺める必要がある。インプットファイルで dos のモードを指定すると、drvmsh.f で cemesr.f が呼ばれ、実軸に沿って等間隔にメッシュがとられる。また、詳細な DOS を描くためには十分たくさんの k 点をとっておかなくてはならない。
通常のアウトプットに加えて、DOS の計算結果が表示される。表示されるエネルギーはフェル
ミレベルから測った値である。Valence の電子のバンドを完全にとらえているか注意する。アウ
トプットから計算値をコピーして適当なソフト (ワークステーション上では gnuplot など) でグラ
フにすると良い。例として、いくつかの代表的な固体についてその DOS を見てみる。
8
7
6

DOS

5
4
3
2
1
0-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

Energy (Ry)
図 7: Be の状態密度

33

0

0.2

0.4

35
30

paramagnetic

25
20
15
10
5
35
30

Spin Up

DOS

25
20
15
10
5
0
5
10
15
20
25
30

Spin Down

35-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

Energy (Ry)

図 8: Fe の状態密度

7.1.1

Beryllium

典型元素で金属となるものから例としてベリリウムをあげる。非磁性の金属であるから計算はnmag のモードを用いた (nk=512, mse=235)。自由電子の状態密度と似ているが、フェルミレベ
ルは状態密度の谷に位置していて、その上下で状態密度が増えている。つまり、原子が集合しそれぞれの軌道が重なり合い、結合軌道と反結合軌道が形成され、そのうちの結合軌道が占有されて結合を生じる。

7.1.2

Iron

典型的な強磁性金属の例として鉄についての計算結果を示す (nk=506, mse=235)。図の一番上は nmag のモードでの計算結果である。フェルミレベルは状態密度の高いところに位置している。mag のモードで計算すると片方のスピンの d バンドはほとんど占有され、フェルミレベルはもう一方のスピンのバンドの中央の谷に位置するようになる。この谷も直観的には結合、反結合
軌道の形成に対応している。

34

30
25

Spin Up

20
15

DOS

10
5
0
5
10
15
20
25

Spin Down

30-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

Energy (Ry)

図 9: Ni の状態密度

7.1.3

Nickelニッケルについての例を示す (nk=505, mse=235)。片方のバンドがほぼ占有され、もう一方のバンドにホールができているのは鉄のばあいと似ているが、フェルミレベル付近での状態密度は高く急激な変化を示している。このような場合フェルミレベルのすこしの変化に対して計算は敏感に反応するので慎重に進まなくてはならない。つまり、多くの k 点をとり、pmix を小さくとるのが望ましい。

7.2

スーパーセル法

結晶中に導入された不純物は不規則に分布しており、そのような系は周期構造を前提とするバ
ンド計算の対称とはならない。Green 関数法では、ホスト結晶の Green 関数とホスト原子、不純物原子の散乱行列から結晶中にただひとつ置かれた不純物の電子状態を求めたり、不規則合金の電子状態を求めたりする方法が確立されているが、局所的な電子構造を反映する物理量を議論するのに重要な格子緩和が取り入れられない。この困難を回避するのに次のような簡便な方法をとる。つまり、複数の単位格子を新しく単位格子と定義して、その中に不純物を導入する。こうして作った単位格子に対してバンド計算を行なう。不純物どうしの相互作用は十分小さいという前提のもとで、結晶中に周期的に不純物を置いてバンド計算を適用したということになる。

7.3

格子定数、格子緩和

密度汎関数法に基づき局所近似を適用して系の全エネルギーを計算しているが、格子定数や、
スーパーセル法でシミュレートされた不純物まわりの格子緩和を全エネルギー極小の条件から見積もることができる。このような計算を行なう場合、最適化するパラメーターに対する全エネル

35Pure System

With Impurities

Super Cell

図 10: Super Cell
ギーの相対的な変化が重要であるから、それぞれの計算においてその他のパラメーターを同じにしておかなくてはならない。とくに緩和を見積もる場合、マフィンティン半径は原子を動かしてもぶつからないよう小さくとっておく必要がある。

8

参考文献

最後に、このノートを書くにあたって特に参考にした文献を列挙しておく。KKR-Green 関数法を用いた計算結果に関する報告や超微細相互作用に関する文献は数多くあるがここにはあげなかった。全体にわたって、1995 年度から始められた赤井教授の南園研究室に於けるゼミのノート、大学院の講義、ソースコードのコメントが参考にされているのはいうまでもない。

• 散乱理論– J. J. Sakurai 現代の量子力学 (吉岡書店)– L. I. Schiff 量子力学 (吉岡書店)– J. Alfken 基礎物理数学 (講談社)– 砂川重信 散乱の量子論 (岩波書店)• KKR 法– W. Kohn, N. Rostoker Phys. Rev. 94(1954)1111.
– K. H. Johnson J. Chem. Phys. 45(1966)3085.
– F. Herman, A. R. Williams, K. H. Johnson J. Chem. Phys. 61(1974)3508.
– R. Podloucky, R. Zeller, P. H. Dederichs Phys. Rev. 22(1980)5777.
– R. Zeller, J. Deutz, P. H. Dederichs Solid State Comm. 44(1982)993.
– H. Akai’s note on Korringa-Kohn-Rostoker Method (1993).
– T. Kotani, H. Akai Phys. Rev. B54(1996)16502.
– J. C. Slater 分子軌道計算 (東京図書)• 密度汎関数法– P. Hohenberg, W. Kohn Phys. Rev. 136(1964)B864.
– W. Kohn, L. J. Sham Phys. Rev. 140(1965)A1133.

36