第一原理計算を支える数値解法

第一原理計算は、変分原理で定義された電子基底状態問題を、離散化と反復解法により数値的に解く枠組みである。計算の成否は、自己無撞着（SCF）の収束制御、巨大固有値問題の反復解法、ブリルアンゾーン積分の近似が支配することが多い。

参考ドキュメント

• VASP Wiki: Self-consistency cycle / Smearing technique / RMM-DIIS https://vasp.at/wiki/Self-consistency_cyclehttps://vasp.at/wiki/Smearing_techniquehttps://vasp.at/wiki/RMM-DIIS
• Quantum ESPRESSO: PWscf User's Guide（反復対角化や混合のパラメタを含む） https://www.quantum-espresso.org/Doc/pw_user_guide/
• （日本語）住友電工 技術論文：電子状態計算におけるJacobi-Davidson法と修正方程式（反復固有値解法の見通し） https://sei.co.jp/technology/tr/bn178/pdf/sei10671.pdf

1. 全体像

典型的な密度汎関数理論（Kohn–Sham DFT）は、次の連立問題に帰着する。

• Kohn–Sham方程式

$${\hat{H}}_{\mathrm{KS}}[\rho ]{\psi }_{n\mathbf{k}}(\mathbf{r})={\epsilon }_{n\mathbf{k}}{\psi }_{n\mathbf{k}}(\mathbf{r})$$

• 電荷密度の再構成

$$\rho (\mathbf{r})=\sum _{n,\mathbf{k}}{w}_{\mathbf{k}}{f}_{n\mathbf{k}}|{\psi }_{n\mathbf{k}}(\mathbf{r}){|}^2$$

ここで本質は、$\rho \mapsto \hat{H}[\rho ]\mapsto \{\psi \}\mapsto {\rho }_{\mathrm{out}}$ という写像の固定点 ${\rho }^{\ast }={\rho }_{\mathrm{out}}({\rho }^{\ast })$ を反復で見つける点にある。

SCFの擬似コードは次で表せる。

• 初期密度 ${\rho }_0$ を与える
• 反復 $m=0,1,2,\dots$：
  1. $\hat{H}[{\rho }_m]$ を構築する
  2. 固有値問題を解き、占有数 $f_{n\mathbf{k}}$ を決める
  3. ${\rho }_{\mathrm{out}}$ を生成する
  4. 混合（mixing）で ${\rho }_{m+1}$ を作る
  5. 収束判定（エネルギー差、残差ノルム、電荷差）を行う

2. 離散化

連続空間の波動関数 $\psi (\mathbf{r})$ を有限次元に落とすため、代表的には次を用いる。

• 平面波基底 + 擬ポテンシャル/PAW
  $\psi$ を平面波で展開し、運動エネルギーでカットオフ $E_{\mathrm{cut}}$ を導入する。 平面波は畳み込みがFFTで高速化でき、周期系で扱いやすい。

• 全電子法（LAPW, APW+lo など）
  核近傍と格子間で基底を切り替え、高精度に全電子密度を扱う。固有値問題は一般化固有値問題になりやすい。

• 実空間差分/有限要素/ウェーブレット
  局所的な格子（メッシュ）で離散化し、境界条件や不均一系に強い設計が可能である。

離散化の観点で重要な実務指標は次である。

• 基底サイズ（未知数の数）$N$ が計算量とメモリを決める
• 固有値問題が密（dense）か疎（sparse/作用素形式）かで、使える反復解法が変わる
• FFTやメッシュ演算の並列化がボトルネックになりやすい

3. 巨大固有値問題

離散化すると、一般に次の固有値問題が現れる。

• 標準固有値問題：$Ax=\epsilon x$
• 一般化固有値問題：$Ax=\epsilon Bx$

ここで $A$ の次元は大きく、全固有値を求める直接法（フル対角化）は高コストである。DFTでは「占有される低エネルギー側の固有対」だけが必要なので、反復固有値解法が主役となる。

3.1 反復対角化

• Davidson法（ブロックDavidson）
  低エネルギー固有対をサブスペース上で近似し、残差を用いてサブスペースを拡張する。ブロック化によりBLAS3（行列積）に寄せられる。

• Lanczos法
  3項漸化式でKrylov部分空間を作る。スペクトル情報を効率よく抽出できるが、再直交化などの工夫が必要になりやすい。

• 共役勾配（CG）型（バンド別の最適化）
  各軌道をエネルギー最小化として更新する。計算は堅牢でメモリが軽い設計に向くが、状況により遅くなることがある。

• Jacobi–Davidson法
  Davidsonを拡張し、補助方程式（修正方程式）を解いてサブスペースを改善する。大規模化で直交化コストが顕在化する点を意識した設計が議論される。

3.2 前処理

固有値反復の収束は、残差 $r=Ax-\epsilon x$ をどれだけ「固有空間に近い方向へ」補正できるかに依存する。そこで

$$\mathrm{\Delta }x\approx M^{-1}r$$

の形で、$M^{-1}$（近似逆）を導入する。平面波では「運動エネルギーが支配する高周波成分」を抑えるような前処理が自然に作りやすい。

4. SCF収束

4.1 固定点反復

最も単純な更新は線形混合である。

$${\rho }_{m+1}=(1-\alpha ){\rho }_m+\alpha {\rho }_{\mathrm{out}}({\rho }_m)$$

しかし金属・表面・長尺セルでは、長波長成分が過度に増幅される電荷のゆらぎ（charge sloshing）が起きやすく、発散・停滞の原因になる。

4.2 混合法

SCFは非線形方程式

$$F(\rho )={\rho }_{\mathrm{out}}(\rho )-\rho =0$$

の求根として見なせる。代表的な加速は次である。

• Pulay（DIIS）
  過去の反復の残差を用い、残差が最小になる線形結合で次の密度（またはポテンシャル）を外挿する。

• Anderson mixing
  DIISと本質的に近い枠組みとして整理されることが多い。

• Broyden法（準Newton型）
  ヤコビアン（またはその逆）の近似を更新し、

  $${\rho }_{m+1}={\rho }_m-B_m^{-1}F({\rho }_m)$$

  の形で高速化する。有限記憶（limited-memory）化により大規模系でも扱える。

4.3 前処理：Kerker型

電荷の長波長成分（小さい $|\mathbf{G}|$）が暴れる系では、混合の前に逆空間で残差にフィルタをかける。

代表例としてKerker前処理は

$$P(\mathbf{G})=\frac{|\mathbf{G}{|}^2}{|\mathbf{G}{|}^2+G_0^2}$$

のように小$|\mathbf{G}|$を抑える。$G_0$は遮蔽長（Thomas–Fermi的なスケール）に結びつけて選ばれることが多い。

4.4 金属の占有数

金属ではフェルミ面近傍の占有が不連続であるため、有限温度（または人工的広がり）で占有数を滑らかにする。

$$f_{n\mathbf{k}}=\frac{1}{1+\exp \left(\frac{{\epsilon }_{n\mathbf{k}}-\mu }{k_{\mathrm{B}}T}\right)}$$

これによりSCFが安定化する一方、積分誤差や自由エネルギーの扱い（$T\to 0$外挿、$E_0$評価など）を一貫させる必要がある。

5. ブリルアンゾーン積分：k点サンプリングと評価法

周期系の物理量は

$$⟨A⟩={\int }_{\mathrm{BZ}}A(\mathbf{k})d\mathbf{k}\approx \sum _{\mathbf{k}}{w}_{\mathbf{k}}A(\mathbf{k})$$

で評価される。ここで誤差は次で決まる。

• k点数と対称性低減
• 金属か絶縁体か（フェルミ面の有無）
• 評価法（テトラヘドロン法、スメアリング法）

5.1 テトラヘドロン法（Blöchl補正を含む）

テトラヘドロン法はk空間を四面体に分割し、帯域の線形内挿で積分を行う。絶縁体の全エネルギーやDOSで精度が高い設計として知られる。金属向けには補正（Blöchl補正など）が提案されている。

5.2 スメアリング法（Gaussian, Methfessel–Paxton など）

スメアリングはフェルミ面の不連続を滑らかにし、少ないk点で収束を得やすくする。手法により「自由エネルギーの誤差次数」や「占有の非単調性」が異なるため、目的（緩和、DOS、微小なエネルギー差）に応じて選ぶのがよい。

6. 力・応力と構造最適化

原子座標 $\mathbf{R}$ に対し、全エネルギー $E(\mathbf{R})$ を最小化する。基本量は力である。

$${\mathbf{F}}_I=-\frac{\partial E}{\partial {\mathbf{R}}_I}$$

理想的にはHellmann–Feynman力が主項であるが、基底が $\mathbf{R}$ に依存する場合はPulay項が現れる。セル最適化では応力（Pulay stressを含む）の制御も必要になる。

構造最適化の数値法は次が代表である。

• 最急降下法：単純だが遅い
• 共役勾配法：メモリ軽量で堅牢
• 準Newton法（BFGS/L-BFGS）：収束が速いことが多い$${\mathbf{R}}_{n+1}={\mathbf{R}}_n-H_n^{-1}\mathrm{\nabla }E({\mathbf{R}}_n)$$ここで $H_n^{-1}$ はヘッセ行列逆の近似である。

7. 線形応答・励起

フォノンや誘電応答などでは、摂動に対する一次応答を求める必要がある。代表的な枠組みは次である。

• 有限差分：超胞を作り原子を微小変位させ、力定数を数値微分で得る
• DFPT（密度汎関数摂動論）：Sternheimer方程式などの線形方程式を解き、直接応答を得る

いずれも、SCFの内側に「追加の反復／線形ソルバ」を持つため、収束設計が二重になる点が特徴である。

8. 計算量・スケーリング

代表的な支配要因は次である。

• 反復対角化：作用素適用（FFT）、直交化、サブスペース更新
• k点並列：k点は独立なので並列化しやすいが、各k点の固有値反復が重い
• 混合法：履歴長（mixing_ndim相当）と保存量がメモリに効く
• 大規模化：$O(N^3)$成分（直交化や密行列操作）が顕在化しやすい

大系では、次の方向性が用いられる。

• 作用素形式の徹底（疎/FFT中心）
• 線形スケーリング（密度行列法、分割統治、Fermi演算子展開など）
• 近似の階層化（粗いSCF→高精度SCF、マルチグリッドや多重解像度）
• GPU/通信隠蔽などの実装最適化（FFTと線形代数の比率調整）

9. 数値課題

数値課題	現れる場所	代表解法	典型的な失敗モード
固定点反復の発散・停滞	SCF	線形混合、DIIS、Anderson、Broyden	charge sloshing、金属での揺れ
長波長成分の暴走	SCF	Kerker前処理、誘電モデル前処理	表面・真空・長尺セルで顕著
低エネルギー固有対の抽出	KS方程式	ブロックDavidson、RMM系、CG、Jacobi–Davidson	直交化コスト、メモリ不足
k点積分誤差	BZ積分	Monkhorst–Pack、テトラヘドロン、スメアリング	金属で収束が遅い
微小エネルギー差	異方性・相対安定性	高密度k点、厳密な収束、評価法の統一	数μeV〜meVの誤差に埋もれる
力・応力の誤差	最適化	収束条件強化、基底依存の注意	形だけ収束、構造が揺れる

まとめ

第一原理計算の数値解法は、(1) 固有値問題を反復で解く設計、(2) SCFの固定点反復を混合と前処理で安定化する設計、(3) ブリルアンゾーン積分の評価法を目的に応じて選ぶ設計の三点に整理できる。これらは互いに独立ではなく、占有数（スメアリング）やk点密度、前処理の選択がSCFと固有値反復の両方に影響するため、全体を一つの反復最適化問題として捉えるのが有効である。