金属の有限温度電子密度計算スクリプト

プログラムの動作

N-integration-metal.py は、金属の有限温度 \(T\) における電子密度 \(N\) を計算するためのPythonスクリプトです。このプログラムの主な目的は、異なる数値積分手法（scipy.integrate.quad と scipy.integrate.romberg）の精度と計算効率を比較することです。

プログラムは以下の主要な機能を提供します。

  • 状態密度関数 (DOS) の定義: 3次元自由電子モデルにおける状態密度を簡略化した \(D(E) = \sqrt{E}\) を使用します。

  • フェルミ・ディラック分布関数 (Fermi-Dirac Distribution Function) の定義: 有限温度における電子のエネルギー占有確率 \(f(E, T, E_F)\) を計算します。

  • 数値積分の実行: 定義された状態密度関数とフェルミ・ディラック分布関数の積を被積分関数とし、Scipyライブラリの integrate.quad (ガウス・ルジャンドル積分) と integrate.romberg (ロンバーグ積分) を用いて数値積分を実行します。

  • 複数積分範囲での比較: フェルミエネルギー \(E_F\) 周辺を含む3つの異なる積分範囲で計算を行い、各手法の結果と実行時間を比較します。

  • パラメータの指定: 温度 \(T\) とフェルミエネルギー \(E_F\) をコマンドライン引数で指定可能です。

このスクリプトは、固体物理学における電子の統計的振る舞いを理解し、数値計算手法の特性を評価するためのツールとして機能します。

原理

このプログラムは、金属中の電子密度 \(N\) が、状態密度関数 \(D(E)\) とフェルミ・ディラック分布関数 \(f(E, T, E_F)\) の積をエネルギー \(E\) で積分することによって求められるという物理的原理に基づいています。

\[N = \int_0^\infty D(E) f(E, T, E_F) dE\]

状態密度関数 (Density of States, DOS)

プログラムでは、3次元自由電子モデルにおける状態密度を簡略化した関数として、エネルギー \(E\) に依存する以下の形式を使用しています。

\[D(E) = \sqrt{E}\]

フェルミ・ディラック分布関数 (Fermi-Dirac Distribution Function)

電子が特定のエネルギー状態を占有する確率を表す関数で、以下の式で与えられます。

\[f(E, T, E_F) = \frac{1}{\exp\left(\frac{E - E_F}{k_B T}\right) + 1}\]

ここで、\(E\) は電子のエネルギー、\(T\) は絶対温度、\(E_F\) はフェルミエネルギー、\(k_B\) はボルツマン定数です。 プログラムでは、特殊なケースとして \(T=0.0\) K の場合は以下の条件分岐が適用されます。

\[\begin{split}f(E, 0, E_F) = \begin{cases} 1.0 & (E < E_F) \\ 0.0 & (E \ge E_F) \end{cases}\end{split}\]

被積分関数

上記の状態密度関数とフェルミ・ディラック分布関数の積 \(D(E) \cdot f(E, T, E_F)\) が、電子密度を計算するための被積分関数となります。

数値積分手法

プログラムでは、以下のScipyライブラリの数値積分関数が使用されます。

  • scipy.integrate.quad: 適応型ガウス・ルジャンドル求積法に基づいています。高い精度と信頼性を提供し、広範囲な関数に対して効果的です。

  • scipy.integrate.romberg: ロンバーグ積分法に基づいています。これは台形公式から導かれる補外法であり、比較的滑らかな被積分関数に対して効率的かつ高精度な結果をもたらします。

物理定数

プログラム内で定義されている主な物理定数は以下の通りです。

  • pi: 円周率 (\(3.14159265358979323846\))

  • h: プランク定数 (\(6.6260755 \times 10^{-34}\) J s)

  • hbar: ディラック定数（換算プランク定数、\(1.05459 \times 10^{-34}\) J s)

  • c: 光速 (\(2.99792458 \times 10^8\) m/s)

  • e: 電気素量 (\(1.60218 \times 10^{-19}\) C)

  • kB: ボルツマン定数 (\(1.380658 \times 10^{-23}\) J K\(^{-1}\)) また、ekBT はエネルギーを \(E / (k_B T)\) の形に変換するための係数として使用されます。

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

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

  • numpy: 高度な数値計算を効率的に行うための基本的なライブラリです。配列操作や数学関数を提供します。

  • scipy: 科学技術計算のためのライブラリで、特に scipy.integrate モジュールが数値積分のために使用されます。

  • matplotlib: (プログラム内ではインポートされていますが、現状では直接的な描画処理は行われていません。将来的な機能拡張のために含まれている可能性があります。) データ可視化のためのライブラリです。

これらのライブラリは、Pythonのパッケージマネージャ pip を使用してインストールできます。

pip install numpy scipy matplotlib

必要な入力ファイル

このプログラムは、実行時に特定の入力ファイルを必要としません。 温度 \(T\) とフェルミエネルギー \(E_F\) の値は、コマンドライン引数として与えることができます。引数が指定されない場合、プログラム内のデフォルト値（\(T = 300.0\) K, \(E_F = 5.0\) eV）が使用されます。

生成される出力ファイル

このプログラムは、計算結果をファイルに保存しません。 すべての計算結果は標準出力（コンソール）に直接表示されます。出力には、設定されたパラメータ、状態密度関数 \(D(E)\) の積分の解析解と数値解、そして被積分関数 \(D(E)f(E,T,E_F)\) の数値積分結果（quad および romberg それぞれの積分値と推定誤差）、および各計算の実行時間が含まれます。

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

N-integration-metal.py は、以下の基本的な形式でコマンドラインから実行されます。

python N-integration-metal.py [温度T] [フェルミエネルギーEF]

  • [温度T] (オプション):

    • 絶対温度をケルビン (K) で指定する浮動小数点数です。

    • この引数を省略した場合、デフォルト値の \(300.0\) K が使用されます。

  • [フェルミエネルギーEF] (オプション):

    • フェルミエネルギーを電子ボルト (eV) で指定する浮動小数点数です。

    • この引数を省略した場合、デフォルト値の \(5.0\) eV が使用されます。

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

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.453559924999298
quad   : ret= (7.453559925000001, 8.275529177579628e-14)
romberg: ret= (7.453559925000001, 7.82823611310573e-15)

(1) Integration range:  0.0  -  5.155093717204992  eV
quad   : ret= (7.472718919782559, 8.296580556209212e-14)   time= 0.054737091064453125  s
romberg: ret= (7.472718919782558, 2.083315750868884e-14)   time= 0.07639193534851074  s

(2) Integration range:  0.0  -  4.844906282795008  eV
Analytical integration for De(E) from  0.0  to  4.844906282795008  =  7.245842103417758
quad   : ret= (7.245842103417758, 8.044005697274094e-14)   time= 0.05436301231384277  s
romberg: ret= (7.245842103417758, 7.86801915474667e-15)   time= 0.07635688781738281  s

(3) Integration range:  4.844906282795008  -  5.155093717204992  eV
quad   : ret= (0.2268768163648003, 2.5186026405230354e-15)   time= 0.05445218086242676  s
romberg: ret= (0.2268768163648003, 1.4878267232304856e-16)   time= 0.07664608955383301  s

この出力では、設定された \(T\) と \(E_F\) の値、状態密度関数 \(D(E)=\sqrt{E}\) の積分結果（解析解と quad/romberg による数値解）、そして3つの異なる積分範囲における被積分関数 \(D(E)f(E,T,E_F)\) の数値積分結果（quad/romberg による積分値と推定誤差）およびそれぞれの計算時間が示されています。ret= の後の最初の数値が積分値、2番目の数値が推定誤差です。

2. 特定の温度とフェルミエネルギーでの実行

例えば、温度を \(100\) K、フェルミエネルギーを \(3.5\) eV に設定して実行する場合。

python N-integration-metal.py 100 3.5

実行結果の例:

T =  100.0 k   EF =  3.5  eV

Analytical integration for De(E) from  0  to  3.5  =  4.136067098499252
quad   : ret= (4.1360670985, 4.591024502391629e-14)
romberg: ret= (4.1360670985, 4.414002497647246e-15)

(1) Integration range:  0.0  -  3.551717887508762  eV
quad   : ret= (4.170682229424619, 4.629392686861512e-14)   time= 0.0543668270111084  s
romberg: ret= (4.170682229424618, 4.414704021237937e-15)   time= 0.07613587379455566  s

(2) Integration range:  0.0  -  3.448282112491238  eV
Analytical integration for De(E) from  0.0  to  3.448282112491238  =  4.101454580062406
quad   : ret= (4.101454580062407, 4.552807897103448e-14)   time= 0.05421590805053711  s
romberg: ret= (4.101454580062407, 4.412850785638543e-15)   time= 0.07611894607543945  s

(3) Integration range:  3.448282112491238  -  3.551717887508762  eV
quad   : ret= (0.06922764936221469, 7.68532402167727e-16)   time= 0.05439305305480957  s
romberg: ret= (0.06922764936221469, 7.502859940733562e-17)   time= 0.07619285583496094  s

この例では、指定された温度とフェルミエネルギーに基づいた計算結果が表示されます。これにより、異なる物理条件下での電子密度の計算結果を比較し、各積分手法の挙動を確認できます。