【第32話】コヒーレント状態ってどんな状態?【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## [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]) )