DFTの実践的基礎理論

2024/10/18にPFNが開催したDFT(Density Functional Theory: 密度汎関数理論)勉強会の資料です。

Transcript

1. DFTの実践的基礎理論 2024/10/18 DFT勉強会 株式会社 Preferred Networks リサーチャー 永井瞭
2. 2 1. DFTの概要と基礎理論 ◦ KS方程式、DFTで計算できるもの 2. 理論と実計算をつなぐ技法 ◦ 基底関数, 擬ポテンシャル, … 3. 交換相関汎関数の詳細 ◦ 汎関数の構成法、自己相互作用, … 4. 高精度化手法の比較 ◦ 汎関数の選択の慣例, +U (Hubbard correction), vdW correction, … Index
3. 3 永井瞭 株式会社 Preferred Networks リサーチャー Matlantisのコア部分の開発や 新規シミュレーション技術の研究を担当 経歴 東京大学理学系研究科物理学専攻　博士課程修了 (物性研究所 杉野研究室) University of Irvine, Kieron Burke group (研究留学) 研究テーマ: 機械学習汎関数 興味 DFTプログラムの実装(分子, 固体どちらも) DFTの各種シミュレーション技法の物理的解釈 機械学習モデルの物質シミュレーションへの埋め込み 自己紹介
4. 4 1. DFTの概要と基礎理論
5. 5 ミクロな現象、電気・熱的性質をシミュレーションするには 原子・電子の挙動が重要。特に、電子は量子力学 例: 分子どうしの化学反応、半導体の電気特性など →電子の量子力学的なエネルギーや波動関数、電子密度を計算することでこれらをシ ミュレーションすることができる • 電子+原子のエネルギーを最小化する方向に原子を動かすorフォノン計算 • 波動関数から電子の動きに関わる量を計算 物性計算における量子力学
6. 6 量子力学ベースの手法 (Born-Oppenheimer近似) • 系の記述: 電子の波動関数 Ψ(r 1 , r 2 ,...r N ), 原子位置は固定 • 支配方程式: Schrödinger方程式 • ポテンシャル: 原子核による背景ポテンシャル+ 電子のクーロン 相互作用 • 数学的性質: 固有方程式 原子レベルのシミュレーション手法 electron + nucleus electron wavefunction nucleus 古典力学ベースの手法 • 系の記述: 質点としての原子の位置 {R i }, 電子は顕に扱わない • 支配方程式: F = m d2R i /dt2 • ポテンシャル: 電子の効果を織り込んだ原子間ポテンシャルを フィッティングして使用 • 数学的性質: 時間発展方程式
7. 7 Schrödinger方程式 電子の波動関数 Ψ(r 1 , r 2 ,...r N ) はN次元空間上の量　扱うのにO(VN)の計算量がかかる 現代のコンピュータでは実用的な電子数のシミュレーションができない →より効率的に解くため密度汎関数理論(DFT)が提唱 電子の子力学量方程式
8. 8 • 密度汎関数理論, Density Functional Theory ◦ 基底状態のすべての物理量が電子密度 n(r) で記述できることを証明 ◦ 電子密度は多体波動関数よりも遥かに小さい空間(空間サイズの1乗)で記述できるの で、数値的にも低コストで扱える可能性 DFTとKohn–Sham equation • Kohn-Sham 方程式: Schrödinger方程式と同じn(r)とEが 多項式計算量で得られる • 解であるψ, nが方程式自体に必要 ◦ 自己無撞着場 (SCF) • 固有値方程式, 計算コストは基底関数の約3乗 • 計算量が少なくなった代わりに, E xc や V xc , といった 厳密に書き下せない項が発生 ◦ “交換相関汎関数, exchange-correlation(XC) functional”. → 近似する必要あり (cf. MD力場)
9. 9 入力 原子位置 (+ セル構造) + 計算パラメータ 出力 エネルギー、電子密度 KS 固有値 & 軌道 エネルギー • 安定性 • 反応障壁 力 • (Ab initio) MD • 力場のフィッティング 電子構造 • 電荷分布 • バンド • 磁性 DFT計算の入出力 The DFT (B3LYP) calculated relative energy barrier for the reductive elimination of difluorinated products from the four coordinate PdII complex [Computational and Theoretical Chemistry 1197, 113162 (2021)] Aligned TDOS for the CH3NH3PbI3 systems under various strains [scientific reports 8, 7760 (2018)]
10. 10 Kohn-Sham equation Walter Kohn, 1998年ノーベル化学賞 W. Kohn and L. J. Sham, “Self-Consistent Equations Including Exchange and Correlation Effects” Phys. Rev. (1965) cited >70000 times Approximated XC functional PBE: cited >170,000 times B3LYP: cited >110,000 times DFTにより量子化学計算が実用レベルに至ったが、 それにより手法の複雑さが跳ね上がった What is DFT? – Impact
11. 11 2. 理論と実計算をつなぐ技法
12. 12 XC汎関数の厳密な形は知られていない。 →様々な近似式が提案 (PBE, B3LYPなど) Exchange-Correlation (XC) functional 精度と単純さ (= 計算コストの軽さ)とのトレードオフ Jacob’s ladder [Chemical And Process Engineering, 1-58 (2013)] (hybrid, N^4) (double hybrid, N^5) Molecular Physics 115, 19 (2017)
13. 13 ψ(r) や n(r) は点ではなく場である。数値的に取り扱うには… → 基底関数で展開 B: 基底関数の数 多くの基底関数を使えばより高精度に表現できるが、KS方程式の計算コストが O(B3)で増 大する 主要な基底関数: Gaussian Type Orbital (GTO) φ =xiyjzkexp(-ar2), 原子中心 例: STO-3G, 6-31G, CC-PVDZ, … cf. Basis set exchange ガウス関数は直交してないので, “基底セット”が予め用意されている Plane Wave (PW) φ = exp(iGr) 使用する平面波の最大の波数をエネルギー指定 e.g. 200eV. 基底関数系
14. 14 DFT プログラム: GTO と PW の比較 GTO PW 代表的なプログラム Gaussian, ORCA, PySCF, GAMESS, MOLPRO,... VASP, Quantum Espresso, GPAW, WIEN2k… 得意 孤立分子など 周期系 不得意 周期系 孤立系, 全電子計算 固有の技術 ガウス積分 ブロッホの波動関数(k点) 擬ポテンシャル 基底関数 6-31G, cc-pvtz, def2-tzvp, … カットオフで指定 高精度WFT CC, CI, CASSCF,... RPA Red: Commercial License, Blue: Free List of quantum chemistry and solid-state physics softwares
15. 15 孤立系では状態の数が離散的: ψ n (r): n=1,2,...,N elec /2 無限に続く周期系では無限の状態が発生 • k = (k x , k y , k z ) はブリルアンゾーン内の任意の実ベクトル • ブリルアンゾーン内で離散的な数値グリッドをとって扱う: k点 • k点の数により計算量が増える。計算速度と精度のトレードオフ • 横軸をk, 縦軸をKS固有値としたプロットがBand structure k点(ブロッホの定理) where u(r) = u(r+R), R: cell vector
16. 16 原子核付近のKS orbitalを滑らかなものに置き換える → 小さい波数の平面波だけで表現できるようになり、cutoffが減らせる →計算量削減 多くのDFTソフトは固有の ppを提供している。 VASPなどはPPの優秀さで有名 cf. Frozen core 近似 (軌道の順番的に内側の軌道を固定) 擬ポテンシャル, Pseudo potential (PP)
17. 17 DFT プログラムでのパラメータ設定 例: Quantum Espresso (PW basis) Calculation type: SCF, MD, Structural optimization,... cell構造 平面波のカットオフ SCFを収束させるためのテクニックの指定 x,y,z 方向のk点の数の指定 擬ポテンシャル(各元素ごとに指定) 原子座標 Gaussian (GTO basis)
18. 18 3. 交換相関汎関数の詳細
19. 19 Hartree–Fock (HF) method 多体波動関数が1電子軌道{φ i }のスレーター行列で書けると仮定し電子間クーロン項を計算 E el-el = E H + E X HF E X HF は (厳密) 交換エネルギーと呼ばれ、波動関数の反対称性が由来 • HFで解くべき方程式はKS方程式のE XC をE X HF に置き換えただけ • 計算コスト: O(N4) ◦ 最適化されたGTO コードでは N3程度に抑えられている ◦ DFT よりもコストが重いのにも関わらず精度が悪いことが多い ▪ 相関が入っていないため Hartree–Fock法と交換相互作用
20. 20 相関 HF近似に含まれないエネルギー • 運動エネルギーとクーロンエネルギーの多体項が由来 Post-HF method HF法で得られた軌道を基底とし、さらに多体の波動関数による効果を計算する • 摂動論: MP2 O(N5) • Coupled Cluster: CCSD(T)...O(N^7) • Many-body perturbation theory: RPA…O(N^4) • Configuration Interaction (CI) ◦ Full CI…O(V^N) Post-HF 法 と相関相互作用
21. 21 KS-DFT におけるXC energy E XC は XC energy density ε XC を用いて定義される • ε XC [n](r) (r における交換相関相互作用の強さ) × n(r) (r における電子数) • E XC は電子密度分布[n]に対して原理上1意に決まるが, {ε XC [n](r)} の決め方は一意ではな い。 • ε XC は ε X とε C に分解される。HF法によりε X の厳密系は知られているが、ε X も近似する • ε XC [n](r) は rのまわりの電子密度記述子を変数とした関数で近似される ◦ eg. n(r), ∇n(r), (n ↑ (r) - n ↓ )/n(r), … • energyそのものではなくenergy densityを近似することで、空間拡張性を担保 (cf. Behler–Parrinello NNP, CNN) DFTにおける交換相関項
22. 22 Semi-local XC functional: - ε xc (r) を 近傍の電子密度の記述子の関数として定義 - 一様電子ガスからの漸近展開など、物理的な条件から 局所密度近似　Local spin density approximation (LDA or LSDA) • 記述子: • 一様電子ガスの XC エネルギーを再現 • 交換項については厳密な形が知られている: • 後述するGGA, metaGGAは次のように定義されることが多い ◦ ε xc = ε xc LDA × (enhancement factor)　 • VWN, PW92 などが有名 一般化勾配近似 Generalized gradient density approximation (GGA) • 追加の記述子:∇n (を無次元化した s) • PBE: 物理条件を組み合わせ単純に作ったもの、固体物理で有名 • BLYP: 分子系でよく使われる Semi-local XC (1)
23. 23 • meta-GGA ◦ 追加の記述子: 運動エネルギー密度 ◦ τにより価電子領域の電子が金属的か、共有結合的か、1電子 的かを判別して補正 ◦ “SCAN” (もしくはr2 SCAN)などが有名 利点 • 計算コストが低い (~ N basis 2) • 数値的に安定で、SCFが収束しやすい 問題点 • 単純なので定量性が低い場合がある ◦ 数kcal/molの違いが重要な分野では結果が反転することもあ るため、より高精度な方法が好まれる • 非局所な効果も重要なはずだが、含まれていない。交換相関ホー ルはほぼデルタ関数 Semi-local XC (2) semi-local 汎関数のXC holeの概 形 [PRL 115, 036402 (2015)]
24. 24 • 1電子しかない場合、電子間クーロン項 E H + E XC は 0であるべき ◦ LDA, GGAなど(semi) localな記述子だけから決定する方法だと、0にならな い ▪ 電子が自分自身と反発することに ▪ 結果、電子は偏在しすぎる ◦ HF 法では完全に打ち消す • 多体系でも自己相互作用の問題は発生する ◦ bandgapや活性化障壁の誤差の主な原因 ◦ 多体系では, HF は逆に自己相互作用を打ち消しすぎる ▪ これにより電子は局在化しすぎる 自己相互作用 repulse itself
25. 25 (Coupling-constant-averaged) XC hole により汎関数のデザインが理解できる “位置 r の電子が位置 r’の電子によりどれくらい強く交換相関効果を受けているか” ほとんどの物質で, 　 の長距離部分は　　の長距離部分と打ち消される =スクリーニング効果 交換相関ホール 典型的な交換、相関 hole の概形 semi-local 汎関数のXC hole
26. 26 4. 高精度化手法の比較
27. 27 • DFT ◦ semi-local XC ▪ LDA ▪ GGA ▪ meta-GGA ◦ hybrid functional ◦ (double hybrid functional) ◦ range-separated hybrid ◦ Hubbard U correction ◦ vdW correction ▪ vdW-D ▪ vdW-DF ◦ (Machine learning based) DFTの精度に関わる理論一覧 (+α) • WFT ◦ HF ◦ MP2 ◦ coupled cluster (CCSD(T)) ◦ RPA(ACFDT) ◦ Quantum Monte Carlo (Diffusion Monte Carlo)
28. 28 E xc = α E X HF + (1-α) E X semi-local + E c semi-local • HFのexchangeをsemilocalに αの割合で混ぜる (0<α<1) • スクリーニング効果を再現するために導入 • PBE0 や B3LYP が有名 Hybrid functional Schematic image of XC hole of exact, HF, semilocal, and hybrid XC. 利点 • 自己相互作用問題が改善 • 多くの場合でsemi-local XCよりも改善 ◦ 反応活性化障壁 ~ 3kcal/mol (過大評価傾向) ◦ バンドギャップも改善 (過大評価傾向) 欠点 • αによる線形混合に物理的根拠はない ◦ これにより追加の物質依存パラメータαが発生 • 計算コストがO(N4)でDFT本体よりも重い。特にPWコードで重い • SCF 収束が悪い場合がある mixing (double hybrid) E xc = α E X HF + (1-α) E X semi-local + β E c MP2 + (1-β)E c semi-local
29. 29 Range-separated hybrid functional hybrid functionalにおいて、E x HF の積分範囲を short-range(SR) と long-range(LR) に分割し どちらかを落とす SR LR HSE (半導体でよく使われる) wB97 (分子系でよく使われる) J.Chem.Phys, 118, 18 J.Chem.Phys, 128, 084106 E HF のLRを落としてSRだけを用いている E HF のSRを落とすor比率を下げている
30. 30 Range-separated hybrid functional 利点 • ωは汎関数の非局所性を簡単に制御できる追加のパラメータのような役割を持 ち、単なるhybrid functionalよりも精度が向上 • 計算コスト・数値的安定性も改善 問題点 • カットオフに使われる erf (erfc) は計算が単純だから採用されているだけで、 物理的な根拠はない • カットオフ距離 1/ω という物質依存の追加パラメータが追加で発生 • 半導体に使われるHSE は LRをcutするのに対して分子系で使われる wB97 はLR を保持する。ユニバーサルではない ◦ bulkでは周囲に電子が詰まっているためスクリーニング効果が強く、XC holeが早く減衰する(=LRをcut) ◦ 孤立分子系では周囲の電子が少ないため減衰しにくい (=LRが残る)
31. 31 自己相互作用の補正を、Hubbard補正で知られる項 　で補正 (n i : 軌道iの占有数) 例: 2つの d軌道 があり (coulomb repulsion) >> (hopping energy gain) だとする Experiment Semi-local XC Semi-local XC + U Hubbard +U 補正 反平行スピンがことなる軌道に局 在 (反強磁性秩序) 自己相互作用により電子が押 し出され、遍歴状態が混ざる + +U 項が遍歴状態にペナルティ を与える
32. 32 Hubbard U 補正 利点: • Uを何かしらの物性値が実験値に合うように決めるだけで、他の物理量も実験と 合うようになることがある • 計算が容易 (semi-local functionalの計算とコストはほぼ変わらない) 問題点: • 単純な対症療法である。本来はもっと複雑な補正が必要なはず • Uは物質依存のパラメータ • 占有数のカウントが擬ポテンシャルの実装にも依存する。DFTプログラムごとに 異なる可能性
33. 33 vdW-D method vdW-DF method 利点: 　長距離相互作用が重要な系の精度が向上 問題点: • vdW-D 法は電子密度を用いず原子位置だけから計算 ◦ MDの力場と類似 ◦ パラメータが目的の物質に十分フィットされているか確認する必要がある • vdW-DF は計算コストが重い ◦ 2重グリッド積分 ◦ 収束性が悪い vdW 補正 PRL 92, 246401
34. 34 孤立分子系: B3LYP, M06-2X, wB97X, Coupled Cluster, … • 分子系向けにチューニングされた汎関数で精度を求める傾向 固体: PBE, PBE+U, SCAN, HSE,... • 物理由来のパラメータレスなものが好まれる • referenceとなる実験データが少ないため、精度よりも安定性が好まれ る 表面系、液体など: vdW-D, vdW-DF, … • 分散力補正が重要に よく使われる汎関数
35. 35 • DFT計算に伴うパラメータ（基底関数、汎関数など）は計算したい物質 によって適切に選択する必要がある ◦ 慣例的に使われている手法も実は合理的な選択に基づいている場合 が多い • DFTは交換相関相互作用を自在に設計することで、正攻法である波動関 数理論に比べて効率よく精度を高めている手法と見ることができる • 交換相関汎関数の設計は以下の観点から理解できる ◦ 電子密度の記述子 ◦ 自己相互作用の補正 ◦ 交換相関ホールの再現方法 まとめ
36. 36 PFN Material Discovery チームではエンジニア・リサーチャーを募集しています！ 主な業務内容 • 物質科学と数値アルゴリズムや機械学習を組み合わせた新規シミュレーション 手法開発 • 最新の手法を用いた物質探索 • Matlantisサービス運用, 新機能の実装 キーワード: MD, NNP, 電子状態理論, 結晶構造探索, 反応経路探索 など 採用募集 Careers page 材料リサーチャー 求人詳細
37. Making the real world computable