【第２５話】量子ビットの作り方３：静電場の与えたときの電子状態【量子ドットコンピュータの原理③】

############################################################################
#　無限に深い量子井戸中の電子に静電場を加えたときの固有状態（ステップ３）
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate
import numpy as np
import numpy.linalg as LA

#図全体
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']

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

######################################
#　物理系の設定
######################################
#量子井戸の幅
L = 1 * 10**-9
#計算区間
x_min = -L / 2.0
x_max = L / 2.0
#状態数
n_max = 10
#行列の要素数
DIM = n_max + 1
#空間分割数
NX = 500
#空間刻み間隔
nm = 1.0 * 10**-9
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#電場の強さ
Ex = 1.0 * 10**10
#電場の強さ
Ex_max = 1.0 * 10**10
#電場強度の分割数
NEx = 10

#静電場無しの固有関数
def verphi(n, x):
	kn = np.pi * (n + 1) / L
	return np.sqrt(2.0 / L) * np.sin(kn * (x + L / 2.0))

#ポテンシャル項
def V(x, Ex):
	return(e * Ex * x)
	
#静電場無しの固有エネルギー
def Energy(n):
	kn = np.pi * (n + 1) / L
	return hbar * hbar * kn * kn / (2.0 * me)

#被積分関数
def integral_matrixElement(x, n1, n2, Ex):
	return verphi(n1 ,x) * V(x, Ex) * verphi(n2, x) / eV


######################################
#　展開係数の計算
######################################
#固有値・固有ベクトルの初期化
eigenvalues = [0] * (NEx + 1)
vectors = [0] * (NEx + 1)
for nEx in range(NEx + 1):
	eigenvalues[nEx] = []
	vectors[nEx] = []




#存在確率分布グラフ描画用の配列初期化
xs = []
phi = [0] * (NEx + 1)
for nEx in range(NEx + 1):
	phi[nEx] = [0] * 2


#中心の電場依存性グラフ描画用の配列初期化
averageX = [0] * 2
for n in range(len(averageX)):
	averageX[n] = [0] * (NEx + 1)


######################################
#　静電場を加えた固有関数
######################################
def verphi_Ex(nEx, n, x):
	phi = 0 + 0j
	for m in range(n_max+1):
		phi += vectors[nEx][n][m] * verphi(m, x)	
	return phi


######################################
#　#静電場強度ごとの固有関数の計算
######################################
for nEx in range(NEx + 1):
	print("電場強度：" + str( nEx * 100 / NEx ) + "%")
	#静電場の強度を設定
	Ex = Ex_max / NEx * nEx
	#エルミート行列（リスト）
	matrix=[]

	###行列要素の計算
	for n1 in range(n_max + 1):
		col=[]
		for n2 in range(n_max + 1):
			#ガウス・ルジャンドル積分
			result = integrate.quad(
				integral_matrixElement, #被積分関数
				x_min, x_max,		        #積分区間の下端と上端
				args=(n1, n2, Ex)			  #被積分関数へ渡す引数
			)
			real = result[0]
			imag = 0j
			#無静電場のエネルギー固有値（対角成分）
			En = Energy(n1)/eV if (n1 == n2) else 0
			#行の要素を追加
			col.append( En + real )
		#行を追加
		matrix.append( col )

	#リスト → 行列
	matrix = np.array( matrix )

	###固有値と固有ベクトルの計算
	result = LA.eig( matrix )
	eig = result[0] #固有値
	vec = result[1] #固有ベクトル

	#小さい順に並べるためのインデックス（配列）
	index = np.argsort( eig )
	#固有値を小さい順に並び替え
	eigenvalues[nEx] = eig[ index ]

	#転置行列
	vec = vec.T
	#固有ベクトルの並び替え
	vectors[nEx] = vec[ index ]

	### 検算：MA-EA=0 ?
	sum = 0
	for i in range(DIM):
		v = matrix @ vectors[nEx][i] - eigenvalues[nEx][i] * vectors[nEx][i]
		for j in range(DIM):
			sum += abs(v[j])**2
	print("|MA-EA| =" + str(sum))


#アニメーション作成用
ims = []
for nEx in range(NEx + 1):
	###固有関数の空間分布
	phi[nEx][0] = abs(verphi_Ex(nEx, 0, x))**2 / (1.0 * 10**9)
	phi[nEx][1] = abs(verphi_Ex(nEx, 1, x))**2 / (1.0 * 10**9)
	#各コマを描画
	img  = plt.plot(x/nm, phi[nEx][0], colors[0], linestyle='solid', linewidth = 5 )
	img += plt.plot(x/nm, phi[nEx][1], colors[1], linestyle='solid', linewidth = 5 )
	time = plt.text(0.4, 4.1, "電場強度：" + str( nEx * 100 / NEx ) + "%" , ha='left', va='center', fontsize=12, fontname="Yu Gothic")
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )


#########################################################################
# グラフの描画（固有関数）
#########################################################################
plt.title( u"無限に深い量子井戸に静電場を加えたときの固有状態", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ |\varphi^{(m)}(x)|^2$", fontsize=30)
#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.12, top = 0.95)
#罫線の描画
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.6, 0.6])
plt.ylim([0, 4.0])

#井戸の概形
plt.vlines([-0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0],  -0.5, 0.5, "black", linestyles='solid', linewidth = 10 )


#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

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