<figure><img src="image_1.png" /></figure>

HAL
open science

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

Fabien Pacaud, Jean-Marc Delaye, Thibault Charpentier, Laurent Cormier, Mathieu Salanne

► To cite this version:

Fabien Pacaud, Jean-Marc Delaye, Thibault Charpentier, Laurent Cormier, Mathieu Salanne. Structural study of Na2O-B2O3-SiO2 glasses from molecular simulations using a polarizable force field. The Journal of Chemical Physics, 2017, 147, pp.161711. 10.1063/1.4992799 . cea-01570528

HAL Id: cea-01570528

https://cea.hal.science/cea-01570528

Submitted on 1 Aug 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Structural study of Na₂O-B₂O₃-SiO₂ glasses from molecular simulations using a polarizable force field

Fabien Pacaud,¹,² Jean-Marc Delaye,² Thibault Charpentier,³ Laurent Cormier,⁴ and Mathieu Salanne¹

¹) Sorbonne Universités, UPMC Univ. Paris 06, CNRS, Laboratoire PHENIX, F-75005 Paris, France

²) CEA, DEN, Laboratoire d'Etude des Matériaux et Procédés Actifs, 30207 Bagnols-sur-Cèze, France

³) NIMBE, CEA, CNRS, Université Paris-Saclay, CEA Saclay, 91191 Gif-sur-Yvette, France

⁴) Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie, Sorbonne Universités - Université Pierre et Marie Curie Paris 6, CNRS UMR 7590, IRD UMR 206, Muséum National d'Histoire Naturelle, F-75005 Paris, France

(Dated: 22 June 2017)

Sodium borosilicate glasses Na₂O-B₂O₃-SiO₂ (NBS) are complex systems from a structural point of view. Three main building units are present: tetrahedral SiO₄ and BO₄ (B<sup>IV</sup>), and triangular BO₃ (B<sup>III</sup>). One of the salient features of these compounds is the change of the B<sup>III</sup>/B<sup>IV</sup> 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<sup>III</sup>/B<sup>IV</sup> with respect to the composition and the temperature.

PACS numbers: 61.43.Bn, 65.60.+a, 81.05.Kf

## I. INTRODUCTION

Sodium borosilicate glasses (Na₂O-B₂O₃-SiO₂, NBS) have a widespread interest in number of fields such as bioactive glasses or the storage of nuclear waste¹. 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²⁻⁴, 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⁵. Indeed, NBS glasses comprise two network formers: silicon with a tetrahedral structure (SiO₄) and boron with a tetrahedral or a triangular structure (BO₄ and BO₃ respectively, noted B<sup>IV</sup> and B<sup>III</sup> 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₄ units. On the other hand, at high alkali oxide contents, they can also act as network

managers when they are located close to a BO₃ unit or a SiO₄ unit, leading to the creation of non-bridging oxygens (NBOs).

<figure><img src="image_3.png" /><figcaption>FIG. 1. Main structural units which coexist in NBS glasses and melts: the tetrahedral units of SiO₄ (a) and BO₄ (b), and the triangular unit of BO₃ (c). Silicon atom is in yellow, boron atoms are in pink and oxygen ones are in red.</figcaption></figure>

Dell and Bray developed a model based on two parameters R = [Na₂O]/[B₂O₃] and K = [SiO₂]/[B₂O₃] to rationalize the variations of boron coordination with the composition². 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<sup>IV</sup> increases linearly with R until R = 0.5 + K/16. In this range reedmergnerite units, which consist of a BO₄ surrounded by four SiO₄, are formed. Some recent studies using Raman spectroscopy⁶⁻⁸ introduced the possibility of danburite units to occur (BO₄ surrounded by three SiO₄ and one BO₄) at this stage. Then, as R increases, the B<sup>IV</sup> 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₃ 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⁹, 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¹⁰⁻¹³. They observed that the amount of B<sup>III</sup> 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¹⁴.

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¹⁵⁻²¹. 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²². 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.²³⁻²⁷. 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²⁸⁻³³. 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<sup>III</sup>/B<sup>IV</sup> ratio with the sodium ion concentration and the temperature. We also study the variations of density and the swelling with temperature.

## II. COMPUTATIONAL METHODS

### A. Polarizable ion model

The PIM force field is composed of four terms: charge, repulsion, dispersion and polarization. The three first terms consist in interactions between pairs of atoms

$$V_{\text{pair}} = \sum_{i,j>i} \left[ \frac{q_i q_j}{r_{ij}} + A^{ij} \exp(-a^{ij} r_{ij}) - f_6^{ij}(r_{ij}) \frac{C_6^{ij}}{r_{ij}^6} - f_8^{ij}(r_{ij}) \frac{C_8^{ij}}{r_{ij}^8} \right] \quad (1)$$

where $f_6^{ij}$ and $f_8^{ij}$ are Tang-Toennies damping functions³⁴ describing the short-range correction to the dispersion. They are defined as follows:

$$f_n^{ij}(r_{ij}) = 1 - \exp(-b_n^{ij} r_{ij}) \sum_{k=0}^{n} \frac{(b_n^{ij} r_{ij})^k}{k!} \quad (2)$$

The polarization term is:

$$V_{\text{polarization}} = \sum_{i} \sum_{j>i} \left[ \frac{q_i \mathbf{r}_{ij} \cdot \boldsymbol{\mu}_j}{r_{ij}^3} g^{ij}(r_{ij}) - \frac{\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij} q_j}{r_{ij}^3} g^{ji}(r_{ij}) \right] \\ + \sum_{i} \sum_{j>i} \left[ \frac{\boldsymbol{\mu}_i \cdot \boldsymbol{\mu}_j}{r_{ij}^3} - \frac{3(\mathbf{r}_{ij} \cdot \boldsymbol{\mu}_i)(\mathbf{r}_{ij} \cdot \boldsymbol{\mu}_j)}{r_{ij}^5} \right] \\ + \sum_{i} \frac{|\boldsymbol{\mu}_i|^2}{2\alpha_i} \quad (3)$$

where $\alpha_i$ is the dipole polarizability of ion i and $\boldsymbol{\mu}_i$ its induced dipole moment. The charge-dipole asymptotic terms are also corrected by Tang-Toennies damping functions for penetration effects at short-range:

$$g^{ij}(r_{ij}) = 1 - c^{ij} \exp(-b^{ij} r_{ij}) \sum_{k=0}^{4} \frac{(b^{ij} r_{ij})^k}{k!} \quad (4)$$

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_{\text{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
<table><caption>TABLE I. Compositions in molar percent of the configurations used for the fitting of the force field. The number of independent configurations for each composition is specified in the second column.</caption><thead><tr><th>System</th><th>Configurations</th><th>%Na<sub>2</sub>O</th><th>%B<sub>2</sub>O<sub>3</sub></th><th>%SiO<sub>2</sub></th><th>%Li<sub>2</sub>O</th></tr></thead><tbody><tr><td>NBS14-18</td><td>7</td><td>14</td><td>18</td><td>68</td><td>-</td></tr><tr><td>NBS15-15</td><td>2</td><td>15.3</td><td>15.3</td><td>69.4</td><td>-</td></tr><tr><td>NBS19-13</td><td>2</td><td>19.1</td><td>12.6</td><td>68.3</td><td>-</td></tr><tr><td>NBS30-15</td><td>3</td><td>30</td><td>15</td><td>55</td><td>-</td></tr><tr><td>SiO<sub>2</sub>-B<sub>2</sub>O<sub>3</sub></td><td>20</td><td>-</td><td>16.7</td><td>83.3</td><td>-</td></tr><tr><td>Li<sub>2</sub>O-B<sub>2</sub>O<sub>3</sub>-1</td><td>1</td><td>-</td><td>90</td><td>-</td><td>10</td></tr><tr><td>Li<sub>2</sub>O-B<sub>2</sub>O<sub>3</sub>-2</td><td>1</td><td>-</td><td>70</td><td>-</td><td>30</td></tr></tbody></table>

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)³⁵ by minimizing the two quantities:

$$ \chi_{\mu}^{2} = \frac{1}{N} \sum_{i=1}^{N} \frac{|\mu_{i}^{\text{DFT}} - \mu_{i}^{\text{PIM}}|^{2}}{|\mu_{i}^{\text{DFT}}|^{2}} \quad (5) $$

$$ \chi_{F}^{2} = \frac{1}{N} \sum_{i=1}^{N} \frac{|\mathbf{F}_{i}^{\text{DFT}} - \mathbf{F}_{i}^{\text{PIM}}|^{2}}{|\mathbf{F}_{i}^{\text{DFT}}|^{2}} \quad (6) $$

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³⁶. The PBE functional was used with a plane-wave cut-off of 100 Ry³⁷. 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)³⁸. A complete theory of electric polarization in crystalline dielectrics has been developed³⁹⁻⁴¹, 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⁴²,⁴³.

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²⁷,³⁵. 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. Table II reports the final parameters which yield values of 0.067 and 0.373 for $\chi_{\mu}^2$ and $\chi_F^2$ respectively. The dispersion parameters were taken from the literature²⁷.

### C. Simulation details

Several experimental results have been reported in the literature for NBS glasses²,¹⁴,⁴⁴⁻⁴⁹. Five NBS composi-

<table><caption>TABLE II. Force field parameters obtained from the fitting process (repulsion and polarization) and from the literature (dispersion)²⁷. All values are in atomic units. Cations interact together only via their strongly repulsive Coulomb interaction. The parameters involving lithium ions are not given since they have not been tested in this work. The polarisabilities of oxide and sodium ions are respectively 9.868 and 0.775 atomic units. The range of the Tang-Toennies correction to the charge-dipole interaction is taken the same within an ion pair, i.e. $b^{ij} = b^{ji}$.</caption><thead><tr><th>Ion pair</th><th>O-O</th><th>Si-O</th><th>B-O</th><th>Na-O</th></tr></thead><tbody><tr><td rowspan="2">repulsion</td><td>A<sup>ij</sup></td><td>201.66</td><td>49.811</td><td>31.209</td><td>201.08</td></tr><tr><td>a<sup>ij</sup></td><td>2.2121</td><td>1.7188</td><td>1.8466</td><td>2.2190</td></tr><tr><td rowspan="3">dispersion</td><td>C<sub>6</sub><sup>ij</sup></td><td>22.186</td><td>2.1793</td><td>-</td><td>2.1793</td></tr><tr><td>C<sub>8</sub><sup>ij</sup></td><td>426.65</td><td>25.305</td><td>-</td><td>25.305</td></tr><tr><td>b<sub>6</sub><sup>ij</sup></td><td>1.4000</td><td>2.2057</td><td>-</td><td>2.2057</td></tr><tr><td rowspan="3">polarization</td><td>b<sub>8</sub><sup>ij</sup></td><td>1.4000</td><td>2.2057</td><td>-</td><td>2.2057</td></tr><tr><td>b<sup>ij</sup></td><td>2.975</td><td>1.9515</td><td>2.0564</td><td>2.1371</td></tr><tr><td>c<sup>ij</sup></td><td>3.3595</td><td>1.4492</td><td>1.2172</td><td>3.3340</td></tr><tr><td></td><td>c<sup>ji</sup></td><td>-</td><td>-</td><td>-</td><td>2.8571</td></tr></tbody></table>

tions (Table III), which have already been studied experimentally, are chosen in this work to compare the simulation data to their experimental counterparts. Several ratios R = [Na₂O]/[B₂O₃] and two ratios K = [SiO₂]/[B₂O₃] are considered in order to estimate the robustness and the transferability of the force field and to study the impact of the composition on the structural properties.

<table><caption>TABLE III. Na₂O-B₂O₃-SiO₂ studied compositions in molar percent.</caption><thead><tr><th>System</th><th>%Na<sub>2</sub>O</th><th>%B<sub>2</sub>O<sub>3</sub></th><th>%SiO<sub>2</sub></th><th>R</th><th>K</th></tr></thead><tbody><tr><td>NBS17-24</td><td>16.5</td><td>23.9</td><td>59.6</td><td>0.69</td><td>2.49</td></tr><tr><td>NBS14-16</td><td>14.4</td><td>15.5</td><td>70.1</td><td>0.93</td><td>4.52</td></tr><tr><td>NBS27-21</td><td>26.8</td><td>20.7</td><td>52.5</td><td>1.29</td><td>2.53</td></tr><tr><td>NBS35-19</td><td>34.5</td><td>18.5</td><td>47.0</td><td>1.86</td><td>2.54</td></tr><tr><td>NBS29-13</td><td>29.0</td><td>13.0</td><td>58.0</td><td>2.23</td><td>4.46</td></tr></tbody></table>

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⁵⁰ and the thermostat⁵¹ 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
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 equilibriations at 3000 K and 300 K are reduced to 100 ps and 1 ns, respectively.

### III. GLASS PREPARATION

<table><caption>TABLE IV. Molar compositions of the experimental glasses.</caption><thead><tr><th>NBS</th><th>17-24</th><th>14-16</th><th>35-19</th><th>29-13</th></tr></thead><tbody><tr><td>%Na<sub>2</sub>O</td><td>11.5</td><td>14.6</td><td>34.4</td><td>28.7</td></tr><tr><td>%B<sub>2</sub>O<sub>3</sub></td><td>25.4</td><td>18.9</td><td>21.2</td><td>15.0</td></tr><tr><td>%SiO<sub>2</sub></td><td>63.1</td><td>65.9</td><td>44.4</td><td>56.3</td></tr></tbody></table>

Analytical grade Na₂CO₃ (purity > 99.95%), SiO₂ (purity 99.8%) and B₂O₃ (enrichment in ¹¹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 larger 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 ($\lambda = 0.72$ Å) that give a q range extending from $0.52$ Å⁻¹ to $15.2$ Å⁻¹ with $q = (4\pi/\lambda)\sin\theta$, where $\lambda$ is the incident wavelength and $2\theta$ is the scattering angle. The banana-like multidetector is constituted of 256 ³He gas giving a 2D diffraction pattern⁵². 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

¹¹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 ZrO₂ rotor, outer diameter (OD) 4 mm) at a sample rotation frequency of 12.5 kHz. To ensure spectrum quantitiveness, 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 BIV units.

### VI. RESULTS AND DISCUSSION

#### A. Structural properties

##### 1. 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⁵⁴:

$$S_{ij}(q) = 1 + \rho \int_{0}^{\infty} 4\pi r^2 [g_{ij}(r) - 1] \frac{\sin(qr)}{qr} dr \quad (7)$$

$$S(q) = \frac{\sum_{i,j} c_i c_j b_i b_j S_{ij}(q)}{\sum_{i,j} c_i c_j b_i b_j} \quad (8)$$

where $c_i$ and $b_i$ are respectively the fractions and the neutron scattering lengths of the atom species i, $\rho$ 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
<figure><img src="image_4.png" /></figure>

FIG. 2. Neutron structure factors of the NBS compositions at room temperature. The systems with around 2600 atoms (red dashed lines) and the mean value of five independent trajectories for the systems with around 330 atoms (black lines) are compared to the experimental data (black circles). No experimental information is available for NBS27-21. A shift in ordinate is made for a better visualization.

et al.¹⁵ 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 Å⁻¹ 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.

## 2. Boron coordination

The structure can also be studied in the real space, which gives information at local range. It is important,

<figure><img src="image_5.png" /></figure>

FIG. 3. The total neutron structure factors (black curve) are split into two contributions: one arising from the partial structure factors involving Na atoms (red) and one involving all the others interactions (blue). The whole set of total neutron structure factors is also given in the right bottom panel.

<figure><img src="image_6.png" /></figure>

FIG. 4. Amount of four-fold coordinated boron atoms during the quenching for different NBS compositions. The stars are the experimental results at room temperature⁴⁴,⁵⁵.

in the NBS systems, to probe the local entities and es-
<table><caption>TABLE V. Amount of %B<sup>IV</sup> in NBS systems calculated at room temperature and compared to the experimental data. The errors are given by the standard deviation of the five independent trajectories for each composition.</caption><tr><td>NBS</td><td>17-24</td><td>14-16</td><td>27-21</td><td>35-19</td><td>29-13</td></tr><tr><td>%B<sup>IV</sup> (simulation)</td><td>57±4</td><td>62±7</td><td>77±3</td><td>78±3</td><td>78±7</td></tr><tr><td>%B<sup>IV</sup> (experiment)</td><td>60</td><td>72</td><td>71</td><td>62</td><td>78</td></tr></table>

especially to analyze the variation of the B<sup>III</sup>/B<sup>IV</sup> 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 <sup>11</sup>B MAS NMR analyses<sup>44</sup>. 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<sub>3</sub> units instead of compensating the BO<sub>4</sub> 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<sub>3</sub> 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<sup>IV</sup> upon melting the glass. A first study made by Michel *et al.*<sup>12</sup> 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<sub>2</sub>O-B<sub>2</sub>O<sub>3</sub><sup>56</sup>, albeit in a different system. Measurements of the X-ray structure factor of molten Na<sub>2</sub>B<sub>4</sub>O<sub>7</sub> 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<sup>IV</sup> at low temperatures is a general feature of such glasses. As described by Alderman *et al.*, structural rearrangement occurs at T > T<sub>g</sub>. A part of the Na compensators creates non-bridging oxygens bounded to Si and B<sup>III</sup> which imply a decrease of the amount of B<sup>IV</sup>.

### 3. Distances and angles between first neighbors

<table><caption>TABLE VI. Average distances between Si-O, B<sup>III</sup>-O and B<sup>IV</sup>-O in the first coordination sphere at 300 K, from simulations and experiments<sup>57,58</sup>.</caption><thead><tr><th>NBS</th><th>17-24</th><th>14-16</th><th>27-21</th><th>35-19</th><th>29-13</th><th>exp.</th></tr></thead><tbody><tr><td>d(Si-O) / Å</td><td>1.63</td><td>1.63</td><td>1.63</td><td>1.63</td><td>1.63</td><td>1.61</td></tr><tr><td>d(B<sup>III</sup>-O) / Å</td><td>1.38</td><td>1.38</td><td>1.38</td><td>1.39</td><td>1.38</td><td>1.37</td></tr><tr><td>d(B<sup>IV</sup>-O) / Å</td><td>1.48</td><td>1.48</td><td>1.48</td><td>1.49</td><td>1.49</td><td>1.48</td></tr></tbody></table>

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. Table VI compares the average distances between Si-O, B<sup>III</sup>-O and B<sup>IV</sup>-O calculated from simulation with experimental values, which were obtained for various crystalline and glassy systems (HBO<sub>2</sub> crystal, vitreous silica, etc)<sup>57-63</sup>. Si-O and B<sup>III</sup>-O bond lengths are slightly overestimated by 0.02 Å and 0.01 Å in our MD simulations, respectively. It is worth noting that DFT simulations involving PBE functional slightly overestimate this distance as well (the Si-O distance is 1.63 Å for NBS system<sup>64</sup> and 1.65 for silicate melt<sup>65</sup>). 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.

<table><caption>TABLE VII. Angles inside and between the silicate and the borate units extracted from simulations and experiments<sup>58-63</sup>.</caption><thead><tr><th>NBS</th><th>17-24</th><th>14-16</th><th>27-21</th><th>35-19</th><th>29-13</th><th>exp.</th></tr></thead><tbody><tr><td>O-Si-O / °</td><td>108.9</td><td>108.8</td><td>108.9</td><td>108.7</td><td>108.7</td><td>109.5-109.7</td></tr><tr><td>O-B<sup>III</sup>-O / °</td><td>117.3</td><td>119.2</td><td>119.4</td><td>119.6</td><td>119.5</td><td>120</td></tr><tr><td>O-B<sup>IV</sup>-O / °</td><td>109.6</td><td>109.4</td><td>109.3</td><td>109.3</td><td>109.2</td><td>109.4</td></tr><tr><td>Si-O-Si / °</td><td>143</td><td>142</td><td>144</td><td>140</td><td>142</td><td>142-144</td></tr><tr><td>Si-O-B / °</td><td>135</td><td>135</td><td>136</td><td>135</td><td>135</td><td>135.6</td></tr></tbody></table>

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<sup>57,61,62</sup>. The silicon tetrahedra have angles slightly lower but close to the expected value<sup>58-60</sup>. The angles centered on oxygen atoms linking the silicate and borate units are also close to the experimental results<sup>59,60,63</sup>. 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.

As discussed previously, the sodium ions can play two roles in sodium borosilicate systems. They can act as charge compensators around B<sup>IV</sup> 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
<table><caption>TABLE VIII. Amount of sodium labelled as charge compensators (Na$_c$) and as network modifiers (Na$_m$), and their average distance to oxygen atoms surrounding them at room temperature.</caption><thead><tr><th>NBS</th><th>17-24</th><th>14-16</th><th>27-21</th><th>35-19</th><th>29-13</th></tr></thead><tbody><tr><td>R</td><td>0.69</td><td>0.93</td><td>1.29</td><td>1.86</td><td>2.23</td></tr><tr><td>%B<sup>IV</sup></td><td>57</td><td>62</td><td>77</td><td>78</td><td>78</td></tr><tr><td>%NBO</td><td>2.0</td><td>4.3</td><td>11.7</td><td>22.3</td><td>19.8</td></tr><tr><td>%Na<sub>c</sub></td><td>87</td><td>70</td><td>58</td><td>41</td><td>37</td></tr><tr><td>%Na<sub>m</sub></td><td>13</td><td>30</td><td>42</td><td>59</td><td>63</td></tr><tr><td>d(Na<sub>c</sub>-O) / Å</td><td>2.42</td><td>2.40</td><td>2.38</td><td>2.38</td><td>2.38</td></tr><tr><td>d(Na<sub>m</sub>-O) / Å</td><td>2.22</td><td>2.22</td><td>2.26</td><td>2.26</td><td>2.26</td></tr></tbody></table>

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<sup>III</sup>/B<sup>IV</sup> 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 ratio decreasing as expected from the Dell and Bray model. When the number of Na ions exceeds the B<sup>IV</sup> concentration, all of them cannot be charge compensator and non-bridging oxygens are formed around BO<sub>3</sub> or SiO<sub>4</sub> 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<sup>+</sup>, 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<sup>III</sup> (they are equally shared between them), while a few of them are incorporated into B<sup>IV</sup> 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<sup>44</sup> at room temperature are also shown. We observe that the density is underestimated for the five compositions, by values ranging between 2 and 7 %.

<figure><img src="image_7.png" /><figcaption>FIG. 5. Number of NBOs close to a network former atom F with F = Si, B<sup>IV</sup> and B<sup>III</sup> (respectively the black, red and green lines) normalized by the number of such atoms. Bottom right panel: percentage of NBOs in all the NBS systems.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. Variation of the density during quenching. The experimental densities at 300 K are shown with stars<sup>44</sup>.</figcaption></figure>

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
<figure><img src="image_9.png" /><figcaption>FIG. 7. Difference between simulated and experimental densities with respect to the SiO₂ concentration.</figcaption></figure>

shown in Figure 7, the underestimation of the density increases with the SiO₂ 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₂ 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⁴⁴,⁶⁶,⁶⁷. 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 Voronoi volumes of all the atoms are calculated. Then the local volume of a FOₙ₋ₘ type entity, i.e. a former F (Si or B) surrounded by a total of n O atoms including m NBOs, is measured as:

$$V(FO_{n-m}) = V(F) + \frac{1}{2} \sum_{i=1}^{n-m} V(O_{BO}) + \sum_{i=1}^{m} V(O_{NBO}) \quad (9)$$

where V(F) is the volume of the former, V(Oₐ) and V(Oₐ) are respectively the volumes of the bridging oxygen and non-bridging oxygen atoms. The volume of the bridging oxygen is divided by two because these O atoms participate in two different local entities.

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ₙ₋ₘ type, and calculated directly for the Na modifiers and charge compensators. The average volumes (Vₐ,ₗ,ᵢ) 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ₐ. 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ₐ). Indeed, if this configuration contains Cᵢ local entities of type Tᵢ, the volume associated to the structure is calculated as:

$$V_{struct}(T) = \sum_{i=1}^{17} C_i(T) V_{avg,i}(T_0) \quad (10)$$

The remaining part may then be associated to the thermal motion:

$$V_{thermal-motion}(T) = V_{conf}(T) - V_{struct}(T) \quad (11)$$

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ₜ), structural swelling (Sₐ) and thermal swelling (Sₐ) are then defined as follow:

$$S_{tot}(T) = \frac{V_{conf}(T) - V_{ref}}{V_{ref}} \quad (12)$$

$$S_{struct}(T) = \frac{V_{struct}(T) - V_{ref}}{V_{ref}} \quad (13)$$

$$S_{therm}(T) = S_{tot}(T) - S_{struct}(T) \quad (14)$$

where Vₐ 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₃ (entities 1 to 4), BO₄ (entities 5 to 9), BO₃+BO₄ (entities 1 to 9) and SiO₄ (entities 10 to 14) with respect to the temperature. One can clearly note the decrease of the amount of BO₄ units for the benefit of the BO₃ units. However, the volumes of the boron reference units are roughly the same (Table IX) which means that the structural swelling involving boron units are essentially due to the NBO creation when
<table><caption>TABLE IX. The 17 reference entities used to characterize the reference configuration structures at ambient temperature and their volumes.</caption><thead><tr><th rowspan="2">Reference entity</th><th rowspan="2">N°</th><th colspan="5">Volumes of each entity AA<sup>3</sup></th></tr><tr><th>NBS17-24</th><th>NBS14-16</th><th>NBS27-21</th><th>NBS35-19</th><th>NBS29-13</th></tr></thead><tbody><tr><td>BO<sub>3</sub> without NBO</td><td>1</td><td>30.34</td><td>31.83</td><td>28.03</td><td>26.91</td><td>29.12</td></tr><tr><td>BO<sub>3</sub> with 1 NBO</td><td>2</td><td>37.91</td><td>39.91</td><td>35.39</td><td>34.28</td><td>35.04</td></tr><tr><td>BO<sub>3</sub> with 2 NBO</td><td>3</td><td>-</td><td>-</td><td>-</td><td>43.27</td><td>43.93</td></tr><tr><td>BO<sub>3</sub> with 3 NBO</td><td>4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>BO<sub>4</sub> without NBO</td><td>5</td><td>31.88</td><td>32.91</td><td>30.12</td><td>29.09</td><td>30.45</td></tr><tr><td>BO<sub>4</sub> with 1 NBO</td><td>6</td><td>38.58</td><td>45.66</td><td>36.03</td><td>35.31</td><td>36.13</td></tr><tr><td>BO<sub>4</sub> with 2 NBO</td><td>7</td><td>-</td><td>-</td><td>-</td><td>41.22</td><td>-</td></tr><tr><td>BO<sub>4</sub> with 3 NBO</td><td>8</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>BO<sub>4</sub> with 4 NBO</td><td>9</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SiO<sub>4</sub> without NBO</td><td>10</td><td>38.96</td><td>41.14</td><td>36.46</td><td>35.35</td><td>37.02</td></tr><tr><td>SiO<sub>4</sub> with 1 NBO</td><td>11</td><td>48.77</td><td>48.94</td><td>44.41</td><td>42.77</td><td>44.40</td></tr><tr><td>SiO<sub>4</sub> with 2 NBO</td><td>12</td><td>-</td><td>60.07</td><td>52.10</td><td>50.03</td><td>50.88</td></tr><tr><td>SiO<sub>4</sub> with 3 NBO</td><td>13</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SiO<sub>4</sub> with 4 NBO</td><td>14</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SiO<sub>5</sub></td><td>15</td><td>39.49</td><td>41.34</td><td>39.09</td><td>38.25</td><td>38.16</td></tr><tr><td>Na modifier</td><td>16</td><td>21.54</td><td>21.25</td><td>18.67</td><td>17.62</td><td>18.23</td></tr><tr><td>Na compensator</td><td>17</td><td>19.41</td><td>20.48</td><td>18.35</td><td>17.58</td><td>18.16</td></tr></tbody></table>

the temperature increases and not the boron structural changes. The NBO creation also explain the increase of the SiO₄ volume units. The structural swelling is therefore mainly impacted by the NBO creation due to the structural changes between BO₃ and BO₄ 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₂O₃ (K = 2.5). This can be explained by a larger number of conversions B<sup>IV</sup> to B<sup>III</sup> 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<sub>g</sub>. This special behavior can be linked to a recent study of the stress corrosion cracking<sup>68</sup>, 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.

<table><caption>TABLE X. Slopes of the structural swelling curves (extracted from Figure 8) at high temperature.</caption><thead><tr><th>NBS</th><th>17-24</th><th>14-16</th><th>27-21</th><th>35-19</th><th>29-13</th></tr></thead><tbody><tr><td>K</td><td>2.49</td><td>4.52</td><td>2.53</td><td>2.54</td><td>4.46</td></tr><tr><td>Slope</td><td>0.0020</td><td>0.0011</td><td>0.0019</td><td>0.0019</td><td>0.0012</td></tr></tbody></table>

## VII. CONCLUSIONS

This study introduces a polarizable force field for the Na₂O-B₂O₃-SiO₂ 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 through MD simulations. A continuous conversion of BO₃ into BO₄ is observed during the quench (at a rate depending on the composition), i.e. the boron speciation widely changes from B<sup>III</sup> to B<sup>IV</sup> 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<sup>III</sup>. 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
<figure><img src="image_10.png" /><figcaption>FIG. 8. Total swelling of each NBS composition (top) and the swelling due to the structural changes with respect to the temperature (bottom).</figcaption></figure>

increases, especially above $T_g$.

## ACKNOWLEDGEMENTS

We are grateful for the computing resources on OCCI-GEN (CINES, French National HPC) obtained through the projects x2015097321 and x2016087684.

## SUPPLEMENTARY DATA

A typical input file including all the force field parameters (and the corresponding output files) is provided for the CP2K code.

<figure><img src="image_11.png" /><figcaption>FIG. 9. Volumes of the all BO₃, BO₄, BO₃ + BO₄ and SiO₄ entities with respect to the temperature in each systemes, display in the left top corner, right top corner, left bottom corner and right bottom corner respectively.</figcaption></figure>

<footnote n="1">¹M. J. Plodinec, Glass Technology **41**, 186–192 (2000).</footnote>

<footnote n="2">²W. J. Dell, P. J. Bray, and S. X. Xiao, Journal of Non-Crystalline Solids **58**, 1 (1983).</footnote>

<footnote n="3">³L.-S. Du and J. F. Stebbins, The Journal of Physical Chemistry B **107**, 10063 (2003).</footnote>

<footnote n="4">⁴F. Angeli, T. Charpentier, D. de Ligny, and C. Cailleteau, Journal of the American Ceramic Society **93**, 2693 (2010).</footnote>

<footnote n="5">⁵L. Deng and J. Du, Journal of Non-Crystalline Solids **453**, 177 (2016).</footnote>

<footnote n="6">⁶B. C. Bunker, D. R. Tallant, R. J. Kirkpatrick, and G. L. Turner, Physics and Chemistry of Glasses **31**, 30 (1990).</footnote>

<footnote n="7">⁷D. Manara, A. Grandjean, and D. Neuville, American Mineralogist **94**, 777 (2009).</footnote>

<footnote n="8">⁸D. Manara, A. Grandjean, and D. Neuville, Journal of Non-Crystalline Solids **355**, 2528 (2009).</footnote>

<footnote n="9">⁹M. M. Smedskjaer, J. C. Mauro, R. E. Yougman, C. L. Hogue, M. Potuzak, and Y. Yue, J. Phys. Chem. B **115**, 12930–12946 (2011).</footnote>

<footnote n="10">¹⁰J. F. Stebbins and S. E. Ellsworth, Journal of the American Ceramic Society **79**, 2247 (1996).</footnote>

<footnote n="11">¹¹S. Sen and J. F. Stebbins, Journal of Non-Crystalline Solids **226**, 29 (1998).</footnote>

<footnote n="12">¹²F. Michel, L. Cormier, P. Lombard, B. Beuneu, L. Galoisy, and G. Calas, Journal of Non-Crystalline Solids **379**, 169 (2013).</footnote>

<footnote n="13">¹³O. L. G. Alderman, M. Liška, J. Macháček, C. J. Benmore, A. Lin, A. Tamalonis, and J. K. R. Weber, J. Phys. Chem. C **120**, 553 (2016).</footnote>

<footnote n="14">¹⁴A. Grandjean, M. Malki, C. Simmonet, D. Manara, and B. Penelon, Physical Review B **75** (2007).</footnote>

<footnote n="15">¹⁵L.-H. Kieu, J.-M. Delaye, L. Cormier, and C. Stolz, Journal of Non-Crystalline Solids **357**, 3313 (2011).</footnote>

<footnote n="16">¹⁶A. J. Connelly, K. P. Travis, R. J. Hand, N. C. Hyatt, and E. Maddrell, Journal of the American Ceramic Society **94**, 151 (2010).</footnote>

<footnote n="17">¹⁷B. Guillot and N. Sator, Geochimica et Cosmochimica Acta **71**, 1249 (2007).</footnote>

<footnote n="18">¹⁸P. Stoch and A. Stoch, Journal of Non-Crystalline Solids **411**, 106 (2015).</footnote>

<footnote n="19">¹⁹M. Bauchy, B. Guillot, M. Micoulaut, and N. Sator, Chemical Geology **346**, 47 (2013).</footnote>

<footnote n="20">²⁰F. Gou, G. Greaves, W. Smith, and R. Winter, Journal of Non-Crystalline Solids **293-295**, 539 (2001).</footnote>

<footnote n="21">²¹J.-M. Delaye and D. Ghaleb, Journal of Non-Crystalline Solids **195**, 239 (1996).</footnote>
²²L.-H. Kieu, J.-M. Delaye, and C. Stolz, Journal of Non-Crystalline Solids **358**, 3268 (2012).

²³P. A. Madden and M. Wilson, Chem. Soc. Rev. **25**, 339 (1996).

²⁴A. Aguado and P. A. Madden, The Journal of Chemical Physics **118**, 5718 (2003).

²⁵S. Jahn, P. A. Madden, and M. Wilson, Phys. Rev. B **74**, 024112 (2006).

²⁶M. Salanne and P. A. Madden, Molecular Physics **109**, 2299 (2011).

²⁷M. Salanne, B. Rotenberg, S. Jahn, R. Vuilleumier, C. Simon, and P. A. Madden, Theor Chem Acc **131** (2012).

²⁸P. A. Madden, R. Heaton, A. Aguado, and S. Jahn, Journal of Molecular Structure: THEOCHEM **771**, 9 (2006).

²⁹J.-P. Piquemal, L. Perera, G. A. Cisneros, P. Ren, L. G. Pedersen, and T. A. Darden, J. Chem. Phys. **125**, 054511 (2006).

³⁰P. D'Angelo and R. Spezia, Chem. Eur. J **18**, 11162 (2012).

³¹A. Marjolin, C. Gourlaouen, C. Clavaguera, P. Y. Ren, J.-P. Piquemal, and J.-P. Dognon, J. Mol. Mod. **20**, 2471 (2014).

³²R. Spezia, Y. Jeanvoine, C. Beuchat, L. Gagliardi, and R. Vuilleumier, Phys. Chem. Chem. Phys. **16**, 5824 (2014).

³³F. Martelli, Y. Jeanvoine, T. Vercouter, C. Beuchat, R. Vuilleumier, and R. Spezia, Phys. Chem. Chem. Phys. **16**, 3693 (2014).

³⁴K. T. Tang and J. P. Toennies, The Journal of Chemical Physics **80**, 3726 (1984).

³⁵S. Tazi, J. J. Molina, B. Rotenberg, P. Turq, R. Vuilleumier, and M. Salanne, The Journal of Chemical Physics **136**, 114507 (2012).

³⁶CPMD version 3.13.2., http://www.cpmd.org/. Copyright IBM Corp 1990-2008, Copyright MPI fur Festkorperforschung Stuttgart 1997-2001.

³⁷J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters **77**, 3865 (1996).

³⁸N. Marzari and D. Vanderbilt, Phys. Rev. B **56**, 12847 (1997).

³⁹R. King-Smith and D. Vanderbilt, Phys. Rev. B **47**, 1651 (1993).

⁴⁰D. Vanderbilt and R. King-Smith, Phys. Rev. B **48**, 4442 (1993).

⁴¹I. Souza, T. Wilkens, and R. M. Martin, Phys. Rev. B **62**, 1666 (2000).

⁴²P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. **82**, 3308 (1999).

⁴³A. Aguado, L. Bernasconi, S. Jahn, and P. A. Madden, Faraday Disc. **124**, 171 (2003).

⁴⁴M. Barlet, A. Kerrache, J.-M. Delaye, and C. L. Rountree, Journal of Non-Crystalline Solids **382**, 32 (2013).

⁴⁵A. Grandjean, M. Malki, and C. Simmonet, Journal of Non-Crystalline Solids **352**, 2731 (2006).

⁴⁶X. Wu and R. Dieckmann, Journal of Non-Crystalline Solids **357**, 2846 (2011).

⁴⁷X. Wu, A. K. Varshneya, and R. Dieckmann, Journal of Non-Crystalline Solids **357**, 3661 (2011).

⁴⁸Y. Kim, Y. Yanaba, and K. Morita, Journal of Non-Crystalline Solids **415**, 1 (2015).

⁴⁹Y. Mirua, H. Kusano, and T. N. ans S. Matsumoto, Journal of Non-Crystalline Solids **290**, 1 (2001).

⁵⁰G. J. Martyna, D. J. Tobias, and M. L. Klein, The Journal of Chemical Physics **101**, 4177 (1994).

⁵¹G. J. Martyna, M. L. Klein, and M. Tuckerman, The Journal of Chemical Physics **97**, 2635 (1992).

⁵²G. J. Cuello, J. Darpentigny, L. Hennet, L. Cormier, J. Dupont, and B. Beuneu, Journal of Physics Condence Series (2016).

⁵³J. Hopf, S. Kerisit, F. Angeli, T. Charpentier, J. Icenhower, B. McGrail, C. Windisch, S. Burton, and E. Pierce, Geochimica et Cosmochimica Acta **181**, 54 (2016).

⁵⁴D. A. Keen, Journal of Applied Crystallography **34**, 172 (2001).

⁵⁵M. Barlet, J.-M. Delaye, T. Charpentier, M. Gennisson, D. Bonamy, T. Rouxel, and C. L. Rountree, Journal of Non-Crystalline Solids **417-418**, 66 (2015).

⁵⁶O. L. G. Alderman, G. Ferlat, A. Baroni, M. Salanne, M. Micoulaut, C. J. Benmore, A. Lin, A. Tamalonis, and J. K. R. Weber, Journal of Physics: Condensed Matter **27** (2015).

⁵⁷W. H. Zachariasen, Acta Crystallographica **16**, 385 (1963).

⁵⁸D. I. Grimley, A. C. Wright, and R. N. Sinclair, Journal of Non-Crystalline Solids **119**, 49 (1990).

⁵⁹H. F. Poulsen, J. Neuefeind, H.-B. Neumann, J. R. Schneider, and M. D. Zeidler, Journal of Non-Crystalline Solids **188**, 63 (1995).

⁶⁰R. L. Mozzi and B. E. Warren, Journal of Applied Crystallography **2**, 164 (1969).

⁶¹R. L. Mozzi and B. E. Warren, Journal of Applied Crystallography **3**, 252 (1970).

⁶²C. T. Prewitt and R. D. Shannon, Acta Crystallographica **B24**, 869 (1968).

⁶³R. T. Downs, H. Yang, R. M. Hazen, L. W. Finger, and C. T. Prewitt, American Mineralogist **84**, 333 (1999).

⁶⁴L. Pedesseau, S. Ispas, and W. Kob, Physical Review B **91** (2015).

⁶⁵R. Vuilleumier, N. Sator, and B. Guillot, Geochimica et Cosmochimica Acta **73**, 6313 (2009).

⁶⁶D. A. Kilymis, J.-M. Delaye, and S. Ispas, The Journal of Chemical Physics **143**, 094503 (2015).

⁶⁷G. Malavasi, M. C. Menziani, A. Pedone, and U. Segre, Journal of Non-Crystalline Solids **352**, 285 (2006).

⁶⁸M. Barlet, J.-M. Delaye, B. Boizot, D. Bonamy, R. Caraballo, S. Peuget, and C. L. Rountree, accepted by Journal of Non-Crystalline Solids (2016).
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /></figure>
<figure><img src="image_19.png" /></figure>

NBS17-24

NBS14-16

NBS27-21

NBS35-19

NBS29-13
<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>