################################################################# ## [032]コヒーレント状態の時間発展 ################################################################# 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'] ################################################################# ## 物理定数 ################################################################# #プランク定数 h = 6.62607015 * 1.0E-34 hbar = h / (2.0 * np.pi) #電子の質量 me = 9.1093837015 * 1.0E-31 #電子ボルト eV = 1.602176634 * 1.0E-19 #ナノメートル nm = 1E-9 #虚数単位 I = 0.0 + 1.0j ################################################################# ## 規格化エルミート多項式 ################################################################# def NormalizedHermitePolynomial( n, x ): H0 = (1.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) H1 = (4.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) * x if(n==0): return H0 if(n==1): return H1 for m in range(2, n+1): H2 = np.sqrt( 2.0 / m ) * x * H1 - np.sqrt( (m - 1) / m ) * H0 H0 = H1 H1 = H2 return H2 ################################################################# ## 物理系の設定 ################################################################# #調和ポテンシャルの強さ omega = 1.0 * 1.0E+15 A = np.sqrt( me * omega / hbar) #固有関数 def verphi(n, x): barX = A * x return np.sqrt(A) * NormalizedHermitePolynomial( n, barX ) #固有エネルギー(eV) def Energy(n): return (n + 1.0/2.0) * hbar * omega ################################################################# ## コヒーレント状態 ################################################################# NC = 400 def CoherentState(v, x, t): _psi = 0.0 + 0.0j for n in range(NC): _psi += ff(n, v) * verphi(n, x) * np.exp( - I * omega * (n + 1.0/2.0) * t ) return _psi * np.exp( -abs(v)**2 / 2.0) def ff(n, v, ln = 1): if( n==0 ): return 1 if( n==1 ): return v * ln return ff( n-1, v, ln * v / np.sqrt(n) ) ######################################################################### # 波動関数 アニメーション ######################################################################### #計算時間の幅 ts = 0 te = 1000 #調和振動子の周期 T0 = 2 * 2.0 * np.pi / omega #時間間隔 dt = T0 / (te - ts + 1) plt.title( u"コヒーレント状態の時間発展", fontsize=20, fontname="Yu Gothic", fontweight=1000) plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=20) plt.ylabel(r"$ \psi_v(x,t) $", fontsize=30) L = 9 * nm #アニメーション作成用 ims = [] #描画範囲 x_min = -L/2 x_max = L/2 #描画区間数 NX = 500 #座標点配列の生成 x = np.linspace(x_min, x_max, NX) #罫線の描画 plt.grid(which = "major", axis = "x", alpha = 0.2, linestyle = "-", linewidth = 1) plt.grid(which = "major", axis = "y", alpha = 0.2, linestyle = "-", linewidth = 1) #描画範囲を設定 plt.xlim([x_min/nm, x_max/nm]) plt.ylim([-1.0, 1.0]) #調和ポテンシャルの概形 yE = 1.0 / 2.0 * me * omega**2 * x**2 /eV /30.0 ims = [] #各時刻における波動関数の計算 for tn in range(ts, te+1): #実時間の取得 t = dt * tn #波動関数の計算 v = 4.0 vx1 = CoherentState(1, x, t) vx4 = CoherentState(4, x, t) vx7 = CoherentState(7, x, t) vx1_real = vx1.real / np.sqrt(A) vx1_imag = vx1.imag / np.sqrt(A) vx4_real = vx4.real / np.sqrt(A) vx4_imag = vx4.imag / np.sqrt(A) vx7_real = vx7.real / np.sqrt(A) vx7_imag = vx7.imag / np.sqrt(A) #波動関数の表示 img = plt.plot( x/nm, yE, linestyle='dotted', color="black", linewidth = 1 ) img += plt.plot( x/nm, vx1_real, colors[0], linestyle='solid', linewidth = 3 ) img += plt.plot( x/nm, vx1_imag, colors[1], linestyle='solid', linewidth = 3 ) img += plt.plot( x/nm, vx4_real, colors[2], linestyle='solid', linewidth = 3 ) img += plt.plot( x/nm, vx4_imag, colors[3], linestyle='solid', linewidth = 3 ) img += plt.plot( x/nm, vx7_real, colors[4], linestyle='solid', linewidth = 3 ) img += plt.plot( x/nm, vx7_imag, colors[5], linestyle='solid', linewidth = 3 ) #アニメーションに追加 ims.append( img ) #余白の調整 #plt.subplots_adjust(left=0.15, right=0.90, bottom=0.1, top=0.99) plt.subplots_adjust(left = 0.12, right = 0.95, bottom = 0.10, top = 0.95) #アニメーションの生成 ani = animation.ArtistAnimation(fig, ims, interval=10) #アニメーションの保存 #ani.save("output.html", writer=animation.HTMLWriter()) #ani.save('anim.gif', writer='pillow') #ani.save("output.mp4", writer="ffmpeg", dpi=300) #グラフの表示 plt.show()
【第31話】調和ポテンシャル中の波動関数の初期配置の与え方【Pythonコピペで量子力学完全攻略マニュアル】
################################################################# ## 調和ポテンシャル中の電子の固有関数 初期値の与え方 ################################################################# import math import numpy as np import scipy.integrate as integrate 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'] ################################################################# ## 物理定数 ################################################################# #プランク定数 h = 6.62607015 * 1.0E-34 hbar = h / (2.0 * np.pi) #電子の質量 me = 9.1093837015 * 1.0E-31 #電子ボルト eV = 1.602176634 * 1.0E-19 #ナノメートル nm = 1E-9 #虚数単位 I = 0.0 + 1.0j ################################################################# ## 規格化エルミート多項式 ################################################################# def NormalizedHermitePolynomial( n, x ): H0 = (1.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) H1 = (4.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) * x if(n==0): return H0 if(n==1): return H1 for m in range(2, n+1): H2 = np.sqrt( 2.0 / m ) * x * H1 - np.sqrt( (m - 1) / m ) * H0 H0 = H1 H1 = H2 return H2 ################################################################# ## 物理系の設定 ################################################################# #量子井戸の幅(m) omega = 1.0 * 1.0E+15 A = np.sqrt( me * omega / hbar) #固有関数 def verphi(n, x): barX = A * x return np.sqrt(A) * NormalizedHermitePolynomial( n, barX ) #固有エネルギー(eV) def Energy(n): return (n + 1.0/2.0) * hbar * omega #状態数 NS = 100 ################################################################# ## 初期分布:ガウス分布 ################################################################# sigma = 1.0 * 10**10 x0 = 1.0 * nm def verphi0(x): return ( sigma**2 / np.pi )**(1.0/4.0) * np.exp( - 1.0/2.0 * sigma**2 * (x - x0)**2 ) #被積分関数 def integral_orthonomality(x, n): return verphi(n, x) * verphi0(x) L = 7 * nm #描画範囲 x_min = -L/2 x_max = L/2 cn = [] for n in range( NS + 1 ): #ガウス・ルジャンドル積分 result = integrate.quad( integral_orthonomality, #被積分関数 x_min, x_max, #積分区間の下端と上端 args = ( n ) #被積分関数へ渡す引数 ) cn.append( result[0] ) #ターミナルへ出力 print( "(" + str(n) + ") " + str( result[0]) ) ################################################################# ## 波動関数 ################################################################# def psi(x, t): _psi = 0.0 + 0.0j for n in range(NS): _psi += cn[n] * verphi(n, x) * np.exp( - I * ( n + 1.0/2.0) * omega * t ) return _psi ######################################################################### # 波動関数 アニメーション ######################################################################### #計算時間の幅 ts = 0 te = 100 #基底状態の周期 T0 = 2.0 * np.pi * hbar / Energy(0) print( "基底状態の周期:" + str(T0) + " [s]" ) #時間間隔 dt = T0 / (te - ts + 1) plt.title( u"調和ポテンシャル中の波動関数", fontsize=20, fontname="Yu Gothic", fontweight=1000) plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=20) L = 7 * nm #アニメーション作成用 ims = [] #描画範囲 x_min = -L/2 x_max = L/2 #描画区間数 NX = 500 #座標点配列の生成 x = np.linspace(x_min, x_max, NX) #罫線の描画 plt.grid(which = "major", axis = "x", alpha = 0.2, linestyle = "-", linewidth = 1) plt.grid(which = "major", axis = "y", alpha = 0.2, linestyle = "-", linewidth = 1) #描画範囲を設定 plt.xlim([-3.5, 3.5]) plt.ylim([-1.5, 6.5]) #調和ポテンシャルの概形 yE = 1.0 / 2.0 * me * omega**2 * x**2 /eV /5.0 ims = [] #各時刻における波動関数の計算 for tn in range(ts, te): #実時間の取得 t = dt * tn #波動関数の計算 ys = psi(x, t).real / np.sqrt(A) * 2+ yE #波動関数の表示 img = plt.plot( x/nm, yE, linestyle='dotted', color="black", linewidth = 1 ) img += plt.plot( x/nm, ys, colors[0], linestyle='solid', linewidth = 3 ) #アニメーションに追加 ims.append( img ) #余白の調整 #plt.subplots_adjust(left=0.15, right=0.90, bottom=0.1, top=0.99) plt.subplots_adjust(left = 0.05, right = 0.95, bottom = 0.10, top = 0.95) #アニメーションの生成 ani = animation.ArtistAnimation(fig, ims, interval=10) #アニメーションの保存 #ani.save("output.html", writer=animation.HTMLWriter()) #ani.save('anim.gif', writer='pillow') #ani.save("output.mp4", writer="ffmpeg", dpi=300) #グラフの表示 plt.show()
【はやくち解説】エルミート多項式の有難みって?【Pythonコピペで量子力学完全攻略マニュアル】
################################################################# ## 規格化エルミート多項式 ################################################################# import math import numpy as np import scipy.integrate as integrate import matplotlib.pyplot as plt #図全体 fig = plt.figure(figsize=(9, 9)) #全体設定 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'] #インデックス最大数 NS = 6 ################################################################# ## エルミート多項式 ################################################################# def NormalizedHermitePolynomial( n, x ): H0 = (1.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) H1 = (4.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) * x if(n==0): return H0 if(n==1): return H1 for m in range(2, n+1): H2 = np.sqrt( 2.0 / m ) * x * H1 - np.sqrt( (m - 1) / m ) * H0 H0 = H1 H1 = H2 return H2 ################################################################# ## 直交性の確認 ################################################################# #被積分関数 def integral_orthonomality(x, n, m): return NormalizedHermitePolynomial( m, x ) * NormalizedHermitePolynomial(n, x) for n in range( NS + 1 ): for m in range( NS + 1 ): #ガウス・ルジャンドル積分 result = integrate.quad( integral_orthonomality, #被積分関数 -10, 10, #積分区間の下端と上端 args=( n, m ) #被積分関数へ渡す引数 ) #ターミナルへ出力 print( "(" + str(n) + ", " + str(m) + ") " + str( result[0]) ) ######################################################################### # グラフの描画(固有関数) ######################################################################### plt.title( u"規格化エルミート多項式", fontsize=20, fontname="Yu Gothic", fontweight=1000) plt.xlabel(r"$x$", fontsize=30) #描画区間数 NX = 500 #描画範囲 x_min = -5 x_max = 5 #座標点配列の生成 x = np.linspace(x_min, x_max, NX) #罫線の描画 plt.grid(which = "major", axis = "x", alpha = 0.6, linestyle = "-", linewidth = 0) plt.grid(which = "major", axis = "y", alpha = 0.6, linestyle = "-", linewidth = 0) #描画範囲を設定 plt.xlim([x_min , x_max]) plt.ylim([-0.5, 10]) plt.yticks(color="None") plt.tick_params(length=0) #グラフの描画 ys = [ 0 ] * (NS + 1) for n in range( NS + 1 ): ys[n] = NormalizedHermitePolynomial( n, x ) + n * 1.5 for n in range( NS + 1 ): plt.plot(x, ys[n], linestyle='solid', linewidth = 5 ) plt.text( 5.1, n * 1.5 - 0.12, r"$H_" + str(n) + "$", fontsize = 25) plt.hlines([n * 1.5], -5, 5, "#000000", linestyles='dotted', linewidth = 1 ) #余白の調整 #plt.subplots_adjust(left=0.15, right=0.90, bottom=0.1, top=0.99) plt.subplots_adjust(left = 0.05, right = 0.92, bottom = 0.12, top = 0.95) #グラフの表示 plt.show()
【第29話】調和ポテンシャル中の電子状態【Pythonコピペで量子力学完全攻略マニュアル】
################################################################# ## 調和ポテンシャル中の電子の固有関数 ################################################################# import math import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation #図全体 fig = plt.figure(figsize=(8, 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.62607015 * 1.0E-34 hbar = h / (2.0 * np.pi) #電子の質量 me = 9.1093837015 * 1.0E-31 #電子ボルト eV = 1.602176634 * 1.0E-19 #ナノメートル nm = 1E-9 #虚数単位 I = 0.0 + 1.0j ################################################################# ## 物理系の設定 ################################################################# #量子井戸の幅(m) omega = 1.0 * 1.0E+15 A = np.sqrt( me * omega / hbar) #エルミート多項式 def HermitePolynomial( n, x ): x0 = 1 x1 = 2.0 * x if(n==0): return x0 if(n==1): return x1 for m in range(2, n+1): x2 = 2.0 * x * x1 - 2.0 * ( m - 1 ) * x0 x0 = x1 x1 = x2 return x2 #固有関数 def verphi(n, x): barX = A * x return ( A**2 / np.pi) **(1.0/4.0) * 1.0 / np.sqrt( 2**n * math.factorial(n) ) * np.exp( - 1.0/2.0 * barX **2 ) * HermitePolynomial( n, barX ) #固有エネルギー(eV) def Energy(n): return (n + 1.0/2.0) * hbar * omega #波動関数 def psi(n, x, t): return verphi(n, x) * np.exp( - I * ( n + 1.0 / 2.0 ) * omega * t ) #状態数 NS = 4 #固有エネルギーの表示(eV) for n in range( NS + 1 ): if( n == 0 ): text = "基底状態" else: text = "第" + str(n) + "励起状態" print( text + ":" + str(Energy( n )/eV) + " [eV]" ) #計算時間の幅 ts = 0 te = 100 #基底状態の周期 T0 = 2.0 * np.pi * hbar / Energy(0) print( "基底状態の周期:" + str(T0) + " [s]" ) ######################################################################### # 波動関数 アニメーション ######################################################################### #時間間隔 dt = T0 / (te - ts + 1) plt.title( u"調和ポテンシャル中の電子の固有関数", fontsize=20, fontname="Yu Gothic", fontweight=1000) plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=20) L = 6 * nm #アニメーション作成用 ims = [] #描画範囲 x_min = -L/2 x_max = L/2 #描画区間数 NX = 500 #座標点配列の生成 x = np.linspace(x_min, x_max, NX) #罫線の描画 plt.grid(which = "major", axis = "x", alpha = 0.2, linestyle = "-", linewidth = 1) plt.grid(which = "major", axis = "y", alpha = 0.2, linestyle = "-", linewidth = 1) #描画範囲を設定 plt.xlim([-3, 3]) plt.ylim([-0.5, 7.5]) NS = 6 #調和ポテンシャルの概形 yE = 1.0 / 2.0 * me * omega**2 * x**2 /eV #各時刻における波動関数の計算 for tn in range(ts, te): #実時間の取得 t = dt * tn img = plt.plot( x/nm, yE, linestyle='dotted', color="black", linewidth = 1 ) #各コマを描画 for n in range( NS + 1 ): #波動関数の計算 y = psi(n, x, t).real / np.sqrt(A) * 0.8 + n + 0.5 #波動関数の表示 img += plt.plot( x/nm, y, colors[n], linestyle='solid', linewidth = 3 ) #アニメーションに追加 ims.append( img ) #E_0~E_6の表示 for n in range( 0, NS + 1 ): plt.text( 3.1, 0.5 + n - 0.12, r"$E_" + str(n) + "$", fontsize = 25) plt.hlines([0.5 + n], -3, 3, "#CCCCCC", linestyles='dotted', linewidth = 1 ) #余白の調整 #plt.subplots_adjust(left=0.15, right=0.90, bottom=0.1, top=0.99) plt.subplots_adjust(left = 0.05, right = 0.92, bottom = 0.10, top = 0.95) #アニメーションの生成 ani = animation.ArtistAnimation(fig, ims, interval=10) #アニメーションの保存 #ani.save("output.html", writer=animation.HTMLWriter()) #ani.save('anim.gif', writer='pillow') ani.save("output.mp4", writer="ffmpeg", dpi=300) #グラフの表示 plt.show()
【第28話】量子ビットの作り方5:1量子ビットの完成!!【量子ドットコンピュータの原理⑤】
############################################################################ # 衝立ありの無限に深い量子井戸に静電場を加えた電子状態 ############################################################################ 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 #全体設定 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.0 * 10**-9 #計算区間 x_min = -L / 2.0 x_max = L / 2.0 #状態数 n_max = 30 #行列の要素数 DIM = n_max + 1 #空間分割数 NX = 500 #空間刻み間隔 nm = 1.0 * 10**-9 #壁の厚さ W = L / 5 #壁の高さの最大値 V_max = 30.0 * eV #電場の強さ Ex_max = 1.0 * 10**8 #電場強度の分割数 NEx = 10 #基底状態と励起状態 N = 2 ###################################### # 衝立なし・静電場無しの解析解 ###################################### #固有関数 def varphi(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): if(abs(x) <= W / 2.0): return e * Ex * x + V_max else: 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 varphi(n1 ,x) * V(x, Ex) * varphi(n2, x) / eV #被積分関数(平均計算用) def average_x(x, a): sum = 0 for n in range(n_max + 1): sum += a[n] * varphi(n, x) return x * sum**2 #固有値・固有ベクトルの初期化 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] * N for n in range( len(phi[nEx]) ): phi[nEx][n] = [0] * (NX + 1) #中心の電場依存性グラフ描画用の配列初期化 averageX = [0] * N for n in range(len(averageX)): averageX[n] = [0] * (NEx + 1) #静電場強度ごとに 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)) ###固有関数の空間分布 for nx in range(NX+1): x = x_min + (x_max - x_min) / NX * nx if(nEx == 0): xs.append( x/nm ) for n in range( len(phi[nEx]) ): for m in range(n_max+1): phi[nEx][n][nx] += vectors[nEx][n][m] * varphi(m, x) #描画用データの整形 phi[nEx][n][nx] = abs(phi[nEx][n][nx])**2 / (1.0 * 10**9) for n in range(len(averageX)): #ガウス・ルジャンドル積分 result = integrate.quad( average_x, #被積分関数 x_min, x_max, #積分区間の下端と上端 args=(vectors[nEx][n]) #被積分関数へ渡す引数 ) #計算結果の取得 averageX[n][nEx] = result[0] * (1.0 * 10**9) ######################################################################### # グラフの描画(エネルギー固有値) ######################################################################### fig1 = plt.figure(figsize=(15, 8)) plt.title( u"無限に深い量子井戸(衝立あり)に静電場を加えたときの固有エネルギー", fontsize=20, fontname="Yu Gothic", fontweight=1000) plt.xlabel( u"静電場の強さ" + r"$\times 10^7 \,[{\rm V/m}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000 ) plt.ylabel( r"$ E\,[{\rm eV}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000 ) #余白の調整 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, 10]) #x軸 exs = range( NEx + 1) #y軸 En_0 = [] En_1 = [] for nEx in range(NEx + 1): En_0.append( eigenvalues[nEx][0] ) En_1.append( eigenvalues[nEx][1] ) #print( str(nV) + " " + str( eigenvalues[nV][0] ) + " " + str( eigenvalues[nV][1] )) #基底状態と第1励起状態のグラフを描画 plt.plot(exs, En_0, marker="o", markersize=15, linewidth = 5) plt.plot(exs, En_1, marker="o", markersize=15, linewidth = 5) #グラフの表示 plt.show() ######################################################################### # グラフの描画(基底状態) ######################################################################### fig2 = plt.figure(figsize=(15, 8)) 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.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 ) #衝立の概形 plt.vlines([-0.1], -1, 3, "black", linestyles='solid', linewidth = 10 ) plt.vlines([0.1], -1, 3, "black", linestyles='solid', linewidth = 10 ) plt.hlines([3], -0.106, 0.106, "black", linestyles='solid', linewidth = 10 ) #描画範囲を設定 plt.xlim([-0.6, 0.6]) plt.ylim([0, 5.0]) #各 for nEx in range(NEx + 1): plt.plot(xs, phi[nEx][0] , linewidth = 5) #グラフの表示 plt.show() ######################################################################### # グラフの描画(第1励起状態) ######################################################################### fig3 = plt.figure(figsize=(15, 8)) plt.title( u"無限に深い量子井戸(衝立あり)に静電場を加えたときの第1励起状態", 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, 5.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 ) #衝立の概形 plt.vlines([-0.1], -1, 3, "black", linestyles='solid', linewidth = 10 ) plt.vlines([0.1], -1, 3, "black", linestyles='solid', linewidth = 10 ) plt.hlines([3], -0.106, 0.106, "black", linestyles='solid', linewidth = 10 ) for nEx in range(NEx + 1): plt.plot(xs, phi[nEx][1] , linewidth = 5) #グラフの表示 plt.show()
【第27話】量子ビットの作り方4:量子井戸の真ん中に衝立を配置したときの電子状態【量子ドットコンピュータの原理④】
############################################################################ # 衝立ありの無限に深い量子井戸中の電子状態 ############################################################################ 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 #全体設定 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.0 * 10**-9 #計算区間 x_min = -L / 2.0 x_max = L / 2.0 #状態数 n_max = 30 #行列の要素数 DIM = n_max + 1 #空間分割数 NX = 500 #空間刻み間隔 nm = 1.0 * 10**-9 #壁の厚さ W = L / 5 #壁の高さの最大値 V_max = 30.0 * eV #壁の高さの分割数 NV = 30 ###################################### # 衝立なしの解析解 ###################################### #固有関数 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, V0): if(abs(x) <= W / 2.0): return V0 else: return 0 #固有エネルギー def Energy(n): kn = np.pi * (n + 1) / L return hbar * hbar * kn * kn / (2.0 * me) ###################################### # 衝立ごとの固有関数の計算 ###################################### #被積分関数(行列要素計算用) def integral_matrixElement(x, n1, n2, V0): return verphi(n1 ,x) * V(x, V0) * verphi(n2, x) / eV #被積分関数(平均計算用) def average_x(x, a): sum = 0 for n in range(n_max + 1): sum += a[n] * verphi(n, x) return x * sum**2 #固有値・固有ベクトルの初期化 eigenvalues = [0] * (NV + 1) vectors = [0] * (NV + 1) for nV in range(NV + 1): eigenvalues[nV] = [] vectors[nV] = [] #存在確率分布グラフ描画用の配列初期化 xs = [] phi = [0] * (NV + 1) for nV in range(NV + 1): phi[nV] = [0] * 2 for n in range( len(phi[nV]) ): phi[nV][n] = [0] * (NX + 1) #中心の電場依存性グラフ描画用の配列初期化 averageX = [0] * 2 for n in range(len(averageX)): averageX[n] = [0] * (NV + 1) #衝立の高さごとに for nV in range(NV + 1): #if(nV == 0): continue print("衝立の高さ:" + str( nV * 100 / NV ) + "%") #衝立を設定 V0 = V_max / NV * nV #エルミート行列(リスト) 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, V0) #被積分関数へ渡す引数 ) real = result[0] #print("(" + str(n1) + "," + str(n2) + ") " + str(real) ) 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[nV] = eig[ index ] #転置行列 vec = vec.T #固有ベクトルの並び替え vectors[nV] = vec[ index ] ### 検算:MA-EA=0 ? sum = 0 for i in range(DIM): v = matrix @ vectors[nV][i] - eigenvalues[nV][i] * vectors[nV][i] for j in range(DIM): sum += abs(v[j])**2 print("|MA-EA| =" + str(sum)) ###固有関数の空間分布 for nx in range(NX+1): x = x_min + (x_max - x_min) / NX * nx if(nV == 0): xs.append( x/nm ) for n in range( len(phi[nV]) ): for m in range(n_max+1): phi[nV][n][nx] += vectors[nV][n][m] * verphi(m, x) #描画用データの整形 phi[nV][n][nx] = abs(phi[nV][n][nx])**2 / (1.0 * 10**9) for n in range(len(averageX)): #ガウス・ルジャンドル積分 result = integrate.quad( average_x, #被積分関数 x_min, x_max, #積分区間の下端と上端 args=(vectors[nV][n]) #被積分関数へ渡す引数 ) #計算結果の取得 averageX[n][nV] = result[0] * (1.0 * 10**9) ######################################################################### # グラフの描画(エネルギー固有値) ######################################################################### fig1 = plt.figure(figsize=(15, 8)) plt.title( u"無限に深い量子井戸に衝立を加えたときの固有エネルギー", fontsize=20, fontname="Yu Gothic", fontweight=1000) plt.xlabel( u"衝立の高さ" + r"$ V\,[{\rm eV}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000 ) plt.ylabel( r"$ E\,[{\rm eV}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000 ) #余白の調整 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, 30]) #x軸 exs = range( NV + 1 ) #y軸 En_0 = [] En_1 = [] for nV in range(NV + 1): En_0.append( eigenvalues[nV][0] ) En_1.append( eigenvalues[nV][1] ) print( str(nV) + " " + str( eigenvalues[nV][0] ) + " " + str( eigenvalues[nV][1] ) ) #基底状態と第1励起状態のグラフを描画 plt.plot(exs, En_0, marker="o", markersize=15, linewidth = 5) plt.plot(exs, En_1, marker="o", markersize=15, linewidth = 5) #グラフの表示 plt.show() ######################################################################### # グラフの描画(基底状態) ######################################################################### fig2 = plt.figure(figsize=(15, 8)) 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, 3.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 ) #衝立の概形 plt.vlines([-0.1], -1, 2, "black", linestyles='solid', linewidth = 10 ) plt.vlines([0.1], -1, 2, "black", linestyles='solid', linewidth = 10 ) plt.hlines([2], -0.106, 0.106, "black", linestyles='solid', linewidth = 10 ) for nV in range(NV + 1): plt.plot(xs, phi[nV][0], linewidth = 5) #グラフの表示 plt.show() ######################################################################### # グラフの描画(第1励起状態) ######################################################################### fig3 = plt.figure(figsize=(15, 8)) plt.title( u"無限に深い量子井戸に静電場を加えたときの第1励起状態", 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, 3.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 ) #衝立の概形 plt.vlines([-0.1], -1, 2, "black", linestyles='solid', linewidth = 10 ) plt.vlines([0.1], -1, 2, "black", linestyles='solid', linewidth = 10 ) plt.hlines([2], -0.106, 0.106, "black", linestyles='solid', linewidth = 10 ) for nV in range(NV + 1): plt.plot(xs, phi[nV][1] , linewidth = 5) #グラフの表示 plt.show()
【第26話】【結果は渋いです】2次元電子波束を斜めに薄膜へ(一応トンネル効果)
################################################################# ## 薄膜へ照射したガウス波束(トンネル効果) ################################################################# import os import numpy as np import matplotlib.pyplot as plt import matplotlib as mat import matplotlib.animation as animation from multiprocessing import Pool flag_calculate = False flag_animation = True #図全体 aspect = 0.9 fig = plt.figure(figsize=(6, 6 * aspect)) #全体設定 plt.rcParams['font.family'] = 'Times New Roman' #フォント plt.rcParams['font.size'] = 16 #フォントサイズ plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント #カラーリストの取得 colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] #並列数 ParallelNumber = 6 #フォルダ指定 dir_path = "output30/" #フォルダ生成 os.makedirs(dir_path, exist_ok = True) ################################################################# ## 物理定数 ################################################################# #プランク定数 h = 6.62607015 * 1.0E-34 hbar = h / (2.0 * np.pi) #電子の質量 me = 9.1093837015 * 1.0E-31 #電子ボルト eV = 1.602176634 * 1.0E-19 #虚数単位 I = 0.0 + 1.0j ###################################### # 物理系の設定 ###################################### #波束の中心エネルギー E0 = 10.0 * eV #入射角度 theta0 = 30.0 / 180 * np.pi #重ね合わせの数 NK = 100 #ガウス分布の幅 sigma = 0.1 / ( 1.25 * 10**-9 ) #波数の間隔 dk = 10.0 * sigma / NK #ガウス分布の中心波数 k0 = np.sqrt( 2.0 * me * E0 / hbar**2) kx0 = k0 * np.cos(theta0) ky0 = k0 * np.sin(theta0) #空間刻み間隔 nm = 1E-9 #空間刻み数 NX = 150 NY = 150 #壁の高さ V1 = 0.0 * eV V2 = 10.5 * eV V3 = 0.0 * eV #壁の幅 d = 1.0 * nm #計算時間の幅 ts = -200 te = 200 #時間間隔 dt = 1.0 * 10**-16 #空間幅 x_min = -30.0 * nm x_max = 30.0 * nm y_min = x_min y_max = x_max ###################################### # 反射係数と透過係数を計算する関数定義 ###################################### def ReflectionCoefficient( MM ): return - MM[1][0] / MM[1][1] def TransmissionCoefficient(MM): return (MM[0][0] * MM[1][1] - MM[0][1] * MM[1][0] ) / MM[1][1] def calculateCoefficient( k1, k2, k3, cos1, cos2, cos3, d): # 転送行列の定義 M21 = np.zeros((2,2), dtype=np.complex) T2 = np.zeros((2,2), dtype=np.complex) M32 = np.zeros((2,2), dtype=np.complex) # 境界 M21[0][0] = (1.0 + 0.0j + ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0 M21[0][1] = (1.0 + 0.0j - ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0 M21[1][0] = (1.0 + 0.0j - ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0 M21[1][1] = (1.0 + 0.0j + ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0 M32[0][0] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0 M32[0][1] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0 M32[1][0] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0 M32[1][1] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0 T2[0][0] = np.exp( I * k2 * cos2 * d ) T2[0][1] = 0 T2[1][0] = 0 T2[1][1] = np.exp( - I * k2 * cos2 * d ) #転送行列の計算 M31 = M32 @ T2 @ M21 #反射係数と透過係数の計算 rc = ReflectionCoefficient(M31) tc = TransmissionCoefficient(M31) return rc, tc, M21 #透過係数と反射係数と転送行列の配列 rcs = np.zeros((NK, NK)) tcs = np.zeros((NK, NK)) rcs = np.array(rcs, dtype=complex) tcs = np.array(tcs, dtype=complex) M21s = [ 0 ] * NK for n1 in range(NK): M21s[ n1 ] = [ np.zeros((2,2), dtype=np.complex) ] * NK for kxn in range(NK): for kyn in range(NK): kx = kx0 + dk * (kxn - NK/2) ky = ky0 + dk * (kyn - NK/2) k = np.sqrt( kx**2 + ky**2 ) 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 )) k3 = np.sqrt(np.complex128( k**2 - 2.0 * me * V3 / hbar**2 )) sin = ky / k cos = kx / k sin1 = sin * k / k1 sin2 = sin * k / k2 sin3 = sin * k / k3 cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 )) cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 )) cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 )) rcs[kxn][kyn], tcs[kxn][kyn], M21s[kxn][kyn] = calculateCoefficient( k1, k2, k3, cos1, cos2, cos3, d) ###################################### # 波動関数の計算 ###################################### #ガウス波束の初期位置 x0 = 0 #-20.0 * nm * np.cos(theta0) y0 = 0 #-20.0 * nm * np.sin(theta0) #ガウス波束の値を計算する関数 def Psi( x, y, t ): #波動関数値の初期化 psi = 0 + 0j #各波数ごとの寄与を足し合わせる for kxn in range(NK): for kyn in range(NK): #各波数を取得 kx = kx0 + dk * (kxn - NK/2) ky = ky0 + dk * (kyn - NK/2) k = np.sqrt( kx**2 + ky**2 ) 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 )) k3 = np.sqrt(np.complex128( k**2 - 2.0 * me * V3 / hbar**2 )) sin = ky / k cos = kx / k sin1 = sin * k / k1 sin2 = sin * k / k2 sin3 = sin * k / k3 cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 )) cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 )) cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 )) k1x = k1 * cos1 k1y = k1 * sin1 k2x = k2 * cos2 k2y = k2 * sin2 k3x = k3 * cos3 k3y = k3 * sin3 #波数から各振動数を取得 omega = hbar / (2.0 * me) * ( kx**2 + ky**2 ) W = np.exp( - I * k1x * x0 - I * k1y * y0) * np.exp( -( (kx - kx0) / (2.0 * sigma) )**2 -( (ky - ky0) / (2.0 * sigma) )**2 ) if( x < 0 ): psi += np.exp( I * ( k1x * x + k1y * y ) - I * omega * t ) * W psi += rcs[kxn][kyn] * np.exp( I * ( - k1x * x + k1y * y ) - I * omega * t ) * W elif( x < d ): A2_p = M21s[kxn][kyn][0][0] + M21s[kxn][kyn][0][1] * rcs[kxn][kyn] A2_m = M21s[kxn][kyn][1][0] + M21s[kxn][kyn][1][1] * rcs[kxn][kyn] psi += A2_p * np.exp( I * ( k2x * x + k2y * y ) - I * omega * t ) * W psi += A2_m * np.exp( I * ( - k2x * x + k2y * y ) - I * omega * t ) * W elif( d <= x ): psi += tcs[kxn][kyn] * np.exp( I * ( k3x * (x - d) + k3y * y ) - I * omega * t ) * W return psi #基準となる振幅を取得 psi_abs = abs(Psi( x0, 0, 0 )) ###################################### # 空間分布の計算 ###################################### #実時間 t = 0 #x軸・y軸のメッシュ生成 xs = np.linspace(x_min, x_max, NX) ys = np.linspace(y_min, y_max, NY) X, Y = np.meshgrid(xs, ys) #並列用ラッパー関数 def wrap_timeEvolution( data ): return timeEvolution( *data ) #時間発展の計算 def timeEvolution( filepath , tn ): #実時間 t = ( ts + tn ) * dt Z = np.sin( np.complex128(X/nm) ) for ix in range(NX): for iy in range(NY): x = xs[ix] y = ys[iy] Z[iy][ix] = Psi( x, y, t ) / psi_abs #ファイルへ保存 np.savetxt( filepath, abs(Z)) #アニメーションフレーム作成用 def update( tn ): #グラフ消去 fig.clear() #両軸位置の設定 axes = fig.add_axes( [0.07, 0.1/aspect, 0.78, 0.78/aspect]) #カラーバー位置の設定 bar_axes = fig.add_axes([0.87, 0.1/aspect, 0.05, 0.78/aspect]) #データ読み込み filepath = dir_path + str(tn + ts) + ".txt" Z = np.loadtxt( filepath ) #透過領域の振幅を定数倍 for ny in range(len(X)): for nx in range(len(X[ny])): if(X[ny][nx] > 1.0 * nm ): Z[ny][nx] = Z[ny][nx] * 1000 Z = Z[:-1, :-1] #2次元カラーグラフ描画 img = axes.pcolormesh( X/nm, Y/nm, Z, vmin=0.0, vmax=1.0, cmap=cmap) #カラーバーの追加 fig.colorbar(img, cax=bar_axes) ####### メイン関数として実行 ######## if __name__ == "__main__": #計算開始 if(flag_calculate): list = [] for tn in range( te - ts + 1 ): filepath = dir_path + str(tn + ts) + ".txt" list.append( (filepath, tn) ) with Pool(processes=ParallelNumber) as p: p.map(wrap_timeEvolution, list) #アニメーション作成 if( flag_animation ): #カラーマップの配色配列 color_list = [] color_list.append( [ 0, "black" ] ) color_list.append( [ 1, "white" ] ) #カラーマップ用オブジェクトの生成 cmap = mat.colors.LinearSegmentedColormap.from_list('cmap', color_list) #アニメーションの設定 ani = animation.FuncAnimation( fig, update, range(te - ts + 1), interval=10) #アニメーションの保存 #ani.save("output.html", writer=animation.HTMLWriter()) #ani.save('anim.gif', writer='pillow') ani.save("output.mp4", writer="ffmpeg", dpi=300) #グラフの表示 #plt.show()
【第25話】量子ビットの作り方3:静電場の与えたときの電子状態【量子ドットコンピュータの原理③】
############################################################################ # 無限に深い量子井戸中の電子に静電場を加えたときの固有状態(ステップ3) ############################################################################ 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()
【第24話】静電場の与え方【量子ドットコンピュータの原理②】
############################################################################ # 無限に深い量子井戸中の電子に静電場を加えたときの固有状態(ステップ1) ############################################################################ import math import cmath import matplotlib.pyplot as plt import matplotlib.animation as animation import scipy.integrate as integrate #図全体 fig1 = plt.figure(figsize=(10, 6)) #全体設定 plt.rcParams['font.family'] = 'Times New Roman' #フォント plt.rcParams['font.size'] = 12 #フォントサイズ #複素数 I = 0.0 + 1.0j ###################################### # 物理定数 ###################################### #プランク定数 h = 6.6260896 * 10**-34 hbar = h / (2.0 * math.pi) #電子の質量 me = 9.10938215 * 10**-31 #電子ボルト eV = 1.60217733 * 10**-19 #電気素量 e = 1.60217733 * 10**-19 ###################################### # 物理系の設定 ###################################### #量子井戸の幅 L = 1.0 * 10**-9 #計算区間 x_min = -L / 2.0 x_max = L / 2.0 #状態数 n_max = 10 #電場の強さ Ex = 1.0*10**10 #固有関数 def verphi(n, x): kn = math.pi * (n + 1) / L return math.sqrt(2.0 / L) * math.sin(kn * (x + L / 2.0)) #ポテンシャル項 def V(x, Ex): return(e * Ex * x) #被積分関数 def integral_matrixElement(x, n1, n2, Ex): return verphi(n1 ,x) * V(x, Ex) * verphi(n2, x) / eV ###行列要素の計算 for n1 in range(n_max + 1): for n2 in range(n_max + 1): #ガウス・ルジャンドル積分 result = integrate.quad( integral_matrixElement, #被積分関数 x_min, x_max, #積分区間の下端と上端 args=(n1,n2, Ex) #被積分関数へ渡す引数 ) real = result[0] imag = 0 #ターミナルへ出力 print( "(" + str(n1) + ", " + str(n2) + ") " + str(real))
【第23話】正規直交系の有難みって?【Pythonコピペで量子力学完全攻略マニュアル】
################################################################# ## 無限に深い量子井戸の固有関数の直交性の確認 ################################################################# import numpy as np import scipy.integrate as integrate ################################################################# ## 物理定数 ################################################################# #プランク定数 h = 6.62607015 * 1.0E-34 hbar = h / (2.0 * np.pi) #電子の質量 me = 9.1093837015 * 1.0E-31 #電子ボルト eV = 1.602176634 * 1.0E-19 ################################################################# ## 物理系の設定 ################################################################# #量子井戸の幅(m) L = 1.0 * 10**-9 #計算区間 x_min = -L / 2.0 x_max = L / 2.0 #固有関数 def verphi(n, x): kn = np.pi * (n + 1) / L return np.sqrt(2.0 / L) * np.sin(kn * (x + L / 2.0)) #被積分関数 def integral_orthonomality(x, n, m): return verphi(m, x) * verphi(n, x) #状態数 NS = 3 for n in range( NS + 1 ): for m in range( NS + 1 ): #ガウス・ルジャンドル積分 result = integrate.quad( integral_orthonomality, #被積分関数 x_min, x_max, #積分区間の下端と上端 args=( n, m ) #被積分関数へ渡す引数 ) #ターミナルへ出力 print( "(" + str(n) + ", " + str(m) + ") " + str( result[0]) )