Transcript 7. シリカガラスの分子動力学シミュレーション

シリカガラスの分子動力学
シミュレーション
福井大学 工学部
葛生 伸
シリカガラスの分子動力学シミュレーション
体積 （密度) の温度依存性
表面の構造
体積の温度依存性
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 , r 2 , , r n ) 
 
 u 2 ( ri , r j ) 
i, j
 ZiZ je 2
u 2 ( ri , r j )  
 r
ij

u 3 ( ri , r j , rk ) 
  
 u 3 ( ri , r j , rk )  
i , j ,k

 r ij 
 r ij 
 erfc 

  A ij exp  
  


L


ij 


h ( r ij , r ik , jik )  h ( r jk , r ji , ijk )  h ( r ki , r kj , ikj )
h ( r ij , r ik , jik ) 
 
i
i

 i exp

C
 r  r i r  rCi
ik
 ij

 (cos  jik  cos  jik c ) 2


C
C


r  r i and r ik  r i
 ij



その他の場合，
０


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
D e n sity ( g /cm
3
IP 10
2
5
Pa
10
AP 10
Pa
1
AP 10
0
0
2000
4000
6000
Temperature ( K )
5
Pa
8000
密度の温度履歴
2.50
D e n s ity / g c m
-3
Isotropic
105 Pa
2.00
1.50
0
Cooling
Heating
2000
4000
6000
Temperature/K
8000
10000
結合角分布の温度依存性 （常圧等方性ポテンシャル)
5
Isotropic
Isotropic
10 Pa
5
10 Pa
9000K
9000K
8000K
8000K
50
100
O-Si-O Bond Angle / degree
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
2000K
150
80
100
120
140
Si-O-Si Bond Angle / degree
160
180
結合角分布の温度依存性 （加圧異方性ポテンシャル)
Anistropic
10
10
Pa
Anisotropic
10
50
100
O-Si-O Bond Angle / degree
9000K
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
2000K
150
80
10
Pa
100
120
140
Si-O-Si Bond Angle / degree
160
180
結合角分布の温度依存性 （常圧異方性ポテンシャル)
Anisotropoc
Abisotropic
5
5
10 Pa
10 Pa
9000K
8000K
7000K
5000K
5000K
2000K
O-Si-O Bond Angle / degree
7000K
6000K
3000K
100
8000K
6000K
4000K
50
9000K
150
4000K
3000K
2000K
80
100
120
140
Si-O-Si Bond Angle / degree
160
180
結合角分布の温度依存性 （加圧等方性ポテンシャル)
Isotropic
10
10
Pa
9000K
10
Isotropic
10
Pa
9000K
8000K
8000K
7000K
6000K
7000K
6000K
5000K
5000K
4000K
4000K
3000K
3000K
2000K
2000K
50
100
O-Si-O Bond Angle / degree
150
80
100
120
140
Si-O-Si Bond Angle / degree
160
180
配位数分布の温度依存性
1.0
Anisotropic
10
10 Pa
1.0
Isotropic
5
Si-(4O)
10 Pa
Si-(4O)
0.5
0.5
Si-(3O)
Si-(3O)
Si-(5O)
F ra c tio n
F ra c tio n
Si-(5O)
Si-(6O)
0
1.0
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
5
10 Pa
10
10
Pa
Si-(4O)
Si-(4O)
0.5
0.5
Si-(5O)
Si-(3O)
Si-(3O)
Si-(6O)
F ra c tio n
F ra c tio n
Si-(5O)
0
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まで急冷前後密度依存性
D e n s ity / g c m
-3
3
2
IP 10
1
0
0
IP 10
5
Pa
10
Pa
Before Quenching
Quenched
Before Quenching
Quenched
2000
4000
6000
8000
Fictive Temperature / K
10000
結合角分布の仮想温度依存性 （常圧等方性ポテンシャル)
Isotropic
Isotropic
5
5
10 Pa
10 Pa
TF
TF
9000K
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
4000K
3000K
2000K
50
100
O-Si-O Bond Angle / degree
150
5000K
4000K
3000K
2000K
80
100
120
140
Si-O-Si Bond Angle / degree
160
180
結合角分布の仮想温度依存性 （加圧異方性ポテンシャル)
Anisotropic
10
10
Anisotropic
Pa
10
10
Pa
TF
9000K
TF
9000K
8000K
8000K
7000K
7000K
6000K
6000K
5000K
5000K
4000K
4000K
3000K
2000K
50
100
O-Si-O Bond Angle / degree
150
3000K
2000K
80
100
120
140
Si-O-Si Bond Angle / degree
160
180
配位数分布の仮想温度依存性
1.0
1.0
Si-(4O)
Isotropic
10
10
Anisotropic
Si-(4O)
Pa
10
0.5
Pa
0.5
Si-(3O)
0
1.0
Si-(6O)
Si-(5O)
F ra c tio n
Si-(5O)
F ra c tio n
10
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 ポテンシャル + クーロンポテンシャル
U Coulomb
r ij   
ij
 

U MS r ij  D 0 e
R0 (Å)
Q iQ j
r ij
 1 r ij / R 0 
 2e


/ 2  1 r ij / R 0
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 Q 1 , , Q N    i 
0
 J r Q
ij
j
平衡条件
ij
j
J ij : クーロン積分
 i : 定　数
0
1   2     i     N
電荷の保存
N
Q
i
 Q tot
i 1
基礎方程式
CQ  D
i0 (eV)
O
8.741
Si
4.168
電荷平衡 (QEq) 法の基礎方程式
1


 J 21  J 11



J  J
11
 N1
1

1
J 22  J 12

J 2 N  J 1N



J N 2  J 12

J NN  J 1N
 Q 1

 Q 1



 Q
 N
  q tot
  0
0



  1
2


 
 0  0
N
  1
クーロン積分
J ij r  
 ni
d
3
r i  
r i d r j  n r i 
3
i
Ai r
n i 1
e
2
 
1


 nj r j
ri  r j  r
 i r
Ai: 規格化因子，ni: 主量子数，zi: 定数
2







クーロン積分
J ij r  
J(r) / eV
O-O
d
3
r i d r j  n r i 
3
i
14..4 / r
O-Si
Si-Si
r/Å
2
 
1


 nj r j
ri  r j  r
2
シミュレーションセルと周期境界条件
z
x, y
粒子数 ６４８
温度履歴
Surface Generation
Stepwise Cooling
100 K / 1 ps
Relaxation
Stepwise Cooling
100 K / 1 ps
Relaxation
Analysis
Ｓｉおよび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
FQ
Debye –Waller Factor / Å2
Debye –Waller Factor / Å2
表面付近のSi-O-Si結合角分布
O
Si
z/Å
z/Å
Deb
Debye-Waller
Factor
lim
t
lim
t
1
N Si
1
NO



r i t   r i 0 
all Si

all O


r i t   r i 0 
表面の欠陥構造
E’ センター