Interatomic Potential for xNa₂O-(1−x)SiO₂ Glasses

Abstract

The interatomic potential for a wide composition range of sodium silicate glasses was proposed by first-principles calculation. Point charge was set for each glass composition on the basis of population analysis of the alkali silicate crystals by using the density functional theory. The potential parameters were obtained from the energy surface of the SiO²⁺ model by using the molecular orbital method. The molecular dynamics simulation using the new interatomic potential showed improved structures of sodium silicate glasses.

1 目的

ケイ酸塩ガラスは元素の混合比を変えることで多様な性質を持つ．これは組成によってガラスの原子構造が変化するからであるとされている．そのため各組成の構造を詳細に解析することはガラスの物性の理解に必要である．分子動力学(Molecular Dynamics, MD) 法はガラスの構造解析によく用いられる．同じ元素から成る系の古典MD 計算には全く同じ原子間相互作用を適用することが多い [1,2]．しかしケイ酸塩ガラスの場合，組成に依存して構造が変化することが知られており，組成に合わせた原子間相互作用でなければ組成毎の構造を正しく表現できない恐れがある．非経験的な計算から組成毎に設定した原子間相互作用を使用することは，組成固有の構造や物性を表現するために重要であると考えている．これまでに，第一原理計算を用いて非経験的なアプローチから，SiO₂ 多形結晶に適用する原子間相互作用を設定する試みが行われていた [3]．しかし，第一原理計算を用いてガラスに適用する原子間相互作用を組成毎に最適化する研究は未だない．そこで本研究では第一原理計算の結果を利用し，組成毎に対応させる原子間相互作用について検討した．

2 方法

本研究では(1) 式に示す佐久間ら [4] が使用した原子間相互作用を用いた．   

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

r_ij は原子i, j間の距離, z_i はi番目のイオンの点電荷，e は素電荷，a_i(b_i) はBorn-Mayer 型の斥力における有効半径(ソフトネスパラメータ)，f₀ は定数(41.865 kJ nm⁻¹ mol⁻¹) ，c_i はvan der Waals 力の係数，D_1ij, D_2ij, β_1ij, β_2ij は共有結合力の係数である．z_i はガラスと同元素からなるアルカリケイ酸塩結晶の電荷解析から求めた．我々はLi, Na 等のアルカリ金属イオンだけではなく，Ca 等のアルカリ土類金属イオンを含むガラスについて解析を行いたいと考えている．そのため，複雑な構造をもつ結晶をはじめ，数多くの結晶の電荷解析を行う必要がある．このことから，結晶の電荷を少ない計算コストで精度を維持した解析手法を用いて検討する必要があると考えている．ガウス型関数等の局在関数基底系は，基底関数のサイズを大きくするとハミルトニアンの対角化が困難になる問題がある．平面波基底はカットオフ波数を増やすことで系統的に基底関数のサイズを拡大でき，また直交性が保障されているため，この問題は生じないとされている [5]．そのため，z_i は平面波基底密度汎関数理論(Density functional theory, DFT) 計算による電荷解析から求めた．この際，ソフトウェアはCASTEP [6] を用いた．対象とした結晶は以前の研究 [7] で用いたリチウムケイ酸塩結晶のほかに，ナトリウムケイ酸塩結晶Na₂Si₂O₅ [8], Na₂SiO₃ [9], Na₄SiO₄ [10], Na₂O [11] を加えた．k 点はそれぞれ，2 × 1 × 3, 3 × 3 × 3, 3 × 3 × 2, 5 × 5 × 5 とした．他の計算条件は以前の研究 [6] に合わせた．電荷解析にはMulliken 法 [12]とHirshfeld 法 [13] を用いた．SiO²⁺ モデルを用いSi-O 間のポテンシャル面をMP2[14]/6–311+g (d) [15] レベルの分子軌道(Molecular Orbital, MO) 計算で求めた．ソフトウェアはGaussian09W [16] を用いた．このSi-O ポテンシャル面に組成毎に定めたz_i を適用した式(1) を，非線形最小二乗法でフィッティングさせ，Si とO のa_i, b_i, c_i, D_1ij, D_2ij, β_1ij, β_2ij を得た．Na のa_i, b_i, c_iは経験的に設定した．xNa₂O-(1−x) SiO₂ (x = 0.3, 0.4, 0.5) ガラスのMD 計算をMXDORTO [17] コードを用いて行った．NPT アンサンブルで計算を行い，刻み時間は1.0 fs とした．はじめに，一辺が約50 Åのセルに粒子約8000 個をランダムに配置させた．次に，冷却速度を0.01 K/step (1.0 × 10¹³ K/s)とし，0.1 MPa で3000 K∼300 K間を合計4270 psで緩和した．300 Kの構造を50 ps 間について解析した．本研究で新たに定めた原子間相互作用(新原子間相互作用)と我々が従来使用してきた，z_i のみを組成毎に変化させた経験的原子間相互作用 [18] (旧原子間相互作用)を使用し，得られた構造と比較した．

3 結果・考察

Figure 1 にMulliken法とHirshfeld 法より得たアルカリケイ酸塩結晶の電荷をサイト毎に示す．Figure 1 (a) Mulliken 法のリチウムおよびナトリウムケイ酸塩結晶のSi の電荷はx が増加すると減少する傾向を示した．(b) Hirshfeld 法のリチウムケイ酸塩結晶のSi の電荷はx の変化に対し一定となった．アルカリケイ酸塩ガラス中のSi の電荷はx の増加に伴い減少することが蛍光X 線の結果 [19,20] から分かっており，Mulliken 法の結果が実験値の傾向と一致した．また，Mulliken 法ではx が増加するとアルカリ金属イオンの電荷は減少する傾向を示すが，Hirshfeld 法では一定となった．以前の我々の研究において，リチウムケイ酸塩結晶のMD 計算では，x 増加に伴いアルカリ金属イオンの電荷も減少するMulliken 法の傾向を用いたほうが，x の広い範囲で構造を改善できることを示した [7]．そこで，ガラスのMD 計算でもMulliken 法の傾向を反映させることにし，各イオンの電荷を式(2)~(4) で与えた．   

$z_{Na}=-0.25x+1.0\left(0<x\le 1\right)$(2)

  

$z_{Si}=-0.50x+2.4\left(0\le x<1\right)$(3)

  

$z_O=-\frac{2x\cdot z_{Na}+(1-x)\cdot z_{Si}}{2-x}$(4)

Figure 1.

 (a) Mulliken and (b) Hirshfeld charges of xA₂O-(1−x) SiO₂ crystals (A = Li (blank) or Na (oblique line))(○:A, □:Si, △:O).

式(4) によってO の電荷は電気的中性を保つ．Figure 2 にx = 0.3 を例にSi-O ポテンシャル面に式(1) をフィッティングさせた結果を示す．全エネルギーはSi-O 間の結合距離が0.146 nm で最小となった．同様の解析により組成x 毎に求めた原子間相互作用のパラメータをTable 1 に示した．新原子間相互作用と旧原子間相互作用を用いた場合の原子間距離を調べた．旧原子間相互作用を用いた場合，Na-O, Si-Si 間距離が結晶構造データとは大きく異なった．特にx = 0.5 のとき，SiO₄ 四面体同士が稜を共有する構造を取った．この構造はZachariasen [21] のガラス形成則と矛盾する．新原子間相互作用を用いた場合，ガラス中の原子間距離は結晶に見られる原子間距離の分布Si-O (0.148–0.167 nm), Na-O (0.220–0.253 nm), O-O (0.245–0.277 nm), Si-Si (0.300–0.322 nm) に収まった．また，ガラスの密度は実測値ではx = 0.3, 0.4, 0.5 でそれぞれ，2.466, 2.532, 2.560 g cm⁻³ である [22]．旧原子間相互作用による密度はそれぞれ，2.056, 2.077, 2.069 g cm⁻³ となり傾向を再現できなかった．新原子間相互作用による密度はそれぞれ，2.044, 2.186, 2.296 g cm⁻³ となり，実測値の組成変化の傾向を再現した．Figure 3 にSi に配位している架橋酸素の数n で分類したSiO₄ 四面体ユニットの存在比Q_nを示す．Figure 3 より新原子間相互作用は旧原子間相互作用の結果よりもNMR [23] のQ_nの増減傾向の再現性が向上しており，x = 0.5 のガラスについても構造が改善できていることが分かった．しかし，Q_nの値を再現するには至らなかった．この原因として，MD シミュレーションは現実の系に比べ，(1)冷却速度が速い，(2)粒子数が少ない，(3)時間スケールが短いことが挙げられる [1]．我々は，将来的には実験による取り扱いが困難な組成のガラスに対しても，大規模粒子系を扱える古典MD 法で構造解析を行いたいと考えている．そのため，計算された値が実験値に限りなく近く，また，その組成変化を幅広いガラス組成で再現できる原子間相互作用を求めている．新原子間相互作用は旧原子間相互作用よりも広範囲のアルカリ組成で構造を改善できることを示唆した．今後は，比熱などの熱力学的値について再現性が向上できているか検討する予定である．

Figure 2.

 Energy surface of SiO²⁺ model calculated by MO and fitting result (x = 0.3).

Table 1.  Interatomic potential parameters used in this study for xNa₂O-(1−x) SiO₂ glasses

Figure 3.

 Q_n distribution of simulated glasses with experimental data (○:Q₄, △:Q₃, □:Q₂, ◇:Q₁, +:Q₀).

4 結論

第一原理計算を用いてガラス組成毎に最適化した原子間相互作用を検討した．各原子の電荷は平面波基底のDFT 計算を用いて結晶の電荷解析をもとに組成毎に設定した．電荷以外のパラメータはMO 計算を用いてSiO²⁺ モデルから求めたポテンシャル面に原子間相互作用をフィッティングさせた．本研究の原子間相互作用は，分子動力学法によるガラスの構造解析を幅広い組成範囲に対してサポートできることが分かった．

参考文献

• 1   J. Du, et al., J. Non-Cryst. Solids, 349, 66 (2004).
• 2   R. N. Mead, G. Mountjoy, J. Phys. Chem. B, 110, 14273 (2006).
• 3   S. Tsuneyuki, M. Tsukada, H. Aoki, Y. Matsui, Phys. Rev. Lett., 61, 869 (1988).
• 4   H. Sakuma, et al., Geochim. Cosmochim. Acta, 73, 4100 (2009).
• 5   H. Kasai, et al., Introduction to Computational Materials Design, Osakadaigakushuppan (2005).
• 6   J. Clark, et al., Z. Kristallogr., 220, 567 (2005).
• 7   Y. Yamamoto, et al., J. Comput. Chem. Jpn., 13, 173 (2014).
• 8   A. K. Pant, et al., Acta Crystallogr., 24, 13 (1968).
• 9   A. Grund, et al., Acta Crystallogr., 5, 837 (1952).
• 10   W. H. Baur, et al., Monatsh. Chem., 117, 793 (1986).
• 11   E. Zintl, et al., Zeit. Elektrochem., 40, 588 (1934).
• 12   R. Mulliken, J. Chem. Phys., 23, 1833 (1955).
• 13   F. L. Hirshfeld, Theor. Chim. Acta, 44, 129 (1977).
• 14   C. Møller, et al., Phys. Rev., 46, 618 (1934).
• 15   L. A. Curtiss, et al., J. Chem. Phys., 103, 6104 (1995).
• 16   M. J. Frisch, et al., Gaussian 09, Revision A. 1, Gaussian Inc., Wallingford CT, (2009).
• 17   K. Kawamura, MXDORTO, JCPE, #29.
• 18   N. Sawaguchi, et al., J. Comput. Chem. Jpn., submitted.
• 19   S. Sakka, et al., J. Non-Cryst. Solids, 22, 57 (1976).
• 20   T. Maekawa, et al., Bull. Chem. Soc. Jpn., 51, 780 (1978).
• 21   W. H. Zachariasen, J. Am. Chem. Soc., 54, 3841 (1932).
• 22   P. B. Narottam, et al., Handbook of Glass Properties, Academic, Orlando, (1986).
• 23   H. Maekawa, et al., J. Non-Cryst. Solids, 127, 53 (1991).