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

プログラムの動作

H1s-HF-LDA.py は、Hartree-Fock-Slater (HFS) 理論に基づき、水素様1s軌道の電子エネルギー準位を計算するPythonプログラムです。特に、SlaterのX-\(\alpha\)交換相関ポテンシャルを考慮し、変分パラメータ \(k_a\) を持つSlater型軌道 (exp(-ka * Z * r)) を用いて全エネルギーを最小化します。

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

  • 電子数 Ne と動径関数パラメータ \(k_a\) の値を指定して、全エネルギー（運動エネルギー、核-電子引力、電子-電子反発、X-\(\alpha\)交換相関エネルギー）を計算します。

  • \(k_a\) を変分パラメータとして、scipy.optimize.minimize を用いて全エネルギーを最小化し、最適な \(k_a\) の値と最低エネルギーを探索します。

  • \(k_a\) または Ne の範囲をスイープし、それぞれのパラメータに対する全エネルギーの変化を標準出力に表示し、matplotlib を用いてグラフをプロットします。

  • デバッグモードでは、動径分布関数 \(R(r)\)、電荷分布 \(Q(r)\)、および各種ポテンシャル（核によるポテンシャル \(U_Z(r)\)、電子-核相互作用ポテンシャル \(U_{eZ}(r)\)、電子-電子相互作用ポテンシャル \(U_{e\rho}(r)\)、X-\(\alpha\)ポテンシャル \(U_{X\alpha}(r)\)）をグラフで可視化します。

このプログラムは、経験的なパラメータ \(k_a\) を用いたSlater型軌道における水素様原子のエネルギー計算を通じて、X-\(\alpha\)交換相関ポテンシャルを取り入れた電子間の相互作用を考慮し、変分原理によりエネルギーを最小化する最適な \(k_a\) を見つけることを目的としています。

原理

本プログラムは、Slater型軌道 (STO) を試行関数として採用し、HFS理論に基づく全エネルギーを計算し、変分原理によってそのエネルギーを最小化します。

1. Slater型軌道 (STO)

1s軌道の動径関数 \(R_{1s}(r)\) は以下のように定義されます。

\[ R_{1s}(r) = \frac{1}{\sqrt{\pi}} (k_a Z)^{3/2} e^{-k_a Z r} \]

ここで、\(k_a\) は変分パラメータ、\(Z\) は原子番号、\(r\) は動径座標です。この動径関数は3次元空間で規格化されており、\(\int_0^\infty 4\pi r^2 R_{1s}(r)^2 dr = 1\) を満たします。 電子密度 \(\rho(r)\) は動径関数の二乗として定義されます。

\[ \rho(r) = R_{1s}(r)^2 \]

2. 全エネルギーの構成要素

全エネルギー \(E_{tot}\) は、運動エネルギー \(T\)、核-電子間引力 \(U_{eZ}\)、電子-電子間反発 \(U_{ee}\)、およびX-\(\alpha\)交換相関エネルギー \(U_{X\alpha}\) の和として表されます。

\[ E_{tot} = T + U_{eZ} + U_{ee} + U_{X\alpha} \]

各項の計算は以下の通りです。

  • 運動エネルギー (\(T\)): 電子の運動エネルギーは、試行関数を用いて以下のように計算されます。

    \[ T = -\frac{1}{2} \int R(r) \left( \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dR(r)}{dr}\right) \right) 4\pi r^2 dr \]

    これはコード内の calT 関数で数値積分されます。

  • 核-電子間引力 (\(U_{eZ}\)): 原子核の電荷 \(Z\) と電子の電荷分布 \(\rho(r)\) の間のクーロン引力です。

    \[ U_{eZ} = -Z \int_0^\infty \frac{\rho(r)}{r} 4\pi r^2 dr \]

    これはコード内の calUeZ 関数で数値積分されます。

  • 電子-電子間反発 (\(U_{ee}\)): 電子の電荷分布 \(\rho(r)\) 間のクーロン反発です。

    \[ U_{ee} = \frac{1}{2} \iint \frac{\rho(\mathbf{r}_1)\rho(\mathbf{r}_2)}{|\mathbf{r}_1-\mathbf{r}_2|} d^3\mathbf{r}_1 d^3\mathbf{r}_2 = \int_0^\infty U_{e\rho}(r) \rho(r) 4\pi r^2 dr \]

    ここで、\(r\) における電子の電荷分布によるポテンシャル \(U_{e\rho}(r)\) は、内部電荷 \(Q(r)\) を用いて以下のように表されます。

    \[ U_{e\rho}(r) = \frac{Q(r)}{r} + \int_r^\infty \frac{\rho(r')}{r'} 4\pi r'^2 dr' \]

    内部電荷 \(Q(r)\) は、\(r\) までの電荷の総和です。

    \[ Q(r) = \int_0^r \rho(r') 4\pi r'^2 dr' \]

    これらの計算は calUrho および calUee 関数で数値積分されます。

  • X-\(\alpha\)交換相関エネルギー (\(U_{X\alpha}\)): SlaterのX-\(\alpha\)交換相関ポテンシャル \(U_{X\alpha}^{local}(r)\) を用いて計算されます。

    \[ U_{X\alpha}^{local}(r) = -3\alpha \left(\frac{3}{4\pi}\rho(r)\right)^{1/3} \]

    ここで、\(\alpha\) は交換係数で、デフォルト値は \(2/3\) です。全交換相関エネルギー \(U_{X\alpha}\) は、この局所ポテンシャルと電子密度を積分して求められます。

    \[ U_{X\alpha} = \int_0^\infty U_{X\alpha}^{local}(r) \rho(r) 4\pi r^2 dr \]

    これらの計算は calLocalUXa および calUXa 関数で数値積分されます。

3. 変分原理

全エネルギー \(E_{tot}\) は、\(k_a\) の関数として表現されます。変分原理に基づき、scipy.optimize.minimize 関数を用いて \(k_a\) を最適化し、\(E_{tot}\) を最小化します。この最小エネルギーを与える \(k_a\) が最適なパラメータとなります。

4. 単位変換

計算されたエネルギーは原子単位系 (Hartree) に対応する形で計算され、最終的に HartreeToeV = 27.2116 を掛けることで電子ボルト (eV) に変換されて出力されます。

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

このプログラムは、以下の非標準Pythonライブラリを使用しています。

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

  • scipy: 科学技術計算ライブラリ。特に数値積分 (integrate.quad)、内挿 (interpolate.interp1d)、最適化 (optimize.minimize) に使用されます。

  • matplotlib: グラフ描画ライブラリ。特に pyplot モジュールが使用されます。

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

pip install numpy scipy matplotlib

必要な入力ファイル

H1s-HF-LDA.py は、コマンドライン引数のみで動作します。特定の入力ファイルは必要ありません。

生成される出力ファイル

  • 計算結果（各エネルギー項、全エネルギーなど）は、実行中に標準出力（コンソール）にテキスト形式で表示されます。

  • グラフ表示オプション (g モード) が有効な場合、matplotlib によって生成されたプロットが画面に表示されます。これらのグラフはプログラムによってファイルとして自動的に保存されることはありませんが、表示されたウィンドウからユーザーが手動で保存することができます。

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

プログラムの実行には、以下の形式でコマンドライン引数を与えます。

Usage: Variables in () are optional
  (i) python H1s-HF-LDA.py mode Z ka Ne
      mode: combination of the following key characters
               d: debug mode, plot fundamental graphs
               v: add variational calculations
               e: output based on energy (default: based on E 1s eigen value)
               k: sweep ka
               n: sweep Ne
               g : plot graph
     ex: python H1s-HF-LDA.py nvg 1.0 1.0 1.0
     ex: python H1s-HF-LDA.py k 1.0 1.0 1.0

  • mode: プログラムの動作モードを指定する1文字以上の文字列。複数のモードを組み合わせることも可能です（例: nvg）。

    • d: デバッグモード。基本的なグラフをプロットし、詳細な情報を出力します。

    • v: 変分計算を追加します。\(k_a\) を最適化して最低エネルギーを探索します。

    • e: 全エネルギー \(E_{tot}\) を基準に出力します（デフォルトでは1s軌道の固有エネルギーに基づきます）。

    • k: \(k_a\) の値をスイープし、エネルギーの変化を調べます。

    • n: Ne (電子数) の値をスイープし、エネルギーの変化を調べます。

    • g: グラフをプロットして表示します。このモードは他のスイープ/デバッグモードと組み合わせて使用します。

  • Z: 原子番号（デフォルト: 1.0）

  • ka: Slater型軌道の指数部係数（デフォルト: 1.0）

  • Ne: 電子数（デフォルト: 1.0）

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

1. デバッグモードで基本情報を表示し、グラフをプロットする

原子番号 \(Z=1.0\)、指数部係数 \(k_a=1.0\)、電子数 \(Ne=1.0\) の水素様原子について、デバッグモードで実行し、動径分布関数やポテンシャルなどのグラフを表示します。

python H1s-HF-LDA.py dg 1.0 1.0 1.0

実行結果の説明: 標準出力には、指定されたパラメータでの各エネルギー項（運動エネルギー、核-電子引力、電子-電子反発、X-\(\alpha\)交換相関エネルギー）と全エネルギーがeV単位で表示されます。また、matplotlib ウィンドウが開き、動径分布関数 \(R(r)\)、電荷分布 \(Q(r)\)、および各種ポテンシャル（\(U(Z)\)、\(-U(Z-\rho)\)、\(U(\rho)\)、\(U(Xa)\)）のグラフが表示されます。プログラムはグラフウィンドウが表示された状態で一時停止し、Enterキーを押すと終了します。

2. \(k_a\) をスイープしてエネルギーの変化を調べ、グラフをプロットする

原子番号 \(Z=1.0\)、電子数 \(Ne=1.0\) の水素様原子について、\(k_a\) の値を広範囲にわたって変化させながら、全エネルギーを計算し、その結果をプロットします。

python H1s-HF-LDA.py kg 1.0 1.0 1.0

実行結果の説明: 標準出力には、\(k_a\) の各値に対する各エネルギー項と全エネルギーが表形式で出力されます。matplotlib ウィンドウには、\(k_a\) の値に対する全エネルギーのグラフが表示され、解析解（水素原子の場合は -13.6 eV）との比較も示されます。プログラムはグラフウィンドウが表示された状態で一時停止し、Enterキーを押すと終了します。

3. 電子数 Ne をスイープし、\(k_a\) を変分最適化しながらエネルギーを調べ、グラフをプロットする

原子番号 \(Z=1.0\)、指数部係数 \(k_a=1.0\) (初期値) の水素様原子について、電子数 Ne を変化させながら、各 Ne の値に対して \(k_a\) を変分的に最適化し、そのときの最低エネルギーと最適化された \(k_a\) の値をプロットします。

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

実行結果の説明: 標準出力には、まず \(k_a\) を最適化しない場合の Ne スイープ結果が表示され、続いて \(k_a\) を最適化した場合の Ne スイープ結果が表形式で出力されます。scipy.optimize.minimize の各イテレーションでの \(k_a\) とエネルギーの値も出力されます。matplotlib ウィンドウには、Ne の値に対する全エネルギー（最適化あり・なし）と、最適化された \(k_a\) の値のグラフが表示されます。プログラムはグラフウィンドウが表示された状態で一時停止し、Enterキーを押すと終了します。