【第16話】周期ポテンシャル（クローニッヒ・ペニーモデル）の波動関数アニメーション【Pythonコピペで量子力学完全攻略マニュアル】

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

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

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

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

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

    #転送行列
    return M , M21T1, M12T2


#################################################################
## 転送行列の固有値
#################################################################
#バンド端エネルギー
BandEdgeEs = [0.158706, 0.266136, 0.615037, 0.849982]
#バンド中心エネルギー
BandCenterEs = [0.205879, 0.7251]
#エネルギー
E = ( BandCenterEs[1] )* eV 
#E = ( BandEdgeEs[3] )* eV 

k = np.sqrt( 2.0 * me * E / hbar ** 2)
M, M21T1, M12T2 = TransferMatrix( k )
print("------------------------")
eig, A = LA.eig(M)
print( eig )
print( A )
print("------------------------")
print( abs( M @ A -  A * eig ) )

LN = 5
LNs = np.linspace(0, LN, LN + 1)

Ap = [ 0 + 0j ] * ( 2 * LN + 1 )
Am = [ 0 + 0j ] * ( 2 * LN + 1 )



#################################################################
## 波動関数の描画
#################################################################

#描画範囲
x_min = 0.0 * nm
x_max = l * LN
#角振動数
omega = E/hbar
#１周期
T = 2.0 * np.pi / omega
#時間間隔
#dt = 0.2 * 10**-16
dt = T / 100
#時間刻み数
NT = 100

#アニメーション作成用
ims=[]

#描画区間数
NX = 100

xs = [0] * LN * 2
xK =np.linspace( x_min,  x_max, 500 )

#座標点配列の生成
for n in range(LN):
    xs[ 2 * n ] = np.linspace( l * n,  l * n + d1, NX )
    xs[ 2 * n + 1 ] = np.linspace( l * n + d1,  l * n + l, NX )

Ap[0] = A[0][0]
Am[0] = A[1][0]
for n in range(LN):
    Ap[2*n+1] = M21T1[0][0] * Ap[2*n] + M21T1[0][1] * Am[2*n]
    Am[2*n+1] = M21T1[1][0] * Ap[2*n] + M21T1[1][1] * Am[2*n]
    Ap[2*n+2] = M12T2[0][0] * Ap[2*n+1] + M12T2[0][1] * Am[2*n+1]
    Am[2*n+2] = M12T2[1][0] * Ap[2*n+1] + M12T2[1][1] * Am[2*n+1]

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 ))

phi_p = [ 0 + 0j ] * ( 2 * LN )
phi_m = [ 0 + 0j ] * ( 2 * LN )
phi = [ 0 + 0j ] * ( 2 * LN )

#各時刻における計算
for tn in range(NT):
    t = dt * tn

    #波動関数の計
    for n in range(LN):
        #領域１
        phi_p[2*n] = Ap[2*n] * np.exp(   I * k1 * ( xs[2*n] - l * n ) - I * omega * t )
        phi_m[2*n] = Am[2*n] * np.exp( - I * k1 * ( xs[2*n] - l * n ) - I * omega * t )
        phi[2*n] = phi_p[2*n] + phi_m[2*n]
        #領域２
        phi_p[2*n+1] = Ap[2*n+1] * np.exp(   I * k2 * ( xs[2*n+1] - ( l * n + d1 ) ) - I * omega * t )
        phi_m[2*n+1] = Am[2*n+1] * np.exp( - I * k2 * ( xs[2*n+1] - ( l * n + d1 ) ) - I * omega * t )
        phi[2*n+1] = phi_p[2*n+1] + phi_m[2*n+1]

        if( n == 0 ):
            img  = plt.plot(xs[2*n]/nm, phi[2*n].real, colors[0], linestyle='solid', linewidth = 5.0)
        else:
            img  += plt.plot(xs[2*n]/nm, phi[2*n].real, colors[0], linestyle='solid', linewidth = 5.0)

        img += plt.plot(xs[2*n]/nm, phi[2*n].imag, colors[1], linestyle='solid', linewidth = 5.0)
        img += plt.plot(xs[2*n+1]/nm, phi[2*n+1].real, colors[0], linestyle='solid', linewidth = 5.0)
        img += plt.plot(xs[2*n+1]/nm, phi[2*n+1].imag, colors[1], linestyle='solid', linewidth = 5.0)


    ims.append( img )


#グラフの描画（波動関数）
plt.title( u"クローニッヒ・ペニーモデルの波動関数（" + r"$l=2.5[{\rm nm}]$" + " , " + r"$d=1.0[{\rm nm}]$"+ " , " + r"$V=-10.0[{\rm eV}]$" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ \psi(x, t) $", 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)


#境界の描画
for n in range(LN):
    plt.hlines([0], n*l/nm,  (n*l+d1)/nm, "black", linestyles='solid', linewidth = 5 )
    plt.vlines([(n*l+d1)/nm], -1, 0, "black", linestyles='solid', linewidth = 5 )
    plt.hlines([-1], (n*l+d1)/nm,  (n*l + l)/nm, "black", linestyles='solid', linewidth = 5 )
    plt.vlines([(n*l + l)/nm], -1, 0, "black", linestyles='solid', linewidth = 5 )



#描画範囲を設定
plt.xlim([x_min/nm, x_max/nm])
plt.ylim([-1.5, 1.5])

#アニメーションの生成
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()