分子動力学冬季講習 ポテンシャル3 共有結合性材料 分子動力学冬季講習 ポテンシャル3 共有結合性材料 東京大学工学系研究科 泉 聡志
概要 近年のナノテクノロジーの進歩により、シリコンプロセスやカーボンナノチューブのMD解析が脚光を浴びている。ここでは、シリコン・炭素系のStillinger-Waber、Tersoff、Brennerポテンシャルをその原理に踏み込んで解説し、応用例をいくつか紹介する。
講習内容 共有結合ポテンシャルの概要 ポテンシャルの紹介・微分形等の実践に必要な知識 ポテンシャルの今後の展望 Stillinger-Waberポテンシャル(Si) Tersoffポテンシャル(Si, C) Brennerポテンシャル(C, H) ポテンシャルの今後の展望
共有結合性材料の特徴 sp2,sp3結合等、方向性を持った結合 C, Siの最外殻はs, px, py, pzの4軌道 sp3混成軌道
分子動力学の適用(共有結合系) Ab-initio 分子動力学法 (CP法) Tight-Binding分子動力学法 →Order-N法 Car-Parinello Tight-Binding分子動力学法 →Order-N法 C.Z.Wang, et al. 経験的分子動力学 Tersoff, Stillinger- Waber, Brenner, Pettifor, …. 近似 速度
Si（C）系の経験的ポテンシャル Stillinger-Waberポテンシャル Tersoffポテンシャル “多体クラスターポテンシャル” ダイヤモンド構造と液体構造にフィッティング Tersoffポテンシャル “ボンドオーダーポテンシャル” 多形態（ダイマ・グラファイト・ダイヤモンド・単純立方・体心立方など）のシリコンの構造にフィッティング 発展形としてBrenner(C-H)、Marty(Si-H)、大平(Si-H)ポテンシャル
多体クラスターポテンシャル ある特定の構造について、結合距離（二体項）、結合角（三体項）、二面角（四体項）を合わせこむ。 合わせ込んだ構造に以外に遷移する場合（化学反応・相転移等）に対応できない！
ボンドオーダーポテンシャル 見かけは二体ポテンシャルbαβに多体項が含まれる。 配位数が大きくなれば、結合を作るための十分な価電子がなくなり、電子が非局在化して結合間で共鳴し、結合を弱める効果を表現出来る。 局所状態密度の二次モーメント近似より導かれる。
講習内容 共有結合ポテンシャルの概要 ポテンシャルの紹介・微分形等の実践に必要な知識 ポテンシャルの今後の展望 Stillinger-Waberポテンシャル(Si) Tersoffポテンシャル(Si, C) Brennerポテンシャル(C, H) ポテンシャルの今後の展望
SWポテンシャルの形 式(1)～(6) 二体項はrαβのみの関数 三体項はrαβ、 rαγ、 cosθの関数 hは1/3。 θ=109.47°でh=cosθ (ダイヤモンド構造が安定) α γ β θ
SWポテンシャルの合わせこみ ダイヤモンド構造の凝集エネルギ・格子定数・融点 液体構造
SWポテンシャルの微分 三体項の微分 力は基本的に二体ポテンシャルと同じヘルマン･ファインマン則 ただし、三体項V3はrαβ、 rαγ、 cosθの関数。よって、これらを経由した偏微分を行う必要がある。一つの三体項に対して、力はα,β,γの三つの原子に作用する。 α γ β θ もちろん、三角定理により、cosθをrαβ、 rαγ 、 rβγの関数として、微分してもよい。
SWポテンシャルの微分 具体的な微分形（一つの三体項） 原子間距離と方向余弦の微分形 テキスト式(16)～式(18)参照
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
SWポテンシャルの物性値２ クラスターの形状は合わない テキスト 図3
SWポテンシャルの物性値 バルク（ダイヤモンド構造）の性質は概ね合っている ダイヤモンド構造以外の他の形態は保証なし クラスターの形状は合わない Si系のポテンシャルの評価はH. Balamaneらによって詳細に行われている。Phys. Rev. B46(1992), 2250 文献[7]
SWポテンシャルの応用 格子間シリコンの平衡濃度・拡散挙動（Brown, Maroudas 1994） 転位の移動度(Bulatov, Yip 1996) Epitaxy結晶成長(Gilmer 1992) インプランテーション(Rubia, Gilmer 1995) SiO2系への拡張(Jiang, Brown 1995)
講習内容 共有結合ポテンシャルの概要 ポテンシャルの紹介・微分形等の実践に必要な知識 ポテンシャルの今後の展望 Stillinger-Waberポテンシャル(Si) Tersoffポテンシャル(Si, C) Brennerポテンシャル(C, H) ポテンシャルの今後の展望
ボンドオーダーポテンシャル 配位数が大きくなれば、結合を作るための十分な価電子がなくなり、電子が非局在化して結合間で共鳴し、結合を弱める効果を表現出来る。 局所状態密度の二次モーメント近似より導かれる。
Tersoffポテンシャル ボンドオーダーに配位数依存（η,δで調整）とともに、角度依存性g(θ)を入れ、シリコンの物性値にフィッティングした。
角度依存性・配位数依存性 120°付近でg(θ)が極小→ボンドオーダー増加 配位数が増大→ボンドオーダー減少 Bond-order g(θ) 配位数増加 g(θ) Bond-order 126° ボンドオーダの角度・配位数依存性 角度依存項
Tersoffポテンシャルの合わせこみ 異なる結晶構造（結合状態）のエネルギ・格子間距離 ダイマ→グラファイト→ダイヤモンド→単純立方→BCC→FCC 異なる二つのバージョンが存在 表面構造合わせこみ（T2）→あまり使われていない 弾性定数合わせこみ（T3)
Tersoffポテンシャルの微分 φ bαβ ζαβ gα(θ） fcαβ VR VA rαβ fcαγ rαγ cosθ 複雑な依存関係に注意して、ヘルマンファイマン則に従って偏微分する。 ② ① ③ ④ ④ ⑤
Tersoffポテンシャルの微分 φ bαβ ζαβ rαβ rαγ cosθ 式(27)～(29) α γ β θ
フローチャート 注意点 ボンドオーダbαβを計算してから力の計算を行う。 三体項はα≠β≠γで計算 配位数が1なら bαβ≡1 微分が発散してしまうので注意!! 二体項は独立に計算する。 図2 参照
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
Tersoffポテンシャルの物性値２ クラスターの形状は合わない テキスト 図3
Tersoffポテンシャルの物性値 バルク（ダイヤモンド構造）の性質は概ね合っている ダイヤモンド構造以外についてもエネルギは保証→相転移にも対応できる クラスターの形状は合わない（T2はクラスターの形状は合うが、エネルギは×）。 Si系のポテンシャルの評価はH. Balamaneらによって詳細に行われている。Phys. Rev. B46(1992), 2250 文献[7]
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)
シリコン系へのTersoffポテンシャルの応用
アモルファスシリコンの固相成長（SPE) ポテンシャル Tersoff 計算条件 NPT, 1600K
SiH4のCVDプロセス AVI 300K SiH4の解離吸着過程
ナノインデンテーション 大阪大学渋谷研究室 シリコン基板へのダイヤモンド圧子の押し込み
シリコン原子のMBE（Molecular Beam Epitaxy） 基板温度８００K
インプランテーション 10keV Ar impacted to the silicon substrate at 300K
講習内容 共有結合ポテンシャルの概要 ポテンシャルの紹介・微分形等の実践に必要な知識 ポテンシャルの今後の展望 Stillinger-Waberポテンシャル(Si) Tersoffポテンシャル(Si, C) Brennerポテンシャル(C, H) ポテンシャルの今後の展望
Brennerが考えたこと～Tersoffポテンシャルは一重、二重、三重結合が混在した場合の記述が貧弱 １、ラジカルの過結合 右の状態は、一重結合とラジカル軌道が形成 →TersoffはV=1/2(Vab+Vba)なので、一重と二重の中間の結合になってしまう。非対称性の考慮が必要 ２、C-C共役結合性 ケクレ構造はπ電子が非局在化して共役結合が生じるが(CH3)2C=C(CH3)2は二重結合が生じ、非共役系である。 →Tersoffはこの二つが区別できない。非局所の効果が必要
Brennerポテンシャル形 炭化水素の様々な効果を補正項（FCC, HCC, HCH）に入れ込む 合わせこみは多種の分子（5０種類程度）、バルク（Tersoffと同様に多形態）の性質で行う。 しかし、パラメータが非常に多い！！
Brennerポテンシャルの物性値 ※表6 色々なバージョンが存在する。 膨大な炭化水素を表現出来る。 Diamond struct. ※表6 色々なバージョンが存在する。 膨大な炭化水素を表現出来る。 Diamond struct. Exp. Tersoff Brenner Lattice constants 3.566[Å] 3.56 Cohesive energy 7.35[eV] 7.37 7.32 Elastic constants C11 1081[GPa] 1090 613(1078*) Elastic constants C12 125[GPa] 120 405(131*) Elastic constants C44 579[GPa] 640 631(680*) Vacancy energy (7.2)[eV] 4.32 7.4 Graphite Bond length 1.42 [Å] 1.46 1.45 7.38[eV] 7.40 7.38 1060[GPa] 1210 1037 180[GPa] -190 155
Brennerポテンシャルの応用 CVDプロセス(Brenner 1990) ダイヤモンド表面の摩擦(Harrison 1992) C60の性質(Brenner 1991) 炭素クラスター(Yamaguchi, Maruyama 1997) 問題点：バルクの弾性定数の再現が悪い
第二世代Brennerポテンシャル バルクの弾性定数の再現のため、2002年に新しいポテンシャルが提案された 変更点 斥力・引力の二体項の関数系を変更 角度依存項を関数ではなく、離散値（合わせこんだ値）の六次のスプライン関数に変更 その他は同じであるが、パラメータは１から再フィッティング
Brennerポテンシャルの応用 ナノインデンテーション(Sinnott 1997) 多結晶ダイヤモンドの破壊(Shenderova 2000) カーボンナノチューブ(Che 1999, Mao 1999, Srivastava 1999) ナノテクノロジーへの応用が目立つ
C60の生成過程 東京大学丸山研究室
CNTの生成過程１ 東京大学丸山研究室
CNTの生成過程２（Niクラスタ触媒） 東京大学丸山研究室
CNTの変形 大阪大学渋谷研究室
講習内容 共有結合ポテンシャルの概要 ポテンシャルの紹介・微分形等の実践に必要な知識 ポテンシャルの今後の展望 Stillinger-Waberポテンシャル(Si) Tersoffポテンシャル(Si, C) Brennerポテンシャル(C, H) ポテンシャルの今後の展望
Tight-Binding分子動力学 経験的ポテンシャルでは、バルク・クラスターなど多種多様な物性を表現するのが難しい。扱えてもパラメータが多くなる。 量子効果を最もシンプルに取り込むTightーBinding法（TB法）が用いられている TB法はマトリクスの対角化が必要になるため、計算時間はO(N2)。
Tight-Bindingの解法の基礎１ n番目の固有状態ψ(n)の波動関数を原子軌道φiα（原子i, 軌道α）の線形結合で表す 波動関数に代入して永年方程式を得る Tight－Binding法ではHiαjβの計算にSlater-Kosterの式を使う。よって、積分不要。
Tight-Bindingの解法の基礎２ ハミルトニアンマトリックスHiα,jβの対角化（固有値問題）を行い、固有値（エネルギ）を求める。 小さいエネルギ準位から電子を二つずつ詰めて行くと、バンドエネルギが得られる。
Order-N Tb法 対角化をなくしてエネルギを求める Bond Order Potential法 Density Matrix法 Fermi Operator法 Global Density of State法 BOP法は、Abell-Tersoffポテンシャルの基盤となる考え方。
Abellポテンシャルの導出 現在提案されてる、EAMポテンシャル・FSポテンシャル・Tersoffポテンシャルは、すべて上の式で表される二次モーメント近似によるポテンシャルである。 →二次モーメント近似の導出
状態密度と結合エネルギー Energy 原子 固体 エネルギの状態密度分布 ntotal(E) 結合は、原子状態より低いエネルギの分子軌道が作られることによって生じる。 → 状態密度の分布形状が結合エネルギを決める。
状態密度の形状を 分布のモーメントで表す E n(E) 二次モーメント近似(軌道数1) E n(E) W/2 -W/2 W/2 -W/2 1/W しかし、これではLOCALな結合の記述が出来ない
局所状態密度の定義 個々の分子軌道からの全体の結合エネルギへの寄与
局所状態密度と結合形態の関係（モーメント定理） 局所状態密度のモーメントは、該当原子から出発し、戻ってくるすべてのパスのN回の“HOP”（ハミルトニアンマトリックスの積）の和によって決まる。 2-Hops 二次モーメント 4-Hops 四次モーメント
二次モーメント近似 局所状態密度の二次モーメントは近接原子数によってのみ決まる。よって、結合エネルギは以下の式のように配位数の1/2乗となる。 Z Ebind
計算をはじめるにあたって Tersoffポテンシャル 応相談 izumi@fml.t.u-tokyo.ac.jp Tight-Bindingポテンシャル L.Colombo, Comp. Mater. Sci, 12(1998),278 ElsevierのHP上にソースコードあり。 ポテンシャルWG（日本材料学会分子動力学部門委員会WG「局所構造を反映した原子間ポテンシャル作成手法の調査」 ） http://www.fml.t.u-tokyo.ac.jp/wg-cmd/