benzene.py ダウンロード/コピー

benzene.py をダウンロード

benzene.py
benzene.py
  1"""
  2概要:
  3    ベンゼン分子のHückel近似による電子構造計算と状態密度（DOS）の可視化を行うスクリプト。
  4詳細説明:
  5    このスクリプトは、Hückel近似法を用いてベンゼン分子（6炭素環）の電子構造を解析します。
  6    Hückel行列を定義し、その固有値問題を解くことで、分子軌道のエネルギー準位と波動関数を計算します。
  7    計算されたエネルギー準位はソートされ、その縮重度が特定されます。
  8    その後、ガウス型ブロードニングを適用して状態密度（DOS）を計算します。
  9    最後に、matplotlibライブラリを使用して、計算された電子準位と状態密度の両方を一つの図にプロットし、
 10    ベンゼン分子の電子状態を視覚的に表現します。
 11関連リンク:
 12    benzene_usage
 13"""
 14
 15import numpy as np
 16from pprint import pprint
 17import matplotlib.pyplot as plt
 18from scipy.ndimage import gaussian_filter1d
 19from collections import Counter
 20from matplotlib.gridspec import GridSpec
 21
 22#===================
 23# Global parameters
 24#===================
 25
 26# Hückel行列の定義（ベンゼンは6炭素環）
 27# α=0, β=1 とした隣接行列形式
 28H_BENZENE = np.array([[0, 1, 0, 0, 0, 1],
 29                      [1, 0, 1, 0, 0, 0],
 30                      [0, 1, 0, 1, 0, 0],
 31                      [0, 0, 1, 0, 1, 0],
 32                      [0, 0, 0, 1, 0, 1],
 33                      [1, 0, 0, 0, 1, 0]])
 34
 35# 状態密度（DOS）パラメータ
 36NUM_POINTS = 1000
 37SIGMA = 0.1  # ブロードニング幅
 38
 39def main():
 40    """
 41    概要:
 42        Hückel行列の対角化を行い、エネルギー準位と状態密度をプロットします。
 43    詳細説明:
 44        この関数は、ベンゼン分子のHückel行列を用いて固有値問題を解き、
 45        分子軌道のエネルギー準位を計算します。
 46        計算されたエネルギー準位はソートされ、その縮重度が特定されます。
 47        その後、ガウス関数を用いたブロードニングによって状態密度（DOS）を計算します。
 48        最後に、matplotlibの機能を用いて電子準位とDOSを一つの図に描画し、表示します。
 49    """
 50    print("\nHückel行列（ベンゼン）:")
 51    pprint(H_BENZENE)
 52
 53    # 固有値・固有ベクトルの計算
 54    eigenvalues, eigenvectors = np.linalg.eig(H_BENZENE)
 55    # 実数部のみ取得しソート
 56    eigenvalues = np.sort(eigenvalues.real)
 57
 58    # エネルギー準位のソートと出力
 59    print("\nエネルギー準位（ソート後）:")
 60    print(eigenvalues)
 61    print()
 62
 63    # 縮重度の計算
 64    # 浮動小数点の微小な差を丸めてカウント
 65    rounded_eigenvalues = np.round(eigenvalues, 8)
 66    degeneracies = Counter(rounded_eigenvalues)
 67
 68    # 状態密度（DOS）の計算
 69    energy_range = np.linspace(min(eigenvalues) - 1, max(eigenvalues) + 1, NUM_POINTS)
 70    dos = np.zeros(NUM_POINTS)
 71
 72    for energy, deg in degeneracies.items():
 73        # ガウス関数によるブロードニング
 74        dos += deg * np.exp(-((energy_range - energy) ** 2) / (2 * SIGMA ** 2))
 75
 76    # 状態密度のさらなる平滑化（オプション）
 77    dos = gaussian_filter1d(dos, SIGMA * (NUM_POINTS / (max(energy_range) - min(energy_range))))
 78
 79    # 図を描画
 80    fig = plt.figure(figsize=(8, 6))
 81    gs = GridSpec(1, 2, width_ratios=[1, 3], wspace=0.4)
 82
 83    # 電子準位図の描画
 84    ax0 = fig.add_subplot(gs[0])
 85    ax0.set_title('Hückel Levels')
 86    for energy, deg in degeneracies.items():
 87        # 縮重している場合は少しずらすか、線の太さを変える等の表現も可能ですが
 88        # オリジナルのプロット形式を維持します
 89        ax0.hlines(energy, 0, 1, colors='k', linestyles='solid')
 90        
 91    ax0.set_ylabel('Energy (units of β)')
 92    ax0.set_xticks([])
 93    ax0.set_xlabel('')
 94
 95    # DOSの描画
 96    ax1 = fig.add_subplot(gs[1], sharey=ax0)
 97    ax1.plot(dos, energy_range, 'k-')
 98    ax1.fill_betweenx(energy_range, 0, dos, color='gray', alpha=0.3)
 99    ax1.set_title('Density of States (DOS)')
100    ax1.set_xlabel('DOS (Arbitrary Units)')
101
102    plt.tight_layout()
103    plt.show()
104
105if __name__ == "__main__":
106    main()