pw1d.py 技術ドキュメント

プログラムの動作

pw1d.py は、一次元周期ポテンシャル中の電子に対する平面波基底セットを用いたバンド構造計算および波動関数可視化を行うPythonプログラムです。このプログラムの主な目的は、固体物理学における電子のエネルギーバンド理論の基本的な概念を、視覚的に理解しやすく提供することです。

主な機能は以下の3つのモードに分かれます。

  1. フーリエ変換 (ft):

    • 指定された一次元周期ポテンシャル（矩形障壁またはガウス型障壁）を定義し、その空間領域での形状とフーリエ変換結果を可視化します。

    • ポテンシャルの周期性とその逆格子空間表現の関係を理解するのに役立ちます。

  2. バンド構造計算 (band):

    • 平面波基底を用いて、k空間における電子のエネルギーバンド構造を計算しプロットします。

    • 自由電子の分散関係と、周期ポテンシャルによって導入されるバンドギャップの形成を比較して示します。

  3. 波動関数計算 (wf):

    • 特定のk点とエネルギー準位に対応する波動関数を計算し、その実部、虚部、および確率密度 (\(|\Psi|^2\)) を空間領域でプロットします。

    • 同時に、元の周期ポテンシャルもプロットし、波動関数とポテンシャルの関係を視覚化します。

このプログラムは、固体中の電子の量子力学的振る舞いを、基本的なモデルを使ってシミュレートし、その結果をグラフで表示することで、教育的ツールとして機能します。

原理

pw1d.py は、一次元周期ポテンシャル \(V(x)\) 中の電子に対するシュレディンガー方程式を、平面波基底展開法（Plane Wave Basis Set method）を用いて解きます。

シュレディンガー方程式は以下の通りです。 $\( \left[ -\frac{\hbar^2}{2m_e} \frac{d^2}{dx^2} + V(x) \right] \psi(x) = E \psi(x) \)\( ここで、\) \hbar \( はディラック定数、\) m_e \( は電子の質量、\) E \( はエネルギー、\) \psi(x) $ は波動関数です。

  1. Blochの定理: 周期ポテンシャル \(V(x) = V(x+a)\)（\(a\) は格子定数）の場合、波動関数はBlochの定理により以下の形で書けます。 $\( \psi_k(x) = e^{ikx} u_k(x) \)\( ここで、\)k\( はBloch波数、\)u_k(x)\( はポテンシャルと同じ周期を持つ関数です。 \)u_k(x)\( はフーリエ級数で展開できるため、波動関数 \) \psi_k(x) \( も以下のように平面波の重ね合わせとして展開できます。 \)\( \psi_k(x) = \sum_G c_{k,G} e^{i(k+G)x} \)\( ここで、\)G = n \frac{2\pi}{a}\( は逆格子ベクトル（\)n$ は整数）です。

  2. ポテンシャルのフーリエ変換: 周期ポテンシャル \(V(x)\) もフーリエ級数で展開できます。 $\( V(x) = \sum_G V_G e^{iGx} \)$ プログラムでは、cal_fft 関数がこのフーリエ変換を numpy.fft.fft を用いて計算します。ポテンシャルは矩形障壁 (rect) またはガウス型障壁 (gauss) のいずれかとして定義されます。

  3. 固有値問題の定式化: 波動関数とポテンシャルのフーリエ展開をシュレディンガー方程式に代入し、各平面波 \(e^{i(k+G')x}\) の係数を比較することで、以下のような固有値方程式が得られます。 $\( \sum_G \left[ \left( \frac{\hbar^2}{2m_e} (k+G')^2 - E \right) \delta_{G',G} + V_{G'-G} \right] c_{k,G} = 0 \)\( これは行列形式 \)H \mathbf{c} = E \mathbf{c}\( で表され、ハミルトニアン行列 \)H\( の要素は以下のようになります。 \)\( H_{G'G} = \frac{\hbar^2}{2m_e} (k+G)^2 \delta_{G'G} + V_{G'-G} \)\( プログラムでは、`free_KE` 関数が自由電子の運動エネルギー項 \) \frac{\hbar^2}{2m_e} (k+G)^2 \( を計算します。`Vftlist` はポテンシャルのフーリエ変換係数 \)V_G\( を保持し、`find_Vft` 関数で必要な \)V_{G'-G}$ を取得します。

  4. ハミルトニアンの構築と対角化: solve_pw 関数は、指定された逆格子ベクトル数 nG に基づいてハミルトニアン行列 \(H\) を構築します。そして、numpy.linalg.eig 関数を用いてこの行列を対角化し、固有値（エネルギー \(E\)）と固有ベクトル（平面波係数 \(c_{k,G}\)）を計算します。固有値はエネルギーバンドを形成し、固有ベクトルは波動関数を構成する各平面波の重みを表します。

  5. 波動関数の計算: cal_wf 関数は、solve_pw で得られた平面波係数 \(c_{k,G}\) とBloch波数 \(k\) を用いて、空間領域での波動関数 \( \psi_k(x) = \sum_G c_{k,G} e^{i(k+G)x} \) を計算します。これにより、実空間での波動関数の形状と確率密度 \(|\Psi|^2\) を可視化できます。

物理定数（h, hbar, c, e, e0, kB, me, R, a0）はプログラム冒頭で定義され、特に運動エネルギーの係数 KE はeVとオングストローム単位に変換されています。

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

pw1d.py の実行には、以下の非標準ライブラリが必要です。

  • numpy: 数値計算（特に配列操作、フーリエ変換、線形代数）

  • matplotlib: グラフのプロット

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

pip install numpy matplotlib

必要な入力ファイル

pw1d.py は、外部の入力ファイルを必要としません。 全てのパラメータ（格子定数、ポテンシャルの種類と高さ、k点の範囲、基底関数の数など）は、ソースコード内で初期設定されるか、プログラム起動時にコマンドライン引数として与えられます。

生成される出力ファイル

pw1d.py は、ディスクにファイルを生成しません。 全ての出力は以下の形式で行われます。

  • 標準出力 (コンソール):

    • プログラムの起動パラメータのサマリー。

    • フーリエ変換されたポテンシャルの詳細（ft モードおよびband, wf モードで基底抽出の際）。

    • 固有値方程式の計算中の詳細（エネルギー固有値、平面波係数など）。

  • matplotlibによるグラフウィンドウ:

    • ft モード: 周期ポテンシャルの実空間での形状と、そのフーリエ変換（生データおよび整理された形式）のグラフが表示されます。

    • band モード: k空間における電子のエネルギーバンド構造のプロットが表示されます。自由電子の分散関係も比較のためにプロットされます。

    • wf モード: 特定のk点、エネルギー準位での波動関数の実部、虚部、および確率密度が、周期ポテンシャルとともに表示されます。

グラフウィンドウは、ユーザーがコンソールでEnterキーを押すまで表示され続けます。

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

pw1d.py は、実行モードと関連するパラメータをコマンドライン引数として受け取ります。

基本構文:

python pw1d.py (mode) (a na pottype bwidth bpot) [mode specific arguments]

各引数の説明:

  • mode: 実行モード (ft, band, wf)。

  • a: 格子定数（単位: オングストローム）。

  • na: FFTのための空間メッシュ分割数（\(2^n\) である必要があります）。

  • pottype: ポテンシャルの種類 (rect または gauss)。

  • bwidth: ポテンシャルの幅（単位: オングストローム）。rect の場合は障壁の幅、gauss の場合はガウス関数の半値幅に相当。

  • bpot: ポテンシャルの高さ（単位: eV）。

モード固有の引数:

  • band モード:

    • nG: 平面波基底の数（通常は奇数）。

    • kmin: k空間の開始点（単位: \( \pi/a \)）。

    • kmax: k空間の終了点（単位: \( \pi/a \)）。

    • nk: k空間のプロット点数。

  • wf モード:

    • nG: 平面波基底の数。

    • kw: 波動関数を計算するk点（単位: \( \pi/a \)）。

    • iLevel: 波動関数を計算するエネルギー準位のインデックス（0から始まる）。

    • xwmin: 波動関数プロットのx軸最小値（単位: オングストローム）。

    • xwmax: 波動関数プロットのx軸最大値（単位: オングストローム）。

    • nxw: 波動関数プロットのx軸点数。

プログラムが引数なしで起動された場合、または不正な引数が与えられた場合は、上記のUsageメッセージが表示されます。

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

ここでは、pw1d.py の各モードにおける具体的な使用例と、その実行結果について説明します。 デフォルト値はソースコード内のグローバル変数で設定されており、引数を省略するとそのデフォルト値が使用されます。 （例: a = 5.4064, na = 64, pottype = 'rect', bwidth = 0.5, bpot = 10.0 など）

1. フーリエ変換モード (ft)

ポテンシャルの形状と、そのフーリエ変換結果を表示します。

実行コマンド:

python pw1d.py ft 5.4064 64 rect 0.5 10.0

実行結果の説明: このコマンドは、格子定数 \(a=5.4064 \text{ Å}\)、FFT分割数 \(na=64\)、矩形ポテンシャル rect、幅 \(bwidth=0.5 \text{ Å}\)、高さ \(bpot=10.0 \text{ eV}\) を持つ一次元周期ポテンシャルのフーリエ変換を実行します。 コンソールには、入力パラメータのサマリーと、ポテンシャルのフーリエ変換結果（各逆格子ベクトル iG に対応するフーリエ係数 ci）が出力されます。

同時に、matplotlibによって以下の3つのグラフが表示されるウィンドウが開きます。

  • 左上: 実空間における周期ポテンシャル \(V(x)\) の形状。

  • 左下: numpy.fft.fft の生データとして得られるフーリエ変換結果の実部と虚部。横軸はFFTのインデックス。

  • 右下: 整理された形式でのフーリエ変換結果の実部、虚部、絶対値。横軸は正規化された逆格子空間の座標（\(1/ \text{Å}\)）。このグラフから、ポテンシャルの周期性に対応する逆格子点のフーリエ係数が見て取れます。

ウィンドウには "Press ENTER to exit>>" と表示され、Enterキーを押すことでプログラムが終了します。

2. バンド構造計算モード (band)

電子のエネルギーバンド構造を計算しプロットします。

実行コマンド:

python pw1d.py band 5.4064 64 rect 0.5 10.0 3 -0.5 0.5 21

実行結果の説明: このコマンドは、上記のポテンシャル設定に加えて、平面波基底数 nG=3、k空間の範囲 \(kmin=-0.5 \pi/a\) から \(kmax=0.5 \pi/a\) までを nk=21 点で計算し、バンド構造を求めます。 コンソールには、入力パラメータ、ポテンシャルのフーリエ変換結果、各k点でのハミルトニアンの対角化結果（固有エネルギー ei と固有ベクトル ci）が出力されます。

matplotlibによって1つのグラフが表示されるウィンドウが開きます。

  • グラフ: 横軸にBloch波数 \(k\)（\( \pi/a \) 単位）、縦軸にエネルギー \(E\)（eV単位）をとったバンド構造がプロットされます。黒い丸はポテンシャルを考慮したバンド構造を示し、赤い実線は自由電子（ポテンシャルがない場合）の分散関係 \(E = \frac{\hbar^2}{2m_e} (k+G)^2\) を示します。周期ポテンシャルによるバンドギャップの形成が確認できます。

ウィンドウには "Press ENTER to exit>>" と表示され、Enterキーを押すことでプログラムが終了します。

3. 波動関数計算モード (wf)

特定のk点、エネルギー準位での波動関数を計算しプロットします。

実行コマンド:

python pw1d.py wf 5.4064 64 rect 0.5 10.0 3 0.0 0 0.0 16.2192 101

実行結果の説明: このコマンドは、上記のポテンシャル設定および nG=3 の基底数を用いて、k点 \(kw=0.0 \pi/a\)（Γ点）における iLevel=0 番目（基底状態）のエネルギー準位の波動関数を計算します。波動関数は \(xwmin=0.0 \text{ Å}\) から \(xwmax=16.2192 \text{ Å}\) の範囲を nxw=101 点でプロットします。16.2192 は \(3 \times a\) に相当します。 コンソールには、入力パラメータ、フーリエ変換結果、ハミルトニアンの対角化結果、選択されたエネルギー準位の値、およびその波動関数の平面波係数が出力されます。

matplotlibによって1つのグラフが表示されるウィンドウが開きます。

  • グラフ: x軸に空間座標 \(x\)（Å単位）、左のy軸にポテンシャル \(V(x)\)（eV単位）、右のy軸に波動関数 \( \Psi(x) \) の値をとります。黒い線は波動関数の実部 \( \text{Re}(\Psi(x)) \)、オレンジ色の線は虚部 \( \text{Im}(\Psi(x)) \)、緑色の線は確率密度 \( |\Psi(x)|^2 \) を示します。青い線は元の周期ポテンシャル \(V(x)\) です。これにより、ポテンシャルと波動関数の空間的な関係が視覚的に理解できます。

ウィンドウには "Press ENTER to exit>>" と表示され、Enterキーを押すことでプログラムが終了します。