第5章: 量子化学計算の実践

5.1 PySCFによる分子計算

📚 PySCFとは

PySCF (Python-based Simulations of Chemistry Framework)は、量子化学計算のためのオープンソースPythonライブラリです。

Hartree-Fock法: RHF, UHF, ROHF
DFT: 多様な交換相関汎関数（LDA, GGA, hybrid）
Post-HF法: MP2, CCSD, CI
励起状態: TD-DFT, EOM-CCSD
固体計算: PBC (Periodic Boundary Conditions)

水分子のHartree-Fock計算

PySCFで水分子（H₂O）のエネルギーと分子軌道を計算します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# PySCF互換のミニマル実装（デモ用）
class WaterMoleculeHF:
    """水分子のHartree-Fock計算（簡易版）"""
    def __init__(self, basis='sto-3g'):
        # 水分子の構造（O-H距離0.96Å, H-O-H角104.5°）
        self.atoms = [
            ('O', [0.0, 0.0, 0.0]),
            ('H', [0.0, 0.757, 0.587]),
            ('H', [0.0, -0.757, 0.587])
        ]
        self.n_electrons = 10  # 8(O) + 1(H) + 1(H)
        self.n_basis = 7      # STO-3G: O(5) + H(1) + H(1)

    def compute_integrals(self):
        """積分行列の計算（簡易版）"""
        # 重なり積分行列（対角化で単位行列に近い形）
        S = np.eye(self.n_basis)
        for i in range(self.n_basis):
            for j in range(i+1, self.n_basis):
                S[i,j] = S[j,i] = 0.2 * np.exp(-0.5 * (i-j)**2)

        # コアハミルトニアン（運動エネルギー + 核-電子引力）
        H_core = -np.diag([20.0, 15.0, 15.0, 15.0, 15.0, 7.0, 7.0])

        # 2電子積分（簡易版）
        g = np.zeros((self.n_basis, self.n_basis, self.n_basis, self.n_basis))
        for i in range(self.n_basis):
            for j in range(self.n_basis):
                g[i,i,j,j] = 10.0 / (i+j+2)
                g[i,j,i,j] = 8.0 / (i+j+3)

        return S, H_core, g

    def build_fock(self, P, H_core, g):
        """Fock行列の構築"""
        F = H_core.copy()
        for i in range(self.n_basis):
            for j in range(self.n_basis):
                for k in range(self.n_basis):
                    for l in range(self.n_basis):
                        F[i,j] += P[k,l] * (g[i,j,k,l] - 0.5 * g[i,k,j,l])
        return F

    def run_scf(self, max_iter=30, conv=1e-6):
        """SCF計算"""
        S, H_core, g = self.compute_integrals()

        # 初期密度行列
        P = np.zeros((self.n_basis, self.n_basis))
        E_old = 0.0
        energies = []

        for iteration in range(max_iter):
            # Fock行列構築
            F = self.build_fock(P, H_core, g)

            # 一般化固有値問題
            epsilon, C = eigh(F, S)

            # 密度行列更新（10電子 = 5占有軌道）
            n_occ = self.n_electrons // 2
            P_new = 2 * C[:, :n_occ] @ C[:, :n_occ].T

            # エネルギー計算
            E_elec = np.sum(P_new * (H_core + F)) / 2
            E_nuc = 9.0  # 核間反発（簡易値）
            E_total = E_elec + E_nuc
            energies.append(E_total)

            # 収束判定
            if np.abs(E_total - E_old) < conv:
                print(f"SCF収束: {iteration+1}反復")
                break

            E_old = E_total
            P = P_new

        return E_total, epsilon, C, energies

# 水分子のHF計算実行
h2o = WaterMoleculeHF()
E_total, epsilon, C, energies = h2o.run_scf()

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# (1) SCF収束
axes[0,0].plot(energies, 'o-', linewidth=2)
axes[0,0].set_xlabel('Iteration', fontsize=12)
axes[0,0].set_ylabel('Energy (Hartree)', fontsize=12)
axes[0,0].set_title('SCF Convergence', fontsize=14, fontweight='bold')
axes[0,0].grid(True, alpha=0.3)

# (2) 軌道エネルギー準位
n_occ = h2o.n_electrons // 2
colors = ['blue' if i < n_occ else 'red' for i in range(h2o.n_basis)]
axes[0,1].barh(range(h2o.n_basis), epsilon, color=colors, alpha=0.7)
axes[0,1].axvline(0, color='black', linestyle='--', alpha=0.5)
axes[0,1].set_xlabel('Energy (Hartree)', fontsize=12)
axes[0,1].set_ylabel('Orbital Index', fontsize=12)
axes[0,1].set_title('Orbital Energy Levels', fontsize=14, fontweight='bold')
axes[0,1].legend(['Occupied', 'Virtual'], loc='best')

# (3) HOMO-LUMO係数
homo_idx = n_occ - 1
lumo_idx = n_occ
x = np.arange(h2o.n_basis)
axes[1,0].bar(x - 0.2, C[:, homo_idx], width=0.4, label='HOMO', alpha=0.8)
axes[1,0].bar(x + 0.2, C[:, lumo_idx], width=0.4, label='LUMO', alpha=0.8)
axes[1,0].set_xlabel('Basis Function', fontsize=12)
axes[1,0].set_ylabel('Coefficient', fontsize=12)
axes[1,0].set_title('HOMO-LUMO Molecular Orbitals', fontsize=14, fontweight='bold')
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3, axis='y')

# (4) 計算結果サマリー
axes[1,1].axis('off')
summary = f"""
H₂O Hartree-Fock Calculation

Total Energy:  {E_total:.6f} Hartree
               {E_total * 27.211:.2f} eV

HOMO Energy:   {epsilon[homo_idx]:.6f} Hartree
LUMO Energy:   {epsilon[lumo_idx]:.6f} Hartree
HOMO-LUMO Gap: {(epsilon[lumo_idx] - epsilon[homo_idx]) * 27.211:.2f} eV

Basis Set:     STO-3G
Electrons:     {h2o.n_electrons}
Basis Funcs:   {h2o.n_basis}
"""
axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace',
                   verticalalignment='center')

plt.tight_layout()
plt.savefig('water_hf_calculation.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"H₂O Total Energy: {E_total:.6f} Hartree")
print(f"HOMO-LUMO Gap: {(epsilon[lumo_idx] - epsilon[homo_idx]) * 27.211:.2f} eV")

DFT計算との比較

同じ水分子をDFT（B3LYP汎関数）で計算し、HF法と比較します。

class WaterMoleculeDFT(WaterMoleculeHF):
    """水分子のDFT計算（B3LYP簡易版）"""
    def __init__(self, functional='b3lyp'):
        super().__init__()
        self.functional = functional

    def xc_energy(self, rho):
        """交換相関エネルギー（B3LYP近似）"""
        # LDA exchange
        C_x = (3/4) * (3/np.pi)**(1/3)
        E_x_lda = -C_x * rho**(4/3)

        # GGA correction (gradient approximation)
        E_x_gga = E_x_lda * (1.0 + 0.1 * rho**(1/3))

        # B3LYP mixing (20% HF + 80% DFT)
        E_xc = 0.8 * E_x_gga

        return E_xc

    def build_fock_dft(self, P, H_core, g):
        """Kohn-Sham Fock行列"""
        F = H_core.copy()

        # Coulomb + Exchange
        for i in range(self.n_basis):
            for j in range(self.n_basis):
                for k in range(self.n_basis):
                    for l in range(self.n_basis):
                        # Coulomb
                        F[i,j] += P[k,l] * g[i,j,k,l]
                        # Exchange (20% HF for B3LYP)
                        F[i,j] -= 0.2 * P[k,l] * g[i,k,j,l]

        # XC potential (simplified)
        rho = np.diag(P)
        V_xc = self.xc_energy(rho + 1e-10)
        F += np.diag(V_xc)

        return F

    def run_dft(self, max_iter=30, conv=1e-6):
        """DFT-SCF計算"""
        S, H_core, g = self.compute_integrals()
        P = np.zeros((self.n_basis, self.n_basis))
        E_old = 0.0
        energies = []

        for iteration in range(max_iter):
            F = self.build_fock_dft(P, H_core, g)
            epsilon, C = eigh(F, S)

            n_occ = self.n_electrons // 2
            P_new = 2 * C[:, :n_occ] @ C[:, :n_occ].T

            E_elec = np.sum(P_new * (H_core + F)) / 2
            E_nuc = 9.0
            E_total = E_elec + E_nuc
            energies.append(E_total)

            if np.abs(E_total - E_old) < conv:
                print(f"DFT-SCF収束: {iteration+1}反復")
                break

            E_old = E_total
            P = P_new

        return E_total, epsilon, C, energies

# DFT計算実行
h2o_dft = WaterMoleculeDFT()
E_dft, eps_dft, C_dft, energies_dft = h2o_dft.run_dft()

# HF vs DFT 比較
print(f"\n=== H₂O Energy Comparison ===")
print(f"HF  Energy: {E_total:.6f} Hartree")
print(f"DFT Energy: {E_dft:.6f} Hartree")
print(f"Difference: {(E_dft - E_total)*1000:.2f} mHartree")

n_occ = h2o.n_electrons // 2
homo_hf = epsilon[n_occ-1]
lumo_hf = epsilon[n_occ]
homo_dft = eps_dft[n_occ-1]
lumo_dft = eps_dft[n_occ]

print(f"\nHOMO-LUMO Gap:")
print(f"HF:  {(lumo_hf - homo_hf)*27.211:.2f} eV")
print(f"DFT: {(lumo_dft - homo_dft)*27.211:.2f} eV")

5.2 構造最適化と振動解析

📚 ポテンシャルエネルギー曲面

構造最適化は、ポテンシャルエネルギー曲面（PES）上で最小エネルギー構造を探索する過程です。

エネルギー勾配:

\[ \mathbf{g}_I = -\frac{\partial E}{\partial \mathbf{R}_I} \]

力が0になる点（\(\mathbf{g}_I = 0\)）が安定構造です。

水素分子の結合距離最適化

H₂分子の結合距離を最適化し、平衡構造とエネルギーを求めます。

from scipy.optimize import minimize

def h2_energy(R):
    """H₂分子のエネルギー（Morse potential近似）"""
    D_e = 4.75  # eV
    alpha = 1.44  # Å⁻¹
    R_e = 0.741  # Å

    E = D_e * ((1 - np.exp(-alpha * (R - R_e)))**2 - 1)
    return E

def h2_gradient(R):
    """エネルギー勾配"""
    D_e = 4.75
    alpha = 1.44
    R_e = 0.741

    exp_term = np.exp(-alpha * (R - R_e))
    dE_dR = 2 * D_e * alpha * (1 - exp_term) * exp_term
    return dE_dR

# 構造最適化（初期値: R=1.0Å）
result = minimize(h2_energy, x0=[1.0], method='BFGS', jac=h2_gradient)

R_opt = result.x[0]
E_opt = result.fun

print(f"最適化結果:")
print(f"  平衡距離: {R_opt:.4f} Å")
print(f"  結合エネルギー: {E_opt:.4f} eV")

# ポテンシャル曲線とPES
R_range = np.linspace(0.3, 3.0, 100)
E_range = [h2_energy(r) for r in R_range]

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# (1) ポテンシャルエネルギー曲線
axes[0,0].plot(R_range, E_range, 'b-', linewidth=2, label='Morse Potential')
axes[0,0].plot(R_opt, E_opt, 'ro', markersize=10, label=f'Minimum: {R_opt:.3f}Å')
axes[0,0].axhline(0, color='k', linestyle='--', alpha=0.3)
axes[0,0].set_xlabel('H-H Distance (Å)', fontsize=12)
axes[0,0].set_ylabel('Energy (eV)', fontsize=12)
axes[0,0].set_title('H₂ Potential Energy Curve', fontsize=14, fontweight='bold')
axes[0,0].legend()
axes[0,0].grid(True, alpha=0.3)

# (2) エネルギー勾配
gradient_range = [h2_gradient(r) for r in R_range]
axes[0,1].plot(R_range, gradient_range, 'g-', linewidth=2)
axes[0,1].axhline(0, color='r', linestyle='--', linewidth=2, label='Zero gradient')
axes[0,1].plot(R_opt, 0, 'ro', markersize=10)
axes[0,1].set_xlabel('H-H Distance (Å)', fontsize=12)
axes[0,1].set_ylabel('dE/dR (eV/Å)', fontsize=12)
axes[0,1].set_title('Energy Gradient', fontsize=14, fontweight='bold')
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)

# (3) 振動エネルギー準位
# Morse振動子の準位
omega = 4401.2  # cm⁻¹ (H₂振動数)
D_e_cm = D_e * 8065.5  # eV → cm⁻¹
n_max = int(D_e_cm / omega)

n_levels = np.arange(0, min(n_max, 15))
E_vib = omega * (n_levels + 0.5) - omega**2 * (n_levels + 0.5)**2 / (4 * D_e_cm)

axes[1,0].barh(n_levels, E_vib / 8065.5, height=0.8, alpha=0.7)
axes[1,0].set_ylabel('Vibrational Quantum Number (n)', fontsize=12)
axes[1,0].set_xlabel('Energy (eV)', fontsize=12)
axes[1,0].set_title('Vibrational Energy Levels', fontsize=14, fontweight='bold')
axes[1,0].grid(True, alpha=0.3, axis='x')

# (4) 最適化サマリー
axes[1,1].axis('off')
summary = f"""
H₂ Geometry Optimization

Equilibrium Distance:  {R_opt:.4f} Å
Bond Energy:           {-E_opt:.4f} eV
                       {-E_opt * 96.485:.1f} kJ/mol

Vibrational Frequency: {omega:.1f} cm⁻¹
Zero-Point Energy:     {E_vib[0] / 8065.5:.4f} eV

Optimization Method:   BFGS
Convergence:           {result.success}
Iterations:            {result.nit}
"""
axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace',
                   verticalalignment='center')

plt.tight_layout()
plt.savefig('h2_geometry_optimization.png', dpi=150, bbox_inches='tight')
plt.show()

5.3 励起状態計算（TD-DFT）

📚 時間依存密度汎関数理論

TD-DFT (Time-Dependent DFT)は、励起状態の電子構造を計算する手法です。

線形応答理論に基づき、基底状態からの励起エネルギーを求めます:

\[ \omega_n = E_n - E_0 \]

吸収スペクトル、発光スペクトル、光化学反応の解析に使われます。

エチレン分子の励起状態

エチレン（C₂H₄）のπ→π*遷移を計算し、紫外吸収スペクトルを予測します。

class EthyleneTDDFT:
    """エチレン分子のTD-DFT計算（簡易版）"""
    def __init__(self):
        self.n_occ = 6  # 占有軌道数
        self.n_virt = 4  # 仮想軌道数

        # 軌道エネルギー（実験値近似、Hartree単位）
        self.epsilon_occ = np.array([-0.8, -0.7, -0.6, -0.5, -0.4, -0.35])
        self.epsilon_virt = np.array([0.1, 0.2, 0.3, 0.4])

    def compute_excitations(self):
        """一電子励起エネルギー"""
        excitations = []

        for i, e_occ in enumerate(self.epsilon_occ):
            for a, e_virt in enumerate(self.epsilon_virt):
                omega = e_virt - e_occ
                # 振動子強度（簡易近似）
                f = 0.1 * np.exp(-2 * (i + a))
                excitations.append({
                    'from': i,
                    'to': a,
                    'energy_hartree': omega,
                    'energy_eV': omega * 27.211,
                    'wavelength_nm': 1239.8 / (omega * 27.211),
                    'oscillator_strength': f
                })

        return sorted(excitations, key=lambda x: x['energy_eV'])

    def absorption_spectrum(self, excitations, lambda_range):
        """吸収スペクトルの計算"""
        spectrum = np.zeros_like(lambda_range)

        for exc in excitations:
            # Lorentzian broadening
            gamma = 20  # nm (半値幅)
            lambda_0 = exc['wavelength_nm']
            f = exc['oscillator_strength']

            spectrum += f * (gamma / 2) / ((lambda_range - lambda_0)**2 + (gamma / 2)**2)

        return spectrum

# TD-DFT計算実行
ethylene = EthyleneTDDFT()
excitations = ethylene.compute_excitations()

# 主要な励起状態（振動子強度 > 0.01）
major_exc = [exc for exc in excitations if exc['oscillator_strength'] > 0.01]

print("主要な励起状態:")
for i, exc in enumerate(major_exc[:5], 1):
    print(f"{i}. {exc['energy_eV']:.2f} eV ({exc['wavelength_nm']:.1f} nm) "
          f"[{exc['from']}→{exc['to']}] f={exc['oscillator_strength']:.3f}")

# スペクトル計算
lambda_range = np.linspace(100, 400, 500)
spectrum = ethylene.absorption_spectrum(excitations, lambda_range)

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# (1) 軌道エネルギーダイアグラム
y_occ = np.arange(ethylene.n_occ)
y_virt = np.arange(ethylene.n_virt) + ethylene.n_occ + 1

axes[0,0].barh(y_occ, ethylene.epsilon_occ * 27.211,
                height=0.6, color='blue', alpha=0.7, label='Occupied')
axes[0,0].barh(y_virt, ethylene.epsilon_virt * 27.211,
                height=0.6, color='red', alpha=0.7, label='Virtual')
axes[0,0].axvline(0, color='k', linestyle='--', alpha=0.3)

# π→π*遷移の矢印
homo = ethylene.n_occ - 1
lumo = ethylene.n_occ + 1
axes[0,0].annotate('', xy=(ethylene.epsilon_virt[0]*27.211, lumo),
                    xytext=(ethylene.epsilon_occ[-1]*27.211, homo),
                    arrowprops=dict(arrowstyle='->', lw=2, color='green'))
axes[0,0].text(-5, (homo+lumo)/2, 'π→π*', fontsize=12, color='green',
                fontweight='bold')

axes[0,0].set_xlabel('Energy (eV)', fontsize=12)
axes[0,0].set_ylabel('Orbital Index', fontsize=12)
axes[0,0].set_title('Molecular Orbital Diagram', fontsize=14, fontweight='bold')
axes[0,0].legend()

# (2) 励起エネルギースペクトル（スティックスペクトル）
energies = [exc['energy_eV'] for exc in major_exc]
intensities = [exc['oscillator_strength'] for exc in major_exc]

axes[0,1].stem(energies, intensities, basefmt=' ')
axes[0,1].set_xlabel('Excitation Energy (eV)', fontsize=12)
axes[0,1].set_ylabel('Oscillator Strength', fontsize=12)
axes[0,1].set_title('Excitation Spectrum', fontsize=14, fontweight='bold')
axes[0,1].grid(True, alpha=0.3)

# (3) UV-Vis吸収スペクトル
axes[1,0].plot(lambda_range, spectrum, 'b-', linewidth=2)
axes[1,0].fill_between(lambda_range, spectrum, alpha=0.3)
axes[1,0].set_xlabel('Wavelength (nm)', fontsize=12)
axes[1,0].set_ylabel('Absorbance (a.u.)', fontsize=12)
axes[1,0].set_title('UV-Vis Absorption Spectrum', fontsize=14, fontweight='bold')
axes[1,0].set_xlim(100, 400)
axes[1,0].grid(True, alpha=0.3)

# 主要ピーク位置を注釈
max_idx = np.argmax(spectrum)
max_lambda = lambda_range[max_idx]
axes[1,0].annotate(f'{max_lambda:.1f} nm',
                    xy=(max_lambda, spectrum[max_idx]),
                    xytext=(max_lambda+50, spectrum[max_idx]*0.8),
                    arrowprops=dict(arrowstyle='->', color='red'),
                    fontsize=11, color='red', fontweight='bold')

# (4) 励起状態サマリー
axes[1,1].axis('off')
summary = f"""
C₂H₄ TD-DFT Excitation Analysis

Lowest Excitation:
  Energy: {major_exc[0]['energy_eV']:.2f} eV
  Wavelength: {major_exc[0]['wavelength_nm']:.1f} nm
  Type: π → π*
  f: {major_exc[0]['oscillator_strength']:.3f}

Max Absorption:
  λmax: {max_lambda:.1f} nm
  Region: {'UV-C' if max_lambda < 280 else 'UV-B' if max_lambda < 315 else 'UV-A'}

Total Excitations: {len(excitations)}
Major (f>0.01): {len(major_exc)}
"""
axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace',
                   verticalalignment='center')

plt.tight_layout()
plt.savefig('ethylene_tddft.png', dpi=150, bbox_inches='tight')
plt.show()

5.4 固体バンド計算

📚 周期境界条件とBloch定理

固体は周期的な結晶構造を持ちます。Bloch定理により、波動関数は:

\[ \psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}} u_{n\mathbf{k}}(\mathbf{r}) \]

ここで\(u_{n\mathbf{k}}\)は格子周期関数です。\(\mathbf{k}\)は波数ベクトルで、バンド構造\(E_n(\mathbf{k})\)を定義します。

1次元周期ポテンシャルのバンド構造

Kronig-Penneyモデルで1次元固体のバンドギャップを計算します。

class KronigPenney1D:
    """1次元周期ポテンシャル（Kronig-Penneyモデル）"""
    def __init__(self, a=1.0, V0=10.0, N_pw=21):
        self.a = a        # 格子定数
        self.V0 = V0      # ポテンシャル深さ
        self.N_pw = N_pw  # 平面波数

    def potential(self, x):
        """周期ポテンシャル（コサイン型）"""
        return -self.V0 * np.cos(2 * np.pi * x / self.a)

    def hamiltonian_matrix(self, k):
        """k点でのハミルトニアン行列"""
        # 平面波基底: exp(i(k + n*G)x), G = 2π/a
        G = 2 * np.pi / self.a
        n_indices = np.arange(-self.N_pw // 2, self.N_pw // 2 + 1)

        H = np.zeros((self.N_pw, self.N_pw), dtype=complex)

        for i, n_i in enumerate(n_indices):
            for j, n_j in enumerate(n_indices):
                k_i = k + n_i * G
                k_j = k + n_j * G

                if i == j:
                    # 運動エネルギー
                    H[i,j] = 0.5 * k_i**2
                elif abs(n_i - n_j) == 1:
                    # ポテンシャル（第一近接のみ）
                    H[i,j] = -self.V0 / 2

        return H

    def compute_bands(self, k_points):
        """バンド構造計算"""
        bands = []

        for k in k_points:
            H = self.hamiltonian_matrix(k)
            eigenvalues = eigh(H, eigvals_only=True)
            bands.append(eigenvalues)

        return np.array(bands)

# バンド構造計算
kp = KronigPenney1D(a=1.0, V0=5.0, N_pw=21)

# 第1ブリルアンゾーン: -π/a ≤ k ≤ π/a
k_points = np.linspace(-np.pi/kp.a, np.pi/kp.a, 100)
bands = kp.compute_bands(k_points)

# バンドギャップ検出
def find_band_gap(bands):
    """直接バンドギャップを検出"""
    valence_band_max = np.max(bands[:, :11])  # 下位バンド
    conduction_band_min = np.min(bands[:, 11:])  # 上位バンド
    gap = conduction_band_min - valence_band_max
    return gap, valence_band_max, conduction_band_min

gap, vbm, cbm = find_band_gap(bands)

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# (1) 周期ポテンシャル
x = np.linspace(0, 3*kp.a, 500)
V = kp.potential(x)
axes[0,0].plot(x, V, 'b-', linewidth=2)
for i in range(4):
    axes[0,0].axvline(i*kp.a, color='gray', linestyle='--', alpha=0.5)
axes[0,0].set_xlabel('Position x/a', fontsize=12)
axes[0,0].set_ylabel('Potential V(x)', fontsize=12)
axes[0,0].set_title('Periodic Potential', fontsize=14, fontweight='bold')
axes[0,0].grid(True, alpha=0.3)

# (2) バンド構造
for n in range(min(10, bands.shape[1])):
    color = 'blue' if n < 5 else 'red'
    axes[0,1].plot(k_points * kp.a / np.pi, bands[:, n], color=color, linewidth=2)

axes[0,1].axhline(vbm, color='blue', linestyle='--', alpha=0.5, label='VBM')
axes[0,1].axhline(cbm, color='red', linestyle='--', alpha=0.5, label='CBM')
axes[0,1].axvline(0, color='black', linestyle='-', alpha=0.3)
axes[0,1].set_xlabel('Wave Vector (k/π/a)', fontsize=12)
axes[0,1].set_ylabel('Energy (a.u.)', fontsize=12)
axes[0,1].set_title('Band Structure', fontsize=14, fontweight='bold')
axes[0,1].set_xlim(-1, 1)
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)

# バンドギャップ領域をハイライト
axes[0,1].fill_between([-1, 1], vbm, cbm, alpha=0.2, color='yellow')

# (3) 状態密度（DOS）
energy_bins = np.linspace(bands.min(), bands.max(), 200)
dos, _ = np.histogram(bands.flatten(), bins=energy_bins, density=True)

axes[1,0].plot(dos, energy_bins[:-1], 'b-', linewidth=2)
axes[1,0].fill_betweenx(energy_bins[:-1], dos, alpha=0.3)
axes[1,0].axhline(vbm, color='blue', linestyle='--')
axes[1,0].axhline(cbm, color='red', linestyle='--')
axes[1,0].set_xlabel('Density of States', fontsize=12)
axes[1,0].set_ylabel('Energy (a.u.)', fontsize=12)
axes[1,0].set_title('Density of States (DOS)', fontsize=14, fontweight='bold')
axes[1,0].grid(True, alpha=0.3)

# (4) 計算結果サマリー
axes[1,1].axis('off')
summary = f"""
1D Periodic Crystal Band Structure

Lattice Constant: {kp.a:.2f}
Potential Depth:  {kp.V0:.2f}
Plane Waves:      {kp.N_pw}

Band Gap:         {gap:.4f} a.u.
                  {gap * 27.211:.2f} eV

VBM:              {vbm:.4f} a.u.
CBM:              {cbm:.4f} a.u.

Material Type:    {'Insulator' if gap > 0.2 else 'Semiconductor'}
"""
axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace',
                   verticalalignment='center')

plt.tight_layout()
plt.savefig('band_structure_1d.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"\nバンドギャップ: {gap * 27.211:.2f} eV")
print(f"材料タイプ: {'絶縁体' if gap > 0.2 else '半導体'}")

🎯 演習問題

メタン（CH₄）のHF計算: 正四面体構造のメタン分子でHF-SCF計算を実装し、全エネルギーと分子軌道を求めよ。
CO₂の構造最適化: CO₂分子のC=O結合距離と∠OCO角度を最適化し、赤外活性な振動モードを計算せよ。
ベンゼンの励起状態: ベンゼンのπ→π*遷移をTD-DFTで計算し、紫外吸収スペクトルを予測せよ。
シリコンのバンド構造: ダイヤモンド構造シリコンのバンドギャップを計算し、実験値（1.1 eV）と比較せよ。
反応経路探索: H₂ + F → HF + H反応の遷移状態を探索し、活性化エネルギーを求めよ。

5.5 材料科学への応用

📚 第一原理計算の材料応用

量子化学計算は、実験に先立つ材料スクリーニング、反応機構解明、物性予測に不可欠です。

触媒設計: 吸着エネルギー、反応障壁の計算
電池材料: 電極電位、イオン拡散係数
半導体: バンドギャップ、キャリア移動度
光学材料: 吸収スペクトル、非線形光学応答
機械学習ポテンシャル: DFTデータで訓練、高速MD

触媒表面吸着エネルギー

金属表面へのCO吸着を計算し、触媒活性を評価します。

class MetalSurfaceAdsorption:
    """金属表面吸着エネルギー計算（簡易モデル）"""
    def __init__(self, metal='Pt'):
        self.metal = metal
        # 金属の仕事関数（eV）
        self.work_functions = {'Pt': 5.65, 'Pd': 5.12, 'Cu': 4.65, 'Au': 5.1}
        self.phi = self.work_functions.get(metal, 5.0)

    def lennard_jones(self, r, epsilon, sigma):
        """Lennard-Jonesポテンシャル"""
        return 4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)

    def adsorption_energy(self, z, molecule='CO'):
        """吸着エネルギーの計算"""
        # LJ パラメータ（金属-CO）
        epsilon = 0.5  # eV
        sigma = 2.5    # Å

        # 吸着エネルギー = LJ + 静電相互作用
        E_vdw = self.lennard_jones(z, epsilon, sigma)
        E_elec = -0.3 * self.phi * np.exp(-z / 2.0)

        return E_vdw + E_elec

    def find_optimal_height(self):
        """最適吸着高さの探索"""
        z_range = np.linspace(1.5, 6.0, 200)
        E_ads = [self.adsorption_energy(z) for z in z_range]

        min_idx = np.argmin(E_ads)
        z_opt = z_range[min_idx]
        E_opt = E_ads[min_idx]

        return z_opt, E_opt, z_range, E_ads

# 複数金属の比較
metals = ['Pt', 'Pd', 'Cu', 'Au']
results = {}

for metal in metals:
    surface = MetalSurfaceAdsorption(metal)
    z_opt, E_opt, z_range, E_ads = surface.find_optimal_height()
    results[metal] = {'z': z_opt, 'E': E_opt, 'curve': (z_range, E_ads)}
    print(f"{metal}: z = {z_opt:.2f} Å, E_ads = {E_opt:.3f} eV")

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# (1) 吸着エネルギー曲線
colors = {'Pt': 'blue', 'Pd': 'green', 'Cu': 'orange', 'Au': 'gold'}
for metal, data in results.items():
    z_range, E_ads = data['curve']
    axes[0,0].plot(z_range, E_ads, color=colors[metal], linewidth=2, label=metal)
    axes[0,0].plot(data['z'], data['E'], 'o', color=colors[metal], markersize=8)

axes[0,0].axhline(0, color='k', linestyle='--', alpha=0.3)
axes[0,0].set_xlabel('Distance from Surface (Å)', fontsize=12)
axes[0,0].set_ylabel('Adsorption Energy (eV)', fontsize=12)
axes[0,0].set_title('CO Adsorption on Metal Surfaces', fontsize=14, fontweight='bold')
axes[0,0].legend()
axes[0,0].grid(True, alpha=0.3)

# (2) 吸着エネルギーの比較
metal_names = list(results.keys())
E_ads_values = [-results[m]['E'] for m in metal_names]  # 正の値で表示
bar_colors = [colors[m] for m in metal_names]

axes[0,1].bar(metal_names, E_ads_values, color=bar_colors, alpha=0.7)
axes[0,1].set_ylabel('|E_ads| (eV)', fontsize=12)
axes[0,1].set_title('Adsorption Energy Comparison', fontsize=14, fontweight='bold')
axes[0,1].grid(True, alpha=0.3, axis='y')

# 最強吸着をハイライト
max_idx = np.argmax(E_ads_values)
axes[0,1].text(max_idx, E_ads_values[max_idx] + 0.05, 'Strongest',
                ha='center', fontsize=10, fontweight='bold')

# (3) 仕事関数 vs 吸着エネルギー
phi_values = [MetalSurfaceAdsorption(m).phi for m in metal_names]
axes[1,0].scatter(phi_values, E_ads_values, s=150, c=bar_colors, alpha=0.7)

for i, metal in enumerate(metal_names):
    axes[1,0].annotate(metal, (phi_values[i], E_ads_values[i]),
                        xytext=(5, 5), textcoords='offset points',
                        fontsize=11, fontweight='bold')

# 線形フィット
z = np.polyfit(phi_values, E_ads_values, 1)
p = np.poly1d(z)
phi_fit = np.linspace(min(phi_values), max(phi_values), 100)
axes[1,0].plot(phi_fit, p(phi_fit), 'r--', alpha=0.5, label=f'y = {z[0]:.2f}x + {z[1]:.2f}')

axes[1,0].set_xlabel('Work Function (eV)', fontsize=12)
axes[1,0].set_ylabel('|E_ads| (eV)', fontsize=12)
axes[1,0].set_title('Work Function vs Adsorption Energy', fontsize=14, fontweight='bold')
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3)

# (4) 触媒活性サマリー
axes[1,1].axis('off')
best_metal = metal_names[max_idx]
summary = f"""
CO Adsorption on Metal Surfaces

Best Catalyst: {best_metal}
  E_ads:       {-results[best_metal]['E']:.3f} eV
  Height:      {results[best_metal]['z']:.2f} Å
  Work Func:   {MetalSurfaceAdsorption(best_metal).phi:.2f} eV

Ranking (by |E_ads|):
"""
for i, m in enumerate(sorted(metal_names, key=lambda x: -results[x]['E']), 1):
    summary += f"  {i}. {m:3s}  {-results[m]['E']:.3f} eV\n"

axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace',
                   verticalalignment='center')

plt.tight_layout()
plt.savefig('catalyst_adsorption.png', dpi=150, bbox_inches='tight')
plt.show()

まとめ

本章では、量子化学計算の実践的手法を学びました:

PySCFによる分子計算: HF法とDFTで水分子の電子状態を計算
構造最適化: エネルギー勾配に基づく平衡構造の探索
振動解析: 調和振動子近似による振動スペクトル予測
TD-DFT: 励起状態計算とUV-Vis吸収スペクトル
バンド計算: 周期系のBloch定理と固体電子構造
材料応用: 触媒、電池、半導体への実用的応用

これらの手法は、Materials Informaticsにおいて、機械学習と組み合わせることで、高速な材料スクリーニングと物性予測を可能にします。第一原理計算で生成された高精度データは、機械学習ポテンシャルの訓練、記述子の設計、予測モデルの構築に活用されます。