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

プログラムの動作

plot_wf.py は、水素原子の波動関数（軌道）を様々な形式で可視化するためのPythonスクリプトです。主に3次元空間における波動関数の振る舞いを静的なプロットやアニメーションとして表現します。このプログラムは、tkWavefunction_H（水素原子波動関数計算）および tkPlot3d（3Dプロットユーティリティ）というカスタムモジュールに強く依存しています。

主な機能は以下の通りです。

水素原子の軌道可視化: 主量子数 $n$、方位量子数 $l$、磁気量子数 $m$ で指定される水素原子の軌道について、波動関数の実数部、虚数部、絶対値、または絶対値の2乗を3次元空間にプロットします。
多様な描画モード:
- rplot / rsurface: 球面調和関数の振幅に基づく3D曲面プロット。
- ranim: 波動関数の位相変化を伴う3D曲面アニメーション。
- dot / dotc / dotanim: モンティカルロ法によるリジェクションサンプリングを用いた3Dドットプロット（確率密度に応じた点の分布）。
- iso3d / iso3danim: 波動関数の等値面（isosurface）プロットとそのアニメーション。
- Rnl: 動径波動関数 $R_{nl}(r)$ と動径確率密度分布 $P_{nl}(r)$ の2Dプロット。
- Ylm / Phi: 球面調和関数の角度依存性（$\theta$ または $\phi$）の2Dプロット。
- normalize: 波動関数の正規化のチェック（モンテカルロ積分による）。
汎用的な3D図形描画: トーラス、球の等値面、指数関数の等高線など、一般的な3D関数プロットもサポートしています。
アニメーション機能: 時間経過による波動関数の位相変化や、モンテカルロサンプリングにおける点の増加をアニメーションとしてGIF形式で保存できます。
柔軟な設定: コマンドライン引数を通じて、描画モード、表示する波動関数の種類、サンプリング数、プロット表示の有無、ファイル保存の有無などを設定できます。

このプログラムは、量子力学における水素原子の波動関数を視覚的に理解するための強力なツールとして機能します。

原理

plot_wf.py は、水素原子のシュレーディンガー方程式の解である波動関数 $\psi_{nlm}(r, \theta, \phi)$ を基に可視化を行います。

1. 水素原子の波動関数

水素原子の波動関数は、量子数 $n$（主量子数）、$l$（方位量子数）、$m$（磁気量子数）によって一意に決定され、動径部分 $R_{nl}(r)$ と角度部分 $Y_{lm}(\theta, \phi)$ の積で表されます。

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

動径波動関数 $R_{nl}(r)$: 原子核からの距離 $r$ に依存し、関連するラグランジュ多項式 (Associated Laguerre polynomials) を用いて表現されます。 $$ R_{nl}(r) = \sqrt{\left(\frac{2Z}{na_0}\right)^3 \frac{(n-l-1)!}{2n[(n+l)!]^3}} e^{-Zr/(na_0)} \left(\frac{2Zr}{na_0}\right)^l L_{n-l-1}^{(2l+1)}\left(\frac{2Zr}{na_0}\right) $$ ここで、$Z$ は原子番号、$a_0$ はボーア半径、$L_p^{(\alpha)}(x)$ は関連するラグランジュ多項式です。
球面調和関数 $Y_{lm}(\theta, \phi)$: 天頂角 $\theta$ と方位角 $\phi$ に依存し、ルジャンドル陪多項式 (Associated Legendre polynomials) と指数関数を用いて表現されます。 $$ Y_{lm}(\theta, \phi) = (-1)^m \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} P_l^m(\cos\theta) e^{im\phi} $$ ここで、$P_l^m(\cos\theta)$ は関連するルジャンドル陪多項式です。

プログラム内部では、scipy.special モジュールの sph_harm (球面調和関数), lpmv (ルジャンドル陪多項式), assoc_laguerre (関連するラグランジュ多項式) などの関数がこれらの計算に使用されています。

2. 可視化手法

3D曲面プロット (rplot, ranim, rsurface): 球面調和関数の絶対値や実数部などを動径方向の半径 $r$ として、球面座標からデカルト座標に変換してプロットします。 $$ x = r \sin\theta \cos\phi $$ y = r \sin\theta \sin\phi $$ z = r \cos\theta $$ `ranim` モードでは、$\psi_{nlm}(r, \theta, \phi, t) = \psi_{nlm}(r, \theta, \phi) e^{-i E_n t / \hbar}$ のように位相因子 $e^{-i E_n t / \hbar}$ を導入し、この位相の時間変化をアニメーションとして表現します。プログラムでは、位相変化を phase_m 引数として渡すことでシミュレートしています。
モンテカルロ・リジェクションサンプリング (dot, dotc, dotanim): 波動関数の絶対値の2乗 $| \psi_{nlm}(x,y,z) |^2$ は、電子が存在する確率密度を表します。dot モードでは、この確率密度に比例して3次元空間に点をサンプリングします。
1. 指定された範囲 ($[x_{min}, x_{max}]$, $[y_{min}, y_{max}]$, $[z_{min}, z_{max}]$) 内で一様な乱数 $(x, y, z)$ を生成します。
2. その点における確率密度 $P(x,y,z) = |\psi_{nlm}(x,y,z)|^2$ を計算します。
3. $P(x,y,z)$ の最大値 $M$ を事前に推定し、$[0, M]$ の範囲で一様乱数 $u$ を生成します。
4. もし $u < P(x,y,z)$ であれば、その点 $(x, y, z)$ を受け入れてプロットに追加します。これを nsamples 回繰り返します。これにより、確率密度の高い領域に多くの点が、低い領域には少ない点が描画され、軌道の形状が確率的に可視化されます。
等値面プロット (iso3d, iso3danim): 波動関数の実数部や絶対値がある特定の閾値と等しくなるような3次元空間上の点の集合を表面として描画します。scikit-image ライブラリの marching_cubes_lewiner のようなアルゴリズムが tkPlot3d モジュール内部で使われている可能性があります（コードからは直接確認できないが、等値面生成の一般的な手法）。
正規化の確認 (normalize): 波動関数は全空間で積分すると1になるように正規化されています。プログラムでは、動径部分について $$ \int_0^\infty 4\pi r^2 |R_{nl}(r)|^2 dr = 1 $$ を `scipy.integrate.quad` で数値積分したり、3次元全体について $$ \iiint |\psi_{nlm}(x,y,z)|^2 dx dy dz = 1 $$ をモンテカルロ積分 (integ_MC3d) で数値的に計算し、正規化が保たれているかを確認します。モンテカルロ積分は、サンプリングされた点の関数値の平均に体積をかけることで近似されます。

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

このプログラムを実行するには、以下の非標準Pythonライブラリが必要です。

numpy: 数値計算（配列操作、数学関数など）
scipy: 科学技術計算（積分、特殊関数など）
matplotlib: グラフ描画（2D/3Dプロット、アニメーションなど）
imageio: GIFアニメーションの保存（imagemagick ライターを使用する場合、内部で利用される可能性があります）

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

pip install numpy scipy matplotlib imageio

カスタムモジュール:

さらに、以下の2つのカスタムモジュールがこのプログラムと同じディレクトリに存在する必要があります。

tkWavefunction_H.py: 水素原子の波動関数に関する計算を行うモジュール。
tkPlot3d.py: 3Dプロットのユーティリティ関数を提供するモジュール。

これらのファイルが提供されていない場合、plot_wf.py は正常に動作しません。

必要な入力ファイル

plot_wf.py 自体は、コマンドライン引数でパラメータを受け取り、直接的な入力データをファイルから読み込むことはありません。

ただし、プログラムの実行には、前述の「必要な非標準ライブラリ」に加えて、以下のカスタムPythonモジュールが plot_wf.py と同じディレクトリに配置されている必要があります。

tkWavefunction_H.py
tkPlot3d.py

これらのモジュールが存在しない場合、ImportError が発生し、プログラムは実行できません。

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

プログラムの実行モードと設定（--fsave 1 の場合）に応じて、以下の種類のファイルが生成されます。

静止画像ファイル (.png)
- plot_wf_Rnl.png: 動径波動関数 $R_{nl}(r)$ と動径確率密度分布 $P_{nl}(r)$ の2Dグラフが複数表示された画像。
- plot_wf_{mode}_{functype}.png: rplot, dot, dotc, iso3d, rsurface, isosurface3d, contour3d などの静的プロットモードで生成される3Dグラフの画像。
  - 例: plot_wf_rplot_2pxa.png (2p軌道x成分絶対値の3Dプロット)
アニメーションファイル (.gif)
- plot_wf_{mode}_{functype}.gif: ranim, dotanim, iso3danim, 3dMC_animation などのアニメーションモードで生成されるGIFアニメーションファイル。
  - 例: plot_wf_ranim_2pz.gif (2p軌道z成分アニメーション)

{mode} には実行時に指定されたモード名が、{functype} には指定された波動関数のタイプ（例: 2pxa）が入ります。

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

plot_wf.py はコマンドライン引数で動作モードと波動関数の種類などを指定します。

基本的な書式は以下の通りです。

python plot_wf.py [mode] [functype] [nsamples] [--fplot <0|1>] [--fsave <0|1>]

mode: 実行する操作モードを指定します。
- 可視化モード: rplot, ranim, dot, dotc, dotanim, iso3d, iso3danim, rsurface, 3dMC_animation, contour3d, torus, isosurface3d
- 関数プロットモード: Rnl, Ylm, Phi
- 解析モード: normalize
functype: 描画対象の波動関数の種類を指定します。
- 形式は [n][l][m][タイプ] の組み合わせです。
  - n: 主量子数 (1, 2, 3...)
  - l: 方位量子数 (s, p, d, f...)
  - m: 磁気量子数 (x, y, z, 0, +-1, +-2...)
  - タイプ:
    - a: 絶対値 $| \psi |$
    - a2: 絶対値の2乗 $| \psi |^2$ (確率密度)
    - r: 実数部 $Re(\psi)$
    - i: 虚数部 $Im(\psi)$
- 例: 1s, 2pxa, 3dza, 2p0r
nsamples: (オプション) モンティカルロサンプリングを使用するモード (dot, dotanim, 3dMC_animation) におけるサンプル数を指定します。デフォルトは 30000 です。
--fplot <0|1>: (オプション) グラフを表示するかどうかを制御します。1 で表示（デフォルト）、0 で非表示。
--fsave <0|1>: (オプション) 生成されたグラフをファイルに保存するかどうかを制御します。1 で保存（デフォルト）、0 で非保存。

使用例のヒント:

General usage:
usage: plot_wf.py [-h] [--fplot {0,1}] [--fsave {0,1}] [mode] [functype] [nsamples]

Plot wave functions of H


positional arguments:
  mode        Mode of operation
  functype    Function type
  nsamples    Number of samples for Monte Carlo sampling

options:
  -h, --help  show this help message and exit
  --fplot {0,1}
              Flag to plot graphs (1: enable, 0: disable)
  --fsave {0,1}
              Flag to save figure file (1: enable, 0: disable)

Examples:
usage: python .\plot_wf.py [rplot|ranim|dot|dotc|dotanim|iso3d|iso3danim|rsurface] [100r|2pxa etc]
usage: python .\plot_wf.py [Rnl|Ylm|Theta|Phi|normalize]
usage: python .\plot_wf.py [3dMC_animation|contour3d|torus]

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

1. 2p軌道x成分の絶対値を3D曲面としてプロット

2pxa は、$n=2$, $l=1$ (p軌道), $m=x$ (磁気量子数1に相当), 絶対値 $| \psi |$ を意味します。

python plot_wf.py rplot 2pxa

実行結果の説明: 画面に2p_x軌道の絶対値 $| \psi |$ を3D曲面として描画したウィンドウが表示されます。波動関数の振幅に応じて色が変化し、軌道のロブ形状が視覚的に表現されます。ファイル保存が有効な場合 (--fsave 1)、plot_wf_rplot_2pxa.png という画像ファイルが生成されます。

2. 1s軌道の確率密度に応じた3Dドットプロットを生成（アニメーションなし、ファイル保存なし）

1s は、$n=1$, $l=0$ (s軌道), $m=0$ を意味します。デフォルトでは絶対値 a が適用されます。

python plot_wf.py dot 1s --nsamples 50000 --fsave 0

実行結果の説明: 画面に1s軌道の確率密度分布をモンテカルロ法でサンプリングした3Dドットプロットが表示されます。--nsamples 50000 により50000個の点が生成され、中心に近いほど点が密に分布し、球対称な形状が確認できます。--fsave 0 のため、画像ファイルは保存されません。

3. 2p軌道z成分の位相変化アニメーションをGIFファイルとして保存

2pz は、$n=2$, $l=1$ (p軌道), $m=z$ (磁気量子数0に相当) を意味します。

python plot_wf.py ranim 2pz --fplot 0 --fsave 1

実行結果の説明: このコマンドは、2p_z軌道の波動関数の位相が時間とともに変化する様子をアニメーションとして計算し、plot_wf_ranim_2pz.gif というGIFファイルを生成します。--fplot 0 により、プロットウィンドウは表示されず、バックグラウンドで処理が進みます。GIFファイルには、波動関数の複素位相変化に伴って3D曲面が回転する様子が収められます。

4. 動径波動関数と動径確率密度分布をプロット

python plot_wf.py Rnl

実行結果の説明: 様々な $n$ と $l$ の組み合わせについて、動径波動関数 $R_{nl}(r)$ と動径確率密度分布 $P_{nl}(r)$ の2Dグラフがそれぞれ9つずつ、合計2枚のウィンドウに表示されます。ファイル保存が有効な場合、plot_wf_Rnl.png と plot_wf_Pnl.png の2つの画像ファイルが生成されます。

5. 球面調和関数の角度依存性（$\theta$）をプロット

python plot_wf.py Ylm

実行結果の説明: $Y_{lm}(\theta, \phi)$ の $\theta$ 依存性（$\phi=0$ 固定）を示すグラフが1枚のウィンドウに表示されます。様々な $l$ と $m$ の組み合わせについて、$\theta$ に対する球面調和関数の値がプロットされ、角度的な形状の違いが確認できます。ファイル保存が有効な場合、plot_wf_Ylm_.png という画像ファイルが生成されます。