schrodinger1d.py 技術ドキュメント

プログラムの動作

schrodinger1d.py は、1次元時間独立シュレーディンガー方程式をシューティング法とVerlet積分を用いて数値的に解くPythonプログラムです。主に量子力学における粒子の波動関数 \(\psi(x)\) とその確率密度 \(|\psi(x)|^2\) を計算し、結果をグラフ表示およびExcelファイルに出力します。

主な機能:

• 1次元シュレーディンガー方程式の解法: 調和振動子ポテンシャル \(V(x) = 0.5 x^2\) (原子単位系) における波動関数を計算します。ポテンシャル関数は V(x) 関数で定義されており、必要に応じて変更可能です。

• シューティング法: ユーザーが与えたエネルギー \(E\) に対して、指定された境界条件から波動関数を数値積分し、その振る舞いを評価します。

• Verlet積分: 波動関数の数値積分にはVerletアルゴリズムを採用しています。

• 境界条件の選択: 波動関数の初期条件として、「漸近条件 (asymptotic)」または「ゼロ条件 (zero)」を選択できます。

• 発散検出: 波動関数が閾値を超えて発散した場合、その時点で計算を停止し、発散した区間までの結果を報告します。

• 結果の可視化: 計算された波動関数 \(\psi(x)\) および確率密度 \(|\psi(x)|^2\) を matplotlib を用いてグラフ表示します。発散した場合は、プロットが対数スケールに切り替わります。

• データ出力: 計算結果（\(x\), \(\psi(x)\), \(|\psi(x)|^2\)）をExcelファイル (.xlsx) 形式で出力します。

解決する課題: このプログラムは、量子系の波動関数を数値的にシミュレートし、束縛状態のエネルギー固有値や波動関数の形を探索するための基本的なツールを提供します。特に、シューティング法を用いることで、適切なエネルギーがどのような振る舞いを導くかを視覚的に理解することができます。

原理

schrodinger1d.py は、1次元時間独立シュレーディンガー方程式を数値的に解きます。

1. 1次元時間独立シュレーディンガー方程式

原子単位系 (\(m=\hbar=e=1\)) における1次元時間独立シュレーディンガー方程式は、以下の形で表されます。

\[-\frac{1}{2}\frac{d^2\psi(x)}{dx^2} + V(x)\psi(x) = E\psi(x)\]

ここで、\(\psi(x)\) は波動関数、\(V(x)\) はポテンシャルエネルギー、\(E\) は全エネルギー固有値です。

プログラムでは、ポテンシャルとして調和振動子のポテンシャル \(V(x) = 0.5 x^2\) をデフォルトで使用しています。

2. Verlet積分

シュレーディンガー方程式を数値的に解くために、Verlet積分（またはStörmer-Verlet法）が使用されます。上記のシュレーディンガー方程式を整理すると、

\[\frac{d^2\psi(x)}{dx^2} = 2(V(x) - E)\psi(x)\]

となります。この2階微分方程式を有限差分近似すると、離散化されたメッシュ点 \(x_i\) における波動関数 \(\psi_i = \psi(x_i)\) に対して、以下の関係が得られます。

\[\frac{\psi_{i+1} - 2\psi_i + \psi_{i-1}}{\Delta x^2} \approx 2(V(x_i) - E)\psi_i\]

ここで \(\Delta x\) はメッシュ間隔 \(x_{i+1} - x_i\) です。この式を \(\psi_{i+1}\) について解くと、

\[\psi_{i+1} = 2\psi_i - \psi_{i-1} + \Delta x^2 \cdot 2(V(x_i) - E)\psi_i\]

この漸化式を用いて、初期の2点 \(\psi_0\) と \(\psi_1\) から順次 \(\psi_2, \psi_3, \dots\) を計算していきます。プログラム内では \(k = 2(V(x_i) - E)\) と定義されており、

\[\psi_{i+1} = 2\psi_i - \psi_{i-1} + \Delta x^2 \cdot k \cdot \psi_i\]

と実装されています。

3. シューティング法

シューティング法は、境界値問題を初期値問題に変換して解く数値手法です。このプログラムでは、領域の左端 \(x=-L\) における初期条件（\(\psi(-L)\) と \(d\psi/dx(-L)\)）を与え、そこから波動関数を積分していきます。

正確なエネルギー固有値 \(E\) が与えられた場合、波動関数は物理的に妥当な振る舞いをします（例えば、束縛状態では遠方で指数関数的に減衰する）。しかし、\(E\) が固有値からずれると、波動関数は遠方で発散する傾向があります。

プログラムは与えられた \(E\) に対して波動関数を積分し、その結果（発散したかどうか）をユーザーに報告します。ユーザーは、波動関数が発散しないような \(E\) を試行錯誤で探すことになります。

4. 境界条件

プログラムは以下の2種類の境界条件をサポートしています。

• bc="zero" (ゼロ条件): \(x=-L\) において、波動関数がゼロに近く、その微分が微小な値を持つと仮定します。 $\(\psi(-L) = 0.0\)\( \)\(\frac{d\psi}{dx}(-L) = \epsilon\)\( 数値積分では、\)\psi_0 = 0.0\( および \)\psi_1 = \epsilon \cdot \Delta x\( と初期設定されます。\)\epsilon$ は --eps 引数で指定します。

• bc="asymptotic" (漸近条件): \(x=-L\) において、波動関数がその領域での漸近的な振る舞いに従うと仮定します。これは、古典的な禁制領域 (\(V(x) > E\)) において波動関数が指数関数的に減衰する振る舞いを模倣します。 $\(\frac{d\psi}{dx}(-L) = \kappa \psi(-L)\)\( ここで \)\kappa = \sqrt{2(V(-L) - E)}\( です。 数値積分では、初期値 \)\psi_0 = \text{psi0}\(（例: \)10^{-6}\(）を与え、\)\psi_1 = \psi_0(1 + \kappa \cdot \Delta x)$ と初期設定されます。

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

schrodinger1d.py を実行するためには、以下のPythonライブラリが必要です。これらは標準ライブラリには含まれていないため、事前にインストールする必要があります。

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

• pandas: データ構造（DataFrameなど）を提供し、Excelファイル出力に使用されます。

• matplotlib: グラフ描画（プロット）に使用されます。

これらのライブラリは、Pythonのパッケージマネージャーである pip を使用してインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy pandas matplotlib

必要な入力ファイル

このプログラムは、実行時に外部の入力ファイルを必要としません。すべての設定はコマンドライン引数として与えられます。

生成される出力ファイル

プログラムは、以下の種類のファイルを生成または表示します。

1. Excelファイル (.xlsx):

   • ファイル名: デフォルトは schrodinger.xlsx。--outfile 引数で変更可能です。

   • 内容: 計算された波動関数に関するデータが保存されます。

     • x: 空間座標の値 (a.u.)

     • psi: その座標における波動関数 \(\psi(x)\) の値

     • psi2: その座標における確率密度 \(|\psi(x)|^2\) の値

2. プロット画像 (GUI表示):

   • ファイル形式: matplotlib によって生成され、実行時にGUIウィンドウで表示されます。ファイルとしては保存されません。

   • 内容:

     • 波動関数 \(\psi(x)\)（発散していない場合）

     • 確率密度 \(|\psi(x)|^2\)

     • 波動関数が発散した場合、発散開始点を示す垂直線がプロットに追加され、y軸は対数スケールに切り替わります。

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

schrodinger1d.py は、以下の形式でコマンドラインから実行します。

python schrodinger1d.py --E <Energy_Guess> [options]

必須引数:

• --E <Energy_Guess>: エネルギー固有値の初期推測値 (Ha)。この値が波動関数の計算に用いられます。

オプション引数:

• --L <Half_Domain_Size>: 計算領域の半分のサイズ。空間は [-L, L] の範囲で計算されます。デフォルト値は 10.0 です。

• --nx <Number_Of_Mesh_Points>: 計算領域内のメッシュ点数。デフォルト値は 2000 です。

• --psi_max <Divergence_Threshold>: 波動関数がこの値を超えた場合、発散したと判断して計算を停止します。デフォルト値は 1e10 です。

• --report_step <Interval>: 途中経過をコンソールに表示する間隔。0 を指定すると表示されません。デフォルト値は 200 です。

• --outfile <Output_Filename>: 出力されるExcelファイルの名前。デフォルト値は schrodinger.xlsx です。

• --bc <Boundary_Condition>: 左端 \(x=-L\) での境界条件。asymptotic (漸近条件) または zero (ゼロ条件) のいずれかを選択します。デフォルト値は asymptotic です。

• --eps <Initial_Derivative>: 境界条件が zero の場合にのみ使用される、初期微分値 \(d\psi/dx\) (\(x=-L\) における)。デフォルト値は 1e-6 です。

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

ここでは、調和振動子ポテンシャル (\(V(x) = 0.5 x^2\)) における波動関数の計算例を示します。調和振動子のエネルギー固有値は \(E_n = (n + 0.5)\) Ha (ここで \(n = 0, 1, 2, \dots\)) です。基底状態は \(E_0 = 0.5\) Ha、第一励起状態は \(E_1 = 1.5\) Ha です。

例1: 調和振動子の基底状態 (\(E=0.5\) Ha)

基底状態のエネルギー \(E=0.5\) Ha を指定し、漸近条件で波動関数を計算します。波動関数は発散せずに積分が完了し、左右対称な基底状態の波動関数が得られます。

python schrodinger1d.py --E 0.5 --bc asymptotic --L 10 --nx 2000 --outfile schrodinger_E0.5.xlsx

期待される実行結果 (コンソール出力の一部):

=== 1D Schrodinger Solver (atomic units) ===
E = 0.5 Ha,  L = 10.0,  nx = 2000
Boundary condition: asymptotic
[INFO] step=200, x=-8.004, psi=2.253e-02
[INFO] step=400, x=-6.008, psi=1.503e-01
...
[INFO] step=1800, x=7.996, psi=2.253e-02
[RESULT] integration completed successfully
psi(end) = 2.253e-02
[INFO] Data saved to schrodinger_E0.5.xlsx

この後、matplotlib によってプロットウィンドウが表示され、基底状態の波動関数（左右対称、中央にピーク）とその確率密度関数が確認できます。schrodinger_E0.5.xlsx ファイルには計算データが保存されます。

例2: 固有値から少しずれたエネルギー (\(E=0.6\) Ha)

固有値 \(E=0.5\) からわずかにずれたエネルギー \(E=0.6\) Ha を指定して計算します。このようなエネルギーでは、波動関数は通常、計算領域の端で発散します。

python schrodinger1d.py --E 0.6 --bc asymptotic --L 10 --nx 2000 --outfile schrodinger_E0.6_diverged.xlsx

期待される実行結果 (コンソール出力の一部):

=== 1D Schrodinger Solver (atomic units) ===
E = 0.6 Ha,  L = 10.0,  nx = 2000
Boundary condition: asymptotic
[INFO] step=200, x=-8.004, psi=2.000e-02
[INFO] step=400, x=-6.008, psi=1.385e-01
...
[INFO] step=1400, x=3.992, psi=4.956e+08
[INFO] step=1600, x=5.988, psi=1.831e+14
[RESULT] diverged at x=5.993
psi(end) = 1.831e+14
[INFO] Data saved to schrodinger_E0.6_diverged.xlsx

この場合、matplotlib のプロットウィンドウには、波動関数が発散していく様子が対数スケールで表示され、発散開始点を示す赤い点線が描画されます。psi(end) の値が非常に大きいことが、発散を示しています。

例3: ゼロ境界条件を使用する例 (\(E=0.5\) Ha)

ゼロ境界条件 (--bc zero) を使用して基底状態のエネルギー \(E=0.5\) Ha を計算します。この場合、初期微分値 --eps の指定が必要です。

python schrodinger1d.py --E 0.5 --bc zero --eps 1e-6 --L 10 --nx 2000 --outfile schrodinger_E0.5_zero_bc.xlsx

期待される実行結果 (コンソール出力の一部):

=== 1D Schrodinger Solver (atomic units) ===
E = 0.5 Ha,  L = 10.0,  nx = 2000
Boundary condition: zero
[INFO] step=200, x=-8.004, psi=1.127e-07
[INFO] step=400, x=-6.008, psi=7.514e-07
...
[RESULT] integration completed successfully
psi(end) = 1.127e-07
[INFO] Data saved to schrodinger_E0.5_zero_bc.xlsx

こちらも基底状態の波動関数が得られますが、--eps の値によっては波動関数が全体的にスケーリングされることに注意してください。これは、シューティング法では線形方程式の解の任意のスケーリング因子が許容されるためです。