Phonopy によるフォノン計算

Phonopyは、第一原理計算で得た原子に働く力から力定数を構築し、フォノン分散・フォノンDOS・熱力学量などを計算するためのポストプロセスである。結晶安定性評価、熱物性推定、相転移や自由エネルギー差の議論を、同一の計算骨格で実行できる点が強みである。

参考ドキュメント

• Phonopy 公式ドキュメント（Overview, Formulation, Interfaces, Features） https://phonopy.github.io/phonopy/
• Togo and Tanaka, First principles phonon calculations in materials science (2015)（概説論文） https://arxiv.org/abs/1506.08498
• （日本語）MateriApps 記事：PhonopyとQuantum ESPRESSOを用いたNaClのフォノン計算 https://ma.issp.u-tokyo.ac.jp/app-post/7146
• （日本語）HPCI ソフトウェア情報：Phonopy https://www.hpci-office.jp/for_users/appli_software/appli_phonopy

1. 位置づけとスコープ

• Phonopyの主対象は調和（harmonic）および準調和（quasi-harmonic）レベルのフォノン物性である
• 代表的なアウトプット
  • フォノン分散、フォノンDOS（全DOS・部分DOS）
  • 格子振動の熱力学量：自由エネルギー、比熱、エントロピー
  • モード・グリュナイゼン定数、群速度、熱振動楕円体（平均二乗変位）など
  • （必要に応じて）非解析項補正（LO-TO分裂）や動的構造因子など
• 格子熱伝導やフォノン寿命（3次以上の非調和）は、通常はphono3py等の別枠に委ねる設計である

2. 原理：力定数・力学行列・固有値問題

2.1 調和近似と力定数

原子変位を小さいと仮定し、全エネルギーを平衡位置まわりで2次まで展開する。

$$E\approx E_0+\frac{1}{2}\sum _{i\alpha ,j\beta }{\mathrm{\Phi }}_{i\alpha ,j\beta }{u}_{i\alpha }{u}_{j\beta }$$

力は $F_{j\beta }=-\partial E/\partial u_{j\beta }$ であるから、力定数は

$${\mathrm{\Phi }}_{i\alpha ,j\beta }=-\frac{\partial F_{j\beta }}{\partial u_{i\alpha }}$$

として定義される。Phonopyは主に有限変位法により、この微分を差分で近似して $\mathrm{\Phi }$（力定数行列）を推定する。

2.2 力学行列とフォノン周波数

結晶の並進対称性を用いて、波数 $\mathbf{q}$ ごとの力学行列 $D(\mathbf{q})$ を構成し、

$$\sum _{j\beta }{D}_{i\alpha ,j\beta }(\mathbf{q})e_{j\beta }^{(\nu )}(\mathbf{q})={\omega }_{\nu }^2(\mathbf{q})e_{i\alpha }^{(\nu )}(\mathbf{q})$$

の固有値問題を解く。固有値 ${\omega }_{\nu }(\mathbf{q})$ がフォノン分散を与える。虚数周波数（${\omega }^2<0$）が現れる場合、構造の動的不安定性や計算条件不備（収束不足、スーパーセル不足、緩和不十分）を疑う入口になる。

3. 力定数の作り方：有限変位スーパーセル法の実務

3.1 基本手順

1. 単位胞（十分に緩和された構造）を用意する
2. 指定したスーパーセルを作り、対称性で非等価な原子変位の集合だけを生成する
3. 各変位構造ごとに第一原理計算を行い、原子に働く力を出す
4. 力と変位から力定数を復元し、音響和則などの拘束を必要に応じて適用する
5. 分散・DOS・熱力学量を計算し、収束と妥当性を点検する

3.2 設計パラメータ

• スーパーセルサイズ
  • 力定数の到達距離を切る操作に等しいため、系が長距離相互作用を持つほど大きいサイズが必要になりやすい
• 変位量（displacement amplitude）
  • 小さすぎると数値誤差に埋もれ、大きすぎると非調和が混入するため、中庸な値を採る（材料と計算精度に依存）
• 力の精度
  • 力定数は力の差分から得るため、電子収束（EDIFF）や基底・k点収束が甘いとフォノンが破綻しやすい
• 対称性の扱い
  • VASP側の対称性で力が“平均化”されると望まない挙動になる場合があるため、Phonopy側の対称性と合わせて設計する

4. Phonopyで得られる代表物性と式の見取り図

4.1 フォノンDOSと熱力学量

フォノン周波数 ${\omega }_{\nu }(\mathbf{q})$ が得られれば、調和近似の振動自由エネルギーは

$$F_{\mathrm{vib}}(T)=\sum _{\mathbf{q},\nu }[\frac{1}{2}\hslash {\omega }_{\nu }(\mathbf{q})+k_{B}T\ln (1-e^{-\hslash {\omega }_{\nu }(\mathbf{q})/k_{B}T})]$$

として計算できる。比熱 $C_V$ やエントロピー $S$ も同様に導かれ、熱容量の温度依存や低温極限の挙動（Debye的な立ち上がりなど）を点検できる。

4.2 準調和近似（QHA）で熱膨張・相安定性へ

体積 $V$ をパラメータとしてフォノンを計算し、$F(V,T)$ を得る。与えた圧力 $p$ のもとで

$$G(T,p)=\underset{V}{min}[F(V,T)+pV]$$

として体積最適化を行えば、熱膨張、体積弾性率の温度依存、相の自由エネルギー差による相安定性評価へ展開できる。

5. 極性結晶の注意点：LO-TO分裂（非解析項補正）

絶縁体・半導体で長距離クーロン相互作用が効く場合、$\mathrm{\Gamma }$ 点近傍でLO-TO分裂が生じる。Phonopyは、Born有効電荷と高周波誘電率を用いて非解析項補正（NAC）を取り込み、分散のギブス振動や不自然な分裂欠落を抑える設計を持つ。NACに必要な量は第一原理側（例：VASPの線形応答）から別途求め、Phonopyへ与える。

6. ワークフロー

6.1 入力と受け渡しファイル

• 入力側（VASP）
  • 緩和済みの構造（POSCAR/CONTCAR）
  • 変位構造ごとの力（OUTCAR等に含まれる）
  • （必要なら）Born有効電荷・誘電率（NAC用）
• 受け渡し（Phonopy）
  • 変位構造群（VASPに投げるPOSCAR群）
  • 力定数を作るための集約ファイル（FORCE_SETSやFORCE_CONSTANTS等）
  • 設定ファイル（スーパーセル、バンド経路、メッシュなど）

6.2 運用方針

• 変位構造ごとのVASP計算は、入力条件を完全に固定して同列に扱う（カットオフ、k点、スミアリング、収束条件など）
• 力のノイズが見える場合は、構造緩和の厳密化、k点密度増加、電子収束強化、スーパーセル拡大を優先して疑う
• 虚数モードが出た場合は、物理的不安定性と数値要因を切り分ける（スーパーセル収束、変位量、力精度、残留応力など）

7. 活用パターン

• 動的安定性スクリーニング：新規結晶候補や置換系で虚数モードの有無を高速に判定する
• 熱物性の一次推定：比熱・エントロピー・自由エネルギーを通じて、温度依存の議論を整える
• QHAによる相安定性：温度で結晶相が入れ替わる可能性や、熱膨張の寄与を見積もる
• 実験分光との接続：分散・DOS・動的構造因子を介して、INS/IXS等の解釈を補助する

8. 注意点

• 構造が十分に緩和されていない
  • 残留力・残留応力が大きいと、音響モードや低周波が崩れやすい
• スーパーセルが小さすぎる
  • 力定数の切断誤差が分散に現れる（特に低周波や分枝の形状）
• 力の数値精度が不足している
  • フォノンは力の差分で決まるため、電子収束が甘いと全体が破綻し得る
• 極性結晶でNACを入れていない
  • $\mathrm{\Gamma }$ 近傍のLO-TO分裂が再現されず、比較が難しくなる
• 結果の比較条件が揃っていない
  • 例えば置換系列でPOTCAR/ENCUT/k点が暗黙に変わると、フォノン差が物理由来か判別不能になる

まとめ

Phonopyは、第一原理計算で得た力から力定数を構築し、調和・準調和レベルのフォノン物性を一貫して算出するための基盤である。実務では、スーパーセル設計と力精度、極性結晶における非解析項補正の有無が結果を支配するため、収束設計と点検項目を手順に組み込んで運用するのが有効である。