酸化物ガラスを対象とした分子動力学計算に用いる
原子間相互作用の設定方法の最適化に関する研究
メタデータ

言語: Japanese
出版者:
公開日: 2019-06-25
キーワード (Ja):
キーワード (En):作成者: 山本, 優也
メールアドレス:所属:

URL

https://doi.org/10.15118/00009905

酸化物ガラスを対象とした分子動力学計算に用いる
原子間相互作用の設定方法の最適化に関する研究

室蘭工業大学大学院
工学専攻

工学研究科

博士後期課程

先端生産システム工学コース
学籍番号 16096004

山本 優也(Yamamoto Yuya)

2018 年 12 月 21 日
室蘭工業大学・学位論文

論文目次
第 1 章 緒言 ------------------------------------------------------------------------------------------------------------- 1§1 – 1. 目的 ------------------------------------------------------------------------------------------------------------ 1§1 – 2. 参考文献 ----------------------------------------------------------------------------------------------------- 14
第 2 章 xNa2O-(1 – x)SiO2 ガラスに適用する原子間相互作用の設定 ---------------------------------- 17§2 – 1. 緒言 ----------------------------------------------------------------------------------------------------------- 17§2 – 2. 原子間相互作用関数 -------------------------------------------------------------------------------------- 20§2 – 3. Coulomb 力項のイオンの電荷の設定 ----------------------------------------------------------------- 22§2 – 4. ポテンシャルカーブの形状の設定 -------------------------------------------------------------------- 26§2 – 5. 古典 MD シミュレーションの計算方法 ------------------------------------------------------------- 29§2 – 6. 結果 ----------------------------------------------------------------------------------------------------------- 31§2 – 7. 考察 ----------------------------------------------------------------------------------------------------------- 38§2 – 8. 結言 ----------------------------------------------------------------------------------------------------------- 42§2 – 9. 参考文献 ----------------------------------------------------------------------------------------------------- 43
第 3 章 0.25{xK2O-(1 – x)Na2O}-0.75SiO2 ガラスシミュレーションの混合アルカリ効果の再現性 ---------------------------------------------------------------------------------------- 46§3 – 1. 緒言 ----------------------------------------------------------------------------------------------------------- 46§3 – 2. 原子間相互作用 -------------------------------------------------------------------------------------------- 48§3 – 3. MD シミュレーションの計算条件 --------------------------------------------------------------------- 51§3 – 4. 結果 ----------------------------------------------------------------------------------------------------------- 52§3 – 5. 考察 ----------------------------------------------------------------------------------------------------------- 61§3 – 6. 結言 ----------------------------------------------------------------------------------------------------------- 80§3 – 7. 参考文献 ----------------------------------------------------------------------------------------------------- 81
第 4 章 MD シミュレーションによる xSiO2 - (1 – x)NaAlO2 ガラスの構造評価 --------------------- 84§4 – 1. 緒言 ----------------------------------------------------------------------------------------------------------- 84§4 – 2. 原子間相互作用 -------------------------------------------------------------------------------------------- 87§4 – 3. MD シミュレーションの計算条件 --------------------------------------------------------------------- 92§4 – 4. 結果 ----------------------------------------------------------------------------------------------------------- 93§4 – 5. 考察 ----------------------------------------------------------------------------------------------------------- 99§4 – 6. 結言 --------------------------------------------------------------------------------------------------------- 105§4 – 7. 参考文献 --------------------------------------------------------------------------------------------------- 106
第 5 章 xNa2O - (1 – x)B2O3 ガラスのホウ酸異常現象の再現性の評価 --------------------------------- 109§5 – 1. 緒言 --------------------------------------------------------------------------------------------------------- 109§5 – 2. 原子間相互作用 ------------------------------------------------------------------------------------------ 112

§5 – 3. MD シミュレーションの計算条件 ------------------------------------------------------------------- 118§5 – 4. 結果 --------------------------------------------------------------------------------------------------------- 119§5 – 5. 考察 --------------------------------------------------------------------------------------------------------- 130§5 – 6. 結言 --------------------------------------------------------------------------------------------------------- 139§5 – 7. 参考文献 --------------------------------------------------------------------------------------------------- 140
第 6 章 総括 ---------------------------------------------------------------------------------------------------------- 144
謝辞 -------------------------------------------------------------------------------------------------------------------- 145

第 1 章 緒言§1 – 1. 目的
これまで、研究者の経験則と多くの物質を合成して特性評価を行う材料探索が行われてきたが、より効率よく材料開発を行うために、これらの実験データを様々な因子で整理し、その相関性を基に機能発現のメカニズムを推定、そして理論モデルを構築するという手法が考えられた。その典型例として、電気陰性度やイオン半径といった単純な説明変数を用いた経験的法則に基づいた
原子・電子スケールからの材料探索が挙げられる。この経験則により、現在に至るまで多くの新規材料や機能が発見されてきた。現在では、結晶性材料において、このような経験則に依ることなく材料探索を行う、第一原理計算とデータマイニングを併せた未知材料のデータベース化も行われている [1]。このように材料探索に関する研究は世界的に重要なテーマの一つとして注目されている。一般的に、結晶の材料科学で取り扱う元素の種類は 80 ~ 90 程度であるが、その組み合
わせは二元系で数千、三元系になると数万以上のバリエーションの材料が生まれる。また、組成比や結晶多形も考慮すると更に数が増える。一方、酸化物ガラスは非周期的な原子配列を有し 4元素以上からなる組成も多く、結晶性材料に比べると組成選択性が非常に高い。また、ガラスは熱履歴によってもネットワーク構造が変化することが知られているため [2, 3]、膨大なバリエー
ションの新規材料が生まれる可能性を秘めている。そのため、ガラスに対し材料探索の研究を進めることが可能となれば、新たな機能性材料の開発に貢献できると考えている。
無機ガラスの一つである酸化物ガラスは工業的に板ガラス、曇りガラス、光学ガラスなどとして製造されており、我々の生活の一部として無くてはならない材料となっている。本論文では、
このような酸化物ガラスについて触れるにあたり、これらのガラスに特有の性質に基づく定義が必要となるとし、（１）原子配列が X 線的に不規則な網目構造をもつ、及び（２）ガラス転移現象を示す固体を、扱うガラスとして定義した [4]。
（１）は構造上非晶質であることを示す。Figure 1 – 1 に結晶とガラスモデルの模式図を示す。
結晶は原子が規則的に配列する。一方、Zachariasen [5] により提唱された、全体にわたり原子が不規則に配列するモデルである不規則構造説は、Warren [6] 、Sun [7] を始め多くのガラス研究者によって支持されている。Zachariasen はガラスが形成されるためには、不規則構造においても内部エネルギーが結晶とほぼ同一であること、また、構造単位が繋がって三次元の連続網目を構成することが必要であると考えた。これに基づいて、Zachariasen は酸化物 AmOn のガラス形成条件として次の四項目を挙げた。
1. 陽イオン A の配位数は小さくなければならない
2. 1 個の酸化物イオンは 3 個以上の陽イオンと結合しない
3. 酸素多面体は隣接多面体と頂点のみを共有し、面や量を共有しない
4. 各酸素多面体は少なくとも、頂点を共有する。
1 – 3 の条件は不規則網目と規則網目の内部エネルギーの差が小さいための条件であり、4 は三次元の網目構造を作るための条件である。これらの四条件を満足する単一酸化物に B2O3, As2O3,
SiO2, GeO2, P2O5 などがあるが、これらの酸化物は実際に急冷しなくてもガラスを形成する。これ
らの単独でガラスを形成する酸化物を Sun は網目形成酸化物と名付けた。これらの単一酸化物にさらに成分数を増やしたガラスに対しては、同一の基準で述べられるとは限らなくなる。あくま
1 / 145でも、Zachariasen は単一酸化物に対して、このような仮説を立てている。Zachariasen は網目形成酸化物にそれ自体ではガラスを形成しない酸化物(網目修飾酸化物, 例えばアルカリ金属酸化物、
アルカリ土類金属酸化物など) を連続網目が壊れない程度に加えても、ガラスが形成されると考えた。添加される網目修飾酸化物のガラス形成に対する影響について、Zachariasen はアルカリ金属またはアルカリ土類金属イオンが大きいほど、陽イオン間の静電反発力が小さいために、ガラ
スと結晶のエネルギー差が小さく、ガラスが形成されやすいと考えている。Figure 1 – 1 (c) は微結晶構造説に基づく不規則構造モデルである。図中の破線で示した領域が微結晶部分であり、微結晶の大きさは 2 nm 以下で、その割合は図のように少ないものから微結晶どうしが接触する程度に多いものまであるとされている。このような微結晶を繋ぎ合わせて固体とするためには、非晶質のマトリックス部分が存在する必要がある。微結晶の向きが不規則であるため、微結晶が小さければ、ガラスの X 線非晶質性、等方性、透明性といった物性を不規則網目構造説と同様に説明できる。微結晶説は Randall らによって提唱され [8]、Porai – Koshits はこの説に修正を加えている [9]。また、SiO2 ガラスなどの X 線的に非晶質の材料に電子顕微鏡暗視野像中に明るいスポッ
トが観察され [10, 11]、これが微結晶存在の証拠であるともされている [11] が、このコントラス
トは不規則網目構造モデルで説明できるとの意見 [10] もあり、このガラスの構造モデルに関する問題の解決は容易ではない。本論文では、Zachariasen の不規則網目構造説を基に酸化物ガラスの議論を行う。
（２）は非晶質の中で特に、ガラスと言えるものを区別する条件である。Figure 1 – 2 にガラス形成液体の体積の温度変化を示す。ガラスを加熱すると膨張するが、ガラス転移温度と呼ばれる温度で、ガラスは過冷却液体に変わり、膨張係数が大きくなる。これがガラス転移現象である。
過冷却液体を冷却するときには、逆に液体がガラス転移温度で固体のガラスに変わる現象がガラス転移である。ガラス転移温度 Tg は融点 Tm の 2/3 に近く、また Tg ではガラスの種類に関わらず粘度が 1013 P に近い値を取ることが知られている [12]。高温の融液を冷却してガラスを作る場
合(A → B → C)、ガラス転移温度付近では粘度が高いため、融液の構造は各温度での平衡構造にまで変化しないまま冷却され、ガラスとなる。この場合、冷却速度が大きいほど高温の平衡構造を取る。ガラスを転移温度領域で再加熱する(C → D) と、ガラスのもつ平衡構造に対応する温度(仮想温度) と異なる温度では、その温度での平衡構造をとろうとする (D →ラスの緩和あるいは安定化と呼ばれている。

2 / 145

E)。この現象はガ

(b)

(a)

(c)

0.2 nmFigure 1 – 1. Structural model of crystal (a), glass (random network structure) (b) and glass (micro
crystalline) (c).Dashed line shows crystal region.

3 / 145

A
Melt

Molar volume, V

Q1 > Q2 > Q3
Glass
Cooling rate Q1

C

DBStabilization
E

Crystallization

Q2
Q3

CrystalGlass transition range

Fictive temperature Tf
Temperature, T
Figure 1 – 2. Schematic diagram of glass transition temperature.

4 / 145

酸化物ガラスは構造緩和に伴い、ガラス転移温度、イオンの自己拡散係数、粘性係数、電気伝
導度、熱膨張係数、体積弾性率、比熱などの物性が変化することが知られており、これは原子ス
ケールの構造やイオンのダイナミクスと関係している。したがって、ガラスの製造及び加工上、原子スケールの構造を知ることは重要であると考えている。一般的な酸化物ガラスでは、Si4+、B3+、
P5+ などの網目形成 (Network Former, NWF) イオンは NWF – O – NWF からなる架橋したネットワーク構造を有するのに対し、Na+ イオンや K+ イオンなどの網目修飾 (Network Modifier, NWM)
イオンを添加すると、酸素架橋構造を切断し、NWF – O1− … NWM1+ 構造となる (Figure 1 – 3)。
例えば、xNa2O − (1 – x)SiO2 系ガラスの熱膨張係数は Na2O/SiO2 比, x の増加に伴い増大することが分かっており、原子スケールで見たときの三次元の網目構造から二次元の平面的な構造へ変化することで説明づけられる [13, 14] (Figure 1 – 4)。

(a) Network forming (NWF) oxide glasses

NWF

O

NWF

(b) Add network modified (NWM) oxide to network forming (NWF) oxide glasses

NWF

O1−

NWM1+

Figure 1 – 3. Network former, (a) and network former/network modifier, (b).

5 / 145(a) 3-D network structure (x = 0)(b) 2-D or 1-D network structure (x = 0.5)

0.5 nmFigure 1 – 4. Conceptual diagram of three-dimensional network structure, (a) and two-dimensional
network structure, (b) (drawn by VESTA [15] ).

結晶性材料の探索を行う場合、結晶系、空間群、原子座標、格子定数、軸角などの格子空間の情報を基に、任意の元素の種類を変更しシミュレーションもしくは、実験的研究を逐次的に行うことで、材料の物性を推定する。一方、ガラスは Zachariasen の理論に基づけば、非周期的な原子配列を有するため、格子という概念がない。そのため、ガラス材料の探索を効率良く行うためには、ガラス中の基本構造を知る必要がある。ガラスの構造は実験的手法と計算科学的手法により調査されてきた。主に、実験的手法では X 線散乱測定及び中性子回折測定による原子間距離の情報 [16 – 20]、核磁気共鳴法 (Nuclear Magnetic Resonance, NMR)による配位数解析 [21 – 26]、赤外分光法及びラマン分光法によるガラス内部のリング構造の解析 [27 – 32] などが行われている。
一方、分子シミュレーションを用いた研究も行われてきた。例えば、乱数によるモンテカルロ
シミュレーション [33 – 36]、物質の電子状態を扱う第一原理計算 [37 – 40] などの手法がある。
前者は、基本的に乱数と確率を用いた構造最適化の手続きを行う。X 線散乱測定や中性子回折測
定から得られた実験的観測による構造情報と計算結果との二乗残差を最小とするまで乱数により試行させ、構造を予測する手法である。実験データに基づきシミュレーションを行うことから、
逆モンテカルロシミュレーションと呼ばれる。この方法は、仮定をほとんど行わずに解析できる
ことが大きな利点であるが、粒子を特徴づけるものが原子散乱因子のみであるため、多成分系への適用が困難な一面がある。後者は基本入力パラメータを原子番号と系の構造のみとし、実験結果を参照せず、系の電子状態を求める計算手法である。そのため、非経験的な結果をシュレーデ
ィンガー方程式に基づいて導くため、一般的に計算結果への信頼性が高いとされているが、原子核と核周辺の電子との相関を計算するため、対応可能な粒子数に限りが出てくる。
ガラスの構造解析を行うシミュレーション手法として、古典的分子動力学 (Molecular Dynamics,
MD) 法が用いることが多い。MD 法はガラスの物性や構造のデータを得ることができる。MD シ
6 / 145ミュレーションは系を構成している粒子の位置と速度などを用いて、与えられた粒子間の相互作用関数と Newton の運動方程式に従い、多数の粒子を同時に運動させ、全体として一つの物質系を構成させる手法である。一般的に、ガラスのシミュレーションは数千個以上の粒子数を必要とし、粘性係数や自己拡散係数などの輸送係数を求める場合は、長時間の解析データが必要となる。
計算コストを抑え、ガラスの物性と構造の情報をシミュレーションで得ることは、ガラス材料の開発を効率よく行うために重要だと考えており、古典的な MD シミュレーションは他の計算科学的手法よりも、ガラスのシミュレーションに適していると考えている。MD シミュレーションが大規模な計算を可能とするのは、粒子間に働く引力と斥力の相互作用 (以下、原子間相互作用)を経験的な関数, Uij(rij) で表現するためである。例えば、従来の MD シミュレーションより用いられる基本的な関数の形として式 (1 – 1) のようなモデルがある [41]。Figure 1 – 5 に式 (1 – 1 ) を示す。

𝑧𝑖 𝑧𝑗 𝑒 2

𝑎𝑖 +𝑎𝑗 −𝑟𝑖𝑗𝑈𝑖𝑗 (𝑟𝑖𝑗 ) = 4𝜋𝜀 𝑟 + 𝑓0 (𝑏𝑖 + 𝑏𝑗 )exp (0 𝑖𝑗

𝑏𝑖 +𝑏𝑗𝑐𝑖 𝑐𝑗) − 𝑟6(1 – 1)𝑖𝑗

rij はイオン i と j 間の距離を示し、右辺第一項目が Coulomb 力項、第二項目が近接反発力項、第三項目が van der Waals 力項からなり、それぞれ、zi, ai, bi, ci は変数を示す。e は素電荷、𝜀0 は真空の誘電率を、f0 は単位を調整するための定数 (4.185 kJ nm-1 mol-1) である。
本研究では、MD シミュレーションを酸化物ガラスの材料設計ツールとして、材料開発の分野に貢献したいと考えている。具体的には、MD シミュレーションであらかじめ計算を行い組成に対する物性値を算出し、材料研究に優先順位をつけることで、材料開発を効率化することができると考えている。また、MD シミュレーションは冷却速度や保持温度を変更させた計算も行うことができるため、ガラス製造のプロセスを効率化させるアイデアを与えるツールになると考えている。他には、この材料設計ツールでシミュレーションしたガラスの構造をデータベースとし構造モデルを提供することが可能となれば、他の高精度及び高コストの計算科学的手法の効率化に繋がると考えている。したがって、古典 MD シミュレーションで多様なガラス構造を網羅的に表現することができれば、ガラスの材料開発に貢献できると考えている。そこで、MD シミュレーションの再現性に大きく関わる要素として、原子間相互作用の設定方法が重要になると考えている。

7 / 145

Uij (rij )Uij (rij )Short repulsive force

van der Waals force

Coulomb force

rij
Figure 1 – 5. Interatomic potential by eq. (1 – 1). Total potential, Coulomb force term, short
repulsive force term and van der Waals force term show red, blue, green and yellow line,
respectively.

8 / 145一般的に MD 計算では、二体間の相互作用のみを適用し、シミュレーションを行うため、実験系研究で報告されている結果と定量的には必ずしも一致しない。しかし、組成変化に伴う物性値の傾向を相対的にシミュレーション可能であり、原子間相互作用の設定を最適化させることで、実験系研究に歩み寄った議論を行いやすくなると考えている。これまでの古典 MD シミュレーシ
ョンの原子間相互作用の設定方法は、対象となる物質の実験結果と一致するように、原子間相互作用の変数を試行錯誤的に設定することが一般的であった。そのため、経験論に頼りすぎることから、計算結果の信頼性に対して疑問視されることが少なからず存在し、また効率よく多くのガラス組成への計算を行うための設定方法を検討する必要があった。
そこで、Tsuneyuki ら [42, 43] は、第一原理計算から原子間相互作用を設定する方法を考えた。
Figure 1 – 6 に Tsuneyuki らの原子間相互作用の設定方法を示す。彼らはガラス中に存在するとされる結晶構造の一部をクラスター分子モデルとして抜き取り、第一原子計算を用いて原子間相互作用のカーブ形状, 𝑈𝑖𝑗 (𝑟𝑖𝑗 )を設定することで、計算結果の信頼性を向上させた。この設定方法による原子間相互作用を適用した MD シミュレーションは、ケイ酸塩結晶/ガラス、アルミノケイ酸
塩結晶及びアルミノリン酸塩結晶の構造を再現できることが報告されている [42 – 45]。しかし、この手法はケイ酸塩系やリン酸塩系ガラスといった酸素多面体がネットワーク構造中で 4 もしく
は 6 配位構造のみを取る組成には対応できるが、例えば、アルカリホウ酸塩ガラスのように、BO3と BO4 ユニットが共存する組成において、数十パターンもの分子モデルを検討する必要がある。
したがって、彼らの設定方法は組成への汎用性が低く、4 元素以上を含む系において使用された例が少ない。そのため、複雑なガラス組成のシミュレーションを行う場合は、試行錯誤で原子間相互作用の変数を調整することが現在においても主流となっている。
式 (1 – 1) の右辺第一項目の Coulomb 力項のイオンの電荷, zi はガラスの構造や物性値の再現性に大きく関わる。SiO2 ガラスにおいて、NWM 酸化物が増加すると Si − O 結合が弱化することが蛍光 X 線を用いた解析により報告されている [46]。このとき、Si イオンと O イオンがもつ電荷が組成比によって変動する。従来の原子間相互作用の設定方法ではこのような実験報告を基に、電荷に対し経験的な値を用いる場合が多い。Noritake ら [47] や Cormack ら [48] は Na2O-SiO2系ガラス/融体において、Na2O/SiO2 比を変数とする関数として Si イオンと O イオンに対し経験的な電荷の式を設定し、用いている。Sawaguchi ら [49] も Li2O-SiO2 系ガラス/融体において、
Li2O/B2O3 比に依存する経験的な式を原子間相互作用中のイオンの電荷に適用している。このよ
うに、組成に対して連続的にイオンの電荷を与える式を導入してきたが、対象とする系によっては実験的観測が困難な場合もあるため、材料設計の視点では最適な設定手法とは言えない。
これまでの原子間相互作用の設定方法の課題として、対象としたガラス組成の物性値を再現可
能なシミュレーションを行える原子間相互作用を試行錯誤的に設定する必要があり、このとき設定を行う使用者により、原子間相互作用の設定基準が何通りも生まれ、決定される原子間相互作
用も設定基準と同数存在していた。これは、既知の組成ならば、シミュレーションで検討したい
構造を調査する際に有効であるが、未知の組成の場合、設定する基準を設けることができない。
そこで、第一原理計算を基準とする設定方法が開発された [42, 43] が、組成への汎用性に課題があることが分かった。

9 / 145

Crystal

Cluster molecule

Potential, Uij (rij )/ Hartree

Interatomic potential

-590

-591

-592
0.10
0.20
0.30
Distance, rSiO / nm

MD simulation for glasses

Figure 1 – 6. A determination method of interatomic potential by first-principle calculation [42, 43].

10 / 145このような、原子間相互作用の設定方法の課題を解決するために、従来とは異なる方法で第一原理計算を利用することで、酸化物ガラスの材料設計ツールとして貢献できる MD シミュレーシ
ョンを行えるのではないかと考えることにした。そこで本研究では、第一原理計算を利用し、多様なガラス組成へ対応可能な原子間相互作用関数を設定する方法を考えることにした。
本研究で対象としたガラス組成を Figure 1 – 7 に示した。Na2O − SiO2 系ガラスは実用ガラスの基本組成であり、実験及び計算データが豊富に存在するため、原子間相互作用の設定について評価を行いやすいと考えた。Na2O − K2O − SiO2 系ガラスは混合アルカリ効果と呼ばれる、アルカリ金属イオンの自己拡散係数が減少する傾向を示すことで知られている。アルカリ金属イオンの低拡散性は実用材料にも利用されている現象の一つであり、MD シミュレーションで再現することができれば、拡散現象を予測する設計ツールとしての活用を見い出せる可能性があると考えている。Na2O − Al2O3 − SiO2 系ガラスは実用材料でも扱われており、マグマの主成分でもある。NWF
酸化物である SiO2 と NWF と NWM の二つの役割を有する中間酸化物の Al2O3 が含まれる系であり、ガラスを形成する酸化物のバリエーションを増やした際の原子間相互作用の汎用性について、評価できる系と考えている。また、Na2O − Al2O3 − SiO2 系ガラスの構造は複雑な酸素多面体
構造やリング構造からなることが報告されており、MD シミュレーションで再現することができ
れば、NWF イオンの配位構造からリング構造まで広範囲の構造ユニットのシミュレーションが期待できる。Na2O − B2O3 系ガラスはホウ酸異常と呼ばれる Na2O/B2O3 比に対し非線形な線熱膨張係数の傾向を示す。このホウ酸異常はこれまでの MD シミュレーションの研究結果でも再現できた報告例がない。また、従来の MD シミュレーションによる報告では、二体間相互作用に加え、角度を補正する三体間相互作用を加えることが多い [50]。本研究では、シンプルな二体間に働く相互作用のみで酸化物ガラスの構造を表現したいと考えている。そのため、二体間相互作用のみを適用した場合、Na2O − B2O3 系ガラスの構造をどこまで再現できるかという視点で記述している。熱膨張係数はガラス製造において、重要な物性値の指標の一つであり、あらかじめ計算科学的手法から熱膨張係数をシミュレーションすることができれば、MD 法の計算精度への信頼性について示すことができると考えている。

11 / 145

12 / 145

Figure 1 – 7. Schematic drawing of this work.

以上に述べた目的に従って、本論文は次のように 6 章から構成されており、以下に各章の要点を述べる。
第 1 章は、本研究を進める上での予備的考察について述べ、本研究の背景及び目的を示した。
第 2 章は Na2O − SiO2 系ガラスを対象に、原子間相互作用の設定方法について記述し、ガラスの構造と物性の再現性について示す。
第 3 章は Na2O − K2O − SiO2 系ガラスを対象に、第 2 章で示した原子間相互作用の設定方法を適用したシミュレーションを行った結果を記述する。
第 4 章は Na2O − Al2O3 − SiO2 系ガラスを対象に、第 2 章で示した原子間相互作用の設定方法を適用したシミュレーションを行った結果を記述する。
第 5 章は Na2O − B2O3 系ガラスを対象に、第 2 章で示した原子間相互作用の設定方法を適用したシミュレーションを行った結果を記述する。
第 6 章は本研究で得られた結果を総括する。

13 / 145§1 – 2. 参考文献
[1] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevier, K. A.
Persson, G. Ceder, Python materials genomics (pymatgen): A robust, open-source python library formaterials analysis, Computational Materials Science 68 (2013) 314 – 319.
[2] C. J. Brinker, R. J. Kirkpatrick, D. R. Tallant, B. C. Bunker, B. Montez, NMR confirmation of strained“Defects” in amorphous silica, J. Non-Cryst. Solids 99 (1988) 418 – 428.
[3] G. K. White, Thermal expansion of reference materials: copper, silica and silicon, J. Phys. D:Appl. Phys.
6 (1973) 2070 – 2078.
[4] S. Sakka, Gaarasu no kagaku to kiso to ouyo, Uchida Rokakuho Publishing Co., Ltd. (1997).
[5] W. H. Zachariasen, The atomic arrangement in glass, J. Am. Chem. Soc. 54 (1932) 3841 – 3851.
[6] B. E. Warren, X-Ray diffraction study of the structure of glass, Chem. Rev., 26 (1940) 237 – 255.
[7] Kuan – Han Sun, Fundamental condition of glass formation, J. Am. Ceram. Soc., 30 (1947) 277 – 281.
[8] J. T. Randall and H. P. Rooksby, The structure of glasses: The Evidence of X-ray diffraction, J. Soc.
Glass Technol., Proc., 14 (1930) 219 – 229.
[9]E. A. Porai-Koshits, The Possibilities and results of X-ray methods for investigation of glassy
substances, in The Structure of Glass, Lebedev, A.A., Ed., New York: Consultants Bureau, (1958) 25
– 35.

[10]

A. Howie, High resolution electron microscopy of amorphous thin films, J. Non-Cryst. Solids 31(1978) 41 – 55.

[11] Y. Bando and K. Ishizuka, Study of the structure of silica glass by high-resolution electron microscopy,
J. Non-Cryst. Solids 33 (1979) 375 – 382.
[12] S. Sakka and J. D. Mackenzie, Relation between apparent glass transition temperature and liquids
temperature for inorganic glasses, J. Non-Cryst. Solids 6 (1971) 145 – 162.
[13] R. Knoche, D.B. Dingwell, F.A. Seifert, S.L. Webb, Non-linear properties of supercooled liquids inthe system Na2O-SiO2, Chem. Geol. 116 (1994) 1 – 16.
[14] J. E. Shelby, Density of vitreous silica, J. Non-Cryst. Solids 349 (2004) 331 – 336.
[15] K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric andmorphology data, J. Appl. Crystallogr. 44 (2011) 1272 – 1276.
[16] R.L. Mozzi and B.E. Warren, The structure of vitreous silica, J. Appl. Cryst. 2 (1969) 164 – 172.
[17] G.S. Henderson, A Si K-edge EXAFS/XANES study of sodium silicate glasses, J. Non-Cryst. Solids
183 (1995) 43 – 50.
[18] N. Zotov, R. G. Delaplane, H. Keppler, Structural changes in sodium tetrasilicate glass around theliquid - glass transition: a neutron diffraction study, Phys. Chem. Minerals 26 (1998) 107 – 110.
[19]

A. C. Wright, B. A. Shakmatkin, N. M. Vedishcheva, The chemical structure of oxide glasses: A
concept consistent with neutron scattering studies ?, Glass Phys. Chem. 27 (2001) 97 – 113.

[20] M. Taylor and G. E. Brown Jr., Structure of mineral glasses-Ⅱ. The SiO2-NaAlSiO4 join, Geochim
Cosmochim. Acta 43 (1979) 1467 – 1473.

14 / 145

[21] J. F. Stebbins, Effects of temperature and composition on silicate glass structure and dynamics: Si-29
NMR results, J. Non-Cryst. Solids 106 (1988) 359 – 369.
[22] H. Maekawa, T. Maekawa, K. Kawamura, T. Yokokawa, The structural groups of alkali silicate glasses
determined from 29Si MAS-NMR, J. Non-Cryst. Solids 127 (1991) 53 – 64.
[23] H. Maekawa, T. Saito, T. Yokokawa, Water in silicate glass: 17O NMR of hydrous silica, albite andNa2Si4O9 glasses, J. Phys. Chem. B 102 (1998) 7253 – 7529.
[24] P. J. Bray and J. G. O’Keefe, Nuclear magnetic resonance investigations of the structure of alkali
borate glasses, Phys. Chem. Glasses 4 (1963) 36 – 46.
[25] S. K. Lee and J. F. Stebbins, Nature of cation mixing and ordering in Na-Ca silicate glasses and melts,
J. Phys. Chem B 107 (2003) 3141 – 3148.
[26] S. K. Lee and J. F. Stebbins, The structure of aluminosilicate glasses: High-resolution 17O and 27Al
MAS and 3QMAS NMR study, J. Phys. Chem B 104 (2000) 4091 – 4100.
[27] E. A. Robinson, Characteristic vibrational frequencies of oxygen compounds of silicon, phosphorus,
and chlorine: correlation of stretching frequencies and force constants with bond lengths and bondorders, Canadian J. Chem. 41 (1963) 3021 – 3033.
[28] P. Mcmillan, Structural studies of silicate glasses and melts – applications and limitations of Raman
spectroscopy, Am. Miner. 69 (1984) 622 – 644.
[29] A. E. Geissberger and F. L. Galeener, Raman studies of vitreous SiO2 versus fictive temperature, Phys.
Rev. B 28 (1983) 3266 – 3271.
[30]

P. Muralidharan, M. Venkateswarlu, N. Satyanarayana, Sol-gel synthesis, characterization and
impedance studies of lithium borosilicate glass, Mat. Res. Bulletin 39 (2004) 1753 – 1762.

[31] B. O. Mysen, D. Virgo, C. M. Scarfe, Relations between the anionic structure and viscosity of silicate
melts- a Raman spectroscopic study, Am. Miner. 65 (1980) 690 – 710.
[32] B. N. Meera and J. Ramakrishna, Raman spectral studies of borate glasses, J. Non-Cryst. Solids 159(1993) 1 – 21.
[33] L. Cormier, G. Calas, B. Beuneu, Structural changes between soda-lime silicate and melt, J. Non-Cryst.
Solids 357 (2011) 926 – 931.
[34] N. Zotoz and H. Keppler, The structure of sodium tetrasilicate glass from neutron diffraction, reverse
Monte Carlo simulations and Raman spectroscopy, Phys. Chem. Minerals 25 (1998) 259 – 267.
[35]

J. Swenson, A Matic, C. Karlsson, L. Börjesson, C. Meneghini, W. S. Howells, Random ion
distribution model: A structural approach to the mixed-alkali effect in glasses, Phys. Rev. B 63 (2001)
132202.

[36] S. Kerisit, J. V. Ryan, E. M. Pierce, Monte Carlo simulations of the corrosion of aluminoborosilicate
glasses, J. Non-Cryst Solids 378 (2013) 273 – 281.
[37] T. Tamura, S. Ishibashi, S. Tanaka, M. Kohyama, M. H. Lee, First-principles analysis of the opticalproperties of structural disorder in SiO2 glass, Phys Rev. B 77 (2008) 085207.
[38] G. Ferlat, T. Charpentier, A. P. Seitsonen. A. Takada, M. Lazzeri, L. Cormier, G. Calas, F. Mauri,Boroxol rings in liquid and vitreous B2O3 from first principles, Phys. Rev. Lett. 101 (2008) 065504.
15 / 145

[39] P. Umari and A. Pasquarello, Fraction of boroxol rings in vitreous boron oxide from a first-principles
analysis of Raman and NMR spectra, Phys. Rev. Lett. 95 (2005) 137401.
[40] T. Ohkubo, E. Tsuchida, M. Gobet, V. S.-Kanian, C. Bessada, Y. Iwadate, First-principles molecular
dynamics simulation and conductivity measurement of a molten xLi2O-(1 – x)B2O3 system, J. Phys.
Chem. B 117 (2013) 5668 – 5674.
[41] F. G. Fumi and M. P. Tosi, Ionic sizes and Born repulsive parameters in the NaCl-type alkali halides.
I., The Huggins-Mayer and Pauling forms, J. Phys. Chem. Solids, 25 (1964) 31 – 43.
[42] S. Tsuneyuki, M. Tsukada, H. Aoki, Y. Matsui, First-principles interatomic potential of silica applied
to molecular dynamics, Phys. Rev. Lett. 61 (1988) 869 – 872.
[43] S. Tsuneyuki, Molecular dynamics simulation of silica with a first-principles interatomic potential,
Mol. Engineer. 6 (1996) 157 – 182.
[44] G.J. Kramer, N.P. Farragher, B.W.H. van Besset, R.A. van Santen, Interatomic force fields for silicas,
aluminophosphates, and zeolites:derivation based on ab initio calculations, Phys. Rev. B 43 (1991)
5068 – 5080.
[45] J. Habsaki and I. Okada, Molecular dynamics simulation of alkali silicates based on the quantum
mechanical potential surfaces, Mol. Sim. 9 (1992) 319 – 326.
[46] S. Sakka and K. Matsusita, Studies on binary silicate glasses based on the SiKα and SiKβ emission Xrays, J. Non-Cryst. Solids 22 (1976) 57 – 66.
[47] F. Noritake, K. Kawamura, T. Yoshino, E. Takahashi, Molecular dynamics simulation and electricalconductivity measurement of Na2O・3SiO2 melt under high pressure; relationship between its
structure and properties, J. Non-Cryst. Solids 358 (2012) 3109 – 3118.
[48] A.N. Cormack, J. Du, T.R. Zeitler, Sodium ion migration mechanisms in silicate glasses probed bymolecular dynamics simulations, J. Non-Cryst. Solids 323 (2003) 147 – 154.
[49] N. Sawaguchi, K. Yamaguchi, M. Sasaki, K. Kawamura, Interatomic potential model for molecular
dynamics simulation of lithium borate melts/glasses, J. Comput. Chem. Jpn. 14 (2015) 139 – 146.
[50]

A. Takada, C. R. A. Catlow, G. D. Price, Computer modelling of B2O3: part 1. New interatomic
potentials, crystalline phases and predicted polymorphs, J. Phys.:Condens. Matter 7 (1995) 8659 –
8692.

16 / 145

第 2 章 xNa2O − (1 − x)SiO2 ガラスに適用する原子間相互作用の設定§2 – 1. 緒言
従来の原子間相互作用の設定方法は試行錯誤的に行うことが一般的であった。しかし、このような経験論に基づく方法は未知の組成の場合、指標とする構造や物性値が分からないため、設定に膨大な検討時間を要する。そのため、原子間相互作用の設定には、実験結果を指標としない方法を考えた。非経験的な結果を導く第一原理計算を利用し、原子間相互作用を設定することにした。対象とするガラスは、実用ガラスの基本組成である xNa2O − (1 − x)SiO2 ガラス(0 ≤ x ≤ 0.5)とした。
xNa2O − (1 − x)SiO2 ガラスにおける Si − O 間の結合は、x が増加するに伴い、弱化することが知られている [1]。このとき、Si と O の電荷は x 毎に変化する。従来の MD シミュレーションに適用する原子間相互作用中の Coulomb 力項の電荷は、シミュレーションのステップ数に関係なく固定した値を設定、もしくはシミュレーションのステップごとに各粒子に対し電荷を解析する
電荷平衡法 [2, 3] などが考えられてきた。本研究ではシンプルな設定方法を追究していることから、ステップ数に関係なく固定した電荷の値を与える方法を用いることにした。x により変化する各イオンの電荷は、ガラスと同組成の結晶中のイオンの電荷を解析することにより、見積もった。これにより、x 変化による Si − O 結合の弱化を表現した。ここで、xNa2O − (1 − x)SiO2 ガラ
スでは、組成, x の違いにより Si − O 間距離が大きく変化するという報告はない [4 − 7]。このことから、Si – O 間の結合の深さは x により変化するが、Si − O 間の平衡原子間距離は変化しないと考えた。そこで、ポテンシャルカーブの形状(曲率)と平衡原子間距離は酸化物ガラス中に多様な構造ユニットが存在する場合でも対応できるように、すべての組成に対しシンプルな二原子分子モデルに働くポテンシャルを第一原理計算で求め、得られたポテンシャル曲線を適用することとした (Figure 2 – 1)。

17 / 145

Potential, USiO (rSiO)

Figure 2 – 1. Schematic daiagram of Si – O potential in xNa2O − (1 − x)SiO2 glasses in this work.

18 / 145

§2 − 2 では MD シミュレーションで使用した原子間相互作用モデルを説明した。§2 − 3 では
アルカリケイ酸塩結晶の電荷解析を行い、その傾向を基に組成を変数とする電荷の式を連続関数として与え、イオンの電荷に適用した。§2 − 4 では SiO2+ モデルの第一原理計算で得た Si − O 間のポテンシャルエネルギー面に原子間相互作用を Na2O/SiO2 比毎にフィッティングさせ、得られた Si − O 間の原子間相互作用の結果を記述した。§2 – 5 では MD シミュレーションの計算条件を記した。§2 – 6 では MD シミュレーションで得られた計算結果を示す。§2 – 7 では§2 – 6 で得られた結果を考察した。§2 – 8 は本章の結言を記した。§2 – 9 は本章で引用した参考文献を示す。

19 / 145§2 – 2. 原子間相互作用関数
式 (1 – 1) は完全イオン型ポテンシャル関数と呼ばれ、NaCl や CsCl といったイオン性が高い結晶に対して主に使用されてきた [8]。しかし、酸化物系では系中のイオンが形式電荷とならない
場合が多く、ケイ酸塩結晶・ガラスのように Si – O 間に方向性をもつ共有結合性を考慮すること
で、実際の系に近い構造を表現できると考えられた。基本的にこのイオン性結合に関する考え方
は Pauling が提唱したイオン結合度 [9] に関連して考えられることが多い。Pauling が提案したイオン結合度の式では電気陰性度の差から Si – O 間のイオン結合性が約 60 %を示す。このことを背景に、従来の MD シミュレーションに適用する Coulomb 力項の Si と O の電荷はそれぞれ、+2.4 及び– 1 . 2 と試行錯誤的に設定することが多い [10]。 そして、残りの約 40 % の共有結合は無視できないとし、原子間相互作用の改善が行われた。しかしながら、あくまでも、Pauling のイオン結合度の式は大まかな目安にしか過ぎないため、原子間相互作用関数の考察は常に行われている [11]。本研究では、酸化物系結晶及びガラスにおいて、構造や物性の再現性が良いことが報告されている式 (2 − 1) を用いた [12]。Figure 2 – 2 に式 (2 – 1) を示す。

𝑈𝑖𝑗 (𝑟𝑖𝑗 ) =𝑧𝑖 𝑧𝑗 𝑒 2
𝑎𝑖 + 𝑎𝑗 − 𝑟𝑖𝑗
𝑐𝑖 𝑐𝑗+ 𝑓0 (𝑏𝑖 + 𝑏𝑗 )exp () − 6 + {𝐷1𝑖𝑗 exp(−𝛽1𝑖𝑗 𝑟𝑖𝑗 ) − 𝐷2𝑖𝑗 exp(−𝛽2𝑖𝑗 𝑟𝑖𝑗 )}4𝜋𝜀0 𝑟𝑖𝑗𝑏𝑖 + 𝑏𝑗𝑟𝑖𝑗(2 − 1)

右辺第四項目は共有結合力項であり D1ij, β1ij, D2ij, β2ij はそれぞれ変数である。

20 / 145

Uij (rij )Uij (rij )Short repulsive force

Covalent force
van der Waals force

Coulomb force

rij
Figure 2 – 2. Interatomic potential by eq. (2 – 1). Total potential, Coulomb force term, short repulsive
force term, van der Waals force term and covalent force term show red, blue, green, yellow and parple line,
respectively.

21 / 145§2 – 3. Coulomb 力項のイオンの電荷の設定
式 (2 – 1) の Coulomb 力項のイオンの電荷を設定するために、ガラスと同組成の結晶のイオンの電荷を解析することにより、見積もった。対象とした結晶は Na2O − SiO2 系に加え、x の変化による電荷の傾向を得るために、Li2O − SiO2 系、K2O − SiO2 系結晶についても計算を行った。イオンの電荷, zi は CASTEP code [13] を用いてアルカリケイ酸塩結晶, xA2O − (1 − x)SiO2 (A = Li, Na及び K) を密度汎関数理論計算 (Density Functional Theory calculation, DFT) を行うことによって設定した。平面波基底とノルム保存型擬ポテンシャルを用いた。ノルム保存型擬ポテンシャルは、Li : 2s1, Na : 2s22p63s1, K : 3s23p64s1, Si : 3s23p2, O : 2s22p4 を使用した [14]。交換相関汎関数はGeneralized Gradient Approximation Perdew − Burke − Ernzerhof (GGA PBE) [15] を適用した。Table 2− 1 に計算対象として選んだ結晶 [16 − 28] と k 点セットを示す。波動関数のカットオフエネル
ギーを 1500 eV とし、SCF 収束値は 5.0×10−7 eV / atom とし構造最適化計算を行った。構造最適化後のすべての結晶の格子定数と軸角は、文献データ [16 − 28] と比較して 2.88 % 以内の誤差に収まった。平面波基底を用いる Mulliken population analysis [29, 30] により、イオンの電荷を解析した。平面波基底を用いる Mulliken population analysis はカットオフ距離領域内に存在する電子を、対象とするイオンの占有軌道電子としてカウントする。このカットオフ距離はすべての結晶で 0.3
nm とした。
Figure 2 − 3 に xA2O − (1 − x)SiO2 (A = Li, Na 及び K) 結晶の電荷解析の結果を示す。x が増加するとともに、Si イオンの電荷の値が減少し、A イオンの電荷も減少傾向を示した。この Si イオ
ンの電荷が x が増加するに伴い減少する傾向を示した。これは xNa2O − (1 − x)SiO2 ガラスにおける蛍光 X 線の報告 [1] と一致している。x = 0 のとき、Si イオンの電荷の値が＋2.4 であり、これは Pauling のイオン結合度の式から導かれる Si – O 間のイオン性 60 % と一致した値を示しているが、相関関係はないと考えている。この Mulliken population analysis の結果を使用して導出した zi を用いた MD シミュレーションを行うと、ケイ酸塩ガラス・結晶の構造の再現性が良いことが分かっている [31, 32]。一方、別の電荷解析手法の一つである Hirshfeld population analysis [33]の結果を使用して導出した zi を用いたアルカリケイ酸塩結晶の MD シミュレーションではその結晶の構造の再現性が悪い [31, 32]。そのため、本研究ではガラスの構造の再現性が良くなると考えられる Mulliken population analysis を使用し、電荷を設定することにした。

22 / 145

23 / 145

Table 2 – 1. The diviation denotes the difference of structural optimization by DFT. Deviation shows in absolute value.

4
Li

Na

K

Mulliken charge

3
2

Si

1
Alkali

0
-1

O

-2
0

0.2

0.4
0.6
0.8
Composition, x

1

Figure 2 – 3. Mulliken charges of ions in several xA2O − (1 – x )SiO2 crystals (listed in Table 2 – 1 ).
Circles, triangles and squares show Li, Na and K ions, respectively.

24 / 145

Mulliken population analysis はオーバーラップポピュレーションを等分に分配するため、結晶の組成変化に対する傾向のみで判断することが多い [29]。そのため、本研究で用いる電荷は従来使用してきた経験的な値 [34 – 36] を考慮し、また、連続したガラス組成に対応するため、x を変数とした数式で示すことにした。以下の式 (2 − 2) − (2 − 4) に zNa, zSi, zO をそれぞれ示す。
𝑧𝑁𝑎 = −0.25𝑥 + 1.0 (0 < 𝑥 ≤ 1)

(2 – 2)𝑧𝑆𝑖 = −0.50𝑥 + 2.4 (0 ≤ 𝑥 < 1)

(2 – 3)𝑧𝑂 = neutrality(2 – 4)

実際の酸化物系を考えた場合、系中の酸化物イオンが最も電子を引き付けやすく、様々な電荷値のパターンを取ると考えられる。本論文において、ステップ毎に電荷を固定して適用する方法を取り入れているため、酸化物イオンの電荷, zO で電気的中性を保つことにした。

25 / 145§2 – 4. ポテンシャルカーブの形状の設定
式 (2 – 1) に用いる Si と O の ai, bi, ci, D1ij, β1ij, D2ij, β2ij は非経験的分子軌道(Molecular Orbital,MO) 計算より求めた。クラスター分子モデルとして、特定の配位構造を表すものではなく、シンプルな SiO2+ モデルを選んだ。これは、単純な二原子分子モデルとすることで、NWF と O イオンの多面体構造にバリエーションが表れても対応できるようにするためである。ソフトウェアはGaussian09 [37]、電子相関手法は Møller − Plesset (MP2) [38] とし、基底関数は 6 − 311 + g(d) [39]を使用した。SiO2+ モデルのポテンシャルエネルギー面を求め、これに§2 – 3 で設定した電荷を組み込んだ式 (2 − 1) をフィッティングさせ、Si − O 間の原子間相互作用を得た。このとき zi はx に依存して設定しているため、x 毎に ai, bi, ci, D1ij, β1ij, D2ij, β2ij を求めている。カーブフィッテ
ィングは MD シミュレーションの結果がアルカリケイ酸塩結晶の Si − O 間距離を再現できる原子間相互作用となるように行った。Figure 2 − 4 に例として、x = 0 のポテンシャルエネルギー面とそのフィッティング結果を示す。SiO2+ モデルの平衡 Si − O 間距離は 0.146 nm だった。Table
2 − 2 に設定した原子間相互作用を示す。ここで、Na – O 間に対して、カーブフィッティングを行っていない。これは、ガラス中の Na – O 間距離は Si – O 間距離に比べ、広範囲の分布を示す
ことが報告されているため [6, 7]、フレキシブルなポテンシャルとする必要があると考えている。
また、Na – O 間にカーブフィッティングを行うと試行変数が少ないことから、最小二乗法が収束しない。原子間相互作用の力の項を追加し、試行する変数を増加させることは、原子間相互作用を設定する上での経験的要素を排除できるかもしれないが、多変数の原子間相互作用はシミュレーションで得られる結果をより複雑化させるとも考えている。そのため、aNa, bNa, cNa のみは経験的に決めた値を使用した [35]。Figure 2 − 5 に本研究で考案した原子間相互作用の設定方法を示す。

Energy difference, ΔE ×1018 / J

0.030
0.025

SiO2+ model
O
Si

MO
MD

0.020
0.015
0.010
0.005
0.000
0.140

0.145
0.150
Distance, rSiO / nm

0.155

Figure 2 – 4. Energy surface determined from MO calculation (MP2 / 6 – 311 + g(d)) of an
SiO2+ model and potential energies calculated by eq. (2 – 1) with fitted parameters for x = 0.
26 / 145

Table 2 – 2. Detemined interatomic potential parameters for the simulation of xNa2O-(1 – x)SiO2 glasses.
(a) IP1x =0

z

a / nmb / nmOSi

neutrality－0.50x+2.4

0.184640.100070.014110.00799−1−1

c / (kJ / mol)1/2 nm30.05605
0.00000−1−1

D 1ij / kJ mol

β 1ij / nm

D 2ij / kJ mol

β 2ij / nm

Si-O

627600.00

52.05

83680.00

34.21

x = 0.1Ozneutralitya / nm0.18421b / nm
0.01401c / (kJ / mol)1/2 nm30.05605

0.09967

0.00798

Si－0.50x+2.4−1−10.00000−1

Si-O

D 1ij / kJ mol
627600.00

β 1ij / nm
52.06

D 2ij / kJ mol
83680.00

β 2ij / nm−1
34.20

x = 0.2
OSizneutrality－0.50x+2.4

Si-O

D 1ij / kJ mol−1627600.00

a / nm0.18397
0.09944
β 1ij / nm−1
52.06

b / nm
0.013960.00797
D 2ij / kJ mol−183680.00

c / (kJ / mol)1/2 nm30.05605
0.00000
β 2ij / nm−1
34.19

x = 0.3za / nmb / nmOSi

neutrality－0.50x+2.4

0.18362
0.09930

0.01391
0.00797c / (kJ / mol)1/2 nm30.05605
0.00000

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1

Si-O

627650.21

52.01

83680.00

34.17

x = 0.4za / nmb / nmO

neutrality

0.18326

0.01386c / (kJ / mol)1/2 nm30.05605Si－0.50x+2.4

0.099140.007960.00000

Si-O

D 1ij / kJ mol−1627696.23

β 1ij / nm−1
51.96

D 2ij / kJ mol−183680.00

β 2ij / nm−134.15

x = 0.5za / nmb / nmO

neutrality

0.18300

0.01383c / (kJ / mol) nm
0.05605Si－0.50x+2.4

0.09902

0.00796

−1

−1

1/2

3

0.00000−1−1

D 1ij / kJ mol

β 1ij / nm

D 2ij / kJ mol

β 2ij / nm

627725.52

51.9283680.00

34.14

z

a / nmb / nmc / (kJ / mol)1/2 nm3－0.25x+1.0

0.132200.01150

0.01227

Oz [35]
neutralitya / nm0.18610b / nm
0.01510c / (kJ / mol)1/2 nm30.05605Si

2.4+1.6x /(3x -4)

0.101200.008300.00000

Na

1.0000D 1ij / kJ mol−1222170.40

0.13220
β 1ij / nm−1
50.00

0.01150
D 2ij / kJ mol−113849.04

0.01227
β 2ij / nm−1
22.40

Si-O

Na
(b) IP2

Si-O

27 / 145

Determination method of interatomic potential§2 – 2

Interatomic potential
U ij (rij ) = ai + a j − rij  ci c j−+ D1ij exp(− 1ij rij ) − D2ij exp(−  2ij rij )
+ f 0 (bi + b j ) exp b + b  rij 6
4 0 rijijzi z j e 2

S tep 1 §2 – 3

4
Li

Na

S tep 2 §2 – 4

K

Mulliken charge

3
2

Si

1
Alkali

0

-1

O

-2
0

0.2

0.4
0.6
0.8
Composition, x

1

Figure 2 – 5. Determination method of interatomic potential in this work.

28 / 145§2 – 5. 古典 MD シミュレーションの計算方法
Na2O − SiO2 系ガラスの MD シミュレーションは MXDORTO システム [40] を使用した。長距
離 Coulomb 力の計算には Ewald 法 [41] を用いた。三次元周期境界条件を用いた。MD シミュレーションでは、熱力学的な状態関数三つを指定したうえで系が取り得るあらゆる状態の集合、すなわちアンサンブルを設定する。従来の研究 [34, 42, 43] は実験によりガラスの密度が既知の場
合、粒子 (N)、体積 (V)及び温度 (T) を一定とするシミュレーションを行うことが多い。一方、本研究では未知の組成や高温、高圧下における構造、熱処理過程のシミュレーションを行うことを想定している。そのため、体積が不明の系や一定圧力下で体積が変化する過程の解析に適してい
る、温度 (T)と圧力 (P) 一定の NPT アンサンブルを使用した。速度ベルレ法を用い、時間刻みを
2.0 fs とした。
Figure 2 − 6 にガラス構造を得るまでの冷却過程を示す。それぞれの組成において、一片が約 4.5
nm の立方体の MD セルの中に疑似乱数を用いて 4950 の粒子数 (N) を配置した。Table 2 – 3 に粒子数の詳細を示す。次に、3000 K (T) において 0.1 MPa (P) 下で 250 000 step かけて構造を緩和した。このとき、緩和が充分であることを内部エネルギーが安定になることで確認した。ガラス構造を得るために冷却速度を 0.01 K/step として、段階的に冷却を行った。3000 K – 2100 K では
250 000 step の緩和を行い、1800 K – 300 K では 500 000 step の緩和計算を行った。1800 K – 300
K では融体からガラスになるとき、イオンの動きが遅くなることを考慮し、構造緩和の時間を増加した。ガラス構造を調査するために、10 000 step の解析計算を 300 K, 600 K で行った。
本研究では二つの異なる原子間相互作用をガラスの MD シミュレーションに適用した。（１）
本研究で第一原理計算に基づき設定した原子間相互作用, IP1 と（２） 従来の研究 [32, 35, 36] で経験論に基づき設定された原子間相互作用, IP2 を適用した。IP1 と IP2 を適用したシミュレーシ
ョンを行い、ガラス構造や物性の再現性について比較し検討を行った。

Simulation time step / step
0

1000000 2000000 3000000 4000000 5000000

3500
Temperature, T / K

3000
2500

2000
1500
1000

500
0

0

2000

4000
6000
8000
Simulation time, t / ps

10000

Figure 2 – 6. Cooling process of xNa2O-(1 – x)SiO2 glass simulation in this work.

29 / 145

Table 2 – 3. Number of particles (N) in xNa2O − (1 – x )SiO2 glass simulations

x
0
0.1
0.2
0.3
0.4
0.5

O
3300
3135
2970
2805
2640
2475

Si
1650
1485
1320
1155
990
825

30 / 145

Na
0
330
660
990
1320
1650

Total (N )
4950
4950
4950
4950
4950
4950

§2 – 6. 結果
Table 2 – 4 に xNa2O-(1 – x)SiO2 ガラスの MD シミュレーションから求めた第一近接原子間距
離 (dSi – O, dNa – O) の解析結果を示す。このとき、第一近接原子間距離は二体相関関数 (Paircorrelation function, PCF) のトップピークの値を使用した。IP1 を適用したシミュレーションの dSiO はいずれも、x = 0, 0.2, 0.3, 0.4, 0.5のとき XRD [4] , EXAFS [5 − 7] 測定結果、結晶構造データ[16,23, 24] の値と一致した。
IP1 を適用したシミュレーションにおいて、
x = 0.1 ~ 0.5 の dNa-O は EXAFS[5 − 7] や結晶構造データ [23, 24] で報告されている値と一致した。すべての組成で IP2 を適用したシミュレーションは IP1 よりも dSi-O, dNa-O が短くなった。

Table 2 – 4. Si – O and Na – O bond lengths in xNa2O − (1 – x )SiO2 glasses estimated from the paircorrelation function simulated at 300 K.
Distance / nm

x
0

IP1

IP2

Referencesa0.162(g) / 0.1608(g)b / 0.159(c)cd Si-O

0.164

0.157

d Si-O

0.164

0.157

d Na-O

0.225

0.221

0.232(g)dd Si-O

0.164

0.157

0.1617(g)bd Na-O

0.226

0.225

d Si-O

0.164

0.157

0.1668(g)b / 0.163(g)e †/ 0.156-0.163(c)f†d Na-O

0.230

0.227

0.235(g) / 0.228-0.237(c)d Si-O

0.164

0.158

0.1586(g)bd Na-O

0.231

0.225

d Si-O

0.164

0.158

0.155-0.167(c)gd Na-O

0.232

0.223

0.224-0.244(c)g

0.1

0.2

0.3

e†

0.4

0.5
(g) : glass (c) : crystalab†x = 0.33

cref. [4] ref . [5] ref . [16] d ref . [6] e ref . [7] f ref . [23] g ref . [24]31 / 145f†

Figure 2 − 7 に xNa2O − (1 − x)SiO2 ガラスのモル体積の結果を示す。IP1 を適用したシミュレー
ションのモル体積は、対象とした組成範囲において実験データ [44, 45] よりも大きい値を取ったが、x の増加に伴う減少傾向は再現した。IP2 を適用したシミュレーションのモル体積の値はす
べての x において、実験データ [44, 45] の値と一致した。

40

IP1
IP2
Takahashi et al.[44]Narottam et al.[45]

Molar volume, V / cm3 mol-1

38
35
33
30

28
25
23
20
0

0.1

0.2
0.3
Composition, x

0.4

0.5

Figure 2 – 7. Molar volumes of xNa2O-(1 – x)SiO2 glasses simulated at 300 K.

32 / 145

xNa2O − (1 − x)SiO2 ガラスの熱膨張係数を 300 K と 600 K のモル体積の値を用いて式 (2 − 5)から算出した。
𝛽=1𝑑𝑉(2 – 5)( )𝑉300K 𝑑𝑇V はモル体積、T は温度を表しており、V300K は 300K のときのモル体積である。Figure 2 − 8に xNa2O − (1 − x)SiO2 ガラスの熱膨張係数を示す。 IP1 と IP2 を適用したシミュレーションのThermal expansion coefficient, β / ppm K-1

300 K − 600 K 間の熱膨張係数は、x が増加するに伴い値が増加する傾向 [46] を再現した。

80
IP1
IP2
Knoche et al.[46]

70
60
50
40
30
20
10
0
0

0.1

0.2

0.3

0.4

0.5

Composition, x

Figure 2 – 8. Thermal expansion coefficients of xNa2O-(1 – x )SiO2 glasses obtained from simulations at 300
and 600 K.

33 / 145

Figure 2 − 9 に xNa2O − (1 − x)SiO2 ガラスシミュレーションから得られた O − Si − O 角 (θO – Si –O)と Si – O – Si 角(θSi – O – Si) の分布を示す。 解析に用いたカットオフ距離は 0.2 nm とした。IP1,

IP2 を適用したすべての x のシミュレーションにおいて、 θO – Si – O は約 109 ° を中心とする分布を示した。また θO – Si – O 分布の範囲は x 増加に伴い小さくなり、これはガラス中の SiO4 四面体の歪みが減少していることを示唆している。IP1 を適用したシミュレーションの θSi – O – Si 分布はすべての組成範囲で 120° ~ 180°であった。一方で、IP2 を適用したシミュレーションの θSi – O – Si
分布はすべての組成で約 110° ~ 180°の範囲となり、IP1 より分布範囲が広かった。また、x = 0 のとき、IP1 を適用したシミュレーションの θSi – O – Si 分布は、約 140°と 160°に二つのピーク分布を示し、x が増加するに伴い 140°のピークは増加し 160°のピークが減少した。

34 / 145

(a) IP1

Distribution

θO-Si-O

θSi-O-Si

x=0
0.1

0.2
0.3
0.4
0.5

80

100

120
140
160
Angle, θ / deg.

180(b) IP2
θSi-O-Si

Distribution

θO-Si-O

x=0
0.1

0.2
0.3
0.4
0.5

80

100

120
140
160
Angle, θ / deg.

180

Figure 2 – 9. Distribution of O – Si – O (θO – Si – O ) and Si – O – Si (θSi – O – Si ) bond angles in xNa2O-(1 –x )SiO2 glasses simulated at 300 K using IP1 (a) and IP2 (b). Dashed and solid lines show θO – Si – O and θ
Si – O – Si, respectively.

35 / 145一般に、SiO4 四面体をユニット中の架橋酸素数の数で分類し全 SiO4 四面体に占める割合を Qn比という [47]。Qn 比の概念図を Figure 2 – 10 に示す。Figure 2 – 11 に MD シミュレーションで得た Qn 比を示す。このとき、解析に用いたカットオフ距離は 0.2 nm とした。IP1, IP2 はともに x= 0, 0.1 において Si が O と 3 配位する構造が 1, 2 個出現したが、これはシミュレーションの時
間・空間的制約により生じるガラス構造の歪みもしくは、シンプルな二体間の相互作用のみを適用していることによる構造の歪みと考え、それらを除外して Qn 比を求めた。Maekawa ら [48] による Na2O – SiO2 系ガラスの 29Si MAS − NMR による報告では 0.2 ≤ x ≤ 0.4 において Q3 と Q2が増加し、Q4 が減少することが報告されており、IP1 と IP2 によるシミュレーションはこの結果と一致した。0.4 ≤ x ≤ 0.5 において IP1 は Q1, Q2 が増加し Q3 と Q4 が減少する実験結果の傾向を再現し、IP2 よりも Qn 比の値に改善が見られた。

Network modifier
Bridging oxygen (BO)Q3

Q2

Non-bridging oxygen (NBO)

Q0
Silicon, Si

6-membered ring

Q4

Q1

4-membered rings

5-membered rings

Figure 2 – 10. Qn ratio of Na2O – SiO2 system glasses.

36 / 145

100
IP1

IP2

NM R[48]

90
Q4

80

Q3
Q2

70
Qn (%)

60
50

Q3
Q1

40

30

Q3

Q4

20
Q2

10

Q0

0
0.0

0.1

0.2
0.3
0.4
Composition, x

0.5

0.6

Figure 2 – 11. Qn ratios in xNa2O-(1 – x)SiO2 glasses simulated at 300 K (circles = Q4, triangles = Q3,
squares = Q2, diamonds = Q1, and plus signs = Q0).

37 / 145§2 – 7. 考察
Figure 2 − 12 に IP1 と IP2 の Si – O と Na – O の Coulomb 力項の zizj 成分を示す。Figure 2 – 3に示したように x の増加に伴いアルカリ金属イオンとケイ素イオンの電荷は減少する。Figure 2 –12 より、いずれの x についても IP1 の Si – O 間の Coulomb 力が IP2 よりも弱く働くため、dSi –Oが長くなったと考えられる。また、Si – O 間と同様に IP1 の Na – O 間の Coulomb 力も IP2 より弱く働くため、dNa – O が長くなったと考えられる。

(a) Si-O
-2.78
-2.80

IP1

zSizO

-2.82
-2.84
-2.86

IP2
-2.88
-2.90
0

0.1

0.2
0.3
Composition, x

0.4

0.5

0.4

0.5

(b) Na-O
-1.00

-1.05
-1.10

IP1

-1.15
zNa zO

-1.20

-1.25
-1.30

IP2

-1.35
-1.40

-1.45
-1.50
0

0.1

0.2
0.3
Composition, x

Figure 2 – 12. Si – O, (a) and Na – O, (b) of zizj in Coulomb force (circles = IP1 and triangles = IP2).

38 / 145

Figure 2 – 4 より、電子相関手法を MP2 としているのは、他の手法を用いると MD シミュレー
ションの結果が収束しなかったためである。現時点で MP2 が MD シミュレーションと相性が良
いのかはわからない。MD の原子間相互作用の Coulomb 力候は原子間距離の逆数に対応して作用するため、他の力の項に比べ原子間の相互作用に影響する度合いが大きい。有機分子による MP2法を用いた解析では Coulomb 相互作用を他の電子相関手法に比べ強く表すことが知られていることから [49]、MP2 法を用いることでガラス構造の再現性が良いと考えている。
Figure 2 – 7 より、IP2 は実験データのモル体積の値 [44, 45] と一致するように経験的に設定した原子間相互作用を適用している。しかし、Table 2 – 4 に示したように IP2 の Si – O 間距離は実際よりも短いため、NPT アンサンブル下でのネットワーク構造 (θSi – O – Si 分布, Qn 比)に歪みが生じている恐れがある。本論文では二体間に働く相互作用のみで、ガラスのシミュレーションを行っている。そのため、シミュレーションで得られた物性値と実験値との間には少なからず差が生じると考えており、絶対量ではなく相対的な x の変化に意味があると解釈している。IP2 によるシミュレーションでは、実験報告の室温のモル体積の値となるように原子間相互作用を設定していたため、ネットワーク構造に歪みが生じる結果を示したと考えている。一方、IP1 はモル体積の値が一致するように原子間相互作用を設定しているのではなく、第一原理計算の結果に基づいて設定していることを考えると、x 毎のモル体積を良く再現していると考えている。
Figure 2 – 8 より、0 < x < 0.3 において、IP1 及び IP2 によるシミュレーションでは実験報告[46]に比べ、熱膨張係数が小さい値となった。これは、実際のガラスよりも硬いガラスとしてシミュ
レーションされていると考えられる。実際の xNa2O-(1 – x )SiO2 ガラスにおいて、SiO2 rich な組成
では体積弾性率が大きく、硬い構造であることが知られている [44]。IP1 及び IP2 によるシミュ
レーションでは、シンプルな原子間相互作用モデルのみで取り扱っていることから、構造の硬さを過剰に表現した形で熱膨張係数に現れたものと解釈している。0.3 ≤ x では、IP1 は IP2 よりも実験報告の熱膨張係数の値と一致していた。これは、IP1 の方が IP2 によるシミュレーションに比べ、ガラス中のネットワーク構造に歪みを生じる割合が少ないことと関係していると考えており、熱膨張係数の傾向にネットワーク構造の影響が現れたと考えている。
Mozzi ら [4] は X 線構造解析の結果から x = 0 における θSi – O – Si が 120° ~ 180°に分布することを報告している。Farnan ら [50] は 0.2K2O – 0.8SiO2 ガラス中の θSi – O – Si は 120° − 130 °～180°に分布すると 17O DAS – NMR 測定結果から報告している。Figure 2 – 9 より IP1 は θSi – O – Si が120° ～ 180° の分布を示しており、IP2 よりも実験報告と近い結果となった。
Figure 2 – 13 に MD シミュレーションで得た Na2O – SiO2 ガラスの Si と O から構成されるリ
ングサイズの分布を示す。リングサイズはリング中の Si の数で定義している。IP1 を適用したシ
ミュレーションでは、x = 0 のとき 7, 8 員環が多く、x が増加するに伴い 4, 5 員環が増加した。これはアルカリ酸化物量の増加によってネットワーク構造が切断されていることを反映していると考えられる。Figure 2 – 9 に示した θSi – O – Si の分布にもネットワーク構造の変化の影響が反映されていると考えられる。140°のピークは 3 ~ 5 員環の小さいリング構造を、160°のピークは 6 ~8 員環の大きいリング構造を構成している角度を示すと考えられる。

39 / 145

(a) IP1

Number of rings per 100 silicon atoms (counts)

120

100

0

0.1

0.2

0.3

0.4

0.5

80

60

40

20

0

0

1

2

3

4

5

6

7

8

5

6

7

8

Ring size(b) IP2

Number of rings per 100 silicon atoms (counts)

120

100

0

0.1

0.2

0.3

0.4

0.5

80

60

40

20

0
0

1

2

3

4

Ring size

Figure 2 – 13. Ring size distributions in xNa2O-(1 – x)SiO2 glasses simulated at 300 K using IP1 (a) andIP2 (b).

40 / 145

Na2O-SiO2 系ガラスの XANES の SiK–edge から得られたエネルギーシフトの傾向 [5] より、x増加に伴い最近接の Si – Si 距離が短くなることが報告されており、θSi – O – Si の分布が高角度側か
ら低角度側へシフトすることが予想される。Figure 2 – 9 からも x の増加に伴い θSi – O – S i のピー
ク分布が低角度側へのシフトがみられる。一方で、IP2 による MD シミュレーションでは x = 0.4,
0.5 のとき 3 員環が IP1 に比べて多くなっており、これはモル体積を実験データに合わせたことによるネットワーク構造の歪みの表れと考えている。
Figure 2 – 11 に示す IP1, IP2 による Qn 比の値は実験値と一致しなかったが、IP1 は IP2 より
も 0 ≤ x ≤ 0.5 の範囲で Qn 比の増減傾向の再現性が改善していることが分かった。これは、IP1が従来の IP2 のようにすべての x に対し同一の原子間相互作用を適用する方法ではなく、x 毎に原子間相互作用を設定していることで、ネットワーク構造の歪みを解消できたことに繋がっていると考えている。

41 / 145§2 – 8. 結言
第一原理計算の結果により定めた原子間相互作用を適用した MD シミュレ―ションを Na2O –SiO2 系ガラスを対象に行った。原子間相互作用に用いる電荷は組成毎にアルカリケイ酸塩結晶のMulliken population analysis の結果を用いて設定した。Si – O 間の原子間相互作用は SiO2+ モデルの分子軌道計算から求めた。この原子間相互作用を適用することにより、原子間距離、熱膨張係数、モル体積、Si – O – Si 結合角、Qn 比の組成変化の実験値の傾向を従来の経験論に基づく設定
方法によるシミュレーションと同等のシミュレーション結果を示すことが分かった。

42 / 145§2 – 9. 参考文献
[1] S. Sakka, K. Matsusita, Studies on binary silicate glasses based on the SiKα and SiKβ emission X-rays,
J. Non-Cryst. Solids 22 (1976) 57 – 66.
[2] J. Gasteiger and M. Marsili, Iterative partial equalization of orbital electronegativity- a rapid access to
atomic charges, Tetrahedron 36 (1980) 3219 – 3288.
[3] A. K. Rappe and W. A. Goddard Ⅲ, Charge equillbration for molecular dynamics simulations, J. Phys.
Chem. 95 (1991) 3358 – 3363.
[4] R.L. Mozzi and B.E. Warren, The structure of vitreous silica, J. Appl. Cryst. 2 (1969) 164 – 172.
[5] G.S. Henderson, A Si K-edge EXAFS/XANES study of sodium silicate glasses, J. Non-Cryst. Solids
183 (1995) 43 – 50.
[6] C. Mazzara, J. Jupille, A.-M. Flank, P. Lagarde, Stereochemical order around sodium in amorphoussilica, J. Appl. Cryst. B 104 (2000) 3438 – 3445.
[7] G.N. Greaves, A. Fontaine, P. Lagarde, D. Raoux, S.J. Gurman, Local structure, of silicate glasses,
Nature 293 (1981) 611 – 616.
[8] F. G. Fumi and M. P. Tosi, Ionic sizes and Born repulsive parameters in the NaCl-type alkali halides. I.,
The Huggins-Mayer and Pauling forms, J. Phys. Chem. Solids, 25 (1964) 31 – 43.
[9] L. Pauling, M. Koizumi (Eds.), The nature of the chemical bond, Cornel University Press, (1960).
[10] K. Kawamura and N. Sawaguchi, Molecular dynamics calculations of mixing of silicate melts, Bull.
Volcanol. Soc. Jpn. 42 (1997) S1 – S11.
[11] Y. Xiang, J. Du, M. M. Smedskjaer, J. C. Mauro, Structure and properties of sodium aluminosilicate
glasses from molecular dynamics simulations, J. Chem. Phys. 139 (2013) 044507.
[12] H. Sakuma, K. Kawamura, Structure and dynamics of water on muscovite mica surfaces, Geochim.
Cosmochim. Acta 73 (2009) 4100 – 4110.
[13] S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.J. Probert, K. Refson, M.C. Payne, First principles
methods using CASTEP, Zeitschrift für Kristallographie, 220 (2005) 567 – 570.
[14] D. R. Hamann, M. Schluter, C. Chiang, Norm-conserving pseudopotentials, Phys. Rev. Lett. 43 (1979)
1494 – 1497.
[15] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev.
Lett. 77 (1996) 3865 – 3868.
[16] A.F. Wright and M.S. Lehmann, The structure of quartz at 25 and 590 ℃ determined by neutron
diffraction, J. Solid State Chem. 36 (1981) 371 – 380.
[17] H. Krüger, V. Kahlenberg, R. Kaindl, Li2Si3O7: Crystal structure and Raman spectroscopy, J. Solid
State Chem. 180 (2007) 922 – 928.
[18] B.H. de Jong, H.T. Supèr, A.L. Spek, N. Veldman, G. Nachtegaal, J.C. Fischer, Mixed alkali systems:
Structure and 29Si MASNMR of Li2Si2O5 and K2Si2O5, Acta Cryst. B54 (1998) 568 – 577.
[19] K.F. Hesse, Refinement of the crystal structure of lithium polysilicate, Acta Cryst. B33 (1977) 901 –
902.
[20] H. Voellenkle, A. Wittmann, H.N. Nowotny, Die kristallstruktur der verbindung Li6(Si2O7), Monatsh.
43 / 145

Chem. 100 (1969) 295 – 303.
[21]

R. Hoppe, K. Bernet, A. Moeller, Einkristall-Synthese hochschmelzender Oxyde bei niederer
Temperatur: Gamma-(Li4(SiO4))-ohne Fehlordnung, isotyp mit Na4(SnO4). (Was heisst 'Isotypie'?)Zeitschrift fuer Anorganische und Allgemenie Chemie, 629 (2003) 1285 – 1293.

[22] R. Hofmann and R. Hoppe, Ein neues Oxogermanat: Li8GeO6 = Li8O2(GeO4) (Mit einer Bemerkung
ueber Li8SiO6 und Li4GeO4), Zeitschrift fuer Anorganische und Allgemeine Chemie, 555 (1987) 118– 128.
[23] A.K. Pant and D.W.J. Cruickshank, The crystal structure of alpha-Na2Si2O5, Acta Cryst. B24 (1968)
13 – 19.
[24] A. Grund and M.M. Pizy, Structure crystalline du metasilicate de sodium anhydre, Na2SiO3, Acta
Cryst. 5 (1952) 837 – 840.
[25] W.H. Baur, E. Halwax, H. Voellenkle, Comparison of the crystal structures of sodium orthosilicate,
Na4SiO4, and sodium orthogermanate Na4GeO4, Monatsh. Chem. 117 (1986) 793 – 797.
[26] M. Jansen, Zur kristallstruktur von kaliumdisilicat, Zeitschrift fuer Kristallographie, 160 (1982) 127 –
133.
[27]

K. Bernet and R. Hoppe, Zur Kristallstruktur von K4(SiO4), Zeitschrift fuer Anorganische undAllgemeine Chemie 589 (1990) 129 – 138.

[28] E. Zintl, A. Harder, B. Dauth, Gitterstruktur der oxide, sulfide, selenide und telluride des lithiums,
natriums und kaliums, Zeitschrift fuer Elektrochemie 40 (1934) 588 – 593.
[29] R.S. Mulliken, Electronic population analysis on LCAO-MO molecular wave functions. I, J. Chem.
Phys. 23 (1955) 1833 – 1840.
[30] M.D. Segall, R. Shah, C.J. Pickard, M.C. Payne, Population analysis of plane-wave electronic structure
calculations of bulk materials, Phys. Rev. B 54 (1996) 16317 – 16320.
[31]

Y. Yamamoto, N. Sawaguchi, M. Sasaki, Interatomic potential for xNa2O-(1–x)SiO2 glasses, J.
Comput. Chem., Jpn. 14 (2015) 63 – 66.

[32] Y. Yamamoto, N. Sawaguchi, M. Sasaki, A new determination method of interatomic potential for
sodium silicate glass simulations, J. Non-Cryst. Solids 466 – 467 (2017) 29 – 36.
[33] F. L. Hirshfeld, Bonded-atom fragments for describing molecular charge densities, Theoret. Chim.
Acta 44 (1977) 129 – 138.
[34] A.N. Cormack, J. Du, T.R. Zeitler, Sodium ion migration mechanisms in silicate glasses probed bymolecular dynamics simulations, J. Non-Cryst. Solids 323 (2003) 147 – 154.
[35] F. Noritake, K. Kawamura, T. Yoshino, E. Takahashi, Molecular dynamics simulation and electricalconductivity measurement of Na2O・3SiO2 melt under high pressure; relationship between its structure
and properties, J. Non-Cryst. Solids 358 (2012) 3109 – 3118.
[36] N. Sawaguchi, K. Yamaguchi, M. Sasaki, K. Kawamura, Interatomic potential model for molecular
dynamics simulation of lithium borate melts/glasses, J. Comput. Chem. Jpn. 14 (2015) 139 – 146.
[37] Gaussian 09, Revision E.01, M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R.
Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li,
44 / 145

H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota,
R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A.
Montgomery, Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N.
Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J.
Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J.
Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski,
R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich,
A.D. Daniels, Ö. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian, Inc., WallingfordCT, 2009.
[38] M. Head-Gordon, J.A. Pople, M.J. Frisch, MP2 energy evaluation by direct methods, Chem. Phys.
Lett. 153 (1988) 503 – 506.
[39] G.A. Petersson, A. Bennett, T.G. Tensfeldt, M.A. Al-Laham, W.A. Shirley, J. Mantzaris, A complete
basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row
atoms, J. Chem. Phys. 89 (1988) 2193 – 2218.
[40] K. Kawamura, MXDORTO, Japan Chemistry Program Exchange P029, (1996).
[41] P. P. Ewald, Evaluation of optical and electrostatic lattice potentials, Ann. Phys. (Leipzig) 64 (1921)
253 – 287.
[42]

J. Du and A.N. Cormack, The medium range structure of sodium silicate glasses: A molecular
dynamics simulation, J. Non-Cryst. Solids 349 (2004) 66 – 79.

[43] R.N. Mead and G. Mountjoy, A molecular dynamics study of the atomic structure of (CaO)x(SiO2)1–
x glasses, J. Phys. Chem. B 110 (2006) 14273 – 14278.
[44] K. Takahashi and A. Osaka, Elastic properties of alkali silicate glasses, Yogyo-Kyokai-Shi, 91 (1983)
22 – 26.
[45] P.B. Narottam, R.H. Doremus, Handbook of glass properties, Acdemic Orlando, (1986).
[46] R. Knoche, D.B. Dingwell, F.A. Seifert, S.L. Webb, Non-linear properties of supercooled liquids inthe system Na2O-SiO2, Chem. Geol. 116 (1994) 1 – 16.
[47] J.F. Stebbins, Identification of multiple structural species in silicate glasses by 29Si NMR, Nature 330(1987) 465 – 467.
[48] H. Maekawa, T. Maekawa, K. Kawamura, T. Yokokawa, The structural groups of alkali silicate glasses
determined from 29Si MAS-NMR, J. Non-Cryst. Solids 127 (1991) 53 – 64.
[49] S. Tsuzuki, Yuukibunshi no bunshikanryoku: Ab initio bunshikidouhou ni yoru bunshikansougosayou
enerugi no kaiseki, University of Tokyo Press, (2015).
[50]

I. Farnan, P.J. Grandinetti, J.H. Baltisberger, J.F. Stebbins, U. Werner, M.A. Eastman, A. Pines,
Quantification of the disorder in network-modified silicate glasses, Nature 358 (1992) 31 – 35.

45 / 145

第 3 章 0.25{xK2O - (1 – x)Na2O}-0.75SiO2 ガラスシミュレーションの混合アルカリ効果の再現性§3 – 1. 緒言
実用ガラスの主成分は SiO2 が多く、その中でも混合アルカリケイ酸塩ガラスは実用的な組成としてよく使用される。混合アルカリガラスは、二種類以上のアルカリ金属酸化物を導入することで、アルカリ金属酸化物イオンの移動度を減少させる現象、いわゆる混合アルカリ効果が起きる。混合アルカリ効果によって、ガラスの粘性や電気抵抗を増加させることが報告されて以来、混合アルカリ効果のメカニズムを理解することは新たなガラス材料の開発に役立つとして、多くの研究が行われてきた [1 – 4]。
放射性同位体を用いた実験やガラス中の内部摩擦を求める実験では、異なるアルカリ金属イオ
ンが協調運動を行うことで混合アルカリ効果を引き起こしていると解釈されている [5 – 11]。X
線吸収 [12, 13]、NMR [14, 15]、赤外分光測定 [16, 17] などからも、混合アルカリ効果に関する研究が行われている。一方、逆モンテカルロシミュレーション [18] や MD シミュレーションを用いた研究も行われてきた。MD シミュレーションでは、diffusion controlled relaxation モデル [19]や jump diffusion モデル [20]、動的構造モデル [21 – 23] が提案されるなど、数多くの研究がなされてきた。また、混合アルカリ効果の起源についての研究も行われており、Balasubramanian と Rao[24]、 Habasaki ら [25] は、Li2O – K2O – SiO2 系ガラスにおいて、混合アルカリ効果はガラス中の異種アルカリ金属イオンのサイズ差のみによってアルカリ金属イオン同士の拡散が阻害される効果として報告している。このように MD シミュレーションはケイ酸塩系の混合アルカリ効果のメカニズムの解釈に活用されてきたが、ガラス中の局所的な構造のシミュレーションの解析によって、これらのことを結論づけている点がある。混合アルカリ効果のメカニズムを解釈するためには、系全体におけるアルカリ金属イオンの拡散挙動から考察していく必要があると考えている。また、従来の MD シミュ

レーションで用いられている原子間相互作用は、実験結果に一致するように試行錯誤的にポテン
シャルパラメータを設定することが多いが [26 – 28]、経験的な要素を多く含むため、計算結果の信頼性の面で疑問視されることが少なからず存在する。Habasaki らは Figure 1 – 6 に示した方法
によって第一原理計算から原子間相互作用を決定しているが [29]、Coulomb 力項の各イオンの電荷もフィッティングの考慮としているため、アルカリ金属イオン種の組成変化による電荷の傾向を、
原子間相互作用に反映していない背景をもつ。そのため、混合アルカリケイ酸塩ガラスに対し、原子間相互作用の設定方法を見直し、再度 MD シミュレーションを行う必要があると考えている。
第 3 章では、0.25{xK2O-(1−x)Na2O}-0.75SiO2 ガラス(0 ≤ x ≤ 1) を対象に、第 2 章で考案した第一原理計算の結果から設定した原子間相互作用を適用したシミュレーションを行い、混合アルカリ効果を再現できるか、また混合アルカリ効果のメカニズムについて新たな知見が得られるかを調査した。混合アルカリ効果はアルカリ金属イオンのサイズ差が開くほど、その効果が大きく働くとされている一方で、サイズ差が小さいとその効果も小さいという報告がある [3, 4]。MD シミュレーショ
ンでは、混合アルカリ効果の影響が明確に表れやすい Li2O – K2O – SiO2 系ガラスのシミュレーションの適用例が多いが [25, 26]、Na2O – K2O – SiO2 系のシミュレーションの事例は少ない。そのため、Na2O– K2O – SiO2 系の計算を行うことで、本研究の MD シミュレーションでも僅かな影響しか示さないイ
オンのダイナミクスの再現性を調査するために適した組成であると考えている。
§3 – 2 では、適用した原子間相互作用の設定方法について記した。§3 – 3 では、MD シミュレー

46 / 145

ションの計算条件について示した。§3 – 4 では MD シミュレーションを行った結果について示した。
§3 – 5 では§3 – 4 で得られた結果について考察した。§3 – 6 では本章の結言を示した。§3 – 7 に本章で引用した参考文献を示す。

47 / 145§3 – 2. 原子間相互作用
原子間相互作用は第 2 章で示した方法により決定した。原子間相互作用は式 (2 – 1) を使用し、電荷は第 2 章§2 – 3 の計算条件を用いて、KxNa2−xSi2O5 結晶 [30 – 35] に対し電荷解析を行った。
対象とした結晶を Table 3 – 1 に示し、Figure 3 – 1 に KxNa2−xSi2O5 結晶の電荷解析結果を示す。ア
ルカリ金属イオンの交換により、すべてのイオンの電荷が一定の値を示した。この電荷解析の傾向を基に y{xK2O-(1−x)Na2O}-(1−y)SiO2 系中のイオンの電荷, zNa, zK, zSi, zO を以下の式 (3 – 1) - (3 – 4)に示すように定式化させた。
𝑧𝑁𝑎 = −0.25𝑦 + 1.0 (0 ≤ 𝑥 < 1, 0 < 𝑦 ≤ 1)

(3 – 1)𝑧𝐾 = −0.25𝑦 + 1.0 (0 < 𝑥 ≤ 1, 0 < 𝑦 ≤ 1)

(3 – 2)𝑧𝑆𝑖 = −0.50𝑦 + 2.4 (0 ≤ 𝑥 ≤ 1, 0 ≤ 𝑦 < 1)

(3 – 3)𝑧𝑂 = neutrality(3 – 4)zO に系の電気的中性を保たせた。式(3 – 1) – (3 – 4) は Figure 2 – 3 の結果も取り入れ、zNa, zK, zSi
は y に対してのみ変化し、アルカリ金属酸化物の混合量, x に依らない式とした。Si – O 間のポテ
ンシャルカーブ形状は第 2 章§2 – 4 と同様のものを使用した。aNa, bNa, cNa, aK, bK 及び cK は経験的に決めた値を使用した。設定した原子間相互作用を Table 3 – 2 に示す。

48 / 145

49 / 145

Table 3 – 1. Percent deviation, which gives the difference in structural optimization by DFT. Deviation shows in absolute value.

4

M ulliken charge

3
2

Si

1

K (◇)
or Na (○)

0
-1

O

-2
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 3 – 1. Mulliken charges of ions in several 0.33{xK2O-(1 – x )Na2O}-0.67SiO2 crystals (listed
in Table 3 – 1). Squares, diamons, circles and triangles show Si, K, Na and O ions.

Table 3 – 2. Interatomic potential parameters for 0.25{xK2O-(1 – x )Na2O}-0.75SiO2 glasses

39

1/2

z

a / nmb / nmO-1.24200.18384

0.01393c / (kJ / mol)3
nm0.05605Si

2.27500.993220.007970.00000

Na

0.93750.132200.01150

0.01227

0.93750.17500

0.01450

0.02864

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1627600

52.06

83680

34.18

K,

133

Si-O

K

50 / 145§3 – 3. MD シミュレーションの計算条件
MD シミュレーションは MXDORTO システム [36] を使用し、時間刻みは 2.0 fs とした。周期境界条件は三次元周期境界を用いた。NPT アンサンブル下でシミュレーションを行い、一辺が約5.0 nm の立方体 MD セル中に疑似乱数を用いて粒子(N = 5400) をランダムに配置させ、圧力(P =0.1 MPa ) および温度(T = 3000 K) の融体のシミュレーションを行い、250 000 step にかけて全エ
ネルギーが一定になるまで構造緩和を行った。Table 3 – 3 に本シミュレーションで扱った粒子数を示す。ガラス状態は冷却速度を 0.01 K/step とし、段階的に冷却することで得た。3000–2100 K
では 250 000 steps 、ガラス状態に近づくほどイオンの拡散が遅くなることを考慮し、1800–300 K
では 500 000 steps の構造緩和を行った。300 K から 1200 K にかけて 300 000 step の 50 step 間隔でデータを記録し解析計算を行った。アルカリ金属イオンの自己拡散係数は以下の式 (3 – 5) を用いた [37]。

1

𝐷 = lim 6𝑡 〈[𝑟(𝑡 + 𝑡0 ) − 𝑟(𝑡0 )]2 〉(3 – 5)𝑡→∞r(t) は時間 t におけるイオンの位置、 t0 は基準時間、<[r(t+t0) − r(t0)]2> は平均二乗変位を示す。
従来の MD シミュレーションの報告では、アルカリ金属イオン間の質量差は混合アルカリ効果に関係はないという報告がある [25 ]。しかし、本研究では、従来の研究とは異なる原子間相互作用の設定方
法を用いていることを考慮し、質量効果についても改めて調査することにした。そこで、Na-39K 系に加え、仮想的な系として Na-133K 系と 39K133KSi3O7 系のシミュレーションも行った。133K に関連する原子間相互作用のポテンシャルパラメータは 39K と同様のものを適用した。

Table 3 – 3. Number of particles (N) in 0.25{xK2O-(1 – x )Na2O}-0.75SiO2 glass simulations.
x
0
0.2
0.4
0.5
0.6
0.8
1

O
3150
3150
3150
3150
3150
3150
3150

Si
1350
1350
1350
1350
1350
1350
1350

K
0
180
360
450
540
720
900

51 / 145

Na
900
720
540
450
360
180
0

Total (N )
5400
5400
5400
5400
5400
5400
5400

§3 – 4. 結果
Figure 3 – 2 に 300 K と 900 K の 0.25{xK2O-(1−x)Na2O}-0.75SiO2 ガラスのモル体積を示す。300
K のモル体積の値は実験値 [38] と異なるが、x の増加に伴ってモル体積が増加する傾向は一致した。

38

M olar volume, V / cm3 mol-1

36This work(900 K)

34
This work(300 K)

32

30
Exp. by Shelby et al. (300 K) [38]
28
26
24
22
0.00

0.20

0.40
0.60
Composition, x

0.80

1.00

Figure 3 – 2. Molar volumes of 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 300 K and 900 K. Dashedline shows ideal mixing. Open circles, open triangles and solid diamonds show this work at 300 K, this work at
900 K and experimental values [38] at 300 K.

52 / 145

Figure 3 – 3 に 900 K における二体相関関数 (Pair Correlation Function, PCF) を示す。第一近接Si –O、O – O、及び Si – Si 間の PCF はすべての x において、それぞれ 0.145 – 0.205 nm、0.230– 0.300 nm 及び 0.270 – 0.360 nm であり、x による変化は見られなかった。第二近接以遠の Si –O、O – O、及び Si – Si 間の PCF に x に対してわずかな構造変化が見られた。300 K における結晶の第一近接 Si – O、O – O、及び Si – Si 間の PCF はそれぞれ、0.154 – 0.165 nm、0.257 – 0.273
nm 及び 0.301 – 0.320 nm であり、XRD による 300 K におけるガラスの Si – O 間距離は 0.162
nm (SiO2 glass) [39]、0.163 nm (Na2Si2O5 glass) [40]、0.1.61 nm (0.125Li2O – 0.125K2O – 0.75SiO2 glass)[41] であり、これらは本研究結果と一致する。一方、第一近接 Na – O、K – O、Na – Na、K – K
及び Na – K 間の PCF は x が増加するに伴い分布に変化が見られた。300 K における結晶の第一近接 Na – O、K – O、Na – Na、K – K 及び Na – K 間の PCF はそれぞれ、0.228 – 0.297 nm、0.256 – 0.334 nm、0.308 – 0.328 nm、0.315 – 0.395 nm 及び 0.357 – 0.425 nm であり、XRD による
300 K におけるガラスの Na – O 間距離は 0.235 nm (Na2Si2O5 glass) [40]、K – O 間距離は 0.279
nm (0.125Li2O – 0.125K2O – 0.75SiO2 glass) [41] と報告されており、これらは本研究結果と一致する。1200 K 以下のすべての PCF は 900 K と同様の結果となった。

53 / 145

x

Na–K

Na–Na

x

x

K–K

PCF

Si–Si

O–O

Na–O
x
x

K–O

Si–O

0.1

0.2

0.3
0.4
Distance, r ij / nm

0.5

0.6

Figure 3 – 3. Pair correlation function for 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 900
K. x = 0, 0.2, 0.4, 0.5, 0.6, 0.8 and 1 shows black dashed, red, yellow, green, blue, parple and blacksolid line, respectively.

54 / 145

積算配位数(Running Coordination Number, RCN) を式 (3 – 6) より求めた。
𝑅𝐶𝑁𝑖𝑗 (𝑟) =1 𝑟∫ 𝑛 (𝑟 ′ )𝑑𝑟 ′𝑁𝑖 0 𝑖𝑗(3 – 6)nij(r’) は距離, r’ におけるイオン, i の周囲に存在するイオン, j の数を示す。Figure 3 – 4 に 900 Kの RCN を示す。Si – O 間の RCN は 4 配位を、O – O 間の RCN は約 5 配位を示した。Si – Si の
RCN は約 3 配位、Na – O と K – O の RCN は 4 ~ 6 配位であった。これは、17O DAS NMR 測定[42] や従来の MD 計算の報告 [43] と一致していた。第一近接の Si – O、O – O、Si – Si、Na – O
及び K – O の RCN は x によって大きな変化は見られない。0.5 nm 以降において Si – O、Si – Siの RCN に混合アルカリに対するわずかな配位数変化が見られた。一方、Na – Na の RCN は x が増加するに伴い減少し、K – K の RCN は x が増加するに伴い増加した。x が増加するにつれ Na– K の RCN は増加し、K – Na の RCN は減少し、これは Na – Cs 系の 133Cs MAS NMR による実験報告と同様の傾向を示した [14]。

55 / 145

56 / 145

0.2

0.2

0.3

(d) Na–K

0.1

Si–O

x

0.4 0.5 0.6
Distance, rij / nm

x

Na–K

x
Si–Si

0.7

0.3
0.4
0.5
Distance, rij / nm

O–O

0.8

0.6

10
9
8
7
6
5
4
3
2
1
0

10
9
8
7
6
5
4
3
2
1
0
0.3

0.2

0.3

(e) K–Na

0.2

x

Na–O

Na–Na

x
0.4 0.5 0.6 0.7
Distance, rij / nm

K–Na

0.4 0.5 0.6 0.7
Distance, rij / nmx(b) Na–O and Na–Na

0.8

0.8

x

dashed, red, yellow, green, blue, parple and black solid line, respectively.

Figure 3 – 4. Running coordination number for 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 900 K.

10
9
8
7
6
5
4
3
2
1
0

RCN

RCN

(a) Si–O, O–O and Si–Si

RCN
RCN

10
9
8
7
6
5
4
3
2
1
0
0.2

x

0.3

K–O

K–K
0.4 0.5 0.6 0.7
Distance, rij / nmxx(c) K–O and K–K

0.8

x = 0, 0.2, 0.4, 0.5, 0.6, 0.8 and 1 shows black

10
9
8
7
6
5
4
3
2
1
0

RCN

一般的に温度が低いほど、平均二乗変位の直線性が不安定になることが経験的に知られている。
0.25{xK2O - (1 − x)Na2O}-0.75SiO2 ガラスの 700 K におけるアルカリの平均二乗変位 (meansquared displacement, msd.)を Figure 3 – 5 – 1, – 2 に示す。700 K における msd. は 200 ps – 500 ps
間で直線性を示したことから、この領域から自己拡散係数を見積もった。700 K 以上の温度に対しても同様の操作を行った。700 K より低い温度に対しては、msd. が直線性を示さなかったため、自己拡散係数は見積もれなかった。Figure 3 – 6 に Na-39K と Na-133K 系における Na+ イオンとK+ イオンの自己拡散係数, DNa 及び DK を示す。Na-39K と Na-133K 系は、いずれの温度においても異種アルカリ金属イオンを混合することで DNa 及び DK が減少した。また、K の質量が増加すると、K+ イオンの自己拡散係数は減少し、対となる Na+ イオンの自己拡散係数も減少した。

57 / 145

(a) x = 0

15
10

5

25

25

20

20

msd. of Na / nm2

20

15

10
5

200
400
Time, Δt / ps

0600(d) x = 0.5200400Time, Δt / ps20

20

msd. of Na / nm2

25

15

10
5

0

25

15

10
5

600

200
400
Time, Δt / ps

600(f) x = 0.8

20
15

10
5

0

0
200
400
Time, Δt / ps5(e) x = 0.6

25

0

10

600

msd. of Na / nm2

0

15

0

0

0

msd. of Na / nm2(c) x = 0.4(b) x = 0.2

msd. of Na / nm2

msd. of Na / nm2

25

0
0

200
400
Time, Δt / ps

600

0

200
400
Time, Δt / ps600

Figure 3 – 5 – 1. Mean squared displacement of Na+ ion in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses
simulated at 700 K. Solid and dashed line show Na+ ion in Na-39K and Na-133K system, respectively.

58 / 145

30

30

25

20
15
10

35
30

25
20
15
10
5

5
200
400
Time, Δt / ps(d) x = 0.6
35

30

30
msd. of K / nm2

35
25

20
15
10

200
400
Time, Δt / ps

Figure 3 – 5 – 2.

600

200
400
Time, Δt / ps

600(f) x = 1.0
35
30

10
0

10

0

(e) x = 0.8

15

0

15

600

20

5
0

200
400
Time, Δt / ps

25

5

20

0
0

600

msd. of K / nm2

0

25

5

0

0

msd. of K / nm2(c) x = 0.5(b) x = 0.4

msd. of K / nm2

35msd. of K / nm2

msd. of K / nm2(a) x = 0.2
35

25

20
15

10
5
0

0

200
400
Time, Δt / ps

600

0

200
400
Time, Δt / ps600

Mean squared displacement of K+ ion in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses

simulated at 700 K. Solid and dashed line show K+ ion in Na-39K and Na-133K system, respectively.

59 / 145

(a) at 700 K, 900 K and 1100 K(b) at 800 K, 1000 K and 1200 K

-4.5

-4.5

1200 K39-5.01100 KK-5.01000 K
Self-diffusion coefficient, log (D / cm2 s-1)

Self-diffusion coefficient, log (D / cm2 s-1)

133K

900 K
-5.5

-6.0
700 K

-6.5

-7.0

-5.5

800 K
-6.0

-6.5

-7.0

Na+

Na+

K+

K+

-7.5

-7.5
0.00

0.20

0.40
0.60
Composition, x

0.80

1.00

0.00

0.20

0.40
0.60
Composition, x

0.80

1.00

Figure 3 – 6. Self-diffusion coefficients of 39K+, 133K+ and Na+ ions in 0.25{xK2O-(1−x)Na2O}-0.75SiO2
glasses at 700, 900, and 1100 K (a), and at 800, 1000, and 1200 K (b). Dashed lines and solid lines show thediffusion coefficient for Na+ ions and K+ ions, respectively. Solid and open symbols show the results for Na-39K
and Na-133K system, respectively.

60 / 145§3 – 5 . 考察
Figure 3 – 2 に示した 300 K におけるモル体積は文献値 [38] よりも大きい値となったが、x が増加するにつれ、モル体積が増加する傾向 [38] を再現できていることから、相対的な x の変化のシミュレーションは行えていると考えている。シミュレーションによる系のモル体積は異種ア
ルカリ金属イオンを混合することで、破線で示した理想混合線よりも値を下回った。この傾向は
文献値の傾向 [38] を再現しており、Na+イオンと K+ イオンの混合によって、パッキング性が良くなっていることを反映していると考えられる。Na – K 系ガラスの電子照射実験からも Figure 3– 2 に示したような体積変化が起きることが報告されており、アルカリ金属イオンのサイズが混合アルカリ効果との関連性が指摘されている [41]。
ここで、Figure 3 – 7 – 1 及び – 2 に 0.25{xK2O-(1−x)Na2O}-0.75SiO2 ガラスのシミュレーションで得られたエンタルピーと混合エンタルピーを示す。Figure 3 – 7 – 1 より、x が増加するに伴い系のエンタルピーは正の方向に増加する傾向を示した。900 K のエンタルピーが 300 K のエンタルピーよりも正に増加したのは、温度が増加すると粒子の運動が激しくなるためである。Figure 3– 7 – 2 より、300 K 及び 900 K の混合エンタルピーはすべての x において負の値を示しており、エンタルピー的に混合することを示す。Takahashi ら [44] による同ガラス組成の実験報告においても混和傾向を示すことが確認されている。この実験系研究による混合アルカリ効果の熱力学的解釈として、異種アルカリ金属イオン間の相互作用が同種のアルカリ金属イオン間よりも増大することから、混合エンタルピーが負の傾向を示していると考えられている [7, 44]。実験的にも Na+
イオンと K+ イオンの混合系は負に混合エンタルピーを示すことが報告 [45, 46] されていることから、本研究は熱力学的な物性に関してもシミュレーションが行えることを示した。
まず、SiO2 ネットワーク構造が混合アルカリ効果と関連性があるかを調査し、その後アルカリ金属イオン間のダイナミクスと混合アルカリ効果との関係性について調査を行うことにした。

61 / 145

-3940

Enthalpy, H / kJ mol−1-3960This work at 900 K

-3980
-4000
-4020
This work at 300 K

-4040
-4060
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 3 – 7 – 1. Enthalpy of 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 300 K and
900 K. Triangles and circles show this work at 300 K and 900 K, respectively.

62 / 145

Mixing enthalpy, ΔHmix / kJ mol−10Exp. [44] at 300 K-1-2This workat 300 K

This workat 900 K

-3
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 3 – 7 – 2. Enthalpy of mixing in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 300
K and 900 K. Triangles, circles and solid circles show this work at 300 K, this work at 900 K andexperimental data at 300 K [44], respectively.

63 / 145

Figure 3 – 3 に示した、Si–O と O–O の PCF は x に対して、短距離構造である SiO4 四面体ユニット構造に変化がないことを示す。ここで、Qn を Figure 3 – 8 に示す。このとき、解析に用いたカットオフ距離は 0.2 nm とした。x が増加するに伴い、わずかに Q3 が増加し、Q2 と Q4 が減少した。Q0 と Q1 の存在量は約 3 % 以下であった。これらは Maekawa ら [47] が NMR で測定した傾向と一致している。アルカリ金属イオンの交換により、Si – O ネットワーク構造にわずかな変化が見られた。アルカリ種の違いによるネットワーク修飾能が中距離的な構造に影響を与えるかを調べるために、Figure 3 – 9 にリング構造を示した。

60

50

Q4
Q3

Qn (%)

40

30

20
Q2
10

Q0

Q1
0
0.00

0.20

0.40
0.60
Composition, x

0.80

1.00

Figure 3 – 8. Qn ratios in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 900 K (crosses =Q4, diamonds = Q3, squares = Q2, triangles = Q1, and circles = Q0).

64 / 145

Counts of rings per 100 Si atoms

30
3
6

25

4
7

5
8

20

15
10
5

0
0

0.2

0.4

0.6

0.8

1

Composition, x

Figure 3 – 9. Ring size distributions in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses simulated at 900 K.
Circles, triangles, squares, diamonds, crosses and plus signs show 3-, 4-, 5-, 6-, 7- and 8-membered rings,
respectively.

65 / 145

x が増加するに伴い 3 ~ 5 員環が増加し、6 ~ 8 員環が減少した。これは x が増加するほど、Na+
イオンよりも修飾能力が高い K＋イオンが増えるため、大きいリング構造が切断され、小さいリングが多くなったためと考えられる。これは PCF (Figure 3 – 3) と RCN (Figure 3 – 4 ) で示した 0.3
nm 以遠の Si – O, O – O 及び Si – Si 間距離がわずかに変化することに対応した結果を示している。リング構造とモル体積の結果から考えると、Na rich 組成では K rich 組成よりも密なネットワークが形成されている。Na – Na の RCN は x が増加するに伴い、0.4 – 0.5nm の領域では、RCNが 3 から 1 へ減少し、K – K では 0.5 から 2.8 へ増加した。Na+ イオンの方が K+ イオンに比べ、わずかに自身の周囲に同種のアルカリを配置する傾向が強いことが分かった。原子間相互作用の設
定の際、近接に働く斥力項の変数, ai はイオン半径に関するパラメータとして取り扱っており、aNa < aK
としてシミュレーションを行っている。このため、K+ イオンは Na+ イオンよりも近接に働く斥力が強
く、すなわちイオン半径の大きさがリング構造などの Si – O ネットワーク構造に影響を与えていると考えている。これらの結果より、アルカリ金属イオンを交換した場合の Si – O ネットワーク構造の変化に関して、混合アルカリ効果の傾向が見られないことから、アルカリ金属イオンの位置情報が混合アルカリ効果と関係していることが考えられる。
Si – O ネットワーク構造中の BO 及び NBO とアルカリ金属イオン間の配位構造を解析した。
このとき解析に使用するカットオフ距離は Na – O と K – O 間の PCF よりそれぞれ 0.30 nm 及び0.40 nm とした。解析した結果を Figure 3 – 10 及び 3 – 11 に示す。Figure 3 – 10 より K – NBO 間
距離の方が Na – NBO 間距離も長距離側へ現れるのは、K+ イオンのイオン半径、すなわちシミュ
レーションでは近接反発力が Na+ イオンよりも大きいためと考えられる。また、両アルカリ金属
イオンにおいて、– NBO 間距離の方が– BO 間距離に比べ短くなった。この結果は形式的に、Oイオンの電荷を – 2 とした場合、BO は – 2 / 2 = − 1 を、NBO では – 2 / 1 = − 2 と考えることができ、より Coulomb 相互作用が強く働く NBO – アルカリ金属イオン間距離が BO – アルカリ金属イオン間距離よりも短距離の距離分布を示すことが予測できる。
Figure 3 – 11 に O イオン周辺に配位する Na+ イオン及び K+ イオンを解析した。両アルカリ金属イオンにおいて、NBO の周辺にアルカリ金属イオンが存在しやすいことが分かった。これは、Figure 3 – 10 より、両アルカリ金属イオンが NBO との距離が近いため、配位数が多くなったと考えられる。NBO 周辺の両アルカリ金属イオンでは、異種アルカリが混合されることにより、配位数が減少しており、これは、対象としている系中のアルカリ金属イオンの存在量が減少している
ことと対応している。BO – K においても同様の理由で、アルカリの交換による配位数の傾向を反映していると考えている。BO – Na の配位数は BO – K と比べると、異種アルカリ金属イオンの混合によって配位数の変化がほぼ一定の値を示した。これは Na+ イオンの方が K+ イオンよりも近接反発力が弱く働くため、周囲のイオンとの相互作用によって遠くへジャンプしづらくなり、
同一の O イオンの近くに存在する確率が多くなるためではないかと考えられる。Figure 3 – 3 のNa – O と K – O の PCF において、x が増加するに伴い短距離側へシフトしたのは、Figure 3 – 10のように、アルカリ金属イオンが負に大きなポテンシャルを持つ NBO との結合を好むようになったためと考えられる。

66 / 145

0.40

Distance of alkali - oxygen, r / nm

0.38
0.36

K - BO

0.34
0.32
K - NBO

0.30
0.28

0.26

Na - BO

0.24
Na - NBO

0.22
0.20
0

0.2

0.6
0.4
Composition, x

0.8

1

Figure 3 – 10. Distance of alkali metal ions – BO and NBO in 0.25{xK2O-(1−x)Na2O}-0.75SiO2
glasses simulated at 900 K. Circles and triangles show distance of Na – oxygen and K – oxygen,
respectively. Solid and broken line show BO and NBO, respectively.

67 / 145

Coordination number of alkali around oxygen (counts)

3.50
3.00
NBO - K
2.50

NBO - Na

2.00
1.50

BO - K

1.00

BO - Na
0.50
0.00
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 3 – 11. Coordination number of BO and NBO – alkali metal ions in 0.25{xK2O(1−x)Na2O}-0.75SiO2 glasses simulated at 900 K. Circles and triangles show coordination
number of oxygen – Na and oxygen – K, respectively. Solid and broken line show BO
and NBO, respectively.

68 / 145

Figure 3 – 6 に示した Na+ イオンと K+ イオンの自己拡散係数は混合アルカリ効果を示した。
133

K+ イオンの自己拡散係数は 39K+ イオンよりも小さい値を示すのは、自身の質量が重くなったことで、拡散しづらくなっていることを示している。一方、Na – 133K 系の Na+ イオンは、Na – 39K
系の Na+ イオンに比べ、わずかに自己拡散係数が減少する傾向を示しており、対となる K+ イオ
ンの質量によって影響が現れた可能性がある。この傾向は、温度が増加するほど明確になった。
Figure 3 – 12 に、700 – 1200 K におけるアルカリ金属イオンの自己拡散係数のアレニウスプロ
ットを示す。また、Figure 3 – 13 に Figure 3 – 12 のアレニウスプロットの傾きから求めたアルカ
リ金属イオンの拡散の活性化エネルギーを示す。Table 3 – 3 にアルカリ金属イオンの自己拡散係数と拡散の活性化エネルギーをまとめた。Na-39K 系において、Na+ イオンの活性化エネルギーは43.8 ～ 56.3 kJ/mol であり、K+ イオンは 36.2 ～ 45.9 kJ/mol である。実験報告にある活性化エ
ネルギーの値 [5] は、Na+ イオンでは 115.5 ～ 65.7 kJ/mol であり、K+ イオンでは 64.4 ～ 106.7
kJ/mol であった。Figure 3 – 13 に示した Na-133K 系の Na+ イオンの活性化エネルギーは Na-39K
系の Na+ イオンの活性化エネルギーと近い値を示している。Na-133K 系において、39K から 133K
イオンへと拡散しづらくなることは、133K イオン自身が重くて動きづらく、133K イオンが Na+ イ
オンの拡散を妨げる、二つの役割を持っているためと考えられる。そのため、Na-39K 系と Na-133K
系のアルカリ金属イオンの活性化エネルギーが類似していると考えられる。
すべての x において、K+ イオンの方が Na+ イオンよりも活性化エネルギーが小さく、これは、K+ イオンの方が Na+ イ
オンよりも拡散しやすいことを示していると考えている。

69 / 145

-11(a) K+ ions of Na-39K system
0.4

0.5

0.6

0.8

0

1

-12

-12

-13

-13
ln(D / cm2 s-1)ln(D / cm2 s-1)

0.2

-14

-15

-16

-16

-17

0.4

0.5

0.6

0.8

-17
0.8

0.9

1.0

1.1
1.2
1.3
T-1 ×1000 / K-1

1.4

1.5

(c) K+ ions of Na-133K system
0.2

0.4

0.5

0.8

0.6

0.8

-12

-13

-13

-15

1.0

0

1

-12

-14

0.9

1.1 1.2 1.3
T -1 ×1000 / K -1

1.4

1.5

0.6

0.8

(d) Na+ ions of Na-133K system-11ln(D / cm2 s-1)ln(D / cm2 s-1)

0.2

-14

-15

-11

(b) Na+ ions of Na-39K system

-11

0.2

0.4

0.5

-14

-15

-16

-16

-17

0.8

0.9

1.0

1.1 1.2 1.3
T ×1000 / K -1
-1

1.4

1.5

-17
0.8

0.9

1.0

1.1 1.2 1.3
T -1 ×1000 / K -1

1.4

1.5

Figure 3 – 12. Arrhenius plot of Na+ and K+ ions for 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses.

70 / 145

135(a) Experimental result [5]60(b) M D simulation resut of this work

55
Activation energy, E / kJ mol-1

Activation energy, E / kJ mol-1

Na-133K: △K (Na-K)

115
105

Na (Na-K)

95
85
75
65
55
45

35
0.00

Na (Na-133K)

Na-39K: ○

125

0.20

0.40
0.60
0.80
Composition, x

1.00

Na (Na-39K)

50

K (Na-39K)

45

40
K (Na-133K)

35
0.00

0.20

0.40
0.60
0.80
Composition, x

1.00

Figure 3 – 13. Activation energy of Na+ and K+ ions in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses.

71 / 145

Table 3 – 3. Self-diffusion coefficients log10(D / cm2 s−1) and activation energies E for 23Na, 39K, and 133K
diffusion in 0.25{xK2O-(1−x)Na2O}-0.75SiO2 glasses.

x

523

623

673

0.00
0.20
0.40
0.50
0.60
0.80
1.00

Temperature, T / K
700
723
800
900
39
Na- K systems
23
Na
-6.27
-5.86 -5.54-6.49-6.06 -5.89-6.64-6.21 -5.75-6.83-6.39 -5.85-6.82-6.55 -6.02-7.12-6.46 -6.04

E / kJ mol-1

1000

1100

1200

-5.25
-5.38
-5.58
-5.66
-5.63
-5.78

-5.10
-5.24
-5.31
-5.40
-5.49
-5.58

-4.90
-5.01
-5.12
-5.19
-5.28
-5.31

43.8
47.9
49.0
53.1
52.5
56.3

-5.89
-5.79
-5.64
-5.74
-5.57
-5.48

-5.60
-5.44
-5.53
-5.46
-5.31
-5.29

-5.35
-5.37
-5.25
-5.24
-5.18
-5.06

-5.19
-5.05
-5.08
-5.08
-4.98
-4.93

45.9
46.0
45.2
43.0
41.2
36.2

Na-133K systems
23
Na-5.86 -5.54-6.13 -5.73-6.31 -5.93-6.40 -6.06-6.47 -6.13-6.63 -6.21

-5.25
-5.49
-5.53
-5.71
-5.77
-5.94

-5.10
-5.22
-5.32
-5.51
-5.53
-5.69

-4.90
-5.06
-5.20
-5.25
-5.35
-5.49

40.1
46.7
43.0
41.9
39.9
37.5

-5.80
-5.63
-5.67
-5.62
-5.59
-5.43

-5.57
-5.39
-5.35
-5.41
-5.37
-5.23

-5.42
-5.20
-5.26
-5.22
-5.22
-5.07

43.8
46.9
52.9
50.5
56.0
49.5

39

K

0.00
0.20
0.40
0.50
0.60
0.80
1.00

-6.55
-6.52
-6.47
-6.36
-6.27
-6.07

0.00
0.20
0.40
0.50
0.60
0.80
1.00

-6.27
-6.49
-6.79
-6.82
-7.10
-7.00

-6.30
-6.12
-6.13
-6.12
-5.89
-5.70

133

K

0.00
0.20
0.40
0.50
0.60
0.80
1.00

-6.65
-6.67
-6.56
-6.55
-6.46
-6.25

-6.34
-6.18
-6.21
-6.12
-6.10
-5.87

-5.98
-5.90
-5.93
-5.85
-5.77
-5.62

references [5]
22
Na
0.00
0.40
0.70
1.00

-9.57
-11.08-11.56-11.64

-8.53
-9.44
-9.81
-9.92

-8.12
-8.81
-9.13
-9.26

-7.77
-8.26
-8.54
-8.68

115.5
127.2
101.3
65.69
42

K

0.00
0.40
0.70
1.00

-12.39 -10.53 -9.80-12.39 -10.34 -9.55-11.49 -9.86 -9.22-10.01 -8.97 -8.55

-9.18
-8.86
-8.68
-8.20

72 / 145

64.4
101.7
108.4
106.7

アルカリ金属イオンの back tracking とは異種のアルカリによって、拡散しようとするサイトへのジャンプが邪魔され拡散が起きにくくなる運動を指す [23]。Figure 3 – 14 にアルカリ金属イオ
ンの back tracking の概念図を示す。過去の MD シミュレーションによる報告 [23, 25] において、混合アルカリ系は互いのアルカリ金属イオンによって、拡散を阻害していることを示した。これ
らの研究報告は、アルカリ金属イオンの back tracking が増加し、アルカリ金属イオンの拡散能を大きく減少させることを指摘した [23, 25]。 しかしながら、このときの解析は一部のピックアッ
プしたアルカリ粒子の解析にとどまっており、系全体の粒子平均のデータとして得る必要があると考えている。本研究では、カットオフ距離, rc を超えないアルカリ金属イオンの運動を backtracking と定義し、rc を超えるジャンプ運動を continuous jump として定義した。系中のそれぞれのアルカリ金属イオンに対して、back tracking と continuous jump は 600 ps 間でカウントした。
back tracking と continuous jump は以下の式を用いた。

𝐴∑𝑁 (𝑏 𝐴 )( 𝑖=1𝑁𝐴 𝑖 )⁄𝐴
𝐵ave. =
600 𝑝𝑠(3 – 7)𝐴∑𝑁 (𝐶 𝐴 )( 𝑖=1 𝐴 𝑖 )𝑁⁄𝐴𝐶ave. =
600 𝑝𝑠(3 – 8)BaveA と CaveA は系中のアルカリ金属イオン種 A の back tracking と continuous jump の平均のカウント数であり、biA と ciA はそれぞれアルカリ金属イオンの back tracking と continuous jump のカ
ウント数を示し、NA はアルカリ金属イオン種 A の粒子数を表す。BaveA と CaveA を解析時間 600ps
当たりの単位として規格化している。

73 / 145

Alkali oxide ion

Back tracking

rc : 0.5 nm

rc : 0.3 nm
rc / 2

SiO4 unit

Figure 3 – 14. Schematic diagram of back tracking motion.

74 / 145

Figure 3 – 15 と Figure 3 – 16 に 900 K における Na+、 39K+、及び 133K+ イオンの BaveA と CaveAを示す。このとき、rc = 0.3 nm、 0.5 nm とした。rc = 0.3 nm は Na – O と K – O の PCF の分布領
域 (Figure 3 – 3) であるため、アルカリが近接の SiO4 四面体から拡散を阻害されるかを調べるために設定した。rc = 0.5 nm は Na – Na, K – K 及び Na – K 間の PCF の分布 (Figure 3 – 3) を目安としており、アルカリ金属イオンと SiO4 四面体によるジャンプの阻害への影響について調査するために設定した。Figure 3 – 15 (a) の rc =0.3 nm の場合、Na – 39K と Na – 133K における Bave.Na+は x が増加するに伴い、徐々に増加した。Na – 39K 系の Bave39K+ と Na – 133K 系の Bave133K+ は xが減少するに伴い、増加した。これは Figure 3 – 15 (b) の rc = 0.5 nm の場合も同様の結果となった。この結果より、異種アルカリが混合されると遠くへジャンプする頻度が減少することが分かった。本シミュレーションは Bunde と Ingram ら [23] の経験的に設定した原子間相互作用を適用した MD 計算の報告と同様に、back tracking motion の存在を示した。また、Figure 3 – 7 – 2 で示した負の混合エンタルピーでは、二種類のアルカリ金属イオンを混合することで、相互作用が増大すると実験系研究において解釈されており [7, 44]、このとき本研究では back tracking motionが増加した。これは、back tracking motion による異種アルカリ金属イオン間の相互作用の増大が、混合アルカリ効果を引き起こす因子の一つであることを示していると考えている。
Figure 3 – 16 (a) rc = 0.3 nm の場合、Na – 39K 及び Na – 133K における Cave.Na+ は x が増加するに伴い、徐々に増加した。Na – 39K 系における Cave39K+ と Na – 133K 系における Cave133K+ は x が減少するに伴い増加した。 これは Figure 3 – 16 (b) rc = 0.5 nm のときも同様の結果となった。カットオフ距離が異なっても、back tracking の傾向が同じであり、自己拡散係数の傾向 (Figure 3 – 6)とも類似していた。Na+ イオンの back tracking はいずれの x においても 39K+ と 133K+ イオンよりも頻度が増えた。 39K+ と 133K+ イオンの continuous jump はいずれの x においても Na+ イオンに比べ頻度が増加した。これは Figure 3 – 16 の K+ イオン の方が Na+ イオンよりも自己拡散係数が大きい傾向に対応しており、K+ イオン の方が Na+ イオンよりも遠くへジャンプする頻度が多いことを示す。また、Figure 3 – 11 に示した、BO – アルカリ金属イオンの配位数の傾向にも対応しており、K+ イオン の方が Na+ イオンよりも近接に働く斥力が大きいため、遠くへジャンプする頻度が多いことを示していると考えている。Na – 39K 系の CaveNa+ は Na – 133K 系よりも、わずかに増加した。rc = 0.3 nm が rc = 0.5 nm よりも continuous jump の頻度が多いのは、アルカリ金属イオンの近接に存在する、SiO4 四面体ネットワーク構造による拡散の阻害が小さいことを示し、SiO4 四面体ネットワーク構造を飛び越えて、遠くへジャンプする確率が多くなると考えられる。

75 / 145

(a) r c = 0.3 nm
29000

Na in 133K
27000
K in 133K

Bave.i

25000
23000

21000
Na in 39K
19000

K in 39K

17000

15000
0

0.2

0.4
0.6
Composition, x

0.81(b) rc = 0.5 nm
31000

Na in 133K

30000

K in 133K

Bave.i

29000
28000
27000

Na in 39K
K in 39K

26000
25000
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 3 – 15. Back tracking motion of Na-39K (circles) and Na-133K (triangles) systems in 0.25{xK2O(1−x)Na2O}-0.75SiO2 glasses by rc = 0.3 nm (a) and rc = 0.5 nm (b) at 900 K (open and solid symbols
show K+ and Na+ ions, respectively).

76 / 145

(a) r c = 0.3 nm
12000

K in 39K
10000

Na in 39K
Cave.i

8000
6000
K in 133K
4000
Na in 133K
2000
0

0

0.2

0.4
0.6
Composition, x

0.81(b) r c = 0.5 nm

4500
4000

K in 39K

3500

Cave.i3000

Na in 39K

2500
2000
1500

K in 133K

1000

500

Na in 133K

0
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 3 – 16. Continuous jump motion of Na-39K (circles) and Na-133K (triangles) systems in 0.25{xK2O(1−x)Na2O}-0.75SiO2 glasses by rc = 0.3 nm (a) and rc = 0.5 nm (b) at 900 K (open and solid symbols show
K+ and Na+ ions, respectively).

77 / 145

アルカリ金属イオンの質量による影響を調査するために、39K133KSi3O7 系の K+ イオンの backtracking についても調べた。 Figure 3 – 17(a), (b) に rc = 0.5 nm の 700 K から 1200 K における、
39

K133KSi3O7、 39K2Si3O7、 及び 133K2Si3O7 ガラスの BaveA と CaveA を示す。133K2Si3O7 系の Bave.133K+

はすべての系と比較し、最も大きくなった。39K2Si3O7 系の Bave39K+ は 133K2Si3O7 系の Bave.133K+ よりも小さくなった。133K+ 自身は重いことから、133K2Si3O7 系の Cave.133K+ はすべての系で最も小さい値を示した。39K2Si3O7 系の Cave39K+ は 133K2Si3O7 系の Cave.133K+ よりも頻度が多くなった。温度が増加するに伴い遠くへジャンプする頻度が増えると考えられるため、back tracking が減少し
continuous jump が増加したと考えられる。39K133KSi3O7 系において、もし、アルカリ金属イオンのサイズのみが拡散に影響を与えるならば、 39K133KSi3O7 系の 39K と 133K の back tracking とcontinuous jump は、39K2Si3O7 系の 39K と 133K2Si3O7 系の 133K がそれぞれ同じ頻度を示すべきであるが、本研究の結果は互いの系で異なる結果となった。このことから、混合アルカリ効果は、
アルカリ金属イオンのサイズ差だけではなく、質量差も関連している現象であることを本研究で示した。

78 / 145

(a) Backtracking jump
31000
29000

BAve.i per 600 ps

27000
25000

23000
21000
133K in 133K Si O
2 3 7
133K in 39K 133KSi O
3 7
39K in 39K 133KSi O
3 7

19000

17000

39

K in 39K 2Si3O 7

15000

600

900

1200
1500
Temperature, T / K

1800(b) Continuous jump

14000
39K in 39K Si O
2 3 7
39K in 39K 133KSi O
3 7

12000

133K in 39K 133KSi O
3 7

CAve.i per 600 ps

10000

133

K in 133K 2Si3O 7

8000
6000
4000

2000
0

600

900

1200
1500
Temperature, T / K

1800

Figure 3 – 17. Back tracking, (a) and continuous jump, (b) motion of K+ ions in 39K, 133K and
39K – 133K systems at 900 K. Solid cicles, open cicles, solid triangles and open triangles show
39K+ ions in 39K, 39K+ ions in 39K – 133K, 133K+ ions in 39K – 133K and 133K+ ions in 133K systems,

respectively.

79 / 145§3 – 6. 結言
0.25{xK2O-(1−x)Na2O}-0.75SiO2 ガラスに適用する原子間相互作用を第 2 章の方法により設定した。
この原子間相互作用を適用することにより、MD シミュレーションによるモル体積、原子間距離、Qn
比、アルカリ金属イオンの自己拡散係数及び活性化エネルギーは実験報告の傾向を再現した。
モル体積の結果より、異種アルカリが混合されると系中のパッキング性が増したことから、アルカ
リ金属イオンのサイズ効果が混合アルカリ効果に影響していることを示した。x が増加すると 3 ~ 5 員環構造が増加し、6 員環以上の構造は減少したことから、Na rich 組成の方が密なネットワーク構造を作っていることが分かった。
アルカリ金属イオンのダイナミクスを調べたところ、異種アルカリが混合されると遠くへジャンプする頻度が減少する back tracking motion が起きた。異種アルカリ金属イオンによる backtrakcking motion が起きることで、混合アルカリ効果が起きる可能性を示した。また、SiO4 四面体
ネットワーク構造が混合アルカリ効果に与える影響は小さいことが示唆された。39K133KSi3O7 系の 39K と 133K の back tracking と continuous jump は、39K2Si3O7 系の 39K と 133K2Si3O7 系の 133K
がそれぞれ同じ頻度を示さなかった。そのため、混合アルカリ効果は、アルカリ金属イオンのサイズ差だけではなく、質量差も関連している現象ということを明らかにした。

80 / 145§3 – 7. 参考文献
[1] K. L. Nagi, The dynamics of ions in glasses: importance of ion-ion interactions, J. Non-Cryst. Solids
323 (2003) 120 – 126.
[2] A. Bunde, K. Funke, M. D. Ingram, Ionic glasses: History and challenges, Solid State Ionics 105 (1998)
1 – 13.
[3] J. O. Isard, The mixed alkali effect in glass, J. Non-Cryst. Solids 1 (1969) 235 – 261.
[4] D. E. Day, Mixed alkali glasses – their properties and uses, J. Non-Cryst. Solids 21 (1976) 343 – 372.
[5] J. W. Fleming Jr. and D. E. Day, Relation of alkali mobility and mechanical relaxation in mixed-alkali
silicate glasses, J. Am. Ceram. Soc. 55 (1972) 186 – 192.
[6]

B. Vessal, G. N. Greaves, P. T. Marten, A. V. Chadwick, R. Mole, S. Houde-Walter, Cation
microsegregation and ionic mobility in mixed alkali glasses, Nature 356 (1992) 504 – 506.

[7] R. Terai, The mixed alkali effect in the Na2O-Cs2O-SiO2 glasses, J. Non-Cryst. Solids 6 (1971) 121 –135.
[8] A. Riviére and L. C.-Arnault, Mechanical relaxations in single and mixed alkali silicate glasses, Mate.
Sci. Engineer. A 521 – 522 (2009) 272 – 275.
[9] C. T. Moynihan, N. S. Saad, D. C. Train, A. V. Lesikar, Mixed-alkali effect in the dilute foreign-alkali
region. Failure of the strong electrolyte/cationic interaction modeⅠ, J. Am. Ceram. Soc. 63 (1980) 458– 464.
[10] J. E. Shelby Jr. and D. E. Day, Mechanical relaxations in mixed alkali silicate glasses: Ⅱ, discussion,
J. Am. Ceram. Soc. 53 (1970) 182 – 187.
[11] G. L. Mcvay and D. E. Day, Diffusion and internal friction in Na-Rb silicate glasses, J. Am. Ceram.
Soc. 53 (1970) 508 – 513..
[12] S. N. H.-Walter, J. M. Inman, A. J. Dent, G. N. Greaves, Sodium and silver environments and ionexchange processes in silicate and aluminosilicate glasses, J. Phys. Chem. 97 (1993) 9330 – 9336.
[13] G. N. Greaves and K. L. Nagi, Reconciling ionic-transport properties with atomic structure in oxide
glasses, Phys. Rev. B 52 (1995) 6358 – 6380.
[14] T. Minami, Y. Tokuda, H. Masai, Y. Ueda, Y. Ono, S. Fujimura, T. Yoko, Structural analysis of alkali
cations in mixed alkali silicate glasses by 23Na and 133Cs MAS NMR, J. As. Ceram. Soc. 2 (2014)
333 – 338.
[15] S. Kaneko, Y. Tokuda, Y. Takahashi, H. Masai, Y. Ueda, J. As. Ceram. Soc, 5 (2017) 7 – 12.
[16] E. I. Kamitsos, A. P. Patsis, G. D. Chryssikos, Infrared reflectance investigation of alkali diborate
glasses, J. Non-Cryst. Solids 152 (1993) 246 – 257.
[17] E. I. Kamitsos, Y. D. Yiannopoulos, H. Jain, W. C. Huang, Far-infrared spectra of alkali germanate
glasses and correlation with electrical conductivity, Phys. Rev. B 54 (1996) 9775 – 9783.
[18]

J. Swenson, A Matic, C. Karlsson, L. Börjesson, C. Meneghini, W. S. Howells, Random ion
distribution model: A structural approach to the mixed-alkali effect in glasses, Phys. Rev. B 63 (2001)
132202.

[19] C. Karlsson, E. Zanghellini, J. Swenson, B. Roling, D. T. Boeron L. Borjesson, Structure of mixed
81 / 145

alkali/alkali-earth silicate glasses from neutron diffraction and vibrational spectroscopy, Phys. Rev. B
72 (2005) 064206.
[20] S. R. Elliott, Frequency-dependent conductivity in ionic glasses: A possible model, Solid State Ionics
27 (1988) 131 – 149.
[21] M. D. Ingram, Towards a theory of ion transport in glass, Physica A 266 (1999) 390 – 399.
[22] P. Maass. A. Bunde, M. D. Ingram, Ion transport anomalies in glasses, Phys. Rev. Lett. 68 (1992) 3064
– 3067 .
[23] A. Bunde, M. D. Ingram, P. Maass, K. L. Nagi, Diffusion with memory: a model for mixed alkali
effects in vitreous ionic conductors, J. Phys. A: Math. Gen. 24 (1991) L881.
[24] S. Balasubramanian and K. J. Rao, A molecular dynamics study of the mixed alkali effect in silicate
glasses, J. Non-Cryst. Solids 181 (1995) 157 – 174.
[25] J. Habasaki, I. Okada, Y. Hiwatari, MD study of mixed alkali effect in a lithium-potassium metasilicate
glass, J. Non-Cryst. Solids 183 (1995) 12 – 21.
[26] A.N. Cormack, J. Du, T.R. Zeitler, Sodium ion migration mechanisms in silicate glasses probed bymolecular dynamics simulations, J. Non-Cryst. Solids 323 (2003) 147 – 154.
[27] N. Sawaguchi, K. Yamaguchi, M. Sasaki, K. Kawamura, Interatomic potential model for molecular
dynamics simulation of lithium borate melts/glasses, J. Comput. Chem. Jpn. 14 (2015) 139 – 146.
[28] K. Kawamura and N. Sawaguchi, Molecular dynamics calculations of mixing of silicate melts, Bull.
Volcanol. Soc. Jpn. 42 (1997) S1 – S11.
[29] J. Habasaki and I. Okada, Molecular dynamics simulation of alkali silicate based on the quantum
mechanical potential surfaces, Mol. Sim. 9 (1992) 319 – 326.
[30] A.K. Pant and D.W.J. Cruickshank, The crystal structure of alpha-Na2Si2O5, Acta Cryst. B24 (1968)
13 – 19.
[31] S. Rakic and V. Kahlenberg, B. C. Schmidt, Hydrothermal synthesis and structural characterization of
kappa-(Na2SiO5) and Na1.84K0.16Si2O5, Solid State Sciences 5 (2003) 473 – 480.
[32] S. Rakic and V. Kahlenberg, The crystal structure of a mixed alkali phyllosilicate with composition
Na1.55K0.45Si2O5, European Journal of Mineralogy 13 (2001) 1215 – 1221.
[33] S. Rakic, V. Kahlenberg, Single crystal structure investigation of twinned NaKSi2O5-a novel single
layer silicate 3 (2001) 659 – 667.
[34] A. L. Spek and B. H. de Jong, Acta cryst. Section E. Structure Reports Online 61 (2005) i188 – i190.
[35] B.H. de Jong, H.T. Supèr, A.L. Spek, N. Veldman, G. Nachtegaal, J.C. Fischer, Mixed alkali systems:
Structure and 29Si MASNMR of Li2Si2O5 and K2Si2O5, Acta Cryst. B54 (1998) 568 – 577.
[36] K. Kawamura, MXDORTO, Japan Chemistry Program Exchange P029, (1996).
[37] J. P. Hansen, I. R. Mcdonald, Theory of Simple Liquids, Academic Press, (1986).
[38] J.E. Shelby, Jr. and D. E. Day, Mechanical relaxations in mixed-alkali silicate glasses: 1, results, J.
Amer. Ceram. Soc. 52 (1969) 169 – 174.
[39] R.L. Mozzi and B.E. Warren, The structure of vitreous silica, J. Appl. Cryst. 2 (1969) 164 – 172.
[40] G.N. Greaves, A. Fontaine, P. Lagarde, D. Raoux, S.J. Gurman, Local structure, of silicate glasses,
82 / 145

Nature 293 (1981) 611 – 616.
[41] O. Gedeon, K. Jurek, I. Drbohlav, Mixed-alkali effect in sodium potassium glasses irradiated withelectrons, J. Non-Cryst Solids. 356 (2010) 456 – 560.
[42] P. Florian, K. E. Vermillion, P. J. Grandietti, I. Farnan, J. F. Stebbins, Cation distribution in mixedalkali disilicate glasses, J. Am. Chem. Soc. 118 (1996) 3493 – 3497.
[43] G. N. Greaves, A. Fontaine, P. Lagarde, D. Raoux, S. J. Gurman, Structural studies of the mixed alkali
effect in disilicate glasses, Solid State Ionics 105 (1998). 243 – 248.
[44] K. Takahashi and T. Yoshio, Heats of mixing in mixed alkali silicate glasses, Zairyou 26 (1976) 785 –
790.
[45] M. Tomozawa, The mixed alkali effect and thermodynamic state of glasses, Solid State Ionics 105(1998) 249 – 255.
[46] J. Lumsden, Thermodynamics of molten mixtures of alkali-metal halides, J. Discussions Faraday Soc.
32 (1961) 138 – 146.
[47] H. Maekawa, T. Maekawa, K. Kawamura, T. Yokokawa, The structural groups of alkali silicate glasses
determined from 29Si MAS-NMR, J. Non-Cryst. Solids 127 (1991) 53 – 64.

83 / 145

第 4 章 MD シミュレーションによる xSiO2 - (1 - x)NaAlO2 ガラスの構造評価§4 – 1. 緒言
アルミノケイ酸塩系ガラスは機械的特性や化学的耐久性に優れていることから、スマートフォンなどのカバーガラスに使われている [1, 2]。これらの物性は原子スケールのガラス中の構造によって決定されると考えられている。
一般的に、Na2O – SiO2 – Al2O3 系ガラスにおいて、Al / Si 及び Al / Na 比 = 1 のとき、Si イオ
ンは O イオンと SiO4 ユニットを作る。Al イオンも O イオンと AlO4 ユニットを作る。このとき、[AlO4]1– の近接に Na+ イオンが存在することで電気的中性を保つと考えられている (Figure 4– 1)。Al / Si 及び Al / Na 比が 1 ではないとき、O イオンの結合について、いくつかの異なる説
や実験報告が存在する。その一つとして、Lacy [3] は、O イオンが 3 つの AlO4 ユニットと結合する酸素トリクラスター構造を形成する説を唱えた (Figure 4 – 1 )。Schmucker ら [4] の NMR による研究は、Al2O3 – SiO2 と Na2O – Al2O3 – SiO2 系ガラスにおいて酸素トリクラスター構造が存在することを報告した。一方、ラマンスペクトル解析では、Al/Na 比が 0.03 ~ 1.61 において、ア
ルミニウムイオンの量比が増加すると、酸素トリクラスター構造は存在せず AlO4 ユニットがAlO6 ユニットの形で存在することを報告している [5]。また、ラマンスペクトル解析では酸素ト
リクラスター構造のラマンスペクトルバンド (418 cm – 1 及び 1200 cm –1) が Na2O – Al2O3 – SiO2系ガラスにおいて、確認されないことが報告されており、たとえすべて AlO4 ユニットとなったとしても酸素トリクラスター構造の存在を否定する指摘もある [5]。系中の O イオンが多くなる
時、Al イオンが 6 配位構造を取る報告も存在する。Okuno ら [6] は X 線散乱、 IR 及びラマン
分光測定により、SiO2 – Al2O3 系ガラスにおいて、Al2O3/SiO2 比が増加するに伴い、 AlO5 もしく
は AlO6 ユニットが存在することを報告した。Morikawa ら [7] は X 線散乱測定により、SiO2 –Al2O3 系ガラスにおいて AlO5 ユニットの存在を報告した。ガラス内の bond order についても議論されており、Lowenstein [8] は Al – O – Al 構造は熱力学的理由から、Si – O – Al と Si – O – Si 構造に比べて存在確率が低いという説を唱えているが、Stebbins と Xu [9] が 3Q MAS NMR によって、NaAlSi3O8 と CaAl2Si2O8 ガラスには Al – O – Al 構造が存在することを報告したことで、bondorder も配位数と同様に複雑な構造を取っていることが予想されている。このように、アルミノケ
イ酸塩系ガラスは実験的観測によっても、構造に様々な解釈が存在する。

84 / 145

Network modifier
Non-bridging oxygen (NBO)

Bridging oxygen (BO)
Q2

SiO4

Q3

x = 0.5
AlO4

Q0

BO

±0

-1

+1

Silicon, Si

Q46-membered ringNa

Network modifier
Bridging oxygen (BO)Q3

Q2

4-membered ringsx ≠ 0.5

Q1

Non-bridging oxygen (NBO)

5-membered rings
aluminum, Al

Q0

Silicon, Si

Q46-membered ringQ1

4-membered rings

tricluster

SiO45-membered rings
AlO4
aluminum, Al

Figure 4 – 1. Charge neutralization of AlO4 units in xSiO2 - (1 – x )NaAlO2 glasses.

85 / 145原子サイズのオーダーで構造解析可能な手段の一つとして、MD シミュレーションを用いた研究も行われてきた [10 – 17]。Montorsi ら [10] は Na2O – Al2O3 – SiO2 ガラスの Al2O3/SiO2 比を変化させたときの MD シミュレーションを行い、酸素トリクラスター構造と AlO4 ユニットが存在することを報告した。Xiang ら [17] も Na2O – Al2O3 – SiO2 系ガラスにおいて、酸素トリクラス
ター構造と AlO4 ユニットが存在することを報告した。これらの研究では、粒子, N、体積, V 及び温度, T を一定とする、NVT アンサンブル下の MD シミュレーションが行われてきたが、NVT ア
ンサンブル下で適用した原子間相互作用を用いて、圧力, P を一定とする NPT アンサンブル下のシミュレーションを行ったところ、歪んだ構造とされている AlO3 ユニットが現れ、実験結果に比べ密度が極端に大きくなることが報告されている [17]。NPT アンサンブル下のシミュレーションも可能とすることは、融体からガラスへ冷却される際の構造についても調査可能となるため、実際の材料開発ツールとして役立つと考えている。
これまで、本研究では NPT アンサンブル下でのシミュレーションを行ってきており、Na2O –SiO2 系ガラス (第 2 章) [18] や Na2O – K2O – SiO2 系ガラス (第 3 章) [19] の構造と物性の再現性が良いことを報告した。また、第 2, 3 章ではネットワーク形成酸化物を SiO2 のみであるとし原子間相互作用の設定を試みてきたが、アルミノケイ酸塩系では中間酸化物である Al2O3 においてもネットワーク構造を形成するため、Si – O と Al – O を同時に考慮する必要があると考えている。
本章では、xSiO2-(1−x)NaAlO2 ガラス(0 ≤ x ≤ 1) に注目した。この組成は、酸化物イオンの数が系中では変化しないと考えられるため、Al イオンの配位数は四配位構造のみの多面体構造を取ることが考えられる。そのため、トリクラスター構造を形成するかどうか、また Al イオンと Si
イオンの bond order に焦点を当てた調査ができると考えた。NPT アンサンブル下での MD シミュレーションを xSiO2-(1−x)NaAlO2 系ガラスで行い、これまでに報告してきた原子間相互作用が二種類のペアポテンシャルを必要とする場合についても対応可能であるかを検討し、アルミノケ
イ酸塩ガラスの構造について調査した。
§4 – 2 は原子間相互作用について記述した。§4 – 3 は MD シミュレーションの計算条件を示した。§4 – 4 は MD シミュレーションで得た計算結果を示す。§4 – 5 は§4 – 4 で得られた結果を基に、配位構造と bond order について考察した。§4 – 6 及び§4 – 7 はそれぞれ、本章の結言と参考文献を記した。

86 / 145§4 – 2. 原子間相互作用
原子間相互作用は第 2 章の方法で決定した。原子間相互作用は式 (2 – 1) を使用し、イオンの電
荷は第 2 章§2 – 3 の計算条件を用いて yA2O-(1−y)Al2O3、及び xSiO2-(1−x)AAlO2 系 [20 – 38] (A =Li, Na 及び K) 結晶に対し電荷解析を行った。LiAlSiO4 と KAlSiO4 結晶に関しては、化学量論比を考慮してそれぞれ 2 × 1 × 1 のスーパーセルを作成し、計算した。使用したノルム保存型擬ポテ
ンシャルは Li, 2s1; Na, 2s22p63s1; K, 3s23p64s1; Al, 3s23p1; Si, 3s23p2; 及び O, 2s22p4 である。Table 4 –1 に結晶構造データ [20 – 38] と計算結果を示す。Figure 4 – 2 に yA2O-(1−y)Al2O3、及び xSiO2(1−x)AAlO2 系結晶の電荷解析結果を示す。Figure 4 – 2 (a) より、Figure 2 – 3 の結果と同様に、
Al2O3 に対しアルカリ金属酸化物を加えると、アルミニウムイオンとアルカリ金属酸化物イオンの電荷が減少した。Figure 4 – 2 (b) より、y が増加するほど、すべての陽イオンの電荷が増加する
傾向を示した。この電荷解析の傾向を基に、xSiO2-(1−x){yNa2O-(1−y)Al2O3} 系中のイオンの電荷を、以下の式 (4 – 1) - (4 – 4) で示した。
𝑧𝑁𝑎 = +0.1𝑥 + (−0.25𝑦 + 1.0) (0 ≤ 𝑥 < 1) (0 < 𝑦 ≤ 1)

(4 – 1)𝑧𝐴𝑙 = +0.15𝑥 + (−0.30𝑦 + 2.1) (0 ≤ 𝑥 < 1) (0 ≤ 𝑦 < 1)

(4 – 2)𝑧𝑆𝑖 = +0.25𝑥 + (−0.50𝑦 + 2.4) (0 < 𝑥 ≤ 1) (0 ≤ 𝑦 ≤ 1)

(4 – 3)𝑧𝑂 = neutrality(4 – 4)zO に系の電気的中性を保たせた。SiO2+ 及び AlO+ クラスターモデルを作成し、第 2 章§2 – 4
で示した計算条件により MO 計算を行った。Figure 4 – 3 に x = 0.5 における Si – O 及び Al – Oのカーブフィッティングの結果を示す。このとき、Si – O 及び Al – O のカーブフィッティングを同時に行い、原子間相互作用を決定した。aNa, bNa, cNa は経験的に決めた値を使用した。Table 4 – 2に各組成におけるポテンシャルパラメータを示す。

87 / 145

88 / 145

Table 4 – 1. The deviation denotes the difference of structural optimization by DFT. Deviation shows in absolute value.

(a) yA2O-(1−y)Al2O3

3

M ulliken charge

2
Al
1
A

0
-1

O

-2
0

3

0.2

0.4
0.6
Composition, y

0.81(b) xSiO2-(1−x)AAlO2
Si

M ulliken charge

2

Al

A

1
0
-1

O

-2

0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 4 – 2. Mulliken charges of ions in yA2O-(1−y)Al2O3, (a) and xSiO2-(1−x)AAlO2 , (b) crystals(A = Li, Na, and K) (listed in Table 4 – 1). Circles, squares, diamonds and triangles are Al, A, Si and O
ions, respectively. White, black and dotted are Li, Na and K, respectively.

89 / 145

Energy difference, ΔE×1018 / J

0.030(a) SiO2+

SiO2+
0.025

Si

O

0.020

MO
MD

0.015
0.010
0.005
0.000
0.140

0.145
0.150
Distance, rSiO / nm

0.155(b) AlO +Energy difference, ΔE×1018 / J

0.016

0.014

Al

AlO+
O

0.012
0.010

MO
MD

0.008
0.006
0.004

0.002
0.000
0.150

0.155
0.160
Distance, rAlO / nm

0.165

Figure 4 – 3. Energy surface determined from MO calculation (MP2/6-311+g(d)) of SiO2+
model, (a) and AlO+ model, (b). Potential energies calculated by eq. (2 – 1) with fitted
parameters for x = 0.5.

90 / 145

91 / 145

Table 4 – 2. Determined interatomic potential parameters for the simulation of xSiO2-(1−x)NaAlO2 glasses.

§4 – 3. MD シミュレーションの計算条件
MD シミュレーションは MXDORTO システム [39] を使用し、時間刻みは 2.0 fs とした。三次元周期境界下で計算を行った。NPT アンサンブルのシミュレーションを行い、一辺が約 5.0 nm の
立方体 MD セル中に疑似乱数を用いて粒子(N = 4950 ～ 5280) をランダムに配置させ、圧力 (P =0.1 MPa ) および温度 (T = 3000 K) の融体のシミュレーションを行い、250 000 step にかけて全エ
ネルギーが一定になるまで構造緩和を行った。Table 4 – 3 に本章のシミュレーションで扱った粒子数を示す。ガラス状態は冷却速度を 0.01 K/step とし、段階的に冷却することで得た。3000–2100
K では 250 000 steps 、ガラス状態に近づくほど、イオンの拡散が遅くなることを考慮し、1800–
300 K では 500 000 steps の構造緩和を行った。300 K における構造をガラス状態として、10 000
steps の解析計算を行った。

Table 4 – 3. Number of particles (N) in xSiO2-(1−x)NaAlO2 glass simulations.

x

O

Si

Al

Na

Total (N )

0

2600

0

1300

1300

5200

0.33

2880

480

960

960

5280

0.41

2924

602

860

860

5246

0.47

2964

702

780

780

5226

0.50

3000

750

750

750

5250

0.60

3100

930

620

620

5270

0.67

3120

1040

520

520

5200

0.75

3200

1200

400

400

5200

0.80

3300

1320

330

330

5280

0.83

3312

1380

276

276

5244

0.86

3360

1440

240

240

5280

1

3300

1650

0

0

4950

92 / 145§4 – 4. 結果
Table 4 – 4 に Al – O と Si – O の第一近接配位数を示す。すべての組成において AlO4 及び SiO4ユニットであった。

Table 4 – 4. The nearest coordination number of Al–O and Si–O in xSiO2-(1−x)NaAlO2 glasses.
Al-O (%)

Si-O (%)

x

3

4

5

6

3

4

5

6

0.00

0.00

99.92

0.08

0.00

ー

ー

ー

ー

0.33

0.10

99.90

0.00

0.00

0.00

100.00

0.00

0.00

0.41

0.47

99.53

0.00

0.00

0.00

100.00

0.00

0.00

0.47

0.38

99.62

0.00

0.00

0.00

100.00

0.00

0.00

0.50

0.27

99.73

0.00

0.00

0.00

100.00

0.00

0.00

0.60

1.61

98.39

0.00

0.00

0.00

100.00

0.00

0.00

0.67

1.35

98.65

0.00

0.00

0.00

100.00

0.00

0.00

0.75

2.50

97.50

0.00

0.00

0.08

99.92

0.00

0.00

0.80

2.42

97.58

0.00

0.00

0.23

99.77

0.00

0.00

0.83

3.26

96.74

0.00

0.00

0.07

99.93

0.00

0.00

0.86

2.92

96.67

0.42

0.00

0.00

100.00

0.00

0.00

1.00

ー

ー

ー

ー

0.42

99.58

0.00

0.00

93 / 145

Figure 4 – 4 に PCF を示す。第一近接 Si – O 及び Al – O の PCF は x が増加するに伴い、短距離側へシフトし、実験結果の傾向と一致した [40]。 x における Si – O 及び Al – O のトップピー
クはそれぞれ、0.164 nm – 0.167 nm 及び 0.174 nm – 0.177 nm であり、実験報告による Si – O 及び
Al – O 間距離である 0.160 – 0.161 nm [41] 及び 0.177 nm [42] と一致している。Na – O の PCF は
約 0.20 nm – 0.30 nm であり、実験報告の 0.26 nm – 0.28 nm [43] と一致している。O – O の PCF
は約 0.25 – 0.35 nm であり、トップピークは x が増加するに伴い短距離側へシフトした。

94 / 145

95 / 145

Figure 4 – 4. Pair correlation function, Si – O / Al – O (a) , Na – O (b) and O – O (c) of xSiO2-(1−x)NaAlO2 glasses at 300 K.

Figure 4 – 5 にモル体積を示す。300 K におけるモル体積は実験値の値 [40, 44 – 46] に比べ大きいが、0.5 ≤ x ≤ 1.0 において、x が増加するに伴いモル体積の値がわずかに減少する傾向は実験値と一致していると考えている。

This work
M. Taylor and G. E. Brown Jr. by exp. [40]F. Seifert et al. by exp. [44]A. Klonkowski [45]P. B. Narottam et al. by exp. [46]Molar volume, V / cm3 mol−1

39
37
35
33
31

29
27
25

0

0.2

0.4

0.6

0.8

1

Composition, x

Figure 4 – 5. Molar volumes of xSiO2-(1−x)NaAlO2 glasses at 300 K.

96 / 145

Figure 4 – 6 にガラスと結晶の NBO 比を示す。結晶はすべての x において、BO 比が 100 %を示した[24, 33, 34, 37]。ガラスの BO 比は x が減少するに伴い 100 % から減少した。NBO 比と酸素トリクラスター構造比は x が減少するに伴い 0 % から増加した。

100
BO(Crystals)

90
80
BO
70

%

60
50
40
30

NBO(Crystals)

20

NBO and tricluster

10
0

0.0

0.2

0.4

0.6

0.8

1.0

Composition, x

Figure 4 – 6. NBO, BO and tricluster ratio in all oxygen of xSiO2-(1 – x)NaAlO2 glasses at 300 K. Openand closed symbols are bridging oxygen of this work and crystals, respectively. Triangles and diamonds
show non-bridging oxygen and tricluster, respectively.

97 / 145

Figure 4 – 7 にリング構造を示す。リングサイズはリング中の Si 及び Al の数で定義した。Al
のみからなる 6 員環は x が増加するに伴い減少した。Si のみからなる 6 員環は x が増加するに伴い増加した。Si と Al からなる 4 員環と 6 員環は Si もしくは Al のみのリングよりも、多く存在しており、SiO2 – CaAl2O4 系のラマン分光測定より算出したリング構造の傾向と一致している [44].

(b) Experimental data [44](a) This work

Percentage of 4- and 6-membered ring (%)

90

Percentage of 4- (Si and Al) and 6- (Si or Al ) membered ring (%)

100

6-rings (Si and Al)

80
70
60

50
6-rings (Si)

6-rings (Al)40
30

4-rings (Si and Al)

20
10
4-rings (Al)4-rings (Si)

0
0.0

0.2

0.4

0.6

0.8

100

4-rings (Si and Al)

90
80
70
60

50

6-rings (Si)40
6-rings (Al)

30
20

10
0
0.0

1.0

0.2

Composition, x

0.4
0.6
Composition, x

0.8

1.0

Figure 4 – 7. Ring size distribution in xSiO2-(1−x)NaAlO2 glasses simulated at 300 K, (a) and in xSiO2(1−x)CaAl2O4 systems for experimental data, (b). Triangles, circles and squares show 6-membered ringsof Al, 6-membered rings of Si and 4-membered rings of Si/Al, respectively.

98 / 145§4 – 5. 考察
Table 4 – 4 において、AlO3 及び AlO5 ユニットがわずかに形成されたが、Xiang ら [17] のNPT アンサンブル下の MD シミュレーションによる報告では、x = 0.5 において全 Al イオンのうちの約 20 % 程が AlO3 ユニット、約 13 % 程が AlO5 ユニットとして形成されており、これらの結果と比較すると、本研究では Al の多面体構造の歪みは少なく、ほぼすべて AlO4 ユニットで形成されていることが考えられる。
Figure 4 – 4 (a) 及び (c) より、Si – O 及び、Al – O 間距離はそれぞれ四面体ユニットを構成する際の値を示した。Si – O、Al – O 及び O – O の PCF は x が増加するに伴い、AlO4 から SiO4
ユニットへ類似した変化を行うことで構造が安定化すると考えられる。Figure 4 – 4 (b) の Na – Oの PCF では 0.2 nm – 0.3 nm の分布を取り、実験報告 [43] で報告されている分布と一致していた。これらのことから、本シミュレーションは四面体ユニット内の短距離的な構造を再現できていると考えられる。
Figure 4 – 5 より、シミュレーションで得られたモル体積は文献値 [40, 44 – 46] よりも大きい値を示したが、x が増加するにつれ、モル体積が減少する傾向を再現できていることから、相対的な x の変化のシミュレーションは行えていると考えている。
Figure 4 – 6 より、結晶 [20 – 38] とガラスの BO 比が異なった。Stebbins と Xu [9, 47] は NMR
により、CaAl2Si2O8 ガラスと NaAlSi3O8 ガラスに NBO が存在することを報告した。彼らは、また、ガラス中に NBO と酸素トリクラスターの存在することを指摘している。Maekawa ら [48]も、29Si MAS NMR を用いることで、Na2O – SiO2 – Al2O3 系ガラスにおいて、Al/Si 比を増加させることで、NBO がわずかに存在することを報告している。Lacourse [49] は Na+ イオンが電荷補償により AlO4 ユニットを形成させ、かつネットワークを修飾する二つの役割があることを指摘している。ここで、Figure 4 – 8 に酸素トリクラスター構造に配位しているカチオンの種類を示す。酸素トリクラスターは 3 つの Al、1 つの Si 及び 2 つの Al の場合に形成されることが多い
ことが分かった。Stebbins と Xu [9, 47] は、3 つの Al、または、1 つの Si 及び 2 つの Al のとき、
酸素トリクラスターが表れることを指摘し、本研究結果と一致している。
本研究の結果より、NBO とトリクラスター構造が存在した。AlO4 ユニットがトリクラスター構造をもつ場合のユニット単位で取り得る電荷を BO と NBO のパターン別に考えると、トリク
ラスター構造をもち、他の三つの O イオンが BO であるとすると、[AlO4]−0.67 となる。一方、トリクラスター構造をもち、他の O イオンすべてが NBO であるとした場合、[AlO4]−3.67 となるが、トリクラスター構造を取り得る場合は、四面体同士が稜共有する構造をとるため、すべて NBOになる場合はほぼ存在しないと考えられる。NBO とトリクラスター構造の存在割合はほぼ同一という本研究の結果から、AlO4 ユニット内で NBO とトリクラスターが 1:1 の割合で共存するとしたとき、ユニット全体の電荷が−1.67 前後となるように帯電し、その AlO4 ユニットの近接に Na+ イオンが 1 ~ 2 個存在することが考えられる。もしくは、SiO4 ユニット中に NBO ができることで、系の電気的中性を保つ。この二つの構造パターンによって系の電気的中性が保たれ、構造が形成されていると考えられる。

99 / 145

100
3Al

Tricluster coordination species (%)

90
80
70
60
50
40

1Si2Al30

2Si1Al

20

3Si

10
0
0.0

0.2

0.4
0.6
Composition, x

0.8

1.0

Figure 4 – 8. Coordination species of oxygen tricluster in xSiO2-(1−x)NaAlO2 glasses simulated at 300 K.
Diamonds, squares, triangles and circles show 3Al, 1Si2Al, 2Si1Al and 3Si, respectively.

100 / 145

Figure 4 – 7 において、Si と Al からなるリングが多く存在しており、bond order が Si – O – Alもしくは Si – O – Si 構造となりやすいことと関連していると考えられている [8]。ここで、Figure
4 – 9 に Si – O – Si、 Al – O – Al 及び Si – O – Al 結合の割合を示す。点線は x = 0.5 のとき、Si –O – Al 結合が理想的に 100 % を取る場合の Si – O – Si 及び Al – O – Al 結合の分布を示す。ネッ
トワーク構造において、Si – O – Al の bond order を avoidance [8] と呼ぶ。Al – O – Al 結合は 0≤ x ≤ 0.8 において、理想的な avoidance よりも存在割合が多くなった。Si – O – Si 結合も 0.5 ≤ x< 0.8 において、理想的な avoidance よりも存在量が多くなる結果を示した。Lee と Stebbins [50,51] は Si – O – Si、Al – O – Al 及び Si – O – Al 結合の分布について 17O 3Q MAS NMR を用いて
Figure 4 – 9 (b)の結果を導いた。 彼らは、NMR の解析結果より、すべての x で BO のみを取る
四面体ユニットを形成すると仮定したモデルを構築した。彼らのモデルは、x = 0.5 において、わずかに Al – O – Al 結合の分布が NMR によって観測されることから、x = 0.5 における Si – O – Al結合の分布量が約 80 % を示すことを報告している [50, 51]。しかし、本研究の Si – O – Al 結合の分布量は彼らのモデルより少なく、また x = 0.67 において Si – O – Al 結合が最大量を示した。
理想的には、x = 0.5 において、Si – O – Al 結合が 100 % の分布を取ることが考えられるが、Figure
4 – 9 (a) の Si – O – Al 結合の割合が、それよりも小さく、x = 0.67 で最大を示した。
ここで Figure 4 – 10 – 1 及び – 2 にシミュレーションで得たエンタルピーと混合エンタルピーを示す。Figure 4 – 10 – 1 より、x が増加するにつれエンタルピーが減少した。Figure 4 – 10 – 2 より、0 < x ≤ 0.5 では混合エンタルピーが正の値を取ることから、エンタルピー的には不混和であることを示した。一方、0.5 < x では混合エンタルピーは負の値を示したことから、混合しやすい
結果を示した。Lee と Stebbins [51] は、17O NMR 測定結果から仮説モデルより導いた bond order(混合)のエンタルピーは、すべての x で負の傾向を示しているに対し、Navrotsky ら [52] は
NaAlSi3O8-SiO2 ガラス(本ガラス組成における、0 .75 ≤ x ≤ 1) において、ほぼゼロに近い値の混合エンタルピーを示すことをカロリメトリー測定より報告しており、混和もしくは不混和の両方を取り得ることを示唆している。本研究では混合エンタルピーが負を示す組成領域において、SiO2
rich 組成である x = 0.67 で Si – O – Al 結合が最も多くなった。本研究において、Lowenstein [8]の avoidance の説のように Si – O – Al 結合と Si – O – Si 結合が存在することは、エンタルピー的に混合しやすいことを表していると考えている。Lee と Stebbins のモデルは NBO と酸素トリク
ラスターの影響を考慮しておらず、これが本研究との差異になったことが原因として挙げられる。
対象としている実際のガラスにおいて、中間酸化物である Al2O3 のみではガラス化しづらいことから、本シミュレーションの混合エンタルピーはネットワーク構造の安定性を反映した結果を示している可能性がある。
x が増加すると、Figure 4 – 6 より BO が多くなり、エネルギー的に安定化する (Figure 4 – 10 –1)。このとき、一つの Na+ が[AlO4]1− と互いに電気的中を保つようにバランスを取りやすくなり、構造を形づくるか、SiO4 ユニット中の NBO が一つできる。一方、x が減少するとき、系の電気的中性を保つ手段として、（１） すべて BO で AlO4 ユニットとなるか、（２）NBO ができトリ
クラスターができない、（３）NBO とトリクラスターができる、パターンが考えられる。x = 0 である NaAlO2 結晶は （１）にあたるが、本研究では x が減少するほど、NBO とトリクラスターができる。ここで、NaAlO2 結晶の MD シミュレーションを行うと、NBO とトリクラスター構造
101 / 145

ができないことは確認している。本研究結果より、ガラスでは、（１）すべて BO で 4 配位構造を形成する、もしくは （３）NBO のときトリクラスターができない場合よりも、（２）NBO と
酸素トリクラスター構造を形成する方が存在する確率が高いことが示された。

(b) Experimental model [50, 51](a) This work
100

100
Al–O–Al9080

Si–O–Si (ideal avoidance)

80

Al–O–Al(ideal avoidance)

70

70

Si–O–Si

60

%

60

%

Si–O–Al

Al–O–Al90

50
Si–O–Al40

50
Si–O–Si

40

30

30

20

20

10

10
0

0
0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Composition, xComposition, x

Figure 4 – 9. Network linkage distribution of xSiO2-(1−x)NaAlO2 glasses, (a) simulated at 300 K, and (b)experimental model [50, 51]. Circles, squares and triangles show Si–O–Si, Al–O–Al and Si–O–Al,
respectively. Dashed line shows the ideal avoidance.

102 / 145

-4600-4650
Enthalpy, H / kJ mol−1

-4700

-4750
-4800

-4850
-4900

-4950
-5000

-5050
-5100

-5150
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 4 – 10 – 1. Enthalpy of xSiO2-(1−x)NaAlO2 glasses simulated at 300 K. Circles show thiswork at 300 K.

103 / 145

Mixing enthalpy, ΔHmix / kJ mol−160This work

40

20
0
-20
Exp. model [51]

-40
-60
-80
-100
-120
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 4 – 10 – 2. Enthalpy of mixing in xSiO2-(1−x)NaAlO2 glasses simulated at 300 K. Opencircles and solid circles show this work at 300 K and experimental model at 300 K [51].

104 / 145§4 – 6. 結言
第一原理計算により定めた原子間相互作用を適用し、xSiO2 - (1 – x)NaAlO2 ガラスのシミュレー
ションを行った。この原子間相互作用を適用することにより、MD シミュレーションによる原子間距
離、モル体積、リング構造は実験報告の組成変化に伴う傾向と一致した。本研究では、NPT アンサン
ブル下でのガラスの構造をシミュレーションで表現できた。
本研究はすべての x で SiO4 及び AlO4 ユニットが形成され、酸素トリクラスター構造と NBO が現れた。酸素トリクラスター構造が現れる場合、3 つの Al、または、1 つの Si 及び 2 つの Al のときに形成されることが多いことが分かった。これは NMR による報告と一致した傾向を示した。このとき、AlO4 ユニット内にトリクラスターが 1 つ、BO が 2 つ、NBO が一つとなる[AlO4]−1.67 にNa+ イオンが 1 ~ 2 つ近接する場合か、AlO4 ユニット内にトリクラスターが一つ、BO が 3 つの[AlO4]−0.67 で SiO4 ユニット内に NBO ができる場合が、
電気的中性を保つ構造として考えられる。
x が減少するとき、系の電気的中性を保つ手段として、NBO と酸素トリクラスター構造を形成することが本研究により示された。

105 / 145§4 – 7. 参考文献
[1] A. Ellison and A. Cornejo, Glass substrates for liquid crystal displays, Int. J. Appl. Glass Sci. 1 (2010)87 – 103.
[2] J. C. Mauro, D. C. Allan, M. Potuzak, Nonequilibrium viscosity of glass, Phys. Rev. B 80 (2009) 094204.
[3] E. D. Lacy, Aluminum in glasses and in melts, Phys. Chem. Glasses 4 (1963) 234 – 238.
[4] M. Schmücker, K. J. D. MacKenzie, H. Schneider, R. Meinhold, NMR studies on rapidly solidified
SiO2-Al2O3 and SiO2-Al2O3-Na2O-glasses, J. Non-Cryst. Solids 217 (1997) 99 – 105.
[5] D. A. McKeown, F. L. Galeener, G. E. Brown, Raman studies of Al coordination in silica-rich sodium
aluminosilicate glasses and some related minerals, J. Non-Cryst. Solids 68 (1984) 361 – 378.
[6] M. Okuno, N. Zotov, M. Schmücker, H. Schneider, Structure of SiO2-Al2O3 glasses: Combined X-ray
diffraction, IR and Raman studies, J. Non-Cryst. Solids 351 (2005) 1032 – 1038
[7] H . Morikawa, S. Miwa, M. Miyake, F. Marumo, T. Sata, Structural analysis of SiO2-Al2O3 glasses,
J. Am. Ceram. Soc. 65 (1982) 78 – 81.
[8] W. Loewenstein, The distribution of aluminum in the tetrahedra of silicates and aluminates, Am. Min.
39 (1954) 92 – 96.
[9] J. F. Stebbins and Z. Xu, NMR evidence for excess non bridging oxygen in aluminosilicate glass, Nature
390 (1997) 60 – 62.
[10] M. Montorsi, M. C. Menziani, C. Leonelli, G. C. Pellacani, A. N. Cormack, Molecular dynamics
simulations of alumina addition in sodium silicate glasses, Mol. Simul. 24 (2000) 157 – 165.
[11] C. Leonelli, G. Lusvardi, M. Montorsi, M. C. Menziani, L. Menabue, P. Mustarelli, Influence of Small
Additions of Al2O3 on the Properties of the Na2O·3SiO2Glass, J. Phys. Chem., B 105 (2001) 919 –
927.
[12] C. A. Angell, P. A. Cheeseman, R. R. Kadiyala, Diffusivity and thermodynamic properties of diopsideand jadeite melts by computer simulation studies, Chem. Geology 62 (1987) 83 – 92.
[13]

D. J. Stein and F. J. Spera, Molecular dynamics simulations of liquids and glasses in the system
NaAlSiO4-SiO2: Methodology and melt structures, Am. Mine. 80 (1995) 417 – 431.

[14] F. Noritake, K. Kawamura, K. N. Matsukage, Elastic anomalies of anorthite: Molecular dynamics
simulations, Physics of the Earth and Planetary Interiors 244 (2015) 32 – 41.
[15] A. Winkler, J. Horbach, W. Kob, K. Binder, Structure and diffusion in amorphous aluminum silicate:a molecular dynamics computer simulation, J. Chem. Phys. 120 (2004) 384 – 393.
[16] Xin Dai, Jian He, Jin Bai, Qing Huang, Xiaodong Wen, Lu Xie, Kan Kuo, Jian Zhang, Wen Li, Shiyu
Du, Ash fusion properties from molecular dynamics simulation: roles of the ratio of silicon andaluminum, Energy & Fuels 30 (2016) 2407 – 2413.
[17] Y. Xiang, J. Du, M. M. Smedskjaer, J. C. Mauro, Structure and properties of sodium aluminosilicate
glasses from molecular dynamics simulations, J. Chem. Phys. 139 (2013) 044507.
[18] Y. Yamamoto, N. Sawaguchi, M. Sasaki, A new determination method of interatomic potential for
sodium silicate glass simulations, J. Non-Cryst. Solids 466 – 467 (2017) 29 – 36.
[19] Y. Yamamoto, N. Sawaguchi, M. Sasaki, A study of mixed alkali effect of Na2O-K2O-SiO2 glasses
106 / 145

by molecular dynamics simulation, J. Comput. Chem. Jpn. 16 (2017) 98 – 101.
[20]

V. L. Cherginets,V. N. Baumer, S. S. Galkin, L. V. Glushkova, T. P. Rebrova, Z. V. Shtitelman,
Solubility of Al2O3 in Some Schloride-Fluoride Melts, Inorg. Chem. 45 (2006) 7367 – 7371.

[21] M. Marezio and J. P. Remeika, High-Pressure Synthesis and Crystal Structure of α-LiAlO2, J. Chem.
Phys. 44 (1966) 3143 – 3145.
[22]

F. Stewner, R. Hoppe, Zur Kristallstruktur von α-Li5AlO4, Zeitschrift fuer Anorganische undAllgemeine Chemie 380 (1971) 241 – 243.

[23] N. J. Zhu, F. Guo, Y. Li, S. X. Shen, L. R. Chen, A. C. Li, Study on crystal structure of diaoyudaoite,
Hua Hsueh Hsueh Pao 50 (1992) 527 – 532.
[24] J. A. Kaduk and S. Y. Pei, The crystal structure of hydrated sodium aluminate, NaAlO2 . 5/4(H2O),
and its dehydration product, J. Solid State Chem. 115 (1995) 126 – 139.
[25] M. G. Barker, P. G. Gadd, M. J. Begley, Preparation and crystal structures of the first alkali-rich sodium
aluminates Na7Al3O8 and Na5AlO4, J. Chem. Soc., Chem. Commun. 0 (1981) 379 – 381.
[26] M. G. Barker P. G. Gadd, S. C. Wallwork, A new sodium aluminate Na17Al5O16, J. Chem. Soc. Chem.
Commun. 1982 (1982) 516 – 517.
[27] F. Griesfeller, J. Koehler, R. Hoppe, Das erste Oligo-Oxo-Aluminat: Na14Al4O13, Zeitschrift fuer
Anorganische und Allgemeine Chemie 507 (1983) 155 – 162.
[28] V. I. Voronin, M. G. Surkova, Sh. G. Shekhtman, N. A. Anurova, V. A. Blatov, Conduction mechanismin the low-temperature phase of KAlO2, Inorg. Mat. 46 (2010) 1234 – 1241.
[29] A. Bon, C. Gleitzer, A. Courtois, Synthese et structure cristalline d'un nouvel aluminate de potassium,
J. Protas, Comptes Rendus Hebdomadaires des Seances de l'Academie des Sciences, Serie C, Sciences
Chimiques 278 (1974) 785 – 788.
[30] K. F. Hesse, Crystal structures of natural and synthetic alpha eucryptite, LiAlSiO4, Zeitschrift fuer
Kristallographie 172 (1985) 147 – 151.
[31]

G. J. Redhammer and G. Roth, Structural variation and crystal chemistry of Li Me(3+) Si2 O6
clinopyroxenes Me(3+) = Al, Ga, Cr, V, Fe, Sc and In. Zeitschrift fuer Kristallographie 219 (2004)
278 – 294.

[32] A. Zemann-Hedlik and J. Zemann, Die Kristallstruktur von Petalit, LiAlSi4O10, Acta Crystallographica
8 (1955) 781 – 787.
[33]

A. Kumar, S. J. Dhoble, D. R. Peshwe, J. Bhatt, Structural and Photoluminescence properties ofnepheline-structure NaAlSiO4 : Dy(3+) nanophosphors, J. Alloy Compd. 609 (2014) 100 – 106.

[34] A. C. McCarthy, R. T. Downs, R. M. Thompson, Compressibility trends of the clinopyroxenes and insitu high-pressure single-crystal X-ray diffraction study of jadeite, Am. Min. 93 (2008) 198 – 209.
[35] M. Gregorkiewitz, Y. Li, T. J. White, R. L. Withers, I. Sobrados, The structure of "orthorhombic" K Al
Si O4 - O1: evidence for Al-Si order from MAS NMR data combined with Rietveld refinement andelectron microscopy, Canadian Mineralogist 46 (2008) 1511 – 1526.
[36] G. D. Getta, N. Rotiroti, T. B. Boffa, A. Pavese, Leucite at high pressure: Elastic behavior, phasestability, and petrological implications, Am. Min. 93 (2008) 1588 – 1596.
107 / 145

[37] A.F. Wright and M.S. Lehmann, The structure of quartz at 25 and 590 ℃ determined by neutron
diffraction, J. Solid State Chem. 36 (1981) 371 – 380.
[38] E. Zintl, A. Harder, B. Dauth, Gitterstruktur der oxide, sulfide, selenide und telluride des lithiums,
natriums und kaliums, Zeitschrift fuer Elektrochemie 40 (1934) 588 – 593.
[39] K. Kawamura, MXDORTO, Japan Chemistry Program Exchange P029, (1996).
[40] M. Taylor and G. E. Brown Jr., Structure of mineral glasses-Ⅱ. The SiO2-NaAlSiO4 join, Geochim.
Cosmochim. Acta 43 (1979) 1467 – 1473.
[41] G. N. Greaves, A. Fontanie, P. Lagarde, D. Raoux, S. J. Gurman, Local structure of silicate glasses
Nature 293 (1981) 611 – 616.
[42] D. A. McKeown, G. A. Waychunas, G. E. Brown, Jr., EXAFS study of the coordination environmentof aluminum in a series of silica-rich glasses and selected minerals within the Na2O-Al2O3-SiO2
system, J. Non-Cryst. Solids 74 (1985) 349 – 371.
[43] S. K. Lee and J. F. Stebbins, The distribution of sodium ions in aluminosilicate glasses: a high-field
Na-23 MAS and 3Q MAS NMR study, Geochim. Cosmochim. Acta 67 (2003) 1699 – 1709.
[44] F. Seifert, B. O. Mysen, D. Virgo, Three-dimensional network structure of quenched melts (glass) in
the systems SiO2-NaAlO2, SiO2-CaAl2O4 and SiO2-MgAl2O4, Am. Min. 67 (1982) 696 – 717.
[45] A. Klonkowski, Band character and properties of the mixed alkali system Na2O-K2O-Al2O3-SiO2, J.
Non-Cryst. Solids 57 (1983) 339 – 353.
[46] P.B. Narottam, R.H. Doremus, Handbook of glass properties, Acdemic Orlando, (1986).
[47] S. K. Lee and J. F. Stebbins, The distribution of sodium ions in aluminosilicate glasses: A high-field
Na-23 MAS and 3Q MAS NMR study, Geochim. Cosmochim. Acta 67 (2003) 1699 – 1709.
[48] H. Maekawa, T. Maekawa, K. Kawamura, T. Yokokawa, 29Si MAS NMR investigation of the Na2OAl2O3-SiO2 glasses, J. Phys. Chem. 95 (1991) 6822 – 6827.
[49] W. C. Lacourse, Structural influences on diffusion in glass – The mixed site effect, J. Non-Cryst. Solids
21 (1976) 431 – 434.
[50] S. K. Lee and J. F. Stebbins, Al-O-Al and Si-O-Si sites in framework aluminosilicate glasses with
Si/Al=1:quantification of framework disorder, J. Non-Cryst. Solids 270 (2000) 260 – 264.
[51] S. K. Lee and J. F. Stebbins, Extent of intermixing among framework units in silicate glasses and melt,
Geochim. Cosmochim. Acta 66 (2002) 303 – 309.
[52] A, Navrotsky, R. Hon, D. F. Weill, D. J. Henry, Thermochemistry of glasses and liquids in the systems
CaMgSi2O6-CaAl2Si2O8-NaAlSi3O8, SiO2-CaAl2Si2O8-NaAlSi3O8 and SiO2-Al2O3-CaO-Na2O,
Geochim. Cosmochim. Acta 44 (1980) 1409 – 1423.

108 / 145

第 5 章 xNa2O - (1 – x)B2O3 ガラスのホウ酸異常現象の再現性の評価§5 – 1. 緒言
アルカリホウ酸塩系ガラスは化学的耐久性や熱衝撃性に優れているため、実用ガラスの基本組成として用いられることが多い。xNa2O - (1 − x)B2O3 ガラス中に存在する全ホウ素に対する四配位ホウ素の割合(B4) はアルカリ金属酸化物量, x の導入により増加する。これは四配位ホウ素が微弱な負の電荷を帯びているため、Na+ イオンがこの四配位ホウ素にトラップされる形で近くに存在すると考えられている。xNa2O - (1 − x)B2O3 ガラスの場合、x = 0 において BO3 ユニット(B3) 構造のみを取り、NWM 量, x を増加させると Figure 1 – 3 に示したような B3 – O1− … Na1+ 構造とはならず、B3 が BO4 ユニット(B4) へと変化し、x ≤ 0.40 まで B4 – O – B3 (or B4) 構造を取り、x> 0.40 になると B4 – O±0 – B3 (or B4) 構造が減少し、B3 – O1− … Na1+ 構造が増加することが 11B
NMR 測定により分かっている (Figure 5 – 1) [1]。

1/2Na2O

B

O

O

O

1/2Na2O

1–

O

O

B

Na+

O

O

O

O

Na+

1–

B

O

B

O

1/2Na2O

O

O

B
O

O

1–

1–

O O
Na+

Na+

Figure 5 – 1. Coordination number change of boron unit with Na2O addition.

109 / 145

O
BO

一方、xNa2O - (1 − x)B2O3 ガラスは約 x = 0.20 で線熱膨張係数が極小値を持つことが報告されており [2, 3]、これをホウ酸異常と呼ぶが、x = 0.4 で極大となる B4 と x = 0.2 で極小となる線熱膨張係数の極値傾向が一致しないことから、より大きな構造ユニットとなる B – O − リング構造
が線熱膨張係数と関係があるとされ、議論されてきた。その一つとして、Krogh-Moe [4] は、x = 0のとき、B イオンと O イオンがそれぞれ三つからなる 3 員環 (ボロキソールリング) のみからな
り、x が増加するにつれ員環数が大きくなるモデルを提案した。Figure 5 – 2 に Krogh-Moe が提案したモデルを示す。

Boroxol ring

Penta-borate

Tetra-borate

Di-borate
O

O
O
O
O

B

B

B
O

O

O

B

B
O

O

O B

O

O

B
O

O

O
O B
O

B

O
B
O
B

O B

O
B
O

O

O

O
O B
O

O

B

O

B

O

O

O

O B

B2O3

O
B

O

B

O

Na2O

Figure 5 – 2 . Krogh-Moe’s ring and ring group model [4].

110 / 145

O

O

B

O

O

Krogh-Moe は Na2O/B2O3 比が増加するに伴い、ボロキソールリング、ペンタボレート、テトラ
ボレート及びダイボレートグループに分けてモデルを示した。これらはすべてボロキソールリン
グを基本構造としたリング構造となっている。 ラマン分光測定 [5 – 9]、XRD 及び中性子回折 [10– 15]、NMR [16 – 18] などの測定により Krogh-Moe の仮説モデルは支持されている。一方で、B2O3
ガラスにおいてボロキソールリングは BO3 平面構造であるため、実測よりもガラスの密度が小さくなることを指摘し、Krogh-Moe のモデルを否定する説もある [19] 。しかし、17O DOR NMR 測
定により、ボロキソールリングが系中に約 70 % あり、少し大きい員環で構成されるルーズなボ
ロキソールリングが約 30 % 存在するという報告 [20] から、ボロキソールリングを基本構造としたリング構造によってホウ酸異常を解釈することが一般的となっている。このとき、BO3 平面構
造との間に van der Waals 力が働くことによりガラス構造が形成されているとして解釈することが一般的となっている [3]。
このリング構造について、MD シミュレーションからの研究も行われており、Xu ら [21] は二体間相互作用のみを用いてアルカリホウ酸塩ガラスの第一近接範囲の構造を良く再現したがボロキソールリングは現れなかった。ボロキソールリングが MD シミュレーションで現れなかったのは、B イオンと O イオンの間に働く共有結合が弱かったからであるとして、Takada ら [22] は二体間相互作用に加え、角度と二面角を補強する三体間相互作用と四体間相互作用を適用した MD
シミュレーションを行ったが、ボロキソールリングは実験結果よりも少ないことが報告されている。他にも様々な原子間相互作用を適用する試みがあったが、ボロキソールリングは形成されな
い、もしくは少量しか現れないことが報告されている [21 – 26]。一方で、これらの MD の研究では線熱膨張係数については言及されていない。従来の MD シミュレーションで適用している原子間相互作用では、線熱膨張係数が再現できなかったことが考えられる。第 2 章において、第一原
理計算の結果を利用し原子間相互作用を設定する方法を考案した。そして、第 3, 4 章においてこの設定方法による原子間相互作用を定め、Na2O-SiO2 系ガラス [27, 28]、Na2O-K2O-SiO2 系ガラス
[29] 及び Na2O-Al2O3-SiO2 系ガラスの MD シミュレーションを行い、ガラスの構造や物性の再現性が良いことを報告した。第 2 – 4 章では、Si – O ネットワーク構造を主とするケイ酸塩系に対し、原子間相互作用の設定方法について考察してきた。本論文では、様々な酸化物ガラス組成へ対応可能な原子間相互作用の設定方法を開発することを目的としている。そのため、NWF が複数存在する酸化物ガラスについても原子間相互作用が適応できているかを検証する必要がある。
そこで、本章では B – O ネットワーク構造を形成する xNa2O-(1 − x)B2O3 ガラス(0 ≤ x ≤ 0.7)を対象とし、第 2 章で考案した設定方法による原子間相互作用を適用した MD シミュレーショ
ンを行い、ガラスの構造、及びホウ酸異常の再現性の評価を行った。
§5 – 2 では、適用した原子間相互作用の設定方法について記した。§5 – 3 では、MD シミュレー
ションの計算条件について示した。§5 – 4 では MD シミュレーションを行った結果について示した。
§5 – 5 では§5 – 4 で得られた結果について考察した。§5 – 6 では本章の結言を示した。§5 – 7 に本章で引用した参考文献を示す。

111 / 145

§5 – 2. 原子間相互作用
本研究は、シンプルに二体間相互作用のみで、ガラスの構造をどこまで再現できるかを調べることを目的としている。二体間相互作用は式 (2 – 1) を使用した。また、従来の Na2O-B2O3 系ガ
ラスの MD シミュレーションの大半は二体間相互作用に加え、角度を補正する三体間相互作用を適用することが多い [ 22]。これは二体間相互作用のみでは共有結合力による方向性をもつ力を表現できないとするという考えから、適用される。本研究では、第 2 章で考案した設定方法により定めた二体間相互作用に、三体間相互作用を加えた際のガラス構造の影響についても調査することにした。本研究で適用した三体間相互作用は Takada ら [22] を参考に、O – B – O 角が 120 °となるように、式(5 – 1) [30] を使用した。Figure 5 – 3 に三体間相互作用を示す。
(5 – 1)

𝑈𝑖𝑗𝑘 = −𝑓[cos{2(𝜃𝑘𝑖𝑗 − 𝜃0 )} − 1]√𝑘𝑖𝑗 𝑘𝑖𝑘
𝑘𝑖𝑗 =
𝑘𝑖𝑘 =

1
exp[𝑔𝑟 (𝑟𝑖𝑗 −𝑟𝑚 )]+1
1
exp[𝑔𝑟 (𝑟𝑖𝑘 −𝑟𝑚 )]+1

f は定数、𝜃𝑘𝑖𝑗 は原子 k-i-j 間に相互作用を働かせる角度、 𝜃0、 𝑔𝑟 及び 𝑟𝑚 は共有結合の角度項に関するパラメータである。

0.1
0
-0.1
-0.2
-0.3

-0.4
90

190
0.23

0.22

170

0.21

0.2

150
0.19

0.18

0.16

0.17

0.15

110

130

210

Figure 5 – 3. Three-body, O – B – O interatomic potential by eq. (5 – 1) in this work.

112 / 145

Potential, UOBO / kJ unit -1

0.2

電荷は第 2 章の計算条件を用いて xA2O-(1 – x)B2O3 結晶 (A = Li, Na 及び K)に対し電荷解析を行った。ノルム保存型擬ポテンシャルは、Li, 2s1; Na, 2s22p63s1; K, 3s23p64s1; B, 2s22p1; 及び O, 2s22p4を使用した。Table 5 – 1 に計算対象として選んだ結晶 [31 – 46] と k 点セットを示す。Figure 5 –4 に xA2O-(1 – x)B2O3 結晶 (A = Li, Na 及び K) の電荷解析結果を示す。x が増加するに伴い B イ
オンと O イオンの電荷が減少する傾向を示した。この電荷解析の傾向を基に、xNa2O-(1 – x)B2O3
系中のイオンの電荷を、以下の式 (5 – 2) - (5 – 4) で示した。
𝑧𝑁𝑎 = −0.25𝑥 + 1.0 (0 < 𝑥 ≤ 1)

(5 – 2)𝑧𝐵 = −0.25𝑥 + 1.5 (0 ≤ 𝑥 < 1)

(5 – 3)𝑧𝑂 = neutrality(5 – 4)zO で系の電気的中性を保たせている。BO+ クラスターモデルを作成し、第 2 章で示した計算条件により MO 計算を行った。Figure 5 – 5 に x = 0 における B – O のカーブフィッティングの結果を示す。aNa, bNa, cNa 及び三体間相互作用のポテンシャルパラメータは経験的に決めた値を使用した。Table 5 – 2 に各組成におけるポテンシャルパラメータを示す。

113 / 145

Table 5 – 1 . Percent deviation, which gives the difference in structural optimization by DFT. Deviation shows in absolute value.

.

114 / 145

2.0
1.5

Alkali, A

Mulliken charge

1.0
0.5

B

0.0
-0.5
-1.0

O

-1.5
-2.0
0

0.2

0.4
0.6
Composition, x

0.8

1

Figure 5 – 4. Mulliken charges of ions in xA2O-(1−x)B2O3 crystals; A = Li (open symbols), Na (solid
symbols) and K (dotted symbols). Boron, alkali and oxygen ions show the triangles, squares andcircles, respectively.

115 / 145

Energy difference, ΔE ×1018 / J

0.016
0.014
0.012

BO+ model
B

O

0.010
0.008
0.006
0.004

MO

0.002

MD

0.000

0.115

0.120
0.125
Distance, rBO / nm

0.130

Figure 5 – 5. Energy surface determined from MO calculation (MP2/6-311+g(d)) of an BO+ model andpotential energies calculated by eq. (2 – 1) with fitted parameters for x = 0.

116 / 145

Table 5 – 2. Interatomic potential parameters for the simulation of xNa2O-(1−x)B2O3 glasses.
(a) x = 0
Atomic parametersOz-1.0000

a / nm0.17920c / (kJ/mol)1/2 nm30.05605b / nm
0.01368B

1.50000.11749

0.01198

0.00000

Pair parameters

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1

B−O83680.00

59.29

33472.00

29.48

c / (kJ/mol)1/2 nm30.05605(b) x = 0.1
Atomic parametersza / nmb / nmO-1.0179

0.179200.01368B

1.47500.117490.119780.00000

Pair parameters

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1

B−O83680.00

59.29

33472.00

29.48

z-1.0385
1.4500D 1ij / kJ mol−183680.00

a / nm0.179200.11749
β 1ij / nm−1
59.29

b / nm
0.01368
0.01198
D 2ij / kJ mol−1
33472.00

c / (kJ/mol)1/2 nm30.05605
0.00000
β 2ij / nm−1
29.48

Atomic parametersza / nmb / nmO
B-1.0625
1.42500.179200.11749

0.01367
0.01197c / (kJ/mol)1/2 nm30.05605
0.00000

Pair parameters

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1

B−O83680.00

59.29

33472.00

29.48

c / (kJ/mol) nm
0.05605(c) x = 0.2
Atomic parameters
O
B
Pair parameters
B−O(d) x = 0.3(e) x = 0.4
Atomic parametersza / nmb / nmO-1.09090.179200.013671/2

B

1.40000.11749

0.011970.00000

Pair parameters

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1

B−O83680.00

59.29

33472.00

29.48

3(f) x = 0.5
Atomic parametersza / nmb / nmc / (kJ/mol)1/2 nm3O-1.12500.17876

0.01337

0.05605B

1.3750Pair parameters

D 1ij / kJ mol

β 1ij / nm

0.11755

D 2ij / kJ mol

β 2ij / nm

B−O83680.00

59.33

33472.00

30.00

c / (kJ/mol) nm
0.05605−1−10.012310.00000−1−1(g) x = 0.6
Atomic parametersza / nmb / nmO-1.16670.17875

0.013351/2

3B

1.35000.11755

0.012300.00000

Pair parameters

D 1ij / kJ mol−1

β 1ij / nm−1

D 2ij / kJ mol−1

β 2ij / nm−1

B−O83680.00

59.33

33472.00

30.00

c / (kJ/mol) nm
0.05605(h) x = 0.7
Atomic parametersza / nmb / nmO-1.21880.17876

0.01333B

1.32500

0.11755

0.01228

−1

−1

1/2

3

0.00000−1

Pair parameters

D 1ij / kJ mol

β 1ij / nm

D 2ij / kJ mol

β 2ij / nm−1

B−O83680.00

59.33

33472.00

30.00

Atomic parametersza / nmb / nmc / (kJ/mol)1/2 nm3(i) Na and O−B−O

Na－0.25x +1.0

0.132200.01150

0.01227

3-body parametersf / 10−22 kJ

θ 0 / deg.

r m / nmg r / nm−1

O−B−O

0.35

120.0

0.1880198.0

117 / 145

§5 – 3. MD シミュレーションの計算条件
MD シミュレーションは MXDORTO システム [47] を使用した。B イオンは非常に軽く小刻みな振動を行うことが予想される。そのため時間刻みは 0.4 fs とした。NPT アンサンブルのシミュレーションを行い、一辺が約 5.0 nm の立方体 MD セル中に疑似乱数を用いて粒子(N = 5400 ~ 5700) をランダムに配置させ、圧力 (P = 0.1 MPa ) および温度 (T = 1800 K) の融体のシミュレーションを行い、500 000 step にかけて全エネルギーが一定になるまで構造緩和を行った。本章で扱う xNa2O-(1 − x)B2O3 ガラス中の各イオンの粒子数を Table 5 – 3 に示す。
ガラス状態は冷却速度を 0.01 K/step とし、段階的に冷却することで得た。1800 – 300 K では
500 000 steps の構造緩和を行った。解析時間はすべての温度で、10 000 steps とした。本章では、二体間相互作用のみ (2bp) と三体間相互作用も適用した場合 (2+3bp) のシミュレーショ
ンを行った。

Table 5 – 3. Number of particles (N) in xNa2O-(1−x)B2O3 glass simulations.

x
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7

O
3240
3150
3120
3000
2915
2700
2700
2400

B
2160
2025
1920
1750
1590
1350
1200
900

118 / 145

Na
0
225
480
750
1060
1350
1800
2100

Total (N )
5400
5400
5520
5500
5565
5400
5700
5400

§5 – 4. 結果
Figure 5 – 6 に 300 K における PCF を示す。2bp 及び 2+3bp によるシミュレーションにおいて、第一近接 B – O 間の PCF は 0.143 – 0.190 nm の分布を示している。 Na – O、O – O、
B – B 間の PCF は 0.20 – 0.30 nm、0.22 – 0.30 nm、及び 0.20 – 0.32 nm の分布を示した。ナト
リウムホウ酸塩結晶 [31, 37 – 42] の第一近接 B – O、Na – O、O – O 及び B – B 間の PCF はそれぞれ、0.126 – 0.158 nm、0.221 – 0.314 nm、0.226 – 0.245 nm 及び 0.227 – 0.266 nm であり、
xNa2O - (1 − x)B2O3 ガラスの XRD による B – O 間距離は 0.137 nm (x = 0) [10]、0.140 nm (x =0.25) [11]、0.145 nm (x = 0.33) [12]、Na – O 間距離は 0.240 nm (x = 0.33) [12] である。2bp 及び 2+3bp によるシミュレーションはこれらと一致している。

119 / 145

(a) 2bpx

Na − Na

0

0.1
B−B

0.2

PCF

0.3
0.4
0.5

O −O

0.6

x
Na − O

0.7

B−O

0

0.1

0.2
0.3
0.4
Distance, rij / nm

0.5

0.6

(b) 2+3bpx
Na − Na
0
0.1
B−B

0.2

PCF

0.3
0.4
0.5

O −O
x

0.6

Na − O

0.7

B−O

0

0.1

0.2
0.3
0.4
Distance, rij / nm

0.5

0.6

Figure 5 – 6. Pair correlation function for xNa2O-(1−x)B2O3 glasses simulated by 2bp, (a) and 2+3bp,(b) at 300 K. x = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and 0.7 show black, red, yellow, green, blue, paple, pink
and yellow green, respectively.

120 / 145

Figure 5 – 7 に 300 K におけるモル体積の結果を示す。2bp 及び 2+3bp によるシミュレー
ションは x が増加するに伴い、モル体積が減少する実験の傾向 [12, 48 – 52] と一致した。2bp
及び 2+3bp のモル体積の値は実験値よりも大きい値を示した。また、2bp 及び 2+3bp のモル
体積は x = 0.5 において極大値を示した。

55
Kunugi et al. [48]
Kodama [49]M olar volume, V / cm3 mol-1

50Inaba et al. [50]K, -D. OH et al. [12]Franz et al. [51]45

Eversteijn et al. [52]
40

2+3bp

35

30
2bp

25
0

0.2

0.4
Composition, x

0.6

0.8

Figure 5 – 7. Molar volumes of xNa2O-(1−x)B2O3 glasses simulated by 2bp and 2+3bp at 300 K.

121 / 145

Figure 5 – 8 に温度変化に伴うモル体積の傾向を示す。2bp 及び 2+3bp のモル体積は温度が増加するに伴いモル体積が増加する実験結果の傾向 [12, 48 – 52] と一致した。ガラス転移温度をモル体積の温度変化に伴う傾向から見積もった。2bp シミュレーションによるガラス転移温度は x = 0 ~ 0.7,でそれぞれ、およそ 700, 800, 800, 700, 700, 600, 600 及び 600 K であった。2+3bp シミュレーションによるガラス転移温度は x = 0 ~ 0.7 でそれぞれ、およそ 600, 800,
800, 700, 700, 600, 600 及び 600 K であった。Shelby [2] は x = 0.27 でガラス転移温度が最大値を取るホウ酸異常現象を示すことが報告しており、本研究の 2bp 及び 2+3bp によるシミュレーションはこの実験報告の傾向と一致した。

122 / 145

(a) 2bp (x = 0, 0.1 and 0.2)

60

60
55

Molar volume, V / cm3 mol-1

Molar volume, V / cm3 mol-1

55(d) 2+3bp (x = 0, 0.1 and 0.2)

50
x=0

45

40

x = 0.1

3530

x = 0.2

25

x=0

50
45

x = 0.1

40
35

x = 0.2
30
25

20

0

500

1000

1500

20

2000

0

500

Temperature, T / K
45

45

43

43

41

41

Molar volume, V / cm3 mol-1

Molar volume, V / cm3 mol-1(b) 2bp (x = 0.3 and 0.4)

39
37
35
33
x = 0.3
31

29

1500

2000

1500

2000

(e) 2+3bp (x = 0.3 and 0.4)

39
37
35
x = 0.333
31

x = 0.429

x = 0.4
27

27

25

25
0

500

1000

1500

2000

0

500

Temperature, T / K

1000

Temperature, T / K(c) 2bp (x = 0.5, 0.6 and 0.7)

45

45

43

(f) 2+3bp (x = 0.5, 0.6 and 0.7)

4341

x = 0.6

39
37
35
33
x = 0.5

31
29

Molar volume, V / cm3 mol-1

Molar volume, V / cm3 mol-1

1000

Temperature, T / K

41
x = 0.6
39
37
35

33
x = 0.53129

x = 0.7

x = 0.7

27

27

25

25
0

500

1000

1500

2000

Temperature, T / K

0

500

1000

1500

2000

Temperature, T / K

Figure 5 – 8. Molar volumes of xNa2O-(1−x)B2O3 glasses/melts simulated by (a) – (c) 2bp and (d) – (f)2+3bp. Circles, triangles, squares, diamonds, crosses, asterisks, pluses and minuse signs show x = 0, 0.1,
0.2, 0.3, 0.4, 0.5, 0.6 and 0.7, respectively.

123 / 145

Figure 5 – 9 に、300 K における B4 比を示す。0 ≤ x < 0.4 において、2bp 及び 2+3bp シミュレーションでは、B4 が増加する傾向 [1] を再現し、x ≤ 0.4 において両シミュレーション
は B4 が減少する傾向 [1] を再現した。NMR による実験値 [1] と比較して、2bp によるシミュレーションは B4 の値が大きくなり、2+3bp によるシミュレーションは B4 が小さくなった。

100

Bray et al. (Li) [1]
Bray et al. (Na) [1]
Bray et al. (K) [1]9080

Theoretical curve

B4 (%)

70
2bp

60
50
40
30
20

2+3bp

10
0
0

0.2

0.4
Composition, x

0.6

0.8

Figure 5 – 9. B4 of xNa2O - (1 − x)B2O3 glasses at 300 K. Open and solid symbols show the MD results
and experimental values (solid circles for Na2O-B2O3 glasses, solid triangles for Li2O-B2O3 glasses,
and solid squares for K2O-B2O3 glasses), respectively. Theoretical curve represents the relation B4 =100x/(1 − x) [1].

124 / 145

Figure 5 – 10 (a) ~ (d) に O – B – O 及び B – O – B 角度の分布を示す。2bp 及び 2+3bp シミュレーションによる O – B – O 角は約 120 °を中心とした分布を示し、x = 0.1 ~ 0.5 における
BO4 ユニットが形成するときに示される約 110 °の分布も見られた。2bp による B – O – B
角の分布は各組成で 80 ≤ θ ≤ 180 °分布を示し、90 °付近にホウ素多面体同士が稜共有構
造を取ることで現れるピークが現れた。一方で、2+3bp では O – B – O 角を 120 °として補正した力を加えているため、稜共有構造が現れなかったと考えられる。2bp 及び 2+3bp による B – O - B 角は x が増加するに伴い、120°から 160° を示す分布が多くなった。

125 / 145

(a) 2bp, O − B − O

Distribution

x = 0.7
0.6
0.5

0.4
0.3

0.2
0.1

0
80

100

120

140

160

Angle, θO−B−O / deg.
(b) 2+3bp, O − B − O

Distribution

x = 0.7
0.6
0.5
0.4

0.3
0.2

0.1
0

80

100

120

140

160

Angle, θO−B−O / deg.

Figure 5 – 10. Distribution of O−B−O (θO−B−O) and B−O−B (θB−O−B) bond angles in xNa2O-(1−x)B2O3 glasses
simulated at 300 K.

126 / 145

(c) 2bp, B − O −B

Distribution

x = 0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

20

40

60

80

100

120

140

160

180

140

160

180

Angle, θB−O−B / deg.

(d) 2+3bp, B − O −B

Distribution

x = 0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

20

40

60

80

100

120

Angle, θB−O−B / deg.

Figure 5 – 10 (continued). Distribution of O−B−O (θO−B−O) and B−O−B (θB−O−B) bond angles in
xNa2O-(1−x)B2O3 glasses simulated at 300 K.

127 / 145

線熱膨張係数, αは
𝛼=

1

d𝑙

∙

(5 – 5)𝑙300 K d𝑇から算出した。l は立方体 MD セルの一辺の長さであり、l300K は 300 K における MD セルの長さを示す。式 (5 – 5) より、温度, T = 300 K と 600 K の MD セル長から見積もった線熱膨張係数の結果を Figure 5 – 11 に示す。2bp 及び 2+3bp を適用したシミュレーションは x =0.2 で最小値を示す実験結果の傾向 [2, 3] を再現した。実験で報告されている 293 K から 573
K の線熱膨張係数の値よりも本研究の値は大きくなった。

Linear thermal expansion coefficient, α / ppm K-1

60

Shelby (373K - 573K) [2]
Klyuev et al. (293K - 573K) [3]

50

40
2bp
30

20

2+3bp

10

0
0

0.2

0.4
Composition, x

0.6

0.8

Figure 5 – 11. Coefficients of linear thermal expansion in xNa2O-(1−x)B2O3 glasses smulated at
300 K and 600 K.

128 / 145

Figure 5 – 12 に B – O リング中のホウ素イオンの数で区別したリング構造の存在量を示す。
2bp 及び 2+3bp シミュレーションは共に、x < 0.4 では 5 ~ 8 員環が多く、3 員環 (ボロキソールリング) 及び 4 員環は x が増加するに伴い x = 0.4 で最大値を取り、0.4 < x で減少した。
x = 0.2 では、中距離構造は 5 ~ 8 員環により構成され、ボロキソールリングはほぼ形成されていない。2bp シミュレーションは 2+3bp に比べ、x = 0.4 において 2 及び 3 員環の数が多かった。

(a) 2bp(b) 2+3bp
18

45

2

40

3

Counts for rings per B 100 particles

Counts for rings per B 100 particles

50

4

35

5

30

6

25

7

20

8

15
10

16

2

14

3
4

12

5
10

6

8

7

6

8

4
2

5
0

0

0

0.2

0.4
0.6
Composition, x

0.8

0

0.2

0.4
0.6
Composition, x

0.8

Figure 5 – 12. Ring size distributions in xNa2O-(1−x)B2O3 glasses simulated at 300 K using (a) 2bp and(b) 2+3bp.

129 / 145

§5 – 5. 考察
Figure 5 – 13 に BO、NBO 及び三配位酸素 (O3) 比を示す。
、x = 0 において、2bp シミュレーションでは NBO と O3 が現れ、2+3bp シミュレーションでは BO のみ形成された。0 ≤ x≤ 0.5 において、2bp 及び 2+3bp によるシミュレーションは、BO が減少し、NBO が増加する傾向を示した。0 ≤ x ≤ 0.33 において結晶では、BO のみしか現れず、ガラスの NMR の報告 [1] においても、BO4 ユニットが形成されることで NBO は出現しないとされている。
2bp 及び 2+3bp シミュレーションでは NBO も出現した。2+3bp によるシミュレーションに
おいて、角度に対して力を加えても 2bp と同様に NBO が現れた。これは、2bp では二体間相互作用のみを与え、各元素に対し固定した電荷を与えていることから、局所的な電荷移動の表現ができないため出現したと解釈している。0 ≤ x ≤ 0.5 において、2bp シミュレーショ
ンでは O3 が最大 20.0 % 現れ、2+3bp シミュレーションでは、最大 1.86% 現れた。Zachariasenの理論に基づいて考えるならば、O3 は多面体が稜共有することで現れるため、本来出現しな
いと考えられている。本シミュレーションで O3 が現れるのは、0 ≤ x ≤ 0.33 に NBO の出
現と同様の理由と解釈している。角度に対する力を加えることで 2+3bp シミュレーションの方が、2bp シミュレーションよりも三配位酸素(O3)を抑制できているが、本研究では、構造が分からない未知の組成への MD シミュレーションを行うことを目的に原子間相互作用の設定
方法を考えてきた。そのため、三体間相互作用を適用する場合、既に結合角度が明らかになっている必要があり、未知の組成への援用がしにくい点がある。このことを踏まえると、2bpは二体間の相互作用のみでも x > 0.33 におけるナトリウムホウ酸塩結晶 [31, 37 – 42, 46] のBO、NBO 比の傾向を再現できていることから、定性的な議論を行えるレベルではシミュレーションができていると考えている。

130 / 145

BO (cry.)NBO (cry.)100
90

BO (2+3bp)

80
BO (2bp)

70

(%)

60
50

NBO (2bp)40

NBO (2+3bp)30

O3 (2bp)20
O3 (2+3bp)
10
0
0

0.2

0.4
Composition, x

0.6

0.8

Figure 5 – 13. Bridging oxygen (BO), non-bridging oxygen (NBO) and three coordination oxygen (O3)in xNa2O-(1−x)B2O3 glasses simulated at 300 K. Triangles, circles and squares show the BO, NBO andO3 ratio, respectively. Solid, dashed and dated lines show the 2bp, 2+3bp and xNa2O-(1−x)B2O3 crystals[31, 37 – 42, 46], respectively.

131 / 145

Figure 5 – 13 及び 5 – 14 に O イオン周辺の配位構造を示す。Figure 5 – 13 に第一近接配位圏の Na – NBO, Na – BO 及び Na – O3 間距離を示す。このとき、解析に用いたカットオフ
距離は 0.3 nm とした。x = 0 のとき、Na+ イオンが存在しないので、グラフに表示していな
い。2bp 及び 2+3bp によるシミュレーションで得られた Na – NBO 間距離は x が増加するに伴い長距離側へ変化する傾向を示し、Na – BO 間距離は x の増加に伴いわずかに低距離側へ変化する傾向を示した。2bp 及び 2+3bp によるシミュレーションで得られた Na – O3 間距離
も Na – BO 間距離と同様の傾向を示した。Na – NBO 及び Na – BO 間距離よりも Na – O3 間
距離の方が長距離側へ現れるのは、O3 が NBO と BO に比べポテンシャル的に負に帯電しづらいことを反映しているためであると考えている。
Figure 5 – 14 に、O イオン周辺に存在する Na+ イオンの配位数を示す。すべての O イオ
ン周辺の Na+ イオンの配位数は x が増加するに伴い、増加する傾向を示しており、これは、系中の Na+ イオンの存在量に対応していると考えている。すべての x において、NBO の周辺に Na+ イオンが存在する割合が多いことが分かった。O3 の周辺に Na+ イオンの存在量が少ないのは、系中に O3 が少量しか存在しないためである。Figure 5 – 13 及び 5 – 14 の結果
より、Figure 3 – 10, 3 – 11 で示したように Na+ イオンが BO よりも負に帯電する NBO に対し、選択的に結合することが分かった。本研究では、Na2O-B2O3 系ガラスにおいても、2bp のようにシンプルな原子間相互作用を適用したシミュレーションにおいても、詳細な構造解析を可能とすることが期待できる。Figure 5 – 6 に示した Na – O 間の PCF では x が増加する
につれ、短距離側へのシフトが見られており、これは x が増加するに伴い NBO の周辺に Na+
イオンが配位する割合が多くなることによって表れている傾向だと考えている。Figure 5 – 6の 0.2 nm 付近における B – B 間の PCF は BO4 ユニットが形成されることでネットワーク
構造が密となり、ホウ素多面体間の距離が短くなることから現れたピークであると考えられる。Na – Na 間の PCF は x が増加するに伴い、短距離側へシフトする傾向を示した。これは、系中の B イオンが少なくなり、ネットワーク構造を形成せず、孤立した BO3 ユニットと Na+
イオンが多くなり、Na+ イオン同士が近接する割合が多くなるためであると考えられる。

132 / 145

0.290Na - O3 (2+3bp)

Distance of Na - oxygen, d / nm

0.280
Na - O3 (2bp)

0.270
Na - BO (2+3bp)

0.260
0.250

Na - BO (2bp)

Na - NBO (2bp)

0.240
0.230
Na - NBO (2+3bp)

0.220
0

0.2

0.4
Composition, x

0.6

0.8

.
Figure 5 – 13. Distance of Na+ ions – BO and NBO in xNa2O-(1−x)B2O3 glasses
simulated at 300 K. Circles and triangles show 2bp and 2+3bp simulations,
respectively. Solid, broken and dotted line show – BO, – NBO and – O3 , respectively

133 / 145

Coordination number of Na around oxygen (counts)

5.0
4.5
NBO - Na (2bp)

NBO - Na (2+3bp)

4.0

BO - Na (2bp)

3.5
3.0

BO - Na (2+3bp)

2.5

O3 - Na (2bp)

2.0
1.5
O3 - Na (2+3bp)

1.0

0.5
0.0
0

0.2

0.4
Composition, x

0.6

0.8

Figure 5 – 14. Coordination number of BO and NBO – alkali metal ions in xNa2O-(1−x)B2O3 glasses
simulated at 300 K. Circles and triangles show 2bp and 2+3bp simulations, respectively. Solid,
broken and dotted line show BO –, NBO – and O3 – , respectively

134 / 145

Figure 5 – 7 より、2bp シミュレーションは 2+3bp シミュレーションよりもモル体積が小さくなった。これは 2+3bp シミュレーションは 2bp よりも B4 が少ないため、系中の密度が増加するためである。2bp 及び 2+3bp を適用したシミュレーションは x が増加するに伴いモル体積が減少する傾向を再現した。x = 0.5 において、2bp 及び 2+3bp を適用したシミュレーシ
ョンのモル体積は極大傾向を示した。これは、Figure 5 – 9 より、0.4 ≤ x ≤ 0.5 では B4 が減少することでモル体積が増加し、0.5 < x では、BO4 ユニットが BO3 ユニットとなり、NBOが増加することで、モル体積が減少する傾向を示したと考えている。2bp によるシミュレーシ
ョンで得られたガラス転移温度は、2+3bp を適用したシミュレーションと同様の組成変化に対する傾向と値を示していることから、x に対して定性的なシミュレーションが行えていると考えている。
Figure 5 – 9 に示した B4 のように、2+3bp によるシミュレーションでは実験値よりも B4が少なくなるため、必ずしも良い結果を与えていると考えてはいない。本章のように、組成変化によって酸素多面体の構造が変化するガラスを対象とする場合、角度の補正を加える組
成とそうではない組成が存在する場合が出てくる。このとき、ある組成には角度を補正する
力を加え、別の組成には二体間相互作用のみを加えることは、材料設計の視点で考えた場合、有効とは言えない。シミュレーションの結果と実験結果では、見かけの物性値が一致しない
背景がそれぞれに存在すると考えている。本研究による、2bp 及び 2+3bp を適用したシミュ
レーションはいずれも、実験報告の B4 の値とは一致はしなかったが、x 変化による傾向は再現できている。特に、2bp によるシミュレーションにおいて、角度の補正を加えずとも、B4 の傾向を定性的な結果として得ることができていることから、シンプルな原子間相互作用でも酸素多面体構造が変化する組成のシミュレーションにも対応できていると考えている。
Figure 5 – 10 において 2bp の O – B – O 角は 100 ~ 120 ° 分布を示した。これは。2+3bpと比較すると、2bp の方が B4 が多いことから、B4 の傾向に対応した分布を示していると考えている。2bp シミュレーションにおける、約 90 °付近の B – O – B 角は O3 の存在と対応している。2bp 及び 2+3bp シミュレーションにおいて、B – O – B 角は x が増加するに伴い小さくなっており、この傾向は Na2O – SiO2 系ガラスにも見られた (Figure 2 – 9)。このとき、Na2O – SiO2 系ガラスの場合、第一原理計算から組成変化に伴う各イオンの電荷を与えており、ガラス中の Si – O 結合の弱化と関係した構造変化を示した。Toyuki と Akagi はアルカリホ
ウ酸塩系ガラスにおいて、B4 が増加すると、B – O – B 結合における BO 上のπ電子密度が増加することを示しており、これは B – O 間結合の弱化と関係するという報告をした [53 –56]。本章の Na2O – B2O3 系ガラスにおいても第 2 章と類似した組成変化を示す電荷を適用していることから、B – O 結合の弱化を反映した構造変化を示したと考えている。
Figure 5 – 11 より、両シミュレーションは実験値よりも線熱膨張係数の値が大きくなったが、B4 の実験結果とシミュレーションの結果の差から考えると、この線熱膨張係数の実験値とシミュレーション結果の差は生じてしまうものとして解釈している。このことを踏まえると、2bp 及び 2+3bp を適用したシミュレーションはいずれも、線熱膨張係数のホウ酸異常の傾向を示していると考えている。
Figure 5 – 15, 16 に x 変化に対応する 2bp 及び 2+3bp シミュレーションで得られたネット
135 / 145

ワーク構造を示す。2bp 及び 2+3bp シミュレーションにおいて、0 ≤ x < 0.2 で 5 ~ 8 員環が存在し、0.2 ≤ x < 0.5 で 3 及び 4 員環構造が増加することからネットワーク構造が密になっていることが確認できる。0.5 ≤ x では、3 及び 4 員環が主に存在し、ネットワーク構造を形作るよりも単体のリング構造や BO3 ユニットとして存在することが多い。Figure 5 – 12 より、2bp 及び 2+3bp シミュレーションにおいて、x = 0 のとき、ボロキソールリング及び 4 員
環は系中でほとんど存在せず、5 ~ 8 員環以上のリング構造が存在することを示した。0 < x <0.2 では、3 ~ 8 員環が増加した。0.2 ≤ x ≤ 0.3 では 5 ~ 8 員環の存在量が最も多くなり、0.3 < x ≤ 0.4 で 5 ~ 8 員環は x の増加に伴い減少し、ボロキソールリング及び 4 員環は増加する傾向を示した。x = 0.4 のときボロキソールリング及び 4 員環の存在量が最も多くなっており、両シミュレーションとも類似した傾向を示した。これらの結果は Krogh-Moe [4] のモ
デルで示されているボロキソールリングを主とした組成変化の傾向に反している。
一方で、本研究では組成変化における線熱膨張係数の非線形な傾向 [2, 3] を再現できており、ホウ酸異常は系中のボロキソールリング量が少なくても現れた。線熱膨張係数が極小となる x = 0.2 では、Figure 5 – 2 に示すようなテトラボレート群とダイボレート群の存在がラマン分光測定により確認されており、これらの構造ユニットの存在により熱膨張係数が低下するという考えがある [8, 9]。本シミュレーション結果においても、これらの構造ユニットの存在が熱膨張係数の低下と関連している可能性も考えられるが、本研究の x = 0.2 では、他の組成に比べ系中で 5 ~ 8 員環の存在量が多く、ボロキソールリングはこれらのリング構造と比較すると少量であった。これはボロキソールリングのように小さいリング構造ではなく、ネットワーク構造を形作る大きいリング構造が、線熱膨張係数のホウ酸異常と関係があると考えられる。このような構造ユニットのモデル化のためには、電荷移動も考慮できる詳細な構造解析結果を導く必要があり、そのためには、第一原理計算と組み合わせたシミュレーシ
ョンを行うことが求められる。本研究では、Na2O-B2O3 系ガラスの構造をシミュレーション出来ていることから、この緩和構造を基に大規模な第一原理計算を行うことでホウ酸異常を解き明かす基礎を作ることができたと考えられる。0.4 < x ではすべてのリングの存在量が減
少しており、B – O − ネットワーク構造が切断されていることが考えられる。

136 / 145

137 / 145

0.5 nm(f) x = 0.5(b) x = 0.1

Na

BO3

(g) x = 0.6(c) x = 0.2BO4(h) x = 0.7(d) x = 0.3

Figure 5 – 15. Snapshot of structure in xNa2O-(1−x)B2O3obtained by 2bp simulation (MD cell was sliced at 1.0 nm) (drawn by VESTA [57]).

(e) x = 0.4(a) x = 0138 / 145

0.5 nm(f) x = 0.5(b) x = 0.1

Na

BO 3

(g) x = 0.6(c) x = 0.2

BO 4(h) x = 0.7(d) x = 0.3

Figure 5 – 16. Snapshot of structure in xNa2O-(1−x)B2O3obtained by 2+3bp simulation (MD cell was sliced at 1.0 nm) (drawn by VESTA [57]).

(e) x = 0.4(a) x = 0§5 – 6. 結言
第一原理計算を利用して設定した原子間相互作用を適用した MD シミュレーションを用いて、xNa2O - (1 − x)B2O3 ガラス (0 ≤ x ≤ 0.7) におけるホウ酸異常を再現できるかを調査した。
xNa2O - (1 − x)B2O3 ガラス中の原子間距離、B4 の組成変化の傾向を再現した。系中のボロキソー
ルリング量は実験で報告されている量よりも少なかったが、ホウ酸異常と呼ばれる組成変化に対する線熱膨張係数の非線形的な傾向も再現した。本研究の x = 0.2 では、他の組成域に比べ員環数のバリエーションが多かった。ボロキソールリングのように小さいリング構造ではなく、ネット
ワーク構造を形作る大きいリング構造が系全体に連なることがホウ酸異常と関係があると考えられる。

139 / 145

§5 – 7. 参考文献
[1] P. J. Bray and J. G. O’Keefe, Nuclear magnetic resonance investigations of the structure of alkali borate
glasses, Phys. Chem. Glasses 4 (1963) 37 – 46.
[2] J. E. Shelby, Thermal expansion of alkali borate glasses, J. Am. Ceram. Soc. 66 (1983) 225 – 227.
[3] V. P. Klyuev and B. Z. Pevzner, The influence of aluminum oxide on the thermal expansion, glasstransition temperature, and viscosity of lithium and sodium aluminoborate glasses, Glass Phys. Chem.
28 (2002) 207 – 220.
[4] J. Krogh-Moe, The structure of vitreous and liquid boron oxide, J. Non-Cryst. Solids 1 (1969) 269 –284.
[5] B. N. Meera and J. Ramakrishna, Raman spectral studies of borate glasses, J. Non-Cryst. Solids 159(1993) 1 – 21.
[6] A. A. Osipov and L. M. Osipova, Structure of glasses and melts in the Na2O – B2O3 system from hightemperature raman spectroscopic data: 1. Influence of temperature on the local structure of glasses andmelts, Glass Phys. Chem. 35 (2009) 121 – 131.
[7] C. F. Windisch Jr. and W. M. Risen Jr., Vibrational spectra of oxygen- and boron-isotopically substituted
B2O3 glasses, J. Non-Cryst. Solids 48 (1982) 307 – 323.
[8] T. Yano, N. Kunimine, S. Shibata, M. Yamane, Structural investigation of sodium borate glasses andmelts by Raman spectroscopy. 1. Quantitative evaluation of structural units, J. Non-Cryst. Solids 321(2003) 137 – 146.
[9] T. Yano, N. Kunimine, S. Shibata. M. Yamane, Structural investigation of sodium borate glasses andmelts by Raman spectroscopy. 2. Conversion between BO4 and BO2O- units at high temperature, J.
Non-Cryst. Solids 321 (2003) 147 – 156.
[10] R. L. Mozzi and B. E. Warren, The structure of vitreous boron oxide, J. Appl. Cryst. 3 (1970) 251 –
257.
[11] K. – D. OH, H. Morikawa, F. Marumo, A. Nukui, Structural analysis of 3B2O3・Na2O glass, YogyoKyokai-Shi 90 (1982) 270 – 275.
[12] K. – D. OH, H. Morikawa, M. Okuno, F. Marumo, Structural analysis of 2B2O3・Na2O glass, YogyoKyokai-Shi 90 (1982) 175 – 180.
[13] Y. Iwadate, K. Igarashi, T. Hattori, S. Nishiyama, K. Fukushima, J. Mochinaga, Short range structureof B2O3-Cs2O glasses analyzed by x-ray diffraction and Raman spectroscopy, J. Chem. Phys. 99(1993) 6890 – 6896.
[14] K. Handa, Y. Kita, S. Kohara, K. Suzuya, T. Fukunaga, M. Misawa, T. Iida, H. Iwasaki, N. Umesaki,
Structure of M2O-B2O3 (M: Na and K) glasses and melts by neutron diffraction, J. Phys. Chem. Solids
60 (1999) 1465 – 1471.
[15]

A. C. Wright, B. A. Shakmatkin, N. M. Vedishcheva, The chemical structure of oxide glasses: aconcept consistent with neutron scattering studies ?, Glass Phys. Cgem. 27 (2001) 97 – 113.

[16] Y. H. Yun and P. J. Bray, Nuclear magnetic resonance studies of the glasses in the system Na2O –B2O3 – SiO2, J. Non-Cryst. Solids 27 (1978) 363 – 380.
140 / 145

[17] S. Sen and J. F. Stebbins, Na-ion transport in borate and germanate glasses and liquids: a 23Na and
11B NMR spin-lattice-relaxation study, Phys. Rev. B 55 (1997) 3512 – 3519.
[18] J. Zhong and P. J. Bray, Change in boron coordination in alkali borate glasses, and mixed alkali effects,
as elucidated by NMR, J. Non-Cryst. Solids 111 (1989) 67 – 76.
[19]

S. R. Elliott, A continuous random network approach to the structure of vitreous boron trioxide,
Philosophical Magazine B 37 (1978) 435 – 446.

[20] R. E. Youngman, Short- and intermediate-range structural ordering in glassy boron oxide, Science 269(1995) 1416 – 1420.
[21] Q. Xu, K. Kawamura, T. Yokokawa, Molecular dynamics calculations for boron oxide and sodium
borate glasses, J. Non-Cryst. Solids 104 (1988) 261 – 272.
[22] A. Takada, C. R. A. Catlow, G. D. Price, Computer modelling of B2O3: part 1. New interatomic
potentials, crystalline phases and predicted polymorphs, J. Phys.:Condens. Matter 7 (1995) 8659 –
8692.
[23] A. Takada, C. R. A. Catlow, G. D. Price, Computer modelling of B2O3: part 2. Molecular dynamics
simulations of vitreous structures, J. Phys.:Condens. Matter 7 (1995) 8693 – 8722.
[24] A. Takada and A. N. Cormack, Computer simulation models of glass structure, Phys. Chem. Glasses:
Eur. J. Glass Sci. Technol. B 49 (2008) 127 – 135.
[25] K. Hirao and N. Soga, Molecular dynamics simulation of Eu3+ -doped sodium borate glasses and theirfluorescence spectra, J. Am. Ceram. Soc. 68 (1985) 515 – 521.
[26] N. Ohtori, K. Takase, I. Akiyama, Y. Suzuki, K. Hamda, I. Sakai, Y. Iwadate, T. Fukunaga, N. Umesaki,
Short-range structure of alkaline-earth borate glasses by pulsed neutron diffraction and molecular
dynamics simulation, J. Non-Cryst. Solids 293 – 295 (2001) 136 – 145.
[27]

Y. Yamamoto, N. Sawaguchi, M. Sasaki, Interatomic potential for xNa2O-(1–x)SiO2 glasses, J.
Comput. Chem., Jpn. 14 (2015) 63 – 66.

[28] Y. Yamamoto, N. Sawaguchi, M. Sasaki, A new determination method of interatomic potential for
sodium silicate glass simulations, J. Non-Cryst. Solids 466 – 467 (2017) 29 – 36.
[29] Y. Yamamoto, N. Sawaguchi, M. Sasaki, A study of mixed alkali effect of Na2O-K2O-SiO2 glasses
by molecular dynamics simulation, J. Comput. Chem. Jpn. 16 (2017) 98 – 101.
[30] H. Sakuma, K. Kawamura, Structure and dynamics of water on muscovite mica surfaces, Geochim.
Cosmochim. Acta 73 (2009) 4100 – 4110. / N. Kumagai, K. Kawamura, T. Yokokawa, An interatomic
potential model for H2O: applications to water and ice polymorphs, Mol. Sim. 12 (1994) 177 – 186.
[31] S. L. Strong and R. Laplow, The structure of crystalline B2O3, Acta Cryst. B24 (1968) 1032 – 1036.
[32] Yu. F. Shepelev, R. S. Bubnova, S. K. Filatov, N. A. Sennova, N. A. Pilneva, LiB3O5 crystal structure
at 20, 227 and 377 ℃, J. Solid State Chem. 178 (2005) 2987 – 2997.
[33] A. Jiang, S. Lei, Q. Huang, T. Chen, D. Ke, Structure of lithium heptaborate, Li3B7O12, Acta Cryst.
C46 (1990) 1999 – 2001.
[34] N. Sennova, R. S. Bubnova, G. Cordier, B. Albert, S. K. Filatov, L. Isaenko, Temperature-dependentchanges of the crystal structure of Li2B4O7, Zeitschrift fuer Anorganische und Allgemeine Chemie
141 / 145

634 (2008) 2601 – 2607.
[35] W. H. Zachariasen, The crystal structure of lithium metaborate, Acta Cryst. 17 (1964) 749 – 751.
[36] F. Stewner, Die Kristallstruktur von alpha-Li3BO3, Acta Cryst. B27 (1971) 904 – 910.
[37] R. S. Bubnova, Ju F. Shepelev, N. A. Sennova, S. K. Filatov, Thermal behaviour of the rigid boronoxygen groups in the alpha-(Na2B8O13) crystal structure, Zeitschrift fuer Kristallographie 217 (2002)
444 – 450.
[38] N. Penin, M. Touboul, G. Nowogrocki, Crystal structure of two new sodium borates Na3B7O12 andNa2Tl2, Journal of Alloys and Compounds 363 (2004) 104 – 111.
[39] J. Krogh-Moe, The crystal structure of sodium diborate Na2O (B2O3)2, Acta Cryst. B30 (1974) 578– 582.
[40] M. Marezio, H. A. Plettinger, W. H. Zachariasen, The bond lengths in the sodium metaborate structure,
Acta Cryst. 16 (1963) 594 – 595.
[41] H. Koenig, R. Hoppe, M. Jansen, Zur Kenntnis von Na4B2O5, Zeitschrift fuer Anorganische undAllgemeine Chemie 449 (1979) 91 – 101.
[42] H. Koenig and R. Hoppe, Zur Kenntnis von Na3BO3, Zeitschrift fuer Anorganische und Allgemeine
Chemie 434 (1977) 225 – 232.
[43] J. Krogh-Moe, The crystal structure of the high-temperature modification of potassium pentaborate,
Acta Cryst. B28 (1972) 168 – 172.
[44] R. S. Bubnova, V. S. Fundamenskii, S. K. Filatov, I. G. Polyakova, Crystal structure and thermal
behavior of KB3O5, Doklady Akad. Nauk 398 (2004) 643 – 647.
[45] J. Krogh-Moe The structure of potassium diboratem K2O2B2O3, Acta Cryst. B28 (1972) 3089 – 3093.
[46] E. Zintl, A. Harder, B. Dauth, Gitterstruktur der oxide, sulfide, selenide und telluride des lithiums,
natriums und kaliums, Zeitschrift fuer Elektrochemie 40 (1934) 588 – 593.
[47] K. Kawamura, MXDORTO, Japan Chemistry Program Exchange P029, (1996).
[48] M. Kunugi, A. Konishi, S. Takeuchi, T. Yamate, Density of alkali borate glasses, Material: technicalnotes 230 (1971) 978 – 981.
[49] M. Kodama, Ultrasonic velocity in sodium borate glasses, J. Mat. Sci. 26 (1991) 4048 – 4053.
[50]

S. Inaba, S. Oda, K. Morinaga, Equation for estimating the thermal diffusivity, specific heat andthermal conductivity of oxide glasses, J. Jpn. Inst. Metals 65 (2001) 680 – 687.

[51] H. Franz, “Borate glasses”, L. D. Pye, V. D. Frèchette and N. J. Kreidl ed., Plenum Press (1977).
[52] F. C. Eversteijn, J. M. Stevens, H. I. Waterman, The density, refractive index and specific refraction ofvitreous boron oxide and of sodium borate glasses as functions of composition, method of preparation,and rate cooling, Phys. Chem. Glasses 1 (1960) 123 – 133.
[53] S. Sakka and K. Matsusita, Studies on binary silicate glasses based on the SiKα and SiKβ emission Xrays, J. Non-Cryst. Solids 22 (1976) 57 – 66.
[54] H. Toyuki and S. Akagi, EPR and optical spectra of the vanadyl ion in alkali borate glasses, Phys.
Chem. Glasses 13 (1972) 15 – 19.
[55] H. Toyuki and S. Akagi, EPR and optical spectra of VO2+ ion in xKF-(1 – x)B2O3 glasses, Phys.
142 / 145

Chem. Glasses 15 (1974) 1 – 5.
[56]

H. Toyuki, Studies on bonding mode of Al atoms in K2O-AlF3-B2O3 glasses based on opticalabsorption of VO2+ and AlKα X-ray emission spectra, Yogyo-Kyokai-Shi 85 (1977) 19 – 22.

[57] K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric andmorphology data, J. Appl. Crystallogr. 44 (2011) 1272 – 1276.

143 / 145

第 6 章 総括
以下に、本研究の総括を述べる。
第 1 章では、本研究の背景及び目的について述べた。本研究は、従来の設定方法の課題であった原子間相互作用に使用する経験的な値を減らし、多種のガラス組成への適用性を併せ持つ、原子間相互作用の設定方法を研究することが目的である。
第 2 章では、実用ガラスの基本組成である xNa2O-(1 – x)SiO2 ガラス (0 ≤ x ≤ 0.5) を対象に、原子間相互作用の設定方法を考案し、その手法により定めた原子間相互作用を適用したシミュレーションを行い、計算結果について検証した。Coulomb 力のイオンの電荷とイオン間に働くポテ
ンシャルエネルギーに対し、二つの第一原理計算の手法を取り入れる方法を提案した。本研究による設定方法を用いた原子間相互作用を適用したシミュレーションは、従来の経験的設定方法と同等のシミュレーションの結果を示すだけでなく、従来に比べガラスのネットワーク構造の表現性が改善することを示した。
第 3 章では、混合アルカリ効果を示す 0.25{xK2O-(1 – x)Na2O}-0.75SiO2 ガラス (0 ≤ x ≤ 1) を対象とした。本研究で考案した原子間相互作用を適用すると混合アルカリ効果を再現可能である
ことを示した。また、従来の MD シミュレーションの研究報告では混合アルカリ効果と関係がないとされていた、アルカリ金属イオンの質量差も影響を及ぼすことを明らかにした。また、混合
アルカリ効果が起きるとき、異種アルカリ金属イオンによって拡散が阻害される back tracking
motion が起きていることをダイナミクスと熱力学的視点から示した。
第 4 章では、xSiO2 - (1 – x)NaAlO2 ガラス (0 ≤ x ≤ 1) を対象とした。アルミノケイ酸塩は実用材料ではスマートフォンなどのカバーガラスに使用され、地球惑星科学の分野ではマグマの主成分として知られている。本研究で行ったシミュレーションは xSiO2 - (1 – x)NaAlO2 ガラス中の原子間距離、リング構造、モル体積の再現性が良いことを明らかにした。網目形成酸化物が二種類混合しているアルミノケイ酸ガラスについても適用可能であることを示した。
第 5 章では、xNa2O - (1 – x)B2O3 ガラス (0 ≤ x ≤ 0.7) を対象に第 2 章で考案した原子間相互作
用の設定方法の組成への汎用性について明らかにした。xNa2O - (1 – x)B2O3 ガラスは、これまで多くの MD シミュレーションが行われてきたが、組成に対し線熱膨張係数が非線形的な傾向を示す、
ホウ酸異常と呼ばれる現象の表現性が低いことが知られていた。熱膨張係数はガラスを製造する上で重要な指標となる物性値の一つであり、この線熱膨張係数を本研究で適用した原子間相互作
用によるシミュレーションで再現できれば、MD シミュレーションがガラス材料の設計ツールとしての活用を期待できると考えた。本研究によるシミュレーションでは、原子間距離、モル体積、BO4 ユニット構造比、線熱膨張係数を表現できることを明らかにした。
第 1 章から第 5 章より、第一原理計算の結果を基に原子間相互作用を設定したことにより、原子間相互作用の変数に使用する経験的な値を減らし、ガラス組成への適用性を向上させたシミュ
レーションを可能にした。本研究で開発した原子間相互作用の設定方法は、多種の酸化物ガラスに対して有効であることが期待できる。本研究は、MD シミュレーションによる酸化物ガラスの材料設計ツールの基盤を作ることができたと考えている。

144 / 145

謝辞
本論文は筆者が、室蘭工業大学大学院工学研究科

工学専攻 博士後期課程

機能材料学研究室に在籍中の研究成果をまとめたものである。機能材料学研究室、澤口 直哉 准教授には本研究の実施の機会を与えて戴き、その遂行にあたって、ご指導戴きました。ここに感謝の意を表します。機能材料学研究室、佐々木 眞 教授には本論文の細部にわたり、ご助言とご指導を戴きました。ここに深謝の意を表します。同専攻、岸本 弘立
ここに感謝の意を表します。東京工業大学

教授には副査としてご助言を戴きました。

河村 雄行 名誉教授には、本論文で取り扱ったシミュレーション手法について、多くのご助言を戴きました。深く感謝の意を表します。機能材料学
研究室の各位には、研究遂行にあたり日頃より有益なご討論ご助言を戴きました。ここに感謝の意を表します。
2018/12/21 山本 優也

145 / 145