8
hp150121 「京」産業利用（実証利用） K Industrial Use
分子動力学法を用いた溶融ガラス中の添加元素の拡散シミュレーション
Simulation of diffusion of additive elements in molten glass using molecular
dynamics method
大塚 順 Jun Otsuka 住友電気工業株式会社
SUMITOMO ELECTRIC INDUSTRIES, LTD.
要旨 溶融ガラス中に微量添加した元素の拡散係数を予測するため, 古典分子動力学シミュレーショ ンにおいて「京」を利用した大規模計算を実施した. 計算規模を拡大することで統計量が向上し, 添加濃度が低い元素でも平均二乗変位が時間に対して単調な増加を示すようになり, 結果, 拡散 係数を求めることが可能になることを実証した. Na2O-SiO2系ガラスに関して Na の添加濃度と拡 散係数の関係についても調べた. ガラスの組成を Na2O : α SiO2と表現したとき, α = 3 ~ 10,000 の 組成を計算して, α = 10,000 についても拡散係数を把握するのに十分な統計量を確保できることを 実証した. ただし, 拡散係数は実測値と乖離しており, 古典分子動力学のポテンシャルなどを見 直すことが今後の課題である. キーワード：工業用ガラス, 分子動力学シミュレーション, 拡散解析, 微量添加, 大規模並列 Abstract
In order to predict the diffusion coefficient of elements added in a small amount in molten glass, a large-scale calculation using K computer was performed in the classical molecular dynamics simulation. By increasing the scale of calculation, the statistics became sufficient, and even for elements with low additive concentration, the mean square displacement showed a monotonic increase over time, and as a result, it was proved that the diffusion coefficient could be obtained. We also investigated the relationship between the Na concentration and the diffusion coefficient for Na2O-SiO2 glasses. When the glass composition was expressed as Na2O : α SiO2, the composition of α = 3 to 10,000 was calculated, and it was demonstrated that even for α = 10,000, sufficient statistics could be guaranteed to understand the diffusion coefficient. However, the diffusion coefficient deviates from the measured value, and it is a future task to improve the potential of classical molecular dynamics.
© 2020 Research Organization for Information Science and Technology All rights reserved. Received: 29 October 2019
Accepted: 24 March 2020 Available online: 30 March 2020
9
Keywords：Industrial Glass, Molecular Dynamics Simulation, Diffusion Analysis, trace-amount doping, Massive Parallel 1. 研究の背景と目的 工業用ガラスでは, SiO2や B2O3などを主成分とし, 必要とする特性を発現させるため何らかの 元素が添加される. 当然ながら, その特性に応じて添加する元素の種類や濃度は異なる. 例えば, 光学用途ではガラスの屈折率を上げるため Pb や Ba などが添加される. また, ガラスの融点や熱 膨張率の制御には Na や Ti などが添加元素として用いられる. ガラスの製造は溶融状態を経て成 形することが多い. この場合, 添加元素の拡散という問題が生じる. 即ち, 目指す性能を実現す るため, 添加元素の拡散を制御し, 適切な元素分布を実現する必要がある. 拡散係数は実験的に 求めることもできるが, 多様な濃度や複数の添加元素が共存する状態での拡散挙動を全て実験で 求めることは非常に困難である. そこで我々は分子動力学シミュレーションを用いた拡散挙動の 解析技術を開発し, コンピュータシミュレーションによる材料設計技術の開発に取り組んでいる. 分子動力学シミュレーションで得られる原子座標の平均二乗変位は時間に対して線形に比例し て増加する関係となり, この傾きから拡散係数が求められる. ただし, 添加元素の濃度が 0.1%に も満たない微量な場合では, 平均二乗変位を計算するのに十分な統計量を確保するには計算モデ ル中の原子数が膨大となり, 通常のクラスターマシンでは能力不足である. そこで本研究では, スーパーコンピュータを利用して大規模並列計算を実行することで, 微量 添加された元素の拡散を高精度に予測できることの実証を目的とした. 2. 計算モデル 古典分子動力学法では, 原子間ポテンシャルに基づいて相互作用を計算し, 運動方程式を数値 積分することで個々の原子の運動をシミュレーションする. 本研究では, Born-Mayer 型ポテンシ ャル： 𝛷𝛷 =4π𝜖𝜖e2 0 ZiZj 𝑟𝑟 + B exp �− 𝑟𝑟 𝜌𝜌� を用いた[1]. 第一項はクーロン力, 第二項は反発力である. ここで, e は電気素量, ε0は誘電率, r は 原子間の距離, Ziは各原子の電荷である. 類似組成の結晶や分子での安定距離や実測データとの 一致などを考慮して, 各パラメータを決定した. クーロン力の計算には長距離相互作用を精度良 く高速に計算できる PPPM (Particle-Particle Particle-Mesh)法を用いた. PPPM 法ではポアソン方程 式の計算に FFT (Fast Fourier Transform)を用いるため, 全 PM (Particle-Mesh)格子上の空間電荷が必 要となる. したがって, 複数ノードを使用して並列計算する場合には, 各ノードで計算された空 間電荷をお互いに送受信する必要がある. FFT 計算には FFTW[2]を使用した.
10 原子数は 3 万~400 万とした. 初期構造の作成手順は, ランダムに原子位置を決めた原子構造を 4,000 K で溶かした後, 所定の温度まで冷却させた. NVT アンサンブルによる計算を行った. 時間 ステップは 1 fs とした. 周期境界条件を課した. 3. 並列計算の方法と効果（性能） 分子動力学計算は LAMMPS[3]を使用して行った. 使用した計算機は「京」である. 並列計算で 使用したノード数は 3, 6, 12, 24, 48, 96, 192, 384, 768, 1536 である. 並列計算は各ノードで 8 コア使 用するフラット MPI で実行した. 時間ステップを 100 ステップ計算したときの実行時間を図 1 に 示す. 原子数に対してノード数が少ないときはノード数に反比例して実行時間が短くなったが, ノード数が増やしていくと実行時間が増加に転じた. これはノード数の増加に伴って通信処理に 要する時間が長くなることが原因である. 図から実行時間が最小になるノード数は 40 万原子で 1,000 ノード程度, 400 万原子で 5,000 ノード程度と予想される. 図 1 LAMMPS で 100 ステップ計算したときの実行時間. 横軸はノード数. 粒子数は 30k, 60k, 120k, 400k, 800k, 1600k, 4000k. 赤い点線は通信などに起因して増加する時間を表す. 通信処理におけるボトルネックを特定するため性能プロファイルを調査したところ, FFT が律 速となっていた. プロファイルの一例を示すと, 3 万原子を 3,072 コア使用して並列計算したとき, PPPM 法のルーチンを処理する時間が全体の約 90%を占めており, その 74%が通信時間だった. このとき, FFTW の実行効率は 1%未満, 演算性能は 32 MFLOPS だった. 最近, クーロン力の計算
11
を効率化する方法として粒子格子多重極法が提唱されており[4, 5], 今後, 本テーマに対する有用 性を検証したい.
4. 研究成果 4.1. 平均二乗変位
平均二乗変位 (Mean Square Displacement, MSD)は 1 𝑁𝑁 �(𝑥𝑥𝑛𝑛(𝑡𝑡) − 𝑥𝑥𝑛𝑛(0))2 𝑁𝑁 𝑛𝑛=1 で定義される. ここで, N は全原子数, 𝑥𝑥𝑛𝑛(0) = 𝑥𝑥0は各原子の参照位置, 𝑥𝑥𝑛𝑛(𝑡𝑡)は時刻 t における各 原子の位置を表す. 本研究では, 汎用的な材料であり, 実験とシミュレーションの報告例もいく つかある Na2O-SiO2系ガラスについて, Na2O の添加濃度を変更したときの拡散シミュレーション を実施した. 計算した濃度は, ガラスの組成を Na2O + α SiO2と表現したとき, α = 3 ~ 10,000 の範 囲とした. 温度は 2,000 K とした. 各濃度における計算モデル中の原子数を表 1 に示す. 表 1 計算モデル中の原子数. 図 2 にシミュレーションで得られた Na の MSD を示す. 図 2 Na の MSD の経時変化. Na2O : SiO2 = 1 : α (α = 3 ~ 10,000). 温度は 2,000 K.
12 α = 3 ~ 100 のとき, MSD は時間に対して直線的に変化を示しており, 拡散モードに到達している ことがわかる. 一方, α = 1,000 ~10,000 では, 統計的なゆらぎが現れているものの, MSD は時間に 対してほぼ比例しており, 拡散係数の予測が可能な統計性は確保できていると考えられる.結果 として, 「京」を利用したシミュレーションにより, Na2O : SiO2 = 1 : 10,000 という微量添加の組成 でも Na の MSD を充分な精度で把握できることを実証できた. 4.2. 拡散係数 3 次元の拡散係数 D は 1 6𝑡𝑡〈(𝑥𝑥 − 𝑥𝑥0)2〉 と定義される. 今回のシミュレーションで得られた D の結果を図 3 に示す. 計算モデル中の原子 数は表 1 にまとめてある. α が大きくなるほど, 即ち, Na 濃度が低くなるほど, Si/O/Na のいずれも D が小さくなる傾向を示した. α = 3 のとき, Si の D は 2 10-8 cm2/s となり Horbach らの報告[6]と合 致した. 一方, Na の D は 7 10-5 cm2/s となり Horbach らの報告[6]に比べて 1 桁大きかった. α = 10,000 のときの Si の D は 2 10-10 cm2/s となり, Horbach ら[7]によって報告された無添加の SiO 2に おける Si の値と比べて 2 桁大きかった. D の乖離の要因として古典分子動力学のポテンシャルの 不備が考えられ, 今後, 改善を試みたい. また, Na の拡散係数が, α = 1,000 と α = 10,000 で逆転し ており, 不自然な傾向を示した. この要因としては古典分子動力学のポテンシャルの不備や統計 量の不足など計算精度の不足が考えられ, 今後の課題である. 図 3 Na2O-SiO2系ガラスの拡散定係数 D の Na2O 添加濃度依存性. Na2O : SiO2 = 1 : α (α = 3 ~ 10,000). 温度は 2,000 K. 5. まとめと今後の課題 微量添加した元素の拡散係数を求めるため, 古典分子動力学法を用いて原子数が 400 万までの モデルを計算した. 計算には「京」を利用して 1,536 ノードまでの並列計算を実行した. 大規模並 列計算においてクーロン力の計算がボトルネックとなったことから, FFT 計算を用いない計算手
13 法の調査が今後の課題である. スーパーコンピュータを利用して計算規模を拡大できたことで微 量添加のガラスについて平均二乗変位を精度良く計算できるようになることを実証した. Na2 O-SiO2系ガラスの拡散係数を添加濃度の関数として評価したところ, 先行研究と整合しない点が確 認された. 添加元素の濃度範囲を広くカバーできる古典分子動力学法のポテンシャルを構築する ことが今後の課題である. 将来的には古典分子動力学法を, 所望の特性を発現する材料を設計す るためのツールとして活用したい. 参考文献
[1] Y. Saito, J. Iihara, K. Yamaguchi, Y. Haruna, and M. Onishi, SEI TECHNICAL REVIEW, 67 (2008) 27-32.
[2] M. Frigo, and S. G. Johnson, Proceedings of the IEEE 93 (2), (2005) 216–231. [3] S. Plimpton, J Comp Phys, 117 (1995) 1-19.
[4] Y. Ohno, HPCI Research Report, hp150111 (2019) (Advance Online Publication). [5] K. Nitadori, arXiv: 1409.5931v2 [astro-ph.IM].
[6] J. Horbach, W. Kob, and K. Binder, Chemical Geology, 174 (2001) 87-101. [7] J. Horbach, and W. Kob, PHYSICAL REVIEW B, 60 (1999) 3169-3181.
