QD3D.py ダウンロード/コピー

QD3D.py をダウンロード

QD3D.py
QD3D.py
  1"""
  2量子ドット (Quantum Dot) におけるエネルギー準位計算および関連関数のプロットスクリプト
  3
  4概要:
  5    本スクリプトは、球対称ポテンシャル井戸モデルを用いた量子ドットのエネルギー準位を計算し、
  6    その結果を表示する機能、および関連する物理関数（球Bessel関数や水素原子の動径波動関数）を
  7    プロットする機能を提供する。
  8
  9詳細説明:
 10    cal モードでは、指定された有効質量と半径を持つ量子ドットの電子のエネルギー準位を、
 11    球Bessel関数の零点を用いて計算します。
 12    plot モードでは、球Bessel関数のグラフを描画し、その零点を視覚化します。
 13    plotH モードでは、水素原子の動径波動関数を描画します。
 14    コマンドライン引数により実行モードを選択できます。
 15
 16関連リンク:
 17    3DQD_usage
 18    量子力学I/球対称井戸型ポテンシャル: https://dora.bk.tsukuba.ac.jp/~takeuchi/?%E9%87%8F%E5%AD%90%E5%8A%9B%E5%AD%A6%E2%85%A0%2F%E7%90%83%E5%AF%BE%E7%A7%B0%E4%BA%95%E6%88%B8%E5%9E%8B%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB
 19    量子力学I/3次元調和振動子: https://dora.bk.tsukuba.ac.jp/~takeuchi/?%E9%87%8F%E5%AD%90%E5%8A%9B%E5%AD%A6%E2%85%A0%2F%EF%BC%93%E6%AC%A1%E5%85%83%E8%AA%BF%E5%92%8C%E6%8C%AF%E5%8B%95%E5%AD%90#fe164113
 20"""
 21
 22import sys
 23import numpy as np
 24from numpy import sqrt
 25from math import factorial
 26from scipy.special import spherical_jn, genlaguerre
 27import scipy.optimize as opt
 28import matplotlib.pyplot as plt
 29
 30
 31# 定数
 32url = "https://dora.bk.tsukuba.ac.jp/~takeuchi/?%E9%87%8F%E5%AD%90%E5%8A%9B%E5%AD%A6%E2%85%A0%2F%E7%90%83%E5%AF%BE%E7%A7%B0%E4%BA%95%E6%88%B8%E5%9E%8B%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB"
 33me = 9.10938356e-31   # kg
 34hbar = 1.0545718e-34    # J·s
 35eV   = 1.60218e-19      # J/eV
 36
 37lstr = ['s', 'p', 'd', 'f', 'g', 'h', 'l']
 38
 39
 40meff = 0.067
 41R = 5e-9 # m
 42Z = 1
 43
 44mode = "cal"
 45
 46
 47if __name__ == "__main__":
 48    nargs = len(sys.argv)
 49    if nargs >= 2: mode = sys.argv[1]
 50    if nargs >= 3: R = 1.0e-9 * float(sys.argv[2]) # m
 51    if nargs >= 4: meff = float(sys.argv[3])
 52    if nargs >= 5: Z = float(sys.argv[4])
 53
 54
 55def get_zeros(func, xmin = 0.0, xmax = 10.0, dx = 0.1, eps = 1.0e-10, nmaxiter = 50, h = 1.0e-10, dump = 1.0, print_level = 0):
 56    """
 57    概要:
 58        関数の零点を計算します。
 59    詳細説明:
 60        与えられた関数 func の指定された区間 [xmin, xmax] 内における零点を探索します。
 61        零点の探索には、関数値の符号反転を検出した後、ニュートン法に似た数値的な手法を使用します。
 62        各ステップで接線近似を用いて次点の推定を行い、指定された収束条件 eps または最大反復回数 nmaxiter に達するまで繰り返します。
 63    引数:
 64        :param func: 零点を探す関数。引数を1つ取るcallableオブジェクト。
 65        :type func: callable
 66        :param xmin: 探索範囲の最小値。
 67        :type xmin: float
 68        :param xmax: 探索範囲の最大値。
 69        :type xmax: float
 70        :param dx: 初期探索におけるステップサイズ。
 71        :type dx: float
 72        :param eps: 零点探索の収束判定に用いる許容誤差。
 73        :type eps: float
 74        :param nmaxiter: ニュートン法系反復の最大回数。
 75        :type nmaxiter: int
 76        :param h: 数値微分の計算に使用する微小な差分。
 77        :type h: float
 78        :param dump: ニュートン法におけるステップサイズの調整係数（ダンプファクター）。1.0は標準的なステップサイズ。
 79        :type dump: float
 80        :param print_level: デバッグ情報の出力レベル。0で非表示、1で表示。
 81        :type print_level: int
 82    戻り値:
 83        :returns: 見つかった零点のリスト。
 84        :rtype: list[float]
 85    """
 86
 87    zeros = []
 88    nnode = 0
 89    prev_r = 0.0
 90    prev_f = 1.0
 91    for r in np.arange(xmin, xmax + dx, dx):
 92        f = func(r)
 93
 94        if f == 0.0:
 95            zeros.append(r)
 96        elif r > 0.0 and f * prev_f < 0.0:
 97            r0 = prev_r
 98            r1 = r
 99            for i in range(nmaxiter):
100                if print_level:
101                    print(f" iter#{i}/{nmaxiter}: nnode={nnode} dr={r1:8.4g} - {r0:8.4g} = {r1-r0:8.4g}")
102                if abs(r1 - r0) < eps: 
103                    zeros.append(r1)
104                    break
105
106                r0 = r1
107                f1 = func(r1)
108
109                rh = r1 + h
110                fh = func(rh)
111                fdiff = (fh - f1) / (rh - r1)
112
113                r1 = r1 - f1 / fdiff / dump
114
115            nnode += 1
116
117        prev_r = r
118        prev_f = f
119
120    return zeros
121
122def get_bessel_zeros(l, xmin = 0.0, xmax = 10.0, dx = 0.1, print_level = 0):
123    """
124    概要:
125        球Bessel関数 j_l(x) の零点を計算します。
126    詳細説明:
127        指定された次数 l の球Bessel関数 scipy.special.spherical_jn(l, x) を対象として、
128        get_zeros 関数を用いて零点を探索します。
129    引数:
130        :param l: 球Bessel関数の次数 (軌道角運動量量子数)。
131        :type l: int
132        :param xmin: 探索範囲の最小値。
133        :type xmin: float
134        :param xmax: 探索範囲の最大値。
135        :type xmax: float
136        :param dx: 初期探索におけるステップサイズ。
137        :type dx: float
138        :param print_level: デバッグ情報の出力レベル。get_zeros 関数に渡されます。
139        :type print_level: int
140    戻り値:
141        :returns: 球Bessel関数の零点のリスト。
142        :rtype: list[float]
143    """
144
145    return get_zeros(lambda x: spherical_jn(l, x), xmin, xmax, dx, print_level = print_level)
146
147
148def energy_level(meff, R, n, l, zeros):
149    """
150    概要:
151        球対称量子井戸（量子ドット）のエネルギー準位を計算します。
152    詳細説明:
153        無限に深い球対称ポテンシャル井戸モデルに基づき、与えられた量子数 n と l に対応する
154        エネルギー準位を計算します。エネルギーは、球Bessel関数の零点 alpha_nl を用いて
155        E_nl = (hbar^2 * k_nl^2) / (2 * m_eff * m_e) の式で求められます。
156        ここで k_nl = alpha_nl / R です。
157    引数:
158        :param meff: 電子の有効質量 (自由電子質量 me に対する比率)。
159        :type meff: float
160        :param R: 量子ドットの半径 (m)。
161        :type R: float
162        :param n: 主量子数 (1から始まる)。球Bessel関数のn番目の零点に対応。
163        :type n: int
164        :param l: 軌道角運動量量子数。
165        :type l: int
166        :param zeros: 球Bessel関数の零点のリスト。zeros[l][n-1] の形式で零点にアクセスします。
167        :type zeros: list[list[float]]
168    戻り値:
169        :returns: 計算されたエネルギー準位 (eV) と対応する球Bessel関数の零点 alpha_nl のタプル。零点が見つからない場合は (None, None) を返します。
170        :rtype: tuple[float, float] or tuple[None, None]
171    """
172    # 球Bessel関数のn番目の零点αnlを求め、波数ベクトルknlを計算
173    if len(zeros[l]) <= n - 1: return None, None
174    aR = zeros[l][n-1]
175
176    knl = aR / R
177    Enl = hbar * hbar * knl * knl / 2.0 / meff / me
178    Enl_eV = Enl / eV
179    
180    return Enl_eV, aR
181
182
183def cal():
184    """
185    概要:
186        量子球のエネルギー準位を計算し、結果を標準出力に表示します。
187    詳細説明:
188        設定された有効質量 meff と半径 R を持つ量子球（量子ドット）に対し、
189        指定された範囲の主量子数 n と軌道角運動量量子数 l について、
190        そのエネルギー準位を計算します。
191        まず get_bessel_zeros を用いて球Bessel関数の零点を取得し、
192        次に energy_level を用いて各準位のエネルギーを計算します。
193        計算結果はエネルギーの低い順にソートされ、各準位の情報（量子数、エネルギー、零点値）が出力されます。
194        
195    戻り値:
196        :returns: なし
197        :rtype: None
198    """
199    print()
200    print("Energy levels for quantum sphere")
201    print(f"  effective mass: {meff} me")
202    print(f"  radius: {R*1.0e9} nm")
203
204    nmax = 5
205    zeros = []
206    print("Zero points:")
207    for l in range(0, nmax):
208        zero_points = get_bessel_zeros(l, xmin = 0.1, xmax = 10.0, dx = 0.1, print_level = 0)
209        zeros.append(zero_points)
210        print(f"l={l}:", [float(v) for v in zero_points])
211
212    level_list = []
213    for n in range(1, nmax):
214       for l in range(0, nmax):
215            Enl, R_zero = energy_level(meff, R, n, l, zeros)
216            if Enl is None: continue
217
218            level_list.append({"E": Enl, "label": f"{n+l}{lstr[l]}", "n": n, "l": l, 
219                                "R0": R_zero})
220
221    print(f"l n l+n  E(orb)")
222    for d in sorted(level_list, key=lambda x: x["E"]):
223        E     = d["E"]
224        label = d["label"]
225        n     = d["n"]
226        l     = d["l"]
227        R0 = d["R0"]
228        print(f"{l} {n} {l+n}    E({label})={E:12.8g} eV  alpha_{n}_{l}={R0:12.8g}  alpha_{n}_{l}^2={R0*R0:12.8g}")
229    print(f"see {url} for the definitions of quantum numbers")
230
231def plot_spherical_bessel(lmax, rmax, rmesh = 500):
232    """
233    概要:
234        球Bessel関数 j_l(x) をプロットします。
235    詳細説明:
236        0から lmax までの次数 l について、球Bessel関数 j_l(x) を
237        0 から rmax までの区間で計算し、matplotlib を用いてグラフを描画します。
238        各関数の零点も get_bessel_zeros を使用して計算し、グラフ上にマークします。
239    引数:
240        :param lmax: プロットする球Bessel関数の最大次数 (l)。
241        :type lmax: int
242        :param rmax: r軸の最大値。
243        :type rmax: float
244        :param rmesh: r軸のデータポイント数。
245        :type rmesh: int
246    戻り値:
247        :returns: なし
248        :rtype: None
249    """
250
251    plt.figure(figsize=(8, 6))
252    plt.title(f"Spherical Bessel Function $j_l(kr)$", fontsize=16)
253    plt.xticks(fontsize=16)
254    plt.yticks(fontsize=16)
255
256    r = np.linspace(0, rmax, rmesh)
257    for l in range(lmax + 1):
258        j_l = spherical_jn(l, r)
259        zero_points = get_bessel_zeros(l, xmin = 0.1, xmax = 10.0, dx = 0.1, print_level = 1)
260        print("zeros: l=", l, zero_points)
261
262        plt.plot(r, j_l, label = f"$j_{l}(kr)$")
263        plt.plot(zero_points, [0.0] * len(zero_points), label = f"zero points of $j_{l}(kr)$", linestyle = '', marker = 'o')
264
265    plt.axhline(y=0, color='red', linestyle='--', linewidth = 0.5)
266    plt.axvline(x=0, color='red', linestyle='--', linewidth = 0.5)
267    plt.xlabel(r"$k$$\cdot$$r$", fontsize=16)
268    plt.ylabel(r"$j_l$($k$$\cdot$$r$)", fontsize=16)
269#    plt.grid(True)
270
271    plt.legend(fontsize=12)
272    
273    plt.show()
274
275def plot_H(nmax, rmax):
276    """
277    概要:
278        水素原子の動径波動関数 R_nl(r) をプロットします。
279    詳細説明:
280        水素原子の動径波動関数 R_nl(r) を、主量子数 n と軌道角運動量量子数 l の
281        組み合わせに対して計算し、matplotlib を用いてグラフを描画します。
282        genlaguerre 関数（一般化されたラゲール多項式）を用いて計算が行われます。
283        プロットされる動径 r の範囲は 0 から 20 まで（ボーア半径 a0 単位）。
284        なお、rmax 引数は現在の実装ではプロット範囲に影響を与えません。
285    引数:
286        :param nmax: プロットする主量子数 n の最大値。n は 1 から nmax まで。
287        :type nmax: int
288        :param rmax: (未使用) r軸の最大値として想定される値（ボーア半径 a0 単位）。
289        :type rmax: float
290    戻り値:
291        :returns: なし
292        :rtype: None
293    """
294    a0 = 1.0
295
296    plt.figure(figsize=(8, 6))
297    plt.title("Radial Wave Function $R_{20}(r)$ for Hydrogen Atom", fontsize=16)
298    plt.xticks(fontsize=16)
299    plt.yticks(fontsize=16)
300
301    r = np.linspace(0, 20, 500)
302    for n in range(nmax + 1):
303        for l in range(n):
304            rho = 2 * Z * r / (n * a0)
305            Lnl = genlaguerre(n - l - 1, 2 * l + 1)(rho)
306            R_nl = (rho ** l) * np.exp(-rho / 2) * Lnl
307#            R_nl = sqrt((2 * Z / (n * a0))**3 \
308#                * factorial(n - l - 1) / (2 * n * factorial(n + l))) \
309#                * R_nl
310            plt.plot(r, R_nl, label = f"$R_{n}$$_{l}$($r$)")
311
312    plt.axhline(y=0, color='red', linestyle='--', linewidth = 0.5)
313    plt.axvline(x=0, color='red', linestyle='--', linewidth = 0.5)
314    plt.xlabel("r (in units of $a_0$)", fontsize=16)
315    plt.ylabel("$R_{nl}$($r$)", fontsize=16)
316#    plt.grid(True)
317
318    plt.legend(fontsize=12)
319    
320    plt.show()
321
322def main():
323    """
324    概要:
325        スクリプトのメイン実行関数。
326    詳細説明:
327        コマンドライン引数 sys.argv[1] の値に基づいて、実行モードを決定します。
328        'cal' モードでは cal() 関数を呼び出し、量子球のエネルギー準位計算と表示を実行します。
329        'plot' モードでは plot_spherical_bessel() 関数を呼び出し、球Bessel関数のプロットを実行します。
330        'plotH' モードでは plot_H() 関数を呼び出し、水素原子の動径波動関数のプロットを実行します。
331        上記以外のモードが指定された場合はエラーメッセージを表示し、スクリプトを終了します。
332        
333    戻り値:
334        :returns: なし
335        :rtype: None
336    """
337    if mode == 'cal':
338        cal()
339    elif mode == 'plot':
340        plot_spherical_bessel(lmax = 3, rmax = 10.0)
341    elif mode == 'plotH':
342        plot_H(nmax = 3, rmax = 10.0)
343    else:
344        print(f"Error: Invalid mode={mode}")
345        exit()
346
347if __name__ == '__main__':
348    main()