自由電子 (空格子) バンド
Ψ ( x) = C exp[i (k + G )] = C exp[i (k + ha )] h = ･･･, -2, -1, 0, 1,2, ･･･
k*h

拡張ゾーン表示

1st BZB

2nd BZB

3

還元ゾーン表示
3BZの併進対称性

22kabs = 1.2

22E (k ) =
kabs
2m
1

kabs = 1

1

Blochの波数: kBloch
実際の波数 : kabs

BZの併進対称性
kabs = 1/2-2.0 -1.5 -1.0 -0.5

00.00.5k (in 2π/a)

1.0

1.5

2.0-0.500kBloch (in 2π/a)

0.5

プログラム: 自由電子バンド
free_electron_band.py
Si の Bravais格子のk点に対してプロット
a = 5.4064 Å
m* = 1.0me
W: (1/2 0 1)
L: (1/2 1/2 1/2)
Γ: (0 0 0)
X: (0 0 1)
W: (1/2 0 1)
K: (3/4 0 3/4)

プログラム: 自由電子バンド

free_electron_band.py

from pprint import pprint #リストを整形して出力
a = 5.4064# angstrom, lattice parameter#逆格子のmetricsを計算
rg = np.zeros([3, 3])
rg[0][0] = 2.0 * pi / arg[1][1] = rg[0][0]rg[2][2] = rg[0][0]# バンド構造をプロットするk点の軌跡:# [kx, ky, kz, k点名称]
klist = [[0.5, 0.0, 1.0, "W"]
, [0.5, 0.5, 0.5, "L"], [0.0, 0.0, 0.0, "$¥Gamma$"], [0.0, 0.0, 1.0, "X"], [0.5, 0.0, 1.0, "W"]
, [0.75, 0.0, 0.75, "K"]]# プロットするバンド構造E(k)のk点数の概数
nk = 101# Ehkl(k)を計算するhkl範囲
hrange = [-3, 3]
krange = [-3, 3]
lrange = [-3, 3]# プロットするエネルギー範囲
Erange = [0.0, 10.0] # eV
# 逆格子のmetrixから、2点のk点間の距離を計算
def cal_kdistance(rg, k0, k1):
dkx = k1[0] - k0[0]dky = k1[1] - k0[1]dkz = k1[2] - k0[2]r2 = rg[0][0] * dkx*dkx + rg[1][1] * dky*dky +rg[2][2] * dkz*dkz
r2 += 2.0 * (rg[0][1] * dkx*dky + rg[1][2] *dky*dkz + rg[2][0] * dky*dkx)
return sqrt(r2)プログラム: 自由電子バンド

free_electron_band.py# k点を与えて自由電子のエネルギーを計算
def cal_E(k, Ghkl):
global rg# k点のリストとhkl範囲を与え、Ehkl(k)を計算
def get_cal_Elist(xkvec, hrange, krange, lrange):yE = []
for i in range(len(xkvec)):
kabs2 = rg[0][0] * (k[0] + Ghkl[0])**2
kx = xkvec[i][0]kabs2 += rg[1][1] * (k[1] + Ghkl[1])**2
ky = xkvec[i][1]kabs2 += rg[2][2] * (k[2] + Ghkl[2])**2
kz = xkvec[i][2]Elist = []return KE * kabs2 # in eV
for ih in range(hrange[0], hrange[1]+1):
for ik in range(krange[0], krange[1]+1):# プロットするk点リスト klistとk点数の概数 nk から、
for il in range(lrange[0], lrange[1]+1):# なるべくk点間隔が等間隔になるように、
E = cal_E([kx, ky, kz], [ih, ik, il])
# 計算するk点などをリストアップする
Elist.append(E)# バンド構造プロットに必要なリストも返す
def get_cal_klist(klist, nk):yE.append(Elist)…
return yE

プログラム: 自由電子バンド

free_electron_band.py# バンド構造をプロット
def plot_band(axis, xk, yE, Erange, ktotallist,ktotal_namelist):# 表示範囲は決め打ち
axis.set_xlim([min(xk), max(xk)])
axis.set_ylim(Erange)# バンド構造をプロット
axis.plot(xk, yE, linestyle = 'none',marker = 'o', markerfacecolor = 'none',markeredgecolor = 'black',
markeredgewidth = 0.5, markersize = 2.0)
# Γ点、BZ境界の縦線を引く
for i in range(1, len(ktotallist)):axis.plot([ktotallist[i-1], ktotallist[i-1]], Erange,linestyle = '-', color = 'black', linewidth =0.5)# k軸の目盛りにk点の名称を表示する# グラフ枠が一つであれば plt.xtics()で設定できる# axisに対しては、.setpでattributeを直接書き換える必要があるらしい
plt.setp(axis, xticks = ktotallist, xticklabels =ktotal_namelist)
axis.set_xlabel("k", fontsize = fontsize)
axis.set_ylabel("E (eV)", fontsize = fontsize)
axis.tick_params(labelsize = fontsize)(拡がった)バンド構造は自由電子として理解できる: Siの例
VASP, PBE96, Conventional cell
10Energy / eV

5G hkl = ha* + kb* + lc*

0-5-10W

自由電子モデル2 2 2(k + G hkl )2E free (k ) =K =
2m2mL

Γ

X WK

Na (BCC) のフェルミ面

G: 0 0 0H ½ -1/2 ½
N00½
P¼¼¼
G000

平面波法
一次結合の基底関数として平面波を使う

φk (r ) = exp(ik ⋅ r )∑ Chkl uhkl (r )uhkl (r ) = exp[iG hkl ⋅ r ]波数 Ghkl の平面波は格子周期の関数の完全基底系:
すべての hkl について和を取れば、完全に正しい解になる=> 実際の計算では |Ghkl| < Gmax の範囲で近似する
H 11 − ES11
H 21 − ES 21

H n1 − ES n1uh 'k 'l ' H uhkl = ∫ e

H 12 − ES12
H 22 − ES ss
H n 2 − ES n 2−i (k + G h'k'l' )⋅r

H 1n − ES1n
H 2 n − ES 2 n=0 H nn − ES nn 2 2
 i (k +G hkl )⋅r∇ + V (r )edr
− 2m2 2= δ hkl ,h 'k ',l 'k + V * (G hkl − G h 'k 'l ' )2m
実際の計算のほとんどがポテンシャルのフーリエ変換=> GPUで高速化が容易

プログラム: 一次元平面波法

http://conf.msl.titech.ac.jp/jsap-crystal/

平面波基底による一次元バンド計算 pw1d.py

Usage: python pw1d.py
python pw1d.py (ft a na pottype bwidth bpot)
python pw1d.py (band a na pottype bwidth bpot nG kmin kmax nk)
python pw1d.py (wf a na pottype bwidth bpot nG kw iLevel xwmin xwmax nxw)
pottype: rect|gauss
実行例: python pw1d.py ft 5.4064 64 rect 0.5 10.0
ポテンシャルのフーリエ変換を表示。
格子定数5.4064Å、単位格子を 26 = 64分割 (FFTのためnaは2n)矩形ポテンシャル 0.5 Å幅、10.0 eV高さ
実行例: python pw1d.py band 5.4064 64 rect 0.5 10.0 3 -0.5 0.5 21
バンド構造を計算。構造、分割数、ポテンシャルは上と同じバンド構造を 逆空間内部座標 [-½ ½] (第一ブリルアンゾーン) で21分割して表示
実行例: python pw1d.py wf 5.4064 64 rect 0.5 10.0 3 0.0 0 0.0 16.2192 101
結晶波動関数を表示。構造、分割数、ポテンシャルは上と同じ
波数ベクトルはΓ点に近い3点を用いる。
k = 0.0 (Γ点), 固有解の 0 番目の準位の波動関数を、0.0 ~ 16.2192 オングストロームの範囲で101分割して表示
Energy levels:(注意: 固有解はエネルギー順にソートしていないので、コンソール出力のEnergy levels:で準位の番号を確認)

0120.624459 eV
6.39666 eV
6.08362 eV

プログラム: 一次元平面波法

pw1d.py
Si の格子定数 a = 5.4064 Å m* = 1.0me
ポテンシャル V(x): 障壁幅0.5 Å 障壁高さ 10.0 eV
python pw1d.py ft 5.4064 64 rect 0.5 10.0 9 -0.5 0.5 21

ポテンシャルのフーリエ変換係数
基底数 64

プログラム: 一次元平面波法

pw1d.py
Si の格子定数 a = 5.4064 Å m* = 1.0me
障壁幅0.5 Å 障壁高さ 10.0 eV python pw1d.py band 5.4064 64 rect 0.5 10.0 9 -0.5 0.5 21
python pw1d.py wf 5.4064 64 rect 0.5 10.0 9 0.0 2 0.0 17 501
k = 0.0, 1番目の準位 (E=0.608eV)の波動関数(注意: 準位はエネルギー順にソートしていないので、コンソール出力のEnergy levels:で準位の番号を確認)
Energy levels:
083.3381 eV
147.3631 eV
2 0.608457 eV
36.91001 eV
422.477 eV
583.3212 eV647.2599 eV
720.7386 eV
85.18836 eV

バンド構造
基底数 9V(x)自由電子V(x) = 0