技術ドキュメント: benzene.py

プログラムの動作

benzene.py は、ベンゼン分子の電子状態をHückel近似モデルを用いて計算し、その結果であるエネルギー準位と状態密度（DOS）を可視化するPythonプログラムです。

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

Hückel行列の定義: ベンゼン分子の\(\pi\)電子系に対応するHückel行列をプログラム内で直接定義します。
固有値・固有ベクトルの計算: 定義されたHückel行列に対して、固有値分解を実行し、分子軌道のエネルギー準位（固有値）と分子軌道の係数（固有ベクトル）を算出します。
エネルギー準位のソートと縮重度計算: 算出したエネルギー準位を昇順にソートし、各準位が持つ縮重度（同じエネルギーを持つ準位の数）を特定します。
状態密度（DOS）の計算とブロードニング: 離散的なエネルギー準位を、ガウス関数でブロードニングすることで連続的な状態密度関数を計算します。これにより、各エネルギー準位の存在確率が連続的な分布として表現されます。
結果の可視化: Matplotlibライブラリを用いて、Hückel近似によるエネルギー準位図と、それに対応する状態密度グラフを一つのウィンドウに並べて表示します。

このプログラムは、簡略化されたHückel近似を通じて、共役系分子の電子構造がどのように離散的なエネルギー準位を持ち、それが連続的な状態密度としてどのように表現されるかを視覚的に理解するのに役立ちます。

原理

benzene.py は、化学におけるHückel近似に基づいています。Hückel近似は、共役系分子（特に\(\pi\)電子系）の電子構造を計算するための簡略化された分子軌道理論です。

Hückel近似の基本:
- \(\pi\)電子系のみを考慮し、\(\sigma\)電子は無視します。
- 原子軌道間の重なり積分は無視されます。
- クーロン積分（原子核と電子間の相互作用、\(H_{ii}\)）は、すべての炭素原子で同じ値 \(\alpha\) と仮定されます。プログラムでは、相対的なエネルギーを示すために \(\alpha=0\) と設定されています。
- 共鳴積分（隣接する原子間の結合、\(H_{ij}\)）は、結合している原子間では同じ値 \(\beta\) と仮定されます。プログラムでは、相対的なエネルギーを示すために \(\beta=1\) と設定されています。
- 非隣接原子間の共鳴積分は0とされます。
Hückel行列: これらの仮定により、Hückel近似は線形代数の固有値問題に帰着されます。N個の原子からなる\(\pi\)電子系の場合、Hückel行列はN×Nの行列となります。 benzene.py では、6つの炭素原子を持つベンゼン環のHückel行列 H が以下のように定義されています。
```
H = [[0, 1, 0, 0, 0, 1],
     [1, 0, 1, 0, 0, 0],
     [0, 1, 0, 1, 0, 0],
     [0, 0, 1, 0, 1, 0],
     [0, 0, 0, 1, 0, 1],
     [1, 0, 0, 0, 1, 0]]
```
この行列の対角要素が \(\alpha\)（0）、隣接する原子に対応する非対角要素が \(\beta\)（1）、それ以外の非対角要素が0です。行列の行と列は、ベンゼン環上の各炭素原子に対応しています。
エネルギー準位と固有ベクトル: Hückel行列 H の固有値問題は、以下のシュレーディンガー方程式の変形として表現されます。

\[ \mathbf{Hc} = E\mathbf{c} \]

ここで、\(\mathbf{H}\) はHückel行列、\(E\) は分子軌道のエネルギー（固有値）、\(\mathbf{c}\) は分子軌道を構成する原子軌道の係数ベクトル（固有ベクトル）です。 numpy.linalg.eig(H) 関数がこの固有値問題を解き、ベンゼンの6つの\(\pi\)電子軌道のエネルギー準位（固有値）と、各軌道の原子軌道係数（固有ベクトル）を算出します。
状態密度 (DOS): 状態密度は、特定のエネルギー範囲内に存在する電子状態の数を表す関数です。離散的なエネルギー準位（固有値）から連続的なDOSを構築するために、各エネルギー準位 \(E_i\) に対してガウス関数を適用し、それらを合計します。

\[ \text{DOS}(E) = \sum_{i} g_i \exp\left(-\frac{(E - E_i)^2}{2\sigma^2}\right) \]

ここで、\(g_i\) はエネルギー準位 \(E_i\) の縮重度（同じエネルギーを持つ状態の数）、\(\sigma\) はブロードニング幅（ガウス関数の広がり）です。この操作により、離散的なエネルギー準位がブロードニングされ、連続的なDOS曲線が得られます。プログラムでは、scipy.ndimage.gaussian_filter1d を使用して、このブロードニング処理を効率的に行っています。

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

benzene.py を実行するためには、以下の非標準Pythonライブラリが必要です。

numpy: 数値計算、特に多次元配列の操作や線形代数計算（固有値・固有ベクトルの計算）に利用されます。
matplotlib: グラフの描画、特にエネルギー準位図と状態密度グラフの作成に使用されます。
scipy: 科学技術計算ライブラリで、このプログラムではガウスフィルター処理 (gaussian_filter1d) に使用されます。
pprint: Pythonのデータ構造を見やすく整形して出力するためのモジュールですが、プログラムの主要機能ではありません。
collections: Pythonのコンテナデータ型を提供するモジュールで、このプログラムでは Counter クラスを使用してエネルギー準位の縮重度を計算します。

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

pip install numpy matplotlib scipy

pprint と collections は標準ライブラリに含まれているため、別途インストールは不要です。

必要な入力ファイル

benzene.py は、外部からの入力ファイルを必要としません。ベンゼン分子のHückel行列はプログラムのソースコード内に直接定義（ハードコード）されています。

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

benzene.py は、ファイルを生成して保存することはありません。

プログラムを実行すると、matplotlib によってグラフィカルなウィンドウが表示されます。このウィンドウには以下の2つのグラフが含まれます。

Hückel Approximation for Benzene (左側):
- ベンゼン分子のHückel近似による電子エネルギー準位を示します。各黒い横線が1つのエネルギー準位を表し、縮重度に応じて同じ高さに複数の線が描かれることがあります。縦軸はエネルギー（単位は「eV」と表示されますが、Hückel近似では無次元化された \(\beta\) 単位と解釈するのが一般的です）。
DOS (右側):
- 計算された状態密度（DOS）を示します。横軸はDOSの値、縦軸はエネルギーで、左側のエネルギー準位図と縦軸が共有されています。DOSグラフのピークは、エネルギー準位が存在する位置に対応しており、ブロードニングによって連続的な曲線となっています。

また、プログラムの実行中に、Hückel行列とソートされたエネルギー準位が標準出力（コンソール）に表示されます。

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

benzene.py は引数を必要としません。Pythonインタープリタを使用して直接実行します。

python benzene.py

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

以下のコマンドを実行します。

python benzene.py

実行結果:

コンソール出力: プログラムはまず、Hückel行列と計算されたエネルギー準位をコンソールに出力します。
```
Hückel行列（ベンゼン）:
array([[0, 1, 0, 0, 0, 1],
       [1, 0, 1, 0, 0, 0],
       [0, 1, 0, 1, 0, 0],
       [0, 0, 1, 0, 1, 0],
       [0, 0, 0, 1, 0, 1],
       [1, 0, 0, 0, 1, 0]])

エネルギー準位（ソート後）:
[-2. -1. -1.  1.  1.  2.]
```
この出力は、Hückel行列が正しく設定されており、固有値（エネルギー準位）がベンゼンのHückel近似の結果である \(-2\beta, -1\beta, -1\beta, 1\beta, 1\beta, 2\beta\) （ここで \(\beta=1\) と仮定）と一致していることを示しています。\(-1\) と \(1\) の準位はそれぞれ2重に縮重していることが分かります。
グラフィカルウィンドウの表示: コンソール出力の後、Matplotlibによるグラフィカルなウィンドウが表示されます。このウィンドウは、左側にベンゼンの6つの電子準位を示す水平線（\(-2\beta, -1\beta, -1\beta, 1\beta, 1\beta, 2\beta\) に対応）が描かれたエネルギー準位図、右側にそれらの準位から計算された状態密度（DOS）の曲線が表示されます。DOS曲線には、各エネルギー準位の位置に対応するピークが見られます。このウィンドウは、ユーザーが閉じるまで表示され続けます。