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 analogues, 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 the present 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 modelling 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 dependance 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,7 . 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,[8][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 abinitio 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 8 , 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 modelling 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-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 36 . 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 Hoover 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 perfomed 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.
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 discret 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). In Na 2 SiO 3 and α−Na 2 Si 2 O 5 the 4-fold coordinated Si atoms are connected to two and three bridging (i.e. double-coordinated) oxygens, respectively 48,49 , 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 Tab. I. The cell constants are predicted within 5% compared to experiment 50,51 . 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 52 .
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  In Tab. 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     17,49 . 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. where N is the number of atoms in the system. The total VDOS has been then computed as: 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): where N α is the number of atoms in the group labelled α 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  heat-capacity calculations and concern the lower end of the spectra 64 , and show an intensity increase with increasing sodium content. For the color version, the reader is referred to the web version of the article.
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. 4a), a feature that originates from the changes presented by both oxygen and silicon partial VDOSs (see Figs. 4b and c) 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. 4b. 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  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 Fig. 7d-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. 7e). 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. 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 cm −1 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: 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: with F Iγ is 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: where α p and τ p are the trace and anisotropy of the susceptibility: The depolarized Raman spectra (VH) have been calculated by multiplying the VV spectra by the depolarization ratio ρ p dep : 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: 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 Fig. S1), 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 Suppl. Mat.), 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. 11b 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 intensitiy 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 66 . 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. 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.
An investigation of the per-atom NBO contributions in the high-frequency band, given in Fig. 14 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  Inset: Total contribution from NBO's on Q 2 and Q 3 species for the three models.

IV. CONCLUSION
We have presented in this work the investigation of the vibrational spectra, VDOS and Raman, for three sodosilicate glasses. The visualization of the vibrational modes has been 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 part 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.