分子動力学法によるケイ酸塩系ガラスの構造シミュレーション
環境無機化学講座 12-127 鈴木 智恵【緒言】
ガラス材料はその透明性、光学特性、耐久性、加工性などの特性を利用して日用品から工業製品まで広く利用されている材料である。しかしながらガラスの原子レベルの構造解析は、ガラスに特有の不規則性のため実験的手法では困難である。そこで実験手法とは異なる別のアプローチとして計算機シミュレーションによる構造解析も研究が進められている。本研究で使用した分子動力学(Molecular Dynamic : MD) 法もガラス材料の構造解析に適したシミュレーション技法の一つである。
MD シミュレーションは Newton の運動方程式を利用したシミュレーション法で、ガラスの溶融、急冷時における原子の挙動を経時変化として観察できる。その結果を利用して動径分布関数を始め、IR や Raman など分光学的性質、熱力学的性質などが求められる。
本研究では代表的な実用ガラスであるケイ酸塩系ガラスをテーマにシミュレーションを行った。SiO2100%のシ
リカガラス中では Si は O を正四面体型に四配位し、O は Si を二配位している。O-Si-O 結合角の分布は 109°付近に集中しているのに対し、SiO4 ユニットをつなぐ Si-O-Si 結合角の分布は 147°を中心に 120°‐180°と広く分布していることが NMR や X 線回折実験の結果から報告されている。シリカガラスに Na2O を添加すると、次の模式図で示すように Si の網目構造を形成している架橋酸素（BO）が切断され非架橋酸素（NBO）が生成する。

Na≡Si－O－Si≡ ＋ Na2O→ ≡Si－O(＊)

O(＊)－Si≡Na

BO と NBO は同じ O 原子ではあるが、配位原子の相違から異なNBOる電荷を持っており、分極率などの化学的性質も異なることが知られている。Fig.1 に XPS による Na2O-SiO2 ガラスの O1s スペクトル密度を反映している。このスペクトルから、配位状態の異なる BOと NBO は異なる電子密度を持つことが分かる。また Na2O の増加とともに BO 及び NBO のピークが低エネルギー側にシフトしていることから、BO、NBO ともに電子密度が増加したことが分かる。これ
らの配位状態や電荷の相違は原子間の結合に大きな影響を与え

Na Auger

Intensity / Arb. Unitsを示す。O1s など内殻軌道の束縛エネルギーはその原子の電子

x=0.50x=0.40x=0.33x=0.30x=0.20ているはずであり、原子間力から原子の挙動を追跡する MD シミ

x=0BOュレーションにもその影響を反映させる必要がある。
これまで MD シミュレーションをガラス材料に適用した研究例は数多く報告されており、尚且つ様々なポテンシャル関数が提案されている。しかしそのほとんどが同種原子には同じ電荷を与え、配位状態、原子間距離など原子の周辺環境を反映していなかっ540535530525

Fig.1 O1s photoelectron spectra ofxNa2O・(1-x)SiO2 glasses[2].

た。本研究で使用した MD プログラムではこれまで広く使用されてきた二原子間の相互作用を表現する二体ポテンシャルに加えて、三原子間の結合角を制御できる三体項を含むポテンシャル関数を導入した。三体項の角度制御によって間接的に配位数を制御している。さらに電荷移動の式を導入し、同種原子でも配位環境によって異なる電荷を与えることが可能となっている。これにより、原子の挙動に電荷の影響を反映させることを試みた。
-1-520Binding energy /eV【計算手法】(i) ポテンシャル関数
本研究で使用した MD プログラムは、ポテンシャル関数として Born-Mayer-Huggins 型の二体ポテンシャ
ル φ ij 及び Anderson[2]によって報告されている三体ポテンシャル関数 φ ijk を用いた。二体ポテンシャルは二原子間の相互作用を表現し、原子間距離を制御している。また結合角の項を含む三体ポテンシャルで結合角を制御し、間接的に配位数も制限している。
rij
Z i Z je 2
φ ij =+ B ij ( r ij )ji
4πε 0 r ij・・・(1)2 体項Zj Z exp{ (σ i + σ j − r ij ) ρ ij }B ij ( r ij ) = Aij 1 + i +Pairnin j 

3 体項{()}

φijk = (ϕ ij + ϕ ik ) exp − γ ijk θ − θ ijk 2・・・(2)ϕ ij = −C ij B ij ( ζ ij ) exp{(ζ ij − rij ) ηij }Zi ：i 番目の原子の電荷
ｒij ：原子 i, j 間の距離
Aij ：反発定数
ｎi ：原子 i の最外殻電子数riijikjσi ：原子 i に関するサイズ･パラメーター
ρiｊ ：i, j 間のソフトネス･パラメーター
θijk ：j, i, k 間の結合角
Cij,ζij,ηij,γijk, θ ：任意定数(ii) 電荷移動
電荷移動を表現する式として以下の(3), (4)式を導入した。
Z i = Z 0 + ∑ ∆q ij・・・(3)

rk

θ ijk
3-Body0.020∆q ij =D el

exp{b ( r ij − X )} + 1・・・(4)∆q：I, j 間で移動する電荷
r：結合距離b, Del , X ：任意のパラメーター

三体ポテンシャル項を計算する際に(3), (4)式を導入し、原子間を移動する電荷量を求める。(4)式から求める電荷の移動量∆q は 結合距離に応じて変化させることができる。また(3)式では配位している原子の数を反映しているので、配位数が異なる同一原子にも異なる電荷を与えることができる。(3), (4)式をTransferred charge ∆qj0.015

0.010

0.005

0.00000.10.20.3

0.40.5Distance /nm

Fig.2 Electronic charge transferred
between Si4+ and O2- in Si-O bond

導入することにより、各原子にそれぞれの周辺環境に応じた電荷を与えることを試みた。
Fig.2 は本研究で実際に用いた電荷移動の関数をプロットしたものである。ここで使用したパラメーター
は SiO2 のα-cristobalite 結晶における最適化計算から求めた。

【結果と考察】MD シミュレーションにより xNa2O・(1-x)SiO2 ガラス(x=0, 0.1, 0.2, 0.3, 0.4, 0.5)の構造モデルを作成したが、今回は x=0.0 及び 0.33 のモデルを中心に発表する。これらのモデルを中性子線動径分布関数で評価し、実測データと比較したものを Fig.3 に示した。x=0 のシリカガラスのモデル(Fig.3[a] ) を実測と比較すると、-2-ピークの位置及び面積、すなわち配位数もよく一致している。また Fig.3[b]の x=0.33(Na2O・2SiO2) ガラスでは近接した原子間距離に実測値と若干のずれはあるものの構造は概ね再現できていると言える。

-22.01.02.01.5
1.00.50.00.00.1

2.5O-Na
Na-Na
Total [Calc.]
Total [Obs.]

223.0

Si-Si
Si-O
Si-NaO-O[b]3.0
2 π r ρ(r) / barns nm-24.0 [a]2 π r ρ(r) / barns nm

3.5Si-Si
Si-O
O-O
Total [Calc.]
Total [Obs.]

0.20.3 0.4 0.5
Radius /nm

0.6

0.10.20.30.4

0.50.6

Raius /nmFig.3 Neutron radial distribution functions for [a] SiO2 and [b] Na2O ・ 2SiO2 glasses.
Experimental data were taken from Refs.3 and 4.

SiO 2[a]Fig.4 には MD シミュレーションから得られたモデル中の[b]

Na2O . 2SiO 2

SiO2Na2O . 2SiO 2Obs. [SiO 2]

O-Si-O 角の狭い分布が確認できる。また Fig.4[b]では Si-O-Siの角度分布を X 線回折実験からFrequency /-分 布 を 示 し た 。 Fig.4[a] で はFrequency /-O-Si-O 及び Si-O-Si 結合角の予想された結果[5]と比較した。
実測の結果よりもやや狭い分布を与えているが、分布の頂点を与える角度はほぼ同じであった。

90100110120Bond Angle /degree130

110 120 130 140 150 160 170 180

Bond Angle /degree

Fig.4 Bond angle distribution of [a] O-Si-O and [b] Si-O-Si.

Fig.5 は得られた構造モデルにおける O 原子の電荷を示している。Na2O の添加に伴い、NBO だけではな
く BO も負電荷が増加している。この傾向は XPS の結果とも一致している。また、MD シミュレーションから得られた電荷の値と併せて DV-Xαクラスター分子軌道計算(MO)により求めた同組成結晶の電荷[6]も示した。分子軌道計算の結果と比較すると、MD の結果の方が Na2O の添加に伴う負電荷の増加が大きい。これは今回のシミュレーションでは Na の電荷を固定した設定にしているためである。この Na イオンの電荷の差が構造モデルの再現性に影響しているかどうかは現時点では明確ではない。いずれにしても、これまでの MD シミュレーションでは反映するのが困難であった、配位状態ごとに異なる電荷を与えることに成功した。

-3-BO [MD]NBO [MD]
Average [MD]

-0.9

SiO2
Na2O .2SiO2Frequency /--1.0

BO [MO]NBO [MO]
Average [MO]Charge-1.1-1.2-1.3-1.4

2 3 4 5 6 7 8 9 10111213
Ring Size /-0.0 0.1 0.2 0.3 0.4 0.5
Na O Content, x / molar ratio
2Fig.5 Charge obtained from MD simulation andDV-Xα MO calculation.

Fig.6

The distribution of Ring Size.

ケイ酸塩系ガラスでは…-Si-O-Si-O-…と続く網目構造が形成されており、構造が Si と O 原子による環で構成されていると捉えることができる。ガラス構造は結合角分布の結果からも分かるように、結合角は一定ではないため、環の大きさを実験的手法で測定することができない。そこでシミュレーションから得られたモデルで環の大きさを観察した。Fig.6 には Ring Size を示したが、Ring Size とは環を形成している Si 原子の数を表している。SiO2 ガラスでは六員環が最も多く形成されている。六員環は SiO2 結晶の多形の中に多く見られる環の数であり、ガラス構造は中・長距離の範囲でも結晶構造と類似していることが示唆される。また
Na2O を添加すると環の分布はサイズの大きい方に広がっている。これは-Si-O-Si-の結合が Na イオンによって切断されたためである。

【結論】○ 電荷移動の項を導入した MD シミュレーションで、各原子に周辺環境を反映した電荷を与えることに成功した。
○ MD シミュレーションにより得られたガラス構造モデルは実測データを非常によく満足するものであった。

《参考文献》[1] 松本, 難波, 三浦, J. Ceram. Soc. Japan., 106(4) (1998) 415
[2]
[3]
[4]
[5]
[6]

D. C. Anderson and J. Kieffer, J. Chem. Phys., 98 (1993) 8979
D.I. Grimley, A.C. Wright and R.N. Sinclair, J. Non-Cryst. Solid., 119 (1990) 49
M. Misawa, D.L. Price and K. Suzuki, J. Non-Cryst. Solid., 37 (1980) 85
R.L. Mozzi and B.E. Warren, J. Appl. Crystallogr. , 2 (1969) 164
萩原達也, 岡山大学修士学位論文 (2001)-4-