第4章：フォノン計算と熱力学特性

格子振動と状態密度のイメージを掴み、熱伝導など物性計算の入り口に立ちます。

💡 補足: フォノンは“原子の合奏”。周波数分布から熱や振動の性格が読み解けます。

学習目標

この章を読むことで、以下を習得できます： - 格子振動（フォノン）の理論的基礎を理解する - 調和近似と動力学行列の概念を説明できる - Phonopyで実際にフォノン計算を実行できる - フォノンバンド構造と状態密度を解釈できる - 熱力学特性（比熱、自由エネルギー、エントロピー）を計算できる

4.1 格子振動の理論

調和振動子モデル

結晶中の原子は平衡位置$\mathbf{R}_I^0$の周りで振動します。変位を$\mathbf{u}_I$とすると：

$$ \mathbf{R}_I = \mathbf{R}_I^0 + \mathbf{u}_I $$

ポテンシャルエネルギーをTaylor展開（調和近似）：

$$ U = U_0 + \sum_I \frac{\partial U}{\partial \mathbf{u}_I}\Bigg|_0 \mathbf{u}_I + \frac{1}{2}\sum_{I,J} \mathbf{u}_I \cdot \frac{\partial^2 U}{\partial \mathbf{u}_I \partial \mathbf{u}_J}\Bigg|_0 \cdot \mathbf{u}_J + O(\mathbf{u}^3) $$

平衡位置では第1項がゼロなので：

$$ U \approx U_0 + \frac{1}{2}\sum_{I,J} \mathbf{u}_I \cdot \mathbf{\Phi}_{IJ} \cdot \mathbf{u}_J $$

$\mathbf{\Phi}_{IJ}$は力定数行列（force constant matrix）：

$$ \Phi_{IJ}^{\alpha\beta} = \frac{\partial^2 U}{\partial u_I^\alpha \partial u_J^\beta}\Bigg|_0 $$

運動方程式

原子$I$の運動方程式：

$$ M_I \frac{d^2 u_I^\alpha}{dt^2} = -\sum_{J,\beta} \Phi_{IJ}^{\alpha\beta} u_J^\beta $$

平面波解を仮定：

$$ u_I^\alpha = \frac{1}{\sqrt{M_I}} e_I^\alpha(\mathbf{k}) e^{i(\mathbf{k}\cdot\mathbf{R}_I - \omega t)} $$

動力学行列（dynamical matrix）を導入：

$$ D_{IJ}^{\alpha\beta}(\mathbf{k}) = \frac{1}{\sqrt{M_I M_J}} \Phi_{IJ}^{\alpha\beta} e^{i\mathbf{k}\cdot(\mathbf{R}_I - \mathbf{R}_J)} $$

固有値問題

$$ \sum_{J,\beta} D_{IJ}^{\alpha\beta}(\mathbf{k}) e_J^\beta(\mathbf{k}) = \omega^2(\mathbf{k}) e_I^\alpha(\mathbf{k}) $$

• $\omega(\mathbf{k})$: フォノンの角振動数（フォノン分散）
• $\mathbf{e}(\mathbf{k})$: フォノンの固有ベクトル（偏極ベクトル）

フォノンバンドの数: $3N_{\text{atom}}$本（単位胞内原子数$N_{\text{atom}}$に対して）

4.2 フォノンの分類

音響モード（Acoustic modes）

• $\omega(\mathbf{k}) \to 0$ as $\mathbf{k} \to 0$
• Γ点（$\mathbf{k}=0$）で周波数がゼロ
• 3本（1本の縦波LA、2本の横波TA）
• 物理的意味: 並進運動

光学モード（Optical modes）

• $\omega(\mathbf{k}) \neq 0$ at $\mathbf{k} = 0$
• Γ点で有限の周波数
• $3(N_{\text{atom}} - 1)$本
• 物理的意味: 単位胞内の相対運動

縦波（Longitudinal）vs 横波（Transverse）

• 縦波（L）: 振動方向が波の進行方向と平行
• 横波（T）: 振動方向が波の進行方向と垂直

単原子結晶（Si, Cuなど）: - 音響モード: 3本（1 LA + 2 TA） - 光学モード: なし

2原子結晶（NaCl, GaAsなど）: - 音響モード: 3本 - 光学モード: 3本（1 LO + 2 TO）

4.3 Phonopyによる実践

インストール

pip install phonopy
# またはcondaで
conda install -c conda-forge phonopy

Example 1: Siのフォノン計算（GPAW使用）

Step 1: 基底状態計算

from ase.build import bulk
from gpaw import GPAW, PW

# Si結晶の作成
si = bulk('Si', 'diamond', a=5.43)

# DFT計算
calc = GPAW(mode=PW(400),
            xc='PBE',
            kpts=(8, 8, 8),
            txt='si_gs.txt')

si.calc = calc
si.get_potential_energy()
calc.write('si_groundstate.gpw')

Step 2: Phonopy用のスーパーセル作成

# phonopy_disp.confを作成
cat > phonopy_disp.conf <<EOF
DIM = 2 2 2
ATOM_NAME = Si
EOF

# Phonopyで変位構造を生成
phonopy -d --dim="2 2 2" --gpaw

これによりsupercell-XXX.pyファイルが生成されます。

Step 3: 力の計算

# supercell-001.pyを実行（各変位構造で）
from gpaw import GPAW

calc = GPAW('si_groundstate.gpw', txt=None)
atoms = calc.get_atoms()
forces = atoms.get_forces()

# FORCE_SETSファイルに書き込み（Phonopyが読み込む形式）
import numpy as np
np.savetxt('forces_001.dat', forces)

全変位構造に対して繰り返す（通常はスクリプトで自動化）

Step 4: フォノン計算

from phonopy import Phonopy
from phonopy.interface.gpaw import read_gpaw
import matplotlib.pyplot as plt

# Phonopyオブジェクトの作成
unitcell, calc_forces = read_gpaw('si_groundstate.gpw')
phonon = Phonopy(unitcell, [[2, 0, 0], [0, 2, 0], [0, 0, 2]])

# 力定数行列の設定（FORCE_SETSから）
phonon.set_displacement_dataset(dataset)
phonon.produce_force_constants()

# バンド構造の計算
path = [[[0, 0, 0], [0.5, 0, 0.5], [0.625, 0.25, 0.625],
         [0.375, 0.375, 0.75], [0, 0, 0], [0.5, 0.5, 0.5]]]
labels = ["$\\Gamma$", "X", "U", "K", "$\\Gamma$", "L"]
qpoints, connections = phonon.get_band_structure_plot_data(path)

phonon.plot_band_structure(path, labels=labels)
plt.ylabel('Frequency (THz)')
plt.savefig('si_phonon_band.png', dpi=150)
plt.show()

# 状態密度（DOS）
phonon.set_mesh([20, 20, 20])
phonon.set_total_DOS()
dos_freq, dos_val = phonon.get_total_DOS()

plt.figure(figsize=(8, 6))
plt.plot(dos_freq, dos_val, linewidth=2)
plt.xlabel('Frequency (THz)', fontsize=12)
plt.ylabel('DOS (states/THz)', fontsize=12)
plt.title('Si Phonon DOS', fontsize=14)
plt.grid(alpha=0.3)
plt.savefig('si_phonon_dos.png', dpi=150)
plt.show()

Example 2: 完全な自動化スクリプト

from ase.build import bulk
from gpaw import GPAW, PW
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
import numpy as np

def calculate_phonon(symbol='Si', a=5.43, dim=(2,2,2)):
    """
    完全自動フォノン計算
    """
    # 1. 基底状態計算
    print("Step 1: Ground state calculation...")
    atoms = bulk(symbol, 'diamond', a=a)
    calc = GPAW(mode=PW(400), xc='PBE', kpts=(8,8,8), txt=None)
    atoms.calc = calc
    atoms.get_potential_energy()

    # 2. Phonopyセットアップ
    print("Step 2: Generate displacements...")
    cell = PhonopyAtoms(symbols=[symbol]*len(atoms),
                       cell=atoms.cell,
                       scaled_positions=atoms.get_scaled_positions())

    phonon = Phonopy(cell, np.diag(dim))
    phonon.generate_displacements(distance=0.01)

    # 3. 力の計算
    print(f"Step 3: Calculate forces for {len(phonon.supercells_with_displacements)} supercells...")
    set_of_forces = []

    for i, scell in enumerate(phonon.supercells_with_displacements):
        supercell = bulk(symbol, 'diamond', a=a).repeat(dim)
        supercell.positions = scell.positions
        supercell.calc = calc
        forces = supercell.get_forces()
        drift_force = forces.sum(axis=0)
        for force in forces:
            force -= drift_force / len(forces)  # ドリフト補正
        set_of_forces.append(forces)

    # 4. フォノン計算
    print("Step 4: Phonon calculation...")
    phonon.produce_force_constants(forces=set_of_forces)

    # バンド構造
    path = [[[0, 0, 0], [0.5, 0, 0.5], [0.625, 0.25, 0.625],
             [0.375, 0.375, 0.75], [0, 0, 0], [0.5, 0.5, 0.5]]]
    labels = ["$\\Gamma$", "X", "U", "K", "$\\Gamma$", "L"]

    phonon.auto_band_structure(plot=True, labels=labels, filename=f'{symbol}_band.png')

    # DOS
    phonon.auto_total_dos(plot=True, filename=f'{symbol}_dos.png')

    print("Done!")
    return phonon

# 実行
si_phonon = calculate_phonon('Si', a=5.43, dim=(2,2,2))

4.4 熱力学特性の計算

自由エネルギー

Helmholtz自由エネルギー（NVT）:

$$ F(T) = U_0 + k_B T \sum_{\mathbf{q},j} \ln\left[2\sinh\left(\frac{\hbar\omega_{\mathbf{q}j}}{2k_B T}\right)\right] $$

• $U_0$: ゼロ点エネルギー
• $\omega_{\mathbf{q}j}$: 波数$\mathbf{q}$、バンド$j$のフォノン周波数

内部エネルギー

$$ U(T) = U_0 + \sum_{\mathbf{q},j} \hbar\omega_{\mathbf{q}j} \left[\frac{1}{2} + n_B(\omega_{\mathbf{q}j}, T)\right] $$

$n_B(\omega, T)$はBose-Einstein分布関数：

$$ n_B(\omega, T) = \frac{1}{e^{\hbar\omega/(k_B T)} - 1} $$

エントロピー

$$ S(T) = -\left(\frac{\partial F}{\partial T}\right)_V = k_B \sum_{\mathbf{q},j} \left[\frac{\hbar\omega_{\mathbf{q}j}}{k_B T} n_B(\omega_{\mathbf{q}j}, T) - \ln(1 - e^{-\hbar\omega_{\mathbf{q}j}/(k_B T)})\right] $$

比熱（定積）

$$ C_V(T) = \left(\frac{\partial U}{\partial T}\right)_V = k_B \sum_{\mathbf{q},j} \left(\frac{\hbar\omega_{\mathbf{q}j}}{k_B T}\right)^2 \frac{e^{\hbar\omega_{\mathbf{q}j}/(k_B T)}}{(e^{\hbar\omega_{\mathbf{q}j}/(k_B T)} - 1)^2} $$

Phonopyでの実装

import numpy as np
import matplotlib.pyplot as plt

# 温度範囲
temperatures = np.arange(0, 1000, 10)

# 熱力学特性を計算
si_phonon.set_mesh([20, 20, 20])
si_phonon.set_thermal_properties(t_step=10, t_max=1000, t_min=0)

tp_dict = si_phonon.get_thermal_properties_dict()
temps = tp_dict['temperatures']
free_energy = tp_dict['free_energy']  # kJ/mol
entropy = tp_dict['entropy']  # J/K/mol
heat_capacity = tp_dict['heat_capacity']  # J/K/mol

# プロット
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 自由エネルギー
axes[0,0].plot(temps, free_energy, linewidth=2)
axes[0,0].set_xlabel('Temperature (K)', fontsize=12)
axes[0,0].set_ylabel('Free Energy (kJ/mol)', fontsize=12)
axes[0,0].set_title('Helmholtz Free Energy', fontsize=14)
axes[0,0].grid(alpha=0.3)

# エントロピー
axes[0,1].plot(temps, entropy, linewidth=2, color='orange')
axes[0,1].set_xlabel('Temperature (K)', fontsize=12)
axes[0,1].set_ylabel('Entropy (J/K/mol)', fontsize=12)
axes[0,1].set_title('Entropy', fontsize=14)
axes[0,1].grid(alpha=0.3)

# 比熱
axes[1,0].plot(temps, heat_capacity, linewidth=2, color='green')
axes[1,0].axhline(3*8.314, color='red', linestyle='--', label='Dulong-Petit (3R)')
axes[1,0].set_xlabel('Temperature (K)', fontsize=12)
axes[1,0].set_ylabel('Heat Capacity (J/K/mol)', fontsize=12)
axes[1,0].set_title('Heat Capacity at Constant Volume', fontsize=14)
axes[1,0].legend()
axes[1,0].grid(alpha=0.3)

# 内部エネルギー（F = U - TSから計算）
internal_energy = free_energy + temps * entropy / 1000  # kJ/mol
axes[1,1].plot(temps, internal_energy, linewidth=2, color='purple')
axes[1,1].set_xlabel('Temperature (K)', fontsize=12)
axes[1,1].set_ylabel('Internal Energy (kJ/mol)', fontsize=12)
axes[1,1].set_title('Internal Energy', fontsize=14)
axes[1,1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('si_thermodynamics.png', dpi=150)
plt.show()

# デバイ温度の計算
# 比熱が3R（Dulong-Petit極限）の63%に達する温度
R = 8.314  # J/K/mol
target_cv = 0.63 * 3 * R
idx = np.argmin(np.abs(heat_capacity - target_cv))
theta_D = temps[idx]
print(f"デバイ温度: {theta_D:.1f} K")
print(f"実験値（Si）: 645 K")

4.5 熱膨張係数

準調和近似（Quasi-harmonic approximation）

体積を変えてフォノン計算を複数回実行：

$$ \alpha(T) = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_P = -\frac{1}{V}\frac{(\partial^2 F/\partial T \partial V)}{(\partial^2 F/\partial V^2)} $$

実装例

import numpy as np
import matplotlib.pyplot as plt

# 異なる格子定数で計算
lattice_constants = np.linspace(5.35, 5.51, 5)  # Å
phonons = []

for a in lattice_constants:
    print(f"Calculating phonon for a = {a:.3f} Å...")
    phonon = calculate_phonon('Si', a=a, dim=(2,2,2))
    phonons.append(phonon)

# 各温度、各体積で自由エネルギーを計算
temps = np.arange(0, 1000, 50)
volumes = (lattice_constants / 5.43)**3  # 規格化体積
free_energies = np.zeros((len(temps), len(volumes)))

for i, phonon in enumerate(phonons):
    phonon.set_mesh([20, 20, 20])
    phonon.set_thermal_properties(t_step=50, t_max=1000, t_min=0)
    tp = phonon.get_thermal_properties_dict()
    free_energies[:, i] = tp['free_energy']

# 各温度で平衡体積を見つける（F最小）
V_eq = np.zeros(len(temps))
for i, T in enumerate(temps):
    # 2次多項式フィット
    coeffs = np.polyfit(volumes, free_energies[i], 2)
    V_eq[i] = -coeffs[1] / (2 * coeffs[0])  # 極値

# 熱膨張係数
alpha = np.gradient(V_eq, temps) / V_eq

# プロット
plt.figure(figsize=(8, 6))
plt.plot(temps, alpha * 1e6, linewidth=2)
plt.xlabel('Temperature (K)', fontsize=12)
plt.ylabel('Thermal expansion coefficient (10$^{-6}$ K$^{-1}$)', fontsize=12)
plt.title('Si Thermal Expansion (QHA)', fontsize=14)
plt.grid(alpha=0.3)
plt.savefig('si_thermal_expansion.png', dpi=150)
plt.show()

print(f"室温（300K）での熱膨張係数: {alpha[6]*1e6:.2f} × 10⁻⁶ K⁻¹")
print(f"実験値: 2.6 × 10⁻⁶ K⁻¹")

4.6 フォノン計算の応用

1. 格子熱伝導率

Boltzmann輸送方程式からの推定（簡易版）:

$$ \kappa = \frac{1}{3} \sum_{\mathbf{q},j} C_{\mathbf{q}j} v_{\mathbf{q}j}^2 \tau_{\mathbf{q}j} $$

• $C_{\mathbf{q}j}$: モードごとの比熱
• $v_{\mathbf{q}j}$: 群速度（$\partial\omega/\partial\mathbf{q}$）
• $\tau_{\mathbf{q}j}$: 緩和時間（フォノン散乱）

PhonoPyでの群速度計算:

# 群速度の計算
si_phonon.set_group_velocity()
group_velocities = si_phonon.get_group_velocity()

print("Group velocities at Gamma point:")
print(group_velocities[0])  # [band, cartesian_direction]

2. 超伝導の臨界温度（Tc）

McMillan式（簡易版）:

$$ T_c = \frac{\omega_{\text{log}}}{1.2} \exp\left[-\frac{1.04(1+\lambda)}{\lambda - \mu^*(1+0.62\lambda)}\right] $$

• $\omega_{\text{log}}$: 対数平均フォノン周波数
• $\lambda$: 電子-フォノン結合定数
• $\mu^*$: クーロン擬ポテンシャル

3. 相転移の検出

虚数フォノンモードが存在すると、結晶構造が不安定：

# 負の振動数（虚数モード）のチェック
frequencies = si_phonon.get_frequencies(q=[0, 0, 0])
if np.any(frequencies < -1e-3):
    print("Warning: Imaginary phonon modes detected!")
    print("Structure may be unstable.")

4.7 本章のまとめ

学んだこと

1. 格子振動の理論 - 調和近似 - 動力学行列 - フォノン分散関係

2. フォノンの分類 - 音響モード vs 光学モード - 縦波 vs 横波

3. Phonopyによる実践 - フォノンバンド構造 - フォノン状態密度 - 完全自動化スクリプト

4. 熱力学特性 - 自由エネルギー - 内部エネルギー - エントロピー - 比熱 - デバイ温度

5. 応用 - 熱膨張係数（準調和近似） - 熱伝導率 - 相転移の検出

重要なポイント

• フォノンは量子化された格子振動
• 調和近似で十分な精度（多くの場合）
• 熱力学特性はフォノンから計算可能
• 実験との良い一致（Si: デバイ温度、熱膨張係数）

次の章へ

第5章では、DFT計算と機械学習を統合した最新手法を学びます。

演習問題

問題1（難易度：easy）

音響モードと光学モードの違いを、物理的意味とともに説明してください。

問題2（難易度：medium）

デバイ温度$\theta_D$の物理的意味と、比熱との関係を説明してください。

問題3（難易度：hard）

準調和近似（QHA）で熱膨張係数を計算する際、なぜ複数の体積でフォノン計算が必要なのか説明してください。

データライセンスと引用

使用したソフトウェア

1. Phonopy - Phonon calculation tool (BSD 3-Clause) - フォノン計算ソフトウェア - URL: https://phonopy.github.io/phonopy/ - 引用: Togo, A., & Tanaka, I. (2015). *Scr. Mater.*, 108, 1-5.

2. GPAW (GPL v3) + ASE (LGPL v2.1+) - DFT力計算バックエンド - 第2章と同様のライセンス

計算パラメータ

• 変位振幅: 0.01 Å（調和近似の標準値）
• Supercellサイズ: 2×2×2以上（収束必須）
• q点メッシュ: 20×20×20以上（DOSとenthalpy計算）

コード再現性チェックリスト

環境構築

conda create -n phonon python=3.11
conda activate phonon
conda install -c conda-forge phonopy gpaw ase matplotlib
pip install spglib  # 対称性解析（推奨）

計算時間の見積もり

トラブルシューティング

問題: 虚数フォノン（負の周波数²） 解決:

# 1. 構造最適化が不十分
# → 力の閾値を0.01 eV/Å以下に
# 2. Supercellが小さい
# → 2×2×2 → 3×3×3に拡大
# 3. 対称性の破れ
# → spglibで対称性をチェック

実践的な落とし穴と対策

1. Supercellサイズ不足

# ❌ 不十分: 1×1×1（収束しない）
phonon = Phonopy(unitcell, [[1, 0, 0], [0, 1, 0], [0, 0, 1]])

# ✅ 推奨: 2×2×2以上（収束テスト必須）
for dim in [(2,2,2), (3,3,3), (4,4,4)]:
    phonon = Phonopy(unitcell, np.diag(dim))
    # フォノン周波数の収束をチェック

2. 構造最適化の不完全

# ❌ 不十分: 粗い最適化 → 虚数フォノン
opt.run(fmax=0.1)  # 力が大きすぎ

# ✅ 推奨: 高精度最適化
opt.run(fmax=0.01)  # 0.01 eV/Å以下

3. q点メッシュ不足（DOS計算）

# ❌ 粗い: 10×10×10（構造を見逃す）
phonon.set_mesh([10, 10, 10])

# ✅ 推奨: 20×20×20以上
phonon.set_mesh([20, 20, 20])

4. 単位換算ミス

# Phonopy内部単位: THz
# 変換:
# - cm⁻¹ = THz × 33.356
# - meV = THz × 4.136
# - K (温度) = THz × 47.992

品質保証チェックリスト

フォノン計算の妥当性

• [ ] Γ点（k=0）の音響モードが0 Hz（誤差 < 0.1 THz）
• [ ] 音響モードが3本（並進自由度）
• [ ] すべての周波数が実数（虚数モードなし）
• [ ] 高対称点でバンドが縮退
• [ ] DOSの総和が3N（自由度数）

熱力学特性の妥当性

• [ ] 高温極限で比熱 → 3R（Dulong-Petit）
• [ ] デバイ温度が文献値の±20%以内
• [ ] 熱膨張係数が正（通常の材料）
• [ ] エントロピーが温度とともに単調増加

参考文献

1. Dove, M. T. (1993). *Introduction to Lattice Dynamics*. Cambridge University Press.

2. Togo, A., & Tanaka, I. (2015). "First principles phonon calculations in materials science." *Scripta Materialia*, 108, 1-5. DOI: 10.1016/j.scriptamat.2015.07.021

3. Phonopy Documentation: https://phonopy.github.io/phonopy/

4. Shulumba, N., et al. (2017). "Temperature-dependent elastic properties of Ti$_x$Zr$_{1-x}$N alloys." Applied Physics Letters, 111, 061901.

著者情報

作成者: MI Knowledge Hub Content Team 作成日: 2025-10-17 バージョン: 1.0 シリーズ: 計算材料科学基礎入門 v1.0

ライセンス: Creative Commons BY-NC-SA 4.0