H1s-HF-LDA.py 技術ドキュメント

プログラムの動作

H1s-HF-LDA.py は、水素様原子の1s軌道における電子のエネルギー準位を計算するためのPythonプログラムです。特に、SlaterのX-alphaポテンシャルを用いた局所密度近似（LDA）の枠組みで、運動エネルギー、原子核-電子間引力、電子-電子間反発（Hartree項）、および交換相関エネルギーを評価します。

このプログラムの主な機能は以下の通りです。

波動関数と電子密度の計算: Slater型の1s軌道波動関数（動径関数）を用いて、電子密度分布を算出します。
エネルギー項の計算: 運動エネルギー、原子核と電子間の静電ポテンシャルエネルギー、電子間の静電ポテンシャルエネルギー（Hartree項）、およびSlaterのX-alpha交換ポテンシャルエネルギーをそれぞれ数値積分により計算します。
全エネルギーの算出: 上記のエネルギー項の和として全エネルギー $E_{tot}$ を求めます。
パラメータ掃引と変分計算:
- 波動関数のスケール因子 $k_a$ や電子数 $N_e$ を変化させてエネルギーの変化を観察する掃引機能を提供します。
- scipy.optimize.minimize を用いて、全エネルギー $E_{tot}$ を最小化する最適な $k_a$ 値を変分法により探索する機能が実装されています。
結果の出力と可視化: 計算結果は標準出力に表形式で表示され、オプションで matplotlib を用いたグラフとして可視化されます。

このプログラムは、多電子原子の電子構造計算の基礎となる要素を理解し、特にLDAに基づく交換相関ポテンシャルの概念とその実装を学習するのに役立ちます。

原理

このプログラムは、水素様1s軌道を対象とした簡易的な電子構造計算に基づいています。電子の記述には、Slater型の単一1s軌道波動関数が用いられ、電子間の相互作用はHartree近似とSlaterのX-alpha交換ポテンシャルによって考慮されます。

1. 波動関数と電子密度

プログラムでは、1s軌道の動径関数 $R(r)$ を以下のSlater型の関数で近似します。 $$ R(r) = C \cdot (k_a Z)^{3/2} e^{-k_a Z r} $$ ここで、$C = 2.0 / \sqrt{4\pi}$ は、全空間で波動関数が正規化されるように調整された定数です。$Z$ は原子核の電荷（原子番号）、$k_a$ は波動関数の指数部の係数（変分パラメータ）、$r$ は原子核からの距離です。

電子密度 $\rho(r)$ は、この波動関数の絶対値の2乗として定義されます。 $$ \rho(r) = |\psi(r)|^2 = |R(r)|^2 $$ ここで、$|\psi(r)|^2$ は pi4 * r * r * rho(r) の積分が1になるように正規化されています。

動径内の電荷 $Q(r)$ は、半径 $r$ までの電子密度を積分したものです。 $$ Q(r) = \int_0^r 4\pi r'^2 \rho(r') dr' $$ ただし、プログラム内の `build_Qr` 関数では、`Ne * Q` として電子数 $N_e$ を考慮した累積電荷が計算され、これが補間関数 qfunc を通じて利用されます。

2. エネルギー項

全電子エネルギー $E_{tot}$ は、以下の4つの主要な項の和として計算されます。 $$ E_{tot} = T + U_{eZ} + U_{ee} + U_{Xa} $$

各項の計算方法は以下の通りです（エネルギーは原子単位（Hartree）で計算され、最終的にeVに変換されます）。

運動エネルギー $T$: 1s軌道の運動エネルギーは、古典的な電子軌道とは異なり量子力学的な寄与を持ちます。プログラムでは、動径関数 $R(r)$ を用いて以下の形で計算されます。 $$ T = \frac{1}{2} (k_a Z)^2 \quad \text{[Hartree]} $$ この式は、プログラム内の `calT` 関数で数値積分により計算される `$ \frac{1}{2} \int_0^\infty 4\pi r^2 |\psi(r)|^2 \left( -\frac{1}{2}\nabla^2 \right) |\psi(r)| dr $` の結果と一致します。
原子核-電子間ポテンシャルエネルギー $U_{eZ}$: 電子と原子核間の引力によるエネルギーです。 $$ U_{eZ} = \int 4\pi r^2 \rho(r) \left(-\frac{Z}{r}\right) dr = -Z (k_a Z) \quad \text{[Hartree]} $$ これは、calUeZ 関数で数値積分により計算されます。
電子-電子間静電ポテンシャルエネルギー $U_{ee}$ (Hartree項): 電子間のクーロン反発によるエネルギーです。プログラムでは、calUrho 関数が電子が感じる静電ポテンシャル $U_{potential}(r)$ を計算し、calUee 関数でそのポテンシャルと電子密度を積分します。 calUrho は以下の形のポテンシャルを計算します。 $$ U_{potential}(r) = N_e \frac{Q_{normalized}(r)}{r} + \int_r^\infty 4\pi r' |\psi(r')|^2 dr' $$ ここで $Q_{normalized}(r) = \int_0^r 4\pi r''^2 |\psi(r'')|^2 dr''$ です。そして `calUee` は `$ \int_0^\infty 4\pi r^2 \rho(r) U_{potential}(r) dr $` を計算します。最終的な $U_{ee}$ 項は、プログラムの mode 引数に応じてスケーリングされます。
- 'e' が mode に含まれる場合 (全エネルギー計算): $U_{ee}$ は $ N_e^2 $ 倍されます。
- 上記以外の場合: $U_{ee}$ は $ N_e $ 倍されます。
Slater X-alpha交換ポテンシャルエネルギー $U_{Xa}$: 電子間の交換相互作用を局所密度近似（LDA）で近似したものです。プログラムでは、calLocalUXa 関数が局所的な交換ポテンシャル $U_{Xa}(r)$ を計算します。 $$ U_{Xa}(r) = -3\alpha \left( \frac{3}{4\pi}\rho(r) \right)^{1/3} $$ ここで $\alpha = 2/3$ はSlaterのX-alpha係数です。 `calUXa` 関数では、$ \int_0^\infty 4\pi r^2 \rho(r) U_{Xa}(r) dr $ が計算されます。最終的な $U_{Xa}$ 項は、プログラムの mode 引数に応じてスケーリングされます。
- 'e' が mode に含まれる場合 (全エネルギー計算): $U_{Xa}$ は $ N_e^{4/3} $ 倍されます。
- 上記以外の場合: $U_{Xa}$ は $ N_e^{1/3} $ 倍されます。

3. 変分原理による最適化

'v' オプションが指定された場合、プログラムは scipy.optimize.minimize を使用して、全エネルギー $E_{tot}$ を最小化する最適な $k_a$ 値を探索します。これは、試行波動関数に含まれるパラメータを最適化することで、基底状態のエネルギーに最も近い値を見つける変分原理に基づいています。

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

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

numpy: 数値計算を効率的に行うための基本的なライブラリです。
scipy: 科学技術計算のためのライブラリで、特に本プログラムでは数値積分 (integrate.quad)、補間 (interpolate.interp1d)、および最適化 (optimize.minimize) の機能を利用しています。
matplotlib: グラフ描画のためのライブラリです。

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

pip install numpy scipy matplotlib

必要な入力ファイル

H1s-HF-LDA.py は、外部の入力ファイルを必要としません。必要な全てのパラメータ（原子番号 $Z$、電子数 $N_e$、スケール因子 $k_a$、モードなど）は、コマンドライン引数として直接プログラムに渡されるか、プログラム内部で初期値として設定されています。

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

本プログラムは、計算結果をファイルとして保存することはありません。計算の数値結果は標準出力（コンソール）に表形式で表示されます。また、グラフ描画オプション ('g' を mode に含む) が指定されている場合、matplotlib によって生成されたグラフがGUIウィンドウとして表示されます。

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

H1s-HF-LDA.py は以下の形式でコマンドラインから実行します。

python H1s-HF-LDA.py mode Z ka Ne

mode: 実行モードを指定する1文字以上のキー文字の組み合わせです。
- d: デバッグモード。基本的なエネルギー計算と、波動関数、電荷分布、各種ポテンシャルのグラフを表示します。
- v: 変分計算を追加します。$k_a$ を最適化して全エネルギーを最小化します。
- e: 全エネルギー ($E_{tot}$) を出力の基準とします（デフォルトでは1s固有値としてのエネルギー $E_{1s}$ が基準）。これにより、$U_{ee}$ や $U_{Xa}$ 項の $N_e$ 依存性が変わります。
- k: $k_a$ パラメータを掃引し、エネルギーの変化を調べます。
- n: $N_e$ (電子数) を掃引し、エネルギーの変化を調べます。
- g: グラフをプロットします。このモードは、d, k, n と組み合わせて使用します。
Z: 原子番号（浮動小数点数）。デフォルトは 1.0。
ka: 1s動径関数の指数部の係数（浮動小数点数）。デフォルトは 1.0。
Ne: 電子数（浮動小数点数）。デフォルトは 1.0。

これらの引数のうち、Z, ka, Ne は省略可能で、省略された場合はプログラム内部のデフォルト値が使用されます。

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

1. 電子数 ($N_e$) を掃引し、変分計算とグラフ表示を行う例

電子数 $N_e$ を変化させながらエネルギーを計算し、さらに $k_a$ を最適化してグラフ表示を行います。

実行コマンド:

python H1s-HF-LDA.py nvg 1.0 1.0 0.5

nvg:
- n: 電子数 ($N_e$) を掃引します。プログラム内部で定義された Nearray に従って $N_e$ が変化します。
- v: 変分計算を有効にします。各 $N_e$ に対して最適な $k_a$ を探索し、エネルギーを最小化します。
- g: 掃引結果をグラフで表示します。
1.0: 原子番号 $Z=1.0$ (水素)。
1.0: 初期の $k_a=1.0$。
0.5: 初期または中心となる $N_e=0.5$（n モードではこれが掃引範囲の中心付近として使われる）。

実行結果の説明: 標準出力には、$k_a$ が最適化されていない場合と、変分法により最適化された場合の各 $N_e$ 値での運動エネルギー ($T$)、核-電子間ポテンシャルエネルギー ($U_{eZ}$)、電子-電子間ポテンシャルエネルギー ($U_{ee}$)、X-alpha交換ポテンシャルエネルギー ($U_{Xa}$)、および全エネルギー ($E_{tot}$ または $E_{1s}$) がeV単位で表形式で出力されます。 ELabel はデフォルトで E 1s です。'e' モードが指定されている場合は Etot となります。同時に、matplotlib により、電子数 $N_e$ に対するエネルギー ($E_{1s}$ または $E_{tot}$) および最適化された $k_a$ の値のプロットが表示されます。エネルギーのグラフには、$N_e$ 周りの放物線近似も合わせて表示されます。

2. スケール因子 ($k_a$) を掃引し、グラフ表示を行う例

$k_a$ を変化させながらエネルギーを計算し、結果をグラフ表示します。