Transcript ポテンシャルの基礎と応用（計算工学会講習会資料）

よくわかる計算工学最前線
「原子間ポテンシャルの
基礎から応用まで 」
2003-01-17
東京大学工学系研究科機械工
学専攻 泉 聡志
日本材料学会分子動力学部門
委員会
• ポテンシャルWG
– 局所構造を反映した原子間ポテンシャル作成手法の
調査
– Tersoff型、 MEAM(modified EAM)型、 Tight-binding
法
– データベース化
• 分子動力学冬季講座
– fcc, bcc系のポテンシャル
– HCP系のポテンシャル
– 共有結合系のポテンシャル
概要
• ポテンシャルの基礎（思想）
– 二体ポテンシャル
– イオンポテンシャル
– 多体ポテンシャル(共有結合系)
• クラスターポテンシャル
• ボンドオーダーポテンシャル
– Tight-Binding法
• ポテンシャルの応用
– シリコンクラスター用Tersoffポテンシャルの作成
– SiH4表面反応用ポテンシャルの作成
原子間ポテンシャルの階層
Density Functional Theory ( GGA, LDA )
CASTEP, Wien97
(Order-N) Tight-Binding
C.Z.Wang, D. G. Pettifor …
Charge-Variable Potential
Empirical Potential （Bond-order）
Vashishita…
Tersoff, Finnis…
Empirical 2-body Potential
Lennard-Jones, Morse
近似 速度
密度汎関数法（DFT）
•
•
•
様々な近似（擬ポテンシャル・交換相関
エネルギ）に基づき、電子状態をセルフコ
ンシストに算出
原子個数 50～100個程度（状態の緩和
＋エネルギ計算）
汎用ソフトが整備されつつある
–
–
VASP
CASTEP
密度汎関数法（DFT）
 1




v
(
r
)

  i (r )   i  i (r )
eff
 2

v eff ( r )  v ( r ) 


dr '
n (r ' )
| r  r '|
v ( r ) : 原子核が電子に作用す

 E xc ( n ( r ))
 n (r )
るポテンシャルエネル
ギ（擬ポテンシャル）
n ( r ) : 電子密度 、 E xc ( n ( r )) : 交換相関エネルギ
※ v eff ( r )が n ( r )に依存するので、求め
電子密度 n ( r )が v eff ( r )に用いた
た  i ( r )から占有状態を考慮し
n ( r )とセルフコンシストに
て得られた
なる必要がある。
VASPによるc-Si,SiO2（β-c）電荷密度の計算
VASP（第一原理計算パッケージ）
• 密度汎関数法によるKS方程式の自己無撞着場計算
• 平面波‐擬ポテンシャル法
• ウルトラソフト型擬ポテンシャル
• 相関交換項の計算 : GGA , LDA
• 固有関数の計算 : 共役勾配法 , RMM-DIIS
• 分子動力学 : 共役勾配法, Verlet法
c-Si(64個)の電荷密度の様子[100]面
c-Si(64個)の電荷密度の様子[110]面 SiO2(24個)の電荷密度の様子[100]面
SiO2(24個)の電荷密度の様子[110]面
分子動力学
• 電子状態の計算を経験的ポテンシャル関
数で近似し、原子系のダイナミクス・統計
集合を解く方法
• ポテンシャルの精度が課題
分子動力学ポテンシャル
1. 二体ポテンシャル
2. イオンポテンシャル
1. 固定電荷型
2. 電荷移動型
3. 多体ポテンシャル（共有結合）
1. クラスターポテンシャル
2. ボンドオーダーポテンシャル
4. Tight-Binding法
1．二体ポテンシャル１
• Lennard-Jonesポテンシャル
– 原子の電荷の偏りによる双極子間の引力相
互作用（van der Waals力）
– 希ガス（Ar, Ne）などに適用
 
 r

C1
 12

C2
6
r

r

repulsive
attractive
C 1  4 
12
, C 2  4 
6
1．二体ポテンシャル２
• Morseポテンシャル
– クラスタの結合のために提案
– 金属等に応用
– Expの関数系は後にEAMに使われる（第一原
理計算と良く合う～Erolessi(1994)）
   B exp  2 A ( r
 r
r


 r0



 r0 )  2 exp  A ( r

 r0 )

でφ=-B (結合エネルギ）
引力項∝h , 斥力項∝|h|2 となるWolfsberg-Helmholtz近似
を満たしている。
1、二体ポテンシャルの問題点
• 弾性定数の記述が貧弱
– 必ずCauchyの関係式C12=C44が成り立つ
• 空孔生成エネルギや表面などback-bondの記述
ができない
• 環境に依存したポテンシャルが必要！！
分子動力学ポテンシャル
1. 二体ポテンシャル
2. イオンポテンシャル
1. 固定電荷型
2. 電荷移動型
3. 多体ポテンシャル（共有結合）
1. クラスターポテンシャル
2. ボンドオーダーポテンシャル
4. Tight-Binding法
2. イオンポテンシャル（固定電
荷）
• BHSポテンシャル（SiO2など）
– 長距離ポテンシャルのため、カットオフ（スク
リーニング）もしくはEwald和が必要
 (r



exp(
r
 i  j)


)
 f

r
bi  b j





   

QiQ j
Coulomb
repulsive
C
6
r

attractive
 covalent  
Qは実効電荷、σはイオン半径、bはイオンの堅さ
長距離力の計算により、計算時間は10倍程度になる
2. Charge-Variable MD
• SiO2, Si3N4の電荷の移動
– バルクに対しては固定電荷で良いが、界面・
表面に対しては、局所的な電荷の移動を考慮
する必要がある。
– 電荷の移動＋原子間力の二つを解く必要が
ある。
Vashishita, 保川等
– 計算時間100倍程度
2, Charge-Variable MD
E es 

Ei (qi ) 
i
1
V

2
ij
( rij ; q i , q j )　静電エネルギ
i j
E i (qi )  E i (0)   i qi 
0
1
2
0
2
J i q i i 
0
E
: 電気陰性度
qi
 E
2
Ji 
0
qiq j
V ij ( rij ; q i , q j ) 
：　原子間の静電相互
rij
静電エネルギを電荷一
 収束計算のため計算
定の拘束条件
qi
2
: ( atomic hardness )
作用エネルギ


  q i  0 をつけて最小化。
 i

時間が 2 桁程度上がる
分子動力学ポテンシャル
1. 二体ポテンシャル
2. イオンポテンシャル
1. 固定電荷型
2. 電荷移動型
3. 多体ポテンシャル（共有結合）
1. クラスターポテンシャル
2. ボンドオーダーポテンシャル
4. Tight-Binding法
共有結合性材料の特徴
• 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
3．Si（C）系の多体ポテンシャル
• Stillinger-Waberポテンシャル
– “多体クラスターポテンシャル”
– ダイヤモンド構造と液体構造にフィッティング
• Tersoffポテンシャル
– “ボンドオーダーポテンシャル”
– 多形態（ダイマ・グラファイト・ダイヤモンド・単純立方・
体心立方など）のシリコンの構造にフィッティング
– 発展形としてBrenner(C-H)、Marty(Si-H)、大平(Si-H)
ポテンシャル
3-1. 多体クラスターポテンシャル
V 
V (R , R )  V (R , R , R )  V (R , R , R , R )  

                      

2b

2 body

3b

3  body


4b




4  body
ある特定の構造について、結合距離（二体項）、結合
角（三体項）、二面角（四体項）を合わせこむ。
合わせ込んだ構造に以外に遷移
する場合（化学反応・相転移等）
に対応できない！
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ポテンシャルの物性値１
テキスト 表5
Diamond struct.
Exp.
SW
T3
Lattice constants
5.429[Å]
5.431
5.432
Cohesive energy
4.63[eV]
4.63
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
(4-5) [eV]
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)
3-2. ボンドオーダーポテンシャル
V 
  A exp  
r
A

 b

B exp   B r



b

: bond  order
Abell (1985 )ポテンシャル
見かけは二体ポテンシャルbαβに多体
項が含まれる。
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(θ)を入れ、シリコンの物性値にフィッティングした。
Abellポテンシャル
bond  order : b   Z 
1 / 2
: Z は配位数
現在提案されてる、EAMポテンシャル・FS
ポテンシャル・Tersoffポテンシャルは、す
べて上の式で表される二次モーメント近
似によるポテンシャルである。
→二次モーメント近似の導出(資料参照)
角度依存性・配位数依存性
126°
角度依存項
Bond-order
g(θ)
• 120°付近でg(θ)が極小→ボンドオーダー増加
• 配位数が増大→ボンドオーダー減少
配位数増加
ボンドオーダの角度・配位数依存性
Tersoffポテンシャルの合わせこみ
• 異なる結晶構造（結合状態）のエネルギ・格子間距
離
ダイマ→グラファイト→ダイヤモンド→単純立方→BCC→FCC
• 異なる二つのバージョンが存在
– 表面構造合わせこみ（T2）→あまり使われていない
– 弾性定数合わせこみ（T3)
Tersoffポテンシャルの物性値１
テキスト 表5
Diamond struct.
Exp.
SW
T3
Lattice constants
5.429[Å]
5.431
5.432
Cohesive energy
4.63[eV]
4.63
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
(4-5) [eV]
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ポテンシャ
ルの応用へ
FSとTersoffポテンシャルの類似性


1
/
2




1
pair


 
Vi    Vi
( r )  B   exp   B r

 
2      



       
 Repulsive
attractive







1
pair


   Vi
( r )  B  exp   B r / 2
2   
 



b



 1 





,





1 / 2





 1   exp   B r  r




(  ,  )

  


     



b
:bond  order








g ( ) exp p ( r
 r )

       

angle
bond - length

    


すべての近接について

の和  配位数 ( Z )
Tersoff の   1,   1 / 2 , g ( )  1, p    B , q  1に対応
q
金属系のボンドオーダーポテン
シャル
※金属結合は不飽和な共有結合
• Tersoffポテンシャル
– 結合角・距離（弱い）と配位数でボンドオーダーを決定
• FSポテンシャル
– 結合距離と配位数でボンドオーダーを決定
• EAMポテンシャル
– 結合距離と配位数でボンドオーダーを決定
– ただし、埋め込み関数（ボンドオーダーの効果に相当
する項）を多項式やlog形にして汎用性・柔軟性を持た
せる → FSと傾向は同じ F (  )  
分子動力学ポテンシャル
1. 二体ポテンシャル
2. イオンポテンシャル
1. 固定電荷型
2. 電荷移動型
3. 多体ポテンシャル（共有結合）
1. クラスターポテンシャル
2. ボンドオーダーポテンシャル
4. Tight-Binding法
4. 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βの対角化（固有値問題）を
行い、固有値（エネルギ）を求める。Order N3
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ポテンシャルの基盤となる考え方。
計算をはじめるにあたって
• 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
「局所構造を反映した原子間ポテンシャル作成手法の調査」 ）
http://www.fml.t.u-tokyo.ac.jp/wg-cmd/
ポテンシャルの開発の試み
角度／距離依存性・配位数依存性でどこま
で高精度なポテンシャルが出来るか？
1. シリコンクラスタ用Tersoffポテンシャル
–
T2，T3→T4の開発
2. SiH4のシリコン表面反応ポテンシャル
–
–
–
ボンドオーダーポテンシャルに基づいた形式
ポテンシャルの合わせこみデータの収集
Saddle-Pointの合わせこみ
Tersoffポテンシャルのクラスタへの応用
• T2, T3ともにクラスターの構造が記述できない
– T3は角度依存性が大きく、表面に向かない
– T2は金属的結合で角度依存性が小さい
Tersoff3の結晶成長計算
表面原子が109°の角度（ダイヤモンド構
造の位置）で付着
クラスタ用ポテンシャルの作成
•
クラスターの構造間エネルギ差
Si3～Si6までのクラスタの形状間エネルギ差をエ
ネルギ構造差定理を使うことによりfitting。
• バルク（グラファイト・ダイヤモンド・単純立方・
BCC）の格子間距離・エネルギを最低条件
→ ボンドオーダーのみを調節して、T3より角度依存
性・配井数依存性を弱くして、T2の性質に近づけ
る。
どの構造がより安定化？を重視
結果（クラスターエネルギ）
Tersoffと比較して、エネルギの
順序の再現が大幅に改善！
考察（Global Min. ）
●S1b
Si4, Si5の最安定構造はAb-initioと一致する。
Si6は三角柱が最安定となり、Ab-initioの八
面体とは異なる。Si10もよりOpenな五角柱と
なり、Ab-initioと傾向が異なる。
●S1c
Si6, Si10の最安定構造の傾向があうポテン
シャルを作成した。しかしながら、バルクの
性質との共存は難しく。バルクの性質を犠
牲にせざるを得なかった。具体的には、ダイ
ヤモンド構造が最安定にならず、BCCが最
安定となってしまった。
結論（クラスタ用ポテンシャル）
• ボンドオーダーを調整することによって、ク
ラスタについては、T2・T3より優れたポテン
シャルが作成できた。
• しかしながら、バルクとクラスタの性質の双
方を満足するポテンシャルは角度と配位数
の調節だけでは不可能であった。
• シリコンクラスタはダングリングボンドの数
が多く、古典MDでは難しい。
実用的なCVD表面反応用Si-Hポテン
シャルの開発
• ポテンシャル作成における課題
1. どのような形状（安定or鞍点？）をあわせ
込むのか？
2. フィッティングデータをどうやって集めるの
か？
3. どのような関数形状を選ぶのか？
4. フィッティングアルゴリズムをどうするの
か？
１、どのような形状（安定or鞍点？）をあ
わせ込むのか？
–
–
–
クラスタ→バルクへの成長を表現するため、
様々な大きさのクラスタを合わせこみ対象に
選定
汎用性を高めるために、同一の大きさでも、
異なる形状のクラスタを合わせこみ対象に
選定
化学反応過程を考え、表面反応の素過程に
重要な現象のエネルギ障壁を合わせる
合わせ込んだクラスタ
合わせ込んだ反応過程
• 平衡値だけの合わせこみで良いのか？
– 結晶成長の表面反応をシミュレーションするた
めには、Saddle-Pointのエネルギを合わせる
必要がある。
※結晶成長に重要な素過程
Hの脱離
SiH4の
解離吸着
２、フィッティングデータをどうやって集
めるのか？
• 実験値の収集は大変・実験によりバラツキ
特に0Kでの結合エネルギの値の収集は困難
– 実験値は用いず、すべてGAUSSIAN98のG2
セットの計算結果を用いた（実験値との平均
誤差1.5％）
G2の計算
レシピ
１
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)
以降エネルギ計算（補正）
３、どのような関数形状を選ぶのか？
V  AF 1
1
bij 
2
N Si 
Si  H
b

Si  H
( N Si ) exp(   A rij )  BF 2
Si  Si
 b ji  G ij
ij
f c ( rij )

f c ( rik )
Si
Si  H 結合の Si の配位数
i 原子 ( Si )の配位数
k
Nj 

f c ( r jl )

( N i , N j , N ij ) : bond  order
jH
Ni 
( N Si ) exp(   B rij ) bij
Tersoff３をベースにし
たEBOPを関数形状
に選定
j 原子 ( Si )の配位数
l
N ij 
Si

f c ( rik ) 
k  Si
※ F1
Si  H

f c ( r jl )
N iと N jの中の Si の数
l  Si
Si  H
, F2
は Si  H 結合のみ、
Si  Si
G ij
は Si  Si 二重結合のみの補正
４、フィッティングアルゴリズムをどうす
るのか？
1.
2.
3.
4.
5.
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-SiH, H-Si-Hの三体項を調整することによって合わせこむ
１、はcrude MCによって、最適値を見つける。２，３，４は１
の中からtry and errorで合わせこむ。５は４で設定していな
い三体項を使って合わせこむ。
6.0%
5.0%
4.0%
deviation
Fitting結果（誤差）
10.0%
9.0%
si7h14
si7h14a
si7h16
si6h6
si6h8
si6h8b
si6h8a
si6h10
si6h10b
si6h10a
si6h12
si6h12b
si6h12a
si6h14
si5h6
si5h8
si5h8b
si5h8a
si5h10
si5h10b
si5h10a
si5h12
si4h6b
si4h6a
si4h8
si4h8b
si4h8a
si4h10
si3h4b
si3h4a
si3h6
si3h6b
si3h6a
si3h8
Hsisih3
si2h6
si2h4
si2h2
sih4
sih3
sih2
sih
0.0%
Marty
This study
8.0%
7.0%
3.0%
2.0%
1.0%
誤差は、ほとんどが2%以下、平均で1.1％(Martyは3.3%)
環状のクラスタは苦手（Tersoff3の特性）
Fitting結果 (Hの解離１)
d1
eV
d2
H解離のポテンシャルマップ
Eb=2.37eV (DFT:2.49)
Si-Hのカットオフ距離と
H-H-Si, H-Si-Hの三体項
Fitting結果 (SiH4の解離吸着)
Eb=+0.56eV (ab-inito 0.5～0.6eV)
Eb
Saddlepoint!
-0.70eV
Si-Siのカットオフ距離と
H-Si-Hの三体項
結論
SiH4の表面反応に使用する分子動力学ポテンシャルを作成
1. どのような形状（安定or鞍点？）をあわせ込むの
か？
→クラスタ形状と表面反応の活性化エネルギ
2. フィッティングデータをどうやって集めるのか？
→GAUSSIAN９８による計算データ
3. どのような関数形状を選ぶのか？
→EBOPに補正項を入れる。パラメータが多ければ
あわせこみの可能性も増える
4. フィッティングアルゴリズムをどうするのか？
→Crude MC+try and error
ご清聴ありがとうございました
注意！（Cutoff関数の問題点）
x
H
Si
Si
4.0Å
Si-Hのcutoff 2.4Å
Cutoffにより、短距離
の間に大きなエネルギ
変化が生じる！！
新しい結合
が出来ると、
急激にエネ
ルギ低下＝ｚ
局所的に強
い引力が働
く！！
結合を切ると
きは逆に局
所的に大き
な障壁
Saddle-point energy計算
• 双方向Gradual Ascent法（GA法）
Bulatovらが開発したGA法（初期状態と最終状態
をバネで結び、剛性を除々に強くすることによって、
saddle-pointを見つける）を双方向に拡張
④
初期状態
GA 1
①
GA 3
③
⑤
GA 4
②
GA 2
GAバネ
初期状態・最終状態ともに、少しずつGAで近
づけてsaddle-pointを探す
最終状態
Fitting手順 (Hの解離2)
d1
eV
d2
H解離のポテンシャルマップ
Eb=2.9eV (DFT:2.9)
Si-Hのカットオフ距離と
H-H-Si, H-Si-Hの三体項
Tight-Bindingの解法
• n番目の固有状態ψ(n)の波動関数を原子軌道φiα（原子i,
軌道α）の線形結合で表す
分子軌道　: 
(n)

 C i  i
(n)
i
ˆ
H
• 波動関数に代入して永年方程式を得る
N

H i  , j  C j  
(n)
(n)
(n)
C i
H
i , j
j
(n)
*
   i  Hˆ  j  dV


(n)

• ハミルトニアンマトリックスHiα,jβの対角化（固有値問題）を
行い、固有値（エネルギ）を求める。
E band 

i , j
z
2  C i C j H i , j
(n)
n 1
( n )*
(n)
状態密度と結合エネルギー
原子
固体
Energy
Hˆ 
(n)


(n)
エネルギの状
態密度分布
(n)
波動関数
ntotal(E)
原子軌道　:  i
分子軌道　: 
(n)

C 
(n)
i
i
i
結合は、原子状態より低いエネルギの分子軌道が作られ
E
ることによって生じる。
U 
En
( E ) dE

f
→ 状態密度の分布形状が結合エネルギを決める。
total
状態密度のモーメント
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
4
 2 W  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

H
i , j
H
j  ,i 
j
2-Hops 二次
モーメント
 
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
局所構造を反映した
原子間ポテンシャル作成手法の調査
主旨
材料を構成する原子構造体の外力に対する応答として，き裂先
端場から射出される転位挙動や結晶粒界近傍に生じる種々の
力学現象などが，分子動力学法により正確に描写できるように
なった．しかしながら，非可逆な変形応答やその極限状態 とし
ての破壊事象は，本質的に不均質な構造を対象とした現象で
あり，これまでの平均的で均質系を仮定したマクロ系の材料特
性のみを表現しうる従来の原子間相互作用の記述では，本質
的に捉えられない． そこで，本WGでは従来からあるポテンシャ
ルの問題点などを整理するとともに，不均 質な局所構造の影
響を反映した原子間ポテンシャル作成のための調査を目的と
する．
活動内容
1.
現状の経験的，半経験的ポテンシャルの問題点の整理
Tersoff型、 MEAM(modified EAM)型、 Tight-binding法
2.
非平衡現象に対応するポテンシャルのフィッテング手法の開発
汎用第一原理ＭＤの利用技術の確立
ポテンシャルの第一原理ＭＤへの繰み込み手法
3.
ポテンシャルのデータベース構築の検討
データベースの作成
大学／企業で必要と考えられるポテンシャルの調査
活動方法
1. 年4回程度のミーティングを実施し，あらかじめ決め
られた報告者が，事前にE-mailで議事内容を連絡す
るとともに当日会議の進行を行う．
2. 報告の内容は，設定された調査の経過報告や文献紹
介など，特に制約は設けない．
3. 学会のオーガナイズドセッションやワークショップ，
フォーラムなどを利用した 活動結果の公表を定期
的に行う．