CMP Technical Report No. 92
「
DOS
の原子軌道への分解」
白井光雲
大阪大学・産業科学研究所
２００６年８月１８日
Department of Computational Nanomaterials Design
ISIR, Osaka University
目 次
1 はじめに 1 2 理論 1 3 GaAsの例 3 3.1 SCF計算 . . . . 3 3.2 DOS計算 ーpwbcdの使い方 . . . . 4 3.3 tetrahedron法 ーpdosdrの使い方 . . . . 6 3.3.1 pwbcd . . . . 7 3.3.2 pdosdr . . . . 9 3.3.3 aydosp . . . . 11 3.4 原子軌道展開の取り方 . . . . 12 3.5 pDOSの表現の仕方 . . . . 14 4 pdosdrでのk点メッシュの取り方 14 5 まとめ 16 A label of m components 17 B colors in aydosp 171
はじめに
波動関数の原子軌道への分解は既にOsaka2kで実現されており、Technical Reportシリー
ズの中で何回か取り上げられた[1, 2, 3]。それぞれを相互に参照しあったりするのが大変で あり、また時代が変遷するたびにいろいろ改良が加えられ、前に書いていたことと矛盾する 部分も出てきて、このあたりで統一的に整理する必要があるように思える。それが本レポー トの目的である。 このレポートで使用する機能はpwmのバージョンで VersionNo = ’1.72’ VersionDate = ’08 Jul 2006’ 以降である。そしてpwbcdのバージョンでは VersionNo = ’0.90’ VersionDate = ’27 Jun 2006’ 以降である。
2
理論
κ種の原子の軌道角運動量がlmの波動関数 ϕκlm(r) = uκl(r)Ylm(θ, φ) (1) は擬ポテンシャルを作成するプログラムで得られる。 固体中の波動関数を原子の軌道関数で分解するという場合、基底系の波動関数は、純粋の 孤立原子の波動関数を取り扱うのは非常に困難であるので、そのかわりそのブロッホ和を 取ったものを用いることとする。これは固体中の原子軌道というものをそういうふうに取っ たということであり、定義である。これがベストという保証は無い。ともかく、そうすると 基本格子でa番目の原子（κ種）のlm軌道は ϕk,alm(r) = X R exp[ik· (R + Ra)] ϕκlm(r− R − Ra) (2) とできる。これを結晶中における原子軌道と見なす。 固体中の波動関数ψk(r)のこの原子軌道成分は、差し当たりそれらの内積を取ることで評 価できる。 hψk|ϕκlmi (3) 式（3）の評価には、実空間表示で実空間上での積分で行うこともできるが、G空間での ベクトル同士の内積と取るのが計算効率からいって得策である。すなわち X G hψk|k + Gihk + G|ϕκlmi = X G ψk∗(k + G)ϕκlm(k + G) (4)で評価する。 平面波展開法では初めから波動関数はψk∗(k + G)として表されているのでこちらは問題 ない。問題は原子軌道を平面波展開することである。 ϕκlm(q) = Z ϕκlm(r)e−iq·rdr = (2l + 1)il Z ∞ 0 uκl(r)jl(qr)r2dr Z Ylm(ˆq)Pl(cos θ)dΩ = 4πil Z ∞ 0 uκl(r)jl(qr)r2drYlm(ˆq) (5) となるので、1 Ilκ(q) = 4πil Z ∞ 0 uκl(r)jl(qr)r2dr (6) とすると、{Ilκ(q)}のセットを求めておくと、結局式（4）の評価は、a番目の原子の構造因 子Sa(q)も含めて hψk|ϕκlmi = X q
ψ∗k(q)Ilκ(q)Ylm(ˆq)Sa(q) (7)
でできる。なおこの式（7）から分かるように結晶中の原子軌道ϕκlmiは元の孤立した原子の 波動関数の動径成分Ilκ(q)、角運動量成分Ylm(ˆq)、そして構造因子Sa(q)からなっているこ とが分かる。 しかしながらこういうふうに原子軌道ϕk,almに分解したときの問題点として、 • もともと原子軌道は直交系でない。そのためこのような分解した後の成分は、当然そ の和をとっても1とならない。 • また平面波展開された原子軌道は規格化さえされていない。もちろん展開された平面 波空間の中で再規格化することはできるが、わざとそうしていない。絶対値を比較し たいときがあるからである。 したがって、その成分はあくまで相対的な意味しか持たないことに注意すべきである。 完全なブロッホ和の組 もう一つの注意として、式（2）の中の波数kについてであるが、 それは平面波展開の時は第一ゾーンの中だけで取れば十分であった。つまり ψk(r) =X G ck+Gei(k+G)·r (8) の展開ではkを第一ゾーンに限定しても、それより高いエネルギーの成分は後のGの和の 中に含まれているので式（8）だけで励起状態も全て記述できる。しかしブロッホ和（2）の 1 本来はYlm(ˆq)にはマイナスがつくはずだが、実数化した球面調和関数を使っているので問題ない。
場合、kをそのように制限してしまうと、表現できる波動関数が限定されてしまう。平面波 展開のときと違い、展開の和の中には励起状態が入っていないので、励起状態を表現できな くなってしまう。つまりkをそのように取ってしまうと、{ϕk,alm(r)}は完全系ではなくな る。ブロッホ和（2）の場合はこうしてkを第一ゾーンに限定するべきではないように思え る。これについては後で実例で調べる。
3
GaAs
の例
GaAsを例に取り、はじめに従来の全DOSの計算例を示し、その後で部分DOSへの分解
の仕方を説明する。
3.1
SCF
計算
SCF計算の条件は以下のように取る。
==================== K-Space Setup ============================== k sampling point set
Nkpts = 2
No NM index in p in c A/gmin i/o star WTK
1 LD -1 -1 -1/ 4 -1 -1 -1/ 4 0.25000 1 4 0.25000
2 XY 1 -1 -1/ 4 -3 1 1/ 4 0.47871 1 12 0.75000
==================== PW_Expansion ============================== Cutoff in the reciprocal space
am = 3.10000 (rel. units)
kcut = 3.15801 (ab^-1) with 2Pi
Ecut = 9.97300 (Ry) UNIT of K 1.01871 (a.u.) Planewave expansion with NHDIM = 169 Name: LD -1 -1 -1/ 4 Nstr= 4 WTK= 0.250000 NPW= 162 Name: XY 1 -1 -1/ 4 Nstr= 12 WTK= 0.750000 NPW= 164 ポテンシャルの情報も記す。 ==================== atomic PPOT ==============================
Ga ca nrl nc pseudopotential read from tape
atom-lda 8-JAN- 2 Improved Troullier - Martins potential
4s( 2.00) rc= 2.50 4p( 1.00) rc= 2.35
As ca nrl nc pseudopotential read from tape
atom-lda 8-JAN- 2 Improved Troullier - Martins potential
4s( 2.00) rc= 2.50 4p( 3.00) rc= 2.50
SCF計算の結果は以下のようになっている。
====================SCF calculation==================== convergence parameters
max iteration : 15
iter Eel deE Xsi nst/bk aglmax
(Ry/cell) (Ry/cell) (Ry^2/cell)
1 -0.0090254029 -6.9760E+00 7.0425E-02 5/ 0 0.624308035 2 -0.4061363678 -3.9711E-01 9.0575E-04 5/ 0 0.237960989 3 -0.4127710828 -6.6347E-03 3.2210E-05 5/ 1 0.021868549 4 -0.4129192995 -1.4822E-04 6.4211E-07 5/ 7 0.003967351 5 -0.4129231853 -3.8858E-06 2.5645E-08 5/ 16 0.000715309 6 -0.4129233887 -2.0347E-07 5.1295E-10 5/ 11 0.000085793 7 -0.4129232409 1.4781E-07 3.4848E-11 3/ 8 0.000021200
CG process is stopped because increase in Eel 1.4781E-07
Etot Eel delta E resid iter
============== ============= ========== ========== === -17.25697453 -0.41292324 1.478E-07 3.485E-11 7 ===== KS levels ===== -0.513307 0.047161 0.283862 0.283862 -0.417444 -0.103171 0.106391 0.193900
3.2
DOS
計算 ー pwbcd の使い方
次にDOS計算のやり方を説明するが、はじめは従来の全DOS計算のやり方を説明する。 これはOsaka2kの基本マニュアルに書かれていることであるので復習する意味合いでここ に述べる。 DOS計算プログラムpwbcdへの入力ファイルbcd.paraは JobType dosInput file name gaas.prim
number of division (nkdiv) 8
number of levels you want to draw (NBUP) usually NEPC 8
scan zone only (iscan) 1
print control (ilp) 1
use symmetry (isymm) 1
energy unit (ienun=0 for Ry, 1 for eV) 0
のようにする。
pwbcdの計算結果はファイルdos gaas.outに出力される。まず
========================= DOSMain ==============================
nkdiv per line : 8
nk3dim : 512
N of bands (NBUP) : 8
dim0eng : 4096
Use of symmetry for diagonalization
Decomposition of DOS ON
using pa
n sample k points (nktot: 47 n total k points (nkfull: 564
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Energy \Ry.\ 10 20 30 40 50 States/Ry. DOS of GaAs
図1: Total DOS of GaAs.
の見出しの後、DOS計算の条件が並べられている。
その最後にフェルミ準位の情報などが記されている。 ========================= FindEf ==============================
Nel0prim = 8
fixocc = 2.00000E+00
N of electrons 8 are occupied in the ascending order in energy. N of occupation processe 189
residual of electron num 0.000000 === Fermi Sort Summary ===
Fermi_Level = 0.424483 (Ry) at 1 th k-point GM 0 0 0/ 8 (p) 0 0 0/ 8 (c) Valence Top = 0.424483 (Ry) at 1 th k-point GM 0 0 0/ 8 (p) 0 0 0/ 8 (c) Conduction Bottom = 0.428412 (Ry) at 42 th k-point X 4 4 0/ 8 (p) 0 0 8/ 8 (c) Energy Gap = 0.003930 (Ry)
with energy resolution for degeneracy 1.00E-03 (Ry)
pwbcdの出力ファイルのうちfort.2を用いてpdosdrに掛ける。まず従来の処理方法で ある、全DOSを描くことを行う。 pdosdrへの入力ファイルpdosdr.inpは gaas. NONM 4 4 4 0 -0.6 1.0 0.001 8 1 12
のようにする。このレベルでのファイルの中身はOsaka2kのマニュアルに説明されている。
pdosdrの出力ファイルの内、fort.13を利用してポストスクリプト図を描くのがaydosp
である。その制御ファイルaydosp.inpは 61 -0.6 1.0 150.0 120.0 10.0 50.0 0.376 1 1 0 DOS of GaAs などどする。一行目は出力ファイル番号である。次がエネルギースケールを与えるもので ある。 最終的なDOS図は図1に示されている。
3.3
tetrahedron
法 ー pdosdr の使い方
もともとpdosdrは全DOSを書くためだけのプログラムではない。部分DOSを描くこと
を目的としたものである。次にこの本来の機能を使うことを行ってみる。この際いくつかの パラメータの引き渡しに注意が必要となるので、まず全体の流れをつかんでもらうため計算 パラメータの受け渡しを図2で示す。 JobType dos ... nkdiv 8 ... OPTION BEGIN decompDOS ON OPTION END ======================== Decomposition of DOS parameters maxcmp : 8 ncmp0dos in pdoddr.inp contracted decomposition by nkat= 2 NO ik at l 1 1 ga 0 2 1 ga 1 5 2 as 0 6 2 as 1 The actual n of lm components calculated = 4 N of Bloch functions with G block 3 gaas. NONMSPIN-ORBIT 4 4 4 8 MAXCMP 4 NCOMP 1 Ga2s 2 Ga4p 5 As4s 6 As4p -0.6 1.0 0.001 8 1 12 61 -0.6 1.0 150.0 120.0 1.0 10.0 0.376 1 4 1 2 3 4 pDOS of GaAs renumbering
3.3.1 pwbcd 部分DOSの計算には、前節のpwbcdへの入力ファイルbcd.paraの最後にオプションを 指定する。 OPTION BEGIN decompDOS ON OPTION END これによりpwbcdはDOSを原子軌道に分解する。初期設定ではしない。それはレポート２ ８でも述べられてることであるが、原子軌道に分解するのにポテンシャルを構築するのと同 じくらいの計算時間がかかるからである。
ファイルdos gaas.outのDOS計算の部分の初めに、分解の条件が書かれている。
========================= Decomposition of DOS
parameters
maxcmp : 8
The array dimension for pdos
Write this number for ncmp0dos in pdoddr.inp contracted decomposition by nkat = 2
NO ik at l
1 1 ga 0
2 1 ga 1
5 2 as 0
6 2 as 1
The actual n of lm components calculated = 4 N of Bloch functions with G block 3
Total number of Bloch functions = 15
1 1 - 1
2 2 - 9
3 10 - 15
Type of atomic orbital: psuedo-wavefunction と、初めに分解した部分DOSの全種類数（maxcmp）が出力される。基本的には各原子ごと に（nkatごとに）、s、p、d、f と４種類に分解され出力される。したがって、デフォルト では maxcmp = nkat× 4 (9) である。ここの部分には、もともとpdosdrが全電子計算に付随するものとして書かれたこ とが反映する。擬ポテンシャルではどの原子もこのようなs、p、d、fの４種類の軌道を使っ ているわけではなく、多くはs、pの二つである。その場合はd、f軌道は形式的に出力する が、値は全て0である。そこで次の行以下に、実際に擬ポテンシャルとして部分DOSを計
算した軌道をリストしてある。つまりGa4s、Ga4p、As4s、As4pの４つの原子軌道に分解
される。その出力ファイルにおける番号が、1、2、5、6となる。分解される順番はここに並 べられた順である。そしてその順でファイルfort.12、fort.13に出力される。ラベルik はそれらの原子の番号である。 ファイルdos gaas.outでは、引き続き、原子軌道の質が記述される。まず原子軌道ϕalm の規格性が試される。つまりϕal0m0|ϕalm ® が計算され、それがどれくらいδl0l,m0mと異なる かが調べられる。
Normalization factor Normality of atom 1
1.36E-01 6.11E-20 5.90E-20 6.10E-20
2.33E-20 2.67E-02 1.34E-19 7.35E-19
2.82E-20 1.34E-19 2.67E-02 1.01E-18
2.69E-20 7.35E-19 1.01E-18 2.67E-02
Normality of atom 2
9.50E-02 1.61E-19 1.63E-19 1.56E-19
4.69E-20 5.17E-02 3.93E-19 1.34E-18
4.83E-19 3.93E-19 5.17E-02 2.11E-18
4.41E-20 1.34E-18 2.09E-18 5.17E-02
は、行列要素ϕal0m0|ϕalm ® をspxpypzの順で並べたものである。見て分かる通り、非対角成 分は０で原子軌道の直交性は非常に良い。しかし対角成分は１となっていない。 次に、GaとAs異種原子間でのspxpypz軌道の重なり積分 ϕal0m0|ϕblm ® が試される。 Orthgonality Combination atom 1 - 2
8.55E-02 9.31E-11 9.31E-11 9.31E-11
9.86E-11 7.66E-03 3.69E-10 3.69E-10
9.86E-11 3.69E-10 7.66E-03 3.69E-10
9.86E-11 3.69E-10 3.69E-10 7.66E-03
異なる原子に属する原子軌道は直交していないはずである。それにも関わらずこの結果は、 非対角成分は非常に良い精度で０となっていることを示している。どういうことだろうか？ 原子軌道の直交性 もともとの原子軌道は異なる原子同士では直交していないはずである。 それにも関わらず上で見たように非対角成分は０となっている。その理由は軌道の対称性に ある。異なった二つの原子に属するそれぞれa、b原子軌道の間の行列要素Eabは強結合近 似の記述ではSlater Koster表現で与えられる。例えば、a、b原子軌道としてs、p軌道であ れば Es,x= αVspσ (10) ここに(α, β, γ)は二つの原子の間の方向余弦である。原子位置を明示的に示すとEs,x(R1, R2) となる。もちろんこれが0になることは一般的には無い。 A B
図 3: (Left) a combination of atom A and B. (Right) atomic orbital representation in a crystal
結晶中では、原子位置RはR = Rl+ Rbと、単位胞の位置Rlとその中での基本構造原 子の相対位置Rbの和で表される。これを便宜上、(lb)、あるいは(lmn, b)と記そう。ここ にl = (lmn)は単位胞のインデックスである。pwbcdの中では原子座標として基本構造原子 Rbの全てを分解している。しかし単位胞の位置Rlは区別せず（これは当たり前である。そ うしないと無限個の組み合わせが必要となる）、そのlに関する和、いわゆるブロッホ和だ けを計算している。zincblendeの場合単位胞に２個の原子を含み、この原子軌道分解ではそ の２個の原子は区別している。しかしそのEs,x(1, 2)と言う場合、Es,x((000, 1), (000, 2))を 計算しているのではない。等価の他の単位胞から来ているものものも全て含む、すなわち Es,x= Es,x((0, 1), (000, 2)) + Es,x((0, 1), (¯1¯10, 2))+ Es,x((0, 1), (¯10¯1, 2)) + Es,x((0, 1), (0¯1¯1, 2)) (11) である。この和のため式（10）の方向余弦は打ち消しあって和は０となる。 このように単位胞の中で見ると、一見独立した原子軌道のように思えるものでも、並進対 称性により隣の単位胞の原子とも結ばれ、結果として原子軌道としての直交性がなくとも結 晶の対称性の理由により直交性が現われることがしばしばある。これを避けて孤立原子のと きの性格を引き出したかったら、大きなスーパーセルを用いて、隣の単位胞の原子との相互 作用を弱くするというアイデアがある。 3.3.2 pdosdr pdosdrの使い方は以下のようになる。入力ファイルpdosdr.inpは、 gaas. NONM 4 4 4 8 4 1 Ga4s 2 Ga4p 5 As4s 6 As4p -0.6 1.0 0.001 8 1 12 などとする。第3.3.1節のpwbcdの出力との対応を見るべきである。この意味は表1のよう にまとめられている。前にも述べたように、もともとpdosdrは全電子計算に付随するもの として書かれている。Osaka2kに移植するとき、入力インターフェースはできるだけ元の ものと同じになるように努めたが、それでもいくらか入力フォーマットを変更せざるを得な かった。まず扱う結晶の名前を最初に入れた。
部分DOSの成分数に関しては上の例のように、MAXCMPとNCOMPの二つを指定しなけれ
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Energy \Ry.\ 1 2 3 4 5 States/Ry.
pDOS of GaAs
Ga4s Ga4p As4s As4p図4: Partial DOS of GaAs.
分DOSへの分解のため取っているディメンジョンは式（9）のようになっている。つまり、
今の例GaAsであれば、Ga、Asそれぞれについてs、p、d、f の４種類、合計８種類であ
る。これがMAXCMPに指定すべき数である。しかし実際に擬ポテンシャルとして部分DOS
を計算するのは、Ga4s、Ga4p、As4s、As4pの４つだけである。それらをこの例のように
明示的に指定し、その数をNCOMPに入れる。
ここに、部分DOS（pDOS）を描くときの原子軌道のラベルを知る必要がある。通常それ
は*.primの中の原子の種類順に行われる。今の例では
Number of atom species 2
No Name Zat Zval
1 ga 31 3
2 as 33 5
Kind of atoms 2
Number of atoms 2
L.L. AND U.U. VALENCE ELEMENT
1 1 3.0000 1 ga 2 2 5.0000 2 as となっているので、Ga、Asの順となる。原子の種類というのはプログラム中ではnkatとい う変数で表されるところの既約原子の種類数を意味する。しかしこれだとスーパーセルでは 非常に大きな数となり、成分分解が大変になるのでnkatがある程度大きくなると2、nspec で表されるところの元素の種類について割り振られる。次に各の原子に関してpからf軌道 の順に並べられる。dos *.outの出力ファイルで確認すると、 2 nkmax0dosによって指定される。その初期値は6である。
表 1: input form of pdosdr.inp. When MAXCMP=0, specification NCOMP is not re-quired.
description
NAME name of the crystal with a period
MAGNET specify the magnetic state
(NONM, MAGN, SPIN)
NX, NY, NZ segments of BZ
MAXCMP dimension of decomposed orbitals
(NCOMP number of selected orbitals
followed by each name of selected orbitals) ESTART, ENEND, DE energy scale
begin at E =ESTART end at E =ENEND with spacing DE
NELEC, NB1, NB2 number of electrons
spanning from band NB1 to NB2
Number of atom species Decomposition of DOS
parameters
maxcmp : 8
contracted decomposition by nkat = 2
NO ik at l 1 1 ga 0 2 1 ga 1 5 2 as 0 6 2 as 1 のように原子軌道に通し番号が割り振られる。pdosdr.inpで指定する番号はこの通し番号 である。
pdosdrの出力ファイルfort.13はDOSをエネルギー順に書き下している。はじめに全
DOSを、次にpdosdr.inpで指定されたpDOSをその順にリストする。スピンがあるとき
はスピンアップに対しての全DOSそして要求があればそれのpDOSをリストする。次にス
ピンダウンに関してその手順を繰り返す。
3.3.3 aydosp
aydospはpdosdrの出力（特にfort.13）を受け、それをPSファイルに図示化する。そ
の入力ファイルaydosp.inpは
61
4 1 2 3 4 DOS of GaAs のようにする。
表 2: input form of aydosp.inp description
IFIL specify file number
EO, EM, minimum and maximum of energy
XM, YM, scale of x- and y- axes (mm) DDOS, DOSM, the tick and full scale of DOS EF, ISPIN Fermi level and spin status NLINE, JLINE(NLINE) number of pDOS lines
followed by each pDOS specification
D title
ここでDOS成分の番号付けについて注意が必要である。今度の番号はpdosdrの出力ファ
イルfort.13に現れる番号である。つまりpdosdr.inpで指定される番号と違うので注意
が必要だ。出力ファイルfort.13ではいつでも最初のDOSは全DOSである。これを0番
目と数える。
NLINE, JLINE(NLINE)の行を例示すると次のようになる。
no SPIN
1 0 -> 1 data; total DOS only
2 0 2 -> 2 data; total DOS and a partial DOS of #2
3 1 3 4 -> 3 data: 3 partial DOSs of #1, 3, and 4
SPIN
2 0 1 -> 2 data: up and down total DOSs
最終的なDOS図は図4に示されている。
3.4
原子軌道展開の取り方
理論のところで述べたが、原子軌道のブロッホ和 ϕk,alm(r) = X R exp[ik· (R + Ra)] ϕκlm(r− R − Ra) (12) をもって結晶中における原子軌道と見なしている。平面波展開ではkは第一ゾーンに限定し てもなんら問題はなかったが、式（12）の中では問題があることを指摘した。実際、図4を 見ると、価電子帯では問題ないが、伝導帯をみるとその部分DOS成分がほとんど無くなっ ている。励起状態が式（12）で展開されるのものには含まれていないからである。この問題-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Energy \Ry.\ 1 2 3 4 5 6 7 8 9 10 States/Ry.
pDOS of GaAs
Ga4s Ga4p As4s As4p図5: Partial DOS of GaAs.
はもっと高い主量子数の原子波動関数ψnlを使うことで解決できるが、擬ポテンシャル法で はそのような波動関数は用意されていない。そこで式（12）のkを第一ゾーンの外まで拡大 して励起状態の波動関数を表現することにした。kの大きなものも含めたときブロッホ和が どのようになるのかを示したものが図6である。 実際のところはk∈ 1st BZに対し、k + Gとしてどこまでとるかが問題となる。必要な Gの組は大きさGの三乗に比例して増加し、（対称性を利用し計算量は減少できるとはい え）それだけ原子軌道分解の組{Ilκ(q)}を用意しなければならず一気に計算量が増加する。 計算に入れるk + Gの組はk + Gの大きさ毎にグループ分けしておき（シェルと呼ぼう）、 どこまでのシェルを取るかを指定する。デフォルトでは第１シェルまで取っている。 リスト１で
The actual n of lm components calculated = 4 N of Bloch functions with G block 3
Total number of Bloch functions = 15
1 1 - 1 2 2 - 9 3 10 - 15 となっているのはそのことである。すなわちk + Gの組は第３シェルまで取る。第一シェ ルの中には１つだけ、第二シェルには８個、第３シェルには６個あることを示している。そ の結果が図5に示されている。図4と見比べてみると、今回のものは明らかに伝導帯にも十 分な部分DOS成分が現われていることがわかる。 しかしその代償として、価電子帯にも上のk + G成分が入り込み、かえってs、p成分の 違いが分からなくなってしまった。これでは元も子もない。
G
k X
図6: Bloch sum of s orbital for excited states and s band in the band structure.
3.5
pDOS
の表現の仕方
ところでpDOSの表現の仕方にはいろいろ任意性がある。例えば式（7）で求められる hψk|ϕκlmiはもちろん複素数なので、そのままでは実数座標にはプロットできない。どうする か。まず考えられるのが、その絶対値自乗|hψk|ϕκlmi| 2を取ることである。これまでのもの は（例えば図4）、実際このように取っている。さらにmを分解しないときは、 X m |hψk|ϕκlmi| 2 (13) をプロットしていることになる。しかしこれを X m |hψk|ϕκlmi| (14) と自乗しない形に選ぶこともできる。解釈の違いだけである。 図7にそれらを比較してある。定義の違いだけでどちらでも本質的に同じである。以降 は、定義として一番自然な自乗和を取っている。4
pdosdr
での
k
点メッシュの取り方
pdosdrはDOSを描くのにTetrahedron法という優れたアルゴリズムを採用しているの
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Energy \Ry.\ 1 2 3 4 5 States/Ry.
pDOS of GaAs sum of abs^2
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Energy \Ry.\ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 States/Ry.
pDOS of GaAs sqrt sum
図7: Comparison of two ways of pdos. The left is the sum of square of the matrix element, while the right is the sum of absolute value of the elements.
と違う点があるため、いろいろトラブルの原因となっている。そこで今回問題をユーザーの 前に提示し、どのように使うべきかを示した。
まず、先に推奨されるpwbcdとpdosdrのk点メッシュの取り方を表4にまとめて示す。
ユーザーは迷ったときはいつもこれに従って入力すべきだ。
表3: Appropriate mesh points between pwbcd and pdosdr
crystal ndiv notes
system pwbcd pdosdr on pdosdr
P n n n should be even F 2n n I 2n n C 2n (n, n, 2n) R 3n n Ｃに関しては、monoclinicとorthorhombicとがあるが、全てはまだチェックしていない。 実際のところ、pdosdrの方では、結晶の型でk点メッシュの分割数にさらに細かな制限が ある。場合の数が多いので、pwbcdとpdosdrとの整合性に関しては、まだ著者も把握しき れていない。k点メッシュで問題が生じたとき、pdosdrの標準出力ファイルで、その整合 性をチェックしてみることだ。
例えばrhombohedral systemに例で示す。pdosdrの標準出力ファイルは入力ファイル
fort.2でのk点情報
1 0 0 0 6 1 GM 1 1 2 21 1 0 0 0 6 1 GM 2 1 2 7 2 0 0 -3 6 1 LD 1 1 4 37 2 0 0 -3 6 1 LD 2 1 4 12 2 0 0 -3 6 1 LD 3 2 4 37 3 0 1 -1 6 1 YP 1 1 6 37 3 0 1 -1 6 1 YP 2 1 6 37 4 1 0 -2 6 1 XP 1 1 6 37 4 1 0 -2 6 1 XP 2 1 6 37 5 0 1 -4 6 1 YP 1 1 6 37 5 0 1 -4 6 1 YP 2 1 6 37 ... 31 0 3 6 6 1 F 4 1 3 37 32 2 1 -8 6 1 GN 1 1 12 37 NBAND= 37 1 37 を打ち出す。pwbcdでは6×6×6でメッシュを切っている。pdosdrの出力は、その後pdosdr 自身のk点メッシュ情報を示す。これは2× 2 × 2で切ったものである。 GENERATED KPOINT NO KX, KY, KZ 1 0 0 0 2 1 2 0 0 1 2 2 3 0 1 0 2 20 4 0 1 1 2 24 5 0 2 0 2 10 6 0 2 1 2 21 7 1 0 0 2 20 8 1 0 1 2 25 9 1 1 0 2 31 10 1 1 1 2 24 11 1 2 0 2 20 12 1 2 1 2 24 13 2 0 0 2 10 14 2 0 1 2 2 15 2 1 0 2 20 16 2 1 1 2 25 17 2 2 0 2 1 18 2 2 1 2 2 となっている。この点のセットが上の入力セットに全て含まれていなければならない。
5
まとめ
pwmで原子軌道分解した後、pdosdrにデータを引き渡す際、元のデーターにどのような 原子軌道分解したデータが含まれているかを常に意識しなければならない。それが、パラ メータMAXCMPとNCOMPの区別である。それは、原子軌道をどこまで（つまりlまでかmま でか）分解するか、あるいは問題とする原子のみを取るかで異なる。 関係するbcd.paraのオプションパラメータは以下のようになる。表4: Options for bcd.para Option Default Description toggle options
decompDOS OFF decomposition to partial DOS mdecomp0dos OFF decomposition to the m level value options
targetatom0dos none atom specification to pDOS nGblck0input 1 n of blocks G in Bloch sums
cutoff0KC modify NHDIM with this factor
A
label of m components
Osaka2kで使われる角運動量波動関数の定義を述べる。部分DOSのm成分はこの表に現 れる順で出てくるので、部分DOSの解釈に当たっては、この表を参照する必要がある。B
colors in aydosp
aydospの中で使われている色と番号の対応を記しておく。 c 1:Black c 2:Red c 3:Green c 4:Blue c 5:Orange c 6:lavender c 7:pink c 8:lime c 9:yellow c 10:bright Blue参考文献
[1] Technical Report No. 28「波動関数の原子軌道への成分表示」
[2] Technical Report No. 61「DOSの原子軌道への分解」
[3] Technical Report No. 83「Osaka2kにおけるpdosdrの使い方 ーk点メッシュの切 り方ー」
表5: definition of real-valued spherical harmonics in Osaka2k. The normalization factor is ignored. l m Ylm 0 0 1 1 -1 x 0 y 1 z 2 -2 xy -1 zx 0 z2− x2− y2 1 yz 2 y2− x2 3 -3 x(3y2− x2) -2 xyz -1 x(4z2− x2− y2) 0 z(2z2− 3(x2+ y2)) 1 y(4z2− x2− y2) 2 z(y2− x2) 3 y(y2− 3x2)
