共有結合性ポテンシャル - 酒井・泉研究室
Download
Report
Transcript 共有結合性ポテンシャル - 酒井・泉研究室
分子動力学冬季講習
ポテンシャル3 共有結合性材料
東京大学工学系研究科
泉 聡志
概要
近年のナノテクノロジーの進歩により、
シリコンプロセスやカーボンナノチュー
ブのMD解析が脚光を浴びている。ここ
では、シリコン・炭素系のStillingerWaber、Tersoff、Brennerポテンシャル
をその原理に踏み込んで解説し、応用
例をいくつか紹介する。
講習内容
1. 共有結合ポテンシャルの概要
2. ポテンシャルの紹介・微分形等の実践に
必要な知識
1. Stillinger-Waberポテンシャル(Si)
2. Tersoffポテンシャル(Si, C)
3. Brennerポテンシャル(C, H)
3. ポテンシャルの今後の展望
共有結合性材料の特徴
sp2,sp3結合等、方向性を持った結合
C, Siの最外殻はs, px, py, pzの4軌道
sp3混成軌道
sp2混成軌道
1 ( sp )
3
1 ( sp )
2
2 ( sp )
2
3 ( sp )
2
1 / 3 ( 2 s )
2 / 3 ( 2 p x )
1 / 3 ( 2 s ) 1 / 6 ( 2 p x )
2 ( sp )
3
1 / 2 ( 2 p y )
1 / 3 ( 2 s ) 1 / 6 ( 2 p x ) 1 / 2 ( 2 p y )
3 ( sp )
3
4 ( sp )
3
1 / 2 ( 2 s ) ( 2 p
1 / 2 ( 2 s ) ( 2 p
1 / 2 ( 2 s ) ( 2 p
)
)
)
1 / 2 (2 s) (2 px ) (2 p y ) (2 p z )
x
) (2 p y ) (2 p z
x
) (2 p y ) (2 pz
x
) (2 p y ) (2 p z
分子動力学の適用(共有結合系)
Ab-initio 分子動力学法 (CP法)
Car-Parinello
Tight-Binding分子動力学法
→Order-N法
C.Z.Wang, et al.
経験的分子動力学
Tersoff, Stillinger- Waber, Brenner, Pettifor, ….
近似 速度
Si（C）系の経験的ポテンシャル
Stillinger-Waberポテンシャル
“多体クラスターポテンシャル”
ダイヤモンド構造と液体構造にフィッティング
Tersoffポテンシャル
“ボンドオーダーポテンシャル”
多形態（ダイマ・グラファイト・ダイヤモンド・単純立方・
体心立方など）のシリコンの構造にフィッティング
発展形としてBrenner(C-H)、Marty(Si-H)、大平(Si-H)
ポテンシャル
多体クラスターポテンシャル
V
V (R , R ) V (R , R , R ) V (R , R , R , R )
2b
2 body
3b
3 body
4b
4 body
ある特定の構造について、結合距離（二体項）、結合
角（三体項）、二面角（四体項）を合わせこむ。
合わせ込んだ構造に以外に遷移
する場合（化学反応・相転移等）
に対応できない！
ボンドオーダーポテンシャル
V
A exp
r
A
b
B exp B r
b
: bond order
Abell (1985 )ポテンシャル
見かけは二体ポテンシャルbαβに多体
項が含まれる。
b Z
1 / 2
: Z は配位数
配位数が大きくなれば、結合を作るための十分な価電子
がなくなり、電子が非局在化して結合間で共鳴し、結合を
弱める効果を表現出来る。
局所状態密度の二次モーメント近似より導かれる。
講習内容
1. 共有結合ポテンシャルの概要
2. ポテンシャルの紹介・微分形等の実践に
必要な知識
1. Stillinger-Waberポテンシャル(Si)
2. Tersoffポテンシャル(Si, C)
3. Brennerポテンシャル(C, H)
3. ポテンシャルの今後の展望
SWポテンシャルの形
cutoff
p
p
2
B r
a f c ( r ) f c ( r ) h cos
angular dependent
2 body
fc (r ) A r
A
B
3 body
式(1)～(6)
•二体項はrαβのみの関数
•三体項はrαβ、 rαγ、 cosθの関数
•hは1/3。 θ=109.47°でh=cosθ
(ダイヤモンド構造が安定)
β
θ
γ
α
SWポテンシャルの合わせこみ
ダイヤモンド構造の凝集エネルギ・格子定
数・融点
液体構造
SWポテンシャルの微分
三体項の微分
力は基本的に二体ポテンシャルと
同じヘルマン･ファインマン則
Fi
V
ri
ただし、三体項V3はrαβ、 rαγ、 cosθの関数。よって、これらを
経由した偏微分を行う必要がある。一つの三体項に対して、
力はα,β,γの三つの原子に作用する。
β
もちろん、三角定理により、cosθをrαβ、 rαγ 、 rβγの関
数として、微分してもよい。
θ
γ
α
SWポテンシャルの微分
具体的な微分形（一つの三体項）
Fi
Fi
Fi
V3
ri
V3
ri
V3
ri
V3 r
r
ri
V3 r
r
V3 r
r
ri
ri
V3 r
r
ri
V3
cos
V3
cos
ri
cos
cos
ri
V3
cos
cos
ri
(10 )
(11 )
(12 )
原子間距離と方向余弦の微分形
テキスト式(16)～式(18)参照
SWポテンシャルの物性値１
テキスト 表5
Diamond struct.
Exp.
SW
T3
Lattice constants
5.429[Å]
5.431
5.432
Cohesive energy
4.63[eV]
4.63
4.63
161.6
142.5
C12 65.0[GPa]
81.6
75.4
C44 81.0[GPa]
60.3
69.0
(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
(4-5) [eV]
5.99
5.86
1.416
1.367
Elastic constants C11 167.0[GPa]
Vacancy energy
Surface energy (100) 1.57[/1×1
cell]
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)
講習内容
1. 共有結合ポテンシャルの概要
2. ポテンシャルの紹介・微分形等の実践に
必要な知識
1. Stillinger-Waberポテンシャル(Si)
2. Tersoffポテンシャル(Si, C)
3. Brennerポテンシャル(C, H)
3. ポテンシャルの今後の展望
ボンドオーダーポテンシャル
V
A exp
A
r
b
B exp B r
b
: bond order
Abell (1985 )ポテンシャル
b
Z
1 / 2
: Z は配位数
配位数が大きくなれば、結合を作るための十分な価電子
がなくなり、電子が非局在化して結合間で共鳴し、結合を
弱める効果を表現出来る。
局所状態密度の二次モーメント近似より導かれる。
Tersoffポテンシャル
V A exp A r
b B exp B r
bond order
b
1
,
g ( ) exp p ( r r
( 23 ) ( 25 )
k
angle
bond - length
すべての近接について
2
2
c
c
g ( ) a 1 2 2
2
d
d
(
h
cos
)
( 21 )
q
の和 配位数 ( Z )
b
Z
1 / 2
ボンドオーダーに配位数依存（η,δで調整）とともに、角度依
存性g(θ)を入れ、シリコンの物性値にフィッティングした。
角度依存性・配位数依存性
g(θ)
120°付近でg(θ)が極小→ボンドオーダー増加
配位数が増大→ボンドオーダー減少
126°
角度依存項
Bond-order
配位数増加
ボンドオーダの角度・配位数依存性
Tersoffポテンシャルの合わせこみ
異なる結晶構造（結合状態）のエネルギ・格子間距
離
ダイマ→グラファイト→ダイヤモンド→単純立方→BCC→FCC
異なる二つのバージョンが存在
表面構造合わせこみ（T2）→あまり使われていない
弾性定数合わせこみ（T3)
Tersoffポテンシャルの微分
φ
複雑な依存関係に
注意して、ヘルマン
ファイマン則に従っ
て偏微分する。
Fi
V
ri
b
αβ
②
ζ
fc
αβ
αγ
fc
④
r
αβ
αγ
VR
r
α
αβ
αβ
g (θ）
r
⑤
③
cosθ
VA
①
r
αγ
④
Tersoffポテンシャルの微分
Fi
Fi
Fi
ri
ri
r
r
r
cos
r
ri
r
ri
ri
ri
, r
, cos
5
1
2
3
4
r
r
ri
ri
r
ri
r
,
r
r
ri
cos
cos ri
b
αβ
r
αβ
αβ
式(27)～(29)
β
φ
cos
cos ri
ζ
θ
γ
α
r
αβ
r
αγ
cosθ
フローチャート
注意点
αβ
ボンドオーダb を計算してか
ら力の計算を行う。
三体項はα≠β≠γで計算
αβ
配位数が1なら b ≡1
αβ
配位数が0なら b ≡0
微分が発散してしまうので注
意!!
二体項は独立に計算する。
図2 参照
Tersoffポテンシャルの物性値１
テキスト 表5
Diamond struct.
Exp.
SW
T3
Lattice constants
5.429[Å]
5.431
5.432
Cohesive energy
4.63[eV]
4.63
4.63
161.6
142.5
C12 65.0[GPa]
81.6
75.4
C44 81.0[GPa]
60.3
69.0
(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
(4-5) [eV]
5.99
5.86
1.416
1.367
Elastic constants C11 167.0[GPa]
Vacancy energy
Surface energy (100) 1.57[/1×1
cell]
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プロセス
300K SiH4の解離吸着過程
AVI
ナノインデンテーション
大阪大学渋谷研究室
シリコン基板へのダイヤモンド圧子の押し込み
シリコン原子のMBE（Molecular
Beam Epitaxy）
基板温度８００K
インプランテーション
10keV Ar impacted to the silicon substrate at 300K
講習内容
1. 共有結合ポテンシャルの概要
2. ポテンシャルの紹介・微分形等の実践に
必要な知識
1. Stillinger-Waberポテンシャル(Si)
2. Tersoffポテンシャル(Si, C)
3. Brennerポテンシャル(C, H)
3. ポテンシャルの今後の展望
Brennerが考えたこと～Tersoffポテンシャルは一
重、二重、三重結合が混在した場合の記述が貧弱
１、ラジカルの過結合
右の状態は、一重結合とラジカル軌道が形成
→TersoffはV=1/2(Vab+Vba)なので、一重と二重
の中間の結合になってしまう。非対称性の考慮が
必要
２、C-C共役結合性
ケクレ構造はπ電子が非局在化して共役結合が生
じるが(CH3)2C=C(CH3)2は二重結合が生じ、非共
役系である。
→Tersoffはこの二つが区別できない。非局所の
効果が必要
Brennerポテンシャル形
炭化水素の様々な効果を補正項（FCC, HCC, HCH）に
入れ込む
合わせこみは多種の分子（5０種類程度）、バルク
（Tersoffと同様に多形態）の性質で行う。
V A exp A r
b B exp B r
bond order
b
1
2
b
b
1
b
Fcc ( N 1 , N 2 , N 3 )
H CC ( CH ) ( N 4 , N 5 )
N 1～ N 5は何らかの近接原子数
しかし、パラメータ
が非常に多い！！
Brennerポテンシャルの物性値
※表6 色々なバージョンが存在する。 膨大な炭化水素を表現出来る。
Diamond struct.
Exp.
Tersoff
Brenner
Lattice constants
3.566[Å]
3.56
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
Exp.
Tersoff
Brenner
Bond length
1.42 [Å]
1.46
1.45
Cohesive energy
7.38[eV]
7.40
7.38
Elastic constants C11 1060[GPa]
1210
1037
Elastic constants C12 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の変形
大阪大学渋谷研究室
講習内容
1. 共有結合ポテンシャルの概要
2. ポテンシャルの紹介・微分形等の実践に
必要な知識
1. Stillinger-Waberポテンシャル(Si)
2. Tersoffポテンシャル(Si, C)
3. Brennerポテンシャル(C, H)
3. ポテンシャルの今後の展望
Tight-Binding分子動力学
経験的ポテンシャルでは、バルク・クラス
ターなど多種多様な物性を表現するのが
難しい。扱えてもパラメータが多くなる。
量子効果を最もシンプルに取り込む
TightーBinding法（TB法）が用いられてい
る
TB法はマトリクスの対角化が必要になるた
め、計算時間はO(N2)。
Tight-Bindingの解法の基礎１
n番目の固有状態ψ(n)の波動関数を原子軌道φiα（原子i,
軌道α）の線形結合で表す
分子軌道 :
(n)
C i i
(n)
i
波動関数に代入して永年方程式を得る Hˆ
N
j
H i , j C j
(n)
(n)
(n)
C i
H
i , j
(n)
*
i Hˆ j dV
(n)
Tight－Binding法ではHiαjβの計算にSlater-Kosterの式を
使う。よって、積分不要。
(n)
Tight-Bindingの解法の基礎２
ハミルトニアンマトリックスHiα,jβの対角化（固有値問題）を
行い、固有値（エネルギ）を求める。
N
H i , j C j
(n)
(n)
(n)
C i
j
小さいエネルギ準位から電子を二つずつ詰めて行くと、
バンドエネルギが得られる。
E band
i , j
z
2 C i C j H i , j
(n)
n 1
( n )*
Order-N Tb法
対角化をなくしてエネルギを求める
Bond Order Potential法
Density Matrix法
Fermi Operator法
Global Density of State法
BOP法は、Abell-Tersoffポテンシャルの基盤となる考え方。
Abellポテンシャルの導出
bond order : b Z
1 / 2
: Z は配位数
現在提案されてる、EAMポテンシャル・FS
ポテンシャル・Tersoffポテンシャルは、す
べて上の式で表される二次モーメント近
似によるポテンシャルである。
→二次モーメント近似の導出
状態密度と結合エネルギー
原子
固体
Energy
Hˆ
(n)
(n)
エネルギの状
態密度分布
(n)
波動関数
ntotal(E)
原子軌道 : i
分子軌道 :
(n)
C i i
(n)
i
結合は、原子状態より低いエネルギの分子軌道が作られ
E
ることによって生じる。
U
En total ( E ) dE
f
→ 状態密度の分布形状が結合エネルギを決める。
状態密度の形状を
分布のモーメントで表す
E
一次モーメント：
1
En ( E ) dE
二次モーメント：
2
E n ( E ) dE 偏差（平均二乗幅）
三次モーメント：
3
E
平均エネルギ
2
3
n ( E ) dE ゆがみ
n(E
)
E
二次モーメント近似(軌道数1)
1 0,
W/2
2
1
W
2
12
0
0
n(E
)
1/W
-W/2
E
Ef
3
2
E2 1
W
En ( E ) dE 2
E
dE 2
W / 2
W
2
W
4
W / 2
E f 0
2
1
2
しかし、これではLOCALな結合の記述が出来ない
局所状態密度の定義
個々の分子軌道からの全体の結合エネル
ギへの寄与
n total ( E )
(E
(n)
)
n
n i ( E )
C
(n)
i
2
(E
U band 2
(n)
):
Ef
En total ( E ) dE
(n)
C i : 分子軌道の係数
n
分子軌道 :
(n)
C i i
(n)
i
局所状態密度と結合形態の関係
（モーメント定理）
局所状態密度のモーメントは、該当原子から出発し、
戻ってくるすべてのパスのN回の“HOP”（ハミルトニ
アンマトリックスの積）の和によって決まる。
2 i
2-Hops 二次
モーメント
H
i , j
H
j ,i
j
p i
H
j1 1 , j p p
i , j1 1
H
j1 1 , j 2 2
H i , j : ハミルトニアンマトリ
4-Hops
四次モーメ
ント
H
i , j
H
j p 1 p 1 , i
ックス
*
i Hˆ j dV
二次モーメント近似
局所状態密度の二次モーメントは近接原
子数によってのみ決まる。よって、結合エ
ネルギは以下の式のように配位数の1/2
乗となる。
binding energy
2 i
Ebind
Z
Z
計算をはじめるにあたって
Tersoffポテンシャル
応相談 [email protected]
Tight-Bindingポテンシャル
L.Colombo, Comp. Mater. Sci, 12(1998),278
ElsevierのHP上にソースコードあり。
ポテンシャルWG（日本材料学会分子動力学部門委員会WG
「局所構造を反映した原子間ポテンシャル作成手法の調査」 ）
http://www.fml.t.u-tokyo.ac.jp/wg-cmd/