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
from typing import Any, NoReturn, Tuple, List # 型ヒントのためにインポート
"""
1Dバンド計算 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(s: str) -> float | None:
"""文字列をfloat型に安全に変換します。
変換に失敗した場合でもプログラムを終了せず、Noneを返します。
:param s: 変換対象の文字列。
:returns: float型に変換された値、または変換失敗時にNone。
"""
try:
return float(s)
except ValueError:
return None
[ドキュメント]
def pint(s: str) -> int | None:
"""文字列をint型に安全に変換します。
変換に失敗した場合でもプログラムを終了せず、Noneを返します。
:param s: 変換対象の文字列。
:returns: int型に変換された値、または変換失敗時にNone。
"""
try:
return int(s)
except ValueError:
return None
[ドキュメント]
def getarg(position: int, defval: Any = None) -> Any:
"""コマンドライン引数を安全に取得します。
指定された位置の引数が存在しない場合でもエラーを発生させず、デフォルト値を返します。
:param position: 取得する引数の位置（インデックス）。
:param defval: 引数が存在しない場合に返すデフォルト値。
:returns: 指定された位置の引数、またはdefval。
"""
try:
return sys.argv[position]
except IndexError:
return defval
[ドキュメント]
def getfloatarg(position: int, defval: Any = None) -> float | None:
"""コマンドライン引数をfloat型に変換して安全に取得します。
指定された位置の引数を取得し、`pfloat`関数を使ってfloat型に変換します。
引数が存在しないか、変換できない場合はデフォルト値を返します。
:param position: 取得する引数の位置（インデックス）。
:param defval: 引数が存在しないか、変換できない場合に返すデフォルト値。
:returns: float型に変換された引数、またはdefval。
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position: int, defval: Any = None) -> int | None:
"""コマンドライン引数をint型に変換して安全に取得します。
指定された位置の引数を取得し、`pint`関数を使ってint型に変換します。
引数が存在しないか、変換できない場合はデフォルト値を返します。
:param position: 取得する引数の位置（インデックス）。
:param defval: 引数が存在しないか、変換できない場合に返すデフォルト値。
:returns: int型に変換された引数、またはdefval。
"""
return pint(getarg(position, defval))
[ドキュメント]
def round01(x: float, a: float) -> Tuple[float, int]:
"""数値を指定された周期で正規化し、周期内の値と周期数を返します。
`x`を`a`で割った余り`x0`と、その周期数`n`を計算します。
`x0`は`0 <= x0 < a`の範囲に正規化されます。
:param x: 正規化対象の数値。
:param a: 周期の幅。
:returns: (x0, n) - 周期内の値 (float) と周期数 (int) のタプル。
"""
if x >= 0.0:
n = int(x / a)
else:
n = int(x / a) - 1
x0 = x - n * a
return x0, n
[ドキュメント]
def usage() -> None:
"""プログラムの正しい使用方法を標準出力に表示します。
コマンドライン引数の形式と、各モード（graph, band, wf）での使用例を示します。
:returns: なし
"""
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: str | None = None) -> NoReturn:
"""オプションのメッセージを表示し、プログラムを終了します。
プログラムの実行を中断し、`usage`関数を呼び出して使用方法を表示した後、
終了コード0でプログラムを終了します。
:param message: 表示するエラーまたは終了メッセージ (オプション)。
:returns: なし (プログラム終了)
"""
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: float) -> float:
"""矩形バリアポテンシャルの値を計算します。
周期境界条件を考慮し、単位格子内のx座標に基づいてポテンシャル値を決定します。
バリア領域では`bpot`、それ以外では0.0を返します。
:param x: ポテンシャルを評価する位置 (A)。
:returns: 指定された位置でのポテンシャル値 (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: float, xstep: float, n: int) -> Tuple[np.ndarray, np.ndarray]:
"""指定された範囲とステップでポテンシャルプロファイルの配列を構築します。
`xmin`から開始し、`n`個の点でポテンシャル関数`pot(x)`を評価し、
x座標とポテンシャル値の配列を返します。
:param xmin: ポテンシャル評価の開始x座標 (A)。
:param xstep: x座標のステップ幅 (A)。
:param n: 評価する点の数。
:returns: (xpot, ypot) - x座標の配列 (numpy.ndarray) と対応するポテンシャル値の配列 (numpy.ndarray) のタプル。
"""
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: float, k: float, w: float, b: float, V0: float) -> float:
"""Kronig-Penneyモデルの特性方程式の左辺（誤差）を計算します。
電子のエネルギー`E`、波数`k`、ウェル幅`w`、バリア幅`b`、バリア高さ`V0`を用いて、
Blochの定理に基づいた特性方程式の値を計算します。この値がゼロに近づくエネルギーが
許容帯のエネルギー準位となります。
:param E: 電子のエネルギー (eV)。
:param k: 波数 (単位はpi/aの分数)。
:param w: ポテンシャルの井戸の幅 (A)。
:param b: ポテンシャルの障壁の幅 (A)。
:param V0: ポテンシャルの障壁の高さ (eV)。
:returns: Kronig-Penneyモデルの特性方程式の左辺の値 (float)。
"""
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: List[complex], kw: float, Ei: float, w: float, b: float, V0: float, eps: float, IsPrint: int = 0) -> None:
"""波動関数の係数`ci`が境界条件を満足するかを検証します（デバッグ用）。
計算された波動関数の係数`ci`が、各領域の境界での接続条件
（波動関数とその導関数の連続性）を満たしているかをチェックします。
`Mij @ ci`がゼロベクトルに近ければ、条件は満たされています。
:param ci: 波動関数の係数リスト [A, B, C, D]。
:param kw: 波数 (単位はpi/aの分数)。
:param Ei: 電子のエネルギー (eV)。
:param w: ポテンシャルの井戸の幅 (A)。
:param b: ポテンシャルの障壁の幅 (A)。
:param V0: ポテンシャルの障壁の高さ (eV)。
:param eps: 許容誤差。
:param IsPrint: デバッグ情報を表示するかどうかのフラグ (0:表示しない, 1:表示する)。
:returns: なし (条件が満たされない場合はプログラムを終了)。
"""
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 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: float, E1: float, nmaxiter: int, eps: float, dump: float, k: float, w: float, b: float, V0: float, IsPrint: int = 0) -> Tuple[float | None, float | None, float | None]:
"""Secant法を用いて、Kronig-Penneyモデルの特性方程式の根であるエネルギー準位を精密化します。
`cal_delta(E)`がゼロに近づくエネルギー`E`を見つけるために、
初期の2つのエネルギー値`E0`, `E1`からSecant法を適用し、
指定された収束条件`eps`または最大反復回数`nmaxiter`まで反復します。
:param E0: Secant法の初期エネルギー推定値1 (eV)。
:param E1: Secant法の初期エネルギー推定値2 (eV)。
:param nmaxiter: Secant法の最大反復回数。
:param eps: 収束判定の許容誤差。
:param dump: 収束を助けるための微小なオフセット値。
:param k: 波数 (単位はpi/aの分数)。
:param w: ポテンシャルの井戸の幅 (A)。
:param b: ポテンシャルの障壁の幅 (A)。
:param V0: ポテンシャルの障壁の高さ (eV)。
:param IsPrint: デバッグ情報を表示するかどうかのフラグ (0:表示しない, 1:表示する)。
:returns: (E2, dE, delta2) - 精密化されたエネルギー (float)、最後のエネルギー変化量 (float)、
最終的なdelta値 (float) のタプル。収束しない場合は(None, None, 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: float, Emax: float, nEsearch: int, k: float, w: float, b: float, V0: float) -> Tuple[List[float], List[List[complex]]]:
"""指定されたエネルギー範囲で、Kronig-Penneyモデルの許容エネルギー準位と
対応する波動関数係数を探索します。
`Emin`から`Emax`までエネルギーを走査し、`cal_delta(E)`の符号が反転する点を
Secant法 (`refine_E`) で精密化してエネルギー準位を特定します。
各準位に対して、波動関数の係数（A, B, C, D）も計算し、リストとして返します。
(注: 現在の実装では、`cal_delta`の最初の符号を基準に根を探索するため、
複数の根がある場合に全てを見つけるわけではありません。)
:param Emin: 探索するエネルギー範囲の下限 (eV)。
:param Emax: 探索するエネルギー範囲の上限 (eV)。
:param nEsearch: エネルギー範囲内の探索点数。
:param k: 波数 (単位はpi/aの分数)。
:param w: ポテンシャルの井戸の幅 (A)。
:param b: ポテンシャルの障壁の幅 (A)。
:param V0: ポテンシャルの障壁の高さ (eV)。
:returns: (Elist, Alist) - 許容エネルギー準位のリスト (list[float]) と、
各準位に対応する波動関数係数のリスト (list[list[complex]]) のタプル。
"""
# nEsearch *= 100
Estep = (Emax - Emin) / (nEsearch - 1)
# print("Estep=", Estep)
d0: float | None = None
iband = 0
Elist: List[float] = []
Alist: List[List[complex]] = []
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))
refined_E, dE, delta0 = refine_E(E - Estep, E, nmaxiter, eps, dump, k, w, b, V0, IsPrint = 0)
# refined_EがNoneの場合は収束しなかったためスキップ
if refined_E is None:
continue
# ibandは表示用であり、実際のバンドインデックスとしては使用されていないため、常に0
print(" E[{}]={:12.6g} eV dE={:12.6g} delta={:12.6g}".format(iband, refined_E, dE, delta0))
Elist.append(refined_E)
# Elist.append(E - 0.5 * Estep)
alpha = sqrt(2.0 * me * refined_E * e) / hbar
beta = sqrt(2.0 * me * (V0 - refined_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_coeff = 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_coeff * Mij[1, 0]
V3i[1, 0] = -A_coeff * Mij[2, 0]
V3i[2, 0] = -A_coeff * Mij[3, 0]
Ai = LA.solve(M3ij, V3i)
ci = [A_coeff, Ai[0, 0], Ai[1, 0], Ai[2, 0]]
check_ci(ci, k, E, w, b, V0, 3.0e-3, IsPrint = 0) # Eはループ変数Emin+iE*Estepの値が使用される
Alist.append(ci)
# E += Estep
return Elist, Alist
[ドキュメント]
def cal_wavefunction(ci: List[complex], x: float, kw: float, Ei: float, w: float, b: float, V0: float) -> complex:
"""Kronig-Penneyモデルにおける波動関数（Bloch関数）の値を計算します。
波動関数の係数`ci`と指定された位置`x`、波数`kw`、エネルギー`Ei`を用いて、
Bloch関数 `psi(x) = exp(ikx) * u(x)` の値を計算します。
ここで`u(x)`は周期関数です。
:param ci: 波動関数の係数リスト [A, B, C, D]。
:param x: 波動関数を評価するx座標 (A)。
:param kw: 波数 (単位はpi/aの分数)。
:param Ei: 電子のエネルギー (eV)。
:param w: ポテンシャルの井戸の幅 (A)。
:param b: ポテンシャルの障壁の幅 (A)。
:param V0: ポテンシャルの障壁の高さ (eV)。
:returns: 指定されたx座標における波動関数の値 (complex)。
"""
# IsPrint = 1
a_lattice = w + b
xmin = -b
xmax = w
x0, n = round01(x, a_lattice) # ここではa_lattice (w+b)を使用
# x0が単位格子 [-b, w) の範囲に正規化されるよう調整
# 元のコードのロジックを維持 (ラウンド01は[0, a)を返すため、調整が必要)
if x0 < -xmin: # x0 < b
x0 += a_lattice
if x0 >= xmax: # x0 >= w
x0 -= a_lattice
# 上記調整後の範囲チェック (デバッグ用)
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_lattice * kw * x0 # ここではa_latticeを使用
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_lattice
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_lattice * kw * x) * u # ここではa_latticeを使用
return f + 0.0j
# デバッグ用: 周期関数部分 u(x) を返す
# return u + 0.0j
[ドキュメント]
def wf() -> None:
"""選択されたエネルギー準位に対する波動関数とポテンシャルをグラフ表示します。
グローバル設定に基づいて、特定の波数`kw`と準位番号`iLevel`に対応する波動関数を
`find_Elist`と`cal_wavefunction`を使って計算します。
計算された波動関数は、その実部、虚部、電荷密度（絶対値の2乗）、およびポテンシャル関数と共にプロットされます。
:returns: なし (グラフ表示後、プログラム終了)
"""
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))
if not (0 <= iLevel < len(Elist)):
print(f"Error: iLevel {iLevel} is out of range for {len(Elist)} energy levels found.")
terminate()
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(r"$\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() -> None:
"""Kronig-Penneyモデルのエネルギーバンド構造をグラフ表示します。
グローバル設定の波数範囲(`kmin`から`kmax`)と点数(`nk`)に基づいて、
各波数`k`における許容エネルギー準位を`find_Elist`を使って計算します。
計算されたエネルギー準位は、`k`に対する`E`のプロットとして表示され、バンド構造を形成します。
:returns: なし (グラフ表示後、プログラム終了)
"""
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(r"$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() -> None:
"""Kronig-Penneyモデルの特性方程式の誤差`delta(E)`をエネルギー`E`の関数としてグラフ表示します。
特定の波数`kg`とエネルギー範囲(`Emin`から`Emax`)および点数(`nE`)に基づいて、
`cal_delta(E)`の値を計算しプロットします。このグラフは、許容エネルギー準位
（`delta(E) = 0`となる点）を視覚的に確認するために使用されます。
:returns: なし (グラフ表示後、プログラム終了)
"""
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: List[float] = []
yD: List[float] = []
for i in range(1, nE): # range(1, nE)なのでi=0はスキップされる
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() -> None:
"""プログラムの主要な実行フローを制御します。
コマンドライン引数から設定された`mode`に基づいて、
`graphview`、`band`、または`wf`のいずれかの関数を呼び出します。
無効な`mode`が指定された場合は、エラーメッセージを表示してプログラムを終了します。
:returns: なし
"""
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()