【第22話】【はやくち解説】トンネル効果ってなに？【Pythonコピペで量子力学完全攻略マニュアル】

※計算時間がかなり掛かります。

#################################################################
## 薄膜へ照射したガウス波束（トンネル効果）
#################################################################
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 = True
flag_animation = True

#図全体
aspect = 0.9
fig = plt.figure(figsize=(8, 4 * 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 = 12
#フォルダ指定
dir_path = "output/"
#フォルダ生成
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 = 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 = 200
NY = 100
#壁の高さ
V1 = 0.0 * eV
V2 = 10.5 * eV
V3 = 0.0 * eV
#壁の幅
d = 1.0 * nm

#計算時間の幅
ts = 0
te = 500
#時間間隔
dt = 1.0 * 10**-16

#空間幅
x_min = -50.0 * nm
x_max =  50.0 * nm
y_min = x_min/2
y_max = x_max/2

######################################
#　反射係数と透過係数を計算する関数定義
######################################
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 = -35.0 * nm

#ガウス波束の値を計算する関数
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 ) * 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] * 5

    Z = Z[:-1, :-1]

    #２次元カラーグラフ描画
    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()