pw1d.py 技術ドキュメント

プログラムの動作

pw1d.py は、1次元周期ポテンシャル中の電子の量子力学的振る舞いを、平面波基底セットを用いた第一原理計算によってシミュレートするPythonスクリプトです。このプログラムは、以下の主要な機能を提供し、1次元結晶における電子状態の解析を可能にします。

ポテンシャルのフーリエ変換 (ftモード): 定義された実空間ポテンシャル（矩形バリアまたはガウスバリア）をフーリエ変換し、その結果を可視化します。これにより、周期ポテンシャルが逆格子空間でどのように表現されるかを確認できます。
バンド構造の計算とプロット (bandモード): 1次元ブロッホ関数と平面波展開法を用いて、電子のエネルギーバンド構造を計算し、波数（$k$）に対するエネルギー（$E$）の関係をプロットします。これにより、結晶中の電子が取りうるエネルギー準位とその分散関係を視覚的に把握できます。
波動関数の計算とプロット (wfモード): 特定の波数 $k$ とエネルギー準位における電子の波動関数を計算し、実空間でのポテンシャルとともにプロットします。波動関数の実部、虚部、および確率密度 $| \Psi |^2$ が表示され、電子の空間的な分布を解析できます。

このプログラムは、固体物理学におけるバンド計算の基本的な原理を理解するための教育ツールとしても利用できます。

原理

本プログラムは、1次元周期ポテンシャル $V(x)$ 中の電子に対するシュレディンガー方程式を、平面波基底セットを用いて解くことに基づいています。

\[ \left( -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + V(x) \right) \psi(x) = E \psi(x) \]

ここで、$\hbar$ は換算プランク定数、$m$ は電子の質量、$\psi(x)$ は波動関数、$E$ は電子のエネルギーです。

Blochの定理: 周期 $a$ を持つポテンシャル $V(x+a) = V(x)$ の場合、波動関数はBlochの定理により次のように表されます。 $$ \psi_k(x) = e^{ikx} u_k(x) $$ ここで、$k$ は波数、$u_k(x)$ は結晶の周期 $a$ を持つ関数 ($u_k(x+a) = u_k(x)$) です。
フーリエ展開: 周期関数 $V(x)$ と $u_k(x)$ はフーリエ級数で展開できます。 $$ V(x) = \sum_G V_G e^{iGx} $$ u_k(x) = \sum_G C_G e^{iGx} $$ ここで、$G = n \frac{2\pi}{a}$ は逆格子ベクトル（$n$ は整数）であり、$V_G$ はポテンシャルのフーリエ係数、$C_G$ は波動関数の係数です。ポテンシャルのフーリエ係数 $V_G$ は、プログラムでは np.fft.fft を用いて計算されます。
固有値問題の導出: これらの展開をシュレディンガー方程式に代入し、整理すると、係数 $C_G$ に関する連立方程式が得られます。これは行列形式で表され、固有値問題として解かれます。 $$ \sum_{G'} H_{G,G'} C_{G'} = E C_G $$ 行列要素 $H_{G,G'}$ は次のように与えられます。 $$ H_{G,G'} = \left( \frac{\hbar^2}{2m} (k+G)^2 \right) \delta_{G,G'} + V_{G-G'} $$ ここで、$\delta_{G,G'}$ はクロネッカーのデルタです。
- 対角要素 ($G = G'$): $H_{G,G} = \frac{\hbar^2}{2m} (k+G)^2 + V_0$ これは、波数 $k+G$ を持つ自由電子の運動エネルギー項と、ポテンシャルの平均値（$G=0$ のフーリエ係数）から構成されます。プログラム内では free_KE(k, iG) が運動エネルギー項を計算し、find_Vft(0, Glist, Vftlist) が $V_0$ を返します。
- 非対角要素 ($G \neq G'$): $H_{G,G'} = V_{G-G'}$ これは、$G$ と $G'$ の間で電子を散乱させるポテンシャルのフーリエ係数 $V_{G-G'}$ に対応します。プログラム内では find_Vft(dij, Glist, Vftlist) が $V_{G-G'}$ を返します。

この行列を対角化することで、固有エネルギー $E$（バンド構造）と、それに対応する波動関数の係数 $C_G$ が得られます。プログラムでは numpy.linalg.eig 関数がこの固有値問題を解くために使用されます。

ポテンシャルの種類: プログラムは以下の2種類の1次元周期ポテンシャルをサポートしています。
- 矩形バリア (pottype = 'rect'): $$ V(x) = \begin{cases} bpot & (0 \le x \le bwidth) \\ 0 & (\text{otherwise}) \end{cases} $$ このポテンシャルは周期 $a$ で繰り返されます。
- ガウスバリア (pottype = 'gauss'): $$ V(x) = bpot \exp\left(-\left(\frac{x - 0.5a}{bwidth}\right)^2\right) $$ これも周期 $a$ で繰り返されます。

必要な非標準ライブラリとインストール方法

このプログラムは以下の非標準Pythonライブラリを使用します。

numpy: 数値計算、特に配列操作、線形代数、フーリエ変換に必要です。
matplotlib: 計算結果をグラフとしてプロットし、可視化するために必要です。

これらのライブラリは pip コマンドを使用してインストールできます。

pip install numpy matplotlib

必要な入力ファイル

pw1d.py は外部の入力ファイルを必要としません。すべての設定パラメータはプログラムコード内にデフォルト値として定義されているか、またはプログラム実行時にコマンドライン引数として直接指定されます。

生成される出力ファイル

pw1d.py はディスク上にファイルとして出力ファイルを生成しません。プログラムの実行結果は、以下の形式で出力されます。

標準出力: 計算の進行状況、入力パラメータ、フーリエ変換されたポテンシャルの係数、計算されたエネルギー準位、波動関数の係数などのテキスト情報が表示されます。
グラフィックウィンドウ: matplotlib ライブラリを使用して、以下のようなプロットが画面にリアルタイムで表示されます。
- ft モード: 実空間ポテンシャル、生のFFT結果、変換後のG空間ポテンシャル。
- band モード: エネルギーバンド構造（E-k図）。
- wf モード: 実空間ポテンシャルと、特定の波動関数の実部、虚部、および確率密度。

これらのプロットは、ユーザーがターミナルで Enter キーを押すまで画面上に保持されます。

コマンドラインでの使用例 (Usage)

プログラムは実行モード（ft, band, wf）とそれに付随するパラメータをコマンドライン引数として受け取ります。

基本的な使用法は以下の通りです。

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)

引数の詳細:

mode: 実行モード (ft, band, wf のいずれか)。
a: 格子定数（Å）。
na: FFTのための空間分割数（2のべき乗である必要があります）。
pottype: ポテンシャルの種類 (rect または gauss)。
bwidth: ポテンシャルの幅または広がり（Å）。
bpot: ポテンシャルの高さ（eV）。
nG: 平面波基底関数の数。
kmin: バンド計算の波数 $k$ の最小値（$\pi/a$ 単位）。
kmax: バンド計算の波数 $k$ の最大値（$\pi/a$ 単位）。
nk: バンド計算の $k$ 点数。
kw: 波動関数計算の波数 $k$（$\pi/a$ 単位）。
iLevel: 波動関数をプロットするエネルギー準位のインデックス（0から始まる）。
xwmin: 波動関数プロットの $x$ 軸最小値（Å）。
xwmax: 波動関数プロットの $x$ 軸最大値（Å）。
nxw: 波動関数プロットの $x$ 点数。

括弧 () で囲まれた引数はオプションであり、省略された場合はプログラム内部のデフォルト値が使用されます。

コマンドラインでの具体的な使用例

以下に、各モードでの具体的な実行例と、その結果として表示される内容について説明します。

ftモード (ポテンシャルのフーリエ変換)

実空間のポテンシャルとそのフーリエ変換結果をプロットします。この例では、格子定数 5.4064 Å、FFT分割数 64、幅 0.5 Å、高さ 10.0 eV の矩形バリアポテンシャルを使用しています。
```
python pw1d.py ft 5.4064 64 rect 0.5 10.0
```
実行結果の説明: グラフィックウィンドウに3つのプロットが表示されます。
- 左上: 実空間における周期ポテンシャル $V(x)$。
- 左下: ポテンシャルのフーリエ変換の生のデータ（np.fft.fftの結果）の実部と虚部。横軸はFFTのインデックス。
- 右下: G空間におけるポテンシャルのフーリエ係数 $V_G$ の実部、虚部、絶対値。横軸は逆格子ベクトル $G$（単位：$1/\text{Å}$）。ポテンシャルの周期性により、特定の $G$ 点に大きな値を持つことが分かります。
bandモード (バンド構造の計算)

電子のエネルギーバンド構造を計算し、プロットします。この例では、上記と同じポテンシャル設定で、7つの平面波基底（nG=7）を使用し、$k$ 軸を $-\frac{\pi}{a}$ から $\frac{\pi}{a}$ まで21点で計算します。
```
python pw1d.py band 5.4064 64 rect 0.5 10.0 7 -0.5 0.5 21
```
実行結果の説明: グラフィックウィンドウにE-k図が表示されます。
- 横軸は波数 $k$（単位：$\pi/a$）、縦軸はエネルギー $E$（単位：eV）。
- 赤い線は、ポテンシャルの影響がない場合の自由電子の分散関係 $E = \frac{\hbar^2}{2m}(k+G)^2$ を示します。
- 黒い丸は、周期ポテンシャルを考慮した計算で得られたバンド構造を表します。自由電子のバンドがブリルアンゾーン境界でギャップを開き、バンド構造が形成される様子が観察できます。
wfモード (波動関数の計算)

特定の波数 $k$ とエネルギー準位に対応する波動関数を計算し、プロットします。この例では、上記と同じポテンシャル設定、7つの平面波基底、波数 $k=0.0$、基底状態（iLevel=0）の波動関数を、実空間 $x$ 軸 0.0 Å から 16.2192 Å (約 $3a$) の範囲で101点計算します。
```
python pw1d.py wf 5.4064 64 rect 0.5 10.0 7 0.0 0 0.0 16.2192 101
```
実行結果の説明: グラフィックウィンドウに実空間ポテンシャルと波動関数が重ねてプロットされます。
- ポテンシャル $V(x)$ が薄い線で表示されます。
- 波動関数の実部 ($\Psi_{real}$)、虚部 ($\Psi_{imaginary}$)、および確率密度 ($|\Psi|^2$) が表示されます。
- $k=0.0$ の場合、波動関数は通常、実数または純虚数になり、確率密度は周期ポテンシャルの影響を受けて局在化または非局在化のパターンを示します。特に、ポテンシャルバリアの内部や外部での波動関数の振る舞いを観察できます。