2023 Volume 87 Issue 1 Pages 1-17
Revealing atomic-scale distributions of energy and stress in defective or complex systems, based on the behavior of electrons, should contribute much to materials science and engineering, while only few practical ab initio methods were developed for this purpose. Thus, we developed computational techniques of local-energy and local-stress calculations within the plane-wave PAW (projector augmented wave)-GGA (generalized gradient approximation) framework. This is natural extension of ab initio energy-density and stress-density schemes, while the inherent gauge dependency is removed by integrating these densities in each local region where the contained gauge-dependent terms are integrated to be zero. In this overview, we explain our scheme with some details and discuss the concepts or physical meanings of local energy and local stress via the comparison with related schemes using similar densities, LCAO (linear combination of atomic orbitals) methods, Green’s function formulation implemented by multiple scattering or TB (tight-binding) methods, or EAM (embedded-atom method) potentials. We present recent applications to metallic surfaces, grain boundaries (GBs) with and without solute segregation, tensile tests of metallic GBs, local elastic constants of microstructures in alloys, and machine-learning based GB-energy prediction, where the local-energy and local-stress analyses provide novel aspects of phenomena, deep insights into the mechanism, and effective data for novel machine-learning based modelling. We discuss unsettled issues and future applications, especially for large-scale metallic systems.
Mater. Trans. 62 (2021) 1-15に掲載
①局所密度近似（local density approximation; LDA）や一般化勾配近似（generalized gradient approximation; GGA）を用いた密度汎関数理論（density-functional theory; DFT）1-3）の確立と，②平面波基底の擬ポテンシャル法など4-7）に基づく大規模スーパーセルの電子状態計算の高速計算技法の開発8-10）により，様々な結晶，欠陥，表面，界面，粒界，合金などの現実的な構造の第一原理計算が可能となった．本稿では，①②を合わせた計算技術を「平面波基底の第一原理計算法」（plane-wave DFT method）と総称する．欠陥や粒界・界面は，大きな結晶の中に欠陥などが存在し，周期的に繰り返す「スーパーセル」として扱われる．こうして，様々な系の構造や物性が電子・原子の振る舞いから高精度に解明され，材料科学・工学に大きく貢献している．材料の機械的性質について，第一原理計算を用いた研究のレビューが出ている11-14）．また，VASP，CASTEP，ABINITなどの平面波基底の第一原理計算法の汎用コード開発も，この研究分野の発展に大きく貢献している．
現在の汎用コードでは，スーパーセル内で積分または平均した量として全エネルギーと応力テンソル15,16）が得られるが，スーパーセル内のエネルギーや応力の局所分布は得ることができない．もし，格子欠陥や粒界・界面を含む大きなスーパーセル内のエネルギーや応力の局所的な分布が得られれば，欠陥や粒界・界面の性質が深く解明できる．また，スーパーセル内の欠陥や粒界・界面の周囲の完全結晶部分のエネルギーや応力を調べることで，スーパーセルのサイズ（結晶部分の厚み）が十分かどうかが明確に判断できる．しかし，スーパーセル内のエネルギーや応力の分布を扱う実用的な手法やコードは，ほとんど開発されていない．
そこで筆者らは，最新の平面波基底の第一原理計算法の枠組みとして，PAW（projector augmented wave）6,7）-GGA法に基づき，スーパーセル内の「局所エネルギー」と「局所応力」を求める実用的な計算技術17,18）を開発し，産業技術総合研究所のパッケージソフトウェアQMAS（quantum materials simulator）19）でコード化した．
まず，これまでの平面波基底の第一原理計算法の枠内で提案されていた「エネルギー密度」と「応力密度」の方法20,21）を，PAW-GGA法の枠組みで再定式化した17,19）．これらの密度分布は，通常の第一原理計算を行った後，占有波動関数や自己無撞着な電子密度分布を用いて実空間メッシュ点のデータとして算出され，スーパーセル内の実空間のエネルギーと応力の分布情報を含んでいる20,21）．これらの密度をスーパーセル全体で積分すると，従来法による全エネルギーや応力テンソルが厳密に得られる．しかし，スーパーセル内のエネルギー密度や応力密度のメッシュ点での値や任意の領域で積分した積分値は，ゲージ依存項（スーパーセル全体で積分すると消える不定性の項）の存在のため，一意的な物理量とは見なせない．筆者らは，この問題を解決するために，「ゲージ依存項が積分してゼロになるような条件」で分割された局所領域毎にエネルギー密度と応力密度を積分することで，各領域（原子や原子層，クラスターなど）の局所エネルギーと局所応力を確固とした物理量として計算することを提案した17,18）．筆者らの局所エネルギー・局所応力計算法は，金属表面17,22,23），溶質偏析のある粒界とない粒界18,24-31），金属粒界の引張試験32,33），合金の微細構造の局所弾性定数34），機械学習と組み合わせたモデリング35,36）など，様々な系や現象に適用され，現象の新しい側面やメカニズムが明らかになり，また，局所的環境vs局所エネルギーの機械学習に基づくモデリングに有効なデータを提供できることが示された．このOverviewでは，筆者らの方法と最近の応用例について，いくつかの詳細を含めて紹介し，議論する．局所エネルギーについては，Trinkleと共同研究者により，PAW-GGA法の枠組みの中で非常によく似たスキームも開発されており37），それとの違いについても論じる．
局所エネルギーと局所応力の概念は，材料科学・工学において非常に重要であるが，その理解や検討は十分ではない．EAM（embedded-atom method）38,39）などの古典的な原子間ポテンシャルを用いれば，原子毎のエネルギーや応力の分布を得ることができるが，原子毎の量の定義は一意的ではなく，概念的な困難が存在する．固体中では，波動関数の本質的な広がり，連続的な電子密度分布，長距離の静電相互作用などがあるためである．一方，量子力学的な応力テンソル場は，量子力学の初期に議論され40,41），関連するエネルギー密度（場）・応力密度（場）が固体に対して提案されたが15,20,21），やはり一意性や観測可能性に関する曖昧さがある．一方，LCAO（linear combination of atomic orbitals）法や多重散乱法（KKR（Korringa-Kohn-Rostoker）法），強結合近似（tight-binding; TB）法などに基づくGreen関数法計算により，固体中の欠陥や表面の電子構造を局所的に計算し，それを通じて局所的なエネルギーを議論する取り組みは，長い歴史を持っている42-45）．このOverviewでは，上記の様々な関連する方法やアプローチとの比較を通して，局所エネルギーと局所応力の歴史や物理的な意味についても論じる．なお，本稿で用いる用語について，局所エネルギー，局所応力とは，ある局所空間で定義され，系全体の積分量または平均値に対する，局所的な寄与として計算される量である．
物質・材料中の応力については，注意が必要である．表面，欠陥，アモルファス構造などにおける「原子レベル応力」（atomic-level stress）は，通常の結晶の弾性歪場による「弾性応力」とは，性質も起源も大きく異なる．表面応力では，配位数の少ない表面原子がバルクと異なる面内原子間距離を好むが，面内周期性はバルクのものに制限され，結果として過剰な表面応力が発生する．より詳細には，表面および表面近傍の原子層における電子の再分配や軌道の再混成が，そういう過剰な応力分布を引き起こす46）．局所的な歪がなくても，表面応力は発生するのである．アモルファス構造における「原子レベル応力」も，古典的な原子間ポテンシャルを用いてEgamiとVitek47-49）が明らかにしたように，配位数，原子間距離，原子体積の局所的な変化（deviation）によって発生するもので，歪場によるものではない．同様の原子レベル応力は，筆者らの第一原理計算による局所応力計算によって金属粒界の中心領域でも観測された18,24,25）．粒界の中心領域の構造は，特異な原子環やクラスターからなるアモルファスや転位芯構造に共通する特徴を持つ．最近，EgamiらのグループもGreen関数法による第一原理局所応力計算法を開発したが50），本Overviewではそれとの比較についても論じる．
一方，従来の弾性応力は，結晶中の弾性歪場によって発生し，長い波長を持つ連続的な原子配列の変位から構成される．このような弾性歪や応力は，連続体モデルで記述されるように原子スケールでは極めてゆっくりと変化する．欠陥やアモルファスの原子レベル応力は，隣接する原子間でさえ負から正へと急激に変化するため，連続体モデルでは扱えない．これは，微視的な応力と巨視的な応力の違いによる．後述するように，筆者らの方法も古典的な原子間ポテンシャル法も，原子レベル応力と従来の弾性応力の両方を扱うことができる．
本Overviewでは，第2章で，筆者らの局所エネルギー・局所応力の第一原理計算法17,18）について，若干の詳細を含めて，あらためて説明する．PAW-GGA法の枠組みでのエネルギー密度，応力密度の各項について説明した後，ゲージに依存しない局所領域の条件を明示し，領域分割法の例を示す．第3章では，上でも触れた，最近の類似手法や関連手法（エネルギー密度の別表現やGreen関数法による局所応力の方法）37,50），以前からのLCAO法やKKR法，TB法，EAMポテンシャルを用いた局所エネルギー・局所応力の計算法などと，筆者らの手法との比較を議論し，局所エネルギーと局所応力の歴史と物理的意味について論ずる．第4章では，様々な系や現象への，筆者らの局所エネルギー・局所応力法の適用例を，機械学習技術との組合せも含めて紹介する．第5章では，理論的・計算技術的な未解決問題や今後の展望について論じる．
スーパーセル内の実空間の分布として与えられたエネルギー密度$\varepsilon_{\textit{tot}}(\vec{r})$20），応力密度$\tau_{\alpha \beta }(\vec{r})$21）を，スーパーセル内で空間分割された各々の局所領域で積分することで，局所エネルギー，局所応力が得られる．エネルギー密度，応力密度は，スーパーセルの周期性を持つ実空間の被積分関数で，スーパーセル全体での積分が平面波基底の第一原理計算法による全エネルギーEtot，応力テンソルσαβ15,16）に等しく，次のように定義される．
| \begin{equation} E_{\textit{tot}} = \int_{\Omega}\varepsilon_{\textit{tot}}(\vec{r})d\vec{r} \end{equation} | (1) |
| \begin{equation} \sigma_{\alpha \beta} = \frac{1}{\Omega}\frac{\partial E_{\textit{tot}}}{\partial \varepsilon_{\alpha \beta}} = \frac{1}{\Omega}\int_{\Omega}\tau_{\alpha \beta}(\vec{r})d\vec{r} \end{equation} | (2) |
これらエネルギー密度，応力密度は，ゲージ依存項が存在するため，任意の場所での密度値や任意の領域での積分値は，確固とした物理量とは見なせない20,21）．筆者らの方法17,18）では，含まれるゲージ依存項が各々の領域で積分されてゼロになるように分割された局所領域Ωiを構築し，その局所領域でエネルギー密度，応力密度を積分することで，ゲージ非依存の局所エネルギー$E_{\textit{tot}}^{i}$と局所応力$\sigma_{\alpha \beta }^{i}$を次式のように求める．
| \begin{equation} E_{\textit{tot}}^{i} = \int_{\Omega_{i}}\varepsilon_{\textit{tot}}(\vec{r})d\vec{r} \end{equation} | (3) |
| \begin{equation} \sigma_{\alpha \beta}^{i} = \frac{1}{\Omega_{i}}\int_{\Omega_{i}}\tau_{\alpha \beta}(\vec{r})d\vec{r} \end{equation} | (4) |
| \begin{equation} E_{\textit{tot}} = \sum\nolimits_{i}E_{\textit{tot}}^{i} \end{equation} | (5) |
| \begin{equation} \sigma_{\alpha \beta} = \frac{1}{\Omega}\sum\nolimits_{i}\Omega_{i}\sigma_{\alpha \beta}^{i} \end{equation} | (6) |
従来のエネルギー密度20），応力密度21）は，平面波基底の第一原理計算法のうち，各々NCPP（norm-conserving pseudopotential method4））-LDA法，USPP（ultra-soft pseudopotential method5））-LDA法の枠組みで開発された．筆者らは，それらをPAW6,7）-GGA法の枠組みで再定義した17-19）．PAW法における全エネルギーは，いくつかの項の和であり7），全エネルギー密度$\varepsilon_{\textit{tot}}(\vec{r})$は，対応する各項の密度形の和として以下のように与えられる．
| \begin{align} \varepsilon_{\textit{tot}}(\vec{r}) &= \varepsilon_{\textit{kin}}(\vec{r}) + \varepsilon_{\textit{xc}}(\vec{r}) + \varepsilon_{H}(\vec{r}) + \varepsilon_{\textit{local}}(\vec{r}) \\ &\quad + \varepsilon_{\textit{non-local}}(\vec{r}) + \varepsilon_{\textit{ovrl-self}}(\vec{r}) \end{align} | (7) |
式（7）において，最初の4項は通常の密度形式で表され，最後の2項$\varepsilon_{\textit{non-local}}(\vec{r})$，$\varepsilon_{\textit{ovrl-self}}(\vec{r})$は各原子サイトに存在する量の和で，原子位置$\vec{R}_{I}$毎のデルタ関数$\delta (\vec{r} - \vec{R}_{I})$で空間的に分布させる．$\delta (\vec{r} - \vec{R}_{I})$はイオン半径内部でGaussianによりブロード化し，通常の密度形式とする20）．こうして，$\varepsilon_{\textit{non-local}}(\vec{r})$，$\varepsilon_{\textit{ovrl-self}}(\vec{r})$と$\varepsilon_{\textit{local}}(\vec{r})$の3項は，イオン半径内部の分布の様子には物理的意味はなく，各イオン半径の内側での積分値のみが意味を持つ．応力密度にも同様の問題がある．したがって，エネルギー密度と応力密度を積分する「局所領域」は，各原子サイトのイオン半径の内側を丸ごと含む必要がある．あるいは，対称性のある構造では，イオン球を平面で分割した領域も可能である20,51）．
第1章で触れたように，Trinkleらは，筆者らと同様にPAW-GGA法の枠内で，よく似た局所エネルギー計算法を開発している37）．彼らの全エネルギー密度は同様の項で構成されるが$\varepsilon_{\textit{local}}(\vec{r})$は含まれない．なぜなら，局所擬ポテンシャルがGaussianによるイオン電荷の静電場と打ち消すように作られるからである．
式（7）の各項の詳細は文献17）にある．ここではゲージ依存問題に関わる項を詳しく説明する．まず，価電子の運動エネルギー密度（エネルギー密度のkinetic項）$\varepsilon_{\textit{kin}}(\vec{r})$には，以下のように非対称形と対称形がある：
| \begin{align} \varepsilon_{\textit{kin}}(\vec{r}) &= - \frac{1}{2}\sum_{i}^{\textit{occ}}f_{i}\psi_{i}^{*}(\vec{r})\nabla^{2}\psi_{i}(\vec{r}) \\ & = \frac{1}{2}\sum\nolimits_{i}^{\textit{occ}}f_{i}\nabla \psi_{i}^{*}(\vec{r})\cdot \nabla \psi_{i}(\vec{r}) - \frac{1}{4}\nabla^{2}\rho_{e}(\vec{r}) \end{align} | (8) |
静電エネルギー密度$\varepsilon_{H}(\vec{r})$は，Maxwell形式で以下のように表される20）．
| \begin{equation} \varepsilon_{H}(\vec{r}) = \frac{1}{8\pi}\nabla V_{H}^{tot}(\vec{r})\cdot \nabla V_{H}^{tot}(\vec{r}) \end{equation} | (9) |
一方，静電エネルギー密度は，Coulomb形式で次のようにも表現される．
| \begin{equation} \varepsilon_{H}(\vec{r}) = - \frac{1}{8\pi}V_{H}^{tot}(\vec{r})\nabla^{2}V_{H}^{tot}(\vec{r}) \end{equation} | (10) |
応力密度も同様にいくつかの項の和で表される．
| \begin{align} \tau_{\alpha \beta}(\vec{r}) &= \tau_{\alpha \beta}^{\textit{kin}}(\vec{r}) + \tau_{\alpha \beta}^{\textit{xc}}(\vec{r}) + \tau_{\alpha \beta}^{H}(\vec{r}) + \tau_{\alpha \beta}^{\textit{local}}(\vec{r}) + \tau_{\alpha \beta}^{\textit{non-local}}(\vec{r}) \\ &\quad+ \tau_{\alpha \beta}^{\textit{ovrl-self}}(\vec{r}) + \tau_{\alpha \beta}^{\textit{overlap}}(\vec{r}) \end{align} | (11) |
| \begin{align} E_{\textit{tot}}^{\varepsilon} &= \int_{\det | I + \varepsilon |\Omega}\varepsilon_{\textit{tot}}((I + \varepsilon)^{- 1}\vec{r})d\vec{r} = \int_{\Omega}\eta^{\varepsilon}(\vec{r}')\det | I + \varepsilon |d\vec{r}' \\ & \approx (1 + Tr(\varepsilon))\int_{\Omega}\eta^{\varepsilon}(\vec{r})d\vec{r} \end{align} | (12) |
| \begin{equation} \begin{split} \psi_{i}(\vec{r}) \to \psi_{i}^{\varepsilon}(\vec{r}) = \det | I + \varepsilon |^{- 1/2}\psi_{i}((I + \varepsilon)^{- 1}\vec{r}) \\ \rho_{e}(\vec{r}) \to \rho_{e}^{\varepsilon}(\vec{r}) = \det | I + \varepsilon |^{- 1}\rho_{e}((I + \varepsilon)^{- 1}\vec{r}) \\ V_{a}(\vec{r} - \vec{R}_{a}) \to V_{a}^{\varepsilon}(\vec{r} - \vec{R}_{a}) = V_{a}(\vec{r} - (I + \varepsilon)\vec{R}_{a}) \end{split} \end{equation} | (13) |
式（7）の各エネルギー密度項に対して，この実空間の手順を適用すると，式（11）の応力密度の各項が得られる．なお，いくつかの項や部分は，こうした実空間の手順よりも，逆空間表現を用いる手順の方がより簡単である．また，$\rho_{e}(\vec{r})$に含まれる各原子位置の補償電荷は，上記式（13）とは異なる扱いをする．エネルギー密度には含まれず，応力密度だけの項として，$\tau_{\alpha \beta }^{\textit{overlap}}(\vec{r})$項がoverlap演算子から導出される．
応力密度の各項の詳細は文献17）にある．応力密度のkinetic項$\tau_{\alpha \beta }^{\textit{kin}}(\vec{r})$は以下のように表される．
| \begin{align} \tau_{\alpha \beta}^{\textit{kin}}(\vec{r}) &= \sum_{i}^{\textit{occ}}f_{i}\psi_{i}^{*}(\vec{r})\nabla_{\alpha}\nabla_{\beta}\psi_{i}(\vec{r}) \\ & = - \sum\nolimits_{i}^{\textit{occ}}f_{i}\nabla_{\alpha}\psi_{i}^{*}(\vec{r})\nabla_{\beta}\psi_{i}(\vec{r}) + \frac{1}{2}\nabla_{\alpha}\nabla_{\beta}\rho_{e}(\vec{r}) \end{align} | (14) |
式（7），式（11）の最終形をスーパーセル全体で積分すると，平面波基底のPAW-GGA法で通常に求めたEtot，σαβと厳密に同じ値になり，密度形の妥当性が証明される．エネルギー密度や応力密度における，kinetic項と静電相互作用項の，各々の2つの形式の問題については，スーパーセル全体の積分では違いがなくなる（ゲージ依存項が積分してゼロになる）．
2.3 ゲージに依存しない局所領域への分割法 2.3.1 局所領域のゲージ依存性を除去する条件上記のように，通常の第一原理計算で求めた自己無撞着な電子密度分布や占有波動関数を式（7），式（11）に入れることで，エネルギー密度や応力密度をスーパーセル内の一様なFFT（高速フーリエ変換）メッシュ点上のデータとして算出することができる．しかし，密度データ自体はゲージ依存性，つまり，式（8）-式（10）や式（14）の関数形の選択に依存する不定性を含んでいる．従来のこうした手法では，表面スラブのスーパーセルに対し，エネルギー密度，応力密度を原子層毎に積分したり，巨視的な平均化を行ったり20,21），点欠陥の解析について原子毎のWigner-Seitzセルで密度を積分したり53），また，スーパーセルを対称面やボロノイ多面体で表面（欠陥）領域とバルク領域に分けて各密度を各領域で積分する20,51,54）などの方策で，実効的なエネルギーや応力の局所分布（局所エネルギー，局所応力）を探ることが行われてきた．しかし，ゲージ依存性（不定性）を厳密に検討することは行われていない．
筆者らの方法は，ゲージ依存項が積分してゼロになるような局所領域にスーパーセルを分割し，各領域でエネルギー密度，応力密度を積分することで，ゲージ依存性を直接に除去する．例えば，式（8），式（14）のエネルギー密度，応力密度のkinetic項のゲージ依存項（対称形と非対称形の差：$\nabla^{2}\rho_{e}(\vec{r})$，$\nabla_{\alpha }\nabla_{\beta }\rho_{e}(\vec{r})$）が，それぞれ局所領域ΩiおよびΩjにおいて積分してゼロになる条件は，以下である．
| \begin{equation} \int_{\Omega_{i}}\nabla^{2}\rho_{e}(\vec{r})d\vec{r} = \int_{S_{i}}\nabla \rho_{e}(\vec{r})\cdot \vec{n}dS = 0 \end{equation} | (15) |
| \begin{align} \int_{\Omega_{j}}\nabla_{\alpha}\nabla_{\beta}\rho_{e}(\vec{r})d\vec{r} &= \int_{\Omega_{j}}[\nabla_{\beta}\rho_{e}(\vec{r})]_{x_{\alpha} = x_{\alpha}^{1}}^{x_{\alpha} = x_{\alpha}^{2}}dx_{\beta}dx_{\gamma} \\ &= \int_{S_{j}}\nabla_{\beta}\rho_{e}(\vec{r})\vec{e}_{\alpha}\cdot \vec{n}dS = 0 \end{align} | (16) |
上記条件を満たす分割は，まず，(111)原子層と真空領域からなるfcc-Al(111)表面スラブのスーパーセルで行われた17）．表面平行方向をx，y，表面垂直方向をzとして，領域境界を(111)原子層の間の平坦なx-y平面に取ることで，条件を満たす原子層領域が構築できる．この場合，式（15），式（16）の面積分の被積分関数自体の値がゼロになるのでなく，面積分がトータルで差し引きゼロになる条件を考える．
まず，式（15）の最終形について，$\vec{n}$はz軸に平行なので，被積分関数の勾配はz成分のみで，以下のようになる．
| \begin{equation} \int_{S}\frac{\partial \rho_{e}(\vec{r})}{\partial z}dxdy = \frac{\partial}{\partial z}\int_{S}\rho_{e}(\vec{r})dxdy = \frac{\partial \overline{\rho_{e}}(z)}{\partial z} = 0 \end{equation} | (17) |
この(111)原子層領域は，式（16）の応力の条件について，(111)面の平行方向と垂直方向の応力成分$\tau_{xx} + \tau_{yy}$，$\tau_{zz}$について各々式（16）を満たす．$\tau_{zz}$についての式（16）は式（17）と同じになり，x-y面の境界とx-y面に垂直な境界ともに条件を満たす．$\tau_{xx} + \tau_{yy}$についての式（16）は，x-y面の境界で内積$\vec{e}_{\alpha }\cdot \vec{n}$がゼロで，x-y面に垂直な境界は，$\vec{n}$と$ - \vec{n}$の存在で，式（16）がゼロとなる．
この分割法に基づく，Al(111)表面の原子層毎の表面応力解析17）については後述する（3.3節，4.1節）．
2.3.3 Bader分割Bader分割（Bader partitioning）55）は，Trinkleらが類似した局所エネルギー計算で最初に採用したものである37）．これは，式（15）の面積分において，被積分関数$\nabla \rho_{e}(\vec{r})\cdot \vec{n}$が直接にゼロになるように，スーパーセル内部の$\rho_{e}(\vec{r})$の勾配を分析し，勾配の垂直方向成分ゼロの曲面（領域境界）Siを探索する．通常の金属などでは，原子間に（原子の周りに）垂直勾配がゼロとなる境界が存在し，原子領域（Bader領域）に分割できる．YuとTrinkle56）は，FFTの均一平行メッシュで与えられた$\rho_{e}(\vec{r})$データに対し，高精度にBader分割する効率的アルゴリズムを開発している．
一方，式（16）の対角和は式（15）と等価であるため，Bader領域は応力テンソルの対角和$\tau_{xx} + \tau_{yy} + \tau_{zz}$，すなわち原子領域の静水圧にも適用可能である．したがって，Bader領域は，ゲージに依存しない原子エネルギーに加え，ゲージに依存しない原子の局所静水圧にも利用できる18）．また，式（5），式（6）のように，原子エネルギーの和や原子応力の体積平均によって，スーパーセル内の適当な原子集団のエネルギーや平均応力を求めることができる．これは後述するように異なる原子種のBader領域間の電子移動に伴う曖昧さを取り除くのに有効である（4.5節）．
Trinkleら37）は，さらに静電エネルギー密度におけるゲージ依存の問題を扱っている．静電エネルギー密度のゲージ依存項（式（9）のMaxwell形式と式（10）のCoulomb形式の差）は2.2.1項で述べたように$\frac{1}{8\pi }\nabla \cdot (V_{H}^{\textit{tot}}(\vec{r})\nabla V_{H}^{\textit{tot}}(\vec{r}))$である．この項の体積分がゼロになるような局所領域の条件は，Gaussの定理から領域境界について$\int_{S_{i}}V_{H}^{\textit{tot}}(\vec{r})\nabla V_{H}^{\textit{tot}}(\vec{r})\cdot \vec{n}dS = 0$である．$V_{H}^{\textit{tot}}(\vec{r})$の勾配の法線方向成分がゼロとなる境界ならば条件を満たす．これは，ちょうど，$\rho_{e}(\vec{r})$のBader領域分割（式（15））と同じやり方で見つけられるはずで，類似した原子領域になる（Bader領域とは少し異なる）．Trinkleらは，エネルギー密度の原子領域毎の積分について，静電エネルギー密度の積分を$V_{H}^{\textit{tot}}(\vec{r})$の勾配ゼロの原子領域で行い，他の成分の積分を通常のBader領域で行う方法を提案している．
筆者らの方法では，エネルギー密度と応力密度の静電項にMaxwell形式を採用し，電子のkinetic項のゲージ依存性を消すために決めた局所領域で全ての項を積分する方法を採用している．エネルギー密度，応力密度のkinetic項の対称形，非対称形の差の問題のみを扱い，静電項については，文献15），文献20），文献21）などで議論されているように，Maxwell形式の方がより適切であると考え，選択する．電子のkinetic項のゲージ依存性は，電子の量子的性質に関連した，より本質的な不定性であるが，静電項は有効表現の選択の問題と考えるわけである．
異なる形式のエネルギー密度がCohenらにより提案されている57）．密度汎関数理論の全エネルギーは，電子の運動エネルギー，電子-電子，電子-原子核（イオン），原子核間（イオン間）の静電相互作用エネルギー，電子間の交換相関エネルギーの和として表される．一方で，次のような表現も可能である：占有されたKSエネルギー（Kohn-Sham方程式の固有エネルギー）の和，その和に含まれる電子-電子の静電相互作用と交換相関エネルギーのダブルカウントの補正，および原子核間（イオン間）静電相互作用エネルギーである．この表現は3.4節でも説明する（式（23））．Cohenらは，この形式でのエネルギー密度表現を扱っている．占有KSエネルギーの和の密度形式$\varepsilon_{\textit{KS}}(\vec{r})$が以下である．
| \begin{equation} \varepsilon_{\textit{KS}}(\vec{r}) = \int_{- \infty}^{E_{F}}ED(\vec{r}, E)dE \end{equation} | (18) |
| \begin{equation} D(\vec{r}, E) = \sum\nolimits_{i}f_{i}| \psi_{i}(\vec{r}) |^{2}\delta (E - E_{i}) \end{equation} | (19) |
$\varepsilon_{\textit{KS}}(\vec{r})$のスーパーセル全体での実空間積分は，占有KSエネルギーの和で，バンド構造エネルギー$E_{\textit{band}} = \sum\nolimits_{i}^{\textit{occ}}f_{i}E_{i} $である．局所領域，例えば原子Iの領域（体積ΩI）で$\varepsilon_{\textit{KS}}(\vec{r})$を積分すると局所のバンド構造エネルギーが
| \begin{equation} E_{\textit{band}}^{I} = \int_{\Omega_{I}}\varepsilon_{\textit{KS}}(\vec{r})d\vec{r} = \int_{- \infty}^{E_{F}}ED_{I}(E)dE \end{equation} | (20) |
| \begin{equation} D_{I}(E) = \int_{\Omega_{I}}D(\vec{r}, E)d\vec{r} = \sum\nolimits_{i}f_{i}\int_{\Omega_{I}}|\psi_{i}(\vec{r})|^{2}d\vec{r}\delta (E - E_{i}) \end{equation} | (21) |
量子力学に従う電子と核（イオン）の集団の系の微視的応力と巨視的応力の概念は，量子力学の初期に導入された40,41）．この問題に対する理解は，NielsenとMartin15,16）によって進展し，さらに，いくつかの理論的研究20,21,58-61）が行われた．NielsenとMartinは，量子力学系における力，応力，ビリアル定理の関係を明らかにした．応力定理では，2.2.2項で説明したように，平面波基底の第一原理計算で，全エネルギーの基底状態表現にスケーリングを導入して，巨視的な応力テンソルの表式の定式化を行った．彼らは，また，原子スケールで変化する微視的応力テンソル場の一般形を導き出した．そこでは，発散（divergence）が力の密度場（force-density field）を生み，kinetic項（運動量フラックス項）とポテンシャル項（配置依存項）をもたらす．ただし，こうした応力場には，（古典応力場と同様に）任意の応力場の回転（curl）の付加について本質的な非一意性がある．ポテンシャル項については，一般的なものとしてKugler形式62）が議論されたが，比較的短距離あるいは局所的に表現できることからMaxwell形式の方が便利であるとされた15）．
NielsenとMartin以降の一連の研究58-61）でも，量子力学系の微視的な応力テンソル場は，発散によって正しい力の密度場を生成するように，あるいは運動量バランスを保持するように，それぞれ異なる方法で導出されたが，kinetic項とポテンシャル項からなる応力テンソル場の形の特徴は，基本的に同様である．さらにどの場合も「一意性」についての問題を有している．ポテンシャル項については，どの場合も一般的なKugler形式よりもMaxwell形式が好まれている．
筆者らの方法は，FilippettiとFiorentini21）による応力密度に基づくもので，ChettyとMartin20）によるエネルギー密度と応力密度の両方の議論を受けて，平面波基底の第一原理計算法（PAW-GGA法）の枠組みで微視的な応力場を定式化したものである．kinetic項については，対称形と非対称形の非一意性があり，静電項についてはMaxwell形を採用している．上記の一連の量子力学的な微視的応力場と同類のものと見なせる．交換相関項も類似している15,58-61）．
3.3 LCAO法による局所エネルギーと局所応力密度汎関数理論に基づいて波動関数を原子軌道様の局在基底の線形結合で表現する第一原理LCAO法では，全エネルギーを各原子や各ボンドの寄与に分割するツールとして，エネルギー密度解析（energy-density analysis; EDA）が提案されている63,64）．この方法では，電子密度分布に対して計算される交換相関エネルギーはボロノイ分割により空間的に各原子に分配し65），その他の項はLCAO表現のMulliken流の原子軌道成分分解により原子や結合に割り振る．
Feilbelmanは第一原理LCAO法による応力テンソルの計算手法を開発した66）．これは，3.2節や2.2.2項で触れた平面波基底の第一原理計算におけるNielsen-Martinの方法15,16）に対応する，LCAO法における定式化である．全エネルギーの歪テンソル微分をLCAO表現で明示したもので，原子軌道基底の不完全性に対するPulay補正も厳密に導入されている．一方，上記の全エネルギーの原子軌道成分分解による局所エネルギーと同様に，交換相関項のボロノイ分割による空間分割に加え，応力テンソルのLCAO表現での原子軌道成分分解により局所応力を求めている66,67）．
Fig. 1は，Al(111)表面について，このLCAO法による原子層応力67）と，筆者らの2.3.2項の各原子層領域での応力密度積分による結果17）との比較である．LDAとGGAの違い，基底関数の違いにもかかわらず，かなりの一致が得られている．両結果とも，最上層と2層目に大きな面内張力と圧縮が存在する．第3層，第4層については，筆者らの結果は小さな面内張力，LCAO法の結果は小さな面内圧縮であるが，面内応力の層毎の振動の様子は類似している．各原子層の面内応力振動は，面内ボンド電荷のFriedel型の振動と相関している17）．4.1節で後述するように，Alのsp電子の顕著な再分布に由来する．
Ab initio local stress of an Al (111) surface modeled by a nine-layer slab17). Surface-parallel components of layer-resolved stresses, calculated by our scheme for both relaxed and unrelaxed surface configurations, are plotted, compared with Feibelman’s results using the first-principles LCAO method67). Positive and negative values mean tensile and compressive stresses in units of energy without division by a local volume in eq. (4), because the local volume is not defined for the surface top layer. For thicker slabs, similar stress oscillations at the surface layers were observed17).
固体中の格子欠陥や表面による局所的なエネルギー変化は，多重散乱（KKR）法やTB法を用いたGreen関数法によるLDOS計算により，局所的な電子構造変化から見積もられてきた42-45,50,68-70）．筆者らの局所エネルギーと局所応力の定式化とGreen関数法による定式化の比較を考える．ここでは，簡略化のため，価電子とイオンの系を局所擬ポテンシャルとLDAで扱うスーパーセルを考える．密度汎関数理論での全エネルギーは以下のように表される．
| \begin{align} E_{\textit{tot}} &= \sum_{i}^{\textit{occ}}f_{i}\int \psi_{i}^{*}(\vec{r})\left(-\frac{1}{2}\nabla^{2}\right)\psi_{i}(\vec{r})d\vec{r} + \int V_{\textit{ion}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad+ \frac{1}{2}\int V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + \int \varepsilon_{\textit{xc}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + E_{I - I} \end{align} | (22) |
| \begin{align} E_{\textit{tot}} &= \sum_{i}^{\textit{occ}}f_{i}E_{i} - \frac{1}{2}\int V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad+ \int(\varepsilon_{\textit{xc}}(\vec{r}) - \mu_{\textit{xc}}(\vec{r}))\rho_{e}(\vec{r})d\vec{r} + E_{I - I} \\ & = \int_{- \infty}^{E_{F}}ED(E)dE - \frac{1}{2}\int V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad+ \int(\varepsilon_{\textit{xc}}(\vec{r}) - \mu_{\textit{xc}}(\vec{r}))\rho_{e}(\vec{r})d\vec{r} + E_{I - I} \end{align} | (23) |
全エネルギーを各原子領域の寄与EIに分解することを考える．$ E_{\textit{tot}} = \sum\nolimits_{I}E_{I} $である．式（22），式（23）の原子領域への分解は，以下の式（24），式（25）になる．
| \begin{align} E_{I} &= \sum_{i}^{\textit{occ}}f_{i}\int_{\Omega_{I}}\psi_{i}^{*}(\vec{r})\left(- \frac{1}{2}\nabla^{2}\right)\psi_{i}(\vec{r})d\vec{r} + \int_{\Omega_{I}}V_{\textit{ion}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad + \frac{1}{2}\int_{\Omega_{I}}V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + \int_{\Omega_{I}}\varepsilon_{\textit{xc}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + e_{I} \end{align} | (24) |
| \begin{align} E_{I} &= \int_{- \infty}^{E_{F}}ED_{I}(E)dE - \frac{1}{2}\int_{\Omega_{I}}V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad+ \int_{\Omega_{I}}(\varepsilon_{\textit{xc}}(\vec{r}) - \mu_{\textit{xc}}(\vec{r}))\rho_{e}(\vec{r})d\vec{r} + e_{I} \end{align} | (25) |
式（25）の表現は，多重散乱法やTB法に基づくGreen関数法での局所エネルギーや全エネルギーの扱いで使用される42-45,50,68-70）．式（24）と式（25）の間の関係は，筆者らの局所エネルギー法とGreen関数法による局所エネルギー法との関係と言える．Crampinら45）は，式（25）を用いてfcc金属の積層欠陥や双晶の各(111)原子層毎の局所エネルギーを，バルク結晶の全エネルギーとの差として摂動論で計算した．式（25）の第1項について，原子層毎に高精度の多重散乱法（layer KKR法）によりKSエネルギーとLDOSを計算し，バルク結晶のものを差し引く．式（25）の残りの項については，バルク結晶の対応する項との差が無視できる．このことは，稠密構造の金属の積層欠陥や双晶では，原子体積や電子密度分布の短範囲秩序がバルク結晶のものと変わりないことから，force theoremにより証明される71,72）．こうして，事実上，式（25）の第1項を計算して完全結晶のものとの差を求めれば，局所エネルギーが与えられる．様々なfcc金属について，原子層毎の（バルク結晶比の）エネルギーは，積層欠陥・双晶のhcp的な積層部分で上昇し，金属種に応じて特徴的なプロファイルを持つことが示された45）．筆者らの方法でも，積層欠陥の原子層の局所エネルギーで，同様のプロファイルが得られている．
上記の多重散乱法を用いたGreen関数法で式（25）を計算するアプローチ45）では，（積層欠陥・双晶の短範囲の原子配列や電子密度がバルク結晶と同じなので）第1項以外の項のバルクとの差を無視した．このような摂動論的な取り扱いは，一般的な欠陥に対しては正当化できない．厳密には式（25）の各項を計算する必要がある．最近，Egamiらのグループ50）が同様の多重散乱法73）で式（25）の各項を原子毎に厳密に求めることを行い，原子領域の局所エネルギーを数値的に求めている．さらに歪テンソルによる実際的な有限歪を導入して原子領域のエネルギー変化を数値計算することで，原子応力を計算する方法を開発している．原子領域は，ボロノイ分割で幾何学的に決めている．
このEgamiらの原子応力計算法は，複数の遷移金属元素からなるfccの高濃度多元合金や高エントロピー合金（HEA）の原子静水圧（圧力）に適用された74-76）．原子応力と電子移動に明らかな相関が観測され，電子を獲得する（失う）原子は圧縮（引張）応力を示している．合金中の原子毎の原子応力変動と合金の降伏強度の関係が議論されている76）．一方，Ishibashiら77）は，筆者らと同様の局所応力計算法（Bader分割による原子領域）を高融点bcc金属からなるHEAに適用し，原子応力と電子移動の間に同様の相関関係を観察している．彼らは，異なる原子種間のリアルな電子密度分布の影響を取り入れてchemicalな特徴を探るには，Bader分割の方がボロノイ分割よりも優れる（首尾一貫性がある）ことを指摘している．Bader分割による各構成原子の電荷，原子体積，原子応力は，明確な相関を示し，それは各原子の価電子濃度（valence electron concentration; VEC）とその第一および第二近傍シェルの様子から説明できる77）．分割法の問題については，第5章であらためて議論する．
3.5 強結合近似（TB）法やワニエ（Wannier）関数による局所エネルギー式（25）に基づき，第1項の局所的なバンド構造エネルギー$E_{\textit{band}}^{I}$をGreen関数法で計算する試みは，TB法を用いて早くから遷移金属や半導体の表面・欠陥で実行されてきた42-44,68-70）．式（23）の第1項のバンド構造エネルギーEbandは，次式のように式（25）の第1項の各原子サイトの寄与の総和で表現できる（on-site TB表現）70）．
| \begin{align} E_{\textit{band}} &= \sum\nolimits_{i}^{\textit{occ}}f_{i} \langle \psi_{i}| H |\psi_{i} \rangle =\sum\nolimits_{i}^{\textit{occ}}f_{i} \langle \psi_{i}|\sum\nolimits_{IL}| \chi_{IL} \rangle \langle \chi_{IL} |H|\psi_{i} \rangle \\ & = \sum\nolimits_{i}^{\textit{occ}}f_{i}E_{i}\sum\nolimits_{IL}| \langle \chi_{IL} |\psi_{i} \rangle |^{2}\\ &=\sum\nolimits_{IL}\int_{- \infty}^{E_{F}}E\sum\nolimits_{i}f_{i}| \langle \chi_{IL} |\psi_{i} \rangle |^{2}\delta (E - E_{i})dE \\ & = \sum\nolimits_{IL}\int_{- \infty}^{E_{F}}ED_{IL}(E)dE = \sum\nolimits_{I}\int_{- \infty}^{E_{F}}ED_{I}(E)dE \end{align} | (26) |
TB法や関連する手法では，LDOSはGreen関数法におけるTBモーメント法やリカージョン法によって効率的に求められ42-44,68-70），式（25）の第1項以外の項は，通常原子間ポテンシャルなどで近似される．こうしたTB法は，fcc金属の積層欠陥・双晶の各原子層のエネルギープロファイルの計算に適用され，3.4節で紹介した第一原理計算結果45）と合致する結果を得ている69）．
一方，式（26）にさらに$\sum\nolimits_{I^{\prime}L^{\prime}}| \chi_{I^{\prime}L^{\prime}} \rangle \langle \chi_{I^{\prime}L^{\prime}} | = 1$を入れると，Ebandが原子間のIL-I′L′の軌道間相互作用の寄与の総和で表現できる（inter-site TB表現）70）．
| \begin{align} E_{\textit{band}} &= \sum_{i}^{\textit{occ}}f_{i} \langle \psi_{i}| H |\psi_{i} \rangle \\ & = \sum\nolimits_{i}^{\textit{occ}}f_{i} \langle \psi_{i}|\sum\nolimits_{IL}| \chi_{IL} \rangle \langle \chi_{IL} |H\sum\nolimits_{I^{\prime}L^{\prime}}| \chi_{I^{\prime}L^{\prime}} \rangle \langle \chi_{I^{\prime}L^{\prime}}|\psi_{i} \rangle \\ & = \sum_{ILI^{\prime}L^{\prime}}H_{ILI^{\prime}L^{\prime}}\int_{- \infty}^{E_{F}}\sum_{i}f_{i} \langle \chi_{IL}| \psi_{i} \rangle^{*} \langle \chi_{I^{\prime}L^{\prime}} |\psi_{i} \rangle \delta (E - E_{i})dE \\ & = \sum\nolimits_{ILI^{\prime}L^{\prime}}H_{ILI^{\prime}L^{\prime}}\int_{- \infty}^{E_{F}}D_{ILI^{\prime}L^{\prime}}(E)dE \end{align} | (27) |
半導体や絶縁体の場合，ワニエ（Wannier）関数の観点から別の局所エネルギー表現が可能である．半導体や絶縁体では，占有された（通常の広がった）固有関数（ブロッホ関数）{ψi}のセットから，ユニタリ変換により原子やボンドに局在するワニエ関数のセット{ϕl}が組み立てられる80,81）．式（23）の第1項のバンド構造エネルギーはワニエ関数を用いた表現に変換できる．
| \begin{align} E_{\textit{band}} &= \sum_{i}^{\textit{occ}}E_{i} = \sum_{i}^{\textit{occ}} \langle \psi_{i}| H |\psi_{i} \rangle \\ & = \sum\nolimits_{l} \langle \phi_{l}|H|\phi_{l} \rangle = \sum\nolimits_{l}E_{l} \end{align} | (28) |
式（28）と式（23），式（25）を組み合わせると，ワニエ関数の立場からの局所エネルギーが考えられる．ワニエ関数ϕlに局在する電子に対し，式（25）の第1項がElで，第2項以降は適当に近似するわけである．四面体半導体の結合軌道モデル（bond orbital model; BOM）82）は，この観点で理解できる．半導体の粒界エネルギー83）は，結合軌道モデルにより，歪の効果も入った各原子間ボンドの結合軌道ϕbのエネルギー期待値の和$\sum\nolimits_{b}\langle \phi_{b} | H| \phi_{b} \rangle $で見積もることができる．これは，スーパーセルの固有状態計算で与えられる粒界の系のEband（式（23）の第1項）が，占有された系の固有関数{ψi}によるハミルトニアンのtraceであり，一方，それは，{ψi}を近似的にユニタリ変換した各ボンドの結合軌道（近似的なワニエ関数）のセット{ϕb}でのtrace $\sum\nolimits_{b}\langle \phi_{b} |H| \phi_{b} \rangle $に，式（28）のように等しくなる82）．大規模計算でEbandを求めなくても，各ボンドの結合軌道の計算だけで得られる83）．
3.6 EAM（embedded-atom method）ポテンシャルによる原子応力第1章で触れたように，EgamiとVitek47-49）はEAMポテンシャルなど古典的原子間ポテンシャルを用いて，アモルファス金属中の原子レベル応力を解析している．EAMでは，固体中の原子iのエネルギーEiは，周囲の情報を入れて$E_{i} = \frac{1}{2}\sum\nolimits_{j}V(r_{ij}) + F(\bar{\rho }_{i})$で与えられ，全系のエネルギーが$E_{\textit{tot}} = \sum\nolimits_{i}E_{i} $である．jは原子iの相互作用する隣接原子である．V(r)はペアポテンシャル，$F(\bar{\rho }_{i})$は埋め込みエネルギーで，原子iの位置での仮想的な電子密度$\bar{\rho }_{i}$の関数である．$\bar{\rho }_{i}$は原子間の密度関数$\rho (r_{ij})$から$\bar{\rho }_{i} = \sum\nolimits_{j}\rho (r_{ij}) $で与えられる．VやF，$\rho (r_{ij})$などの関数，パラメータは，その金属の結晶構造や様々な特性を再現するように決められる．ここで，応力テンソルを考える．EAMで表したEtotの歪微分から，$\sigma_{\alpha \beta } = \frac{1}{\Omega }\frac{\partial E_{\textit{tot}}}{\partial \varepsilon_{\alpha \beta }}$である．2.2.2項で議論したNielsen-Martinの方法と共通である15,16,49）．式（6）を仮定し，全体の平均応力への各原子の寄与として，原子応力が
| \begin{align} \sigma_{\alpha \beta}^{i} &= \frac{1}{\varOmega_{i}}\frac{\partial E_{i}}{\partial \varepsilon_{\alpha \beta}} \\ &= \frac{1}{2\varOmega_{i}}\sum\nolimits_{j}V^{\prime}(r_{ij})\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}} + \frac{1}{\varOmega_{i}}F^{\prime}(\overline{\rho_{i}})\sum\nolimits_{j}\rho^{\prime}(r_{ij})\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}} \end{align} | (29) |
さて，文献18）では，Fig. 2とFig. 3に示すように，CuとAlの対応粒界（CSL粒界）について，式（29）によるEAMの原子応力と筆者らの第一原理の原子応力（Bader領域での積分）との比較を行い，さらに各原子領域のエネルギーについても同様の比較を行った．文献86），文献87）のCuとAlのEAMを用いた計算との比較である．原子応力（応力テンソルの対角和）は，Cu粒界では第一原理とEAMが比較的一致するが，Al粒界では，粒界の中心部分の原子でEAMと第一原理との差が顕著になる（Fig. 2）．原子エネルギーについても，Al粒界よりCu粒界の方が一致が良いが，両粒界とも中心領域では顕著な違いが見られる（Fig. 3）．Al粒界は，原子の応力とエネルギーの両方で，EAMと第一原理との一致がCu粒界に比べて良くない．このことは，Alのsp電子の性質に起因する．粒界の中心部分の構造乱れ（原子間距離，配位数，方位など）に対応して，Alのsp電子は（Cuのsd電子に比して）顕著な再分布をする．その効果がEAMでは適切に表現できないと言える．しかし，粒界領域全体の原子のエネルギーの和は，EAMと第一原理で差は大きくない．一方，Cu粒界，Al粒界ともに，一般に原子応力の方が原子エネルギーよりEAMと第一原理との違いが小さい．これは，原子応力は原子体積で割った「示強性（intensive）の量」（温度，圧力他）であるが，原子エネルギーは体積や分量に比例する「示量性（extensive）の量」（体積他）であることが関係する．後者の量は，領域体積や境界の位置，電子移動に強く影響されるためである．
Ab initio local stresses in the {221}Σ9 tilt GBs in (a) Cu and (b) Al18). Both the Cu and Al GBs have similar atomic configurations. Diagonal sums of local-stress tensors for respective atoms in the 42-atom supercell are plotted in units of GPa, compared with those by EAM potentials86,87). Atoms from No. 8 (29) to No. 15 (36) belong to the GB-core region. Atoms of Nos. 10 (31) and 13 (34) are looser sites, and those of Nos. 11 (32) and 12 (33) are tighter sites. The numbers in the parentheses indicate atoms in the other symmetric interface.
金属の表面応力は，表面再構成，表面成長，メゾスコピック構造の形成など，様々な現象を支配している．金属の表面応力の起源には，①電子の再分布，②軌道の再混成という2つのモデルがある46）．この2つの要因は完全には分離しにくく，また各々の金属の結合の性質に依存する．3.3節で説明したAl(111)表面の原子層毎の面内応力$\tau_{xx} + \tau_{yy}$の第一原理計算結果17）（Fig. 1）では，表面近く数層にわたって，面内応力が振動している．この応力振動は，面内の原子間結合の電子密度値の層毎の増減の振動と相関している．面内の原子間電子密度の増加（減少）が面内引張（圧縮）応力を発生させる．特に最上原子層では，Al原子の配位数が減り，上の原子層とのボンドを担う電子が面内のAl-Alボンドに溜まり，（面内原子間距離は完全結晶と同じなので）強い引張応力が生まれている．Alのsp電子は，表面や欠陥など構造乱れや配位数乱れに敏感で，顕著に再分布する．fcc遷移金属の(111)表面では，こうした電子密度の振動による表面応力の振動は観測されていない22,23）．
一方，一連の4dと5dの遷移金属の表面応力22,23）については，筆者らの方法による第一原理による表面応力は，dバンドのTB法の2次モーメント近似（Friedelモデル）44）を式（25），式（26）に適用して求めた表面原子エネルギーの微分からの表面応力と近いものであった．後者のFriedelモデルは，表面原子の配位数減少で，原子間d-d軌道混成で形成されるdバンドの幅が減り，高エネルギーになる描像に基づく．筆者らの第一原理による表面層の応力とFriedelモデルに基づく応力（Friedel応力と呼ぶ）が一致することは，実際に4d，5d遷移金属の表面応力の起源が，電子密度の再分布よりも表面原子の配位数減によるd-d軌道混成の変化であることを示す．
4.2 金属やSiの結晶粒界における局所エネルギーと局所応力CuやAlの対応粒界（Fig. 2，Fig. 3，{122}Σ9粒界）18），bcc Feの対応粒界24,25）に筆者らの局所エネルギー，局所応力法を適用し，各原子（Bader領域）のエネルギー，応力（応力テンソルの対角和）を求めた．結晶粒界は，バルクと異なる原子環やアモルファス的なクラスターが完全結晶に挟まれた構造で，中心部の原子は顕著な構造乱れ（配位数，ボンド長，原子体積，方位他の乱れ）を持つ．図に示されているように，どの金属でも，そういう中心部の原子は原子エネルギーと原子応力が完全結晶と大きく異なるが，界面から数原子層離れるとエネルギーは完全結晶のものに近くなり，原子応力も消える．エネルギーや応力の変化の広がりから粒界層の厚みが判定できる．スーパーセル内のバルク領域のエネルギーが結晶のものと同じで，応力も消えていれば，スーパーセルのサイズ（粒界の垂直方向の結晶部分の厚み）が十分と言える．式（5）のように，粒界の各原子の局所エネルギーの総和から完全結晶のエネルギーを引いたものが粒界エネルギーである．また，式（6）より，粒界の原子応力の体積平均が粒界のスーパーセルの平均応力で，通常はほぼゼロとなる．
一方，どの金属種でも粒界には，原子体積の小さい，原子間の接近したtighter siteと，原子体積の大きい，原子間の伸びたlooser siteの両方のsite（原子位置）が存在する．多くの場合，tighter（looser）siteで圧縮（引張）応力が観察される．bcc Feの粒界では，Stonerモデルで説明されるように，tighter（looser）siteでスピン偏極が減少（増大）する．このような対照的なsiteは，溶質や不純物の偏析挙動に大きな影響を与える26,27,30,31）．
Fig. 3（b）のAl粒界の原子（Bader領域）のエネルギーについて，原子間距離の短いtighter siteには，電子が集中し，バルクに比べて負のエネルギーになっている18）．Alのsp電子の典型的振る舞いとして，表面や格子欠陥の配位数の減った原子間に集まり，原子間距離の短い，強い共有結合的な結合を作る88-90）．粒界のtighter siteで配位数が減って短いボンドになり，電子が集中してバルクよりも低いエネルギーとなっていると言える．これは，3.3節と4.1節で触れたAlの特異な表面応力とも関連している．Al粒界で，粒界原子の周囲の環境により，バルクよりも低い原子エネルギーが起こりうることは，4.6節で説明する機械学習を用いた予測技術でも扱われている35）．
一方，Fig. 4は，実験で頻繁に観察されるSiの2つの<110>傾角粒界における原子応力の計算結果を示す28,29）．筆者らの手法では，原子毎のエネルギーや応力（静水圧，対角和）を，Bader領域でのエネルギー密度，応力密度の積分で求める．Bader領域（Bader分割）を決める計算技術（アルゴリズム）は，Siなど四面体半導体では，原子サイトから離れた場所に電子密度の瘤があるため，うまく働かないことがある．そこで，Si粒界では，ボロノイ分割65）を採用して，原子応力を求めた．3.6節で議論したように局所応力は領域体積で割った示強性の量で，領域体積や領域境界の細部に大きくは依存しないため，ボロノイ分割でも結果の信頼性は高いはずである．Fig. 4の各粒界のSi原子は全て四配位のボンドを保っているが，回転角の大きい{114}Σ9粒界では，<110>軸に沿った長い距離の再構成ボンドの導入が避けられず83），粒界の中心部のボンドに大きな引張応力が生じている（Fig. 4（a））．一方，回転角の小さい{221}Σ9粒界では，<110>軸方向の再構成ボンドなしに四配位の安定構造が形成できるので，引張応力はそれほど大きくない（Fig. 4（b））．このことは，アトムプローブ・トモグラフィーと透過型電子顕微鏡を組み合わせた実験観察で，{221}Σ9粒界に比べ，{114}Σ9粒界で酸素原子の偏析密度が高いこと28,29）を説明する．
Ab initio local stresses in the (a) {114}Σ9 (θ = 141.06 deg.) and (b) {221}Σ9 (θ = 38.94 deg.) tilt GBs in Si28,29). Stress of each bond, estimated by the average of atomic hydrostatic stresses of both-end atoms, is plotted on each atomic configuration obtained by ab initio calculations, viewed along the [110] axis. Only substantial tensile stresses are shown in units of GPa. The atomic configurations are in good agreement with electron microscopy observations28,29).
金属材料では，添加した溶質原子が粒界に高濃度で集まる「粒界偏析」がしばしば発生し，粒内の溶質濃度の低下，粒界脆化，粒界での異相析出などを引き起こす．そのため，各溶質元素の粒界偏析エネルギーEsegとその物理的起源を明らかにすることが重要である91）．置換型溶質原子の場合，Esegは，溶質原子がバルク結晶中にあり，粒界は清浄である「偏析前の系」と，バルク結晶中にあった溶質が粒界のサイトの母相金属原子と入れ替わった「偏析系」の間の全エネルギー変化である．粒界を含むスーパーセルを組み立てて，偏析前の系では溶質をバルク結晶（置換位置）に置き，偏析系では，溶質原子と粒界サイトの母相原子を入れ替える．この2つのスーパーセルの各原子の局所エネルギーを考えると，対応する原子の間の局所エネルギー（Bader領域のエネルギー）の変化の総和がEsegである．そこで，スーパーセル内の原子を以下の4つのグループに分けて，各グループの原子で偏析前と後での局所エネルギーの変化の和をT1∼T4とすれば，EsegがT1∼T4に分解される26,27）．
T1は，偏析前の系で粒界の置換サイトにあった母相原子が，偏析で溶質と入れ替わってバルク結晶に移動する際の，その母相原子自体の局所エネルギー変化（粒界にあったときとバルクの位置を占めたときの差）である．−T1の値は，バルク状態に対する，偏析前の清浄な粒界の偏析サイトの局所エネルギー上昇を意味し，Fig. 3のデータと同様の量である．T2は，溶質原子がバルク結晶部分に置換している状態と粒界の偏析サイトに置換している状態での溶質原子自体の局所エネルギー（Bader領域のエネルギー）の変化である．T3は，粒界の偏析サイトに溶質原子が来て，母相原子と入れ替わった際に，粒界の周囲の母相原子の局所エネルギー変化の総和．溶質原子が粒界のサイトに存在している際の周囲の母相原子の局所エネルギーの和から，清浄粒界のときの偏析サイトの周囲の同じ母相原子の局所エネルギーの和を差し引いたもの．T4は，バルク結晶中に溶質原子があるときの周囲の母相原子の局所エネルギーの和について，溶質原子を母相原子に替えた時の変化の総和である．−T4の値は，逆にバルク結晶について，母相原子が溶質原子に変わった際の周囲の母相原子の局所エネルギーの変化の総和である．これらの項の詳細については，文献26, 27）に図と式で説明されている．なお，偏析前の系と偏析後の系で，同じ原子であってもBader領域の形が少し変わる場合もあるので，T1∼T4の各値自体は，不確定な要素もあるが，4項の総和は厳密に正しい．これら4項を解析することで，粒界偏析の起源について貴重な情報が得られる．
bcc Feの対応粒界（対称傾角粒界，{332}Σ11と{111}Σ3）について，一連の3d遷移金属（transition metal; TM）溶質の優先的な偏析サイトと偏析エネルギーEsegの第一原理計算を行った．さらに局所エネルギー計算から，T1∼T4の各項を求めた30）．周期表の左のearly-TM元素（Sc，Ti，V）はlooser site（原子体積大）の方に，周期表の右のlate-TM元素（Co，Ni，Cu）はtighter site（原子体積小）の方に優先的に偏析し，Feに近いmiddle-TM元素（Cr，Mn）はsite依存性があまりない．Esegの値は，周期表の両端の3d遷移金属元素（Sc，Ti，Ni，Cu）で大きく（絶対値の大きな負値），中央に近づくと小さくなる（負値が上昇してゼロに近づく）が，再び中央付近のMnで少し大きくなる（下がってカスプになる）．周期表に対してM字型のEsegの変化は，3d遷移金属溶質とbcc-Fe中らせん転位の相互作用の第一原理計算結果92）と類似している．
こうしたEsegの溶質元素依存性の起源を探るため，Fig. 5に示すように，EsegのT1∼T4への局所エネルギー分解を行い，またLDOS解析と組み合わせて分析した30）．LDOS分析では，母相のFeでminority-spinのd電子が主に周囲のFeや溶質と相互作用する点が重要である．Fig. 5から，まず，溶質全般にT1とT2は互いに相殺される傾向がわかる．early-TM（Ti，V）では，T3がlooser siteで負であることでlooser siteでのEsegの利得が生まれている．early-TMの溶質の3d電子は，母相のFeより数が少なく，looser siteでanti-ferro的な配置で（溶質の多数派スピンの3d電子が）周囲のFeのminor-spinのd電子と強く混成することで，周囲のFe原子を（置換前のFe原子があるときよりも）安定化させる（T3が負）．late-TM（Ni, Cu）では，T4の大きな負の値がEsegの利得を支配している．上記の−T4値の意味からlate-TM溶質原子はバルク結晶中で周囲のFe原子を不安定化させる．T3の正の値は，late-TM溶質が粒界のサイトにあるときにも周囲のFe原子を不安定化させることを意味する．Fig. 5から−T4 > T3なので，周囲のFe原子を「不安定化させる効果」がバルク中の方が粒界サイトより大きく，これがlate-TM溶質の偏析の起源と言える．late-TM原子はFeより多数のd電子を持ち，増加したd電子は周囲のFeのminority-spinのd電子とあまり混成せず（early-TMと逆），避けあう傾向にあることがLDOSからわかる．
Variations of the four terms, (a) T1, (b) T2, (c) T3 and (d) T4, in the local-energy decomposition of Eseg (eV/atom) for the segregation at looser and tighter sites (blue and red triangles) of each 3d-TM solute in the Σ11(332) (solid triangles) and Σ3(111) (empty triangles) GBs in bcc Fe30). The sum of the four terms is equal to Eseg. The meaning of each term is explained in the text.
middle-TM元素のCrは，Coも含めてFeに近いのでEsegが小さいが，Mnのみが例外的にEsegが負でカスプになる．その原因は，Fig. 5からT4の負の値が主因である．late-TMのNi，Cuに似た特徴がある．Fe結晶中に単独で置換して存在するMnのLDOSを見ると，anti-ferro的な電子配置で，Mn内では多数派スピンのd電子（Feとの関係ではminority-spin側）は，周囲のFeのminority-spinのd電子と混成せず，強い局在性を示している．Mn内のminority-spinの方（Feとの関係ではmajority-spin側）も，周囲のmajority-spinのFeのd電子との混成が少ない（この点はlate-TMとも異なる）．Mn原子のd電子は，フント則で説明されるように局在した高スピン状態として存在する傾向があり，そのせいで周囲のFeとのd-d混成が顕著に弱く，周囲のFe原子を不安定化する．その傾向が，結晶中に比べ粒界では少し緩和されるので，（late-TM溶質と同様に）−T4 > T3から偏析利得が生まれる．Mnでカスプを持つ現象は，3d-TM溶質と点欠陥や転位との相互作用の第一原理計算でも観察されており92），同様の起源と考えられる．
一方，Fe粒界のsp元素溶質，Al，Si，P，Sについての局所エネルギー解析も行われた31）．sp元素の場合もEsegは，T3と−T4の関係，すなわち，T3 + T4の値が支配すると言える．sp元素が，周囲のFe原子を，①粒界でバルク中と比べてどれくらい安定化させるか，あるいは②不安定化させる効果が粒界でバルクに比べて小さいか，が偏析を支配する．さらに，sp元素が周囲のFe原子の安定化（不安定化）に及ぼす効果は，「Fe-Fe間のd-d混成」に比しての「溶質-Fe間のsp-d混成」の強度が決めるので，原子間距離とsp元素原子の価電子原子軌道の空間的広がりに支配されると言える．3.5節で論じたように，文献31）では，COHP解析（式（27））で様々な偏析サイトやバルク置換位置でのAl-Fe，Si-Fe，P-Fe，S-Feの各々の相互作用を分析し，周囲のFe原子の安定化・不安定化の観点から各元素のEsegや偏析挙動（site preference）が説明できることが示された．特に原子軌道の広がりの短いPとSは，Fe結晶中のバルク置換位置では弱いsp-d混成のため近接Fe原子を不安定化するが，Fe粒界の原子間の短いtighter siteでのみ，近傍Fe原子を安定化したり，不安定化の程度が下がるので，大きな偏析利得が生まれる．
3d-TM溶質の場合もsp元素溶質の場合も，溶質原子自体の局所エネルギー変化T2よりも，粒界やバルク中の溶質原子周囲の母相Fe原子の局所エネルギー変化T3，T4の方がEsegを支配すると言える．これは，T2は溶質-Fe相互作用のバルク中と粒界サイトでの変化だが（溶質原子から見た周囲のFe配位数やボンド長の変化），後者のT3，T4は，Fe-Fe相互作用がFe-溶質相互作用に入れ替わる（Feの近接原子がFeから溶質に元素種が替わる）際のFeの局所エネルギー変化なので，後者の方が定量的に大きな効果を持つためである．
4.4 金属粒界の第一原理引張試験第一原理引張試験（first-principles tensile test; FPTT）は，原子や電子の自然な挙動に従って，粒界や異相境界，あるいは欠陥を持つ系の変形・破壊過程や本性的強度を探る強力な手段である12）．スーパーセル構造（格子と内部座標）に小さなステップで引張歪を加え，その歪下で第一原理計算による原子に働く力に従って原子緩和を行い，その構造に再び引張歪を加え，緩和する過程を繰り返し，自然に変形や破壊を発生させる．この引張変形過程での電子やボンドの挙動を探るとともに，局所エネルギー・局所応力の変化を調べることが興味深い．
筆者らは，このスキームをAlとCuの対応粒界（{122}Σ9粒界）に適用した32）．Fig. 2，Fig. 3は，FPTT前の初期緩和構造の原子応力，原子エネルギーの様子である．FPTTの過程で，これらの局所エネルギー，局所応力が変化する．Al粒界とCu粒界で変化の様子がかなり異なることが判明した．Al粒界では，tighter siteで原子間結合が強化され，looser siteでは引き伸ばされて弱化している．この不均一性が，原子エネルギーと応力の複雑な変化を引き起こす．最初，tighter siteの強化されたボンドのバックボンドが主に伸ばされ，tighter siteの局所応力とエネルギーが急激に上昇する．一方，強化されたボンドは短いが，ある程度伸びると急激に弱化し，それをきっかけに全ての界面結合が破断する．各歪段階で，界面および背後の原子の原子エネルギー，原子応力の変動は大きく，Al sp電子に特徴的な，原子配列やボンドの変化に即応して局所電子密度が再分布する現象と関係している．一方，Cu粒界では，積層欠陥エネルギーが小さいために引張変形過程でスーパーセルのバルク領域に積層欠陥が導入されたが，粒界部分の局所エネルギーと応力の変化はむしろ単純で均一的であった．
なお，通常のFPTTでは，スーパーセルの界面平行方向のサイズは大きくは取れないので，粒界に付随するクラックや転位の挙動は十分には扱えない．厳密には，クラックや転位の核生成，注入，移動の自由度が十分に扱える大きなスーパーセルが望ましい．一方，Alのねじり粒界のFPTTについての局所エネルギーの第一原理計算から，粒界の界面層の局所歪vs局所エネルギーのuniversalな曲線が求まり，メゾスコピックな破壊力学における粒界層の変形・破壊用の埋め込み関数として使えることが提案されている33）．
4.5 合金の微細構造の局所弾性定数微細構造や組織を持つ材料の弾性的あるいは力学的特性を理解するために，材料中の局所領域の弾性定数を，「局所応力の1次微分」あるいは「局所エネルギーの2次微分」を介して求めることができるはずである．今，微細構造を含む大きなスーパーセルを考える．セル全体の体積弾性率は，セル全体の等方膨張（+δ），等方圧縮（−δ）でのスーパーセルの静水圧$(\sigma_{xx} + \sigma_{yy} + \sigma_{yy})/3$の差をセル全体の圧縮，膨張の体積変化の割合で割れば求められる．この等方圧縮，等方膨張の過程で，各原子のBader領域の静水圧変化，Bader体積変化を求めれば，各原子の局所体積弾性率Biが以下のように得られる．
| \begin{equation} B_{i} = \frac{\sigma_{i}(+ \delta) - \sigma_{i}(- \delta)}{(\Omega_{i}(+ \delta) - \Omega_{i}(- \delta))/\Omega_{i}^{0}} \end{equation} | (30) |
このスキームを用いて，FeリッチFe-Si合金のSi濃度依存の体積弾性率変化の機構が明らかにされた34）．実験的には，bcc Fe合金試料の体積弾性率は，置換Si濃度の増加に対して10 at％Si程度までは減少するが，この濃度を超えると急激に回復する93）．この機構を探るため，大きなbcc Feスーパーセル内の体心サイトに置換Si原子を配置し，体心位置のSiと隣接bcc Fe原子によるSiFe8立方体クラスターがFe中に適当に配列した合金モデルを扱う．Si濃度に応じてSiFe8クラスターの分布と割合が変化する．Si同士は隣接しないと考えられるので，濃度が増えるとSiFe8クラスター同士がFeを共有して接するようになる．
Fig. 6は，各Si濃度のスーパーセル内のSiFe8クラスター部分と残りの純Fe領域のそれぞれの局所体積弾性率，およびスーパーセル全体の平均の体積弾性率である．全体の体積弾性率は11.11 at％Siまで減少し，このSi濃度を超えると回復する．このSi濃度依存性は実験と合致する93）．各領域について，SiFe8クラスター領域の平均は6.25-11.11 at％Siまで非常に軟らかくなり，それを超えると回復する．残りの純Fe領域の変動は小さい．6.25 at％SiまではSiFe8クラスターが別々に存在しているが，6.25-11.11 at％SiまでのSi濃度では，SiFe8立方クラスターの対角のFe原子を共有するSiFe8クラスターの接触，1次元連鎖が不可避的に発生する．こういう系では，2つのSiFe8クラスターに共有される，角のFe原子部分で局所体積弾性率が顕著に減少している．一方，11.11 at％Siを超える濃度では，SiFe8立方クラスターのエッジが共有され，SiFe8クラスターが3次元あるいは2次元に「より密着して」配列するようになる（そうしないと高濃度Siを収容できない）．そうすると，Fig. 6に示すように，クラスター部分はSi-Feのsp-d混成が強くなり，高い体積弾性率を示すようになる．
Calculated local bulk moduli of (a) pure-Fe and (b) SiFe8-cluster regions as well as (c) the cell-averaged ones for the FeSi alloy supercell models with various Si concentrations are plotted in units of GPa34). At 12.5 at％Si, partial D03 models, I and II, contain two-dimensional and three-dimensional arrangements of edge-shared SiFe8 clusters, respectively. A number on each vertical bar indicates the volume ratio of each region in the FeSi alloy supercell.
このようにSi濃度に依存してSiFe8クラスターの配列や接続の仕方が変化し，それがSiFe8クラスター領域の結合状態の変化，局所弾性率の変化を生み，全体の体積弾性率のSi濃度依存性が生じる．このことは，平均化された電子構造93）を調べるのではなく，筆者らのスキームを用いてそれぞれの領域の局所的な体積弾性率を調べることによってのみ明らかにできる．異なる原子種間のBader分割は，領域の境界を越えた電子移動に関連して，原子個々の局所量の解釈に複雑な問題をもたらすが，SiFe8クラスター領域で和や平均を取ることで，そういう問題が回避される．
4.6 機械学習を通じたモデリング最近の材料科学・工学では，機械学習の応用が注目されている94）．筆者らの局所エネルギー・局所応力の第一原理計算と機械学習の組合せは，この分野で新しいフロンティアを開くことが期待できる．機械学習の目的の1つは，実験や第一原理計算のデータから，新しい予測や洞察を効率的に導き出すことである．局所エネルギーや局所応力を含むようにデータを拡張することで，機械学習に基づく予測の対象や範囲を拡大することが期待できる．表面や粒界，力学的挙動，溶質の挙動，触媒など，局所的な結合が重要な系において有用なはずである．
Tamuraらは，Alの様々な対応粒界の局所構造と局所エネルギーの関係を機械学習で処理することで，粒界エネルギーを予測する技術を構築した35）．まず，比較的小さいΣ値の対応粒界を選択し，筆者らの局所エネルギーの第一原理計算法により，粒界の各原子（各Bader領域）の局所エネルギーを求め，各原子の局所構造の幾何学的データとともに蓄積する．このtraining dataに線形回帰モデルLASSO（least absolute shrinkage and selection operator）95）を適用し，局所環境と原子エネルギーの関係を探索する．局所環境データは，SOAP（smoothed overlap of atomic positions）96）など，各原子の局所的幾何学的特徴の記述子を様々に検討し，予測精度を最適化したLASSOモデルを構築した．これらを新たな対応粒界構造に適用し，各原子のエネルギーを予測し，それらの和で粒界エネルギーを求め，精度を検証した．粒界エネルギー自体もデータに組み込むことで，さらに予測精度が向上している．
Tamuraと共同研究者は，このアプローチをbcc Feの一般粒界，すなわち，対応粒界のような規則構造を持たない，大規模な乱れた構造の粒界に拡張した36）．EAMポテンシャルなどを用いて溶融体から凝固させた大規模多結晶構造での一般粒界のエネルギーの高精度計算を考える．上述のAl粒界でのアプローチでは，スーパーセルで扱える対応粒界の第一原理計算から機械学習のデータを構築した．一般粒界の場合，構造が大規模すぎてスーパーセルによる第一原理計算が実行できない．そこで，一般粒界の構造から適当にピックアップした局所の原子について，それを中心に近傍構造をボックスに切り出し，人為的なスーパーセルに並べて第一原理計算を実行し，中心原子の局所エネルギーを求める．この過程を一般粒界の様々な局所原子で繰り返すことで，上述のAl粒界の場合と同様にbcc Feの一般粒界の粒界原子の局所環境vs局所エネルギー（原子エネルギー）のtraining dataが蓄積でき，Al粒界の場合と同様に線形回帰モデルが構築できる．
比較的小さいスーパーセルで局所構造を切り出し，中心原子の局所エネルギーを求めることがキーポイントである．bcc Feでは，金属電子による遮蔽と，構造乱れによるd-d混成変化の及ぶ範囲が比較的短いため24,25），比較的短い範囲の局所構造の切り出しで，中心原子の電子構造や局所エネルギーが元の大規模構造中のものと同じように再現できるはずである．効率的に乱れた構造の原子をピックアップし，データを蓄積していけば，予測精度が向上していく．なお，局所環境と原子エネルギーの相関を見出すことは，原子間ポテンシャルの開発と同じ意味を持つ．機械学習による原子間ポテンシャルに関しては，Trinkleらの第一原理局所エネルギー法37）の改良版97）がアモルファスSiのニューラルネットワークに基づくポテンシャル構築に使われている98）．
局所応力についての機械学習の応用では，3.4節で触れたように，高融点金属による高エントロピー合金（HEA）の原子応力と局所化学環境の相関が扱われている77）．各原子とその近傍シェルの価電子濃度（VEC）の様子が原子応力を支配するようである．今後，局所エネルギーや局所応力と機械学習の組合せは，さらに様々な現象への適用が期待される．
筆者らの局所エネルギー・局所応力の第一原理計算法は，従来の平面波基底の第一原理計算では未利用であった，スーパーセル内のエネルギーや応力の分布情報を取り出し，活用する手法である．ただ，エネルギー密度や応力密度の積分のためにスーパーセルを局所領域に分割する仕方は，普遍的なものは決めにくく，局所エネルギーや局所応力の値自体は分割の仕方に依存する．しかし，分割された局所エネルギーの総和と分割された局所応力の体積平均は，従来法で得られるスーパーセルの全エネルギーとセル応力に厳密に等しい．局所エネルギー・局所応力の分布は，スーパーセル内の局所的な原子・電子の挙動を反映したもので，複雑な構造や現象の本質や機構を理解するための有用な情報である．
Bader分割は，上で説明した応用例やTrinkleらのグループの応用例37,99,100）も含めて示されているように，エネルギーや応力の局所分布を統一的に探るための自然かつ実用的な分割法である．もちろん，局所エネルギーと局所応力には，3.6節で論じたように示量性（extensive）と示強性（intensive）の問題があり，前者（示量性の局所エネルギー）は領域の境界や大きさに本質的に依存し，解釈に困難が生じる場合もある．一方，2.3.3項や4.5節で議論したように，原子領域の量は，和や体積平均を取ることで，クラスターや原子集団の領域の量に拡張でき，そこでは原子毎の領域境界の問題は低減化される．LDOSやCOHPなどの電子的挙動に関する補足的解析を組み合わせることも有効である30,31）．
第3章で議論したように，類似した局所のエネルギーや応力を扱う理論的手法として，エネルギー密度の別表現，LCAO法，多重散乱法やTB法を用いたGreen関数法，ワニエ関数法，あるいはEAMポテンシャルを用いる方法などがある．筆者らの手法とこれらの手法との比較は，まだ限られているが，3.3節，3.4節，3.6節で説明したように，金属表面，金属中の積層欠陥，合金，結晶粒界の計算結果に関して，基本的に矛盾せず，違いがある場合もリーズナブルなものであった．従来法との比較は，局所エネルギー，局所応力の概念や物理的意味を明確にするために重要である．もちろん，局所エネルギー，局所応力，局所領域の定義やゲージ依存性問題の扱いには自由度があり，各方式の優劣を単純に評価することはできない．しかし，筆者らの方法は，定式化における曖昧さや人為的操作を極力減らしており，計算手法の簡便性・柔軟性・精度，従来の平面波基底の第一原理計算と同じ電子構造を共有できる点など，多くの優位性がある．筆者らの手法は，LCAOやGreen関数法などの局所の電子構造を扱う手法と，広がった波動関数を平面波基底で扱う手法との橋渡しをするものである．
筆者らの手法では，スーパーセルを原子，原子層，クラスターなどの領域に，あるいは欠陥周囲とバルク領域に分割することができる．それは，価電子密度分布に従って自然にゲージ依存性を取り除くように実行できる．一方，LCAO法による局所エネルギー63,64）や局所応力66,67）は，原子軌道基底に依存し，最適な原子軌道基底を曖昧さなく定義する一般的な方法はない．多重散乱法を用いるGreen関数法では，原子球近似（atomic sphere approximation; ASA）45,50）による原子球の構築は，細密充填構造以外の系では，何らかの人為的な操作が必要になり，また，原子のエネルギーと応力は，通常ボロノイ分割を用いて計算される50,74-76）．しかし，3.4節で触れたように，HEAへの応用では，原子配置に基づくボロノイ分割に比べ，電子密度分布の法線勾配をゼロに保つ境界で分割するBader分割の方が，物理的に意味のある量を提供でき，適切であると言える77）．また，単体固体の場合，Bader分割は局所的な電荷中性に近い原子領域を提供するはずで，強結合近似法などにおける前提と合致する68）．
理論面や計算技術上で，いくつかの検討すべき問題が残っている．第一に，Bader領域で区切った場合も各領域は完全に中性ではなく，電荷移動が存在する．それらの局所エネルギー，局所応力に及ぼす影響について一般的に検討されるべきである．局所エネルギーのような示量性の量の場合，影響が大きいはずである．電荷移動が本質的に伴うと思われるイオン性固体，化合物，合金の場合，局所エネルギーと局所応力を，カチオン，アニオンの各々の値で分析するのではなく，隣接するイオンペアまたはクラスターの和または平均として分析する方が有効と言える．そうした新しい分析技術を検討することも望まれる．第二に，上記の点に関連して，局所領域を原子中心領域ではなく，結合領域など，新たな領域分割法を検討することも望まれる．4.2節でも触れたように，通常の原子中心のBader分割のアルゴリズムは，電子密度の勾配や分布に基づく探索のため，電子密度の上昇が原子位置以外で顕著な四面体半導体の系では，必ずしも安定ではない．第三に，エネルギー密度および応力密度における静電項の表現，例えばMaxwellとCoulomb形式，MaxwellとKugler形式の違いによる影響を，2.2節，2.3節，3.2節でも議論したようにさらに検討する必要がある．第四に，通常求めている静水圧応力（応力テンソルの対角和）ではなく，各応力テンソル成分について，式（16）を満たす局所領域について検討する必要がある．第五に，筆者らの方法における局所応力と原子に働く力の関係は，3.2節で触れた量子力学的応力と力場に関する一般的な議論15,20,21,58-61）とは異なり，あまり明確ではない．また，連続体モデルや原子モデルの古典的な応力場との比較84,85,101）も興味深いがあまり検討されていない．
最後に，局所エネルギーの大規模システムへの応用の将来性を強調したい．文献36）では，bcc Feの大規模構造の全エネルギー（完全結晶比のエネルギー上昇）が，各構成原子の局所環境と原子エネルギーの関係を機械学習によってモデル化することで見積もられた．この機械学習のデータは，元の大規模構造から効率的に原子をピックアップし，その原子と周囲の構造を大規模構造から切り出し，人工的に小さなスーパーセルを構築して第一原理計算を行い，局所エネルギー法を適用することで求められた．これは，原子と周囲の構造を切り出した「小さなスーパーセル」の大きさが，スーパーセルの中心原子の電子構造や局所エネルギーを元の大規模構造中でのそれらと同じに再現できるものであればよい．金属では長距離静電相互作用のスクリーニングや平均化により，それほど大きく取らなくても可能のはずである．ここで，もし大規模構造の全ての原子について，上記のように順に周囲を切り出して小さなスーパーセルにして第一原理計算を行い，局所エネルギー法を適用すれば，大規模構造の全原子のBader領域のエネルギーがシームレスに求まるはずである．これは，分割統治（divide and conquer; DC）型のオーダーN法102,103）に相当する．また，例えば，乱れの大きい部分には小さいスーパーセルを介したDC型第一原理計算を適用し，その他の部分には機械学習ベースの予測や経験的ポテンシャルを適用することで，大規模金属系の高信頼性のマルチスケールあるいはハイブリッド計算が可能なはずである．こうしたアプローチ，すなわち局所構造を切り出した小さなスーパーセルでの局所エネルギー計算と機械学習との組合せは，オーダーN法やマルチスケール計算として大きな展開が期待できる．今のところ，与えられた構造に対するエネルギーの計算であるが，原子に働く力や応力を含めて機械学習で扱うなど，発展性がある．
本Overviewの結論は以下のようにまとめられる．筆者らの局所エネルギー・局所応力の第一原理計算法は，平面波基底PAW-GGA法の枠組みの中で，ゲージに依存しない各局所領域でのエネルギー密度，応力密度の積分により，大きなスーパーセル内部のエネルギー分布，応力分布を明らかにする．これまで，類似するエネルギー密度法，LCAO法による方法，多重散乱法やTB法をGreen関数法に組み合わせる方法，EAMポテンシャルを用いる方法などで，類似した局所の物理量を計算する試みが行われている．これらと比較することで，筆者らの方法による局所エネルギー，局所応力の物理的な意味が明らかになり，筆者らの方法の優位性，すなわち，曖昧さの低減，柔軟性と精度，従来の平面波基底の第一原理計算法とのデータ共有や互換性などが示された．なお，局所エネルギーについては，（総和は厳密には正しいが）示量性の量として局所領域の分割法への依存性は無視できず，物理的解釈には注意が必要である．状況に応じて局所領域を原子からクラスターに拡張したり，LDOSやCOHPなどの電子挙動解析を組み合わせたりすることが有効である．最近の金属表面，結晶粒界，粒界への溶質偏析，金属粒界の引張試験，合金の微細構造の機械的性質などへの応用では，局所エネルギー・局所応力計算によって現象の新しい側面が明らかになり，従来の第一原理計算では到達できなかったメカニズムの解明が可能となった．また，局所環境と局所エネルギーの関係を機械学習でモデリングすることで，金属系の粒界エネルギーや大規模な乱れた構造のエネルギーの高精度予測が可能となった．今後の幅広い応用に向けた未解決の課題についても議論した．将来展望として，特に，機械学習技術とDC型の切り出した構造の局所エネルギー・局所応力計算とを組み合わせるアプローチは，大規模金属系の高精度のマルチスケール計算技術への発展が期待される．以上から，局所エネルギー・局所応力の第一原理計算は，材料科学・工学に新たなフロンティアを拓くと言える．
本研究は，文部科学省「ポスト京プロジェクト（次世代の産業を支える新機能デバイス・高性能材料の創出；CDMSI）」，文部科学省「構造材料元素戦略研究拠点（ESISM）」，新エネルギー・産業技術総合開発機構（NEDO）「革新型蓄電池実用化促進基盤技術開発（RISING2）」（JPNP16001）の支援を受けた．本研究の一部は，JST産学共創基礎基盤研究プログラム「ハミルトニアンからの材料強度設計」，科研費新学術領域研究「バルクナノメタル」，科研費補助金（19H05177，20K05696）の支援を受けた．QMASコードと本手法の開発に貢献された石橋章司博士，粒界や大規模システムにおける機械学習との組合せを進められた田村友幸准教授，各種応用を進めていただいたHao Wang博士，Somesh Kr. Bhattacharya博士，Zhuo Xu博士に感謝いたします．また，大野裕准教授には，Si粒界についての共同研究で，伊藤一真博士（日本製鉄㈱），澤田英明博士（日本製鉄㈱），尾方成信教授，陳迎教授，毛利哲夫名誉教授には，Fe粒界およびFe合金に関する共同研究でお世話になりました．感謝致します．
