kronig_penney.py ダウンロード/コピー
kronig_penney.py
kronig_penney.py
1import sys
2import numpy as np
3from numpy import sqrt, exp, sin, cos, tan, cosh, sinh
4import numpy.linalg as LA
5from pprint import pprint
6import csv
7from matplotlib import pyplot as plt
8
9"""
101D band calculation by Kronig-Penney model
11
12Kronig-Penneyモデルによる1次元バンド計算を行うスクリプト。
13半導体の電子バンド構造を、矩形ポテンシャルを用いたKronig-Penneyモデルで計算し、
14その結果をグラフ表示、バンド図プロット、波動関数プロットのいずれかのモードで出力する。
15
16関連リンク: :doc:`kronig_penney_usage`
17"""
18
19#===================================
20# physical constants
21#===================================
22pi = 3.14159265358979323846
23pi2 = 2.0 * pi
24h = 6.6260755e-34 # Js";
25hbar = 1.05459e-34 # "Js";
26c = 2.99792458e8 # m/s";
27e = 1.60218e-19 # C";
28e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
29kB = 1.380658e-23 # JK<sup>-1</sup>";
30me = 9.1093897e-31 # kg";
31R = 8.314462618 # J/K/mol
32a0 = 5.29177e-11 # m";
33
34
35#========================
36# global configuration
37#========================
38mode = 'graph' # graph|band|wf
39
40#========================
41# Crystal definition
42#========================
43# Si
44a = 5.4064 # angstrom, lattice parameter
45
46#========================
47# Potential
48#========================
49bwidth = 0.5 # A, barrier width
50bpot = 10.0 # eV, barrier height
51
52#=====================================
53# 解を走査するグラフ表示
54#=====================================
55kg = 0.0 # k point to be plotted
56# 解を走査するエネルギー範囲
57Emin = 0.0
58Emax = 9.5
59# グラフを表示するエネルギー点数
60nE = 51
61# 解を走査するエネルギー点数
62nEsearch = nE
63# Secant法パラメータ
64eps = 1.0e-8
65nmaxiter = 100
66dump = 0.0
67
68#========================
69# Band
70#========================
71kmin = -0.5 # in pi/a
72kmax = 0.5 # in pi/a
73nk = 21
74
75# プロットするエネルギー範囲
76Erange = [0.0, 10.0] # eV
77
78# リストに保存する準位最大数
79nMaxLevel = 15
80
81#========================
82# Wave function
83#========================
84#波動関数を描画するx範囲
85xwmin = 0.0 # A
86xwmax = 3.0 * a # A
87nxw = 101
88#描画する波動関数の波数
89kw = 0.0
90#描画する波動関数の準位番号
91iLevel = 0
92
93#===================================
94# figure configuration
95#===================================
96figsize = (6, 8)
97fontsize = 12
98legend_fontsize = 8
99
100
101#==============================================
102# fundamental functions
103#==============================================
104def pfloat(str):
105 """
106 概要: 文字列をfloat型に安全に変換する。
107 詳細説明: 変換できない場合はNoneを返す。
108 :param str: str 変換する文字列。
109 :returns: float または None 変換されたfloat値、またはNone。
110 """
111 try:
112 return float(str)
113 except:
114 return None
115
116def pint(str):
117 """
118 概要: 文字列をint型に安全に変換する。
119 詳細説明: 変換できない場合はNoneを返す。
120 :param str: str 変換する文字列。
121 :returns: int または None 変換されたint値、またはNone。
122 """
123 try:
124 return int(str)
125 except:
126 return None
127
128def getarg(position, defval = None):
129 """
130 概要: コマンドライン引数を安全に取得する。
131 詳細説明: 指定された位置の引数が存在しない場合、デフォルト値を返す。
132 :param position: int 取得する引数の位置（インデックス）。
133 :param defval: Any 引数が存在しない場合に返すデフォルト値。
134 :returns: Any 指定された位置の引数、またはデフォルト値。
135 """
136 try:
137 return sys.argv[position]
138 except:
139 return defval
140
141def getfloatarg(position, defval = None):
142 """
143 概要: コマンドライン引数をfloat型に変換して安全に取得する。
144 詳細説明: 指定された位置の引数が存在しない場合やfloatに変換できない場合、デフォルト値を返す。
145 :param position: int 取得する引数の位置（インデックス）。
146 :param defval: float または None 引数が存在しない場合や変換できない場合に返すデフォルト値。
147 :returns: float または None 変換されたfloat値の引数、またはデフォルト値。
148 """
149 return pfloat(getarg(position, defval))
150
151def getintarg(position, defval = None):
152 """
153 概要: コマンドライン引数をint型に変換して安全に取得する。
154 詳細説明: 指定された位置の引数が存在しない場合やintに変換できない場合、デフォルト値を返す。
155 :param position: int 取得する引数の位置（インデックス）。
156 :param defval: int または None 引数が存在しない場合や変換できない場合に返すデフォルト値。
157 :returns: int または None 変換されたint値の引数、またはデフォルト値。
158 """
159 return pint(getarg(position, defval))
160
161def round01(x, a):
162 """
163 概要: 数値を区間 [0, a) にマッピングし、周期的なオフセットを計算する。
164 詳細説明: x = n * a + x0 となるような x0 と整数 n を計算する。x0 は 0 <= x0 < a の範囲に収まる。
165 :param x: float マッピング対象の数値。
166 :param a: float 周期の幅。
167 :returns: tuple (float, int) x を a で割った余り x0 と整数 n のタプル。
168 """
169 if x >= 0.0:
170 n = int(x / a)
171 else:
172 n = int(x / a) - 1
173 x0 = x - n * a
174 return x0, n
175
176def usage():
177 """
178 概要: スクリプトの正しい使用方法を標準出力に表示する。
179 詳細説明: コマンドライン引数の形式といくつかの使用例を示す。
180 """
181 print("")
182 print("Usage: Variables in () are optional")
183 print(" python {}".format(sys.argv[0]))
184 print(" python {} (graph a bwidth bpot k Emin Emax nE)".format(sys.argv[0]))
185 print(" python {} (band a bwidth bpot nG kmin kmax nk)".format(sys.argv[0]))
186 print(" python {} (wf a bwidth bpot kw iLevel xwmin xwmax nxw)".format(sys.argv[0]))
187 print(" ex: python {} {} {} {} {} {} {} {} {}"
188 .format(sys.argv[0], 'graph', a, bwidth, bpot, kg, Emin, Emax, nE))
189 print(" ex: python {} {} {} {} {} {} {} {}"
190 .format(sys.argv[0], 'band', a, bwidth, bpot, kmin, kmax, nk))
191 print(" ex: python {} {} {} {} {} {} {} {} {} {}"
192 .format(sys.argv[0], 'wf', a, bwidth, bpot, kw, iLevel, xwmin, xwmax, nxw))
193
194def terminate(message = None):
195 """
196 概要: エラーメッセージを表示してプログラムを終了する。
197 詳細説明: usage() 関数を呼び出した後、指定されたメッセージ（任意）を表示してプログラムを終了する。
198 :param message: str または None 表示するエラーメッセージ（オプション）。
199 """
200 print("")
201 if message is not None:
202 print("")
203 print(message)
204 print("")
205
206 usage()
207 print("")
208 exit()
209
210#==============================================
211# update default values by startup arguments
212#==============================================
213argv = sys.argv
214#if len(argv) == 1:
215# terminate()
216
217mode = getarg (1, mode)
218a = getfloatarg(2, a)
219bwidth = getfloatarg(3, bwidth)
220bpot = getfloatarg(4, bpot)
221
222if mode == 'graph':
223 kg = getfloatarg(5, kg)
224 Emin = getfloatarg(6, Emin)
225 Emax = getfloatarg(7, Emax)
226 nE = getintarg (8, nE)
227elif mode == 'band':
228 kmin = getfloatarg( 5, kmin)
229 kmax = getfloatarg( 6, kmax)
230 nk = getintarg ( 7, nk)
231elif mode == 'wf':
232 kw = getfloatarg( 5, kw)
233 iLevel = getintarg ( 6, iLevel)
234 xwmin = getfloatarg( 7, xwmin)
235 xwmax = getfloatarg( 8, xwmax)
236 nxw = getintarg (9, nxw)
237
238
239def pot(x):
240 """
241 概要: 矩形ポテンシャル関数を定義する。
242 詳細説明: 0 <= x < a - bwidth の範囲ではポテンシャルは0、a - bwidth <= x < a の範囲では bpot のポテンシャルを返す。周期 a を持つ。
243 :param x: float 位置 (A)。
244 :returns: float 指定された位置 x におけるポテンシャル値 (eV)。
245 """
246 global a
247 global bwidth, bpot
248
249 xred, nred = round01(x, a)
250
251 if a - bwidth <= xred < a:
252 return bpot
253 return 0.0
254
255def build_potential(xmin, xstep, n):
256 """
257 概要: 指定された範囲でポテンシャルプロファイルを作成する。
258 詳細説明: xmin から開始し、n 個の点に対して xstep 間隔でポテンシャル値を計算する。
259 :param xmin: float 最初のx座標 (A)。
260 :param xstep: float x座標のステップ幅 (A)。
261 :param n: int 計算する点の数。
262 :returns: tuple (numpy.ndarray, numpy.ndarray) x座標の配列と対応するポテンシャル値の配列のタプル。
263 """
264 xpot = np.empty(n)
265 ypot = np.empty(n)
266 for i in range(n):
267 xx = xmin + i * xstep
268 xpot[i] = xx
269 ypot[i] = pot(xx)
270 return xpot, ypot
271
272def cal_delta(E, k, w, b, V0):
273 """
274 概要: Kronig-Penneyモデルのエネルギー方程式の左辺 delta(E) を計算する。
275 詳細説明: エネルギー E と波数 k、ポテンシャルパラメータ w、b、V0 を用いて、
276 周期的な結晶における電子のエネルギー準位を決定する方程式の値を計算する。
277 delta(E) = cos(ka) となる E が許容エネルギーとなる。
278 :param E: float 電子のエネルギー (eV)。
279 :param k: float ブロッホ波数 (単位: pi/a)。
280 :param w: float ポテンシャル井戸の幅 (A)。
281 :param b: float ポテンシャル障壁の幅 (A)。
282 :param V0: float ポテンシャル障壁の高さ (eV)。
283 :returns: float Kronig-Penneyモデル方程式の左辺の値。
284 """
285 alpha = sqrt(2.0 * me * E * e) / hbar
286 beta = sqrt(2.0 * me * (V0 - E) * e) / hbar
287 ka = k * pi2
288 alphaw = alpha * w * 1.0e-10
289 betab = beta * b * 1.0e-10
290 delta = (beta*beta - alpha*alpha)/2.0/alpha/beta * sin(alphaw) * sinh(betab) \
291 + cos(alphaw) * cosh(betab) \
292 - cos(ka)
293# print("a=", E, ka, alphaw, betab, delta)
294 return delta
295
296def check_ci(ci, kw, Ei, w, b, V0, eps, IsPrint = 0):
297 """
298 概要: 波束係数 ci がKronig-Penneyモデルの境界条件を満たしているか検証する (デバッグ用)。
299 詳細説明: 4つの境界条件式に ci を代入し、その結果が eps より小さいかを確認する。
300 満たさない場合はエラーで終了する。
301 :param ci: list of complex 波動関数の係数リスト。
302 :param kw: float ブロッホ波数 (単位: pi/a)。
303 :param Ei: float 電子のエネルギー (eV)。
304 :param w: float ポテンシャル井戸の幅 (A)。
305 :param b: float ポテンシャル障壁の幅 (A)。
306 :param V0: float ポテンシャル障壁の高さ (eV)。
307 :param eps: float 許容誤差。
308 :param IsPrint: int 詳細を出力するかどうかのフラグ (0: なし, 1: あり)。
309 """
310 alpha = sqrt(2.0 * me * Ei * e) / hbar
311 beta = sqrt(2.0 * me * (V0 - Ei) * e) / hbar
312 ka = kw * pi2
313 lambda_ = exp(1.0j * ka)
314 alphaw = alpha * w * 1.0e-10
315 betab = beta * b * 1.0e-10
316 alpha *= 1.0e-10
317 beta *= 1.0e-10
318
319 Passed = 1
320 vmax = 0.0
321 if 1:
322 Mij = np.empty([4, 4], dtype = complex)
323 M3ij = np.empty([3, 3], dtype = complex)
324 V3i = np.empty([3, 1], dtype = complex)
325 Mij[0, 0] = Mij[0, 1] = 1.0
326 Mij[0, 2] = Mij[0, 3] = -1.0
327 Mij[1, 0] = 1.0j * alpha
328 Mij[1, 1] = -1.0j * alpha
329 Mij[1, 2] = -beta
330 Mij[1, 3] = beta
331 Mij[2, 0] = exp( 1.0j * alphaw)
332 Mij[2, 1] = exp(-1.0j * alphaw)
333 Mij[2, 2] = -lambda_ * exp(-betab)
334 Mij[2, 3] = -lambda_ * exp( betab)
335 Mij[3, 0] = 1.0j * alpha * exp( 1.0j * alphaw)
336 Mij[3, 1] = -1.0j * alpha * exp(-1.0j * alphaw)
337 Mij[3, 2] = -lambda_ * beta * exp(-betab)
338 Mij[3, 3] = lambda_ * beta * exp( betab)
339 if IsPrint:
340 for i in range(4):
341 print(" ci[{}] = {:12.4g}+j{:12.4g}".format(i, ci[i].real, ci[i].imag))
342 for i in range(4):
343 v = Mij[i, 0] * ci[0] + Mij[i, 1] * ci[1] + Mij[i, 2] * ci[2] + Mij[i, 3] * ci[3]
344 v = abs(v)
345 if IsPrint:
346 print(" abs(Mij@ci[{}]) = {}".format(i, v), eps)
347 if v > eps:
348 Passed = 0
349 if v > vmax:
350 vmax = v
351
352 if not Passed:
353 print("Error: Mij @ ci is not zero: abs(Mij@ci)={} > eps={}".format(vmax, eps))
354 exit()
355
356
357def refine_E(E0, E1, nmaxiter, eps, dump, k, w, b, V0, IsPrint = 0):
358 """
359 概要: Secant法を用いてKronig-Penneyモデルの許容エネルギー E を高精度化する。
360 詳細説明: cal_delta(E)=0 となる E をSecant法で探索し、収束したエネルギー値を返す。
361 :param E0: float 探索範囲の下限エネルギー (eV)。
362 :param E1: float 探索範囲の上限エネルギー (eV)。
363 :param nmaxiter: int 最大繰り返し回数。
364 :param eps: float 収束判定の許容誤差。
365 :param dump: float 収束を助けるための微小な値。
366 :param k: float ブロッホ波数 (単位: pi/a)。
367 :param w: float ポテンシャル井戸の幅 (A)。
368 :param b: float ポテンシャル障壁の幅 (A)。
369 :param V0: float ポテンシャル障壁の高さ (eV)。
370 :param IsPrint: int 詳細を出力するかどうかのフラグ (0: なし, 1: あり)。
371 :returns: tuple (float, float, float) または (None, None, None)
372 収束したエネルギー (float)、最終的なエネルギー変化 (float)、最終的なdelta値 (float) のタプル。
373 収束しない場合はNoneのタプル。
374 """
375 delta0 = cal_delta(E0, k, w, b, V0)
376 delta1 = cal_delta(E1, k, w, b, V0)
377 for i in range(nmaxiter):
378 diff = (delta1 - delta0) / (E1 - E0)
379 if diff >= 0.0:
380 diff += dump
381 else:
382 diff = -(abs(diff) + dump)
383
384 dE = -delta1 / diff
385 E2 = E1 + dE
386 delta2 = cal_delta(E2, k, w, b, V0)
387
388 if abs(dE) < eps:
389 if IsPrint:
390 print(" converged at E = {:12.6g} with dE = {:12.6g} delta = {:12.6g}"
391 .format(E2, dE, delta2))
392 return E2, dE, delta2
393 else:
394 E0 = E1
395 E1 = E2
396 delta0 = delta1
397 delta1 = delta2
398 continue
399 else:
400 print(" Not converged for {} iterations.".format(nmaxiter))
401 print(" E = {:12.6g} with dE = {:12.6g} delta = {:12.6g}".format(E2, dE, delta2))
402 return None, None, None
403
404def find_Elist(Emin, Emax, nEsearch, k, w, b, V0):
405 """
406 概要: 指定されたエネルギー範囲でKronig-Penneyモデルの許容エネルギー準位と対応する係数を見つける。
407 詳細説明: Emin から Emax までを nEsearch 分割して cal_delta(E) の符号反転を探索し、
408 Secant法で厳密なエネルギー値を特定する。同時に波動関数の係数 ci も計算する。
409 :param Emin: float 探索するエネルギー範囲の最小値 (eV)。
410 :param Emax: float 探索するエネルギー範囲の最大値 (eV)。
411 :param nEsearch: int エネルギー範囲の探索ステップ数。
412 :param k: float ブロッホ波数 (単位: pi/a)。
413 :param w: float ポテンシャル井戸の幅 (A)。
414 :param b: float ポテンシャル障壁の幅 (A)。
415 :param V0: float ポテンシャル障壁の高さ (eV)。
416 :returns: tuple (list of float, list of list of complex)
417 許容エネルギー準位のリストと、各準位に対応する波動関数係数のリストのタプル。
418 """
419# nEsearch *= 100
420 Estep = (Emax - Emin) / (nEsearch - 1)
421# print("Estep=", Estep)
422 d0 = None
423 iband = 0
424 Elist = []
425 Alist = []
426 for iE in range(nEsearch):
427 E = Emin + iE * Estep
428 if E == 0.0:
429 continue
430 if V0 <= E:
431 break
432
433 delta = cal_delta(E, k, w, b, V0)
434
435 if d0 is None:
436 d0 = delta
437 continue
438 if d0 * delta < 0.0:
439 d0 = delta
440
441# print(" E[{}]={:12.6g} eV delta={:8.4g}".format(iband, E, delta))
442 E, dE, delta0 = refine_E(E - Estep, E, nmaxiter, eps, dump, k, w, b, V0, IsPrint = 0)
443 print(" E[{}]={:12.6g} eV dE={:12.6g} delta={:12.6g}".format(iband, E, dE, delta0))
444
445 Elist.append(E)
446# Elist.append(E - 0.5 * Estep)
447
448 alpha = sqrt(2.0 * me * E * e) / hbar
449 beta = sqrt(2.0 * me * (V0 - E) * e) / hbar
450 ka = k * pi2
451 lambda_ = exp(1.0j * ka)
452 alphaw = alpha * w * 1.0e-10
453 betab = beta * b * 1.0e-10
454 alpha *= 1.0e-10
455 beta *= 1.0e-10
456
457 Mij = np.empty([4, 4], dtype = complex)
458 M3ij = np.empty([3, 3], dtype = complex)
459 V3i = np.empty([3, 1], dtype = complex)
460 Mij[0, 0] = Mij[0, 1] = 1.0
461 Mij[0, 2] = Mij[0, 3] = -1.0
462 Mij[1, 0] = 1.0j * alpha
463 Mij[1, 1] = -1.0j * alpha
464 Mij[1, 2] = -beta
465 Mij[1, 3] = beta
466 Mij[2, 0] = exp( 1.0j * alphaw)
467 Mij[2, 1] = exp(-1.0j * alphaw)
468 Mij[2, 2] = -lambda_ * exp(-betab)
469 Mij[2, 3] = -lambda_ * exp( betab)
470 Mij[3, 0] = 1.0j * alpha * exp( 1.0j * alphaw)
471 Mij[3, 1] = -1.0j * alpha * exp(-1.0j * alphaw)
472 Mij[3, 2] = -lambda_ * beta * exp(-betab)
473 Mij[3, 3] = lambda_ * beta * exp( betab)
474
475 A = 1.0
476 M3ij[0, 0] = Mij[1, 1]
477 M3ij[0, 1] = Mij[1, 2]
478 M3ij[0, 2] = Mij[1, 3]
479 M3ij[1, 0] = Mij[2, 1]
480 M3ij[1, 1] = Mij[2, 2]
481 M3ij[1, 2] = Mij[2, 3]
482 M3ij[2, 0] = Mij[3, 1]
483 M3ij[2, 1] = Mij[3, 2]
484 M3ij[2, 2] = Mij[3, 3]
485 V3i[0, 0] = -A * Mij[1, 0]
486 V3i[1, 0] = -A * Mij[2, 0]
487 V3i[2, 0] = -A * Mij[3, 0]
488
489 Ai = LA.solve(M3ij, V3i)
490
491 ci = [A, Ai[0, 0], Ai[1, 0], Ai[2, 0]]
492# check_ci(ci, k, E, w, b, V0, 3.0e-3, IsPrint = 0)
493 Alist.append(ci)
494
495 E += Estep
496
497 return Elist, Alist
498
499def cal_wavefunction(ci, x, kw, Ei, w, b, V0):
500 """
501 概要: 指定された係数 ci を用いて、Kronig-Penneyモデルの波動関数 Psi(x) を計算する。
502 詳細説明: ブロッホの定理に基づき Psi(x) = exp(ikx) * u(x) の形式で波動関数を構築する。
503 u(x) は周期関数で、ポテンシャル障壁内外で異なる形を取る。
504 :param ci: list of complex 波動関数の係数リスト。
505 :param x: float 波動関数を評価する位置 (A)。
506 :param kw: float ブロッホ波数 (単位: pi/a)。
507 :param Ei: float 電子のエネルギー (eV)。
508 :param w: float ポテンシャル井戸の幅 (A)。
509 :param b: float ポテンシャル障壁の幅 (A)。
510 :param V0: float ポテンシャル障壁の高さ (eV)。
511 :returns: complex 位置 x における複素波動関数の値。
512 """
513 IsPrint = 1
514
515 a = w + b
516
517 xmin = -b
518 xmax = w
519 x0, n = round01(x, a)
520 if x0 < -xmin:
521 x0 += a
522 if x0 >= xmax:
523 x0 -= a
524 if not xmin <= x0 < xmax:
525 print("Error: x0 out of range: x={:8.4g} {} x0={:8.4g} w={:8.4g} b={:8.4g}".format(x, n, x0, w, b))
526 exit()
527
528# if IsPrint:
529# print("x={:8.4g} {} x0={:8.4g} w={:8.4g} b={:8.4g}".format(x, n, x0, w, b))
530
531# check_ci(ci, kw, Ei, w, b, V0, 3.0e-3)
532
533 alpha = sqrt(2.0 * me * Ei * e) / hbar
534 beta = sqrt(2.0 * me * (V0 - Ei) * e) / hbar
535 alpha *= 1.0e-10
536 beta *= 1.0e-10
537 phase0 = pi2 / a * kw * x0
538 kph0 = exp(1.0j * phase0)
539
540# Calculate the periodic function u(x) from phi(x) in -b <= x < w
541 if xmin <= x0 < 0.0: # in barrier, defined in -b <= x < 0, w <= x < a
542 f = ci[2] * exp(beta * x0) + ci[3] * exp(-beta * x0)
543 u = f / kph0
544 else: # in well, defined in 0 <= x < w
545 f = ci[0] * exp(1.0j * alpha * x0) + ci[1] * exp(-1.0j * alpha * x0)
546 u = f / kph0
547
548# Calculate Bloch function phi(x) = exp(ikx) * u(x)
549 f = exp(1.0j * pi2 / a * kw * x) * u
550
551 return f + 0.0j
552# デバッグ用: 周期関数部分 u(x) を返す
553# return u + 0.0j
554
555
556def wf():
557 """
558 概要: Kronig-Penneyモデルにおける波動関数を計算し、プロットする。
559 詳細説明: 特定の波数 kw と準位番号 iLevel に対応するエネルギー準位を計算し、
560 その波動関数を xwmin から xwmax の範囲で描画する。
561 ポテンシャルプロファイルも重ねて表示する。
562 """
563 global mode
564 global a
565 global bwidth, bpot
566 global nEsearch, nMaxLevel
567 global kw, iLevel
568 global xwmin, xwmax, nxw
569
570 xwstep = (xwmax - xwmin) / (nxw - 1)
571 Estep = bpot / (nEsearch - 1)
572
573 print("")
574 print("=== Input parameterss ===")
575 print("mode:", mode)
576 print("a=", a, "A")
577 print("Wave function to be plotted: k = {} iLevel = {}".format(kw, iLevel))
578 print("x range: {} - {} at {} step, {} points".format(xwmin, xwmax, xwstep, nxw))
579 print("potential: w={} A h={} eV".format(bwidth, bpot))
580
581 print("")
582 V0 = bpot
583 b = bwidth
584 w = a - b
585
586 print("")
587 print("at k={:8.4g}".format(kw))
588
589 Elist, Alist = find_Elist(0.0, V0, nEsearch, kw, w, b, V0)
590
591 xplot, yplot = build_potential(xwmin, xwstep, nxw)
592
593 print("")
594 print("=== Calculate wave function ===")
595 print("Energy levels:", Elist, "eV")
596 print("at k = {}".format(kw))
597 print("{}-th energy level".format(iLevel))
598 Ei = Elist[iLevel]
599 ci = Alist[iLevel]
600
601 print(" E = {:12.6g} eV".format(Elist[iLevel]))
602 print(" A = {:12.4g}+j{:12.4g}".format(ci[0].real, ci[0].imag))
603 print(" B = {:12.4g}+j{:12.4g}".format(ci[1].real, ci[1].imag))
604 print(" C = {:12.4g}+j{:12.4g}".format(ci[2].real, ci[2].imag))
605 print(" D = {:12.4g}+j{:12.4g}".format(ci[3].real, ci[3].imag))
606 sumci = abs(ci[0] + ci[1] - ci[2] - ci[3])
607 print(" sum(ci) = {:12.4e}".format(sumci))
608 alpha = sqrt(2.0 * me * Ei * e) / hbar * 1.0e-10
609 beta = sqrt(2.0 * me * (V0 - Ei) * e) / hbar * 1.0e-10
610 print(" alpha = {:12.6g} A^-1".format(alpha))
611 print(" beta = {:12.6g} A^-1".format(beta))
612
613 print("")
614 print("Normalization")
615 nxintg = int(a / xwstep + 1.0001)
616 xintgstep = a / (nxintg - 1)
617 chg = 0.0
618 for i in range(nxintg):
619 x = 0.0 + i * xintgstep
620 yval = cal_wavefunction(ci, x, kw, Ei, w, b, V0)
621 chg += yval * yval.conjugate()
622
623 chg = chg.real * xintgstep
624 kywf = 1.0 / sqrt(chg)
625 print("integ(|psi(x)|^2) = ", chg)
626 print("Normalization coefficient = ", kywf)
627 for i in range(4):
628 ci[i] *= kywf
629 print(" A = {:12.4g}+j{:12.4g}".format(ci[0].real, ci[0].imag))
630 print(" B = {:12.4g}+j{:12.4g}".format(ci[1].real, ci[1].imag))
631 print(" C = {:12.4g}+j{:12.4g}".format(ci[2].real, ci[2].imag))
632 print(" D = {:12.4g}+j{:12.4g}".format(ci[3].real, ci[3].imag))
633
634 ywf = np.empty(nxw, dtype = complex)
635 for i in range(nxw):
636 x = xwmin + i * xwstep
637 ywf[i] = cal_wavefunction(ci, x, kw, Ei, w, b, V0)
638
639 charge = [(ywf[i] * ywf[i].conjugate()).real for i in range(nxw)]
640
641 fig = plt.figure(figsize = (16, 4)) #figsize)
642 ax2 = fig.add_subplot(1, 1, 1)
643# ax2 = fig.add_subplot(2, 1, 2)
644 ax1 = ax2.twinx()
645
646 ax1.set_xlim([xwmin, xwmax])
647 ax1.plot(xplot, yplot, linewidth = 0.5, label = 'U(x)')
648 ax1.plot(ax1.get_xlim(), [0.0, 0.0], color = 'r', linestyle = 'dashed', linewidth = 0.5)
649 ax2.set_xlim([xwmin, xwmax])
650 ax2.plot(xplot, ywf.real, color = 'r', linewidth = 1.5, label = "real")
651 ax2.plot(xplot, ywf.imag, color = 'b', linewidth = 1.5, label = "imaginary")
652 ax2.plot(xplot, charge, color = 'black', linewidth = 0.5, label = "charge")
653 ax2.plot(ax1.get_xlim(), [0.0, 0.0], color = 'r', linestyle = 'dashed', linewidth = 0.5)
654 ax1.set_xlabel("x (A)", fontsize = fontsize)
655 ax1.set_ylabel("U(x)", fontsize = fontsize)
656 ax2.set_xlabel("x (A)", fontsize = fontsize)
657 ax2.set_ylabel("$\Psi$($x$)", fontsize = fontsize)
658
659 handler1, label1 = ax1.get_legend_handles_labels()
660 handler2, label2 = ax2.get_legend_handles_labels()
661 ax2.legend(handler1 + handler2, label1 + label2, loc = 2, borderaxespad = 0.0, fontsize = legend_fontsize)
662# ax2.legend(fontsize = legend_fontsize)
663
664 ax1.tick_params(labelsize = fontsize)
665 ax2.tick_params(labelsize = fontsize)
666 plt.tight_layout()
667
668 plt.pause(0.1)
669 print("Press ENTER to exit>>", end = '')
670 input()
671
672 terminate()
673
674def band():
675 """
676 概要: Kronig-Penneyモデルのエネルギーバンド構造を計算し、プロットする。
677 詳細説明: 指定された波数範囲 (kmin から kmax) で複数の k 点に対して許容エネルギー準位を計算し、
678 E-k図（バンド構造）をプロットする。
679 """
680 global mode
681 global a
682 global bwidth, bpot
683 global kmin, kmax, nk
684 global nEsearch, nMaxLevel
685
686 kstep = (kmax - kmin) / (nk - 1)
687
688 print("")
689 print("=== Input parameterss ===")
690 print("mode:", mode)
691 print("a=", a, "A")
692 print("potential: w={} A h={} eV".format(bwidth, bpot))
693 print("k range: {} - {} at {} step, {} points".format(kmin, kmax, kstep, nk))
694 print("")
695
696 print("")
697 V0 = bpot
698 b = bwidth
699 w = a - b
700
701 xk = [kmin + i * kstep for i in range(nk)]
702 yE = np.zeros([nMaxLevel, nk])
703 nMaxBand = 0
704 for ik in range(nk):
705 k = kmin + ik * kstep
706 print("at k={:8.4g}".format(k))
707
708 Elist, Alist = find_Elist(0.0, V0, nEsearch, k, w, b, V0)
709 n = len(Elist)
710 if n > nMaxBand:
711 nMaxBand = n
712 for iband in range(min(n, nMaxLevel)):
713 yE[iband][ik] = Elist[iband]
714
715
716 fig = plt.figure(figsize = figsize)
717 ax1 = fig.add_subplot(1, 1, 1)
718
719 ax1.set_xlim([-0.5, 0.5])
720 ax1.set_ylim(Erange)
721# ax1.set_ylim([0.0, ax1.get_ylim()[1]])
722
723 for iL in range(nMaxBand):
724 ax1.plot(xk, yE[iL], linestyle = '', marker = 'o', markersize = 5.0,
725 markerfacecolor = 'none', markeredgecolor = 'black', markeredgewidth = 0.5)
726 ax1.set_xlabel("$k$ $(\pi$$/a)$", fontsize = fontsize)
727 ax1.set_ylabel("E (eV)", fontsize = fontsize)
728 ax1.legend(fontsize = legend_fontsize)
729 plt.tight_layout()
730
731 plt.pause(0.1)
732
733 print("Press ENTER to exit>>", end = '')
734 input()
735
736 terminate()
737
738def graphview():
739 """
740 概要: Kronig-Penneyモデルの許容エネルギー方程式の delta(E) 関数をプロットする。
741 詳細説明: 特定の波数 kg に対して、エネルギー E の関数としての cal_delta(E) の値を計算し、
742 グラフ表示する。delta(E)=0 となる点が許容エネルギー準位を示す。
743 """
744 global mode
745 global a
746 global bwidth, bpot
747 global Emin, Emax, nE
748
749 Estep = (Emax - Emin) / (nE - 1)
750
751 V0 = bpot
752 b = bwidth
753 w = a - b
754
755 print("")
756 print("=== Input parameterss ===")
757 print("mode:", mode)
758 print("a=", a, "A")
759 print(" barrier: w={} A h={} eV".format(b, V0))
760 print(" well : w={} A h={} eV".format(w, 0.0))
761 print("Energy range: {} - {}, {} eV step {} points".format(Emin, Emax, Estep, nE))
762 print("at k = {}".format(kg))
763 print("")
764
765 xE = []
766 yD = []
767 for i in range(1, nE):
768 E = Emin + i * Estep
769 if V0 <= E:
770 break
771
772 delta = cal_delta(E, kg, w, b, V0)
773
774 xE.append(E)
775 yD.append(delta)
776
777 fig = plt.figure(figsize = figsize)
778 ax1 = fig.add_subplot(1, 1, 1)
779
780 ax1.plot(xE, yD)
781 ax1.set_xlim([Emin, Emax])
782 ax1.plot([Emin, Emax], [0.0, 0.0], linestyle = 'dashed', color = 'r', linewidth = 0.5)
783 ax1.set_xlabel("E (eV)", fontsize = fontsize)
784 ax1.set_ylabel("delta", fontsize = fontsize)
785# ax1.legend(fontsize = legend_fontsize)
786 ax1.tick_params(labelsize = fontsize)
787 plt.tight_layout()
788
789 plt.pause(0.1)
790 print("Press ENTER to exit>>", end = '')
791 input()
792
793 terminate()
794
795def main():
796 """
797 概要: プログラムのエントリポイント。
798 詳細説明: コマンドライン引数で指定された mode に応じて、graphview、band、または wf 関数を呼び出す。
799 無効なモードが指定された場合はエラーメッセージを表示して終了する。
800 """
801 global mode
802
803 if mode == 'graph':
804 graphview()
805 elif mode == 'band':
806 band()
807 elif mode == 'wf':
808 wf()
809 else:
810 terminate("Error: Invalid mode [{}]".format(mode))
811
812
813if __name__ == "__main__":
814 main()