Document

シリカガラスの分子動力学
シミュレーション
福井大学 工学部
葛生 伸
シリカガラスの分子動力学シミュレーション
体積 (密度) の温度依存性
表面の構造
体積の温度依存性
R.Brückner, J.Non-Cryst.Solids 5, 123, 1970
一般的のガラス
シリカガラス
以前の研究
D. R. Perchak and J. M. O’Reilly, J.Non-Cryst.Solids 167, 211 (1994)
ポテンシャル
B.P.Feuston, S.H.Garofalini, J.Chem.Phys. 89, 9, 1988
 

 
U (r1, r2 ,, rn )   u2 (ri , r j ) 
i, j
2
 Zi Z j e
u2 (ri ,r j )  
 rij
u3 (ri ,r j ,rk ) 
  
 u3 (ri , r j , rk )  
i , j ,k

 rij 
rij 

 erfc    Aij exp   
  


L


 ij 

h(rij , rik , jik )  h(r jk , r ji , ijk )  h(rki , rkj , ikj )
h(rij , rik , jik ) 
 i


c 2
i

i exp

(cos


cos

jik
jik )
C
C 
 rij  r i rik  r i   r  r C and r  r C 
ik
i 
 ij i
その他の場合, 0 


u2中のパラメタ
Si – Si
Si – O
O–O
IP
AP
IP
AP
ρij /Å
0.29
0.29
0.29
0.29
1.88×10-9
1.88×10-9
2.96×10-9
3.00×10-9
βij /Å
2.50
2.29
2.50
2.34
IP
AP
0.29
0.29
0.725×10-9
1.10×10-9
2.50
2.34
O – Si – O
Si – O – Si
Aij /erg
u3中のパラメタ
ric /Å
λi /erg
3.0
19×10-11
2.6
0.3×10-11
IP: 等方ポテンシャル (u2のみ)
AP: 異方性ポテンシャル (u2+u3)
ri /Å
2.8
2.0
温 度 履 歴
1 step = 0.1 ns
Temperature (K)
8000
6000
4000
2000
0
1
2
Time (ns)
3
Si 粒子数
216
O 粒子数
432
Δt
10-15s
4
各条件での密度の温度依存性
4
IP 10 10 Pa
3
)
5
Density ( g/cm
3
IP 10 Pa
2
AP 10
10
Pa
1
AP 10 5 Pa
0
0
2000
4000
6000
Temperature ( K )
8000
密度の温度履歴
2.50
Density / g cm
-3
Isotropic
105 Pa
2.00
1.50
0
Cooling
Heating
2000 4000 6000 8000 10000
Temperature/K
結合角分布の温度依存性 (常圧等方性ポテンシャル)
105 Pa
Isotropic
Isotropic
105 Pa
9000K
8000K
50
100
9000K
8000K
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
2000K
150
O-Si-O Bond Angle / degree
80
100
120
140
160
Si-O-Si Bond Angle / degree
180
結合角分布の温度依存性 (加圧異方性ポテンシャル)
Anistropic
1010 Pa
Anisotropic
10
50
100
10
9000K
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
2000K
150
O-Si-O Bond Angle / degree
80
Pa
100
120
140
160
Si-O-Si Bond Angle / degree
180
結合角分布の温度依存性 (常圧異方性ポテンシャル)
Anisotropoc
Abisotropic
105 Pa
105 Pa
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
50
100
9000K
150
O-Si-O Bond Angle / degree
2000K
80
100
120
140
160
Si-O-Si Bond Angle / degree
180
結合角分布の温度依存性 (加圧等方性ポテンシャル)
Isotropic
1010 Pa
Isotropic
9000K
9000K
8000K
7000K
6000K
1010 Pa
8000K
7000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
2000K
50
100
150
O-Si-O Bond Angle / degree
80
100
120
140
160
Si-O-Si Bond Angle / degree
180
配位数分布の温度依存性
1.0
Isotropic
Si-(4O)
1.0
105 Pa
0.5
Si-(4O)
0.5
Si-(3O)
Si-(3O)
Si-(5O)
Si-(6O)
Fraction
Fraction
Anisotropic
1010 Pa
0
1.0
Si-(5O)
0
1.0
O-(2Si)
0.5
O-(2Si)
0.5
O-(1Si)
O-(1Si)
O-(3Si)
0
O-(3Si)
2000
4000
6000
Temperature / K
8000
0
2000
4000
6000
Temperature / K
8000
配位数分布の温度依存性
1.0
Anisotropic
1.0 Isotropic
105 Pa
1010 Pa
Si-(4O)
Si-(4O)
0.5
0.5
Si-(5O)
Si-(3O)
Si-(3O)
Si-(5O)
0
Fraction
Fraction
Si-(6O)
1.0
0
1.0
O-(2Si)
O-(2Si)
0.5
0.5
O-(3Si)
O-(1Si)
O-(3Si)
0
2000
4000
6000
Temperature / K
8000
O-(1Si)
0
2000
4000
6000
Temperature / K
8000
密度の仮想温度依存性温度依存性
各温度から300 Kまで急冷前後密度依存性
Density / g cm
-3
3
2
IP 105 Pa
1
0
0
10
IP 10
Pa
Before Quenching
Quenched
Before Quenching
Quenched
2000 4000 6000 8000 10000
Fictive Temperature / K
結合角分布の仮想温度依存性 (常圧等方性ポテンシャル)
Isotropic
Isotropic
105 Pa
105 Pa
TF
TF
9000K
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
50
100
5000K
4000K
4000K
3000K
3000K
2000K
2000K
150
O-Si-O Bond Angle / degree
80
100
120
140
160
Si-O-Si Bond Angle / degree
180
結合角分布の仮想温度依存性 (加圧異方性ポテンシャル)
Anisotropic
Anisotropic
1010 Pa
1010 Pa
TF
9000K
TF
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
2000K
50
100
150
O-Si-O Bond Angle / degree
3000K
2000K
80
100
120
140
160
Si-O-Si Bond Angle / degree
180
配位数分布の仮想温度依存性
1.0
1.0
Si-(4O)
Anisotropic
Isotropic
1010 Pa
Si-(4O)
0.5
0.5
Si-(5O)
Si-(3O)
0
1.0
Si-(6O)
Si-(5O)
Fraction
Fraction
1010 Pa
Si-(3O)
0
1.0
O-(2Si)
O-(2Si)
0.5
0.5
O-(3Si)
O-(3Si)
O-(1Si)
0
O-(1Si)
0
2000
4000
6000
Fictive Temperature / K
8000
2000
4000
6000
Fictive Temperature / K
8000
シリカガラス表面の分子動力学シミュレーション
表面では欠陥構造が多数
⇒ 原子間の電荷の移動
⇒ 電荷平衡法により考慮
シリカ表面の分子動力学シミュレーション
ポテンシャル
Morse-Stretch ポテンシャル + クーロンポテンシャル
  
UCoulomb rij 
 

UMS rij  D0 e
R0 (Å)
ij
Qi Q j
rij
 1rij / R0 
 2e
 / 21rij / R0 
D0 (kcal/mol)

γ
O-O
3.7835
0.5363
10.4112
Si-Si
3.4103
0.2956
11.7139
Si-O
1.6148
45.9970
8.8022
原子間の電荷移動効果
電荷平衡 (QEq)法
原子 i の化学ポテンシャル
 i Q1,,QN    i0 
J r Q
ij
j
平衡条件
ij
j
J ij : クーロン積分
 i0 : 定 数
1   2     i     N
電荷の保存
N
Q
i
 Qtot
i 1
基礎方程式
CQ  D
i0 (eV)
O
8.741
Si
4.168
電荷平衡 (QEq) 法の基礎方程式
1


 J 21  J11



J  J
 N1 11

 J 2N

 J NN
1
J 22  J12

J N 2  J12
1
 Q1   q tot 

  0

0
 J1N  Q1   1   2 









 

0
0
 J1N  QN   1   N 
クーロン積分
J ij r   d r i d r j ni ri 

3
3
2
 
1
   n j r j
ri  r j  r
ni ri   Ai r ni 1e i r
Ai: 規格化因子,ni: 主量子数,zi: 定数
2
クーロン積分
J ij r   d r i d r j ni ri 

J(r) / eV
O-O
3
3
14..4 / r
O-Si
Si-Si
r/Å
2
 
1
   n j r j
ri  r j  r
2
シミュレーションセルと周期境界条件
z
x, y
粒子数 648
温度履歴
Surface Generation
Stepwise Cooling
100 K / 1 ps
Relaxation
Stepwise Cooling
100 K / 1 ps
Relaxation
Analysis
SiおよびOの相対数密度分布 (バルクの総数密度で規格化)
qEQ
FQ
O
n
n
O
Si
Si
z/Å
z/Å
密度分布
Present
z/Å
表面付近の動径分布関数
FQ
g(r)
g(r)
QEq
z/Å
z/Å
表面付近のSi-O-Si結合角分布
Tetrahedral
qEQ
Planer
FQ
配位数分布
4-coord. Si
BO
FQ
FQ
qEQ
qEQ
qEQ
qEQ
FQ
3-coord. Si
z/Å
Deb
NBO
FQ
z/Å
配位数分布
Si
4-coord. Si
charge / e
3-coord. Si
Deb
Total Charge density
NBO
BO
z/Å
qEq
O
Si
Debye –Waller Factor / Å2
Debye –Waller Factor / Å2
表面付近のSi-O-Si結合角分布
FQ
O
Si
z/Å
Deb
Debye-Waller
Factor
z/Å
1
lim
t  N
Si
1
lim
t  N
O



ri t   ri 0
all Si

all O


ri t   ri 0
表面の欠陥構造
E’ センター