Dimitrios Kilymis

Simona Ispas
email: simona.ispas@umontpellier.fr

Bernard Hehlen
email: bernard.hehlen@umontpellier.fr

Sylvain Peuget

Jean-Marc Delaye

Vibrational properties of sodosilicate glasses from first-principles calculations

The vibrational properties of three sodosilicate glasses have been investigated in the framework of density functional theory. The pure vibrational density of states has been calculated for all systems and the different vibrational modes have been assigned to specific atoms or structural units. It is shown that the Na content affects several vibrational features as the position and intensity of the R band or the mixing of the rocking and bending atomic motions of the Si-O-Si bridges. The calculated Raman spectra have been found to agree with experimental observations and their decomposition indicated the dominant character of the non-bridging oxygen contribution on the spectra, in particular for the high-frequency band, above 800 cm -1 . The decomposition of the high-frequency Raman feature into vibrations of the depolymerized tetrahedra (i.e. Q n -units) has revealed spectral shapes of the partial contributions that cannot be accounted for by simple Gaussians as frequently assumed in the treatment of experimentally obtained Raman spectra.

I. INTRODUCTION

Raman and infra-red (IR) spectroscopies are well-established characterization tools of the silicate glasses and melts in both research laboratories and industries. Although the Raman and IR spectra are easy to collect, the spectral responses are complex, combining broad and overlapping peaks due to the inherent structural disorder of these materials. Consequently, the decomposition and extraction of structural information from Raman and IR data are often based on parallelisms to well-known crystalline analogs, and result in a phenomenological and mostly qualitative approach 1,2 . A good and quantitative understanding of the structural significance of the features of the Raman spectra, as well as that of the two other vibrational spectra (IR and inelastic neutron scattering spectroscopies) can be achieved theoretically using atomistic calculations, in particular ab initio ones. In this context, the major aim of this study is to use atomistic simulation approaches in order to improve our ability to connect measured vibrational properties of binary silicates to their structural features.

The great advantage of the atomistic modeling is that the positions of the atoms in the system under study are continuously monitored and the structure is exactly known.

Therefore, correlations can be extracted between the atomic structure and the vibrational properties in order to arrive to a more quantitative interpretation of the vibrational properties and their dependence on the chemical composition, pressure, quench rates etc. The latter is not directly accessible in experiments, but derived through functions describing how the modes with a given frequency are coupled to the probing radiation (photons, neutrons, etc.) [3][4][5] .

The vibrational density of states (VDOS) can be in practice accessed either through the diagonalization of the dynamical matrix of the glass structure (in a local potential minimum), or the Fourier transform of the velocity autocorrelation function. Within a classical approach, most of the empirical potentials existing in the literature are quite reliable from the point of view of the structure but they fail to reproduce the experimental VDOS, even for systems as simple as pure silica or germania 6,[START_REF] Huang | Molecular Dynamics Simulations of Disordered Materials[END_REF] . Ab initio methods are much more accurate in the prediction of the vibrational properties of materials but they have been scarcely used in this scope, due to the high computational cost associated with the treatment of glass systems 4,6,[START_REF] Sarnthein | [END_REF][9][10][11][12][13][14][15][16][17][18] . More specifically, the fact that one has to model a system which lacks symmetry makes the calculations more time-consuming and finally results to systems that at best contain a few hundreds of atoms.

The calculations of the Raman spectra of a glass have in general three stages: (i) generation of a structural model; (ii) the calculation of its vibrational properties; and (iii) the calculation of derivatives of the polarizability tensor with respect to atomic displacements (i.e., the Raman susceptibilities). For the latter, two frameworks can be considered, either a first-principles one 19,20 or an empirical one, called bond polarizability model 21,22 . Within the bond polarizability approximation, calculations for silica and binary sodo-silicate glasses have achieved satisfying agreement with experimental data [23][24][25][26] . We note equally that Raman intensities can be equally computed from the dynamical autocorrelation functions of the dielectric tensor in a molecular dynamics (MD) simulation 5,27 approach which can be applied to liquids as well as to anharmonic solids, but is very demanding, requiring the calculation of the dielectric tensor at each MD step.

Concerning the vibrational spectra of silicate glasses, most work carried out within an ab initio approach has focused on the simplest case of vitreous SiO 2 . The early work of Sarnthein et al. elucidated the nature of the high-frequency modes [START_REF] Sarnthein | [END_REF] , while more recent works further elaborated on the VDOS of silica, as well as its IR and Raman spectra 4,9,10,18 . For more complex cases, theoretical calculations and analysis of the vibrational spectra have been carried out for sodosilicates 25 , borosilicates 17 , aluminosilicates 13,28 , phosphosilicates 14,29,30 and magnesium-silicates 16 .

In this work, we report on the vibrational properties of three binary sodosilicate glasses whose structures have been prepared by combining classical and ab initio MD simulations.

We present their VDOS and its decompositions according to types of vibrations and contributions from specific structural units. At the same time we report results on the total and partial Raman spectra, aiming to help resolve the origins of different vibrational modes and better interpret experimental results. As Raman spectroscopy probes the short-and medium-range order in glasses, this modeling approach may supply quantitative criteria for the spectral analysis.

II. SIMULATION METHODS

The selected sodo-silicate glasses of this study contain 20.0%, 25.0%, and 33.3% Na 2 O, hereafter called NS4, NS3 and NS2, respectively. Two samples per composition have been prepared, containing 180, 204, and 207 atoms for NS4, NS3, and NS2, respectively. In order to circumvent the costly preparation within an ab initio molecular dynamics (MD) simulation, we have generated our structural models by employing classical MD approach during the melt-and-quench cycle followed by a rather short ab initio MD run at room temperature.

As it has been shown in past studies conducted using this combined approach 4,13,[31][32][33][34] the use of the ab initio approach mostly serves as a way to refine the short-range structure, while the medium range remains practically unchanged. The switch from a classical to an ab initio approach can have a dramatic effect on the calculated vibrational properties 6 . Similarly, Giacomazzi et al. extracted the vibrational properties of silica glass prepared with both methods with the results suggesting that they generally produce similar results, with some discrepancies arising due to the slightly different Si-O-Si angle distributions 4 .

The initial preparation of the NS2 and NS3 glass structures, as well as that of one NS4 structure, has been performed using classical MD by employing the Guillot-Sator potential, with a cutoff set at 6 Å35 , starting from a random atomic distribution in a cubic box reproducing the experimental density [START_REF] Bansal | Handbook of Glass Properties[END_REF] . The systems were then heated to 4000 K and equilibrated for 20 ps using the constant volume-constant temperature (NVT) ensemble, before being quenched stepwise to 300 K with a rate of 10 11 K/s using NVT. Due to an inaccurate description of the medium-range order in the case of the NS2, a third system has been prepared this time using a 10 10 K/s quench rate. Two final equilibrations of 20 ps and 5 ps were finally performed, using the constant pressure-constant temperature (NPT) ensemble and the constant volume-constant energy (NVE) ensemble, respectively. Due to slight contractions of the glasses during the NPT phase, some of the final structures were rescaled in order to match the experimental density values, and were then re-equilibrated for 5 ps using the NVE ensemble. The classical MD simulations were carried out with Berendsen thermostats and barostats, using an integration step equal to 1 fs. They were carried out using the DL POLY 4.01 software 37 . As a second structure for the NS4 composition, we have used one previously generated, similarly through combined classical and ab initio MD simulations, albeit with a different classical potential 38 .

In a second step, we used the above classical structural models as input configurations for ab initio MD simulations carried out using the Vienna ab initio package (VASP) 39,40 . The electronic degrees of freedom were treated in the framework of density functional theory (DFT), using the general gradient approximation (GGA) and the PBEsol functional 41,42 .

We started with an NVT run at room temperature for 2 ps, followed by an NVE run of same length. Complete computational details can be found in previous publications 17 . For the second NS4 structure, the ab initio MD simulations have been performed using Car-Parrinello MD as implemented in CPMD code (see simulation details in Ref. 38 ).

Further, in a third step, we have calculated the vibrational frequencies and eigenmodes, as well as the Raman tensors, using a DFT scheme as provided in the QUANTUM ESPRESSO (QE) package (version 5.4) 43 . The atomic positions have been initially relaxed using the BFGS minimizer with a fixed simulation box size setting the convergence criteria for the total energy and forces at 4 × 10 -8 Ry and 10 -5 Ry/a.u., respectively. In case of a pressure higher than 1 GPa after the relaxation, the box was slightly suitably rescaled, keeping the final densities within 1.5% of the experimental value. The calculations were carried out using the GGA-PBE exchange-correlation functional 41 with Troullier-Martins norm-conserving pseudopotentials 44 . The performance of the pseudopotentials was assessed by investigating the structure and vibrational properties of two crystalline phyllosilicates, the sodium metasilicate (Na 2 SiO 3 ) and natrosilite (α-Na 2 Si 2 O 5 ). Calculations on the crystalline compounds were carried out using a kinetic energy cutoff of 100 Ry and Monkhorst-Pack k-point grids with a 4 × 4 × 4 and 3 × 4 × 3 resolution for the metasilicate and natrosilite, respectively.

For the glasses, the same kinetic energy cutoff was used and calculations were carried out at the Γ-point. We should note at this point that, although two different exchange-correlation functionals have been used in the ab initio MD runs, the final atomic relaxation has been carried out using PBE in all cases. Therefore, we consider that the final sets of positions used for the vibrational properties are consistent. As documented by Ganster et al. 13 , the ab initio runs at the later stages of the preparation of a glass model only serve as a way to refine the short-range order and do not affect the medium-range one.

The vibrational properties (eigenfrequencies and eigenmodes) were calculated using the density-functional perturbation theory method (DFPT) as implemented in QE code 43,45 . All spectra for the NS3 and NS4 systems presented in the Sec. III have been averaged over the two models for the corresponding composition, after broadening the discrete spectra using Gaussians with a 20 cm -1 FWHM, unless explicitly stated. In the case of NS2, the results only concern the system prepared with the slower quench rate (i.e., 10 10 K/s).

III. RESULTS AND DISCUSSION

A. Crystalline compounds

In Na 2 SiO 3 and α-Na 2 Si 2 O 5 the fourfold-fold coordinated Si atoms are connected to two and three bridging (i.e., double-coordinated) oxygens, respectively 46,47 , usually called Q 2 and Q 3 units (we recall that Q n species are tetrahedral environments with n = 4, 3, 2, 1, 0 BOs). Unit cells containing 12 and 36 atoms for Na 2 SiO 3 and α-Na 2 Si 2 O 5 were relaxed, with the results presented in Table I. The cell constants are predicted within 5% compared to experiment 48,49 . The intra-tetrahedral structural parameters are more accurately predicted than the inter-tetrahedral ones, whereas the largest deviations from the experimental values are observed in the case of Na-O bond lengths. The overestimation of these bond lengths is probably connected to the use of GGA/PBE approximation, as already reported in the literature 50 . The Raman spectra for the two systems are given in Fig. 1. As there is a constant underestimation of the frequencies of the calculated spectra, the presented values have been multiplied by a factor of 1.04. For both crystal structures, we observe a good agreement of the band positions with respect to available experimental data [51][52][53][54] , except from the weak Si-O stretching band of Na 2 SiO 3 at 1063 cm -1 which is not present in our calculations. However, the overall good description of the spectra suggests that the chosen pseudopotentials and simulation parameters are adequate for the calculation of the spectra for the glasses. When switching from classical to ab initio MD, we have observed changes of some structural features such as the shift to smaller angles shown by the SiOSi bond angle distribution (BAD) (see Fig. 2) or the shifts of the first peaks of SiO and NaO radial distribution functions (RDF) plotted in Fig. 3 for NS3. For the SiOSi BADs plotted in Fig. 2, their average value does not strongly depend on the sodium content (see Table II), but their maximum does, as its location shifts to lower values with increasing Na content. This trend is consistent with experimental observations and is particularly apparent for the ab initio data 38 .

For the SiO RDF, the first peak shifts to larger values together with a broadening, which may reflect the deformation of SiO 4 tetrahedra bearing NBOs, known to have shorter bond lengths with Si atoms compared to Si-BO bonds 17,32 . Concerning the Na-O RDF, Fig. 3 shows a decreasing of NaO bond lengths, with an average NaO bond length moving from 2.73 Å during the classical simulation to 2.50 Å for the ab initio trajectory. This variation improves the comparison with experimental values (see Table II), but overall it seems that, as for the crystal structures, these bond lengths are overestimated, reflecting the already mentioned trend of ab initio calculations using the GGA approximation. In Table II we summarize some structural parameters of the three glass models used for the vibrational calculations (i.e., relaxed to 0 K). The amount of NBO in the structures was found to be in very close agreement to available NMR results [55][56][57] . The percentages of the NBOs are accurately reproduced in all glasses, while the agreement of the Q n species with experimental results is better for the well-polymerized NS4 glass than for NS2. Overall, one can state that the short-range structural features of our models compare reasonably well with experiments, but deviations exist in terms of distribution of the Q n species, and thus for the medium-range structure. This discrepancy can be attributed to the choice of the empirical potential during the initial preparation of the glass models, since the ab initio equilibration at room temperature and the subsequent energy minimization produced negligible changes in the medium-range structure. It can be also linked to the employed high quench rate from the liquid.

It is also interesting to analyze the intertetrahedral Si-O-Si angles with respect to the Q n speciation of the two silicons. In this scope, we have averaged the corresponding values over all five glass models, wherever sufficient statistics were available, and found the largest well related with a steric effect caused by the increased amount of Na around the Q n entities when increasing the number of NBOs. We have hence calculated the average number of Na atoms around a Q n unit (i.e., the Na atoms around the four oxygens and within a sphere of radius equal to the first O-Na coordination shell), and we have found increasing values when going from Q 4 to Q 2 , for a given composition, and when increasing the soda content.

angles being formed by Q 4 -O-Q 4 (144.2 • ), followed by Q 4 -O-Q 3 (143.8 • ), Q 4 -O-Q 2 (141.5 • ), Q 3 -O-Q 3 (140.1 • ), and finally Q 3 -O-Q 2 (137.3 • ). This change in the Q m -O-Q n angle could be
For example, for Q 3 species, we have found ≈ 3.8 sodiums for NS4, ≈ 4.4 sodiums for NS3 and finally ≈ 5.8 for NS2 composition. Further, by calculating the average number of Na around each NBO, we have found values equal to 2.9, 3.2, and 3.4 for NS4, NS3, and NS2, respectively, similarly to previous ab initio calculations on NS4 and a sodium borosilicate glass 17,47 . Of course, this implies that more than one alkali atom is found in the vicinity of each NBO, therefore deviating from the traditional theoretical model for the structure of sodosilicates, which assumes a one-to-one correspondence.

C. Vibrational density of states

As stated in Sec. II, the dynamical matrix of each glass model was computed within the framework of DFPT, and then its diagonalization yielded the 3N eigenfrequencies ω p and their corresponding normalized 3N -component eigenvectors e(ω p ) (p = 1, 2, 3, . . . , 3N ),

where N is the number of atoms in the system. The total VDOS has been then computed as:

g (ω) = 1 3N -3 3N p=4 δ (ω -ω p ) (1)
The knowledge of the eigenfrequencies and eigenvectors makes possible further calculations, starting with the so-called partial VDOS corresponding to the contributions of different groups of atoms (as, for example, the species):

g α (ω) = 1 3N -3 3N p=4 Nα I=1 3 k=1 |e I,k (ω p ) | 2 δ (ω -ω p ) (2)
where N α is the number of atoms in the group labeled α and e I,k (ω p ) are the 3-components of the eigenvector e(ω p ) with the displacement of the particle I. The total VDOSs of NS4, NS3 and NS2 glasses are shown in Fig. 4 alongside the partial contributions of the three species α = Si, O, Na. Firstly we can identify three main bands with features (position, shape and intensity) changing when the Na content varies. The presence of three bands as well as qualitatively similar decompositions with respect to the contributions from the three species have been already reported in previous theoretical studies for NS4 glass using classical and ab initio 23,25 , as well as by fully ab initio studies of more complex silicate glasses 17,63 . For the low frequency range below 300 cm -1 , we observe an intensity increase when going from NS4 to NS2, which is directly correlated to the intensity increase presented by the Na partial VDOSs when the Na sodium content increases [see The change in the Na amount also gives rise to modifications of the total VDOSs for the frequency ranges above 400 cm -1 . We notice an intensity decrease with increasing Na concentration in the range 400 -500 cm -1 [see Fig. and can be correlated to changes in the glass polymerization due to the modifier atoms.

Further on, there is a peak pointing at ≈ 740 cm -1 which is strongly affected by the Na content: its intensity decreases with increasing Na amount, alongside a slight shift to lower frequencies. Again, we can link these variations to changes shown by both Si and O partial VDOS, and especially the ones of the silicons with a decrease of the maximum intensity of 40% between NS4 and NS2 samples, as shown in Fig. 4 . This band is of course reminiscent of the so-called highfrequency band in pure silica, and its softening due to the silicate network depolymerization has been previously reported in ab initio studies for other silicates 16,17,63 . The assignment of this band and the relation with the network connectivity will be discussed in details in the next subsections.

We can initially proceed to a visual inspection of the vibrational modes (videos are provided as Supplemental Material 65 ). The nature of the vibrations was found to be qualitatively similar in all three glasses under study. Figure 5 depicts some characteristic snapshots taken from the NS4 models and visualized using XCrySDen 66 . For the low-end of the VDOS we observe mainly motions of rigid tetrahedra together with neighboring sodium atoms. we mostly observe the stretchings of both Si-BO and Si-NBO bonds, as previously reported for a NS4 glass 23 within a classical approach or for magnesio-silicate glasses 16 within a DFT calculation.

After obtaining a qualitative overview of the vibrational modes of the NSX glasses, we can proceed to a more quantitative analysis. As mentioned earlier, the most intense contribution to the VDOS stems from the O atoms. In Fig. 6, we present the decomposition of the O partial VDOS with respect to BOs and NBOs contributions. The results naturally show the relative decrease of the BO contribution in the entire spectrum, as the glasses become more depolymerized, and the simultaneous increase of the NBO contribution. We notice that around 700 cm -1 the BO contribution remains dominant, even in the case of the NS2 glass, and, as we will see below, the vibrational modes in this region of the spectra result from a superposition of rockings/bendings of the Si-O-Si bridgings.

We have further analysed the vibrational modes of the Si-O-Si bridges according to the decomposition of the BO motions into rocking, bending, and stretching components. For each mode, we have calculated the relative displacement of the BO with respect to its two neighboring Si atoms and projected it over three orthogonal directions 67 the relative contributions of rocking, bending, and stretching motions, respectively, with respect to the total VDOS. The data plotted in Fig. 7 show a mix of rocking and bending motions in the low-and mid-frequency range, with the former being more intense in the region between 100-400 cm -1 , whereas the latter are relatively more pronounced between 500-800 cm -1 . On the contrary, the high frequency band contains essentially stretching modes of BO. Regarding the effect of the glass composition, the increase of Na amount when going from NS4 to NS2 leads to an increase of the contributions from rocking and bending modes at the lower part of the VDOS (below 700 cm -1 ). At the same time, their intensities for the band around 750 cm -1 decrease, and also shift to lower frequencies (the effect being more pronounced for the relative contributions, as shown in Figs. 7(d)-7(e). The increase of Na has also a limited effect on the intensity of the high-frequency band which only shifts to lower frequencies [see Fig. 7(e)]. Similar decompositions of the oxygen motions have been previously carried out for pure silica described using an effective pair potential 26,67 or within an ab initio framework 10,68 , as well as for the NS4 glass using a harmonic potential 23 .

These studies have revealed the same strong dependence of BO motions with respect to the frequency, and a qualitative agreement with respect to the relative ratios between the three contributions. inverse effect is present in the high-frequency band, where a decrease of the Si-O-Si angle produces a downshift of the average contribution to the VDOS. In this case, the downshift can be straightforwardly linked to stretching in Q n species since, as will be discussed in Section III B, smaller Si-O-Si angles are expected as the amount of NBOs on the Si atoms increases.

Finally, we have investigated the contribution from the NBOs by decomposing it with respect to the type of Q n species that they are linked to. Figure 9 shows the contributions to the VDOS from NBOs on Q 2 and Q 3 entities. We remark that NBOs on Q 2 demonstrate a significantly higher intensity, compared to those on Q 3 for all three glasses. However, the most important difference between the two cases can be identified in the high-frequency band, where in all systems we observe a bimodal contribution for NBOs on Q 2 species and an overlapping peak for Q 3 . In the former case, the peak centers were found to be at the same frequencies, 920±10 cm -1 and 1020±10 cm -1 , whereas the maximum of the Q 3 contribution shifts from 1060 to 1010 cm -1 when going from NS4 to NS2. We should also note that we do not comment on the contribution to the VDOS from NBOs on Q 1 entities, since only one exists in the NS2 system and therefore good statistics are not ensured.

D. Raman spectra

The Raman spectra have been obtained by applying uniform electric fields and calculating the second order derivatives of the electron density matrix 19 . The Raman susceptibility for each mode p is given by:

R p ij = √ V N I=1 3 k=1 ∂χ ij ∂R I,k e I,k (ω p ) √ M I ( 3
)
where V is the system volume, i, j, k = 1, 2, 3, M I the atom I mass, while e I,k (ω p ) is the 3-components vector defined in Eq. ( 2). In Eq. ( 3), the Raman tensor ∂χ ij ∂R I,k is given by:

∂χ ij ∂R I,k = 1 V ∂ 2 F I,k ∂E i E j E=0 (4)
with F Iγ the force exerted on atom I and E the applied electric field. For the polarized Raman spectra (VV), the intensity for each mode p is given by:

I p V V = α 2 p + 4τ 2 p 45 (5)
where α p and τ p are the trace and anisotropy of the susceptibility:

α p = R p 11 + R p 22 + R p 33 3 , (6)
τ 2 p = (R p 11 -R p 22 ) 2 + (R p 22 -R p 33 ) 2 + (R p 33 -R p 11 ) 2 + 6 (R p 12 ) 2 (R p 13 ) 2 + (R p 23 ) 2 2 , (7)
The depolarized Raman spectra (VH) have been calculated by multiplying the VV spectra by the depolarization ratio ρ p dep :

I p V H = ρ p dep • I p V V with ρ p dep = 3τ 2 p 45α 2 p + 4τ 2 p ( 8
)
In order to decompose the Raman spectra into partial contributions of a specific group of atoms, we have written the intensity of mode p as a sum of three terms:

I p V V = α I p V V,α + I p V V,comp + I p V V,over (9)
where only the atoms of type α are considered in Eq. ( 3) prior to the calculation of the partial intensities I p V V,α using Eqs. ( 5)-( 7) (i.e., α = Si, O, Na or α = BO, NBO). The second and third terms are respectively the contributions from the remaining atoms in the system (i.e. the ones of Si and Na if α = BO, NBO) and the overlapping (or interference) terms, due to the squared terms in Eqs. ( 5)- (7). As already pointed out by Umari and Pasquarello for the case of pure SiO 2 68 , even if the sum of the partial Raman spectra does not recover the total spectrum, as in the case of the VDOS, it is reasonable to assume that these partial spectra give a measure of the relative ones.

The VV and VH calculated Raman spectra of the glasses under study are given in Fig. 10 and compared to recent experimental results 69 . The calculated spectra result directly from the Gaussian broadening of the discrete intensities [see Eqs. ( 5) and ( 8)], while the experimental ones have been obtained by multiplying the experimentally measured intensities by ω/ [n (ω) + 1], where n (ω) = 1/ [exp (ω/k B T ) -1)] is the Bose factor, k B the Boltzmann constant, and T = 300 K. In the Supplemental Material (see Fig. S1) 65 , we also provide a comparison between the calculated spectra and the experimental ones, using an alternate definition for the reduced intensities. In that case, the latter have been multiplied by 1/[ω(n(ω) + 1)], which is equivalent to dividing the calculated spectra in Fig. 10 by ω 2 .

For either VV and VH spectra, the positions of the main bands are well reproduced by the calculated spectra but there exist discrepancies in terms of band intensities. If we consider the band between 500-600 cm -1 in the VV spectra, its form and intensity is well reproduced for the highly polymerized (high silica content) NS4 composition. As it will be discussed later, this band contains the signal from Si-O-Si vibrations, i.e. the BO atoms.

The calculated band is broader with respect to the experimental one and it is not narrowing with increasing sodium content. This indicates that there are differences in the local and medium range structures present in our glass models and those in the real glasses. These differences are likely due to the small size of our glass models, which in the case of NS2 contains considerably less BOs and therefore poorer statistics compared to NS4. Further, the discrepancies for the high-frequency band (900-1200 cm -1 ), corresponding to the response from the Q n species are due to the fact that observed Q n distributions are not exactly reproduced in our models. For the VH spectra, we see a broad band below 600 cm -1 (see Fig. S1 in Supplemental Material 65 ), which was assigned to two types of alkali motion 69 , one at ≈175 cm -1 and the second one at ≈350 cm -1 . In spite of the small size of our models, we can identify in the calculated VH spectra the presence of small peak at ≈350 cm -1 showing an increasing intensity with increasing sodium content. As for the pure VDOS, we have initially decomposed the VV spectra with respect to the contributions of the three species, Si, O, and Na and present the results in Fig. 11. The partial contributions stemming from the O atoms are by far the most significant, as it was already noticed for pure silica 68 , as well as in the partial VDOS of our three sodosilicates (see Fig. 4). This trend is far more substantial in the polarized Raman spectra. Consequently, there is strong similarity between the form of the oxygen contributions plotted in Fig. 11(b) and the total VV Raman spectra shown in Fig. 10. With increasing sodium content, it is the oxygen contribution exhibiting the most noticeable changes, namely, a shift to the right of the band between 400-700 cm -1 and a small shift to the left of the high-frequency band.

Due to the prominent contribution of oxygen vibrations to the Raman spectra, we have further considered three decompositions: according to their species (BO or NBO), the bond angle of the Si-O-Si bridgings, and the contribution of the NBOs to the high frequency band.

The decomposition according to the contributions from BOs and NBOs is given in Fig. 12.

We initially point out the strong contribution of the NBOs to the high-frequency, or Q n , band. This happens in spite of their smaller concentration with respect to the BOs in any of the three compositions, suggesting that they are the most Raman active in this band. The intensity of these modes largely dominates the one of the oxygens in Si-O-Si bridges, that one finds in the pure silica glass. The vibrations of NBOs then define the shape and intensity of this Raman band, and we expect that they are predominantly of stretching type 23 . For the rest of the spectra (below 800 cm -1 ) the NBOs present a distinct signal whose intensity increases with increasing sodium content, and surpasses the one of BOs in the case of NS2.

As for BOs, they present a significant band between 400-600 cm -1 , particularly for low sodium content. Taking into account the decomposition of the pure VDOS with respect to rocking, bending and stretching motions of the Si-O-Si bridges (see Fig. 7), the Raman BO vibrations could be assigned either to rocking or bending motions, or both. However, for pure SiO 2 , Umari and Pasquarello 68 have shown that the VV spectrum below 900 cm -1 is dominated by bending motions, while the rocking contribution is almost suppressed, in spite of their large weight in the pure VDOS. This conclusion likely holds also for sodosilicate compositions. It is supported by the fact that the relative contributions of the rocking, bending, and stretching vibrations [see Figs. 7(d)-7(f)] do not seem to show a strong dependence on sodium content, on one side, while their frequency dependence also presents a close resemblance to the one reported for a-SiO 2 (see Ref. 68 ). Finally, if we consider the contribution of the BOs to the Q n band, our calculations show that the presence of sodium leads to the suppression of the double-peak feature, well known for pure SiO 2 glass. An investigation of the per-atom NBO contributions in the high-frequency band, given in Fig. 14, shows that the Q n contributions are significantly overlapping. Regarding the NBOs on Q 1 species, the poor statistics caused by the single case observed in the NS2 sample makes us refrain from further comments. In the case of NBOs on Q 2 species, we observe a welldefined doublet structure with maxima around 920 and 1030 cm -1 , whose position does not seem to be strongly affected by composition. The most probable reason for the appearance of the doublet is the presence of neighboring effects from different Q n species linked to Q 2 ones. However, our analysis with respect to the populations of Q n neighbors did not reveal a clear correlation. As for the NBOs on Q 3 species, they exhibit a unimodal distribution whose maximum shifts from 1080 cm -1 in NS4 to 1030 cm It could be argued that the bimodal distribution of the Q 2 contribution to the Raman spectra results from the relatively low total number of these units which lead to poor statistics. In order to clarify this point, we present in Fig. 15 the total polarized Raman spectra for the two NS2 systems prepared with a quench rate of 10 Inset: total contribution from NBO's on Q and Q 3 species for the three models.

NBO Q 3 NBO Q 2 NBO Q 1 BO a) b) c) 930
used to make an initial analysis in order to show that the experimentally predicted modes are reproduced in the calculations. The per-species decomposition of the VDOS revealed the relatively high contribution from O atoms in the whole frequency range, and the contribution of the Na atoms at low frequencies, typically below 400 cm -1 . The mode analysis of the Si-O-Si bridges has shown that the vibrations up to 800 cm -1 are of rocking and bending type, whereas above that frequency they are essentially stretchings. The average contributions of the bridgings to the VDOS seems to be affected by the angle formed by the three atoms, with a slight upshift around 500 cm -1 and a downshift at the high-frequency band as the angle becomes smaller. Concerning the high-frequency band, the results reveal important overlaps between NBO contributions in Q 2 and Q 3 species, as well as with those of the BOs.

The Raman spectra reproduce well the experimental band positions, whereas the calculated intensities show some deviations. We have shown that this is due to the limited system size and the prediction of the Q n species distribution by the classical potential which was used to quench the glasses from the melt. The calculation of the Raman response by Si-O-Si bridges has shown that large angles have a particularly flat contribution to the lower and middle parts of the spectra while, similarly to the VDOS, relatively small angles exhibit an upshift around 500 cm -1 . The calculation of the BO and NBO partials to the Raman spectra clearly shows that the form of the Q n band is defined by the latter. The further analysis of the NBOs in terms of their local environment shows a bimodal contribution with well-defined maxima for those connected to Q 2 and a highly overlapping Q 3 peak. The fact that the BO partial is also overlapping and spans the entire width of the Q n band indicates that one has to be particularly prudent when deconvoluting experimental spectra with the use of simple functions.

FIG. 1 .

 1 FIG. 1. Polarized Raman spectra for a) Na 2 SiO 3 and b) α-Na 2 Si 2 O 5 . Experimental spectra are from Refs. 51,52 .

B.

  FIG. 2. SiOSi BADs of the simulated NS4, NS3, and NS2 glassy samples at 300 K: (a) (Top) Extracted from VASP ab initio MD runs; (b) (Bottom) Extracted from DL -POLY classical MD runs.

FIG. 3 .

 3 FIG. 3. Radial distribution functions for SiO and NaO pairs computed for the classical and ab initio runs, at 300K and for NS3 glass structures.

FIG. 4 .

 4 FIG. 4. (a) Total VDOSs of the NSX glasses and partial contributions from (b) Si,(c) O and (d) Na atoms. The total VDOSs were normalized to unity, while the partial VDOS for a given composition were normalized so that their sum gives the corresponding total VDOS.

Fig. 4 (FIG. 5 .

 45 Fig.4(d)]. We note that the Na contributions are situated only in the bottom part of the spectra below 400 cm -1 , as already pointed out in previous calculations for silicates containing sodium17,25,63 . In terms of relative intensity of the three partials, it is evident that the O contribution dominates the VDOS [see Figs.4(b)-4(d)]. To our knowledge, the only available experimental data in order to compare the VDOS come from heat-capacity calculations and concern the lower end of the spectra 64 , and show an intensity increase with increasing sodium content.

  4(a)], a feature that originates from the changes presented by both oxygen and silicon partial VDOSs [see Figs. 4(b) and 4(c)]

  (b). Finally, the three VDOSs are characterized by a broad band approximatively above 850 cm -1 and up to 1250 cm -1 , with contributions from Si and O atoms which shifts to lower frequencies with increasing Na content [see Figs. 4(a)-4(c)]

Further, for increasing

  frequencies up to 350 cm -1 , we first see dangling and then stretching motions of the NBO-Na bonds. At the same time, the Si-O-Si bridges exhibit a combined rocking/bending motion [Fig.5(a)]. At higher frequencies we can also observe characteristic collective motions like the breathing modes of 4-member and 3-member silicate rings, shown in Figs.5(b) and 5(c). In the region 700-800 cm -1 [see Fig.5(d)], we still have the mixed rocking/bending motions of Si-O-Si bridges but we equally notice an increased contribution stemming from Si atoms, which corroborates the above discussed features of the Si partial VDOS. At a lower extent, we can still identify the presence of a mixed motion on NBOs combining a stretching of the Si-NBO bond together with a motion orthogonal to the Si-NBO direction. Finally, in the high frequency band above 850 cm -1 [see Fig.5(e)],

FIG. 6 .

 6 FIG. 6. Contribution to the VDOS of BO (solid lines) and NBO (dashed lines) for (a) NS4, (b) NS3, and (c) NS2.

  : the rocking direction orthogonal to the SiOSi plane, the bending direction along the bisector of the Si-O-Si angle and finally the stretching direction orthogonal to the two previous ones. The corresponding partial VDOS are presented in Figs. 7(a)-7(c), while Figs. 7(d)-7(f) shows

FIG. 7 .

 7 FIG. 7. Frequency and composition dependencies of the (a) rocking, (b) bending, and (c) stretching components of BO modes in the three NSX glasses (X=4,3,2). (d)-(f) show the relative contributions of the rocking, bending and stretching motions, respectively, with respect to the total VDOS.

FIG. 8 .FIG. 9 .

 89 FIG. 8. Decomposition of the BOs contribution with respect to their Si-BO-Si angle for (a) NS4, (b) NS3, and (c) NS2.

FIG. 11 .

 11 FIG. 11. Partial contributions of three constituent species: silicons (a), oxygens (b), and sodiums (c).

FIG. 12 .

 12 FIG. 12. Partial contributions to the VV spectra of BO (solid lines) and NBO (dashed lines) atoms for a) NS4, b) NS3 and c) NS2.

FIG. 13 .

 13 FIG. 13. Contributions to the polarized Raman spectra of the BO atoms taking into account the value of the Si-O-Si bond angle for a) NS4, b) NS3 and c) NS2.

FIG. 14 .

 14 FIG. 14. Per-atom contribution to the Raman spectra of NBO's on different Q n -species for (a) NS4, (b) NS3, and (c) NS2. The per-atom contribution of BO's is also given for comparison. The maxima of the contributions for Q 2 and Q 3 species are noted on the plots.

Fig. 15 ) 3 FIG. 15 .

 15315 FIG.15. Polarized Raman spectra of the three NS2 samples obtained with different quench rates.

FIG. S1.

  FIG. S1. a) Reduced VV (top panel) and VH (bottom panel) Raman spectra for the three glasses (solid lines), alongside experimental results (dashed lines) 69 . The experimental reduced spectra were multiplied by 1/[ω (n (ω) + 1)], where n (ω) = 1/ (exp (ω/k B T ) -1)) is the Bose factor, k B the Boltzmann constant and T = 300 K. The calculated as well as the experimental spectra are normalized to 1.

TABLE I .

 I Cell size, average values of the bond lengths (Si-BO, Si-NBO and Na-O) and angles (Si-O-Si and O-Si-O) formed by the unit cell atoms for Na 2 SiO 3 and α-Na 2 Si 2 O 5 . Average percent errors for the angles and bond lengths are given in parentheses.

		Na 2 SiO 3	Exp. a	α-Na 2 Si 2 O 5	Exp. b
	a ( Å)	10.92	10.48	6.65	6.409
	b ( Å)	6.34	6.07	15.83	15.422
	c ( Å)	4.97	4.82	5.02	4.896
	Si-O-Si ( • )	139.2	133.7	162.1	160.0
				144.2	138.9
		(4.1)		(2.6)
	O-Si-O ( • )	103.9	103.1	108.6	109.1
		110.7	110.8	107.5	107.5
		116.5	116.9	112.9	113.2
		107.0	107.1	106.0	105.4
				108.6	108.2
				113.0	113.0
		(0.1)		(0.0)
	Si-BO ( Å)	1.689	1.677	1.621	1.609
		1.675	1.668	1.651	1.643
				1.646	1.638
		(0.6)		(0.6)
	Si-NBO ( Å)	1.606	1.592	1.590	1.578
		(0.9)		(0.8)
	Na-O ( Å)	2.436	2.282	2.536	2.386
		2.450	2.303	2.468	2.338
		2.609	2.549	2.480	2.373
		2.483	2.370	2.444	2.290
		2.550	2.404	2.629	2.600
		(5.3)		(4.8)

a Ref.

48

; b Ref.

49

TABLE II .

 II Density, NBO percentage, and Q n species percentages for the NSX glasses, as well average values of some bond lengths and SiOSi angle for the glass models used for vibrational , 50.0 e 43.1 61.0 b , 60.5 c , 62.3 e , 34.9 f 58.7 79.0 b , 81.0 f

	calculations.
			NS4		NS3		NS2
		Sim.	Exp.	Sim.	Exp.	Sim.	Exp.
	ρ (g/cm 3 )	2.38	2.38 a	2.40	2.43 a	2.48	2.49 a
	% NBO	21.8	23.1 b	28.6	27.0 b , 28.0 c ,28.0 d	40.0	40.0 b
	% Q 1					2.2	-
	% Q 2	8.3	3.0 b , 0.8 e	11.8	1.0 b , 2.3 c 2.4 e , 1.4 f	17.4	10.0 b , 8.4 f
	% Q 3 32.3 48.0 b % Q 4 59.4 50.0 b 49.2 e	45.1	38.0 b 37.2 c 35.3 e	21.7	11.0 b 10.3 f
	Si-O ( Å)	1.627	1.617 g	1.628	1.62 g	1.632 1.61-1.62 h , 1.631 i
	Si-BO ( Å)	1.635		1.638		1.649
	Si-NBO ( Å) 1.574		1.581		1.587
	Na-O ( Å)	2.522	2.30 h	2.512		2.551	2.30-2.36 h
	SiOSi ( • )	143.08		142.45		143.49

a Ref.

[START_REF] Bansal | Handbook of Glass Properties[END_REF]

; b Ref.

55

; c Ref.

56

; d Ref.

57

; e Ref.

58

; f Ref.

59

g Ref.

60

; h Ref.

61

, i Ref.

62

  4,70 ). For the high frequency band, these two main contributions show a stronger overlap with respect to the ones found in the pure VDOS. The contribution of BOs with angles between 140 • -160 • is rather asymmetric and seems to become stronger with increasing Na content for the highest frequencies, with respect to the one of BOs with angles smaller than 140 • . Finally we note that the contribution of the very large angles, above 160 • , is almost negligible below 850 cm -1 , and seems to slightly increase with increasing frequency in the Q n -band.

		8e-05
	Intensity (Arb. units)	4e-05 6e-05
		2e-05
		0e+00	800	900	1000	1100	1200	800	900	1000	1100	1200	800	900	1000	1100	1200
					Raman shift (cm -1 )				Raman shift (cm -1 )				Raman shift (cm -1 )

  -1 in NS2. The particular form of the NBO contributions for Q 2 and Q 3 species has been first pointed out by Zotov et al. for NS4 23 and is found to be in remarkable agreement to the more recent experimental works by Malfait et al. for sodosilicates 59 and Woelffel et al. for soda-lime silicates 71 . The Q 4 speciesgive rise to a very broad band ranging from 900 to 1200 cm -1 , underneath the Q 2 , Q 3

contributions. The strong overlaps between the different Q n entities, as well as the form of the individual contributions, indicate that the quite common decomposition of experimental data with gaussian peaks should be considered cautiously.

E. Comparison of the NS2 samples

  11 K/s (henceforth named model I and model II), alongside the one with 10 10 K/s which has already been discussed. Models I and II have very similar percentages of Q 2 and Q 3 species, 26.1% Q 2 and 28.3% Q 3 for model I, 24.4% Q 2 and 31.1% Q 3 for model II, which reflects to an almost identical structure of the Q n band. On the contrary, the NS2 glass prepared with the 10 10 K/s quench rate has a considerably lower amount of Q 2 entities, hence the lower signal around 950 cm -1 . Furthermore, the partial contributions from NBOs on Q 2 (inset of

  V. ACKNOWLEDGEMENTS D.K. would like to thank CEA for the financial support. This work was granted access to the HPC resources of TGCC/CINES/IDRIS under the allocation x2016097572 attributed by GENCI (Grand Equipement National de Calcul Intensif).

		0.003
		0.0025			NS4 exp. NS4 calc.				NS3 exp. NS3 calc.	NS2 exp. NS2 calc.
	Reduced Intensity (Arb. units)	0.001 0.0015 0.002			VV spectra				VV spectra	VV spectra
		0.0005
		0 0	300 600 900 1200	0	300 600 900 1200	0	300 600 900 1200
			Raman shift (cm	-1 )		Raman shift (cm	-1 )	Raman shift (cm	-1 )
	Reduced intensity (Arb. units)	0.001 0.002 0.003 0.004 0.005	VH spectra	NS4 exp. NS4 calc.			NS3 exp. NS3 calc. VH spectra	NS2 exp. NS2 calc. VH spectra
		0 0	300	600	900 1200 0	300	600	900 1200 0	300	600	900 1200
			Raman shift (cm	-1 )		Raman shift (cm	-1 )	Raman shift (cm	-1 )

SUPPLEMENTARY MATERIAL

D. Kilymis, S. Ispas, B. Hehlen, S. Peuget and J.-M. Delaye List of some videos corresponding to atomic vibrations at selected frequencies:

• Video 1: Motions of rigid tetrahedra and dangling-like motions of Na-O bonds (ω = 108 cm -1 ).

• Video 2: Breathing motions of oxygens atoms (highlighted in blue) in a 4-member ring, as well as bending motions of Si-O-Si bridges (ω = 477 cm -1 ).

• Video 3: Breathing motions of oxygens atoms (highlighted in blue) in a 3-member ring (ω = 610 cm -1 ).

• Video 4: Streching-like motions of oxygen atoms in high-frequency band (ω = 1130 cm -1 ) These videos were produced using Avogadro version 1.1.1 [START_REF]Avogadro: an open-source molecular builder and visualization tool[END_REF][START_REF] Hanwell | [END_REF] , and are available upon request from the authors.