(Creation) 2017 年 8 月 31 日

ABCAP を用いた光学定数の計算
浜田典昭
東京理科大学理工学部物理学科〒278-8510 野田市山崎 2641

全電子状態計算プログラム「All electron Band structure CAlculaion Package (ABCAP) 」abcap1707の使用法について説明する。プログラムについてのやや詳しい解説は「計算機マテリアルデザイン入門」[1] を見ていただきたい。ここでは、ZnO を例にバンド計算の手順を簡単に述べた後、バンド間遷移による光学定数の計算方法と計算結果を紹介する。バンド計算で得られた波動関数とバンド構造を用いて誘電関数を計算し、複素屈折率、吸収係数などを評価する。ドーピングはリジッドバンドの範囲で取り扱う。

1 準備
1.1 ディレクトリ構成
プログラムは、ディレクトリ abcap1707 以下に収められ、つぎのようなディレクトリがある。

src/

ソースファイル

data/

データファイル

samples/

計算例

bin/

実行ファイルの格納

lib/

オブジェクトモジュール (.a)
モジュールファイル (.mod)

include/インクルードファイル (.h)

これを使うにはいくつかの環境変数の設定が必要である。

1.2 環境変数の設定例
また、コンパイラの名前などを次の例に習ってそれぞれの環境変数に入れる。
カレントディレクトリ (.) を path に入れ、環境変数 ABCAP にプログラムパッケージのあFC_TYPE=’gen’るディレクトリをセットする。
(bash の場合)

FC=’ifort’
FO=

PATH=.:$PATH

MPIFC=’mpif90’

ABCAP=/home/band/abcap1707

MPIFO=

環境変数 ABCAP には、実際にパッケージがある
ここで FO とあるのは、fortran コンパイラオプ

絶対ディレクトリを入れる。

ションであるが、特に値を入れなくてもとりあえず走ると期待される。
常にこの環境が設定される必要があるので、.bashrc か .bash_profile が呼ばれたときにこの環境が設定されるようにする。
環境設定は、ディレクトリ abcap1707 直下のファイル Envs_sh (または Envs_csh) を参考にしていただきたい。

11.3 コンパイルとリンク
ディレクトリ abcap1707 直下で make を実行すると、必要なプログラムはほとんどすべてコンパイル・リンクされ、実行モジュールがディレクトリ bin の下に格納される。

2 バンド計算の手順
2.1 計算用ディレクトリの設定
次のようにすると便利である。
まず、計算用ディレクトリを作る：

cd

$home

mkdir abc (任意の名前)
cd abccp $abcap1707/samples/Setnew.sh .
次に、一つの計算物質 (ZnO) に付き一つのディレクトリを作る。

cd {$home}/abc
mkdir ZnO
cdZnO../Setnew.sh
コマンド Setnew.sh が、例題のディレクトリ $ABCAP/samples/LaMnO3c_f_6/ および data ディレクト
リ から、計算に必要なファイル (.data, .sh) を計算物質ディレクトリにコピーする。入力ファイル (.data)を修正して目的の物質を計算するが、修正せずに実行すると例題「強磁性立方晶 LaMnO3 」が走る。

2.2 実行順序
計算用ディレクトリの下でHとコマンドを打つと、以下のような計算手順が画面に出る。

-----------------------------------------------------------------(0) (ab_prp.data, atom.data)

ab_prp.sh(1) (ab_input.data)

ab_in.sh-----------------------------------------------------------------(2) Hu (for +U calculation)-----------------------------------------------------------------(3) (ab_input.data)

fl06.sh or fl06b.shcheck.sh-----------------------------------------------------------------(4a) (bn_atps.data)bn_atps.sh(p3_atps.data)p3_atps.sh------------------------------------------ghostview plot.ps -----(4b) <noSOi>(bnpl.data)

bnpl.sh<noSOi>(angle.data)brzone.sh --> brzone.ps2

<SOi>

cd bndraw_so(bndraw.data)

bndraw.sh------------------------------------------ghostview plot?.ps -----(4c) (bndraw.data)

bndraw.sh------------------------------------------ghostview plot.ps -----(4d) (bn_pdos.data)

bn_pdos.sh(p2_dos.data)

p2_dos.sh------------------------------------------ghostview plot?.ps -----番号の順に実行する。ファイル名.data の内容を編集し、ファイル名.sh を実行する。a, b, c, ... で区別されたものの順番は任意である。なお、コマンドスクリプトファイルには拡張子.sh が付いている。
以下、ウルツ鉱型の結晶構造を持つ ZnO の例に沿って、それぞれのジョブを説明する。

2.3 ab prp.sh
入力ファイル ab_prp.data----------------------------------------------------------abcap-ab_prp.data
0!jpr

ZnO (wurtzite)
lattice parameter -2----*----3----*----4----*----5----*----6----*----7
3.249 3.249 5.207 90.0 90.0 120.0 !a,b,c[A], alpha,beta,gamma[degree]
space group
3-2----*----3----*----4----*----5----*----6----*----70202

0 1

0 1

1 2!igen,jgen(2,3)

19

0 1

0 1

1 2

!igen,jgen(2,3)

kinds of atoms!idim, il(r,h,s,f,b,oz,ox,oy),ngen,inv-2----*----3----*----4----*----5----*----6----*----7

2!# of kinds

1 0.6666666666666667 0.3333333333333333 0.0

Zn!jpos,position,name

1 0.6666666666666667 0.3333333333333333 0.345 O!jpos,position,namemagnetic state-2----*----3----*----4----*----5----*----6----*----7

0!jmag(0,1,2)

totally symmetric basis set -3----*----4----*----5----*----6----*----7
24.0

6!cut-off energy[Hr],Lmax

k-points (# of division)
6

6

4!nx,ny,nz

iteration
4

60.10---3----*----4----*----5----*----6----*----7-2----*----3----*----4----*----5----*----6----*----7
0.10!method, n-method, pmix, amix!---*----1----*----2----*----3----*----4----*----5----*----6----*----7----------------------------------------------------------• １行目はファイルのヘッダー（abcap-に続いてファイル名）• ４行目以降の、lattice, space, kinds, magnetic, totally, k-points などのキーワードはデータの属性を表し、最初の４文字がプログラムにより読まれ、判断される。これらのキーワードに続く一固まりのデータ（データブロック）はこのような形式に書かれていなければならない。データブロッ3クは互いに順番を換えても良い。
デ ー タ フ ァ イ ル と し て atom.data が 必 要 で あ る が 、こ れ は デ ー タ ベ ー ス と し て 与 え ら れ て い る

data/atom0.data から初めにコピーされている。たとえば、Zn に対しては、------------------------------Zn
30.0865.31.8

100 200 210 300 310 324 405 415
1.0 1.0 3.0 1.0 3.0 5.0 1.0 0.0
1.0 1.0 3.0 1.0 3.0 5.0 1.0 0.0------------------------------という内容を持っている。2 行目最後の数字が MT 球の半径 [Bohr] である。3 行目の三桁の数字を nlp と書くと、n が主量子数、l が角運動量、p が動径波動関数設定条件パラメータである。p = 0 はそれをコアとして扱うことを意味する。p ̸= 0 の shell は価電子としてバンド計算される。最後の 2 行が up-spin とdown-spin の占有数である。
ab prp.sh を実行すると、以下の計算で使う ab input.data ファイルが作られる。標準的なことをするのであれば、データファイル ab input.data は編集しなくてよい。

2.4 ab in.sh
ab in.sh では、次のプログラムが走る。
ab_in

全対称基底関数を作り、a_spw.dta と c_ssw.dta に書く。

:

ab_inch:原子の電子状態を計算し、その重ね合わせとして初期電子密度を作る。

ab_kpgn:

計算すべきｋ点を生成し、ファイル a_kp0.dta, と a_kp2.dta に書く。

ab_size:

計算に必要なサイズを計算し、ファイル f_size0.dta, f_size1.dta, para.inc, para1.inc. に書く。

MT 球間の全対称基底関数の重なり積分を計算し、f_ospw.dtb に書く。

ab_ospw:

ここで、.dta はアスキーファイル、.dtb はバイナリーファイルを表す。

2.5 fl06.sh
fl06.sh を用いて、自己無道着計算のための繰り返し計算を実行する。
fl06.sh繰り返しの回数はこのファイルの中の変数 ITER_MAIN に指定する。たとえば、set ITER_MAIN = 16 で 16 回の繰り返し計算をする。この中では次のプログラムが走る。

fl_pot :密度汎関数理論に基づきポテンシャルを計算する。

fl_bnd :内殻電子および価電子の固有状態を計算する。

fl_chg :電子密度を計算する。

fl_ptuj:+U 計算をするとき、軌道に依存するポテンシャルを計算する。

fl_pot :全エネルギーを計算する。

fl_mx5:

入力と出力の電子密度の差を計算し、次の入力電子密度を決める。

4

check.sh コマンド
fl06.sh を実行した後、check.sh を実行すると、繰り返し計算の最後の数回の結果が画面に表示される。画面に現れたものの中で、whole cell とあるのが、繰り返し一回毎の入出力の電子密度分布の差 (の２乗平均のルート) [electrons/a.u.3 ] であり、10−5 程度になるとかなり収束しており、10−7 程度になると全エネ
ルギーの 10−6 [Hr] の桁が動かなくなる。必要な物理量の収束が達成されるまで、fl06.sh を繰り返す。

f ef.dta ファイル：--------------------------------------abcap-ef[Hr]:Fermi-level, dos, vale, band-E (/spin)

0.1265509385972843E+00
18

18

134

0.0000000000E+00

34# of fully-occ. bands, # of occ. bands

0.1800000000E+02 -0.1840833730E+01 spin=1--------------------------------------2 行目はフェルミ準位である。3 行目は、18 番目までのバンドが full に詰まっていることを示している。

2.6 結晶構造の描画
bn atps.sh を実行することにより中間ファイル p3 atps.data と p3 atps.dta が得られ、これらに基づいて p3 atps.sh を実行することにより結晶構造の図がファイル plot.ps にポストスクリプト言語で書かれる。ghostscript(gs) コマンド や ポストスクリプトプリンタにより描画することができる。

2.7 バンド構造の描画
bnpl.sh を実行することにより、バンド構造（E-k カーブ）が得られる。出力はポストスクリプトファイ
ル plot.ps である。（磁性体の場合には plot1.ps と plot2.ps に出力される。）

2.8 状態密度の描画
bn pdos.sh に続いて p2 dos.sh を実行する。出力はポストスクリプトファイル plot.ps である。（磁性体の場合には plot1.ps と plot2.ps に出力される。）
描画例を図 1 に示す。

3 光学定数の計算
FLAPW バンド計算によって得られた波動関数 ψkn を用いて誘電関数 ϵij (ω) を計算する。(ここで、i, j = x, y, z である。) バンド間遷移の寄与として次式を計算する（プログラム df01）。
[]e2
1ϵij (ω)/ϵ0 = δij +χ(ω)−nijϵ0 mω 2 m∑[]cv∗χij (ω) = 2
χcv
ij (ω + iη) + χij (−ω + iη)
cv

1
χcv
ij (z) = −N⟨⟩⟨⟩
∑ fkv (1 − fkc ) v pi c c pj v
z − (εkc − εkv )k5(1)(2)(a) 結晶構造(b) バンド構造(c) 状態密度

図 1 ZnO金属の場合には、ドルーデの寄与として次式を計算し、誘電関数に加える。

i σij (0)−1 σij (0)/τ=ω 1 − iωτω ω + i/τ
2eσij (0)/τ =⟨vi vj ⟩F D(εF )ϵ0電子の速度 v は、FLAPW バンド計算によって得られた波動関数 ψkn を用いて⟨⟩pv kn = ψkn
ψkn ,m(3)により計算する。p は運動量、m は電子質量であり、k は結晶波数、n はバンドの番号である。フェルミ面上における速度テンソルは∑ Ωc ∫dSF⟨vi vj ⟩F D(εF ) =v kn v kn
3(2π)
ℏ|v
kn |
n(4)で計算される。フェルミ面における面積分はテトラヘドロン法を用いて比較的簡単に実行できる。ドルーデモデルにおける緩和時間 τ は簡単には評価できない。計算プログラムでは、ℏを入力として与える。
τ3.1 光学定数の計算手順
バンド計算を実行したディレクトリ ZnO の一つ上のディレクトリ（ls コマンドで ZnO が表示されるディレクトリ）で、$ABCAP/samples/Setopt.sh ZnOとコマンドを打つと、新しいディレクトリ ZnO s が作成され、./ZnO/ と$ABCAP/samples/ZnO_5_s0/から光学定数の計算に必要なファイルが ZnO s の下にコピーされる。この ZnO s の下に移動し、Hとコマンドを打つと、以下の計算手順が画面に出る。

6

------------------------------------------------(0) (bn_scis.data)

bn_scis.sh(do this just once)

bngap.sh!optical gap------------------------------------------------(1) (bn_dope.data)

bn_dope.sh!bands at the fermi energy(bnpl.data)

bnpl.sh!E-k curve(2) (df01.data)

df01.sh!dielectric function (interband epsilon(omega))(3) (wpw.data)

wpw.sh!<w|p|w>bn_efp2.sh!Fermi velocity <v_i*v_j>(df_total.data) df_total.sh !dielectric function (+Drude term)------------------------------------------------(4) (p2_dielec.data)

p2_dielec.sh!total dielectric function (epsilon)(p2_optnk.data)

p2_optnk.sh!n, k(p2_absorb.data)

p2_absorb.sh!abcorption coefficient (alpha)optra.sh!T, R, A : no ps-filep2_optra.sh!T, R, A

(5) (optra.data)(p2_optra.data)------------------------------------------------- (plot.ps) ----ここで、(0) においてバンドギャップを調節 (scissors operation) する。(1) においてリジッドバンドを仮定して、ドーピングを行い、(2) においてバンド間遷移による誘電関数を計算し、(3) においてフェルミ面上の電子の速度を求め、ドルーデ項を計算して、誘電関数に加える。(4) において誘電関数、複素屈折率、吸収係数の描画を行う。(5) において、薄膜における反射率、吸収率、透過率を評価し、それらを描画する。

3.1.1

手順 (0)密度汎関数理論によるバンド計算は現在のところバンドギャップを正しく計算できないことが多いので、バンドギャップを実験に合わせて補正する必要がある。そのための入力ファイルは

bn_scis.data------------------------------------------------------19

2.3

i-conduction-band, scisors[eV] (for upspin)

19

2.3

i-conduction-band, scisors[eV] (for down spin)-------------------------------------------------------である。LDA による ZnO のバンドギャップの計算値はおよそ0.95eV である。実験値は温度に依存するがこの値より約 2.3eV大きい。上記の入力は伝導電子バンド（19 番目のバンド以上）を 2.3eV だけ平行移動することを意味する。スピン毎に行えるが非磁性の場合には上段のデータを用いる。

fundamental gap と optical gap がどのように修正されたか見るために、bngap.sh を走らせる。その結果はファイル bn-gap.log に書かれる。また、Γ点のバンドのエネルギー値をフェルミ準位を基準に eV 単位でファイル bngap.txt に書き出す。
バンドギャップ修正後のバンド構造

7

3.1.2

手順 (1)

bn dope.sh はリジッドバンドを仮定してドーピングを行う。
bn_dope.data ファイル-------------------------------------------------------00.001

jdope(0,1), dope(0:/cell, 1:/cm**3)-------------------------------------------------------または-------------------------------------------------------1

1.e19

jdope(0,1), dope(0:/cell, 1:/cm**3)-------------------------------------------------------加える電子数を、基本単位胞あたり、または、立方 cm あたりで与える。正の数が電子ドープ、負の数がホールドープとなる。（ドープしない場合は 0.0 ）フェルミ準位の位置を確認するために、bnpl.sh でバンド構造を描く。

3.1.3

手順 (2)

バンド間遷移による誘電関数を計算する。

df01.data ファイル------------------------------------------------------0jpr

200

ne0.010.0# meshemin,emax [eV] (photon energy)

0.000005 1 0

eshift[Hr],ipath,matopt

318

nv1,nv2: valence bands

1921

nw1,nw2: conduction bands------------------------------------------------------２行目と３行目で (光の) エネルギーメッシュを指定する。４行目の最初の数は計算に用いるエネルギーの虚数部分である。５行目は詰まっているバンド、６行目は空いているバンドを指定する。（bngap.txt が参考になる。）計算量は、詰まっているバンドの数と空いているバンドの数の積に比例する。
磁性体の場合には、最後の２行はスピン毎に指定する。

df01.data ファイルの最後の２行（強磁性金属の例）------------------------------------------------------3 19

3 16

nv1,nv2(up), nv1,nv2(down) : valence bands

18 26

15 25

nw1,nw2(up), nw1,nw2(down) : conduction bands-------------------------------------------------------

3.1.4

手順 (3)

wpw.sh によりバンドの各点 (kn) における速度の計算をする。そのための入力ファイルがwpw.data ファイル8

-------------------------------------------------------33

bnbgn+, bnend+(spin up)-33

bnbgn+, bnend+(spin down)------------------------------------------------------であり、フェルミ準位より下３番目のバンドから上３番目のバンドまで速度を計算することを示す。フェルミ準位を横切るバンドについては必ず計算するので、ドルーデ項を評価するには００とすれば十分であるが、それほど時間のかかる計算ではないので、上記のデータで問題はない。

bn efp2.sh によりフェルミ速度テンソルを計算する。
これを用いて、ドルーデ項も含めた誘電関数を計算する。

df_total.data ファイル------------------------------------------------------1

model0.02

damp[eV]------------------------------------------------------ドルーデモデルの緩和時間の逆数を eV 単位で与える。

3.1.5

手順 (4)

テンソル成分に次の順で番号を付け (1 ∼ 12)、番号を指定することで、テンソル成分またはそれらの平均値を描画をする。

------------------------------------------------------Re(xx), Im(xx),Re(xy), Im(xy),Re(xz), Im(xz)Re(yy), Im(yy),Re(yz), Im(yz)Re(zz), Im(zz)------------------------------------------------------誘電関数 の描画：

p2_dielec.data ファイル------------------------------------------------------00jpr, kpaper

2iconv(2:Hr->eV)

0.0

10.0

1.0

emin, emax, de-1.0

5.0

2.0

dmin, dmax, dd(integer)

8.0

5.0

ZnO

xe, yd(eV)(mm/u, mm/u)

Dielectric function

6ncurve(Imag. part) xx,yy,zz
3

2 8 12nm(i),(im(j,i),j=1,nm(i))(Real part) xx,yy,zz
3

1 7 11nm(i),(im(j,i),j=1,nm(i))(Imag. part) zz

9112(Realnm(i),(im(j,i),j=1,nm(i))part) zz1

11nm(i),(im(j,i),j=1,nm(i))(Imag. part) xx
12(Realnm(i),(im(j,i),j=1,nm(i))part) xx1

1nm(i),(im(j,i),j=1,nm(i))------------------------------------------------------光学定数 の描画：

p2_optnk.data ファイル------------------------------------------------------00jpr, kpaper2

iconv(2:Hr->eV)

0.0

6.0

1.0

emin, emax, de0.0

2.5

1.0

dmin, dmax, dd(integer)

12.0

30.0xe, yd(eV)(mm/u, mm/u)

ZnO
3

ncurveaverage0ndata, jdata(ndata)

zz
1

11

ndata, jdata(ndata)(xx+yy)/2
2

1

7

ndata, jdata(ndata)------------------------------------------------------吸収係数 の描画：

p2_absorb.data ファイル------------------------------------------------------00jpr, kpaper2

iconv(2:Hr->eV)

0.0

14.0

1.0

emin, emax, de0.0

1.5

0.5

dmin, dmax, dd

8.0

50.0

xe, yd(mm/u, mm/u)

ZnO
3ncurveaverage0ndata, jdata(ndata)zz
111ndata, jdata(ndata)xx+yy10(eV)

2

17ndata, jdata(ndata)------------------------------------------------吸収係数は 106 /cm を単位として表示される。

3.1.6

手順 (5)

薄膜に光を垂直入射したときの反射率、吸収率、透過率、および光学定数（屈折率と消衰係数）をファイル出力する：

optra.data ファイル------------------------------------------------------6

ij_epsilon(1:xx,4:yy,6:zz)

200.0

ddd [nm]

1.0d0

1.53d0n-incident, n-substrate-----------------------------------------------------１行目は入射光の電場の向き、２行目に薄膜の厚さ [nm]、３行目は入射側の物質の屈折率と出射側の物質の屈折率である。

薄膜の反射率、吸収率、透過率 の描画：

p2_optra.data ファイル------------------------------------------------------1

0jpr, kpaper(0,1)

1iscale(0,1)

ZnO (s)n$de=10$u2$u1 [cm$u-$u3],
Photon energy

00.0 10.0

1.0’F’

T, R, A[%]
20.000.0

100.0

10.60.1d=200[nm],n$dg=1.53[eV]
iautox,xmin,xmax,dx,FEIx(F,E,I)’F’ymin,ymax,dy,FEIy(F,E,I)
isize,ux,uy (cm/unit)

1 3

ifmt(0:sut format,1:gnu format),# of lines’f_TRAeV.dta’,1 2 1’Transmittance’

iix,iiy,type of line

’f_TRAeV.dta’,1 3 2’Reflectance’

iix,iiy,type of line

’f_TRAeV.dta’,1 4 3’Absorptance’

iix,iiy,type of line-------------------------------------------------------

3.2 計算結果
図 2 に ZnO の計算結果を示した。

11

(a) 誘電関数(b) 屈折率と消衰係数
図 2 ZnO

参考文献[1] 浜田典昭, 「笠井・赤井・吉田編：計算機マテリアルデザイン入門 (大阪大学出版会, 2005)」第７章（基礎編および応用編）.

12