free_electron_band.py ダウンロード/コピー
free_electron_band.py
free_electron_band.py
1import sys
2from pprint import pprint
3import numpy as np
4from numpy import sqrt, exp, sin, cos, tan, pi
5import numpy.linalg as LA
6import csv
7from matplotlib import pyplot as plt
8
9"""
10概要:
11 自由電子モデルを用いたバンド計算を行うモジュール。
12
13関連リンク:
14 free_electron_band_usage
15"""
16
17#===================================
18# physical constants
19#===================================
20pi = 3.14159265358979323846
21pi2 = 2.0 * pi
22h = 6.6260755e-34 # Js";
23hbar = 1.05459e-34 # "Js";
24c = 2.99792458e8 # m/s";
25e = 1.60218e-19 # C";
26e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
27kB = 1.380658e-23 # JK<sup>-1</sup>";
28me = 9.1093897e-31 # kg";
29R = 8.314462618 # J/K/mol
30a0 = 5.29177e-11 # m";
31
32
33#========================
34# Crystal definition
35#========================
36# Si
37a = 5.4064 # angstrom, lattice parameter
38rg = np.zeros([3, 3])
39rg[0][0] = pow(2.0 * pi / a, 2) # reciprocal space metric in angstrom^-2
40rg[1][1] = rg[0][0]
41rg[2][2] = rg[0][0]
42
43# 有効質量
44meff = 1.0 # in me
45
46# E(k) = KE * k*k, E(k) in eV, k in angstrom^1
47KE = hbar * hbar / 2.0 / (meff * me) * 1.0e20 / e
48print("KE = ", KE)
49
50#========================
51# Band plot parameters
52#========================
53# バンド構造をプロットするk点の軌跡: [kx, ky, kz, k点名称]
54klist = [
55# [-0.5, 0.0, 0.0, "-X"]
56# , [0.0, 0.0, 0.0, r"$\Gamma$"]
57# , [ 0.5, 0.0, 0.0, "X"]
58 [0.5, 0.0, 1.0, "W"]
59 , [0.5, 0.5, 0.5, "L"]
60 , [0.0, 0.0, 0.0, r"$\Gamma$"]
61 , [0.0, 0.0, 1.0, "X"]
62 , [0.5, 0.0, 1.0, "W"]
63 , [0.75, 0.0, 0.75, "K"]
64 ]
65# プロットするバンド構造E(k)のk点数の概数
66nk = 101
67
68# Ehkl(k)を計算するhkl範囲
69hrange = [-3, 3]
70krange = [-3, 3]
71lrange = [-3, 3]
72#hrange = [0, 0]
73#krange = [0, 0]
74#lrange = [0, 0]
75
76# プロットするエネルギー範囲
77Erange = [0.0, 10.0] # eV
78
79
80#===================================
81# figure configuration
82#===================================
83figsize = (6, 8)
84#figsize = (6, 4)
85fontsize = 16
86legend_fontsize = 8
87
88
89#==============================================
90# fundamental functions
91#==============================================
92def pfloat(str):
93 """
94 概要:
95 文字列を実数値に変換する。
96 詳細説明:
97 実数値に変換できない文字列をfloat()で変換するとエラーになってプログラムが終了するため、
98 この関数は変換できなかったらNoneを返すが、プログラムは終了させない。
99 引数:
100 :param str: 変換元の文字列
101 :type str: str
102 戻り値:
103 :returns: 変換された実数値、変換できない場合はNone
104 :rtype: float or None
105 """
106 try:
107 return float(str)
108 except:
109 return None
110
111def pint(str):
112 """
113 概要:
114 文字列を整数値に変換する。
115 詳細説明:
116 pfloat()と同様に、変換できない文字列が渡されてもエラー終了させずNoneを返す。
117 引数:
118 :param str: 変換元の文字列
119 :type str: str
120 戻り値:
121 :returns: 変換された整数値、変換できない場合はNone
122 :rtype: int or None
123 """
124 try:
125 return int(str)
126 except:
127 return None
128
129def getarg(position, defval = None):
130 """
131 概要:
132 起動時のコマンドライン引数を安全に取得する。
133 詳細説明:
134 sys.argvに対して範囲外のインデックスを指定してもエラー終了せず、デフォルト値を返す。
135 引数:
136 :param position: 取得したい引数のインデックス
137 :type position: int
138 :param defval: 範囲外アクセス時に返すデフォルト値
139 :type defval: Any
140 戻り値:
141 :returns: 取得した引数の文字列、またはデフォルト値
142 :rtype: str or Any
143 """
144 try:
145 return sys.argv[position]
146 except:
147 return defval
148
149def getfloatarg(position, defval = None):
150 """
151 概要:
152 起動時引数を取得し、実数に変換して返す。
153 引数:
154 :param position: 取得したい引数のインデックス
155 :type position: int
156 :param defval: 範囲外アクセス時や変換失敗時に返すデフォルト値
157 :type defval: Any
158 戻り値:
159 :returns: 変換された実数値、またはデフォルト値
160 :rtype: float or Any
161 """
162 return pfloat(getarg(position, defval))
163
164def getintarg(position, defval = None):
165 """
166 概要:
167 起動時引数を取得し、整数値に変換して返す。
168 引数:
169 :param position: 取得したい引数のインデックス
170 :type position: int
171 :param defval: 範囲外アクセス時や変換失敗時に返すデフォルト値
172 :type defval: Any
173 戻り値:
174 :returns: 変換された整数値、またはデフォルト値
175 :rtype: int or Any
176 """
177 return pint(getarg(position, defval))
178
179def usage():
180 """
181 概要:
182 プログラムの使用方法を表示する。
183 詳細説明:
184 コンソール上に使い方のメッセージを出力する。
185 戻り値:
186 :returns: なし
187 :rtype: None
188 """
189 print("")
190 print("Usage:")
191# print(" python {}".format(sys.argv[0]))
192
193def terminate(message = None):
194 """
195 概要:
196 メッセージを表示し、プログラムを安全に終了する。
197 詳細説明:
198 任意のメッセージと使用方法を表示した上で exit() を呼び出す。
199 引数:
200 :param message: 終了前に表示する文字列（Noneの場合は表示しない）
201 :type message: str or None
202 戻り値:
203 :returns: なし
204 :rtype: None
205 """
206 print("")
207 if message is not None:
208 print("")
209 print(message)
210 print("")
211
212 usage()
213 print("")
214 exit()
215
216def cal_kdistance(rg, k0, k1):
217 """
218 概要:
219 逆格子の計量テンソルを用いて、2つのk点間の距離を計算する。
220 引数:
221 :param rg: 逆格子の計量テンソル(3x3行列)
222 :type rg: list or numpy.ndarray
223 :param k0: 始点のk点座標 [kx, ky, kz]
224 :type k0: list
225 :param k1: 終点のk点座標 [kx, ky, kz]
226 :type k1: list
227 戻り値:
228 :returns: 2点間の距離
229 :rtype: float
230 """
231 dkx = k1[0] - k0[0]
232 dky = k1[1] - k0[1]
233 dkz = k1[2] - k0[2]
234 r2 = rg[0][0] * dkx*dkx + rg[1][1] * dky*dky + rg[2][2] * dkz*dkz
235 r2 += 2.0 * (rg[0][1] * dkx*dky + rg[1][2] * dky*dkz + rg[2][0] * dky*dkx)
236
237 return sqrt(r2)
238
239def cal_E(k, Ghkl):
240 """
241 概要:
242 k点座標と逆格子ベクトルから自由電子のエネルギーを計算する(単位: eV)。
243 詳細説明:
244 kおよびGhklは内部座標系で与える。
245 引数:
246 :param k: k点の内部座標 [kx, ky, kz]
247 :type k: list
248 :param Ghkl: 逆格子ベクトル [h, k, l]
249 :type Ghkl: list
250 戻り値:
251 :returns: 計算されたエネルギー値 (eV)
252 :rtype: float
253 """
254 global rg
255
256 kabs2 = rg[0][0] * (k[0] + Ghkl[0])**2 # in angstrom^-2
257 kabs2 += rg[1][1] * (k[1] + Ghkl[1])**2
258 kabs2 += rg[2][2] * (k[2] + Ghkl[2])**2
259 E = KE * kabs2 # in eV
260 return E
261
262def get_dklist(klist, nk):
263 """
264 概要:
265 プロットするk点軌跡の各距離と、最初のk点からの累積距離リストを計算する。
266 詳細説明:
267 バンド構造プロット用に、k点の名称リストも合わせて生成する。
268 引数:
269 :param klist: [kx, ky, kz, 名称] の形式を持つk点の軌跡リスト
270 :type klist: list
271 :param nk: プロットするk点の概算数
272 :type nk: int
273 戻り値:
274 :returns: dklist, ktotal_list, ktotal_namelist, ktotal のタプル。
275 dklistは各k点間の距離リスト。
276 ktotal_listは最初のk点からの累積距離リスト。
277 ktotal_namelistはk点の名称リスト。
278 ktotalは累積距離の合計。
279 :rtype: tuple
280 """
281 print("")
282
283# list of k distances from the first k point of each k region
284 dklist = []
285# list of k distance from the first k point of the k list
286 ktotal = 0.0
287 ktotal_list = []
288 ktotal_namelist = []
289 ktotal_list.append(0.0)
290 ktotal_namelist.append(klist[0][3])
291 for i in range(1, len(klist)):
292 print("k [{:<10s}: ({:6.4f}, {:6.4f}, {:6.4f}] to [{:<10s}:({:6.4f}, {:6.4f}, {:6.4f}]"
293 .format(
294 klist[i-1][3], klist[i-1][0], klist[i-1][1], klist[i-1][2],
295 klist[i][3], klist[i][0], klist[i][1], klist[i][2]))
296
297 dk_ = cal_kdistance(rg,
298 [klist[i-1][0], klist[i-1][1], klist[i-1][2]],
299 [klist[i][0], klist[i][1], klist[i][2]] )
300 ktotal += dk_
301 dklist.append(dk_)
302 ktotal_list.append(ktotal)
303 ktotal_namelist.append(klist[i][3])
304# print(" dk={} ktotal={}".format(dk_, ktotal))
305
306 return dklist, ktotal_list, ktotal_namelist, ktotal
307
308def get_cal_klist(klist, nk):
309 """
310 概要:
311 指定されたk点軌跡から、等間隔になるように計算用のk点リストを生成する。
312 詳細説明:
313 プロットに必要な累積距離リストやk点名称などのデータもあわせて作成する。
314 引数:
315 :param klist: バンド構造をプロットするk点の軌跡リスト
316 :type klist: list
317 :param nk: プロットするk点数の概数
318 :type nk: int
319 戻り値:
320 :returns: xk, xkvec, ktotallist, ktotal_namelist, res のタプル。
321 xkは各計算k点の累積距離リスト。
322 xkvecは各計算k点の座標リスト。
323 ktotallistは区切りとなるk点の累積距離リスト。
324 ktotal_namelistは区切りとなるk点の名称リスト。
325 resは距離情報や分割数情報を持つ辞書。
326 :rtype: tuple
327 """
328 dklist, ktotallist, ktotal_namelist, ktotal = get_dklist(klist, nk)
329 kstep = ktotal / nk
330
331 nklist = []
332 xk = []
333 xkvec = []
334 ktotal_ = 0.0
335 for i in range(1, len(klist)):
336 nk_ = int(dklist[i-1] / kstep + 1.00001)
337 nklist.append(nk_)
338 if i == len(klist) - 1:
339 ndiv = nk_ - 1
340 else:
341 ndiv = nk_
342 kstepx_ = (klist[i][0] - klist[i-1][0]) / ndiv
343 kstepy_ = (klist[i][1] - klist[i-1][1]) / ndiv
344 kstepz_ = (klist[i][2] - klist[i-1][2]) / ndiv
345
346 print("k: {:<10s} - {:<10s}".format(klist[i-1][3], klist[i][3]), end = '')
347 print(": nk={:3d} ktotal={:6.4f} kstep=({:6.4f}, {:6.4f}, {:6.4f})"
348 .format(nklist[i-1], ktotal_, kstepx_, kstepy_, kstepz_))
349
350 dk_ = cal_kdistance(rg, [0.0, 0.0, 0.0], [kstepx_,kstepy_, kstepz_])
351 for j in range(nklist[i-1]):
352 kx = klist[i-1][0] + j * kstepx_
353 ky = klist[i-1][1] + j * kstepy_
354 kz = klist[i-1][2] + j * kstepz_
355
356 xk.append(ktotal_)
357 xkvec.append([kx, ky, kz])
358 ktotal_ += dk_
359
360 res = {"dklist": dklist, "nklist": nklist, "ktotal": ktotal, "kstep": kstep}
361
362 return xk, xkvec, ktotallist, ktotal_namelist, res
363
364def get_cal_Elist(xkvec, hrange, krange, lrange):
365 """
366 概要:
367 計算用k点リストと逆格子ベクトルの範囲から、各k点でのエネルギーリストを計算する。
368 引数:
369 :param xkvec: 計算を行うk点座標 [kx, ky, kz] のリスト
370 :type xkvec: list
371 :param hrange: 逆格子ベクトルhの計算範囲 [最小値, 最大値]
372 :type hrange: list
373 :param krange: 逆格子ベクトルkの計算範囲 [最小値, 最大値]
374 :type krange: list
375 :param lrange: 逆格子ベクトルlの計算範囲 [最小値, 最大値]
376 :type lrange: list
377 戻り値:
378 :returns: 各k点におけるエネルギー群が格納された入れ子リスト
379 :rtype: list
380 """
381 yE = []
382 for i in range(len(xkvec)):
383 kx = xkvec[i][0]
384 ky = xkvec[i][1]
385 kz = xkvec[i][2]
386 Elist = []
387 for ih in range(hrange[0], hrange[1]+1):
388 for ik in range(krange[0], krange[1]+1):
389 for il in range(lrange[0], lrange[1]+1):
390 E = cal_E([kx, ky, kz], [ih, ik, il])
391 Elist.append(E)
392
393 yE.append(Elist)
394
395 return yE
396
397def plot_band(axis, xk, yE, Erange, ktotallist, ktotal_namelist):
398 """
399 概要:
400 計算した自由電子のバンド構造をグラフにプロットする。
401 引数:
402 :param axis: matplotlibのAxesオブジェクト
403 :type axis: matplotlib.axes.Axes
404 :param xk: プロットするk点の累積距離のリスト
405 :type xk: list
406 :param yE: 各k点に対応するエネルギーのリスト(入れ子でも可)
407 :type yE: list
408 :param Erange: プロットするエネルギーの表示範囲 [最小値, 最大値]
409 :type Erange: list
410 :param ktotallist: k点境界における、最初のk点からの距離の和のリスト
411 :type ktotallist: list
412 :param ktotal_namelist: k点境界における、k点の名称のリスト
413 :type ktotal_namelist: list
414 戻り値:
415 :returns: なし
416 :rtype: None
417 """
418# 表示範囲は決め打ち
419 axis.set_xlim([min(xk), max(xk)])
420 axis.set_ylim(Erange)
421
422# バンド構造をプロット
423 axis.plot(xk, yE, linestyle = 'none',
424 marker = 'o', markerfacecolor = 'none', markeredgecolor = 'black',
425 markeredgewidth = 0.5, markersize = 2.0)
426
427# Γ点、BZ境界の縦線を引く
428 for i in range(1, len(ktotallist)-1):
429 axis.plot([ktotallist[i], ktotallist[i]], Erange,
430 linestyle = '-', color = 'black', linewidth = 0.5)
431
432# k軸の目盛りにk点の名称を表示する
433# グラフ枠が一つであれば plt.xtics()で設定できる
434# axisに対しては、.setpでattributeを直接書き換える必要があるらしい
435 plt.setp(axis, xticks = ktotallist, xticklabels = ktotal_namelist)
436 axis.set_xlabel("k", fontsize = fontsize)
437 axis.set_ylabel("E (eV)", fontsize = fontsize)
438 axis.tick_params(labelsize = fontsize)
439
440
441def main():
442 """
443 概要:
444 自由電子バンド計算のメインルーチン。
445 詳細説明:
446 格子定数やプロットパラメータを設定し、k点のリストとエネルギーを計算して
447 matplotlibを用いてバンド図を描画する。
448 戻り値:
449 :returns: なし
450 :rtype: None
451 """
452 global a, rg
453 global nk, klist
454
455 print("")
456 print("Lattice parameters: ({}) A".format(a))
457 print("Effective mass: {} me".format(meff))
458 print("Reciprocal lattice metric [A^-2]:")
459 pprint(rg)
460 print("khl range: h=[{}, {}] k=[{}, {}] l=[{}, {}]"
461 .format(*hrange, *krange, *lrange))
462 print("Plot E range: {} - {} eV".format(*Erange))
463
464 xk, xkvec, ktotallist, ktotal_namelist, res = get_cal_klist(klist, nk)
465
466 print("")
467 print("k vectors")
468 print(" k_total: {} A^-1".format(res["ktotal"]))
469 print(" nk : ", nk)
470 print(" kstep : {}".format(res["kstep"]))
471 print(" dklist")
472 pprint(res["dklist"])
473 print(" ktotallist:")
474 pprint(ktotallist)
475 print(" nk_list:", res["nklist"])
476# print("xk")
477# pprint(xk)
478
479 yE = get_cal_Elist(xkvec, hrange, krange, lrange)
480
481 print("")
482 print("plot")
483
484 fig = plt.figure(figsize = figsize)
485 ax1 = fig.add_subplot(1, 1, 1)
486
487 plot_band(ax1, xk, yE, Erange, ktotallist, ktotal_namelist)
488
489 plt.tight_layout()
490
491 plt.pause(0.1)
492 print("Press ENTER to exit>>", end = '')
493 input()
494
495 terminate()
496
497
498if __name__ == "__main__":
499 main()