import sys
import numpy as np
from numpy import sqrt, exp, sin, cos, tan, cosh, sinh
import numpy.linalg as LA
from pprint import pprint
import csv
from matplotlib import pyplot as plt
"""
1D band calculation by Kronig-Penney model
Kronig-Penneyモデルによる1次元バンド計算を行うスクリプト。
半導体の電子バンド構造を、矩形ポテンシャルを用いたKronig-Penneyモデルで計算し、
その結果をグラフ表示、バンド図プロット、波動関数プロットのいずれかのモードで出力する。
関連リンク: :doc:`kronig_penney_usage`
"""
#===================================
# physical constants
#===================================
pi = 3.14159265358979323846
pi2 = 2.0 * pi
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
kB = 1.380658e-23 # JK<sup>-1</sup>";
me = 9.1093897e-31 # kg";
R = 8.314462618 # J/K/mol
a0 = 5.29177e-11 # m";
#========================
# global configuration
#========================
mode = 'graph' # graph|band|wf
#========================
# Crystal definition
#========================
# Si
a = 5.4064 # angstrom, lattice parameter
#========================
# Potential
#========================
bwidth = 0.5 # A, barrier width
bpot = 10.0 # eV, barrier height
#=====================================
# 解を走査するグラフ表示
#=====================================
kg = 0.0 # k point to be plotted
# 解を走査するエネルギー範囲
Emin = 0.0
Emax = 9.5
# グラフを表示するエネルギー点数
nE = 51
# 解を走査するエネルギー点数
nEsearch = nE
# Secant法パラメータ
eps = 1.0e-8
nmaxiter = 100
dump = 0.0
#========================
# Band
#========================
kmin = -0.5 # in pi/a
kmax = 0.5 # in pi/a
nk = 21
# プロットするエネルギー範囲
Erange = [0.0, 10.0] # eV
# リストに保存する準位最大数
nMaxLevel = 15
#========================
# Wave function
#========================
#波動関数を描画するx範囲
xwmin = 0.0 # A
xwmax = 3.0 * a # A
nxw = 101
#描画する波動関数の波数
kw = 0.0
#描画する波動関数の準位番号
iLevel = 0
#===================================
# figure configuration
#===================================
figsize = (6, 8)
fontsize = 12
legend_fontsize = 8
#==============================================
# fundamental functions
#==============================================
[ドキュメント]
def pfloat(str):
"""
概要: 文字列をfloat型に安全に変換する。
詳細説明: 変換できない場合はNoneを返す。
:param str: str 変換する文字列。
:returns: float または None 変換されたfloat値、またはNone。
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
概要: 文字列をint型に安全に変換する。
詳細説明: 変換できない場合はNoneを返す。
:param str: str 変換する文字列。
:returns: int または None 変換されたint値、またはNone。
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
概要: コマンドライン引数を安全に取得する。
詳細説明: 指定された位置の引数が存在しない場合、デフォルト値を返す。
:param position: int 取得する引数の位置（インデックス）。
:param defval: Any 引数が存在しない場合に返すデフォルト値。
:returns: Any 指定された位置の引数、またはデフォルト値。
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
概要: コマンドライン引数をfloat型に変換して安全に取得する。
詳細説明: 指定された位置の引数が存在しない場合やfloatに変換できない場合、デフォルト値を返す。
:param position: int 取得する引数の位置（インデックス）。
:param defval: float または None 引数が存在しない場合や変換できない場合に返すデフォルト値。
:returns: float または None 変換されたfloat値の引数、またはデフォルト値。
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
概要: コマンドライン引数をint型に変換して安全に取得する。
詳細説明: 指定された位置の引数が存在しない場合やintに変換できない場合、デフォルト値を返す。
:param position: int 取得する引数の位置（インデックス）。
:param defval: int または None 引数が存在しない場合や変換できない場合に返すデフォルト値。
:returns: int または None 変換されたint値の引数、またはデフォルト値。
"""
return pint(getarg(position, defval))
[ドキュメント]
def round01(x, a):
"""
概要: 数値を区間 [0, a) にマッピングし、周期的なオフセットを計算する。
詳細説明: x = n * a + x0 となるような x0 と整数 n を計算する。x0 は 0 <= x0 < a の範囲に収まる。
:param x: float マッピング対象の数値。
:param a: float 周期の幅。
:returns: tuple (float, int) x を a で割った余り x0 と整数 n のタプル。
"""
if x >= 0.0:
n = int(x / a)
else:
n = int(x / a) - 1
x0 = x - n * a
return x0, n
[ドキュメント]
def usage():
"""
概要: スクリプトの正しい使用方法を標準出力に表示する。
詳細説明: コマンドライン引数の形式といくつかの使用例を示す。
"""
print("")
print("Usage: Variables in () are optional")
print(" python {}".format(sys.argv[0]))
print(" python {} (graph a bwidth bpot k Emin Emax nE)".format(sys.argv[0]))
print(" python {} (band a bwidth bpot nG kmin kmax nk)".format(sys.argv[0]))
print(" python {} (wf a bwidth bpot kw iLevel xwmin xwmax nxw)".format(sys.argv[0]))
print(" ex: python {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'graph', a, bwidth, bpot, kg, Emin, Emax, nE))
print(" ex: python {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'band', a, bwidth, bpot, kmin, kmax, nk))
print(" ex: python {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'wf', a, bwidth, bpot, kw, iLevel, xwmin, xwmax, nxw))
[ドキュメント]
def terminate(message = None):
"""
概要: エラーメッセージを表示してプログラムを終了する。
詳細説明: usage() 関数を呼び出した後、指定されたメッセージ（任意）を表示してプログラムを終了する。
:param message: str または None 表示するエラーメッセージ（オプション）。
"""
print("")
if message is not None:
print("")
print(message)
print("")
usage()
print("")
exit()
#==============================================
# update default values by startup arguments
#==============================================
argv = sys.argv
#if len(argv) == 1:
# terminate()
mode = getarg (1, mode)
a = getfloatarg(2, a)
bwidth = getfloatarg(3, bwidth)
bpot = getfloatarg(4, bpot)
if mode == 'graph':
kg = getfloatarg(5, kg)
Emin = getfloatarg(6, Emin)
Emax = getfloatarg(7, Emax)
nE = getintarg (8, nE)
elif mode == 'band':
kmin = getfloatarg( 5, kmin)
kmax = getfloatarg( 6, kmax)
nk = getintarg ( 7, nk)
elif mode == 'wf':
kw = getfloatarg( 5, kw)
iLevel = getintarg ( 6, iLevel)
xwmin = getfloatarg( 7, xwmin)
xwmax = getfloatarg( 8, xwmax)
nxw = getintarg (9, nxw)
[ドキュメント]
def pot(x):
"""
概要: 矩形ポテンシャル関数を定義する。
詳細説明: 0 <= x < a - bwidth の範囲ではポテンシャルは0、a - bwidth <= x < a の範囲では bpot のポテンシャルを返す。周期 a を持つ。
:param x: float 位置 (A)。
:returns: float 指定された位置 x におけるポテンシャル値 (eV)。
"""
global a
global bwidth, bpot
xred, nred = round01(x, a)
if a - bwidth <= xred < a:
return bpot
return 0.0
[ドキュメント]
def build_potential(xmin, xstep, n):
"""
概要: 指定された範囲でポテンシャルプロファイルを作成する。
詳細説明: xmin から開始し、n 個の点に対して xstep 間隔でポテンシャル値を計算する。
:param xmin: float 最初のx座標 (A)。
:param xstep: float x座標のステップ幅 (A)。
:param n: int 計算する点の数。
:returns: tuple (numpy.ndarray, numpy.ndarray) x座標の配列と対応するポテンシャル値の配列のタプル。
"""
xpot = np.empty(n)
ypot = np.empty(n)
for i in range(n):
xx = xmin + i * xstep
xpot[i] = xx
ypot[i] = pot(xx)
return xpot, ypot
[ドキュメント]
def cal_delta(E, k, w, b, V0):
"""
概要: Kronig-Penneyモデルのエネルギー方程式の左辺 delta(E) を計算する。
詳細説明: エネルギー E と波数 k、ポテンシャルパラメータ w、b、V0 を用いて、
周期的な結晶における電子のエネルギー準位を決定する方程式の値を計算する。
delta(E) = cos(ka) となる E が許容エネルギーとなる。
:param E: float 電子のエネルギー (eV)。
:param k: float ブロッホ波数 (単位: pi/a)。
:param w: float ポテンシャル井戸の幅 (A)。
:param b: float ポテンシャル障壁の幅 (A)。
:param V0: float ポテンシャル障壁の高さ (eV)。
:returns: float Kronig-Penneyモデル方程式の左辺の値。
"""
alpha = sqrt(2.0 * me * E * e) / hbar
beta = sqrt(2.0 * me * (V0 - E) * e) / hbar
ka = k * pi2
alphaw = alpha * w * 1.0e-10
betab = beta * b * 1.0e-10
delta = (beta*beta - alpha*alpha)/2.0/alpha/beta * sin(alphaw) * sinh(betab) \
+ cos(alphaw) * cosh(betab) \
- cos(ka)
# print("a=", E, ka, alphaw, betab, delta)
return delta
[ドキュメント]
def check_ci(ci, kw, Ei, w, b, V0, eps, IsPrint = 0):
"""
概要: 波束係数 ci がKronig-Penneyモデルの境界条件を満たしているか検証する (デバッグ用)。
詳細説明: 4つの境界条件式に ci を代入し、その結果が eps より小さいかを確認する。
満たさない場合はエラーで終了する。
:param ci: list of complex 波動関数の係数リスト。
:param kw: float ブロッホ波数 (単位: pi/a)。
:param Ei: float 電子のエネルギー (eV)。
:param w: float ポテンシャル井戸の幅 (A)。
:param b: float ポテンシャル障壁の幅 (A)。
:param V0: float ポテンシャル障壁の高さ (eV)。
:param eps: float 許容誤差。
:param IsPrint: int 詳細を出力するかどうかのフラグ (0: なし, 1: あり)。
"""
alpha = sqrt(2.0 * me * Ei * e) / hbar
beta = sqrt(2.0 * me * (V0 - Ei) * e) / hbar
ka = kw * pi2
lambda_ = exp(1.0j * ka)
alphaw = alpha * w * 1.0e-10
betab = beta * b * 1.0e-10
alpha *= 1.0e-10
beta *= 1.0e-10
Passed = 1
vmax = 0.0
if 1:
Mij = np.empty([4, 4], dtype = complex)
M3ij = np.empty([3, 3], dtype = complex)
V3i = np.empty([3, 1], dtype = complex)
Mij[0, 0] = Mij[0, 1] = 1.0
Mij[0, 2] = Mij[0, 3] = -1.0
Mij[1, 0] = 1.0j * alpha
Mij[1, 1] = -1.0j * alpha
Mij[1, 2] = -beta
Mij[1, 3] = beta
Mij[2, 0] = exp( 1.0j * alphaw)
Mij[2, 1] = exp(-1.0j * alphaw)
Mij[2, 2] = -lambda_ * exp(-betab)
Mij[2, 3] = -lambda_ * exp( betab)
Mij[3, 0] = 1.0j * alpha * exp( 1.0j * alphaw)
Mij[3, 1] = -1.0j * alpha * exp(-1.0j * alphaw)
Mij[3, 2] = -lambda_ * beta * exp(-betab)
Mij[3, 3] = lambda_ * beta * exp( betab)
if IsPrint:
for i in range(4):
print(" ci[{}] = {:12.4g}+j{:12.4g}".format(i, ci[i].real, ci[i].imag))
for i in range(4):
v = Mij[i, 0] * ci[0] + Mij[i, 1] * ci[1] + Mij[i, 2] * ci[2] + Mij[i, 3] * ci[3]
v = abs(v)
if IsPrint:
print(" abs(Mij@ci[{}]) = {}".format(i, v), eps)
if v > eps:
Passed = 0
if v > vmax:
vmax = v
if not Passed:
print("Error: Mij @ ci is not zero: abs(Mij@ci)={} > eps={}".format(vmax, eps))
exit()
[ドキュメント]
def refine_E(E0, E1, nmaxiter, eps, dump, k, w, b, V0, IsPrint = 0):
"""
概要: Secant法を用いてKronig-Penneyモデルの許容エネルギー E を高精度化する。
詳細説明: cal_delta(E)=0 となる E をSecant法で探索し、収束したエネルギー値を返す。
:param E0: float 探索範囲の下限エネルギー (eV)。
:param E1: float 探索範囲の上限エネルギー (eV)。
:param nmaxiter: int 最大繰り返し回数。
:param eps: float 収束判定の許容誤差。
:param dump: float 収束を助けるための微小な値。
:param k: float ブロッホ波数 (単位: pi/a)。
:param w: float ポテンシャル井戸の幅 (A)。
:param b: float ポテンシャル障壁の幅 (A)。
:param V0: float ポテンシャル障壁の高さ (eV)。
:param IsPrint: int 詳細を出力するかどうかのフラグ (0: なし, 1: あり)。
:returns: tuple (float, float, float) または (None, None, None)
収束したエネルギー (float)、最終的なエネルギー変化 (float)、最終的なdelta値 (float) のタプル。
収束しない場合はNoneのタプル。
"""
delta0 = cal_delta(E0, k, w, b, V0)
delta1 = cal_delta(E1, k, w, b, V0)
for i in range(nmaxiter):
diff = (delta1 - delta0) / (E1 - E0)
if diff >= 0.0:
diff += dump
else:
diff = -(abs(diff) + dump)
dE = -delta1 / diff
E2 = E1 + dE
delta2 = cal_delta(E2, k, w, b, V0)
if abs(dE) < eps:
if IsPrint:
print(" converged at E = {:12.6g} with dE = {:12.6g} delta = {:12.6g}"
.format(E2, dE, delta2))
return E2, dE, delta2
else:
E0 = E1
E1 = E2
delta0 = delta1
delta1 = delta2
continue
else:
print(" Not converged for {} iterations.".format(nmaxiter))
print(" E = {:12.6g} with dE = {:12.6g} delta = {:12.6g}".format(E2, dE, delta2))
return None, None, None
[ドキュメント]
def find_Elist(Emin, Emax, nEsearch, k, w, b, V0):
"""
概要: 指定されたエネルギー範囲でKronig-Penneyモデルの許容エネルギー準位と対応する係数を見つける。
詳細説明: Emin から Emax までを nEsearch 分割して cal_delta(E) の符号反転を探索し、
Secant法で厳密なエネルギー値を特定する。同時に波動関数の係数 ci も計算する。
:param Emin: float 探索するエネルギー範囲の最小値 (eV)。
:param Emax: float 探索するエネルギー範囲の最大値 (eV)。
:param nEsearch: int エネルギー範囲の探索ステップ数。
:param k: float ブロッホ波数 (単位: pi/a)。
:param w: float ポテンシャル井戸の幅 (A)。
:param b: float ポテンシャル障壁の幅 (A)。
:param V0: float ポテンシャル障壁の高さ (eV)。
:returns: tuple (list of float, list of list of complex)
許容エネルギー準位のリストと、各準位に対応する波動関数係数のリストのタプル。
"""
# nEsearch *= 100
Estep = (Emax - Emin) / (nEsearch - 1)
# print("Estep=", Estep)
d0 = None
iband = 0
Elist = []
Alist = []
for iE in range(nEsearch):
E = Emin + iE * Estep
if E == 0.0:
continue
if V0 <= E:
break
delta = cal_delta(E, k, w, b, V0)
if d0 is None:
d0 = delta
continue
if d0 * delta < 0.0:
d0 = delta
# print(" E[{}]={:12.6g} eV delta={:8.4g}".format(iband, E, delta))
E, dE, delta0 = refine_E(E - Estep, E, nmaxiter, eps, dump, k, w, b, V0, IsPrint = 0)
print(" E[{}]={:12.6g} eV dE={:12.6g} delta={:12.6g}".format(iband, E, dE, delta0))
Elist.append(E)
# Elist.append(E - 0.5 * Estep)
alpha = sqrt(2.0 * me * E * e) / hbar
beta = sqrt(2.0 * me * (V0 - E) * e) / hbar
ka = k * pi2
lambda_ = exp(1.0j * ka)
alphaw = alpha * w * 1.0e-10
betab = beta * b * 1.0e-10
alpha *= 1.0e-10
beta *= 1.0e-10
Mij = np.empty([4, 4], dtype = complex)
M3ij = np.empty([3, 3], dtype = complex)
V3i = np.empty([3, 1], dtype = complex)
Mij[0, 0] = Mij[0, 1] = 1.0
Mij[0, 2] = Mij[0, 3] = -1.0
Mij[1, 0] = 1.0j * alpha
Mij[1, 1] = -1.0j * alpha
Mij[1, 2] = -beta
Mij[1, 3] = beta
Mij[2, 0] = exp( 1.0j * alphaw)
Mij[2, 1] = exp(-1.0j * alphaw)
Mij[2, 2] = -lambda_ * exp(-betab)
Mij[2, 3] = -lambda_ * exp( betab)
Mij[3, 0] = 1.0j * alpha * exp( 1.0j * alphaw)
Mij[3, 1] = -1.0j * alpha * exp(-1.0j * alphaw)
Mij[3, 2] = -lambda_ * beta * exp(-betab)
Mij[3, 3] = lambda_ * beta * exp( betab)
A = 1.0
M3ij[0, 0] = Mij[1, 1]
M3ij[0, 1] = Mij[1, 2]
M3ij[0, 2] = Mij[1, 3]
M3ij[1, 0] = Mij[2, 1]
M3ij[1, 1] = Mij[2, 2]
M3ij[1, 2] = Mij[2, 3]
M3ij[2, 0] = Mij[3, 1]
M3ij[2, 1] = Mij[3, 2]
M3ij[2, 2] = Mij[3, 3]
V3i[0, 0] = -A * Mij[1, 0]
V3i[1, 0] = -A * Mij[2, 0]
V3i[2, 0] = -A * Mij[3, 0]
Ai = LA.solve(M3ij, V3i)
ci = [A, Ai[0, 0], Ai[1, 0], Ai[2, 0]]
# check_ci(ci, k, E, w, b, V0, 3.0e-3, IsPrint = 0)
Alist.append(ci)
E += Estep
return Elist, Alist
[ドキュメント]
def cal_wavefunction(ci, x, kw, Ei, w, b, V0):
"""
概要: 指定された係数 ci を用いて、Kronig-Penneyモデルの波動関数 Psi(x) を計算する。
詳細説明: ブロッホの定理に基づき Psi(x) = exp(ikx) * u(x) の形式で波動関数を構築する。
u(x) は周期関数で、ポテンシャル障壁内外で異なる形を取る。
:param ci: list of complex 波動関数の係数リスト。
:param x: float 波動関数を評価する位置 (A)。
:param kw: float ブロッホ波数 (単位: pi/a)。
:param Ei: float 電子のエネルギー (eV)。
:param w: float ポテンシャル井戸の幅 (A)。
:param b: float ポテンシャル障壁の幅 (A)。
:param V0: float ポテンシャル障壁の高さ (eV)。
:returns: complex 位置 x における複素波動関数の値。
"""
IsPrint = 1
a = w + b
xmin = -b
xmax = w
x0, n = round01(x, a)
if x0 < -xmin:
x0 += a
if x0 >= xmax:
x0 -= a
if not xmin <= x0 < xmax:
print("Error: x0 out of range: x={:8.4g} {} x0={:8.4g} w={:8.4g} b={:8.4g}".format(x, n, x0, w, b))
exit()
# if IsPrint:
# print("x={:8.4g} {} x0={:8.4g} w={:8.4g} b={:8.4g}".format(x, n, x0, w, b))
# check_ci(ci, kw, Ei, w, b, V0, 3.0e-3)
alpha = sqrt(2.0 * me * Ei * e) / hbar
beta = sqrt(2.0 * me * (V0 - Ei) * e) / hbar
alpha *= 1.0e-10
beta *= 1.0e-10
phase0 = pi2 / a * kw * x0
kph0 = exp(1.0j * phase0)
# Calculate the periodic function u(x) from phi(x) in -b <= x < w
if xmin <= x0 < 0.0: # in barrier, defined in -b <= x < 0, w <= x < a
f = ci[2] * exp(beta * x0) + ci[3] * exp(-beta * x0)
u = f / kph0
else: # in well, defined in 0 <= x < w
f = ci[0] * exp(1.0j * alpha * x0) + ci[1] * exp(-1.0j * alpha * x0)
u = f / kph0
# Calculate Bloch function phi(x) = exp(ikx) * u(x)
f = exp(1.0j * pi2 / a * kw * x) * u
return f + 0.0j
# デバッグ用: 周期関数部分 u(x) を返す
# return u + 0.0j
[ドキュメント]
def wf():
"""
概要: Kronig-Penneyモデルにおける波動関数を計算し、プロットする。
詳細説明: 特定の波数 kw と準位番号 iLevel に対応するエネルギー準位を計算し、
その波動関数を xwmin から xwmax の範囲で描画する。
ポテンシャルプロファイルも重ねて表示する。
"""
global mode
global a
global bwidth, bpot
global nEsearch, nMaxLevel
global kw, iLevel
global xwmin, xwmax, nxw
xwstep = (xwmax - xwmin) / (nxw - 1)
Estep = bpot / (nEsearch - 1)
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("a=", a, "A")
print("Wave function to be plotted: k = {} iLevel = {}".format(kw, iLevel))
print("x range: {} - {} at {} step, {} points".format(xwmin, xwmax, xwstep, nxw))
print("potential: w={} A h={} eV".format(bwidth, bpot))
print("")
V0 = bpot
b = bwidth
w = a - b
print("")
print("at k={:8.4g}".format(kw))
Elist, Alist = find_Elist(0.0, V0, nEsearch, kw, w, b, V0)
xplot, yplot = build_potential(xwmin, xwstep, nxw)
print("")
print("=== Calculate wave function ===")
print("Energy levels:", Elist, "eV")
print("at k = {}".format(kw))
print("{}-th energy level".format(iLevel))
Ei = Elist[iLevel]
ci = Alist[iLevel]
print(" E = {:12.6g} eV".format(Elist[iLevel]))
print(" A = {:12.4g}+j{:12.4g}".format(ci[0].real, ci[0].imag))
print(" B = {:12.4g}+j{:12.4g}".format(ci[1].real, ci[1].imag))
print(" C = {:12.4g}+j{:12.4g}".format(ci[2].real, ci[2].imag))
print(" D = {:12.4g}+j{:12.4g}".format(ci[3].real, ci[3].imag))
sumci = abs(ci[0] + ci[1] - ci[2] - ci[3])
print(" sum(ci) = {:12.4e}".format(sumci))
alpha = sqrt(2.0 * me * Ei * e) / hbar * 1.0e-10
beta = sqrt(2.0 * me * (V0 - Ei) * e) / hbar * 1.0e-10
print(" alpha = {:12.6g} A^-1".format(alpha))
print(" beta = {:12.6g} A^-1".format(beta))
print("")
print("Normalization")
nxintg = int(a / xwstep + 1.0001)
xintgstep = a / (nxintg - 1)
chg = 0.0
for i in range(nxintg):
x = 0.0 + i * xintgstep
yval = cal_wavefunction(ci, x, kw, Ei, w, b, V0)
chg += yval * yval.conjugate()
chg = chg.real * xintgstep
kywf = 1.0 / sqrt(chg)
print("integ(|psi(x)|^2) = ", chg)
print("Normalization coefficient = ", kywf)
for i in range(4):
ci[i] *= kywf
print(" A = {:12.4g}+j{:12.4g}".format(ci[0].real, ci[0].imag))
print(" B = {:12.4g}+j{:12.4g}".format(ci[1].real, ci[1].imag))
print(" C = {:12.4g}+j{:12.4g}".format(ci[2].real, ci[2].imag))
print(" D = {:12.4g}+j{:12.4g}".format(ci[3].real, ci[3].imag))
ywf = np.empty(nxw, dtype = complex)
for i in range(nxw):
x = xwmin + i * xwstep
ywf[i] = cal_wavefunction(ci, x, kw, Ei, w, b, V0)
charge = [(ywf[i] * ywf[i].conjugate()).real for i in range(nxw)]
fig = plt.figure(figsize = (16, 4)) #figsize)
ax2 = fig.add_subplot(1, 1, 1)
# ax2 = fig.add_subplot(2, 1, 2)
ax1 = ax2.twinx()
ax1.set_xlim([xwmin, xwmax])
ax1.plot(xplot, yplot, linewidth = 0.5, label = 'U(x)')
ax1.plot(ax1.get_xlim(), [0.0, 0.0], color = 'r', linestyle = 'dashed', linewidth = 0.5)
ax2.set_xlim([xwmin, xwmax])
ax2.plot(xplot, ywf.real, color = 'r', linewidth = 1.5, label = "real")
ax2.plot(xplot, ywf.imag, color = 'b', linewidth = 1.5, label = "imaginary")
ax2.plot(xplot, charge, color = 'black', linewidth = 0.5, label = "charge")
ax2.plot(ax1.get_xlim(), [0.0, 0.0], color = 'r', linestyle = 'dashed', linewidth = 0.5)
ax1.set_xlabel("x (A)", fontsize = fontsize)
ax1.set_ylabel("U(x)", fontsize = fontsize)
ax2.set_xlabel("x (A)", fontsize = fontsize)
ax2.set_ylabel("$\Psi$($x$)", fontsize = fontsize)
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax2.legend(handler1 + handler2, label1 + label2, loc = 2, borderaxespad = 0.0, fontsize = legend_fontsize)
# ax2.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def band():
"""
概要: Kronig-Penneyモデルのエネルギーバンド構造を計算し、プロットする。
詳細説明: 指定された波数範囲 (kmin から kmax) で複数の k 点に対して許容エネルギー準位を計算し、
E-k図（バンド構造）をプロットする。
"""
global mode
global a
global bwidth, bpot
global kmin, kmax, nk
global nEsearch, nMaxLevel
kstep = (kmax - kmin) / (nk - 1)
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("a=", a, "A")
print("potential: w={} A h={} eV".format(bwidth, bpot))
print("k range: {} - {} at {} step, {} points".format(kmin, kmax, kstep, nk))
print("")
print("")
V0 = bpot
b = bwidth
w = a - b
xk = [kmin + i * kstep for i in range(nk)]
yE = np.zeros([nMaxLevel, nk])
nMaxBand = 0
for ik in range(nk):
k = kmin + ik * kstep
print("at k={:8.4g}".format(k))
Elist, Alist = find_Elist(0.0, V0, nEsearch, k, w, b, V0)
n = len(Elist)
if n > nMaxBand:
nMaxBand = n
for iband in range(min(n, nMaxLevel)):
yE[iband][ik] = Elist[iband]
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(1, 1, 1)
ax1.set_xlim([-0.5, 0.5])
ax1.set_ylim(Erange)
# ax1.set_ylim([0.0, ax1.get_ylim()[1]])
for iL in range(nMaxBand):
ax1.plot(xk, yE[iL], linestyle = '', marker = 'o', markersize = 5.0,
markerfacecolor = 'none', markeredgecolor = 'black', markeredgewidth = 0.5)
ax1.set_xlabel("$k$ $(\pi$$/a)$", fontsize = fontsize)
ax1.set_ylabel("E (eV)", fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def graphview():
"""
概要: Kronig-Penneyモデルの許容エネルギー方程式の delta(E) 関数をプロットする。
詳細説明: 特定の波数 kg に対して、エネルギー E の関数としての cal_delta(E) の値を計算し、
グラフ表示する。delta(E)=0 となる点が許容エネルギー準位を示す。
"""
global mode
global a
global bwidth, bpot
global Emin, Emax, nE
Estep = (Emax - Emin) / (nE - 1)
V0 = bpot
b = bwidth
w = a - b
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("a=", a, "A")
print(" barrier: w={} A h={} eV".format(b, V0))
print(" well : w={} A h={} eV".format(w, 0.0))
print("Energy range: {} - {}, {} eV step {} points".format(Emin, Emax, Estep, nE))
print("at k = {}".format(kg))
print("")
xE = []
yD = []
for i in range(1, nE):
E = Emin + i * Estep
if V0 <= E:
break
delta = cal_delta(E, kg, w, b, V0)
xE.append(E)
yD.append(delta)
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(1, 1, 1)
ax1.plot(xE, yD)
ax1.set_xlim([Emin, Emax])
ax1.plot([Emin, Emax], [0.0, 0.0], linestyle = 'dashed', color = 'r', linewidth = 0.5)
ax1.set_xlabel("E (eV)", fontsize = fontsize)
ax1.set_ylabel("delta", fontsize = fontsize)
# ax1.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def main():
"""
概要: プログラムのエントリポイント。
詳細説明: コマンドライン引数で指定された mode に応じて、graphview、band、または wf 関数を呼び出す。
無効なモードが指定された場合はエラーメッセージを表示して終了する。
"""
global mode
if mode == 'graph':
graphview()
elif mode == 'band':
band()
elif mode == 'wf':
wf()
else:
terminate("Error: Invalid mode [{}]".format(mode))
if __name__ == "__main__":
main()