【第15話】周期ポテンシャル（クローニッヒ・ペニーモデル）のエネルギーバンドの計算方法【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## クローニッヒ・ペニーモデルのエネルギーバンド
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
#虚数単位
I = 0.0 + 1.0j

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19

######################################
#　物理系の設定
######################################
#空間刻み間隔
nm = 1E-9
#井戸の形状
d1 = 1.5 * nm
d2 = 1.0 * nm
#周期
l = d1 + d2
#井戸の深さ
V1 = 0 * eV
V2 = -10.0 * eV

######################################
#　転送行列の計算
######################################

def TransferMatrix( k, n = 0 ):
	#各領域の波数
	k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
	k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))
	# 転送行列の定義
	T1  = np.zeros((2,2), dtype=np.complex)
	T2  = np.zeros((2,2), dtype=np.complex)
	M21 = np.zeros((2,2), dtype=np.complex)
	M12 = np.zeros((2,2), dtype=np.complex)
	# 境界
	M21[0][0] = (1.0 + 0.0j + k1 / k2) / 2.0
	M21[0][1] = (1.0 + 0.0j - k1 / k2) / 2.0
	M21[1][0] = (1.0 + 0.0j - k1 / k2) / 2.0
	M21[1][1] = (1.0 + 0.0j + k1 / k2) / 2.0
	M12[0][0] = (1.0 + 0.0j + k2 / k1) / 2.0
	M12[0][1] = (1.0 + 0.0j - k2 / k1) / 2.0
	M12[1][0] = (1.0 + 0.0j - k2 / k1) / 2.0
	M12[1][1] = (1.0 + 0.0j + k2 / k1) / 2.0
	T1[0][0] = np.exp( I * k1 * d1 )
	T1[0][1] = 0
	T1[1][0] = 0
	T1[1][1] = np.exp( - I * k1 * d1 )
	T2[0][0] = np.exp( I * k2 * d2 )
	T2[0][1] = 0
	T2[1][0] = 0
	T2[1][1] = np.exp( - I * k2 * d2 )
	#転送行列の計算
	M = M12 @ T2 @ M21 @ T1 

	for _n in range(n):
		M = M @ M

	#転送行列
	return M

#################################################################
## クローニッヒ・ペニーのモデルの分散関係
#################################################################
#計算点数
NK = 1000000
#エネルギーの計算範囲
E_min = 1.0/NK
E_max = 4.0
#エネルギーの行列
Es = np.linspace(E_min, E_max, NK)
#y値配列
ys = []
#Kl値配列と対応するエネルギー配列
Kls = []
EBs = []
#バンド番号
bn = 0
#Kl値配列と対応するエネルギー配列
BandEdgeEs = []
BandCenterEs = []
#バンド番号
BandEdgeNumber = 0
BandCenterNumber = 0
#バンド内フラグ
flag_out = True
flag_center = False
sign = 0

for E in Es:
	E = E * eV
	#波数の計算
	k = np.sqrt(np.complex128( 2.0 * me * E / hbar**2 ))
	M = TransferMatrix( k )
	y = ( M[0][0] + M[1][1]) / 2.0
	ys.append( y.real )

	if( abs(y.real) <= 1.0 ):
		if( flag_out == True ):
			flag_out = False
			BandEdgeEs.append( E/eV )
			Kls.append([])
			EBs.append([])
		Kl = np.arccos( y.real )
		Kls[ bn ].append( Kl / np.pi )
		EBs[ bn ].append( E / eV )

		if( sign * y.real/abs(y.real) < 0 ):
			BandCenterEs.append( E/eV )
		sign = y.real/abs(y.real)
	else:
		if( flag_out == False ):
			flag_out = True
			bn += 1
			BandEdgeEs.append( E/eV )

print(BandEdgeEs)
print(BandCenterEs)

#################################################################
## 転送行列のトレース
#################################################################
plt.title( r"$y=\frac{1}{2}{\rm Tr} (M)$" + u"（" + r"$l=2.5[{\rm nm}]$" + " , " + r"$d=1.0[{\rm nm}]$"+ " , " + r"$V=-10.0[{\rm eV}]$" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$E$" + r"$[{\rm eV}]$" , fontsize=20)
plt.ylabel(r"$y(E)$",  fontsize=20)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)
#描画範囲を設定
plt.xlim([0, 4])
plt.ylim([-2.5, 2.5])
#壁の描画
plt.hlines([-1], 0, 10, "blue", linestyles='dashed', linewidth = 3 )
plt.hlines([1],  0, 10, "blue", linestyles='dashed', linewidth = 3 )
#トレース値のグラフの描画
plt.plot(Es, ys, linestyle='solid', linewidth = 5)

#################################################################
## クローニッヒ・ペニーのモデルのエネルギーバンド
#################################################################
fig2 = plt.figure(figsize=(15, 8))
plt.title( u"クローニッヒ・ペニーモデルのエネルギーバンド（" + r"$l=2.5[{\rm nm}]$" + " , " + r"$d=1.0[{\rm nm}]$"+ " , " + r"$V=-10.0[{\rm eV}]$" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$E$" + r"$[{\rm eV}]$" , fontsize=20)
plt.xlabel(r"$Kl/\pi$",  fontsize=20)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([0, 1])
plt.ylim([0, 4])

#分散関係の描画
for n in range(bn):
	plt.plot(Kls[n], EBs[n], linewidth = 5, color = colors[0])

#グラフの表示
plt.show()