N-integration-metal.py の技術ドキュメント

プログラムの動作

N-integration-metal.py は、金属の有限温度 $T$ における電子のエネルギー分布と総電子密度に関する数値積分の精度と計算効率を検証するPythonプログラムです。特に、自由電子ガスモデルにおける状態密度関数 $D(E) = \sqrt{E}$ とフェルミ・ディラック分布関数 $f(E)$ の積 $D(E)f(E)$ の積分に焦点を当てています。

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

定数定義: 物理定数（$\pi$, プランク定数 $h$, 換算プランク定数 $\hbar$, 光速 $c$, 素電荷 $e$, ボルツマン定数 $k_B$）を設定します。
パラメータ設定: デフォルトで温度 $T = 300.0$ K、フェルミエネルギー $E_F = 5.0$ eV を設定します。これらはコマンドライン引数で変更可能です。
関数定義:
- De(E): エネルギー $E$ に対する状態密度関数（$\sqrt{E}$ に比例）。
- fe(E, T, EF): エネルギー $E$ に対するフェルミ・ディラック分布関数。 $T=0$ K の特殊ケースも処理します。
- Defe(E, T, EF): 状態密度関数とフェルミ・ディラック分布関数の積。
数値積分: scipy.integrate モジュールから quad (適応型求積法) および romberg (ロンバーグ積分法) の2種類の数値積分手法を使用します。
精度と効率の比較: 複数の異なる積分範囲（$0 \sim E_F+dE$, $0 \sim E_F-dE$, $E_F-dE \sim E_F+dE$ など）において、各積分手法の結果と計算時間を比較します。一部のケースでは、解析解との比較も行い、数値積分の精度を確認します。
出力: 計算された積分結果、誤差推定（quad の場合）、および各積分手法の計算時間を標準出力に表示します。

このプログラムは、金属物理学や固体物理学における電子の振る舞いを理解し、数値計算手法の適用性を評価するための教育的・研究的なツールとして機能します。

原理

このプログラムは、金属中の電子のエネルギー分布を記述する量子統計力学と、その数値計算手法に基づいています。

状態密度関数 (Density of States) 自由電子ガスモデルにおいて、エネルギー $E$ における状態密度 $D(E)$ は、電子が占有可能なエネルギー準位の密度を表し、以下のようにエネルギーの平方根に比例します。 $$D(E) \propto \sqrt{E}$$ プログラム中の `De(E)` 関数はこの関係を実装しており、$E^{1/2}$ を返します。
フェルミ・ディラック分布関数 (Fermi-Dirac Distribution Function) 絶対零度以上の有限温度 $T$ において、電子がエネルギー準位 $E$ を占有する確率 $f(E, T, E_F)$ は、フェルミ・ディラック統計によって記述されます。 $$f(E, T, E_F) = \frac{1}{\exp\left(\frac{E - E_F}{k_B T}\right) + 1}$$ ここで、$E_F$ はフェルミエネルギー、$k_B$ はボルツマン定数、$T$ は絶対温度です。プログラム中の `fe(E, T, EF)` 関数はこの式を実装しており、$T=0$ K の場合（$E < E_F$ で $f(E)=1$、$E > E_F$ で $f(E)=0$）も特別に扱います。`ekBT` は $e / (k_B T)$ を表し、エネルギーを $k_B T$ で無次元化するために使用されます。
被積分関数と全電子密度 ある温度 $T$ における全電子密度 $N$ は、状態密度 $D(E)$ とフェルミ・ディラック分布関数 $f(E, T, E_F)$ の積を全エネルギー範囲で積分することによって得られます。 $$N = \int_0^\infty D(E) f(E, T, E_F) dE$$ プログラム中の `Defe(E, T, EF)` 関数は、この被積分関数 $D(E) \cdot f(E, T, E_F)$ を計算します。この積分は、温度によって変化する電子のエネルギー分布を反映した総電子数を表します。
数値積分手法 このプログラムでは、scipy.integrate モジュールから以下の2つの数値積分関数を使用して、上記の積分を実行します。
- scipy.integrate.quad: 適応型求積法（通常はGauss-Kronrod法）に基づき、高い精度と信頼性で積分を実行します。結果として積分の値と絶対誤差推定値を返します。
- scipy.integrate.romberg: ロンバーグ積分法は、複合台形公式やシンプソン公式の結果からリチャードソン外挿を適用して、積分の精度を向上させる手法です。
プログラムは、これらの手法を用いて上記の積分を異なる積分範囲で実行し、その結果の精度（既知の解析解との比較や、異なる手法間の比較）と計算時間を比較することで、数値積分の性能を評価します。

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

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

numpy: 高度な数値計算（特に配列操作）を効率的に行うための基盤ライブラリ。
scipy: 科学技術計算のためのライブラリ。特にscipy.integrateモジュールが数値積分のために使用されます。
matplotlib: グラフ描画などデータ可視化のためのライブラリ。ただし、このプログラムの現在のバージョンではインポートされていますが、直接使用はされていません。

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

pip install numpy scipy matplotlib

必要な入力ファイル

N-integration-metal.py プログラムは、実行時に特定の入力ファイルを必要としません。

プログラムの動作に必要なすべてのパラメータ（温度 $T$ とフェルミエネルギー $E_F$）は、プログラムコード内にデフォルト値が設定されており、必要に応じてコマンドライン引数として渡すことができます。

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

N-integration-metal.py プログラムは、実行時にいかなるファイルも生成したり保存したりしません。

すべての計算結果と実行情報は、標準出力（コンソール）に直接表示されます。出力される主な内容は以下の通りです。

設定された温度 ($T$) とフェルミエネルギー ($E_F$) の値。
状態密度関数 De(E) の解析的な積分結果（$T=0$ K の場合のフェルミエネルギーまでの積分）と、quad および romberg による数値積分結果。
異なる積分範囲（例：$0 \sim E_F + dE$, $0 \sim E_F - dE$, $E_F - dE \sim E_F + dE$）における、被積分関数 Defe(E, T, EF) の quad および romberg による数値積分結果。
各数値積分手法の計算にかかった時間（ncycle 回繰り返した場合の合計時間）。
quad 関数によって返される積分結果と、その誤差推定値（タプル形式で表示されます）。

これらの出力は、異なる数値積分手法の精度と計算効率を比較するための情報を提供します。

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

N-integration-metal.py は、以下の形式でコマンドラインから実行できます。プログラムの実行に必要な引数はオプションであり、指定しない場合はデフォルト値が使用されます。

基本形式:

python N-integration-metal.py [T] [EF]

[T]: 計算に使用する温度をケルビン (K) で指定します。これは浮動小数点数である必要があります。
- この引数を省略した場合、プログラムのデフォルト値である $300.0$ K が使用されます。
[EF]: 計算に使用するフェルミエネルギーを電子ボルト (eV) で指定します。これも浮動小数点数である必要があります。
- この引数を省略した場合、プログラムのデフォルト値である $5.0$ eV が使用されます。
- [T] を指定し、[EF] を省略することはできません。[T] を指定する場合は、[EF] も合わせて指定する必要があります。

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

例1: デフォルト値での実行

温度 $T = 300.0$ K、フェルミエネルギー $E_F = 5.0$ eV のデフォルト値を使用してプログラムを実行します。

python N-integration-metal.py

実行結果の例:

T =  300.0 k   EF =  5.0  eV

Analytical integration for De(E) from  0  to  5.0  =  7.453559924999299
quad   : ret= (7.453559924999298, 8.27552914170364e-14)
romberg: ret= 7.453559924999298

(1) Integration range:  0.0  -  5.155099999999999  eV
quad   : ret= (7.453559924999298, 8.27552914170364e-14)  time= 0.05267786979675293  s
romberg: ret= 7.453559924999298  time= 0.027003765106201172  s

(2) Integration range:  0.0  -  4.8449  eV
Analytical integration for De(E) from  0.0  to  4.8449  =  7.098808167733221
quad   : ret= (7.09880816773322, 7.881775836894046e-14)  time= 0.052824974060058594  s
romberg: ret= 7.09880816773322  time= 0.02699899673461914  s

(3) Integration range:  4.8449  -  5.155099999999999  eV
quad   : ret= (0.015024765660682259, 1.6682977539045766e-16)  time= 0.05432605743408203  s
romberg: ret= 0.015024765660682259  time= 0.0279998779296875  s

この出力は、設定された $T$ と $E_F$ を示し、De(E) 関数の解析解と数値積分結果を比較します。その後、Defe(E, T, EF) 関数を異なる積分範囲で quad と romberg を用いて積分した結果と、それぞれの計算時間を表示します。quad の結果は、積分の値とその誤差推定値のペア（タプル）として表示されます。

例2: 特定の温度とフェルミエネルギーを指定して実行

温度 $T = 10.0$ K、フェルミエネルギー $E_F = 2.5$ eV を指定してプログラムを実行します。

python N-integration-metal.py 10.0 2.5

実行結果の例:

T =  10.0 k   EF =  2.5  eV

Analytical integration for De(E) from  0  to  2.5  =  2.6352313626786937
quad   : ret= (2.6352313626786937, 2.9250020478142273e-14)
romberg: ret= 2.6352313626786937

(1) Integration range:  0.0  -  2.505175  eV
quad   : ret= (2.6352313626786937, 2.9250020478142273e-14)  time= 0.05199694633483887  s
romberg: ret= 2.6352313626786937  time= 0.027000904083251953  s

(2) Integration range:  0.0  -  2.494825  eV
Analytical integration for De(E) from  0.0  to  2.494825  =  2.624796075191427
quad   : ret= (2.624796075191427, 2.9135767690460394e-14)  time= 0.05200386047363281  s
romberg: ret= 2.624796075191427  time= 0.026999950408935547  s

(3) Integration range:  2.494825  -  2.505175  eV
quad   : ret= (0.0001099859253457581, 1.2210875954605912e-18)  time= 0.05200505256652832  s
romberg: ret= 0.0001099859253457581  time= 0.02699899673461914  s

この例では、コマンドラインで指定された $T=10.0$ K と $E_F=2.5$ eV を用いて計算が実行されます。デフォルト値の例と比較すると、フェルミエネルギーの変更に伴い、状態密度の解析積分結果や各数値積分結果が変化していることがわかります。特に低温では、フェルミ・ディラック分布関数の形状がシャープになるため、$E_F$ 近傍の積分値が大幅に小さくなっていることが見て取れます。