Simulation of borosilicate glasses with non-constant force field molecular dynamics 4 D O I: 1 0. 15 82 6/ ch im te ch .2 01 9. 6. 1. 01 Raskovalov A. A. Chimica Techno Acta. 2019. Vol. 6, No. 1. P. 4–11. ISSN 2409–5613 A. A. Raskovalov The Institute of High-Temperature Electrochemistry, Ural Branch of the Russian Academy of Sciences, 20 Akademicheskaya St., Yekaterinburg, 620137, Russian Federation e-mail: other@e1.ru Simulation of borosilicate glasses with non-constant force field molecular dynamics In this study the simulation of microscopical behavior of borosilicate glasses was conducted with non-constant force field molecular dynamics. The suggested model consists of classical pair potentials in the Buckingham form, long range Coulomb interaction, intramolecular bonded interactions and possibility of bond breaking and formation. The latter effects are accompanied by changes in the types of the bond-forming particles. The simulated system corresponds to the structure of borosilicate glasses with predominantly four-coordinated boron atoms. Different structure groups are formed due to the dissociation / formation of intramolecular bonds, and the processes of the glass network rearrangement intensifies with temperature increasing. Keywords: molecular dynamics; non-constant force field; glass; borate; silicate. Received: 11.12.2018. Accepted: 28.12.2018. Published: 29.03.2019. © Raskovalov A. A., 2019 Introduction Borosilicate glasses are materials with a wide range of applications. They can be used for nuclear waste utilization [1], as implants for tissue repair [2], in elec- trochemical devices as  ionic conductors [3–5], as sealants for solid oxide fuel cells [6] and in  many other areas of  human activity. From this point of  view, micro- scopical behavior of these glasses comes to researchers’ attention. Glasses consisting of  several network former oxides (B2O3, SiO2, P2O5 and some others) can contain a  colossal amount of  possible structure groups and their combinations. For this reason numerical simulations of  mixed glasses are very promoting because it can “view” glass structure at molecular level. To date, a molecular dynamics (MD) simula- tion of silicate, borate and mixed borosili- cate glasses was conducted in series of pa- pers [1–5, 7–19]. Most of these works deal with classical MD, except researches [2, 17] performed with Car-Parinello (MD) [20] and reax force field [21] (ReaxFF), correspondingly. The Car-Parinello meth- od brings quantum mechanical calcula- tion into MD simulation. The ReaxFF is constucted to provide chemical reactivity during MD simulation. The last two me- thods are used to include effects of chemi- cal bonds dissociation and recombination to the simulation. These effects are needed to  be considered for  correct simulation of glass behavior. To perform such calcu- lations with classical MD is an ambiguous and complicated, if not impossible, task. 5 In the present paper I demonstrate how one can perform MD simulation of boro- silicate glass with non-constant force field methodology. This approximation consists of possibility of bond breaking and forma- tion (including changing in particle types of bond forming atoms) runtime. Simulation details MD simulations were performed with an original program written in C in canon- ical (NVT) ensemble. Newton’s equations of motion were integrated by Velocity Ver- let algorithm [22] with timestep of 0.5 fs during 100’000 steps. Equilibration time was 0.5 ps (1000 timesteps). Electrostatic interactions were accounted using the Ewald summation. Nosé-Hoover thermo- stat [23] with relaxation time of  0.02 ps was used for maintaining the temperature around desired temperature. There was tested three temperatures: 298, 500 and 1000 K.  The considered system consists of 94 Si atoms, 80 B atoms, 429 O atoms and 242 alkaline metal (Me) atoms (the to- tal number of particles is 845). The oxygen atoms were divided into bridging (Oc) and non-bridging (Ot) species. Bridging oxygen connects with two atoms of glass-forming element (Si / B) and non-bridging — with only one. Initial configuration of the sys- tem was generated with own special code. The box was cubic with the edge length of 20.4 Å to match an approximate density of glassy silicates. Van der Waals interac- tions in the system were given by Bucking- ham pair potential: U A r C rij ij ij ij ij ij � � � � �� � � �� �exp , � 6 (1) where Uij is the potential energy between the i-th and the j-th atoms, rij is the dis- tance between them, and Aij, ρij and Cij are empirical parameters. Besides the Van der Waals short-range interactions, covalent bonds between oxygen and silicon / boron were set in the form of an intramolecular harmonic potential: U k r rij ij ij ij� �� �2 0 2 , (2) where kij is the spring constant and r0ij is the equilibrium distance between parti- cles. Unlike interactions provided by Eq. (1), this potential is applied only to  co- valent bonded pairs (S – Ot, Si – Oc, B – Ot and B – Oc) listed separately. In addition, a  three-body potential energy term was applied to maintain a feasible valent an- gle for Si – Oc – Si, B – Oc – B and Si – Oc – B bonds: U k ijk tb ijk� �� �2 0 2 cos cos ,� � (3) where ktb is the spring constant of the three- body potential, 100.0 eV, θijk is the angle be- tween ij and ik bonds, θ0 is the equilibrium angle, 148.3° (taken from [15] for Si – O – Si angle); i, j and k are the indexes of Oc atom and its covalent-bonded neighbors, cor- respondingly. The values of Aij, ρij and Cij were taken from papers [1, 11, 12, 15, 19]. Simulations were performed with different combina- tion of these parameters. The more suitable for the glass structure description set of the parameters is summarized in Table 1 with corresponding references. The parameters of the valent bond potential (2) are given in Table 2. The atomic charges were sug- gested to be found as 4δ, 3δ, a, — 2δ and -δ-a for Si, B, Me, Oc and Ot, correspond- ingly. The values of δ and a are set to 0.4e and 1.0e. Our MD program allows deleting 6 (adding) valent bonds from (to) the corre- sponding list to mimic processes of bonds dissociation and glass forming: … ≡ Si – Oc – R ↔ ≡Si+ Ot – R (4) and … ≡ B – Oc – R ↔ –B + Ot – R. (5) If the bonds automatically break (form) then the interatomic distance is higher (lower) than a maximal bond distance pa- rameter (rm). The values of this parameter are given in Table 2. The further dissocia- tion with participation of Si+ or B+ species is not allowed. The charges of Si+ and B+ species are chosen to keep electroneutrality during reactions (4) and (5). Table 2 Valent bonds parameters (including the maximal bond distance, rm) used in this study pair k, eV r0, Å rm, Å Si — Ot 10 –4 1.625 — Si — Oc 10 –4 1.650 2.0 B — Ot 10 –5 1.400 — B — Oc 10 –5 1.430 1.9 Results and Discussion The obtained structure of  the system is presented in Fig. 1. One can see that the simulated glass consists of branchy chains of Si — O / B — O bonds which can form various combinations and loops. In prin- ciple, there are two possible boron coor- dinations in borate glasses, triangular and tetrahedral [24]. In our case almost all bo- ron in the glass is four-coordinated, Fig. 1. Radial distribution functions (RDFs) are presented in  Fig.  2. All RDF-curves consist of  one sharp and several broad maximums. This is typical of  liquids or amorphous systems and indicates the presence of  the short-range and the ab- sence of the long-range orders. There are no series of well-resolved maxima, specific for crystalline solids. The first maximum of  RDF-curve corresponds to  the most probable distance in the near-neighbor co- ordination. According to Fig 2(a), the dis- tance between Me and Ot is less than that between Me and Oc. It can be explained by  more negative charge of  Ot species (Simulation Details section) and, there- fore, stronger Coulomb interaction. RDFs almost do not depend on temperature, the differences are observed only for the RDFs Table 1 The Buckingham pair potential parameters used in this study pair A, eV ρ, Å C, eVÅ6 Reference Si – Ot / Oc 13702.905 0.193817 54.681 [19] B – Ot / Oc 206941.81 0.124 35.0018 [19] Me – Ot / Oc 4383.7555 0.243838 30.700 [5] Ot / Oc – Ot / Oc 352.56 0.35 0 [1] Si – B 337.584 0.29 0 [1] Si – Me 861.744 0.29 0 [1] Si – Si 836.16 0.29 0 [1] B – B 121.056 0.29 0 [1] Me – Me 9500.0 0.23 0 [12] 7 with the sharpest maxima, for  example, Si – O. For these RDFs the value of the first maximum decreases with temperature in- creasing and the peak becomes broader, Fig. 2(b). This process occurs, most likely, due to the thermal motion. The first ma- ximum for B – Ot / Oc and Si – Ot / Oc pairs locates at  distance of  1.375–1.425 and 1.575–1.675 Å, correspondingly, which is close to B – O и Si – O interatomic dis- tances in borosilicate glasses [15]. As the system initially does not have crystalline lattice and a  corresponding look of RDF-curves, it is necessary to find another criterion of  transition between liquid and solid states. Such a  criterion can be, for  example, mean square dis- placement (MSD), which characterizes Fig. 2. Some radial distribution functions, g(r) of the simulated glass: (a) for 298 K; (b) for pair Si – Ot at 298 and 1000 K Fig. 1. The structure of the simulated system after 30’000 timesteps at 298 K in (a) stick-and-ball and (b) polyhedral representation. Blue polyhedra correspond to boron central atom, gray — to silicon 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 Me – Ot Me – Si Ot – Oc g( r) r, Å Me – Oc 1 2 3 4 5 6 0 2 4 6 8 10 12 14 16 18 20 22 298 K g( r) r, Å 1000 K a b 8 a  mobility of  ions. In the simulated bo- rosilicate glass, MSDs of  most ions tend to  reach a  plateau with some oscillation at temperature of 298 K, except one for Me ions increased linearly, Fig. 3(a). The time derivative of  MSD grows with tempera- ture increasing, Fig. 3(b). This means that all species start to diffuse and the system behavior becomes liquid-like. Analogous effects can be seen on trajectories of the ion motion. At temperature of 298 K most ions vibrate around one point (although with big enough amplitude), Fig. 4(a). An exception is the Me ion, its trajectory con- sists of a series of hops between positions with long enough oscillation time. With temperature increasing, Me ion spends less time in one position, more often hops and covers greater distance, Fig. 4(b). Its move- ment character becomes more and more liquid-like. Fig. 5 demonstrates the number of Si and B species as a function of time. The changes in these quantities are provided by dissociation of B – Oc and Si – Oc bonds according to equations (4) and (5). Since dissociation in  the suggested model oc- curs upon reaching the determined dis- tance (Table 2), the observed fluctuations of the particles numbers are related with Fig. 3. The mean square displacement (MSD) of the simulated borosilicate glass: (a) for different ions at 298 K; (b) for Me and Si ions at different temperatures Fig. 4. Trajectories of motion in the simulated borosilicate glass for: (a) different ions at 298 K; (b) the same Me ion at different temperatures 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 Ot Oc B Si M SD , Å 2 time, ps Me 5 10 15 20 25 30 35 40 45 50 0 10 20 30 40 50 60 Si (500 K) Me (298 K) Me (500 K) Me (1000 K) Si (1000 K) M SD , Å 2 time, ps Si (298 K) a b 7 8 9 10 11 12 13 14 15 8 9 10 11 12 13 BOc O c Y , Å X, Å Me 10 11 12 13 14 15 16 17 18 19 20 21 4 5 6 7 8 9 10 11 12 500 K 298 K Y , Å X, Å 1000 K a b 9 changing in  the B / Si – Oc distances. Ac- cording to Fig. 5, the suggested pair po- tentials provide stronger B — Oc bonding, then Si – Oc one, as the amplitude of the Si number oscillation is bigger. Also, Fig. 5 demonstrates that the amplitudes of these oscillations increase with temperature. This is a result of thermal motion which pro- motes overcoming of forces of interatomic attraction, as noted above in the discussion of  the radial distribution functions. The bonds breaking and formation realized by methodology of non-constant force field allow glass structure to be “dynamical,” i.e. to destroy ones atomic groups and to cre- ate others. A variety of formed molecular groups can see in Fig. 1. Conclusions In the paper reference classic pair po- tentials have been tested for molecular dy- namic simulation of borosilicate glasses. The suggested combination of the poten- tials and ionic charges describes interato- mic distances in borosilicate glasses rea- sonably well. Oxygen atoms were divided into bridging and non-bridging by  the value of electrical charge. The transition between these states is possible due to bond formation / dissociation. These phenomena are implemented by  non-constant force field molecular dynamics. It has been shown that the degree of the bond disso- ciation increases with temperature growth. Acknowledgments The reported study was funded by Russian Foundation for Basic Research (RFBR), according to the research project No. 16-33-60095 mol_а_dk. References 1. Delaye JM, Ghaleb D. Molecular dynamics simulation of a nuclear waste glass matrix. Materials Science and Engineering B. 1996;37:232–6. 2. Tilocca A. Sodium migration pathways in multicomponent silicate glasses: Car — Parrinello molecular dynamics simulations. J Chem Phys. 2010;133:014701. DOI:10.1063 / 1.3456712. 3. Verhoef AH, Den Hartog HW. Molecular Dynamics simulations of borate glasses. Rad Effect Defect Solids. 1991;119–121(2):493–8. DOI:10.1080 / 10420159108220770. 4. Varsamis C-PE, Vegiri A, Kamitsos EI. Molecular dynamics investigation of lithium borate glasses: Local structure and ion dynamics. Phys Rev B. 2002;65:104203. DOI:10.1103 / PhysRevB.65.104203. Fig. 5. The number of Si and B species as a function of simulation time 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 B (1000 K) Si (1000 K) B (298 K) N time, ps Si (298 K) 10 5. Cormack AN, Du J, Zeitler TR. Alkali ion migration mechanisms in silicate glasses probed by molecular dynamics simulations. Phys Chem Chem Phys. 2002;4:3193–7. DOI:10.1039 / b201721k. 6. Smeacetto F, Miranda A, Chrysanthou A, et al. Novel Glass-Ceramic Composition as Sealant for SOFCs. J Am Ceram Soc. 2014;97:3835−42. DOI:10.1111 / jace.13219. 7. Inoue H, Aoki N, Yasui I. Molecular Dynamics Simulation of the Structure of Borate Glasses. J Am Ceram Soc. 1987; 70(9):622–7. DOI:10.1111 / j.1151–2916.1987.tb05729.x. 8. Xu Q, Kawamura K, Yokokawa  T.  Molecular dynamics calculations for  bo- ron oxide and sodium borate glasses. J Non-Cryst Solids. 1988;104(2-3):261–72. DOI:10.1016 / 0022–3093(88)90397–3. 9. Vessal B, Amini M, Leslie M, Catlow CRA. Potentials for Molecular Dynamics Simula- tion of Silicate Glasses. Mol Sim. 1990;5(1-2):1–7. DOI:10.1080 / 08927029008022407. 10. Cormack AN, Cao Yu. Molecular Dynamics Simulation of Silicate Glasses. Mol Eng. 1996;6:183. DOI:10.1007 / BF00161727. 11. Rossano S, Ramos A, Delaye J-M, Creux S, Filipponi A, Brouder Ch, Calas G. EXAFS and Molecular Dynamics combined study of CaO-FeO-2SiO2 glass. New insight into site significance in silicate glasses. Europhys Lett. 2000;49(5):597–602. 12. Cormak AN, Du J. Molecular dynamics simulation of soda-lime-silicate glasses. J Non-Cryst Solids. 2001;293–295:283–9. 13. Gou F, Greaves GN, Smith W, Winter R. Molecular dynamics simulation of sodium borosilicate glasses. J Non-Cryst Solids. 2001;293(1):539–46. DOI:10.1016 / S0022-3093(01)00775-X. 14. Sawaguchi N, Yamaguchi K, Sasaki M, Kawamura  K.  Interatomic Potential Model for  Molecular Dynamics Simulation of  Lithium Borate Melts / Glasses. J Comp Chem. 2015;14(4):139–46. DOI:10.2477 / jccj.2015–0017. 15. Scherer C. Molecular dynamics simulations of silicate and borate glasses and melts: Structure, diffusion dynamics and vibrational properties [dissertation]. Mainz (Ger- many): Fachbereich Physik, Mathematik und Informatik der Johannes Gutenberg- Universitat; 2015. 189 p. 16. Li X, Song W, Yang K, Anoop Krishnan NM, Wang B, Smedskjaer MM, Mauro JC, Sant G, Balonis M, Bauchy M. Cooling rate effects in sodium silicate glasses: Bridging the gap between molecular dynamics simulations and experiments. J Chem Phys. 2017; 147:074501. DOI:10.1063 / 1.4998611. 17. Yu Y, Wang B, Wang M, Sant G, Bauchy M. Reactive molecular dynamics simulations of sodium silicate glasses — Toward an improved understanding of the structure. Int J Appl Glass Sci. 2017;8:276–84. DOI:10.1111 / ijag.12248. 18. Stevensson B, Yu Ya, Ed´en M. Structure-Composition Trends in Multicomponent Borosilicate-Based Glasses Deduced from Molecular Dynamics Simulations with Improved B — O and P — O Force Fields. Phys Chem Chem Phys. 2018;20:8192–209. DOI:10.1039 / C7CP08593A. 19. Wang M, Anoop Krishnana NM, Wang B, Smedskjaer MM, Mauro JC, Bauchya M. A new transferable interatomic potential for  molecular dynamics simulations 11 of borosilicate glasses. J Non-Cryst Solids. 2018;498:294–304. DOI:10.1016 / j.jnoncrysol.2018.04.063. 20. Car R, Parrinello M. Unified Approach for Molecular Dynamics and Density-Func- tional Theory. Phys Rev Lett. 1985;55:2471–4. DOI:10.1103 / PhysRevLett.55.2471. 21. van Duin ACT, Dasgupta S, Lorant F, Goddard WA. Reax FF: A Reactive Force Field for Hydrocarbons. J Phys Chem A. 2001;105:9396–409. DOI:10.1021 / jp004368u. 22. Verlet  L.  Computer “Experiments” on  Classical Fluids. I.  Thermodynamical Properties of Lennard-Jones Molecules. Phys Rev. 1967;159(1):98–103. DOI:10.1103 / PhysRev.159.98. 23. Martyna GJ, Tuckerman ME, Tobias DJ, Klein ML. Explicit Reversible Integrators for Extended Systems Dynamics. Mol Phys. 1996;87(5):1117–57. DOI:10.1080 / 00268979600100761. 24. Wu J, Stebbins JF. Cation Field Strength Effects on Boron Coordination in Binary Borate Glasses. J Amer Ceram Soc. 2014; 97(9):2794–801. DOI:10.1111 / jace.13100.