分子動力学法による B2O3 ガラスの構造シミュレーション
無機材料化学研究室 大 09-137

榮 秀和

1. 緒言
今日まで、望む機能を持つ材料を設計する「材料設計」は、幾度も実験を繰り返し試行錯誤しつつ、機能と組成との関連について相関性を見出すことが一般的であったが、効率よく「材料設計」を行うべく、近年、コンピューターシミュレーションを構造‑物性‑機能間の相関性の予測手段、もしくはその解析手法として用いる動きが活発となってきた。材料化学の分野で用いられるコンピューターシミュレーションの計算手法としては、多電子系を取り扱う分子軌道法(MO)や、多原子系を取り扱う分子シミュレーション法がある。その中でも多原子系を取り扱う分子動力学法(molecular dynamics method; MD)は、与えられた原子間相互作用を用いてそれぞれの原子に対して働く力を求め、古典的なニュートンの運動方程式を順次解析的に解くことによって全ての原子の位置と速度を時間の関数として得るシミュレーション法であり、広い範囲の温度や圧力における系の構造・物性が予測可能である。さらに、非現実的な環境下における挙動や実験では観察し得ない情報も得られるなど、材料設計を行う上で非常に有用な研究手法であるといえる。
本研究で対象とした B2O3 ガラスの構造に関する研究はこれまでにも多くの研究者によって行われてきた。当初、Zachariasen のガラス形成則により B2O3 ガラスは BO3 平面三角形が頂点を共有し合ったランダムな網目構造[Fig.1(a)]から形成されているとされていたが、後の回折法 1),2)、および分光法 3)-6)による研究により、3 つの BO3 平面三角形が同一平面内で重合したボロキソール環[Fig.1(b)]と呼ばれる特異な中距離構造グループが存在するとの考え方が主流となっている。しかし、MD 法による B2O3 ガラスの構造に関する過去の研究によると、ボロキソール環を再現したという報告は論者の知る限り Inoue ら 7)、Takada ら 8)によるものだけであり、ボロキソール環を再現することは非常に困難であるとされてきた。このような背景からもシミュレーションによるボロキソール環の再現については計算化学に残された大きな課題であると言える。
またホウ酸塩ガラスにはアルカリ金属酸化物等の(a)修飾酸化物の添加に伴いボロキソール環を含む様々な特異構造が存在することが示唆されており、同時に組成に対して密度、熱膨張係数、屈折率等の物性が非直線的に変化する「ホウ酸異常」と呼ばれる現象がよく知られている。この特異構造と物性の組成依存性には密接な相関関係があり、シミュレーショ
B3O6ンによりこの特異構造の再現に挑戦することは、シ(b)（ボロキソール環ール環）
ミュレーションによるガラスの材料設計の実現においても大きな意義がある。
したがって、本研究では B2O3 ガラスの実測データを満足するガラス構造モデルを得るべく、3 体ポテンシャルを導入した MD 法を用いてホウ酸塩結晶の構造・物性を維持するポテンシャルを探索し、得Fig.1 random network model(a)られたポテンシャルを用いて B2O3 ガラスの再現を試
and boroxol ring model(b).
1み、B2O3 ガラス構造モデルについて検討を加えた。また、ここでは MD 法の汎用性をより増すために、従来使用していた MD プログラムの改良および機能拡張を行った。
2. 計算方計算方法
法
MD プログラムは、我々のグループで開発され論者が改良および機能拡張したものを使用した。
ポテンシャル関数は以下に示す Born-Mayer-Huggins 型の 2 体ポテンシャルに Anderson10)により報告されている 3 体ポテンシャルを付加したものを用いた。

f ij =
2 体項

3 体項

Zi Z je2
4pe 0+ Bij (rij )æ ZZj ö÷ exp{(s i + s j - rij )r ij }Bij (rij ) = Aij ç1 + i +ç nn j ÷øiè{ (f ijk = (j ij + j ik ) exp - g ijk q - q ijk)}
2j ij = -C ij Bij (z ij ) exp{(z ij - rij )h ij }ここで、Zi は i 番目の原子の電荷、e2/4pe0 は定数(230.7×10-30Jm)、rij は原子 i, j 間の距離、Aijは反発定数、ni は原子 i の最外殻電子数、si は原子 i に関するサイズパラメーター、ri は原子i に関するソフトネスパラメーター、qijk は原子 i, j, k 間の結合角、Cij, zij, hij, gijk, q は任意定数である。
3 配位ホウ素[B(3)]、4 配位ホウ素[B(4)]を安定化させるために、3 体項は O-B(3)-O、O-B(4)-O、B(3)-O-B(3)、B(3)-O-B(4)のユニットに対して 4 つ適用した。これにより、2 体間の直線的な相互作用だけでなく共有結合性による結合の方向性も決めることができる。また 3 体項を適用するにあたって、結合判定を行う打切距離 rcut と適用配位数の制限を設けており、これにより系の密度および B と O の配位数を制御することができるようにした。
B2O3 ガラス構造モデルを得るまでの本研究のフローを Fig.2 に示す。本研究で取り扱うポテンシャルパラメーターは 30〜40 個と多く、従来の試行錯誤法では最適なパラメーターを得ることは困難である。したがって、結晶構造を維持するポテンシャルはガラスにも適用できると考え、まず結晶にてポテンシャルパラメーターの最適化を行った。ボレート結晶において 300K,0Pa で 2500step 計算を行い、計算後の結晶構造（原子間距離、密度）が初期結晶構造を維持できるように非線形最小二乗法を用いてポテンシャルパラメーターを更新させ、このループを繰り返すことにより結晶構造を維持するポテンシャルパラメーターを得た。従来は各パラメーターを線形独立にするために電荷など幾つかの重要なパラメーターを固定する必要があったが、このパラメーターの導出法は全てのパラメーターを最適化の対象とすることができるため、非常に汎用性の高い手法である
Fig.2 Experimental flow of this study.と言える。
2結晶において最適化を行い得られたポテンシャルパラメーターを用いて B2O3 ガラスの構造シミュレーションを行った。粒子数は合計 500 個であるとし、初期座標は乱数で与えた。
まず、温度 6000K にて NVT アンサンブルを用いて乱数構造の初期緩和を 6000step 行った。
その後 6000K から実際の融体状態である 1300K まで 100K ずつ系の温度を降下させ、1300K
NPT アンサンブル (P=0Pa)から 300K までは 50K ずつ降下させた。
各温度ステップで 5000step、で計算を行った。300K にて NPT アンサンブルで 20000step 構造緩和を行い、B2O3 ガラス構造モデルを得た。
動径分布関数(radial distribution function : RDF)は pair-function 法 11)により導出した。
Pij (l ) =()1 Qmax a i a jexp - bQ 2 cos QldQ
22 ò0gexp(-bQ2)
a は原子の散乱能、g は sharpenning factor、l はピーク位置 rij からの距離(l=r-rij)、
ここで、は減衰因子、Q は波数ベクトルである。a の原子散乱能には X 線回折では原子散乱因子 f 、中性子線回折では散乱長 b が用いられる。ｇは中性子線回折の場合 g=1、X 線回折では一電子あたりの平均散乱因子 fe=Σfi/ΣZi を用いる。RDF の縦軸に 2p2rr(r)を与えた場合、距離 rijに位置する Nij 個の i-j 原子対は NijPij(l)/rij のピークとして表せる。
振動スペクトルは Brawer12)の示した自己相関関数をフーリエ変換することによって得た。
I (w ) = ò-¥ C (t ) × C (0) exp(- iwt )dt¥

Raman : C (t ) = r&i (t )å
jrij (t )Infrared : C (t ) = m (t ) = qi ri (t )rij (t )ここで、wは角振動数、t は時間、r は位置ベクトル(rij=rj-ri)、mは双極子モーメントである。
Raman スペクトルについては、本研究では B を中心とした BOn (n=1〜4)ユニットにおける振動モードと、O を中心とした OBm (m=1〜2)の振動モードの両者を導出し評価を行った。
これらの評価は 2000step、NVT アンサンブルで統計計算を行い導出した。

frequency / -3. 計算結果
計算結果
3-1. ポテンシャルパラメーターの最適化結果
本研究では幾つかのボレート結晶を用いてポテンシャルパラメーターの最適化を行っているが、発表では Li2B4O7 結晶、B2O3-I 結晶を用いた結果を示す。Li2B4O7 結晶は BO3 三角形と BO4 四面体ユニットが 1:1 で存在するダイボレート構造から構成される結晶であり、B2O3I 結晶は BO3 ユニットが鎖状構造を形成している結晶でボロキソール環は存在しない。
Li2B4O7 結晶にてポテンシャルパラメーターの最適化を行ったが、当初 BO3 ユニットと BO4ユニットにおける B-O 結合距離の違いを再現すcal. without charge transfer
cal. with charge transfer
ることができなかった。そこで、配位数と結合
B-O pairLi2B4O 7 crystal
距離に応じて形式電荷を変化させるという考え方を新たに導入し、BO3 と BO4 ユニットにおける B-O 間の電荷移動量に差を与えた。電荷移動量は、以下に示す式で与えた。
Del電荷移動量 =exp{b(r - X )} + 1
0.100.120.140.160.180.20

r : 結合距離,b, Del , X : 任意のパラメーター
3r / nmFig.3 Pair correlation function of B-O pair.

22p rrr (r) / x10 barns nm
2p-222p rr
r (r) / x10 barns nm
2p-2

その結果、B(3):+1.171, B(4):+1.069, O:-0.927 のcal.B-B
neutroncal.B-O10
ように電荷を与えることができ、Fig.3 のよう
cal.B-Li
cal.O-O
8cal.Li-Oに B-O 結合距離の分布をよく維持することがで
cal.Li-Li
cal.total6
きた。この電荷移動も考慮して最適化を行った
cryst.total
4
結 果 、 Fig.4 の中 性 子線 RDF に 示 すよ う に
22500step 計算後もよく構造を維持するポテンシ0ャルが得られた。また振動スペクトルの導出も0.00.1
0.2
0.3
0.40.5行い、振動特性も良く再現できたためこの得ら
r / nmれたポテンシャルを B2O3 ガラスの構造シミュ
Fig.4 Neutron RDFs of Li2B4O7 calculated
レーションに適用した。しかし、Fig.4 中のfrom the initial and resultant structure.
r=0.42nm 付近の再現性については良好とは言い
12cal.B-B
cal.B-O
cal.O-O
cal.total
neutron
難い。これは Li2B4O7 結晶は正方晶(c/a=1.08)でcryst.total
10あるのに対し、この計算で用いた MD プログラ
8ムの仕様上、立方体セルで計算を行ったためで6ある。そこで直方体セルも取り扱える
4MXDORTO9)プログラムに 2 体および 3 体ポテ2ンシャルを導入し、さらに非線形最小二乗プロ0グラムを導入することにより、直方体セルも扱
0.00.1
0.2
0.30.4r / nmえ、従来のプログラムと同様に計算を行えるプ
Fig.5 Neutron RDFs of B2O3-I calculated
from the initial and resultant structure.
ログラム(MXD)を作成した。残念ながら、Li2B4O7結晶については MXD プログラムを用いてポテンシャルパラメーターを最適化させることはできなかったが、B2O3-I 結晶では Fig.5 に示すように良く構造を維持できるポテンシャルパラメーターを得ることができた。この B2O3-I 結晶から得られたポテンシャルも B2O3 ガラスの構造シミュレーションに適用した。

Number of boron / -3-2. B2O3 ガラスの構造シミュレーション
120
Li2B4O7 結晶から得られたポテンシャルを用いたシミュレーションによるガラス構造モデルを
100 B in boroxol rings: 19.5%
MD1 モデルと呼ぶ。モデル中の B の O 配位数80はこれまで知られてきたようにほぼ全て 3 配位60であった。また密度は実測値が 1.82×103kg·m-3 で40あるのに対し、1.78×103kg·m-3 と非常に良く再現していた。本研究にて問題の一つとするボロキ20ソール環の有無について調べるために MD1 モデ
005101520ル中のリングサイズの分布について調べたとこ
Ring size / ろ[Fig.6]、ボロキソール環の存在が認められた。
Fig.6 Ring size distribution in MD1 model.
ボロキソール環中に存在するホウ素原子は全てのホウ素原子の 19.5%であった。Fig.7 に MD1 モデルの中性子線 RDF を示す。0.25〜0.34nmの距離域が B-Opair により明確に再現されているが、0.36nm 付近の特徴的なピークを再現することができなかった。また、振動スペクトルを導出したところ IR スペクトルは実測スペクトルと非常に良い対応を示したが、Raman スペクトルについては実測スペクトルの特徴的な 800cm-1 付近に現れるスペクトルを再現できなかった。しかし、Fig.8 に示すようにボロキ4ソール環を構成している BO3 ユニットのみ抽出すると、800cm-1 付近のスペクトルを再現することができた。
MD1. B-B
MD1. O-O
obs. Ref. [2]neutron

MD1. B-O
MD1. total

Raman

MD1. B-O(BO 3 in B3O 6 units)
obs. Ref. [4]Intensity / arb.units

302520151022p rrr (r) / x10 barns nm
2p-2

35

500.0

0.1

0.2

r / nm

0.3

0.4

400

0.5

Fig.7 Observed and calculated(MD1) neutron
RDFs of B2O3 glass.

800
1200-1wave number / cm

1600

Fig.8 Observed and calculated(MD1) Raman
spectra of B2O3 glass.

22p rrr (r) / x10 barns nm
2p-222p rrr (r) / x10 barns nm
2p-2

B2O3-I 結晶から得られたポテンシャルを用い
MD2.B-B
MD2.B-O30
neutron
MD2.O-O
MD2.total
たシミュレーションにはプログラム MXD を用obs. Ref. [2]25
い、得られたガラス構造モデルを MD2 モデル20と呼ぶ。MD2 モデルは前述の結果と異なり、ボ15
ロキソール環の局所的な生成は認められなかっ10た。Fig.9 に示すように 0.25〜0.34nm の距離域
5を B-Opair は再現できなかった。Johnson ら 2)は0中性子線回折から、Fig.10 に示すように 0.25〜
0.0
0.1
0.2
0.3
0.40.5r /nm
0.36nm 付近のピークはボロキソール環の生成を
Fig.9 Observed and calculated(MD2) neutron示すものであると提案しており、本研究にてボRDFs of B2O3 glass.
ロキソール環の得られた MD1 モデルはその提案に対して非常によく対応しているのに対し、Johnson's boroxol model. Ref.[2]30MD2 モデルはボロキソール環が存在しないた
B-O pairMD1 model
MD2 model
25
め Johnson らのモデルとは一致しなかった。
これらの結果から、従来示唆されてきたよ20(c)うに B2O3 ガラス中にはボロキソール環が存在(b)(a)
15-1
し、Raman スペクトルにおける 800cm 付近の10(c)特徴的なピークはボロキソール環の存在による(b)5(a)ものであることを支持する結果が得られた。過0去の研究 1)-3)によると B2O3 ガラス中には 40〜
0.1
0.2
0.3
0.40.5r / nm
90%のボロキソール環中のホウ素原子が存在し、Fig.10 Neutron RDFs of B-O pair.
3 つの BO3 ユニットが同一平面内に存在することによって RDF における 0.36nm 付近のピークが現れるとされている。本研究によって得られた B2O3 ガラス構造モデルはボロキソール環が少なく、それらにおける BO3 ユニットが同一平面内にないため Raman スペクトルおよび動径分布関数をうまく再現することができなかったものと思われ、実際の B2O3 ガラスではより高分率の平面ボロキソール環の存在が予想される。

5

4.総4.総括○結晶を用いてポテンシャルパラメーターの最適化を行い、B2O3 ガラスの構造シミュレーションを行った結果、従来困難であったボロキソール環の再現に成功した。しかしながら実測データを十分満足するまでには至っておらず、実際の B2O3 ガラス中にはより高分率な平面ボロキソール環が存在することが示唆された。
○本研究で開発した分子動力学計算プログラムは、既知の結晶構造から非線形最小二乗法によりポテンシャルパラメーターの自動最適化を行うことが可能であり、結晶およびガラス材料へのアプローチにおいて非常に有効な手法であると言える。
【参考文献】参考文献】
1) R. L. Mozzi and B. E. Warren, J. Appl. Cryst., 3 (1970) 251.
2) P. A. V. Johnson and A. C. Wright, J. Non-Cryst. Solids, 50 (1982) 281.
3) R. E. Youngman, S. T. Haubrich, J. W. Zwanziger, M. T. Janicke, B. F. Chmelka, SCIENCE, 269(1995) 1416.
4) C. F. Windisch, Jr. and W. M. Risen, Jr., J. Non-Cryst. Solids, 48 (1982) 307.
5) B. N. Meera and J. Ramakrishna, J. Non-Cryst. Solids, 159 (1993) 1.
6) E. I. Kamitsos, A. P. Patsis, M. A. Karakassides and G. D. Chryssikos, J. Non-Cryst. Solids, 126(1990) 52.
7) H. Inoue, N. Aoki and I. Yasui, J. Am. Ceram. Soc., 70 (1987) 622.
8) A. Takada, C. R. A. Catlow and G. D. Price, J. Phys. : Condens. Matter., 7 (1995) 8693.
9) 平尾一之、河村雄行、「パソコンによる材料設計」、裳華房 (1994).

10)D. C. Anderson and J. Kieffer, J. Chem. Phys., 98 (1993) 8978.
11)B. E. Warren, X-RAY DIFFRACTION, Dover Publications, Inc., (1990)
12) S. Brawer, J. Chem. Phys., 79 (1983) 4539.

6