Structural study of Na2O–B2O3–SiO2 glasses from molecular simulations using a polarizable force field

Sodium borosilicate glasses Na 2 O-B 2 O 3 -SiO 2 (NBS) are complex systems from a structural point of view. Three main building units are present: tetrahedral SiO 4 and BO 4 (B IV ), and triangular BO 3 (B III ). One of the salient features of these compounds is the change of the B III /B IV ratio with the alkali concentration, which is very diﬃcult to capture in force ﬁelds-based molecular dynamics simulations. In this work we develop a polarizable force ﬁeld able to reproduce the boron coordination and more generally the structure of several NBS systems in the glass and in the melt. The parameters of the potential are ﬁtted from density functional theory calculations only, in contrast with the existing empirical potentials for NBS systems. This ensures a strong improvement on the transferability of the parameters from one composition to another. Using this new force ﬁeld, structure of NBS systems is validated against neutron diﬀraction and nuclear magnetic resonance (NMR) experiments. A special focus is given to the distribution of B III /B IV with respect to the composition and the temperature.

Structural study of Na 2 O-B 2 O 3 -SiO 2 glasses from molecular simulations using a polarizable force field Three main building units are present: tetrahedral SiO 4 and BO 4 (B IV ), and triangular BO 3 (B III ). One of the salient features of these compounds is the change of the B III /B IV ratio with the alkali concentration, which is very difficult to capture in force fields-based molecular dynamics simulations. In this work we develop a polarizable force field able to reproduce the boron coordination and more generally the structure of several NBS systems in the glass and in the melt. The parameters of the potential are fitted from density functional theory calculations only, in contrast with the existing empirical potentials for NBS systems. This ensures a strong improvement on the transferability of the parameters from one composition to another. Using this new force field, structure of NBS systems is validated against neutron diffraction and nuclear magnetic resonance (NMR) experiments. A special focus is given to the distribution of B III /B IV with respect to the composition and the temperature. Sodium borosilicate glasses (Na 2 O-B 2 O 3 -SiO 2 , NBS) have a widespread interest in number of fields such as bioactive glasses or the storage of nuclear waste 1 . They are chemically inert, mechanically strong and resistant to thermal shocks; these properties explain their widespread usage in the glass industry. The structural study of sodium borosilicate glasses by diffraction techniques is rather complex due to the presence of several components. In contrast NMR has shown to be a powerful probe of the glass structure since the seminal work of Dell and Bray 2-4 , but it also has its limits in the local view of the glass structure. It is therefore of interest to complement such experimental studies by classical molecular dynamics simulations. However, the modeling of these systems is rather complex 5 . Indeed, NBS glasses comprise two network formers: silicon with a tetrahedral structure (SiO 4 ) and boron with a tetrahedral or a triangular structure (BO 4 and BO 3 respectively, noted B IV and B III afterwards) (Figure 1). These two kinds of boron coordinations are an important characteristic of the NBS ternary systems, which display a non-linear structural change with respect to the composition, known as the boron anomaly. The amount of alkaline atoms has a direct impact on the boron coordination number. On the one hand, they play the role of charge compensator for the negatively charged BO 4 units. On the other hand, at high alkali oxide contents, they can also act as network modifiers when they are located close to a BO 3 unit or a SiO 4 unit, leading to the creation of non-bridging oxygens (NBOs). Dell and Bray developed a model based on two pa- to rationalize the variations of boron coordination with the composition 2 . This model was parameterized on data extracted from NMR experiments. Basically, at low R, each sodium ion is a charge compensator and the amount of B IV increases linearly with R until R = 0.5 + K/16. In this range reedmergnerite units, which consist of a BO 4 surrounded by four SiO 4 , are formed. Some recent studies using Raman spectroscopy 6-8 introduced the possibility of danburite units to occur (BO 4 surrounded by three SiO 4 and one BO 4 ) at this stage. Then, as R increases, the B IV concentration reaches a plateau and NBOs start to be created. Finally, if the Na + ion concentration is further increased, the amount of four-fold coordinated boron slowly drops as BO 3 units bearing NBO form. Since then, several improvements to this model were proposed, such as the one of Smedskjaer et al. which is based on topological principles 9 , yielding results in good agreement with the available experimental data.
Further studies based on neutron and X-ray diffraction have shown that temperature strongly impacts on the boron coordination [10][11][12][13] . They observed that the amount of B III largely increases when the glass is heated. This means that when the glass is formed, tetrahedral boron units are created during quenching. This specific feature of sodium borosilicate glasses has a direct effect on dynamic properties like the viscosity and the electrical conductivity 14 .
The main goal of this work is to follow these structural changes along the quenching for several NBS compositions using molecular dynamics. Various interaction potential sets have been developed in order to explore silicate and borosilicate systems by classical molecular dynamics [15][16][17][18][19][20][21] . They traditionally use the Born-Mayer-Huggins (BMH) analytical form and manage to reproduce the structure, the boron coordination at room temperature with respect to the composition and the mechanical properties, in good agreement with the experiments. In most cases they were fitted from experimental data. The main weakness of such interaction potentials is there low transferability. Indeed, some of their parameters (the partial charges in the Coulomb term and the parameters of the short-range repulsion term for B-O interactions) depend on the composition, which means that independent fits have to be made for each studied NBS system 22 . In this work, we propose a different approach. The force field includes a single set of parameters, which is used for a wide range of NBS compositions. Our model uses fixed charges which are set to the formal valence charges. In order to account for many-body polarization effects, we use the Polarizable Ion Model (PIM) developed by Madden et al. [23][24][25][26][27] . This term is an effective representation of the deformation of the electronic density in response to electric fields. As a consequence, the PIM allows a higher flexibility for the local structure than conventional non-polarizable force fields. The parameters are fitted to a large body of reference quantities (individual forces and dipoles) determined for several representative atomic configurations using Density Functional Theory (DFT) calculations. Similar interaction potentials have already been developed and used for many systems. They were shown to reproduce accurately the structural and dynamic properties in many situations for which multivalent ions play a key role [28][29][30][31][32][33] . Here the structure of the NBS systems is validated against neutron diffraction and NMR experiments. This allows us to shed some light on the evolution of the B III /B IV ratio with the sodium ion concentration and the temperature. We also study the variations of density and the swelling with temperature.

A. Polarizable ion model
The PIM force field is composed of four terms : chargecharge, repulsion, dispersion and polarization. The three first terms consist in interactions between pairs of atoms where f ij 6 and f ij 8 are Tang-Toennies damping functions 34 describing the short-range correction to the dispersion. They are defined as follow: The polarization term is: where α i is the dipole polarizability of ion i and µ i its induced dipole moment. The charge-dipole asymptotic terms are also corrected by Tang-Toennies damping functions for penetration effects at short-range: This interaction is many-body in essence since the induced dipoles are additional degrees of freedom which are calculated at each time step (or for each individual configuration) by minimization of V polarization .

B. Fitting procedure
In order to fit the parameters, we perform reference DFT calculations on a series of representative configurations, generated using empirical force fields or extracted from DFT-based MD simulations. To develop a force field able to deal with a large range of NBS compositions and the complexity of these systems, seven compositions were used for the fitting, with a total of thirty six configurations ( Table I). The parameters are fitted on the forces and the dipoles on each atom (obtained by DFT calculations) 35 by minimizing the two quantities: where the sum rules over all configurations. In a first step, the dipoles and the forces applying on each atom are calculated by DFT using the CPMD code 36 . The PBE functional was used with a plane-wave cut-off of 100 Ry 37 . Martins-Troullier's pseudo-potentials were used to account for core electrons. The forces are a direct output from the simulation. For the dipoles, it is necessary to calculate the maximally localized Wannier functions (MLWF) 38 . A complete theory of electric polarization in crystalline dielectrics has been developed [39][40][41] , which validates the calculation of the dipole moments of single ions or molecules from the center of charge of the subset of MLWF which are localized in their vicinity 42,43 . Then, iterations on the PIM parameters are performed until the forces and the dipoles calculated by classical molecular dynamics are close enough to the ab initio ones 27,35 . It is important to remind that the charges q i used in the Coulomb term are not fitted and chosen as -2, +4, +3, +1 and +1 for O, Si, B, Na and Li respectively. The fitting process only takes into account the repulsion and the polarization parameters. To analyze the structural changes and especially the boron speciation occurring during the temperature decrease, two series of simulations were performed. All these simulations used a time-step of 1 fs. The cutoff radii for the short-range interactions is set equal to half of the simulation box length. The first simulation series consisted in using five initial configurations per composition to gather more statistics on the results. Each of them contains approximately 330 atoms. First of all, the systems are equilibrated in the NPT ensemble at 3000 K for 1 ns following isobars at 0 GPa (the barostat 50 and the thermostat 51 relaxation times being respectively fixed to 5 ps and 20 ps). Then a quenching in steps of 10 K is applied to reach 300 K (still in the NPT ensemble), at a quenching rate of 1 K/ps. Once the room temperature is reached, a NVT equilibration is made for 10 ns 4 which allows to study the systems at room temperature, in the glassy state. The second series of simulations is close to the previous one except that only one initial configuration is used per composition which is composed of around 2600 atoms. The other difference is that the equilibrations at 3000 K and 300 K are reduced to 100 ps and 1 ns, respectively. Analytical grade Na 2 CO 3 (purity > 99.95%), SiO 2 (purity 99.8%) and B 2 O 3 (enrichment in 11 B > 99%) were taken as starting materials for the glass preparation. The glass powders were prepared in two steps. The material is first heated at respectively 1250 • C, 1275 • C, 1100 • C and 1250 • C for the NBS17-24, NBS14-16, NBS35-19 and NBS29-13 glasses then naturally cooled at ambient temperature. The glasses are then heated again at the same temperature but during the second cooling, a quench rate of 10 K/h is applied between T g +50 K and T g -50 K (T g corresponds to the glass transition temperature). The T g are respectively equal to 543 • C, 588 • C, 467 • C and 505 • C for the NBS17-14, NBS14-16, NBS35-19 and NBS29-13 glasses. The composition of the final glass powders have been analyzed by ICP-AES (except for the NBS14-16 glass which has been analyzed by MEB-EDS). The results are given in Table IV, showing slight deviation from the targer compositions.

IV. NEUTRON DIFFRACTION
Neutron diffraction measurements have been performed on the 7C2 diffractometer at the Orphée reactor of the Laboratoire Léon Brillouin (Saclay, France), using hot neutrons (λ = 0.72Å) that give a q range extending from 0.52Å −1 to 15.2Å −1 with q = (4π/λ) sin θ, where λ is the incident wavelength and 2θ is the scattering angle. The banana-like multidetector is constituted of 256 3 He gas giving a 2D diffraction pattern 52 . NBS17-24, NBS14-16, NBS35-19 and NBS29-13 glasses have been measured as powder contained in a cylindrical vanadium cell. Each acquisition lasted 12 hours to get a good signal to noise ratio. Several corrections have been applied to the data after conversion from 2D to 1D pattern: detector efficiency, background from the container, multiple scattering attenuation and inelastic scattering. V. NUCLEAR MAGNETIC RESONANCE 11 B MAS NMR spectra were collected on a Bruker Avance II 500WB spectrometer operating at a magnetic field of 11.72 T, using a Bruker CPMAS probe (with ZrO2 rotor, outer diameter (OD) 4 mm) at a sample rotation frequency of 12.5 kHz. To ensure spectrum quantitativeness, a short pulse of 1 µs (tip angle of π/12) and a recycle delay of 2 s were used. Chemical shifts were referenced to an external sample of 1 M boric acid solution (19.6 ppm relative to the boron trifluoride etherate). Data processing occurred via an in-house code (for details see references 4,53) to extract the relative population of B IV units.

VI. RESULTS AND DISCUSSION
A. Structural properties

Structure factors
A first step to evaluate the ability of the potential to correctly reproduce the experiments consists in calculating the structure factor. The structure is then probed at both the short and intermediate range. The partial structure factors S ij are calculated from the partial radial distribution functions g ij and then used to determine the total structure factor as follows 54 : sin(qr) qr dr (7) where c i and b i are respectively the fractions and the neutron scattering lengths of the atom species i, ρ being the number density.
The neutron structure factors of the NBS systems at room temperature are plotted in Figure 2. The results for the small simulation cells (around 330 atoms) are averaged over five independent trajectories for each composition, while only one trajectory was used for the large ones. Both results are compared to the experimental data. First of all, the simulation results yield the correct positions of all the peaks. The force field is therefore able to reproduce the general structure of the NBS glasses. The experimental data displays a single first peak for NBS17-24 and NBS14-16, which is splitted in two sub-peaks for the compositions with a higher amount of sodium. These double-peaks are also present in the simulation results but with slight differences in amplitude. Despite this discrepancy, there is a noticeable improvement compared to the empirical potentials of Kieu et al. 15 which do not reproduce this splitting (only one global peak is observed). The analysis of the partial structure factors can give more information on the origin of the double-peak and on the difference in amplitude between simulations and experiments. We separate the total structure factors in two parts corresponding to the contribution of all the partial structure factors involving interactions with Na on the one hand and the contribution of all the partial structure factors without the Na ones on the other hand ( Figure 3). The second peak located around 2.15Å −1 clearly originates from the correlations involving sodium ions. This effect is more important in NBS27-21, NBS35-19 and NBS29-13 systems due to the larger amount of sodium they contain. In addition, we also observe a decrease of the intensity of the low-q peak that is due to the presence of Na.
A comparison can also be done between the simulation results, obtained from the small systems and the larger ones. The structure factor obtained in both cases are very similar (Figure 2). It means that the structural data determined from both quenches give equivalent results and that no significant size effect occurs. In the following, only the data calculated from the small systems will be used, the results with the larger systems being roughly the same.

Boron coordination
The structure can also be studied in the real space, which gives information at local range. It is important, in the NBS systems, to probe the local entities and es- pecially to analyze the variation of the B III /B IV ratio with respect to the composition and the temperature. In this aim, the change of boron coordination has been studied during the quenching and at room temperature. This property was calculated from the integral of the first peak of the B-O radial distribution function with a cut-off equal to 1.91Å (which corresponds to the first minimum of the corresponding RDF). The percentages of four-fold coordinated boron atoms are plotted in Figure 4. At room temperature (Table V), the simulation data are compared with the experimental ones, obtained from 11 B MAS NMR analyses 44 . The predicted amount of four-fold coordinated boron is in good agreement with the experiments for NBS17-24, NBS14-16, NBS27-21 and NBS29-13. One can clearly note the impact of the composition and especially the amount of sodium on the boron local structure. The force field used is able to adapt the boron coordination with the composition in the same way that the experimental measurements, except for NBS35-19. For this composition, experiments show that a large amount of sodium creates non-bridging oxygens and BO 3 units instead of compensating the BO 4 units. But in the simulations made in NBS35-19, the force field seems to still favor the tetrahedral structure of boron atoms.
When the temperature increases, the systems start to depolymerize and the triangular BO 3 entity is the most stable configuration for boron atom in every composition. In agreement with our simulations (albeit for different compositions), experiments predict a decrease of the amount of B IV upon melting the glass. A first study made by Michel et al. 12 on NBS glasses using neutron diffraction coupled with Reverse Monte Carlo revealed boron coordinations equal to 3.7 at room temperature and 3.3 at 1273 K, respectively. A similar phenomenon has also been observed experimentally by Alderman et al. in Na 2 O-B 2 O 3 56 , albeit in a different system. Measurements of the X-ray structure factor of molten Na 2 B 4 O 7 coupled with a thermodynamic model showed that this coordination number also decreased from 3.42 at room temperature to 3.17 at 1600 K. This suggests that the larger stability at B IV at low temperatures is a general feature of such glasses. As described by Alderman et al., structural rearrangement occurs at T > T g . A part of the Na compensators creates non-bridging oxygens bounded to Si and B III which imply a decrease of the amount of B IV . Now, we focus on the local environment of the silicon and boron atoms by analyzing the distances and the angles in the first coordination sphere at room temperature.  65 ). This could explain the difference we observe between our theoretical and the experimental bond lengths, the potential introduced in this work being adjusted from DFT calculations.  It is also important to analyze the angles, to assess our structural models. They are given in Table VII. The angles within the tetrahedral and triangle structural units involving boron atoms are in good agreement with experiments 57,61,62 . The silicon tetrahedra have angles slightly lower but close to the expected value [58][59][60] . The angles centered on oxygen atoms linking the silicate and borate units are also close to the experimental results 59,60,63 . These distances and angles analyses show that the force field developed is able to reproduce the local structure around boron and silicon atoms in good agreement with the experiments.

Distances and angles between first neighbors
As discussed previously, the sodium ions can play two roles in sodium borosilicate sytems. They can act as charge compensators around B IV units or as network modifiers. They are separated using the following procedure. Firstly, all the NBOs are identified. An oxygen atom is considered as non-bridging if there is no more than one network-former (Si or B) in its first neighbor shell. To detect the Si and B first neighbors, distances cut-off are determined from the first minima of the corresponding radial distribution function, i.e. 2.18Å for Si-O pairs and 1.91Å for B-O pairs. In a second step, the closest Na atom around each NBO is found and labelled as a modifier. All the other Na atoms are then charge compensators. As for the B III /B IV ratio, the relative amounts of modifiers and compensators strongly change with the composition of the system (Table VIII). The percentage of charge modifiers decreases with the R ratio decreasing as expected from the Dell and Bray model. When the number of Na ions exceeds the B IV concentration, all of them cannot be charge compensator and non-bridging oxygens are formed around BO 3 or SiO 4 units. The distances between oxygen and sodium ions can also be analyzed by distinguishing their modifier and charge compensator nature. As shown in Table VIII, the network modifiers are closer to the oxygen atoms than the compensators. One can also remark that these distances seem to depend on the compensator/modifier ratio. The difference between them becomes smaller when the amount of modifier increases.
The relative populations of NBOs and bridging oxygens (BOs) also change during the quenching, as shown on Figure 5. Indeed, the amount of NBOs arises from a complex interplay between the amount of Na + , the boron coordination number, as well as the change of role of sodium ions (modifier to charge compensator). They are preferentially located in the vicinity of Si and of B III (they are equally shared between them), while a few of them are incorporated into B IV units.

B. Density and swelling
The procedure used in this work for the quenching (in the NPT ensemble) allows the volume to fluctuate. It is thus possible to follow the density variations with respect to the temperature, as shown on Figure 6. The experimental values 44 at room temperature are also shown. We observe that the density is underestimated for the five compositions, by values ranging between 2 and 7 %.  Several origins can be proposed for this discrepancy. Of course, there are intrinsic errors linked to the force field, but the fast quenching rate also has a strong impact. As shown in Figure 7, the underestimation of the density increases with the SiO 2 concentration. Since the latter is a strong glass former, it is likely increasing substantially the viscosity of the melt. This result therefore points to the necessity of performing longer quenches for the melts with higher SiO 2 concentrations. During the quenching, the density progressively increases and the density variation is governed by both the thermal motion and the structural changes. In order to separate these two contributions, we performed the following analysis 44,66,67 . Firstly, a set of ten configurations has been considered at ambient temperature for each composition, taken from the final relaxation at 300 K. These configurations are called the reference configurations hereafter. They are characterized by the numbers of local entities they contain. Such local entities are defined according to the list given in Table IX. These 17 entities are called the reference entities.
Once the population of the local entities contained in a reference configuration have been listed, their local volume is estimated. To do this, the individual Voronoï volumes of all the atoms are calculated. Then the local volume of a FO n−m type entity, i.e. a former F (Si or B) surrounded by a total of n O atoms including m NBOs, is measured as: Then, the total volume of a reference configuration is given by the summation of the volumes of all its local entities given by the equation 9 for FO n−m type, and calculated directly for the Na modifiers and charge compensators. The average volumes (V avg,i ) of each reference entity i are then estimated on the basis of the ten reference configurations.
For each temperature and thus at different density, a configuration is characterized by its volume V conf . Knowing the amount of each entity present in this configuration, it is straightforward to estimate the part of the volume due to the structural entities (V struct ). Indeed, if this configuration contains C i local entities of type T i , the volume associated to the structure is calculated as: The remaining part may then be associated to the thermal motion: This model allows to determine, for any configuration, the relative parts of the volume variation associated to the structural changes and to the thermal motion. The total swelling (S tot ), structural swelling (S struct ) and thermal swelling (S therm ) are then defined as follow : where V ref is the mean volume of the reference configurations. Figure 8 displays the swelling of each NBS system versus the temperature and the contribution due to the structural changes. The latter does not exceed 2.5 %, which means that whatever the temperature, the thermal motion is the main origin of the swelling and that structural changes are a secondary effect.
The main structural contribution to the density changes seems to arise from the NBO formation when the temperature increases. Figure 9 displays the changes of the volumes of all entities BO 3 (entities 1 to 4), BO 4 (entities 5 to 9), BO 3 +BO 4 (entities 1 to 9) and SiO 4 (entities 10 to 14) with respect to the temperature. One can clearly note the decrease of the amount of BO 4 units for the benefit of the BO 3 units. However, the volumes of the boron reference units are roughly the same ( the temperature increases and not the boron structural changes. The NBO creation also explain the increase of the SiO 4 volume units. The structural swelling is therefore mainly impacted by the NBO creation due to the structural changes between BO 3 and BO 4 and the depolymerization of the network at high temperature. On the right panel of Figure 8, one can clearly remark the difference between systems with ratios K = 2.5 or K = 4.5. For temperatures higher than the glass transition, the slope of the structural swelling (Table X) is more important for systems having a larger amount of B 2 O 3 (K = 2.5). This can be explained by a larger number of conversions B IV to B III in these systems, which involves the creation of a larger amount of NBO.
NBS35-19 seems to have a particular behavior. Its structural modifications have more impact on the swelling than the other systems at low temperature. It can be explained by its large amount of NBOs, even at room temperature. The plasticity associated to the depolymerisation in this glass explains why reorganizations can occur at a temperature lower than T g . This special behavior can be linked to a recent study of the stress corrosion cracking 68 , where it was observed that NBS35-19 has a special behavior in contact with water. Indeed, in this system, it is easier for water to penetrate into the glassy network which would imply a larger depolymerization than the other systems.

VII. CONCLUSIONS
This study introduces a polarizable force field for the Na 2 O-B 2 O 3 -SiO 2 systems, with parameters fitted to DFT reference data and without any use of empirical information. This force field offers a good transferability, since only one set of parameters is used for all the compositions. Usable for a large range of NBS compositions, it yields structural properties of NBS glasses that agree well with experimental data. New insights into the impact of the composition on the local structure has been obtained though MD simulations. A continuous conversion of BO 3 into BO 4 is observed during the quench (at a rate depending on the composition), i.e. the boron speciation widely changes from B III to B IV when the temperature decreases. Moreover, the impact of sodium ions on the boron coordination at room temperature is in good agreement with experimental observations. For higher R ratios, NBO amount increases and are mainly linked equally to Si and B III . The thermal motion and the structural changes contributions have been separated to explain the swelling. Even if the thermal motion remains the main origin of the swelling, the network depolymerization plays an increasing role when the temperature