Ab initio molecular dynamics simulation of Na-doped aluminosilicate glasses and glass-water interaction

The structure and properties of sodium aluminosilicate (NAS) glasses are investigated using ab initio molecular dynamics and density functional calculations. Four NAS glass models of about 700 atoms with composition (SiO2)0.6(Al2O3)0.4-x(Na2O)x with Na/Al ratio R = 0.0, 0.5, 1.0 and 1.5 are constructed corresponding to x = 0, 0.135, 0.20 and 0.24. Detailed information on network coordination, electronic structure, interatomic bonding and partial charge distribution, mechanical and optical properties of these models are presented and fully analyzed. The structural details for each R are discussed in terms of short- and intermediate-range order manifested in the coordination number, atomic pair and bond angle distributions. It is shown that the mechanical strength of NAS glasses decreases with increasing Na content, indicating that pure aluminosilicate glass is stronger than the alkali-doped glasses. We use the novel concept of total bond order density to characterize the internal cohesion of the NAS glasses. In the case of R = 1 NAS model, 12 water molecules are added to investigate the solvation effect and hydrolysis in NAS glass.


I. INTRODUCTION
Alkali aluminosilicate is an important segment of inorganic glass with numerous industrial and technological applications. They are found in many geological rock-forming materials and their melts. Due to their excellent chemical durability, some of these glasses are used as a medium for the immobilization of radioactive wastes. [1][2][3] Incorporating alumina in silica can provide additional resistance to the glass by enhancing its mechanical properties. 4 Sodium aluminosilicate (NAS) is a modified silicate glass widely used in the chemically strengthened glass after ion exchange phenomena as demonstrated in the famous Gorilla Glass. [5][6][7] They exhibit a great variation in physical properties such as density, molar volume, viscosity, refractive index, thermal expansion and thermodynamic properties when Na and Al contents are varied. They show a drastic change at the composition ratio R close to 1. [8][9][10][11][12] At this composition, viscosity shows a maximum while activation energy for electrical conductivity goes to a minimum. 9,11,13 The causes for this anomaly are mainly explained on the basis of geometrical description of the changes in the coordination number of Al, 12,14,15 the existence of tri-cluster O 16,17 and the degree of network polymerization. [18][19][20][21] Detailed investigation to extract the structural features and internal bonding associated with this peculiar behavior is still lacking. A precise atomic-scale understanding of the structural and physiochemical properties of such multi-component glasses is extremely important for their specific applications.
The binary aluminosilicate glasses are generally composed of well-coordinated Si, and to a lesser extent Al framework connected by the bridging oxygen, whereas the alkali-doped silicate glasses contain depolymerized network due to the presence of non-bridging oxygen (NBO). When both Na and Al are present in silica, the ternary glass network or NAS becomes far more complicated. The Si and Al generally still form tetrahedral units, but the Na ions either act as a charge compensator for Al tetrahedron [AlO 4 ]or break the network to produce NBOs. The roles played by Na are usually determined by the Na/Al ratio R and its effect are discussed by dividing into three regions: peraluminous (R < 1), meta-aluminous (R=1) and peralkaline (R>1). In the R < 1 region, all Na ions are used to compensate the charge in Al tetrahedrons. Since there are not enough Na ions to neutralize the Al tetrahedrons, it leads to the formation of O tri-cluster where three tetrahedral units share one O. 13,16,22 Alternatively, Al acts as a network modifier by being in a higher coordinated state (5 or 6). 8,23,24 Some previous studies have shown the existence of both tetrahedral and octahedral coordinated Al in region R < 1. 12,14,24 Other authors proposed that only four-coordinated Al exits over the entire compositional range including R < 1 along with a formation of three coordinated bridging O. 22,25 In the R = 1 region, sufficient Na ions are present to compensate the charge in [AlO 4 ]units so none of them will produce NBO, indicating that the glass structure is fully polymerized. So, adding equal amount of Al to Na in sodium silicate glass could eliminate NBOs. This notion is supported by some experimentalists. [26][27][28] The observed extremum in physical properties at R = 1 mentioned above supports the structural model of full polymerization. However, spectroscopic studies disagree with this conclusion. The O (1s) level in XPS indicated that some [AlO 6 ] 3exist even at R = 1. 29 Other experimental studies have also shown the existence of NBOs for R = 1. 8,30 All these raise the questions about the validity of the assertion that at R = 1, all Na ions in NAS are used for charge compensation. In the R > 1 region, all Al tetrahedrons can be neutralized by Na and the surplus of Na then act as network modifier in breaking the Si-O-Si bond. In other words, for R > 1, Na acts as both a charge compensator and a network modifier.
The presence of NBO in alkali-doped aluminosilicate glass has a significant influence on their physical properties including transport and thermodynamics properties, chemical durability, hardness and glass transition temperature. This is an important issue for commercial glass products. Apart from the ideal picture that NAS glass consists of both tetrahedrally coordinated Si and Al linked by bridging O and Na behaves as a charge compensator and network modifier, the existence of higher coordinated Al and NBOs in real glasses makes it an extremely complicated and challenging problem. A comprehensive and accurate study of NAS glass covering a wide range of Na/Al ratio is essential to elucidate the short-and intermediate-range order in these complex glasses.
The backbone arrangement of the network former ions in silicate glasses is the fundamental concept central to their structure and properties. Addition of Al and Na to silicate glass perturbs the local structures in both short and intermediate ranges leading to changes in glass properties. Experimental techniques such as X-ray diffraction, 19 x-ray photoelectron spectroscopy, 31,32 nuclear magnetic resonance (NMR), 28,33-40 X-ray radial distribution, 41 extended X-ray absorption fine structure (EXAFS), 20,21 Raman spectroscopy 18 etc. carried out over the entire range of R have provided a wealth of structural information of aluminosilicate glasses. Despite numerous investigations, the exact interaction of Na with Al tetrahedrons in the silica network is still under debate due to lack of precise information on the local atomic-scale structure and their distribution. In this regard, computer simulation is a powerful alternative to extract the structural details, some of which are not accessible via experimental probes. Calculation of physical properties based on the simulated structure can be related to the measured macroscopic properties. Comparison between these two different approaches can validate the simulation results on one side and provide new insights and directions for experiment on the other sides. This creates a win-win situation. Most of the simulations on glasses use classical molecular dynamics (MD) which suffer from a lack of accurate interatomic potential when the system becomes multicomponent and complex such as in NAS. Clearly, this drawback makes it very difficult for classical MD to describe bond formation, bond dissociation and chemical reaction occurring at specific sites. Hence, only a few classical MD simulations on NAS glasses exist. 22,[42][43][44][45][46] Ab initio MD based on quantum mechanics is the best alternative to meet the required high degree of accuracy. The price paid is that they are computationally demanding and are limited to smaller models and shorter simulation time.
In this paper we extract the details of structural, physical and electronic properties of NAS glasses. The goal is to analyze the evolution of glass properties and provide new insights when the composition changes with varying Na/Al ratio. As a first step, our simulation for NAS glasses are limited to bulk models which can provide the basic information and necessary correlations that can be applied to more realistic surface models with solvation effect. This paper is organized as follows: after the introduction section, we describe the methodology used for model construction and properties calculation in Section II. Section III presents the results from our simulation study and discussion including a sub-section III E on the reverse simulation with the removal of all Na and Al oxides from the NAS models. In section IV, we discuss briefly on the solvation and hydrolysis in the NAS model for R = 1. In the last section, we give a brief conclusion and our plan in extending the simulation to surface models of solvated alkali and mixed alkali aluminosilicate glass.

A. Simulation of glass models
Ab initio molecular dynamics (AIMD) based on density functional theory (DFT) is used to simulate the glass structures. This avoids the use of empirical interatomic potential used in classical MD which suffers from the lack of sufficient accuracy and transferability for complex multicomponent glasses. 47 AIMD provides more accurate results and reliable predictability. We have simulated four NAS models (SiO 2 ) 0.6 (Al 2 O 3 ) 0.4-x (Na 2 O)x for x = 0, 0.135, 0.20 and 0.24 corresponding to ratio R of Na/Al = 0.0, 0.5, 1.0 and 1.5 respectively. The AIMD simulations of these models with size of about 700 atoms each are carried out using the Vienna ab initio simulation package (VASP). 48,49 We used a generalized gradient approximation (GGA) PBE potential 50 for exchange correlation functional. The interaction between the ions and electrons is described by projector augmented wave (PAW) method. 51,52 All simulations are carried out under NVT ensemble where the volume of the simulation box, temperature of the system and number of particles remain constant. The NVT ensemble has the advantage that the glass structure with desired density can be simulated.
A Nose thermostat 53,54 is used to control the heat bath temperature. A melt and quench technique is adopted to mimic the thermal process in experiment. The initial model for the glass is generated from the random distribution of atoms in a cubic box with periodic boundary conditions. The cell parameters are ARTICLE scitation.org/journal/adv determined based on the available experimental density. 55, 56 We used a sufficiently high cutoff energy of 600 eV for the plane wave expansion with the convergence limits for ground state energy and force set at 10 -5 eV and 10 -3 eV/Å respectively. A single Γ-point sampling is used for Brillouin zone. The Newton's equations of motion are integrated under Verlet's algorithm with ionic time step of 1 fs. At each MD step, the solutions of instantaneous Kohn-Sham equation provide the wave-function for the residual minimization scheme with direct inversion in the iterative subspace (RMM-DIIS) method. 57 This method is very efficient for diagonalizing the Kohn-Sham Hamiltonian for larger systems such as in our models. All models are first heated from the initial random configuration to a hypothetical temperature of 4000 K, high enough to liquefy the glass within 2 ps in the canonical ensemble. The systems are held at this elevated temperature for 4 ps to equilibrate the melt and erase any memory effect from previous configuration. Next, these models follow a series of discrete reduction in temperature in step of 500 K and equilibration at each temperature until they reached at 500 K and finally to 300 K with a total simulation time of about 80 ps. The models are simulated for 4 ps at each temperature for the additional relaxation of the system to ensure no correlation with the structure from the previous temperature. This relatively short time is not meant for investigation of liquid dynamics but is reasonable for producing satisfactory glass models. The near linear diffusion region after 3 ps in the mean square displacement curve for heaviest ions Al and Si around glass transition region in NAS-0.5 ( Figure S1 in supplementary material) supports this relaxation. Similar results were recorded for other systems. During the simulation, we monitored the temperature and pressure of system at each stage to make sure that the fluctuation was under full control in the specified time period. In the glassforming region below 2000 K, the cooling rate is slowed down to a nominal quenching rate of ∼40 K/ps. This rate is still considerably faster than the actual cooling rate in the laboratory but is quite reasonable for modern AIMD. [58][59][60] The final AIMD relaxed structure at 300 K have converged temperature with a fluctuation of ±10 K around 300 K and pressure with a fluctuation of ±7 kB around 0 kB. These structures are further optimized with no constraint on volume to release the pressure to 0 kB. The resulting atomic positions and cell parameters, total energy and density are recorded. These final optimized structures are the representative models used in our structural analysis and properties calculation. AIMD simulations of large supercell is very time-consuming. In most cases, such simulations are carried out on models of a few hundreds of atoms at most 59,61,62 and rarely exceed one thousand atoms. 63 Our four models of around 700 atoms is considered sufficiently large for the present study.

B. Calculation of physical properties
Similar to our previous studies on other silicate glasses, 64-66 we used the VASP-optimized structure to calculate the elastic properties based on the strain-stress analysis approach. 67 A small strain of 0.5 % is applied to the independent strain element in the relaxed structure at a fixed volume. The elastic coefficients Cij are obtained by solving a set of linear stress vs strain equations. From the values of Cij, we can obtain the bulk mechanical parameters: Young's modulus E, bulk modulus K, shear modulus G and Poisson's ratio η using Voight-Reuss-Hill scheme for polycrystalline approximations. [68][69][70] (See supplementary material).
The electronic structure and bonding are calculated using the first-principles orthogonalized linear combination of atomic orbitals (OLCAO) method. 71 The OLCAO is an all-electron method based on local density approximation (LDA) of DFT using atomic orbitals for basis expansion. This method is extremely efficient for electronic structure calculation for large and complex systems such as inorganic, 64,72-74 organic 75 and biomolecular 76,77 materials containing a large number of atoms in the simulation cell. The centerpiece is the evaluation of bond order (BO) values ρ αβ between all pairs of atoms and the total bond order density (TBOD) as well as the effective charge Q * on each atom based on the Mulliken scheme (See supplementary material). The optical properties of the four NAS models are calculated in the form of frequency-dependent complex dielec- . The imaginary part ε 2 ( ̵ hω) is obtained first and the real part ε 1 ( ̵ hω) is obtained from ε 2 ( ̵ hω) using Kramer-Kronig transformation. 78 Finally, the square root of the real part ε 1 ( ̵ hω) in the zero-frequency limit gives the refractive index n = √ ε 1 (0).

III. RESULTS AND DISCUSSIONS
A. Structural analysis of simulated models The four optimized NAS glass models are sketched in Figure 1 showing polyhedral network formed by Si and Al atoms, and distribution of Na atoms along with bridging, non-bridging and tricluster O. Figure 1(a) shows the structure for binary aluminosilicate glass with no Na whereas Figures 1(b), 1(c) and 1(d) are for ternary systems with R = 0.5, 1.0 and 1.5 respectively. The NAS-1 model is highly polymerized while NAS-1.5 has a more disruptive network due to the presence of a large number of NBO. These figures are just 2D sketches and comprehensive analysis should be discussed in terms of pair distribution function, density, coordination number (CN), bond length (BL) and bond angles (BA) distributions etc. Table I lists the composition, structural parameters, volume and density of the four NAS models.
The macroscopic description of atoms in glass structures is embodied in the pair distribution function (PDF) which reveals the disorderness and short-range order. The calculated partial PDF (PPDF) or simply PDF for all possible atomic pairs in the four simulated NAS models are shown in Figure 2. The strongest first peak from the Si-O pairs (1.61-1.64 Å) merges with the peak from Al-O pairs (1.75-1.77 Å) resulting an asymmetric first peak. The average tetrahedral distance, (Si, Al)-O for NAS-0 is found to be 1.69 Å. Okunu et al. in an X-ray diffraction, IR and Raman studies of aluminosilicate glass with 38 mol. % of alumina reported this distance to be 1.69 Å. 79 Figure 2 shows the intensity of peak changes when R varies. However, there is no    81 It is noteworthy to point out that our result on Al-O-Al BA distribution is in contrast to the rather old Lowenstein's Al avoidance principle which states that when two tetrahedrons are linked by one bridging O, only one central cation can be Al. 82 The coordination number (CN) of atoms in the NAS models provides additional useful insights on local structural environment around network former and modifier in the glass. The CN are calculated based on somewhat arbitrarily defined cutoff length usually estimated from the first minimum in the PDF. The average CN for Al, Na and Si in the four NAS models are listed in Table II. The CN for Si in all models is close to 4, or most of Si is associated with SiO 4 tetrahedron. For NAS-0, the Al coordination is slightly higher than 4 as some five-and six-fold Al exist but for the other 3 models, they are all close to 4. An X-ray, IR and Raman studies of aluminosilicate glass containing 38 mol% of Al 2 O 3 also supported the existence of 5-or 6-fold Al atoms. 79 Sen and Youngman in a NMR study of low alumina containing aluminosilicate glasses showed the presence of 4-, 5-and 6-fold coordinated Al sites. 83 For NAS-0.5 with more Al than Na, our data suggest that a few Al are in higher coordination state. This result simply shows that some of the O atoms become threefold-coordinated tri-cluster O and this changes the network topology rather than simply forming higher coordinated Al. Some of the previous results also have supported this notion. 16,22 Presence of tri-cluster O is consistent with NMR and X-ray diffraction studies of NAS glasses. 40,79 The co-existence of O tri-clusters and tetrahedral aluminum provides a much more consistent interpretation for experimental data than simply attribute to higher coordinated Al. An interesting observation is that the CN for Na listed in Table II is higher than that of Si and Al and all three ions have lowest CN for R = 1, indicating more polymerized structure for NAS-1. Figure 4 summarizes the Al coordination in histogram form for the four NAS models.

B. Electronic structure and interatomic bonding
The electronic structure and bonding in a solid are fundamental properties based on which many other physical properties can be explained. In particular, internal cohesion related to interatomic bonding is a relatively new concept to interpret mechanical properties. The calculated total density of states (DOS) for the 4 simulated glasses are shown in Figure S2. The zero of the energy is set at the top of the occupied valence band (VB). The VB DOS for all glass models consists of two segments separated by a gap of about 5 eV: a lower segment below -15 eV and an upper segment spreading within 0 to -10 eV. The lower segment DOS has a single prominent peak in NAS-0 glass while in other models it consists of two peaks, a small peak or a shoulder along with the prominent peak. The upper segment of VB DOS has one major peak and one noticeable minor peak in NAS-0 glass while in other models, they have two or more minor peaks due to the presence of Na. The conduction band (CB) DOS spectra have similar features for the four models. They also have higher number of peaks near the CB edges in sodium containing glasses than in pure aluminosilicate glass NAS-0. The band gap value (listed in Table II) decreases significantly with increasing R. The trend is shown in Figure S3(a). It is well known that DFT calculation generally underestimates the band gap. The real band gap could be larger, but the same trend is expected.
The BO and effective charge obtained from the quantum mechanical wave functions are central to describe the bonding characteristics of atoms in glasses. The BO quantifies the relative strength of the bond between a pair of atoms. In general, BO values scale inversely with BL. However, BO can be greatly influenced by the nearby atoms and the local environment of bonded atoms. It is superior to using the BL to quantify the strength of bonding as practiced in classical calculations. This point is particularly important for the alkali-doped aluminosilicate glasses because of the complicated and controversial interpretation of the structures in the complex multicomponent glasses. Figure 5 shows the scattered plot of the BO vs BL distributions in the four NAS models. The Si-O is the strongest bond with a high BO value. They follow the order of Si-O > Al-O > Na-O. Therefore Al-O is the second strongest bond pair. In NAS-1 and NAS-1.5 models, they have relatively fewer Al atoms having higher coordination and all other Al exist in tetrahedrons with charge neutralized by Na, hence the Al-O BOs are clustered with less spread than in NAS-0 and NAS-0.5 glasses. The pie chart in the inset shows the contribution in bonding from different pairs. In NAS-0 with no Na ions, the contribution due to Al-O pairs is slightly higher than from Si-O pairs even though Al-O are weaker than Si-O bonds. In Na containing models, the contribution from Si-O bonds goes up while those from Al-O goes down. The contribution from the weaker Na-O bonds cannot compete with the steady increase in the contribution from strong Si-O bonds. The sum of all BO values gives the total bond order (TBO) and when normalized with cell volume gives the total bond order density (TBOD) which is presented in Table II along with the pairwise resolved partial BOD (PBOD). TBOD is a single quantum mechanical metric to assess the internal cohesion in a material. 84 Figure S3(b) shows the variation of TBOD with R. Pure aluminosilicate is found to be a strongest glass with a high TBOD value than other glasses. With the doping of Na ions, R decreases reaching a minimum at R = 1 which is consistent with the low-density value for this glass. The effective charge Q * for all atoms in the NAS models are calculated using Mulliken scheme. 85,86 Figure 6a shows the distribution of calculated partial charge (PC) for each atom in the four NAS models and their average values. The PC of an atom is the difference of effective charge Q * from a neutral charge of that atom Q 0 and expressed as ∆Q * = Q 0 -Q * . The PC provides the information on the charge transfer between atoms. In all models, the PC for Si, Al and Na are positive while O is negative, this implies Si, Al and Na all loose charges to O. The average PC for Si, Al and Na decreases slightly while it increases for O when moving from R = 0 to R = 1.5. The PC for O has a wider range of spread due to the presence of different type of O: bridging, non-bridging and tri-cluster O. Figure 6b shows the histogram distribution of the PC for each atom in the four NAS models. For R = 1.5, Si has a wider spread compared to other three models. This reflects the different connection of Si to bonded O and NBOs. In spite of the existence of different types of O, its PC distribution is almost similar and exhibits a Gaussian type distribution mainly because O is the only negatively charged ion in the glass.

C. Mechanical properties
The mechanical properties of a glass are central to the quality of commercial glass products. A thorough understanding of mechanical properties in relation to its composition and structure is important for its technological advancement. Table III lists   parameters with R is shown in Figure 7. Despite the isotropic nature of amorphous materials in general, a small degree of anisotropy is seen in the elastic coefficients since the calculation is limited to the cell of finite size. The elastic constants decrease with increasing Na/Al ratio as shown in Figure S4. The axial elastic constants C11, C22, C33 show some minor variations with R. They decrease as R increases since addition of Na generally makes the glass softer. The Shear elastic constants C55 and C66 show virtually no variation with R. The minor difference can be accounted for as statistical noises between models of finite size. These elastic constants are mainly related to the short-range order of the glass and they dictate the inter-atomic forces and packing density. 87 Figure 7 Figure 7a shows our calculated values of E and G for NAS-0 are underestimated while the value of K match perfectly with experimental data. It also shows that for NAS-1 glass, our calculated value of E is slightly overestimated than experimental result. The Young's modulus (E), bulk modulus (K) and shear modulus (G) in the series decrease with increasing R indicating the glass becomes softer. Similar results were reported by Livshits et al. 88 in an experimental and by Xiang et al. 44 in a theoretical study of NAS glasses of slightly different composition. The rate of decrease is faster from NAS-0 to NAS-0.5 and then decreases more or less linearly to NAS-1.5. The decrease in mechanical properties with increasing Na content is consistent with the production of NBO and de-polymerization of the network. Figure 7(b) shows the Poisson's ratio (η) and the Pugh's modulus ratio (G/K) which are inverse to each other. Based on the Pugh's definition of critical value separating the ductility and brittleness of materials 89 NAS-0 glass having G/K < 0.57 is more ductile whereas all other three glasses having G/K > 0.57 or more brittle. This is somewhat surprising since it indicates that the addition of alkali ions may make the glass more brittle. A possible strategy for making it more ductile is the use of mixed alkali effect where a larger alkali ion say K is mixed with Na that may ameliorate the brittleness problem especially when the solvent effect is also considered. We will briefly discuss such strategies in the later section.

D. Optical properties
Refractive index of glasses always dominates their optical properties in all aspects of their applications. We have calculated the interband optical properties in the form of frequencydependent complex dielectric functions and the results are presented in Figure 8(a) along with the energy loss function (ELF) in Figure 8 after overcoming the threshold energy determined by the band gap. They exhibit a broader peak or a plateau from 9 to 13 eV as indicated by the two vertical arrows. They remain relatively constant thereafter and show large variations with R. The real part of the dielectric function ε 1 is calculated from ε 2 and the broad peak follows a similar pattern but at a lower energy. Figure 8(b) shows the ELF calculated from the inverse of the complex dielectric function. They all have very smooth spectra with a well-defined peak, Plasmon peak, with frequency ωp at 20.92 eV for R = 0 but shifts slightly to the lower energy as R increases. The Plasmon frequency in a solid is the frequency of collective excitation of electrons and can be easily characterized experimentally. Table II lists the refractive indices for the simulated glasses obtained from n = √ ε 1 (0). Figure S3 (c) shows the variation of refractive index n with R value. n decreases first with increasing R, attains a minimum at R = 1 and then increases again at R = 1.5.

E. Exploration of possible silicate network recovery
In the present simulation for NAS glasses using AIMD, Na and Al oxides merge well with silica and adjust the final network topology as described in section II A above. An interesting question here is what will happen if Na and Al are removed from the well-constructed and fully optimized NAS models. Can the models return to the pure silica network and to what extent? A detailed study of this subject in aqueous environment is instrumental to the understanding of vulnerability of silica network to chemical attack when it is disrupted and can provide insights on the leaching effect and chemical durability of oxide glasses. As a first step to address this issue, we remove all the Na and Al oxide from the four NAS models and perform the AIMD simulation in the reverse direction. The starting models for this simulation, after removing Na 2 O and Al 2 O 3 from the NAS glasses have different initial configurations while the final targeted structures of pure amorphous SiO 2 (a-SiO 2 ) are the same.
After the removal of Na and Al oxides, the initial volume is much larger than the volume of pure silica glass. The AIMD simulation under NPT ensemble is then used to adjust the volume using Parrinello-Rahman technique with Langevin thermostat. 90 Each sample now contains only 360 atoms (120 SiO 2 molecule). This cell size is still sufficiently large, so a single gamma point sampling is used. After the initial shake up under NPT ensemble for volume adjustment, the simulation is then followed by NVT ensemble until the temperature and pressure fluctuation is under good control. The final AIMD structures of the pure silica glass at room temperature are further optimized without volume constraint similar to the construction of the NAS models. The same OLCAO method is now applied to analyze the recovered a-SiO 2 models. Figure S5 shows the scattered plots of BO vs BL for the four recovered silica glasses. The figure in the inset on the top right panel shows the BO vs BL plot for a separate independently simulated a-SiO 2 glass (384 atoms) using AIMD. Comparison of this plot for pure a-SiO 2 with the four silica glass models obtained from NAS models shows they are almost identical except its TBOD value is slightly higher than those four recovered a-SiO 2 glass models. This clearly indicates that after the removal of Na and Al from the NAS models, they can recover to pure a-SiO 2 structure. Figures S6 and S7 show the TDOS and dielectric functions for the four recovered a-SO 2 models. They show no discernable differences.

IV. HYDROLYSIS OF NAS BULK GLASS
Glass-water interaction is an important issue for long-term durability of silicate glass. Its properties can be greatly altered due to corrosive effect by water and ions contained therein either in the ARTICLE scitation.org/journal/adv interior or on the surface of the glass. When water reacts with pure silica or alkali silicate glasses, the network is depolymerized due to formation of hydroxyl group OH from dissociated water. It breaks the Si-O-Si bond forming the silanol group Si-OH. [91][92][93] In the case of alkali aluminosilicate glass this reaction is much more complex and different water dissolution mechanism have been proposed. Earlier studies on hydrous NAS glasses have shown dissociated water causes the de-polymerization of the glass network by the rupture of T-O-T framework (T = Si, Al). 94,95 In contrast, Kohn et al. in a NMR study argued that water does not depolymerize the aluminosilicate framework, rather Si-O-Al linkage gets protonated forming the bridging hydroxyl Si-OH-Al and Na-OH complexes. [96][97][98] Many studies have criticized this mechanism and supported the depolymerization model along with the existence of Si-OH and Al-OH groups. [99][100][101][102][103] Many experimental and theoretical works have been conducted to explore the reaction mechanism between alkali aluminosilicate glass and water. 101,104-108 Despite much effort, a definitive conclusion and well-accepted model for hydrated alkali aluminosilicate glass is still lacking. Rigorous atomistic calculation with information on interatomic bonding, distribution of different cations, un-dissociated water molecules and hydroxyl group is vital to understand this complex water dissolution mechanism. This motivated us to conduct a preliminary study presented in this subsection.
We choose the NAS-1 model with R = 1 for this study. The first step is to determine the appropriate number of water molecules to be included in the simulation cell. We constructed six glass models from NAS-1 with addition of 6, 8, 10, 12, 15 and 20 water molecules. The initial volume of the cell is expanded by about 5-8 % and then H 2 O molecules are placed at interstitial positions in a random fashion making sure that the water molecules do not overlap with the other atoms. These initial models are first relaxed with atomic positions fixed to adjust the volume. The model is then fully relaxed without any constraints. The final six optimized solvated models are then used as input for OLCAO calculation. The interatomic bonding and TBOD for each model with different water contents are analyzed. It turns out that the model with 12 H 2 O molecules has the highest TBOD (see Table S1) and is chosen as the representative model for the solvated NAS glass. Figure 9 shows the BO vs BL plot for hydrated NAS-1 glass showing the complex interplay of different types of bonding including hydrogen bonding (HB). It is clear that water dissolution in NAS glass leads to a complex bonding pattern with different BO values. The dissociated water breaks the glass framework in forming the strong covalent Si-OH and Al-OH bonds having a BL of ∼1 Å with high BO values. The intramolecular covalent O-H bonds in H 2 O have similar BL and bonding strength. The weak HB (O⋯H) between different hydroxyl group and H 2 O molecules cover a wide range having BL ∼1.5 -2.5 Å. The strength of these HBs are less than covalent Si-O and Al-O bonds but are comparable to Na-O ionic bonds. Similar to the un-hydrated model, the contribution to bonding follows the order of Si-O > Al-O > Na-O > O⋯H. A slight decrease in TBOD value is observed in the hydrated NAS-1 compared to dry model. Bartholomew et al. in their infrared spectroscopy study showed water in glass exist as molecular water H 2 O and dissociated water in the form of OH group. 109 Stuke and coworkers further showed that increasing the water content tends to reduce the dissociation. 110 Our result also shows the presence of both H 2 O and OH with later as the dominating group.
Our preliminary exploration on hydrolysis of NAS glass is aimed at providing the atomistic details based on rigorous DFT calculations. This can facilitate better and deeper interpretation of the experimentally measured data in many laboratories. The glasswater interaction with mixed alkali ions can be even more complex and intriguing. Fundamental understanding of this reaction is of paramount importance in ion exchange and chemically strengthening of many novel inorganic glasses. The present work on this hydrated NAS-1 glass is a first step for such studies.

V. CONCLUSIONS
In this work, we have systematically investigated the network structure, electronic, mechanical and optical properties of sodium aluminosilicate glasses as a function of Na/Al ratio R. Detailed analysis on changes in structure and properties are carried out in composition regions of R < 1, R = 1, and R > 1 in comparison with experiments. The use of accurate AIMD simulation enables us to clarify many confusing and contradictory issues regarding the structure and properties of this complex glass system. Both the mechanical and optical properties of these NAS glasses are calculated and discussed in comparison with limited existing data so far.
We find that in NAS glasses, Al prefers tetrahedral network connection sharing the tri-cluster O rather than being in a higher coordinated (5 or 6) state. These results clarify many conflicting reports in the literature and refute the longstanding Lowenstein's rule of Al avoidance principle. A novel concept of TBOD is used to describe the complex interplay of different types of bonding and to characterize the internal cohesion on the simulated glasses.
In the AIMD simulation carried in reverse direction by removing all Na and Al from NAS glasses, we found that silica network is ARTICLE scitation.org/journal/adv recovered again but with a weaker internal cohesion. We have also attempted a preliminary study on the solvation effect and hydrolysis in NAS glass with R = 1. The results are very encouraging and revealing and set the stage for similar investigations on ion exchange and chemically strengthening in many novel inorganic glasses.

SUPPLEMENTARY MATERIAL
The supplementary material file contains description of: 1-Computational Methods. 2-Additional results in the form of Tables and Figures. 3-References.