wavefunction2D.py 技術ドキュメント

プログラムの動作

wavefunction2D.py は、様々な物理モデルにおける波動関数を2次元平面で計算し、その実部、虚部、および確率密度を可視化するPythonスクリプトです。主に、以下の3つのモードをサポートしています。

  1. pw (Plain Wave - 自由電子): Blochの定理に基づく平面波モデルの波動関数を計算します。

  2. qbox (Quantum Box - 量子井戸/量子ドット): 無限深の量子井戸（量子ドット）モデルにおける定在波的な波動関数にBloch波の位相項を組み合わせた波動関数を計算します。

  3. H (Hydrogen-like - 水素様原子): 水素様原子の動径波動関数と球面調和関数を組み合わせた波動関数にBloch波の位相項を組み合わせた波動関数を計算します。

プログラムは、X軸に沿った1次元プロット、Y軸に沿った1次元プロット、およびXY平面における2次元（3D表面・コンター）プロットを生成します。これらのプロットは matplotlib ライブラリを使用してインタラクティブに表示されます。また、計算された波動関数の一部の値を標準出力に出力します。

原理

wavefunction2D.py は、Blochの定理に基づいて、周期的なポテンシャル場における波動関数をモデル化しています。波動関数 \(\Psi(\mathbf{r})\) は一般に、周期関数 \(u_{\mathbf{k}}(\mathbf{r})\) と平面波成分 \(e^{i \mathbf{k} \cdot \mathbf{r}}\) の積として表されます。

\[\Psi(\mathbf{r}) = u_{\mathbf{k}}(\mathbf{r}) e^{i \mathbf{k} \cdot \mathbf{r}}\]

ここで \(\mathbf{k} = (k_x, k_y, k_z)\) はBloch波ベクトルです。本プログラムでは、\(\mathbf{r}=(x,y,z)\) の座標が単位格子内 \([-0.5, 0.5)\) にマッピングされた \(x_0, y_0, z_0\) を用いて \(u_{\mathbf{k}}(\mathbf{r})\) を構築し、Bloch波項は \(\phi(x,y,z) = e^{i 2\pi (k_x x/a_x + k_y y/a_y + k_z z/a_z)}\) として計算されます。ここで \(a_x, a_y, a_z\) は格子定数です。

プログラムの中心的な機能は、以下の波動関数計算モードにあります。

1. pw (Plain Wave) モード

自由電子の波動関数を計算します。このモードでは、\(u_{\mathbf{k}}(\mathbf{r})\) 部分も平面波として扱われます。したがって、量子数 \(n_1, n_2, n_3\) は Bloch 波ベクトルの整数倍として機能します。

\[\Psi_{\text{pw}}(x,y,z) = e^{i 2\pi \left( (n_1+k_x) \frac{x}{a_x} + (n_2+k_y) \frac{y}{a_y} + (n_3+k_z) \frac{z}{a_z} \right)}\]

ここで \(n_1, n_2, n_3\) はコマンドライン引数で与えられる「量子数」であり、Bloch波ベクトルの定在波成分を表します。\(k_x, k_y, k_z\) は Bloch 波ベクトルの繰り込み帯域内の成分です。

2. qbox (Quantum Box) モード

無限深の量子井戸（量子ドット）モデルの波動関数にBloch波の位相項を組み合わせます。単位格子内の座標 \(x_0, y_0, z_0\) は reduce2unitcell 関数によって \([-0.5, 0.5)\) の範囲に変換され、さらに kqbox 係数で空間的な広がりが調整されます。 1次元の波動関数成分は、量子数 \(n_i\) に応じて \(\sin(n_i \pi x_0 / a_x)\) または \(\cos(n_i \pi x_0 / a_x)\) の形を取ります。

\[F(x_0,y_0,z_0; n_1,n_2,n_3) = f_x(x_0; n_1) f_y(y_0; n_2) f_z(z_0; n_3)\]

ここで \(f_x(x_0; n_1) = \begin{cases} \cos(n_1 \pi x_0 / a_x \cdot \text{kqbox}[n_1-1]) & \text{if } n_1 \text{ is odd} \\ \sin(n_1 \pi x_0 / a_x \cdot \text{kqbox}[n_1-1]) & \text{if } n_1 \text{ is even} \end{cases}\) \(f_y\) も同様に定義されますが、\(f_z\) はコード内で常に \(1.0\) と定義されているため、Z方向には定在波的な依存性がありません。

最終的な波動関数は、中心の単位セルと隣接する単位セルからの複数のBloch波項の重ね合わせとして計算されます。

\[\Psi_{\text{qbox}}(x,y,z) = \sum_{v_x, v_y, v_z \in \{-1,0,1\}} F(x_0,y_0,z_0; n_1,n_2,n_3) \cdot e^{i 2\pi \left( (k_x+v_x) \frac{x}{a_x} + (k_y+v_y) \frac{y}{a_y} + (k_z+v_z) \frac{z}{a_z} \right)}\]

この和は、周期的な境界条件下で局在した波動関数を表現するための近似的なアプローチと考えられます。

3. H (Hydrogen-like) モード

水素様原子の波動関数にBloch波の位相項を組み合わせます。波動関数は動径波動関数 \(R_{nl}(r)\) と球面調和関数 \(Y_{lm}(\theta, \phi)\) の積で表されます。

\[R_{nl}(r) Y_{lm}(\theta, \phi)\]

動径波動関数 \(R_{nl}(r)\) は主量子数 \(n\) と方位量子数 \(l\) に依存し、球面調和関数 \(Y_{lm}(\theta, \phi)\) は方位量子数 \(l\) と磁気量子数 \(m\) に依存します。プログラムでは、特定の \(n, l, m\) の組み合わせに対応する具体的な数式が実装されています。動径部分のスケールは kR 係数と原子番号 \(Z\)、ボーア半径 \(a_B\) によって調整されます。球面調和関数は、実数値で扱えるように虚数部分が分離されています（例：\(Y_{1,1} \propto \sin\theta \cos\phi\), \(Y_{1,-1} \propto \sin\theta \sin\phi\)）。 座標は reduce2unitcell で \([-0.5, 0.5)\) に変換された \(x_0, y_0, z_0\) を使用して極座標 \(r_0, \theta_0, \phi_0\) が計算されます。

\[\Psi_{\text{H}}(x,y,z) = \sum_{v_x, v_y, v_z \in \{-1,0,1\}} R_{nl}(r_0) Y_{lm}(\theta_0, \phi_0) \cdot e^{i 2\pi \left( (k_x+v_x) \frac{x}{a_x} + (k_y+v_y) \frac{y}{a_y} + (k_z+v_z) \frac{z}{a_z} \right)}\]

これも qbox モードと同様に複数のBloch波項の重ね合わせで表現されます。

共通の処理

  • reduce2unitcell(x): 座標 \(x\) を中心を \(0.0\) とする単位格子内 \([-0.5, 0.5)\) の範囲にマッピングします。これにより、格子定数 \(a_x, a_y, a_z\) で定義される周期的な系において、単位セル内の波動関数を考えることができます。

  • phi(x, y, z, ...): Bloch波の位相項 \(e^{i 2\pi (k_x x/a_x + k_y y/a_y + k_z z/a_z)}\) を計算します。これはすべてのモードで共通に使用されます。

  • 計算は \(z=0.0\) のXY平面上で行われ、3次元メッシュ数 nmeshz は使用されません。プロット範囲 nlxrange と nlyrange はデフォルトで \([-1.5, 1.5]\) に設定されており、これは中心の単位セルと周囲の単位セルをカバーすることで、周期的な挙動を視覚化できるようにしています。

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

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

  • numpy: 数値計算を効率的に行うためのライブラリ。

  • matplotlib: グラフ描画のためのライブラリ。mpl_toolkits.mplot3d は matplotlib の一部です。

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

pip install numpy matplotlib

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。 すべての設定（モード、量子数、Bloch波ベクトル、メッシュ数など）は、スクリプト内の初期値として定義されるか、コマンドライン引数として与えられます。

生成される出力ファイル

このプログラムは、いかなるファイルも生成しません。 計算結果は標準出力にテキスト形式で出力されるほか、matplotlib を用いてグラフウィンドウにインタラクティブに表示されます。

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

プログラムの基本的な実行コマンドと引数は以下の通りです。

python wavefunction2D.py MODE N1 N2 N3 [KX KY KZ] [NMESHX NMESHY NMESHZ]

引数の説明：

  • MODE: 必須。計算モードを指定します。以下のいずれかを指定します。

    • pw: 自由電子（平面波）モデル。

    • qbox: 量子ドットモデル（無限深量子井戸）。

    • H: 水素様原子モデル。

  • N1, N2, N3: 必須。選択されたモードに応じた量子数を指定します。

    • pw モードの場合: kx0, ky0, kz0 (Bloch波ベクトルの整数部分 \(n_x, n_y, n_z\) に対応)

    • qbox モードの場合: nx, ny, nz (量子数)

    • H モードの場合: n, l, m (量子数)

  • KX, KY, KZ: オプション。Bloch波ベクトルの繰り込み帯域内の成分 (\(k_x, k_y, k_z\)) を指定します。単位は \(\pi/a\) です。デフォルト値は 0.0, 0.0, 0.0 です。

  • NMESHX, NMESHY, NMESHZ: オプション。波動関数の計算に使用するX, Y, Z方向のメッシュ数を整数で指定します。デフォルト値は 120, 120, 120 です。NMESHZ は2Dプロットのため実際には使用されません。

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

例1: 自由電子（pwモード）

X方向に波数 \(0.5 \pi/a\) を持つ平面波（定在波成分なし）の波動関数を計算します。

python wavefunction2D.py pw 0 0 0 0.5 0.0 0.0 120 120 120

実行結果の説明:

  • 標準出力には、X軸上、Y軸上 (\(y=0.0, z=0.0\))、およびXY平面 (\(z=0\)) 上で計算された波動関数の実部 (\(\Psi_r\))、虚部 (\(\Psi_i\))、および確率密度 (\(|\Psi|^2 = \Psi_r^2 + \Psi_i^2\)) の値が表示されます。Bloch波項 \(\phi(x)\) の実部と虚部も表示されます。

  • matplotlib ウィンドウが表示され、以下のグラフが表示されます。

    • 上段左: XY平面における波動関数実部のコンタープロット。

    • 上段中央: XY平面における波動関数虚部のコンタープロット。

    • 上段右: XY平面における波動関数実部の3D表面プロット。

    • 下段左: X軸に沿った波動関数実部と虚部の1Dプロット。

    • 下段中央: Y軸に沿った波動関数実部と虚部の1Dプロット。

    • 下段右: XY平面における波動関数虚部の3D表面プロット。

  • この例では、\(\Psi(x,y,z) = e^{i 2\pi (0.5 x/a_x)} = \cos(\pi x/a_x) + i \sin(\pi x/a_x)\) となるため、X方向に周期的な振動を示す波が描画されます。実部はX軸に関して偶関数、虚部は奇関数となります。

例2: 量子ドット（qboxモード）

X, Y方向に量子数 \(n_x=1, n_y=1\) を持ち、Bloch波項を持たない (\(k_x=k_y=k_z=0\)) 量子ドットの波動関数を計算します。

python wavefunction2D.py qbox 1 1 0 0.0 0.0 0.0 120 120 120

実行結果の説明:

  • 標準出力には、計算された波動関数の数値データが表示されます。

  • matplotlib ウィンドウが表示され、qbox モデルにおける波動関数の実部と虚部の2D/3Dプロット、および1Dスライスプロットが描画されます。

  • この例では、量子数 \((1,1,0)\) の定在波的な波動関数が描画されます。Bloch波項がないため、波動関数は単位セルに対して実数となり、原点付近で最大値を取るコサイン形状の波形が見られるでしょう。XY平面の等高線図では、格子中心 (\(x=0, y=0\)) にピークを持つ定在波パターンが確認できます。

例3: 水素様原子（Hモード）

\(1s\) 軌道 (\(n=1, l=0, m=0\)) の水素様原子波動関数にBloch波項を考慮しない (\(k_x=k_y=k_z=0\)) 場合を計算します。

python wavefunction2D.py H 1 0 0 0.0 0.0 0.0 120 120 120

実行結果の説明:

  • 標準出力には、計算された波動関数の数値データが表示されます。

  • matplotlib ウィンドウが表示され、H モデルにおける波動関数の実部と虚部の2D/3Dプロット、および1Dスライスプロットが描画されます。

  • この例では、1s軌道の波動関数は動径方向のみに依存する球対称な形状となり、実部のみが存在します（球面調和関数 \(Y_{00}\) は実数）。Bloch波項がないため、原点を中心に最大値を持つガウス関数的な形状が見られるでしょう。グラフタイトルには mode:H n=(1,0,0)1s k=(0.0,0.0,0.0) のように表示されます。