pw1d.py 技術ドキュメント
このドキュメントは、Pythonプログラム pw1d.py の技術的な側面、機能、および使用方法について説明します。
プログラムの動作
pw1d.py は、平面波基底セット法を用いて1次元周期ポテンシャル中の電子の挙動を計算するPythonスクリプトです。主に以下の3つのモードで動作します。
ft(Fourier Transform) モード: 設定された周期ポテンシャルの実空間形状と、そのフーリエ変換の結果を可視化します。これにより、周期ポテンシャルが逆格子空間でどのように表現されるかを理解できます。band(Band Structure) モード: 1次元結晶中の電子のエネルギーバンド構造（エネルギー分散関係 \(E(k)\)）を計算し、グラフとしてプロットします。自由電子のエネルギーバンドも同時に表示され、ポテンシャルの影響を比較できます。wf(Wave Function) モード: 特定の結晶波数 \(k\) とエネルギー準位iLevelに対応する電子の波動関数を計算し、実空間でプロットします。ポテンシャルと波動関数の実部、虚部、および確率密度 \(|\Psi(x)|^2\) が表示されます。
このプログラムは、固体物理学におけるバンド計算の基本的な概念、特に平面波基底セット法の適用と周期ポテンシャルの影響を理解するための教育ツールとして機能します。
原理
pw1d.py は、周期ポテンシャル中の電子に対する1次元のシュレーディンガー方程式を、平面波基底セット法（Plane Wave Basis Set Method）を用いて数値的に解くことで、エネルギーバンド構造や波動関数を求めます。
周期ポテンシャル プログラムは、格子定数 \(a\) を持つ1次元の周期ポテンシャル \(V(x)\) を用います。ポテンシャルは単位胞 \(0 \le x < a\) の範囲で定義され、周期的に繰り返されます。
矩形バリア (
rect): $\(V(x) = \begin{cases} b_{\text{pot}} & \text{if } 0 \le x < b_{\text{width}} \\ 0 & \text{otherwise} \end{cases}\)\( ここで \)x\( は \)0 \le x < a$ に正規化された座標です。ガウス型 (
gauss): $\(V(x) = b_{\text{pot}} \exp\left(-\left(\frac{x - 0.5a}{b_{\text{width}}}\right)^2\right)\)$ これらのポテンシャルはpot(x)関数によって実装されています。
ポテンシャルのフーリエ変換 周期ポテンシャル \(V(x)\) はフーリエ級数展開されます。 $\(V(x) = \sum_{G} V_G e^{i G x}\)\( ここで \)G = \frac{2\pi}{a} m\( は逆格子ベクトル（\)m\( は整数）です。フーリエ係数 \)V_G\( は以下の積分で与えられます。 \)\(V_G = \frac{1}{a} \int_0^a V(x) e^{-i G x} dx\)\( プログラムでは、`cal_fft` 関数がFFT（高速フーリエ変換）を用いて \)V_G$ を数値的に計算します。
シュレーディンガー方程式と平面波基底 1次元の結晶中の電子の波動関数は、ブロッホの定理により \(\Psi_k(x) = u_k(x) e^{i k x}\) の形をとります。ここで \(u_k(x)\) は実空間の周期 \(a\) を持つ関数、\(k\) は結晶波数です。 \(u_k(x)\) を平面波基底で展開します。 $\(u_k(x) = \sum_{G} c_G e^{i G x}\)\( これをシュレーディンガー方程式に代入し、各平面波基底に対する係数 \)c_G\( を求める固有値問題を解きます。その結果、以下の行列方程式が得られます。 \)\( \sum_{G'} \left( \frac{\hbar^2 (k+G)^2}{2m_e} \delta_{G,G'} + V_{G-G'} \right) c_{G'} = E c_G \)\( これは行列形式で \)\mathbf{H} \mathbf{c} = E \mathbf{c}\( と書くことができ、行列要素 \)H_{G,G'}\( は以下のようになります。 \)\(H_{G,G'} = \frac{\hbar^2 (k+G)^2}{2m_e} \delta_{G,G'} + V_{G-G'}\)\( ここで、\)\delta_{G,G'}$ はクロネッカーのデルタです。
対角要素 (\(G=G'\)): 自由電子の運動エネルギー項とポテンシャルの \(G=0\) 成分 (\(V_0\))。 $\(H_{G,G} = \frac{\hbar^2 (k+G)^2}{2m_e} + V_0\)$
非対角要素 (\(G \ne G'\)): ポテンシャルのフーリエ成分 \(V_{G-G'}\)。 $\(H_{G,G'} = V_{G-G'}\)\( プログラム中の `solve_pw` 関数は、この行列 \)\mathbf{H}\( を構築し、`numpy.linalg.eig` を用いて固有値（エネルギー \)E\(）と固有ベクトル（平面波の係数 \)c_G$）を計算します。
物理定数 プログラムでは、以下の物理定数が定義されています。
pi,hbar(換算プランク定数),me(電子質量),e(電気素量) など。 運動エネルギーの係数KEは、hbar * hbar / (2.0 * me * e) * 1.0e20として計算され、エネルギーが eV 単位、長さが Å (オングストローム) 単位になるように調整されています。 具体的には、運動エネルギー項の係数は以下のように定義されます。 $\(KE_{\text{coeff}} = \frac{\hbar^2}{2m_e e} \times 10^{20} \quad [\text{eV} \cdot \AA^2]\)\( そして、実際に利用される係数 \)kK\( は、結晶波数 \)k\( が \)\pi/a\( 単位で入力されることを考慮して計算されます。 \)\(kK = KE_{\text{coeff}} \left(\frac{2\pi}{a}\right)^2 \quad [\text{eV}]\)\( これにより、自由電子の運動エネルギー \)E_{\text{free}}(k, G)\( は、入力された \)k\( (単位 \)\pi/a\() と逆格子ベクトル \)G\( (単位 \)\pi/a\() を用いて、eV 単位で計算されます。 \)\(E_{\text{free}}(k, G) = kK (k+G)^2\)$
必要な非標準ライブラリとインストール方法
pw1d.py を実行するためには、以下の非標準Pythonライブラリが必要です。
numpy: 数値計算（配列操作、高速フーリエ変換、線形代数演算など）のために使用されます。matplotlib: グラフ描画のために使用されます。
これらのライブラリは、Pythonのパッケージマネージャーである pip を使ってインストールできます。
pip install numpy matplotlib
必要な入力ファイル
pw1d.py は、外部の入力ファイルを必要としません。
全ての計算パラメータはプログラムのソースコード内に初期値として定義されており、コマンドライン引数を通じて実行時に変更することが可能です。
生成される出力ファイル
pw1d.py は、いかなるファイルもディスクに保存しません。
プログラムの実行結果は以下の形式で出力されます。
標準出力 (コンソール): プログラムの動作状況、入力パラメータ、計算途中の詳細な数値データ（例：フーリエ変換されたポテンシャル係数、固有エネルギー、波動関数の係数など）が出力されます。
グラフィカルウィンドウ:
matplotlibを用いて、計算結果がグラフとしてリアルタイムで表示されます。ftモード: ポテンシャルとそのフーリエ変換。bandモード: エネルギーバンド構造。wfモード: 波動関数とポテンシャル。 これらのグラフウィンドウは、ユーザーがコンソールでENTERキーを押すまで開いたままになります。
コマンドラインでの使用例 (Usage)
基本的な使用方法は以下の通りです。丸括弧 () 内の引数は省略可能です。
Usage: Variables in () are optional
python pw1d.py
python pw1d.py (ft a na pottype bwidth bpot)
python pw1d.py (band a na pottype bwidth bpot nG kmin kmax nk)
python pw1d.py (wf a na pottype bwidth bpot nG kw iLevel xwmin xwmax nxw)
pottype: rect|gauss
各引数の意味は以下の通りです。
mode: 実行モード (ft,band,wf)。省略時はftになります。a: 格子定数 (\(\AA\))。na: FFT計算のための実空間の分割数 (2のN乗である必要があります)。pottype: ポテンシャルの種類 (rect(矩形バリア) またはgauss(ガウス型))。bwidth: ポテンシャルの幅 (\(\AA\))。rectの場合はバリア幅、gaussの場合はガウス関数の幅を表すパラメータ。bpot: ポテンシャルの高さ (eV)。nG: 平面波基底の数。kmin: バンド計算の開始波数（単位は \(\pi/a\)）。kmax: バンド計算の終了波数（単位は \(\pi/a\)）。nk: バンド計算の波数点の数。kw: 波動関数を計算する結晶波数（単位は \(\pi/a\)）。iLevel: 波動関数をプロットするエネルギー準位のインデックス (0から始まる)。xwmin: 波動関数プロットの開始X座標 (\(\AA\))。xwmax: 波動関数プロットの終了X座標 (\(\AA\))。nxw: 波動関数プロットのX軸点の数。
コマンドラインでの具体的な使用例
1. ポテンシャルのフーリエ変換表示 (ft モード)
以下のコマンドは、格子定数 \(a=5.4064 \AA\)、FFT分割数 \(na=64\)、矩形バリアポテンシャル（幅 \(0.5 \AA\)、高さ \(10.0 \text{ eV}\)）を設定し、そのポテンシャルとフーリエ変換結果をプロットします。
python pw1d.py ft 5.4064 64 rect 0.5 10.0
実行結果の説明: 実行すると、グラフィカルウィンドウが開き、4つのサブプロットが表示されます。
左上のプロット: 実空間における周期ポテンシャル \(V(x)\) の形状。
左下のプロット: FFTの生の結果（虚軸が配列インデックス）。
右下のプロット: 正規化された逆格子空間におけるポテンシャルのフーリエ変換 \(|V_G|\) および実部、虚部。 コンソールには、ポテンシャル設定とフーリエ変換の数値結果の一部が表示されます。
2. バンド構造計算 (band モード)
以下のコマンドは、上記と同じポテンシャル設定に加え、7個の平面波基底を用いて、結晶波数 \(k\) を \(-\frac{0.5\pi}{a}\) から \(\frac{0.5\pi}{a}\) まで21点計算し、エネルギーバンド構造をプロットします。
python pw1d.py band 5.4064 64 rect 0.5 10.0 7 -0.5 0.5 21
実行結果の説明: 実行すると、グラフィカルウィンドウが開き、エネルギーバンド構造が表示されます。横軸は結晶波数 \(k\)（単位 \(\pi/a\)）、縦軸はエネルギー (eV) です。黒い丸でプロットされた点が計算されたバンド構造であり、赤い線は比較のための自由電子のエネルギーバンドです。コンソールには、計算に使用されたパラメータ、フーリエ変換されたポテンシャル係数、各 \(k\) 点での固有エネルギーと固有ベクトルの一部が出力されます。
3. 波動関数計算 (wf モード)
以下のコマンドは、上記と同じポテンシャル設定と7個の平面波基底を用い、結晶波数 \(k=0.0\) で、基底状態（0番目のエネルギー準位）の波動関数を計算し、プロットします。プロット範囲は \(x=0.0 \AA\) から \(x=16.2192 \AA\) (\(3a\)) で、101点を使用します。
python pw1d.py wf 5.4064 64 rect 0.5 10.0 7 0.0 0 0.0 16.2192 101
実行結果の説明: 実行すると、グラフィカルウィンドウが開き、ポテンシャルと波動関数が表示されます。横軸は実空間座標 \(x\) (\(\AA\)) です。青線は波動関数の実部、橙線は虚部、緑線は確率密度 \(| \Psi(x) |^2\) を示し、黒線はポテンシャル \(V(x)\) です。コンソールには、計算に使用されたパラメータ、フーリエ変換されたポテンシャル係数、各 \(k\) 点での固有エネルギーと固有ベクトル、選択されたエネルギー準位の波動関数係数などが出力されます。