20160419 先駆的科学計算に関するフォーラム@九州大学箱崎キャンパス

大規模第一原理計算による金属ガラスのナノスケール変形機構の解明
豊田工業大学大学院工学研究科
椎原 良典

SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

1	

第一原理計算と材料の機械的性質
第一原理計算： 量子力学に基づき電子状態を記述する波動関数を求め，エネルギーや原子間力等を計算する手法
波動関数→電子密度

理想強度

セル応力
ヤング率
原子に働く力系のエネルギー

MgZn2合金の原子構造（周期境界条件）
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

ポアソン比
線膨張係数

結晶の応力ひずみ線図2	

金属ガラスのナノスケール変形機構
金属ガラス： 優れた機械的特性，一方，変形機構の理解は不十分
微視的な力学応答：
弾性体

金属ガラス

低ひずみ状態から塑性的挙動が起こる
CuZr系金属ガラス

微視的視点：電子・原子構造と力学状態の相関を明らかとしたい
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

3	原子応力による応力緩和機構の解明
荷重下の構造緩和と応力分布変化の関係を直接捉える

弾性体

原子応力解析で位置

金属ガラスを特定，機構を解明

大規模計算を通じて
CuZr系金属ガラス

構造変化の広がりを見る

目的：第一原理原子応力解析法を開発，大規模原子系へ適用
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

4	

成果

1. 金属ガラスの原子応力計算
2. 1500原子規模金属ガラスの大規模第一原理計算

SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

5	

成果

1. 金属ガラスの原子応力計算
2. 1500原子規模金属ガラスの大規模第一原理計算

SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

6	

第一原理計算：支配方程式
密度汎関数法による近似Schrödinger方程式
電子密度Ĥ[ ρ ]ψi = ε iψi
Hamiltonian

一電子波動関数 固有エネルギー• 全ての電子・スピンについて解く（方程式の数は電子数ｘスピン）• Hamiltonianは電子密度依存：非線形性• 波動関数は規格直交条件を満たす：条件付き
電子密度

ρ = ∑ fiψii2

規格直交条件*∫ ψ ψ dr = δijij

条件付き連立非線形（一般化）固有値問題を解く
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

7	

第一原理原子応力計算
第一原理計算の基本的なアウトプット：
全エネルギー，力，セル応力，etc.

セル応力
原子応力

∂Etot
Etot = Etot + ∑εαβ +!
αβ ∂εαβε

エネルギーのセルひずみ微分

第一原理原子応力計算法を提案
N atom

N atom∂Etot, iEtot = ∑ Etot, i + ∑ ∑εαβ +!
ii αβ ∂εαβε

原子あたりエネルギーのセルひずみ微分

第一原理計算ソフトOpenMXにセル応力・原子応力を実装
例：Mg合金内の介在物クラスター近傍応力分布
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

8	金属ガラスの応力緩和
対象：Zr114Cu86 (分子動力学でメルト・急冷後にセル緩和して作成)
ひずみ0の構造を一方向に引張り，構造緩和する
応力緩和

参考：VASPによるZr112Cu88（別の構造）の結果

原子構造緩和により応力緩和が発生
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

9	原子応力分布
原子応力テンソルのグリフ表示：主応力テンソルの図示

ZrCu

ひずみ0でも内部に応力分布は存在．Zrは圧縮，Cuは引張SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

10	セルひずみと原子応力分布
ひずみを印加しても原子応力分布にはほとんど変化が見られず
原子種固有の応力が大きいため原子応力の差異で議論する必要

ひずみ0.0%

ひずみ1.0%(緩和)

ひずみ2.5%(緩和)SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

ひずみ5.0%(緩和)11	

応力緩和と原子応力
x方向ひずみ下における構造緩和前後の応力変化

原子応力計算法により，応力緩和をもたらす原子を特定可能
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

12	

成果

1. 金属ガラスの原子応力計算
2. 1500原子規模金属ガラスの大規模第一原理計算

SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

13	

第一原理計算：大規模系での問題
3つのループ：波動関数の求解，原子構造緩和，セル構造緩和Ĥ[ ρ ]ψi = ε iψiρ = ∑ fiψi 2

密度汎関数法による近似 i
Schrödinger方程式

波動関数→電子密度

Fj = 0( j = 1,..., N atom )σ ij = 0セルに生じる応力が0

全原子に働く力が0

大規模原子系における問題：１．メモリ

セル応力

２．計算量が原子数の３乗
３．反復回数の増大
原子に働く力

Nオーダー法による系のエネルギー解決を模索

SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

14	

Nオーダー法とそのパラメータN：原子数• 一般的な第一原理計算はNの3乗の計算量

Ĥ[ ρ ]ψi = ε iψi
行列のサイズがNの2乗，電子軌道数がNに比例，計3乗• 行列を疎行列化できれば計算オーダーを低減できる
Nオーダー法のアイデア：局所性を利用した分割統治法• Krylov部分空間法（Ozaki, 2006）でのパラメーター
クラスターサイズ，部分空間の次元
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

調整が必要
15	

パラメーターの決定“クラスターサイズ”と”Krylov部分空間の次元”を従来計算法との誤差から決定した（200原子の計算）

1200次元117原子クラスターサイズとエネルギー精度

部分空間の次元とエネルギー精度• milliHartreeレベルの誤差• 部分空間の次元の方が精度に影響大
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

16	

大規模計算：収束の問題
電子密度から生じる非線形の反復計算（SCFループ）長波長成分が関与⇒大規模系で収束性悪化

1500 atoms

200 atoms

SCF反復間におけるエネルギー変化量

パラメータの調整により収束解を得たが，反復回数は増大
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

17	

1500原子規模金属ガラスの原子応力分布

Fig: Atomic stress tensor [eV] visualized in tensor glyph.

• 電子状態，原子応力の計算に成功• 力・セル応力が収束した結果を得るには至らす
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE

18	

成果のまとめ
第一原理原子応力計算法を開発し，CuZr系金属ガラスに適用した．また，Nオーダー法による1500金属ガラス計算を実施し，SCF収束解と原子応力分布を得た
・第一原理原子応力計算
第一原理計算ソフトOpenMXに第一原理局所応力計算を実装した．
CuZr系金属ガラスの一方向変形に際して応力緩和をもたらす原子を特定した．
今後，原子構造・電子構造変化との関係を調査し，変形機構の解明を目指す．

・大規模原子応力計算
1500原子規模金属ガラスのSCF収束解と原子応力を得たが，力とセル応力を収束させるに至らなかった．パラメータの更なる調整，高精度な分子動力学を前処理として行う等の改善を通じて所望の大規模計算を実現する．
SOLID MECHANICS LAB., TOYOTA TECHNOLOGICAL INSTITUTE19