On Molecular Dynamics Simulations of Silicate Melt/Glass−Molecular Simulations of Inorganic Compounds−

Abstract

無機化合物のMD計算を念頭において，高分子様の構造を持つ無機凝集体について構造などの概要を示し，古典分子動力学法による2元系ケイ酸塩溶融体/ガラスの計算における構造緩和の問題を述べる．また無機凝集体のMD計算の困難さ，意義，および必要性について述べる．

Translated Abstract

Structural features of inorganic polymer, oxoate compounds like silicate, are described. Then molecular dynamics simulations of sodium silicate melt / glasses were carried out using an empirical inter atomic potential model. The behavior of the energy and structural relaxation were investigated. The Q-species (the most important characters describing structure of silicate) change during an (NPT)-MD, show no obvious systematic change even in the clear energy relaxation. The difficulty of discussing the relation between energy-structure relaxation was shown. The necessity of computational methods of inorganic solids are discussed.

1 はじめに

有機分子集合体・高分子・生体分子などの分子シミュレーションは近年ますます盛んになっている．その様な物質系の重要性を考えると当然のことであろう．一方で，無機物質，特に固体，液体， 非晶質の分子シミュレーションによる研究は，学会発表の数で見る限り，増加しているとは言えないどころか，減少しているように思える．無機凝集体の分子シミュレーションを必要とする研究対象は種々多数あると考えられるにもかかわらずそのような状況にあると私は考えている．

このような研究が必要とされる物質系の1つが酸化物溶融体/ガラスである．実際にMD計算を行ってみると，以下に示すように大きな困難性があることが明らかとなってきた．

加えて，経験原子間モデルに基づく古典MD計算の意義についても考えておく必要がある．計算機の発展とともに非経験電子状態計算が多くの系で可能になってきている．それでも非常に大規模な系を用いないと解明できない物理量は多くある．近年メソスケール現象の観測が必要とされることが多くなってきている．古典分子シミュレーションはそのような要請に応えられる有力な手法の第一である．

2 ケイ酸塩溶融体/ガラスのMD計算

酸化物の中でケイ酸塩化合物は地球表層における固体構成物の大部分を占め，また工業的にも，原料や製品として重要であると考えられる．ケイ酸塩を主成分の1つとするガラスは高機能化などのため研究開発が続けられており，金属精錬においてもスラグとして重要な役割をはたしている．地球構成物質としてのマグマや鉱物は地球や惑星の構成・発達・進化や物質移動を解明するために不可欠な研究対象である．

ケイ酸塩の結晶固体・非晶質・ガラス・溶融体の特徴を述べるために，最も簡単な系としてA₂O-SiO₂あるいはBO-SiO₂系を考えることが適切であろう．Aはアルカリ金属元素でありBはアルカリ土類元素である．このような系では，化学式mA₂O·nSiO₂あるいはmBO·nSiO₂でm/n≦<2において，化学式中の，したがって構造中のすべての酸素原子がSiと結合していることになる．

Siは常圧下，室温から融点程度の温度においては酸素原子について4配位で，[SiO₄]^4-四面体を形成しているが，化学組成により，重合して≡Si-O-Si≡架橋して(この酸素原子を架橋酸素と呼ぶ)，高分子的な様々な構造を形成している．[SiO₄]^4-が独立して存在するネソケイ酸塩，-Si-O-の連結の仕方が異なる，ソロ，サイクロなどm/n比に応じて多様なケイ酸塩構造があり，m=0の酸化ケイ素のテクト構造までがある(石英，クリストバライトなど)．

したがって，ガラスや溶融体の構造は，Si-O種の分布により特徴づける．すなわち，Siに配位している4個の酸素原子の内の架橋酸素の数を用いて，Si種をQ₀からQ₄までに分類し，その相におけるQ種の頻度分布でガラス・溶融体の構造を特徴づけることができる．このQ種分布はガラスについては²⁹Si-NMRにより解析可能である．しかし，ガラス転移点以上の高温ではQ種間の交換速度が速くて，Q種の分別定量は不可能であることが知られており [1]．高温ラマン分光法により溶融体でもQ種の分別定量が可能である [2]．

同様な様相を有する系は他にホウ酸塩があり，さらにはホウケイ酸塩とアルミノケイ酸塩などがある．

2.1 MD計算

計算にはMXDORTO/MXDTRICL [3]を用いた．静電相互作用の計算にはEWALD法を用い，運動方程式の差分方程式にはVerletの方法を用いた．積分時間間隔は2fsとした．

原子間相互作用モデルは次のものを用いた:   

$\begin{array}{c}{u}_{ij}\left(r_{ij}\right)=\frac{z_i{z}_j{e}^{\text{2}}}{\text{4}\pi {\epsilon }_0{r}_{ij}}-\frac{c_i{c}_j}{r_{ij}^{\text{6}}}+f_{\text{0}}\left(b_i+b_j\right)\text{exp}\left(\frac{a_i+a_j-r_{ij}}{b_i+b_j}\right)\\ +D_{\text{1}ij}\text{exp}\left(-{\beta }_{\text{1}ij}{r}_{ij}\right)+D_{\text{2}ij}\text{exp}\left(-{\beta }_{\text{2}ij}{r}_{ij}\right)\end{array}$

すなわち，クーロン項，van der Waals項，近接反発項，および共有結合を表す2つの指数関数項である．a，b，c，D，およびβのパラメータは固定であらかじめ結晶構造や物性などを用いて決まられたものである．イオン電荷(z)はNaは+1，Siは+2.222で定数であるがOの電荷は組成により変化する．Na₂O．2SiO₂ではz (O)=-1.2889である．

ガラス化する酸化物融体は，融体からの冷却過程において，

[高温] 液体 → (過冷却液体) → ガラス転移 → ガラス [低温・室温]という過程を経て，ガラス状態となる．MD計算はNa₂O·2SiO₂組成 (Na₂Si₂O₅)を対象とし，乱数発生による構造を用い，高温過ぎずに迅速に平衡化できる温度と考えられる温度から始めた．ここでは2273.15 K (2000°C)において平衡化を行い，100Kごとに段階的に温度降下し，それらの各温度でエネルギーと構造の緩和を調べた． ある温度から次の温度降下の開始は，後でしめすように，その温度に到達した時点から始まる初期の顕著なエネルギー緩和の終了時点とした．おおむね各温度での(NPT) MD計算の開始から10∼20 ns後である．したがって冷却速度は10⁸ K/s程度かそれより大きいことになる．

計算した系は4系で，単位セルあたりの原子数と1273.15 Kにおける基本セル(ほぼ立方体)長を示す(Table 1)．

Table 1. The number of atoms and cell edge length of each system

		Basic cell
System	No. atoms	Edge length/Å(1273.15 K)
N5400	5400	43.13
N9000	9000	51.19
N18000	18000	64.48
N36000	36000	81.43

N18000系の冷却過程におけるエンタルピーと密度，および自己拡散係数の変化をFigure 1に示す．高温(2273.15 K)から温度降下によりエンタルピーはほぼ直線的に減少し，密度も直線的に増大しているが，1600 K付近でそれぞれ別の直線に乗り移っている．高温側と低温側の2直線の交点の温度はこのMD計算における極めて急冷でのガラス転移に関する仮想温度と考えられる．あとで示すように，エンタルピーと密度はこの屈曲付近の各温度における長時間の緩和によりそれぞれ減少と増大し(密度の増大が上矢印であらわされている)，新たな交点は1200K付近まで下降しているようである．このモデル系のガラス転移温度はこの温度よりさらに低いものと考えられる．

Figure 1.

 Enthalpy (●) and density (◆) as function of temperature of the N18000 system of Na2O·2SiO2 melt/glass. The arrow indicate the direction of relaxation with time in (NPT)-MD calculations.

1273·15Kにおける(NPT)-MD計算によるエンタルピーと密度，およびQ種とO種の存在割合の時間変化を，N5400，N9000，N18000およびN36000系のそれぞれについて，Figure 2，Figure 3，Figure 4およびFigure 5に示す．この温度ではSi-O結合の組み換えによる構造緩和がMD計算の時間の範囲内で観測可能であると考えられる．N5400系とN9000系ではプロットの1点が100万ステップ×2fs=2nsであり，N18000系とN36000系では50万ステップ×2fs=1nsである．したがって，N5400系は全体として250ns，から最も大きなN36000系で47nsの計算となっている．エンタルピーの変化から，N5400系では150ns程度で，N18000系では55ns程度で緩和が終わっているように見えるが，N9000系とN36000系ではこの一連の計算では緩和が終わっていないように思える．また，N5400系とN9000系では，エネルギー緩和は明らかに起こっているが，より大きな2系に比べて，その揺らぎが顕著に大きい．

Figure 2.

 Results of (NPT)-MD calculations at 1273. 15 K and 0.1 MPa of N5400 system. Enthalpy and density are plotted in upper graph, and Q-species and O^2- ratio (×) are plotted in lower graph.

Figure 3.

 Results of (NPT)-MD calculations at 1273.15 K and 0.1 MPa of N9000 system. Enthalpy and density are plotted in upper graph, and Q-species and O^2- ratio (×) are plotted in lower graph.

Figure 4.

 Results of (NPT)-MD calculations at 1273.15 K and 0.1 MPa of N18000 system. Enthalpy and density are plotted in upper graph, and Q-species and O^2- ratio (×) are plotted in lower graph.

Figure 5.

 Results of (NPT)-MD calculations at 1273.15 K and 0.1 MPa of N36000 system. Enthalpy and density are plotted in upper graph, and Q-species and O^2- ratio (×) are plotted in lower graph.

このようなエネルギーと密度の変化に対して，構造の詳細がどのように対応して変化しているかを調べるために，各系のQ種分布とO^2-の割合の変化を調べた(Figure 2からFigure 5の下図)．

温度降下に対する構造緩和は，主なものとして，Q種分布とO種分布の変化に現れると考えられる．他にNa原子の配置もあるが，前2者は，Si-O結合の開裂・再結合とネットワーク構成原子(SiとO)の拡散が関わり，後者は大きな自己拡散係数を持っているネットワーク修飾原子(Na)の移動による緩和である．

2元系ケイ酸塩(M₂O-SiO₂)溶融体・ガラスの場合，網目形成元素は Si原子のみで，網目修飾元素は M原子(1価)のみである．Na₂O·2SiO₂=Na₂Si₂O₅結晶では全てのOはSiに結合しており，Si種はQ₃のみからなっている(SiO₂ではQ₄のみ，Na₂SiO₃ではQ₂のみ)．溶融体ではQ_0，Q_1，Q₂， Q_3，Q₄ が混在しており，その存在比は温度とともに変化することが知られている [2] (Figure 6)．すなわち，温度低下によりエントロピーの寄与が相対的に小さくなり，次の均化反応が進むとされる．    Q_n-1 + Q_n+1 -> 2Q_n   

Figure 6.

 Q_n distribution of Na₂O·2SiO₂ melt/glass as a function against temperature. MD calculated molar fractions of Q_n species of the N18000 system are displayed.

Na₂O·2SiO₂の場合にはn=3，すなわちQ₃の方向へSi種が寄っていくことが期待される．ただしMD計算で作成したこの組成の溶融体では，O種はO⁰(Si-O-SiのO:架橋酸素)，O⁻(Si-O．．NaのO:非架橋酸素)に加えて，O^2-(Siに結合していない酸素)が少量存在している．したがって，次の反応も温度低下と緩和により起こっていることが考えられる:    2Q₄ + O^2- → 2Q₃   

またO種のみの間に関しても，上式にも含まれているが次の様な均化緩和が考えられる:    O⁰ + O^2- → 2O⁻   

Q種分布は1273.15Kの定温過程においてあまり変化していない．F Maeharaらの実験測定 [2]によるQn種のモル分率で、1000 °CでQ₃が60%，Q₂とQ₄がそれぞれ20%程度になっているが，MD計算ではいずれの系でもそれに比べてQ₃は約20%少なく，Q₄は10%多くなっている．N18000系においてはO^2-種の存在比が時間とともに減少していることが明らかである(Figure 4下)が他の系では減少しているとは認められない．

N18000系の1773·15Kから1073·15KのQ種分布はTable 2のようになってる．この冷却過程においてQ種分布はほとんど変化していないと言える．

Table 2.  Ratio of Q species in the N18000 system.

T/K	Q0	Q1	Q2	Q3	Q4
1773.15	0.24	3.93	21.52	44.40	29.90
1673.15	0.38	4.16	21.03	43.73	30.69
1573.15	0.24	3.56	20.96	45.96	29.28
1473.15	0.33	3.25	21.88	43.88	30.66
1373.15	0.17	3.89	20.19	45.92	29.83
1273.15	0.16	3.81	20.14	45.79	30.09
1173.15	0.16	3.27	20.60	46.00	29.98
1073.15	0.19	3.45	20.08	45.67	30.62

このようなMD計算によるQ種分布はFigure 6における1300°C = 1573.15Kを超える温度のQ種分布に相当するものである [2]．すなわちMD計算では温度降下に対する構造緩和が十分には行われていないと考えられる。

以上見てきたように，ケイ酸塩溶融体のMD計算は，1万原子以上のなるべく多くの基本セル内原子数を必要とし，かつ溶融体の作成においてもかなり長い計算ステップ(100 ns以上)を必要としており，それでも構造緩和が十分に行われていないように思われる．熔融体からガラス状態への緩和、すなわちガラス転移に関してはさらに極めて長い計算ステップが必要であろう．このように大規模計算が不可欠であるので、原子間相互作用モデルが妥当であるかの検討もさらに必要に思える．ここで用いている原子間相互作用モデルでは1100 K以下ではSi-O結合の開裂・再結合の頻度が非常に小さく，構造緩和がMD計算の時間では極めて起こりにくい．

このような困難があっても，酸化物溶融体/ガラスの本質を解明するためには，NMR測定ではわからない構造の詳細やXRD，X線や中性子線小角散乱(SAX)，赤外分光・ラマン散乱などの測定結果の解析にMD計算が有効であり必要であると考えられる．

3 無機凝集体のMD計算

酸化物結晶についても古典MD計算が適切な研究手法と考えられる興味ある様々な系がある． その1つは 固溶体についてである．無機材料物質や鉱物結晶は多くの場合複雑な固溶体となっている．固溶体の挙動が簡単な固溶モデル，すなわち正則溶液モデルなどで説明できる場合もあるが，一般には物理・化学的に意味のあるモデルで説明することは難しいことが多い．MD計算を用いて平衡状態の原子・イオン配置を求めることはほぼ不可能であろう．しばしば結晶内拡散は極めて遅く，また各種欠陥構造の存在も関係するからである． しかし種々の原子配置を生成し，MD計算で検証すること事は比較的容易である．また，実際には出現しない非平衡構造を調べることも可能であり，構造-物性の関係を理解するために重要なこと事である．

いくつかの無機結晶固体は誘電体や固体電解質となる．BaTiO₃は高温では等方体だが，室温では強誘電体であり、対称心のない正方晶となって自発分極する．このような挙動を再現できる古典MD計算は確認していないが容易ではないであろう．特に多結晶体で(結晶粒界を含めて)，このような挙動を調べることができれば，材料科学などで非常に意義があると考える．固体電解質も同様で，ZrO₂は高温(立方晶)でOの結晶内拡散はMD計算で実現できるようであるが，室温の単斜晶への転移を，予定調和的ではない原子間相互作用モデルで再現した例は存在しないように思う．これらの相転移はいずれも変位型相転移で，古典MDでも再現できそうに考えられる．固体電解質・イオン伝導結晶固体でも同様である．

構造相転移(特に圧力誘起相転移)は変位型と再構成型に分けることができる．変位型は強誘電体や固体電解質への相転移にしばしば現れる．再構成型は，骨格結合が開裂・再結合するものと言え，特に圧力誘起相転移では配位数が変化するものが多い．このような相転移には古典MD計算でも実現できるものもある．

これらの現象には多かれ少なかれ結晶の欠陥構造も関係している，点欠陥(化合物では欠陥対など)， 線欠陥(転位)，面欠陥(積層欠陥など)，さらに亜粒界，や粒界・表面などを含んだ分子シミュレーションは計算規模が大きくなることもあり，まず古典的なモデルでやってみることが薦められる．

4 おわりに

無機凝集体の物理化学の発展に分子シミュレーションが寄与できると考える．メルト/ガラスや結晶固体などの物性化学は，本質的に構造と性質を体系化するものであり，現代無機化学の主題である．地球・惑星構成物質の化学も同様である．無機化学は有機物ではない物質の化学を扱うという事だと思うが，フランス語でChimie inorganiqueという語も使われるが，Chimie mine'raleとも呼ばれる．

References

• [1]S． Shimokawa，， H． Maekawa． E． Yamada，Y． Nakamura， and T． Yokokawa， Chem． Lett．， 617 (1990/4)
• [2]   T． Maehara， J. Non-Cryst. Solids， 351， 3685 (2005)． doi:10.1016/j.jnoncrysol.2005.10.003
• [3]K． Kawamura， MXDORTO， Japan Chemistry Program Exchange， #029 (1996)