よくわかる計算工学最前線 「原子間ポテンシャルの 基礎から応用まで 」

Presentation on theme: "よくわかる計算工学最前線 「原子間ポテンシャルの 基礎から応用まで 」"— Presentation transcript:

1 よくわかる計算工学最前線 「原子間ポテンシャルの 基礎から応用まで 」
よくわかる計算工学最前線 「原子間ポテンシャルの 基礎から応用まで 」 東京大学工学系研究科機械工学専攻　泉　聡志

2 日本材料学会分子動力学部門委員会 ポテンシャルWG 局所構造を反映した原子間ポテンシャル作成手法の調査
Tersoff型、 MEAM(modified EAM)型、 Tight-binding法 データベース化 分子動力学冬季講座 fcc, bcc系のポテンシャル HCP系のポテンシャル 共有結合系のポテンシャル

3 概要 ポテンシャルの基礎（思想） ポテンシャルの応用 二体ポテンシャル イオンポテンシャル 多体ポテンシャル(共有結合系)
クラスターポテンシャル ボンドオーダーポテンシャル Tight-Binding法 ポテンシャルの応用 シリコンクラスター用Tersoffポテンシャルの作成 SiH4表面反応用ポテンシャルの作成

4 原子間ポテンシャルの階層 Density Functional Theory ( GGA, LDA )
CASTEP, Wien97 C.Z.Wang, D. G. Pettifor … (Order-N) Tight-Binding Tersoff, Finnis… Empirical Potential　（Bond-order） Vashishita… Charge-Variable Potential Lennard-Jones, Morse Empirical 2-body Potential 近似 速度

5 密度汎関数法（DFT） 様々な近似（擬ポテンシャル・交換相関エネルギ）に基づき、電子状態をセルフコンシストに算出
原子個数　50～100個程度（状態の緩和＋エネルギ計算） 汎用ソフトが整備されつつある VASP CASTEP

6 密度汎関数法（DFT）

7 VASPによるc-Si,SiO2（β-c）電荷密度の計算
密度汎関数法によるKS方程式の自己無撞着場計算 平面波‐擬ポテンシャル法 ウルトラソフト型擬ポテンシャル 相関交換項の計算 : GGA , LDA 固有関数の計算 : 共役勾配法 , RMM-DIIS 分子動力学 : 共役勾配法, Verlet法 c-Si(64個)の電荷密度の様子[100]面 c-Si(64個)の電荷密度の様子[110]面 SiO2(24個)の電荷密度の様子[100]面 SiO2(24個)の電荷密度の様子[110]面

8 分子動力学 電子状態の計算を経験的ポテンシャル関数で近似し、原子系のダイナミクス・統計集合を解く方法 ポテンシャルの精度が課題

9 分子動力学ポテンシャル 二体ポテンシャル イオンポテンシャル 多体ポテンシャル（共有結合） Tight-Binding法 固定電荷型
電荷移動型 多体ポテンシャル（共有結合） クラスターポテンシャル ボンドオーダーポテンシャル Tight-Binding法

10 1．二体ポテンシャル１ Lennard-Jonesポテンシャル 原子の電荷の偏りによる双極子間の引力相互作用（van der Waals力）
希ガス（Ar, Ne）などに適用

11 1．二体ポテンシャル２ Morseポテンシャル クラスタの結合のために提案 金属等に応用
Expの関数系は後にEAMに使われる（第一原理計算と良く合う～Erolessi(1994)） でφ=-B (結合エネルギ） 引力項∝h , 斥力項∝|h|2　となるWolfsberg-Helmholtz近似を満たしている。

12 1、二体ポテンシャルの問題点 弾性定数の記述が貧弱 空孔生成エネルギや表面などback-bondの記述ができない
必ずCauchyの関係式C12=C44が成り立つ 空孔生成エネルギや表面などback-bondの記述ができない 環境に依存したポテンシャルが必要！！

13 分子動力学ポテンシャル 二体ポテンシャル イオンポテンシャル 多体ポテンシャル（共有結合） Tight-Binding法 固定電荷型
電荷移動型 多体ポテンシャル（共有結合） クラスターポテンシャル ボンドオーダーポテンシャル Tight-Binding法

14 2. イオンポテンシャル（固定電荷） BHSポテンシャル（SiO2など）
長距離ポテンシャルのため、カットオフ（スクリーニング）もしくはEwald和が必要 Qは実効電荷、σはイオン半径、bはイオンの堅さ 長距離力の計算により、計算時間は10倍程度になる

15 2. Charge-Variable MD SiO2, Si3N4の電荷の移動
バルクに対しては固定電荷で良いが、界面・表面に対しては、局所的な電荷の移動を考慮する必要がある。 電荷の移動＋原子間力の二つを解く必要がある。 Vashishita, 保川等 計算時間100倍程度

16 2, Charge-Variable MD

17 分子動力学ポテンシャル 二体ポテンシャル イオンポテンシャル 多体ポテンシャル（共有結合） Tight-Binding法 固定電荷型
電荷移動型 多体ポテンシャル（共有結合） クラスターポテンシャル ボンドオーダーポテンシャル Tight-Binding法

18 共有結合性材料の特徴 sp2,sp3結合等、方向性を持った結合 C, Siの最外殻はs, px, py, pzの4軌道 sp3混成軌道

19 3．Si（C）系の多体ポテンシャル Stillinger-Waberポテンシャル Tersoffポテンシャル “多体クラスターポテンシャル”
ダイヤモンド構造と液体構造にフィッティング Tersoffポテンシャル “ボンドオーダーポテンシャル” 多形態（ダイマ・グラファイト・ダイヤモンド・単純立方・体心立方など）のシリコンの構造にフィッティング 発展形としてBrenner(C-H)、Marty(Si-H)、大平(Si-H)ポテンシャル

20 3-1. 多体クラスターポテンシャル ある特定の構造について、結合距離（二体項）、結合角（三体項）、二面角（四体項）を合わせこむ。
合わせ込んだ構造に以外に遷移する場合（化学反応・相転移等）に対応できない！

21 SWポテンシャルの形 式(1)～(6) 二体項はrαβのみの関数 三体項はrαβ、 rαγ、 cosθの関数
hは1/3。　θ=109.47°でh=cosθ (ダイヤモンド構造が安定) α γ β θ

22 SWポテンシャルの合わせこみ ダイヤモンド構造の凝集エネルギ・格子定数・融点 液体構造

23 SWポテンシャルの物性値１ テキスト 表5 Diamond struct. Exp. SW T3 Lattice constants
5.429[Å] 5.431 5.432 Cohesive energy 4.63[eV] 4.63 Elastic constants C11 167.0[GPa] 161.6 142.5 C12 65.0[GPa] 81.6 75.4 C44 81.0[GPa] 60.3 69.0 Vacancy energy (3-4) [eV] 2.82 3.70 Interstitial energy Td (5-6) [eV] 5.25 3.45 H (4-5) [eV] 6.95 4.61 Bond-centered 5.99 5.86 Surface energy (100) 1.57[/1×1 cell] 1.416 1.367 Surface energy (111) 1.39[/1×1 cell] 1.158 0.957

24 SWポテンシャルの物性値２ クラスターの形状は合わない テキスト 図3

25 SWポテンシャルの物性値 バルク（ダイヤモンド構造）の性質は概ね合っている ダイヤモンド構造以外の他の形態は保証なし
クラスターの形状は合わない Si系のポテンシャルの評価はH. Balamaneらによって詳細に行われている。Phys. Rev. B46(1992), 2250 文献[7]

26 SWポテンシャルの応用 格子間シリコンの平衡濃度・拡散挙動（Brown, Maroudas 1994）
転位の移動度(Bulatov, Yip 1996) Epitaxy結晶成長(Gilmer 1992) インプランテーション(Rubia, Gilmer 1995) SiO2系への拡張(Jiang, Brown 1995)

27 3-2. ボンドオーダーポテンシャル 見かけは二体ポテンシャルbαβに多体項が含まれる。 配位数が大きくなれば、結合を作るための十分な価電子がなくなり、電子が非局在化して結合間で共鳴し、結合を弱める効果を表現出来る。 局所状態密度の二次モーメント近似より導かれる。

28 Tersoffポテンシャル ボンドオーダーに配位数依存（η,δで調整）とともに、角度依存性g(θ)を入れ、シリコンの物性値にフィッティングした。

29 Abellポテンシャル 現在提案されてる、EAMポテンシャル・FSポテンシャル・Tersoffポテンシャルは、すべて上の式で表される二次モーメント近似によるポテンシャルである。 →二次モーメント近似の導出(資料参照)

30 角度依存性・配位数依存性 120°付近でg(θ)が極小→ボンドオーダー増加 配位数が増大→ボンドオーダー減少 Bond-order g(θ)
配位数増加 g(θ) Bond-order 126° ボンドオーダの角度・配位数依存性 角度依存項

31 Tersoffポテンシャルの合わせこみ 異なる結晶構造（結合状態）のエネルギ・格子間距離
ダイマ→グラファイト→ダイヤモンド→単純立方→BCC→FCC 異なる二つのバージョンが存在 表面構造合わせこみ（T2）→あまり使われていない 弾性定数合わせこみ（T3)

32 Tersoffポテンシャルの物性値１ テキスト 表5 Diamond struct. Exp. SW T3
Lattice constants 5.429[Å] 5.431 5.432 Cohesive energy 4.63[eV] 4.63 Elastic constants C11 167.0[GPa] 161.6 142.5 C12 65.0[GPa] 81.6 75.4 C44 81.0[GPa] 60.3 69.0 Vacancy energy (3-4) [eV] 2.82 3.70 Interstitial energy Td (5-6) [eV] 5.25 3.45 H (4-5) [eV] 6.95 4.61 Bond-centered 5.99 5.86 Surface energy (100) 1.57[/1×1 cell] 1.416 1.367 Surface energy (111) 1.39[/1×1 cell] 1.158 0.957

33 Tersoffポテンシャルの物性値２ クラスターの形状は合わない テキスト 図3

34 Tersoffポテンシャルの物性値 バルク（ダイヤモンド構造）の性質は概ね合っている
ダイヤモンド構造以外についてもエネルギは保証→相転移にも対応できる クラスターの形状は合わない（T2はクラスターの形状は合うが、エネルギは×）。 Si系のポテンシャルの評価はH. Balamaneらによって詳細に行われている。Phys. Rev. B46(1992), 2250 文献[7]

35 Tersoffポテンシャルの応用 照射欠陥(Kitabatake 1993)
アモルファス化・固相成長(Motooka 1993, 2000) 圧力相変態(Mizushima, Yip 1994) インプランテーション（Nordlund 1998） アモルファス化（Bond-defect）(Marques 2001) 水素系への応用 Si-H Marty(1995), 大平(1994), 泉(2002) C-H Brenner(1990), (2002)

36 シリコン系へのTersoffポテンシャルの応用へ

37 FSとTersoffポテンシャルの類似性

38 金属系のボンドオーダーポテンシャル Tersoffポテンシャル FSポテンシャル EAMポテンシャル ※金属結合は不飽和な共有結合
結合角・距離（弱い）と配位数でボンドオーダーを決定 FSポテンシャル 結合距離と配位数でボンドオーダーを決定 EAMポテンシャル ただし、埋め込み関数（ボンドオーダーの効果に相当する項）を多項式やlog形にして汎用性・柔軟性を持たせる　→　FSと傾向は同じ　　

39 分子動力学ポテンシャル 二体ポテンシャル イオンポテンシャル 多体ポテンシャル（共有結合） Tight-Binding法 固定電荷型
電荷移動型 多体ポテンシャル（共有結合） クラスターポテンシャル ボンドオーダーポテンシャル Tight-Binding法

40 4. Tight-Binding分子動力学 経験的ポテンシャルでは、バルク・クラスターなど多種多様な物性を表現するのが難しい。扱えてもパラメータが多くなる。 量子効果を最もシンプルに取り込むTightーBinding法（TB法）が用いられている TB法はマトリクスの対角化が必要になるため、計算時間はO(N2)。

41 Tight-Bindingの解法の基礎１ n番目の固有状態ψ(n)の波動関数を原子軌道φiα（原子i, 軌道α）の線形結合で表す
波動関数に代入して永年方程式を得る Tight－Binding法ではHiαjβの計算にSlater-Kosterの式を使う。よって、積分不要。

42 Tight-Bindingの解法の基礎２ ハミルトニアンマトリックスHiα,jβの対角化（固有値問題）を行い、固有値（エネルギ）を求める。Order N3 小さいエネルギ準位から電子を二つずつ詰めて行くと、バンドエネルギが得られる。

43 Order-N Tb法 対角化をなくしてエネルギを求める Bond Order Potential法 Density Matrix法
Fermi Operator法 Global Density of State法 BOP法は、Abell-Tersoffポテンシャルの基盤となる考え方。

44 計算をはじめるにあたって Tersoffポテンシャル http://www.fml.t.u-tokyo.ac.jp/~izumi/MD/
　にてソースコード＋マニュアル配布 Tight-Bindingポテンシャル L.Colombo, Comp. Mater. Sci, 12(1998),278 ElsevierのHP上にソースコードあり。 ポテンシャルWG（日本材料学会分子動力学部門委員会WG「局所構造を反映した原子間ポテンシャル作成手法の調査」 ）

45 ポテンシャルの開発の試み

46 角度／距離依存性・配位数依存性でどこまで高精度なポテンシャルが出来るか？
シリコンクラスタ用Tersoffポテンシャル T2，T3→T4の開発 SiH4のシリコン表面反応ポテンシャル ボンドオーダーポテンシャルに基づいた形式 ポテンシャルの合わせこみデータの収集 Saddle-Pointの合わせこみ

47 Tersoffポテンシャルのクラスタへの応用
T2, T3ともにクラスターの構造が記述できない T3は角度依存性が大きく、表面に向かない T2は金属的結合で角度依存性が小さい Tersoff3の結晶成長計算 表面原子が109°の角度（ダイヤモンド構造の位置）で付着

48 クラスタ用ポテンシャルの作成 クラスターの構造間エネルギ差
Si3～Si6までのクラスタの形状間エネルギ差をエネルギ構造差定理を使うことによりfitting。 バルク（グラファイト・ダイヤモンド・単純立方・BCC）の格子間距離・エネルギを最低条件 →　ボンドオーダーのみを調節して、T3より角度依存性・配井数依存性を弱くして、T2の性質に近づける。 どの構造がより安定化？を重視

49 結果（クラスターエネルギ） Tersoffと比較して、エネルギの順序の再現が大幅に改善！

50 考察（Global Min. ） ●S1b Si4, Si5の最安定構造はAb-initioと一致する。Si6は三角柱が最安定となり、Ab-initioの八面体とは異なる。Si10もよりOpenな五角柱となり、Ab-initioと傾向が異なる。 ●S1c Si6, Si10の最安定構造の傾向があうポテンシャルを作成した。しかしながら、バルクの性質との共存は難しく。バルクの性質を犠牲にせざるを得なかった。具体的には、ダイヤモンド構造が最安定にならず、BCCが最安定となってしまった。

51 結論（クラスタ用ポテンシャル） ボンドオーダーを調整することによって、クラスタについては、T2・T3より優れたポテンシャルが作成できた。
しかしながら、バルクとクラスタの性質の双方を満足するポテンシャルは角度と配位数の調節だけでは不可能であった。 シリコンクラスタはダングリングボンドの数が多く、古典MDでは難しい。

52 実用的なCVD表面反応用Si-Hポテンシャルの開発
ポテンシャル作成における課題 どのような形状（安定or鞍点？）をあわせ込むのか？ フィッティングデータをどうやって集めるのか？ どのような関数形状を選ぶのか？ フィッティングアルゴリズムをどうするのか？

53 １、どのような形状（安定or鞍点？）をあわせ込むのか？
クラスタ→バルクへの成長を表現するため、様々な大きさのクラスタを合わせこみ対象に選定 汎用性を高めるために、同一の大きさでも、異なる形状のクラスタを合わせこみ対象に選定 化学反応過程を考え、表面反応の素過程に重要な現象のエネルギ障壁を合わせる

54 合わせ込んだクラスタ

55 合わせ込んだ反応過程 平衡値だけの合わせこみで良いのか？
結晶成長の表面反応をシミュレーションするためには、Saddle-Pointのエネルギを合わせる必要がある。 ※結晶成長に重要な素過程 SiH4の　解離吸着 Hの脱離

56 ２、フィッティングデータをどうやって集めるのか？
実験値の収集は大変・実験によりバラツキ 特に0Kでの結合エネルギの値の収集は困難 実験値は用いず、すべてGAUSSIAN98のG2セットの計算結果を用いた（実験値との平均誤差1.5％） １ HF/6-31G(d) Opt Freq 予備計算 ２ MP2(full)/6-31G(d) Opt 形状最適化 ３ MP4/6-311G(d,p) 以降エネルギ計算（補正） ４ MP4/6-311+G(d,p) ５ MP4/6-311G(2df,p) ６ QCISD(T)/6-311G(d,p) ７ MP2/6-311+G(3df,2p) G2の計算レシピ

57 ３、どのような関数形状を選ぶのか？ Tersoff３をベースにしたEBOPを関数形状に選定

58 ４、フィッティングアルゴリズムをどうするのか？
SiH4のエネルギ・結合距離・振動数を合わせこみ ( A, B, λA, λB, c, d ) SiH,SiH2,SiH3のエネルギ・結合距離を合わせこみ（ F1Si-H, F2Si-H で補正） Si2H6, Si2H4, Si2H2, HSiSiH3のエネルギと結合間距離を合わせこみ ( GijSi-Siで補正) Si3Hx～Si7Hxのエネルギがよく一致するように（特にエネルギの順序）合わせこみ（ GijSi-Siで補正） 4まで出来た後に、水素脱離とSiH4解離吸着をH-H-Si, H-Si-H, H-Si-Hの三体項を調整することによって合わせこむ １、はcrude MCによって、最適値を見つける。２，３，４は１の中からtry and errorで合わせこむ。５は４で設定していない三体項を使って合わせこむ。

59 誤差は、ほとんどが2%以下、平均で1.1％(Martyは3.3%) 環状のクラスタは苦手（Tersoff3の特性）
Fitting結果（誤差） 誤差は、ほとんどが2%以下、平均で1.1％(Martyは3.3%)　環状のクラスタは苦手（Tersoff3の特性）

60 Fitting結果 (Hの解離１) Eb=2.37eV (DFT:2.49) d1 d2 H解離のポテンシャルマップ
Si-Hのカットオフ距離と H-H-Si, H-Si-Hの三体項

61 Fitting結果 (SiH4の解離吸着) -0.70eV Eb=+0.56eV (ab-inito 0.5～0.6eV)
Saddle-point! Eb Si-Siのカットオフ距離と H-Si-Hの三体項

62 結論 どのような形状（安定or鞍点？）をあわせ込むのか？ →クラスタ形状と表面反応の活性化エネルギ
SiH4の表面反応に使用する分子動力学ポテンシャルを作成 どのような形状（安定or鞍点？）をあわせ込むのか？ →クラスタ形状と表面反応の活性化エネルギ フィッティングデータをどうやって集めるのか？ →GAUSSIAN９８による計算データ どのような関数形状を選ぶのか？ →EBOPに補正項を入れる。パラメータが多ければあわせこみの可能性も増える フィッティングアルゴリズムをどうするのか？ →Crude MC+try and error

63 ご清聴ありがとうございました

64

65 注意！（Cutoff関数の問題点） x 新しい結合が出来ると、急激にエネルギ低下＝ｚ局所的に強い引力が働く！！ H Si 4.0Å
Si-Hのcutoff 2.4Å Si H 新しい結合が出来ると、急激にエネルギ低下＝ｚ局所的に強い引力が働く！！ 結合を切るときは逆に局所的に大きな障壁 Cutoffにより、短距離の間に大きなエネルギ変化が生じる！！

66 Saddle-point energy計算
双方向Gradual　Ascent法（GA法） Bulatovらが開発したGA法（初期状態と最終状態をバネで結び、剛性を除々に強くすることによって、saddle-pointを見つける）を双方向に拡張 ② ④ 初期状態 GA 1 GA 3 GA 4 ① GA 2 ③ ⑤ GAバネ 最終状態 初期状態・最終状態ともに、少しずつGAで近づけてsaddle-pointを探す

67 Fitting手順 (Hの解離2) Eb=2.9eV (DFT:2.9) d1 d2 H解離のポテンシャルマップ Si-Hのカットオフ距離と
H-H-Si, H-Si-Hの三体項

68 Tight-Bindingの解法 n番目の固有状態ψ(n)の波動関数を原子軌道φiα（原子i,軌道α）の線形結合で表す
波動関数に代入して永年方程式を得る ハミルトニアンマトリックスHiα,jβの対角化（固有値問題）を行い、固有値（エネルギ）を求める。

69 状態密度と結合エネルギー Energy 原子 固体 エネルギの状態密度分布 ntotal(E)
結合は、原子状態より低いエネルギの分子軌道が作られることによって生じる。 →　状態密度の分布形状が結合エネルギを決める。

70 状態密度のモーメント E n(E) 二次モーメント近似(軌道数1) E n(E) W/2 -W/2
W/2 -W/2 1/W しかし、これではLOCALな結合の記述が出来ない

71 局所状態密度の定義 個々の分子軌道からの全体の結合エネルギへの寄与

72 局所状態密度と結合形態の関係 モーメント定理
局所状態密度のモーメントは、該当原子から出発し、戻ってくるすべてのパスのN回の“HOP”（ハミルトニアンマトリックスの積）の和によって決まる。 4-Hops 四次モーメント 2-Hops 二次モーメント

73 二次モーメント近似 局所状態密度の二次モーメントは近接原子数によってのみ決まる。よって、結合エネルギは以下の式のように配位数の1/2乗となる。→二次モーメント近似 Z Ebind

74

75 局所構造を反映した 原子間ポテンシャル作成手法の調査
　　　　　 主旨 　　材料を構成する原子構造体の外力に対する応答として，き裂先端場から射出される転位挙動や結晶粒界近傍に生じる種々の力学現象などが，分子動力学法により正確に描写できるようになった．しかしながら，非可逆な変形応答やその極限状態 としての破壊事象は，本質的に不均質な構造を対象とした現象であり，これまでの平均的で均質系を仮定したマクロ系の材料特性のみを表現しうる従来の原子間相互作用の記述では，本質的に捉えられない． そこで，本WGでは従来からあるポテンシャルの問題点などを整理するとともに，不均 質な局所構造の影響を反映した原子間ポテンシャル作成のための調査を目的とする．

76 活動内容 現状の経験的，半経験的ポテンシャルの問題点の整理
　　　 Tersoff型、 MEAM(modified EAM)型、 Tight-binding法 非平衡現象に対応するポテンシャルのフィッテング手法の開発 汎用第一原理ＭＤの利用技術の確立 ポテンシャルの第一原理ＭＤへの繰み込み手法 ポテンシャルのデータベース構築の検討 データベースの作成 大学／企業で必要と考えられるポテンシャルの調査

77 活動方法 1. 年4回程度のミーティングを実施し，あらかじめ決められた報告者が，事前に で議事内容を連絡するとともに当日会議の進行を行う． 2. 報告の内容は，設定された調査の経過報告や文献紹介など，特に制約は設けない． 3. 学会のオーガナイズドセッションやワークショップ，フォーラムなどを利用した 活動結果の公表を定期的に行う．