Zhen Zhang

Simona Ispas
email: simona.ispas@umontpellier.fr

Walter Kob
email: walter.kob@umontpellier.fr

The critical role of the interaction potential and simulation protocol for the structural and mechanical properties of sodosilicate glasses

Keywords: silica, sodosilicate glass, computer simulation, structure, mechanical properties, finite size effects

We compare the ability of various interaction potentials to predict the structural and mechanical properties of silica and sodium silicate glasses. While most structural quantities show a relatively mild dependence on the potential used, the mechanical properties such as the failure stress and strain as well as the elastic moduli depend very strongly on the potential, once finite size effects have been taken into account. We find that to avoid such finite size effects, samples of at least 75,000 atoms are needed. Finally we probe how the simulation ensemble influences the fracture properties of the glasses and conclude that fracture simulations should be carried out in the constant pressure ensemble.

Introduction

Computer simulations have by now been established to be a highly valuable tool to get insight into the properties of complex systems such as liquids and glasses [START_REF] Kob | Computer simulations of supercooled liquids and glasses[END_REF][START_REF] Massobrio | Molecular dynamics simulations of disordered materials[END_REF][START_REF] Kob | First-principles simulations of glass-formers[END_REF]. The most important ingredient for such simulations is the interaction potential between the particles since its accuracy is crucial for obtaining a reliable description of the properties of the material on the micro-and macroscopic scale. There are two possibilities to describe these interactions: ab initio or using effective potentials (also called "classical simulations"). In the first approach the interactions are calculated directly from the electronic structure of the system, a procedure that is very accurate (although not exact) but computationally very expensive [START_REF] Marx | Ab initio molecular dynamics: basic theory and advanced methods[END_REF]. As a consequence only relatively small systems can be simulated, i.e., typically less than 10 3 particles, over a short time (less than 10 ns). In contrast to this, effective potentials allow to access significantly larger systems and longer times: 10 7 atoms and tens of micro-seconds. This advantage comes at the cost of a deteriorated accuracy of the potentials and hence the predicted material properties are less reliable. Despite this drawback many numerical investigations are done with effective potentials since for many studies one needs to have systems that are relatively large, e.g. to probe mechanical properties, and cooling rates that are not too high. As a consequence significant efforts have been taken to obtain effective potentials that are reliable and for systems like silica one can find in the literature dozens of potentials, see Ref. [START_REF] Sundararaman | New optimization scheme to obtain interaction potentials for oxide glasses[END_REF] for an exhaustive list. (In recent times force fields are sometimes obtained from approaches called "machine learning" [START_REF] Behler | Perspective: Machine learning potentials for atomistic simulations[END_REF][START_REF] Liu | Machine learning for glass science and engineering: A review[END_REF]. Since at present it is not clear to what extent these approaches give transferable potentials we will here not discuss them further.)

It is obviously most important to know which ones of these potentials are reliable and which ones are not. Therefore one can find several studies in the literature in which the performance of various potentials are compared [START_REF] Hemmatti | IR absorption of silicate glasses studied by ion dynamics computer simulation. I. IR spectra of SiO 2 glass in the rigid ion model approximation[END_REF][START_REF] Hemmatti | Comparison of pair potential models for the simulation of liquid SiO 2 : Thermodynamic, angular-distribution, and diffusional properties[END_REF][START_REF] Soules | Silica molecular dynamic force fields -A practical assessment[END_REF][START_REF] Cowen | On force fields for molecular dynamics simulations of crystalline silica[END_REF][START_REF] Yuan | Local structures of MD-modeled vitreous silica and sodium silicate glasses[END_REF][START_REF] Bauchy | Structural, vibrational, and elastic properties of a calcium aluminosilicate glass from molecular dynamics simulations: The role of the potential[END_REF]. Roughly speaking one finds that the structural properties, such as the static structure factor, are relatively independent of the potential considered, a result that is not surprising since the parameters of many potentials have been optimized to reproduce the glass structure determined in experiments. Dynamical quantities, such as the diffusion constant, show a much stronger dependence on the potential used [START_REF] Hemmatti | Comparison of pair potential models for the simulation of liquid SiO 2 : Thermodynamic, angular-distribution, and diffusional properties[END_REF]. This result is related to the fact that dynamics depends not only on the interaction energy close to the local minima of the potential energy landscape, which governs the structural properties of the glass, but also on the local barriers, i.e., quantities that are usually not taken into account when the parameters of a potential are optimized. As an example we mention here the work by Hemmatti and Angell who probed how the diffusion constant of Si in SiO 2 depends on the potential and found that the predicted values varied by orders of magnitude [START_REF] Hemmatti | Comparison of pair potential models for the simulation of liquid SiO 2 : Thermodynamic, angular-distribution, and diffusional properties[END_REF].

As mentioned above, the investigation of mechanical properties of glasses are one of the important applications for simulations with effective potentials. Whether the potentials that are currently used to simulate such systems are indeed reliable to probe these properties has, however, not been tested so far. The goal of the present work is therefore to fill this gap by investigating the fracture behavior of glasses. In parallel we will also study how the ensemble used for the simulation, constant pressure or constant volume, influences the fracture.

For the present study we focus on the case of pure silica and the sodo-silicate system Na 2 O-3SiO 2 , in the following denoted as NS3, since these are two glass-forming systems that are representatives for many oxide-glasses. The paper is organized as follows: In Sec. 2 we give the details of the investigated potentials and the simulation procedure. In Sec. 3 we present the results on the static, dynamical and mechanical properties of the various glass samples. The paper is concluded with a summary and a discussion in the last section.

Simulation details 2.1 Interatomic potentials

For the simulation of SiO 2 and Na 2 O-3SiO 2 (NS3) we considered four pair potentials, all of which have the same functional form given by

V (r ij ) = q i q j e 2 4π 0 r ij + A ij e -r ij /B ij - C ij r 6 ij , (1)
i.e., they are given by the sum of a Coulomb and Buckingham term. (Here r ij is the distance between atoms i and j.) Thus the difference between the potentials are just the values for the various parameters of the potential, i.e., q i , A ij etc. In particular we will consider the potentials proposed by Habasaki and Okada (HO) [START_REF] Habasaki | Molecular dynamics simulation of alkali silicates based on the quantum mechanical potential surfaces[END_REF], by Teter et al. [START_REF] Cormack | Alkali ion migration mechanisms in silicate glasses probed by molecular dynamics simulations[END_REF], by Guillot and Sator (GS) [START_REF] Guillot | A computer simulation study of natural silicate melts. Part I: Low pressure properties[END_REF],

and by Sundararaman et al. (SHIK) [START_REF] Sundararaman | New optimization scheme to obtain interaction potentials for oxide glasses[END_REF][START_REF] Sundararaman | New interaction potentials for alkali and alkaline-earth aluminosilicate glasses[END_REF]. For the case of silica we also considered the potential by van Beest et al. (BKS) [START_REF] Van Beest | Force fields for silicas and aluminophosphates based on ab initio calculations[END_REF] since it has been used in many previous studies and found to be able to reproduce quite well many properties of real silica [START_REF] Vollmayr | Cooling-rate effects in amorphous silica: A computersimulation study[END_REF][START_REF] Horbach | Static and dynamic properties of a viscous silica melt[END_REF][START_REF] Luo | Molecular dynamics study of correlations between IR peak position and bond parameters of silica and silicate glasses: Effects of temperature and stress[END_REF]. The parameters of the various potentials used in the present work are given in Table 1.

The presence of the Coulomb term in Eq. ( 1) makes the use of such potentials computationally expensive since they have to be evaluated by means of approaches like the Ewald summation [START_REF] Allen | Computer simulation of liquids[END_REF]. One possibility to avoid this problem is to use the method proposed by Wolf et al in which the Coulomb term is replaced by

q i q j e 2 4π 0 r ij →      q i q j e 2 4π 0 1 r ij - 1 r c + r ij -r c r 2 c r < r c 0 r ≥ r c (2)
where r c is a cutoff distance [START_REF] Wolf | Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r -1 summation[END_REF][START_REF] Carré | Amorphous silica modeled with truncated and screened Coulomb interactions: A molecular dynamics simulation study[END_REF]. In this form the potential becomes thus short ranged and hence computationally much more efficient. For the case of the SHIK potential we have compared the structural, dynamical, and mechanical properties of the glass if the Coulomb term is replaced by the expression given by Eq. ( 2) and in the following we will show that there is no significant difference. Hereafter we will denote the potential with Eq. ( 2) as SHIK and the one with the full Coulomb interaction as SHIKc.

Due to the van der Waals term in Eq. ( 1), the potentials have a singularity at short distances.

To prevent that particles fuse together in a unphysical manner, we added a short range repulsive term [START_REF] Sundararaman | New interaction potentials for alkali and alkaline-earth aluminosilicate glasses[END_REF]. This modification does not affect at all the properties of the system at intermediate and low temperatures and hence can be considered to be just a computational trick to avoid this singularity.

Finally we mention that the charge of the oxygen atoms for the HO and SHIK potentials depends on composition in order to maintain charge neutrality of the system when the sodium concentration is varied q O = (1 -y) q Si + 2y q Na 2 -y , where y is the Na 2 O mole concentration [START_REF] Habasaki | Molecular dynamics simulation of alkali silicates based on the quantum mechanical potential surfaces[END_REF][START_REF] Sundararaman | New interaction potentials for alkali and alkaline-earth aluminosilicate glasses[END_REF]. Thus the oxygen charges reported in Table 1 have been calculated for y = 0.25 which corresponds to the sodosilicate composition studied in the present work.

In Fig. 1 we plot the potential energy and the forces for the potentials considered. The oxygen charge used for silica when using the SHIK potential is -q Si /2. Also included are the cutoff distances R cut for the short range/long range parts of the potentials. the potential used. Note that, in order to show that we are indeed plotting the relevant range in distance, we also include the radial distribution function g αβ (r), with α, β ∈ {Si, O} (green lines with symbols) in the graphs for the force. (These g αβ (r) are for the SHIK potential, at 300 K, but the other potentials predict radial distribution functions that are very similar.)

Simulation procedure

The systems we consider are pure silica and a sodium silicate with composition Na 2 O-3SiO 2 corresponding to a Na 2 O molar concentration of 25 %. The glass samples were produced by using the conventional melt-quench method. We used cubic boxes (periodic boundary conditions) that contained between 5000 and 600,000 particles, which corresponds to sizes between 4 and 20 nm at room temperature. For the long-range Coulomb interaction, the Wolf truncation method (see Eq. ( 2)) was employed when using the SHIK potential while for the other potentials this interaction was evaluated with the particle-particle particle-mesh (PPPM) solver algorithm with an accuracy of 5 × 10 -5 . As mentioned above, we also carried out simulations for silica using the SHIK parameters for the short-range part of the potential, and Coulomb interactions evaluated using PPPM algorithm, and the corresponding data will be labeled SHIKc.

The samples were first equilibrated at a high temperature in the canonical ensemble (N V T ) using a fixed volume that corresponds to the experimental value of the density of the glass at room temperature [START_REF] Bansal | Handbook of glass properties[END_REF]. These N V T runs were done at 3600 K for silica and 3000 K for NS3, both for about 300 ps, a time that is sufficiently long to equilibrate the samples completely.

These liquids were subsequently equilibrated in the N P T ensemble (constant number of atoms, pressure, and temperature) at the same temperatures and at zero pressure. The lengths of these N P T runs depended on the potential considered and were sufficiently long to equilibrate the samples, i.e. the mean squared displacement of the particles was more than 100 Å2 (see below).

For the GS potential, the N P T equilibration of the NS3 liquid was done at 2100 K since for higher temperatures the samples became unstable because they were above the boiling point of the GS potential.

All simulations were carried out using the LAMMPS software [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] with a time step of 1.6 fs and using a Nosé-Hoover thermostat and barostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF][START_REF] Hoover | Constant-pressure equations of motion[END_REF]. After equilibration, the liquid samples were quenched to 300 K in the N P T ensemble at zero pressure. The cooling rate was 0.25 K/ps, i.e. a value that is relatively small for MD simulations of sodium silicate glasses with a comparable system size. Previous simulation studies have shown that this cooling rate is small enough so that the properties of the system do not depend on the cooling rate in a significant manner [START_REF] Vollmayr | Cooling-rate effects in amorphous silica: A computersimulation study[END_REF][START_REF] Lane | Cooling rate and stress relaxation in silica melts and glasses via microsecond molecular dynamics[END_REF][START_REF] Li | Cooling rate effects in sodium silicate glasses: Bridging the gap between molecular dynamics simulations and experiments[END_REF]. The glass samples at 300 K were then annealed in the N P T ensemble for 160 ps before we started the mechanical tests. Since room temperature is well below the glass transition temperature, there is no need to consider a longer thermalization run as the atomic configurations are practically frozen, and neither the density nor the structural properties change noticeably.

The results presented in the following sections have been obtained by using for each potential only one sample per composition. For the structural and dynamical properties discussed in Sec.

3.1, the samples contained 36,480 and 38,400 atoms for SiO 2 and NS3, respectively, corresponding to box sizes around 8 nm. These system sizes are sufficiently large to make sample-to-sample fluctuations small.

For a given composition and a given potential, the mechanical properties presented in Subsec. 3.2 are the average obtained from one sample put under uniaxial tensile strain in the three independent directions, i.e. x, y, and z directions, and the error bars correspond to standard deviation. To calculate the stress-strain curve we increased the dimension of the box in one direction linearly in time, L(t). The strain is then given by

(t) = ˙ t ( 3
)
where ˙ is the strain rate. The stress tensor is obtained from the usual expression [START_REF] Thompson | General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions[END_REF]

σ = 1 V N i=1 [m i v i ⊗ v i + r i ⊗ f i ] (4)
where V and N are the volume and the total number of atoms of the simulation box, respectively, while m i is the mass of atom i, and v i , r i and f i are the velocity, position and force vector of atom i, respectively.

3 Results

Structural and dynamical properties in the liquid and glassy states

The results in the present subsection are all for a cubic system with the size of the edge given by L = 8 nm which corresponds to N = 36480 and N = 38400 atoms for SiO 2 and NS3, respectively.

This size is sufficiently large to avoid that dynamic and static quantities are affected by finite size effects.

To investigate the dependence of the dynamics on the potential we have calculated the mean squared displacement (MSD) of a tagged particle:

r 2 α (t) = 1 N α Nα j=1 |r j (t) -r j (0)| 2 , (5)
where α ∈ {Si, O, Na} and N α is the number of particles of species α.

In Fig. 2 we show the time dependence of the MSD for Si, i.e., the species that moves the slowest. The temperatures are those at which we have equilibrated the samples in the N P T ensemble and the curves correspond to the different potentials. One sees that, at these T s, the dynamics is already somewhat glassy in that at intermediate times the MSD has a plateau [START_REF] Binder | Glassy materials and disordered solids: An introduction to their statistical mechanics[END_REF].

These graphs also show that at long times the MSD is a linear function of time, i.e., that the particles have reached the diffusive regime, indicating that the system is equilibrated.

In agreement with previous results [START_REF] Hemmatti | Comparison of pair potential models for the simulation of liquid SiO 2 : Thermodynamic, angular-distribution, and diffusional properties[END_REF], we find that the MSD at long times, and hence the diffusion constant, shows a very strong dependence on the potential, more than a factor of ten, although the structure of the glass does not vary that much (see below). Since it can be expected that the activation energy for the diffusion constant also depends on the potential considered, the diffusion constant at lower temperatures will differ even more, as previously reported in

Ref. [START_REF] Hemmatti | Comparison of pair potential models for the simulation of liquid SiO 2 : Thermodynamic, angular-distribution, and diffusional properties[END_REF].

In Fig. 3 we show the temperature dependence of the mass density as predicted by the various potentials. These curves were obtained by cooling the samples at zero pressure from high temperatures to 0 K, i.e. at high T they correspond to the equilibrium density of the system whereas at low T to the density of the glass phase. The data demonstrates that the T -dependence of the density depends significantly on the potential: Not only the absolute the temperature is 3600 K. For NS3 the temperature is 3000 K except for the GS curve which is for 2100 K.

values differ but also the slopes, i.e. the thermal expansion coefficients. For NS3 these slopes are larger than the ones for SiO 2 a result that matches the experimental findings [START_REF] Bansal | Handbook of glass properties[END_REF]. Also included in the graphs are the experimental values of the density at room temperature (left most green diamond) and higher. For SiO 2 , one sees that three potentials, namely Teter, HO, and GS, predict a density that is significantly too high, whereas the BKS and SHIK potentials are able to predict this quantity quite well. For NS3, the HO and SHIK potentials make an accurate prediction of the experimental density (green diamonds) whereas the ones by Teter and GS predict values that differ somewhat more from the experimental one. These graphs show thus that there are potentials which are able to give a quantitative good prediction of the density at room temperature. Note that the densities at room temperature are influenced by the cooling rate of the sample [START_REF] Vollmayr | Cooling-rate effects in amorphous silica: A computersimulation study[END_REF]. However, this dependence is relatively mild since the density is directly related to the fictive temperature at which the sample falls out of equilibrium and this temperature depends only logarithmically on the cooling rate. Thus if the discrepancy between the predicted density and the experimental value is too large, this cannot be rationalized by a too high cooling rate but must instead be considered as a flaw of the used potential.

We have also included in Fig. 3 the experimental values of the density at higher temperatures, calculated using the experimental values of the thermal expansion coefficient [START_REF] Bansal | Handbook of glass properties[END_REF]. One recognizes that, for silica, the slope of this data matches very well the one predicted by the SHIK and BKS potentials while the three other potentials predict a larger expansion coefficient. For NS3 the graph shows that the HO potential predicts a slope which is in excellent agreement with the experimental data, the SHIK and Teter potentials are in fair agreement, while the GS potential is far off from the reality. (b) NS3 Experimental data for silica are taken from Ref. [START_REF] Susman | Temperature dependence of the first sharp diffraction peak in vitreous silica[END_REF] and for NS3 from Ref. [START_REF] Pöhlmann | Structure and dynamics of hydrous silica(tes) as seen by molecular dynamics computer simulations and neutron scattering[END_REF].

We now probe how the properties of the glass depend on the potential. To start we show in Fig. 4 the static structure factor as seen in a neutron scattering experiment which is the weighted sum of the partial structure factors [START_REF] Binder | Glassy materials and disordered solids: An introduction to their statistical mechanics[END_REF]:

S tot (q) = N α=Si,O,Na N α b 2 α α,β=Si,O,Na b α b β S αβ (q) . ( 6
)
Here q is the wave-vector, b α is the neutron scattering length for species α and S αβ (q) is the partial structure factor:

S αβ (q) = f αβ N Nα j=1 N β k=1 exp(iq.(r j -r k )) , α, β = Si, O, Na, with f αβ = 1
for α = β and f αβ = 1/2 otherwise, and N the total number of atoms. From the figure we recognize that for the case of silica all considered potentials agree very well with the experimental data. This is not that surprising since in most cases the parameters of these potentials have been optimized to reproduce the structure of the glass. Qualitatively the same conclusion can be drawn for the case of NS3, panel (b). There is, however, one exception: The GS potential shows a strong increase of the signal at small q. This behavior indicates the presence of a phase separation, in this case the formation of large domains of Na atoms, and a visual inspection of the sample shows that this is indeed the case. (We note that this defect of the potential is not readily seen in the radial distribution functions and seems to have gone unnoticed so far.) A further useful quantity to characterize the structure of silicate glasses is the distribution of the bond angles. Of particular interest is the angle SiOSi since it gives information about the relative orientation of two neighboring tetrahedra. For the case of SiO 2 this distribution is presented in Fig. 5a and we recognize that there is not much dependence on the potential considered and that all of them predict a position of the maximum of the distribution that is compatible with the experimental estimate for that angle [START_REF] Malfait | 29 Si NMR spectroscopy of silica glass: T 1 relaxation and constraints on the SiOSi bond angle distribution[END_REF]. Also for NS3 the distribution of this angle is basically independent of the potential used, Fig. 5b. The only exception is the data obtained from the Teter potential which predicts a peak at significantly larger angles. The calculated mean SiOSi angles are 143.7 GS, HO and Teter potential, respectively. When compared to the experimental mean value of 141.7 • extracted from 29 Si MAS NMR measurements [START_REF] Angeli | Insight into sodium silicate glass structural organization by multinuclear NMR combined with first-principles calculations[END_REF], the SHIK data is in good agreement, while GS and HO present a reasonable one. In addition, one can compare these distributions with the ones from ab initio simulations [START_REF] Kilymis | Vibrational properties of sodosilicate glasses from first-principles calculations[END_REF][START_REF] Tilocca | Structural and electronic properties of modified sodium and soda-lime silicate glasses by Car-Parrinello molecular dynamics[END_REF] of very similar compositions, which predicted peak positions close to 140 • , i.e. values that are compatible with the one predicted by the SHIK, HO, and GS potentials, but not with the prediction by the Teter potential.

Finally we discuss two local structural quantities that probe the local environment of an atom, namely the oxygen speciation and the distribution of the Q n tetrahedral species in NS3.

For the former we have used the first minimum in the radial distribution function of the Si-O pair to determine the connectivity of a given oxygen atom and thus to determine whether the atom is free (FO), non-bridging (NB), bridging (BO), or three-fold coordinated (TBO). In Fig. 6a we show the corresponding probability for the different potentials and one recognizes that the considered force fields all give the same distribution and it agrees very well with the experimental measurements [START_REF] Nesbitt | Bridging, non-bridging and free (O2) oxygen in Na 2 O-SiO 2 glasses: An X-ray Photoelectron Spectroscopic (XPS) and Nuclear Magnetic Resonance (NMR) study[END_REF]. Thus one can conclude that this quantity is very robust or in other terms the oxygen speciation is not very useful indicator for testing the quality of a potential. Things are different for the Q n species, i.e. the probability that a silicon atom is connected to exactly n bridging oxygens, n = 0, 1, 2, 3, 4. These probabilities are shown in Fig. 6b and one finds that, e.g., the frequency of the Q 3 structure depends significantly on the potential. In particular one observes that the GS potential gives a probability that is rather low with respect to the other potentials and also to the experimental data. Thus the observation that the fraction of Q 3 units depends significantly on the used potential indicates that this quantity can be used as an indicator for evaluating the quality of a potential. In addition we mention that the network depolymerization depends also on the cooling rate: For sodium silicate glasses, it has been shown that the percentage of Q 3 increases with decreasing quench rate while the Q 4 one is decreasing, and hence improve the agreement with the experimental data [START_REF] Li | Cooling rate effects in sodium silicate glasses: Bridging the gap between molecular dynamics simulations and experiments[END_REF].

Mechanical properties of the glass

To determine the mechanical properties of the annealed glass samples we strained them in one direction using a strain rate of 0.5/ns. As we will show below this value is sufficiently small to obtain results that do not depend in a significant manner on the rate. We considered two ensembles for this transformation: A constant pressure ensemble (N P T ) in which the pressures in the directions orthogonal to the strain were set to zero and a constant volume ensemble (N V T ) in which the two orthogonal directions were not allowed to change, i.e. the cross section was constant. From this transformation we can thus directly obtain the stress-strain curve in the pulling direction discussed in the following.

System size and cooling rate dependence

Before starting the discussion of the physical properties of the systems, we probe how these properties depend on the system size. Usually the properties of glasses can be obtained from simulations with relatively moderate system size, say 10 3 -10 4 atoms. However, elasticity and fracture are associated to non-local processes and hence finite size effects can be important. That this is indeed the case is demonstrated in Fig. 7 which shows the stress-strain curve for different (cubic) box sizes L. (These results are for the SHIK potential, but for the other force fields a similar behavior has been found [START_REF] Zhang | Fracture of oxide glasses: New insights from atomistic simulations[END_REF].) Error bars have been estimated by considering three different samples. For the case of silica, panel (a), one finds that the elastic regime is basically (b) NS3

Figure 7: Dependence of stress-strain curve on the system size for the SHIK potential at 300 K. (a) Silica and (b) NS3. The system size ranges from around 5000 to around 600,000 atoms, corresponding to the increase of the size of the simulation box size from 4 nm to 20 nm at room temperature density (at zero strain). The error bars were estimated from the standard deviation for the three samples.

independent of L. However, once the failure point (i.e. the strain at which the stress has a maximum) has been attained there are very strong finite size effects in that the small systems break in a much gentler manner than the large ones. Only if L has reached 10 nm, which corresponds to around 75,000 atoms, the stress-strain curve becomes basically independent of the system size (see Refs. [START_REF] Yuan | Molecular dynamics simulation of amorphous silica under uniaxial tension: From bulk to nanowire[END_REF] and [START_REF] Pedone | Dynamics of Fracture in Silica and Soda-Silicate Glasses: From Bulk Materials to Nanowires[END_REF] for related studies).

For the case of NS3 the system size dependence is more pronounced in that one has to use systems of about L = 16 nm, corresponding to about 300,000 atoms, before the stress-strain curve converges, see panel (b) in Fig. 7. These stronger finite size effects are likely related to the presence of the Na atoms which make a more heterogeneous structure in the NS3 glass compared to that of SiO 2 [START_REF] Greaves | EXAFS and the structure of glass[END_REF][START_REF] Horbach | Dynamics of sodium in sodium disilicate: Channel relaxation and sodium diffusion[END_REF]. However, for strains smaller than the failure point the curves for the different system sizes superimpose quite nicely and therefore also for this composition a system size with around 36,000 atoms, i.e. box size 8 nm, is sufficient to study the elastic regime.

In order to make a fair comparison of the behavior predicted by the different potentials we have used the same system size L = 8 nm, which corresponds to N = 36480 and N = 38400 atoms for SiO 2 and NS3, respectively. Although for this system size one still can observe finite size effects, they are minor and hence do not preclude to understand which potentials give rise to a realistic fracture behavior and which ones do not. Another important parameter for the properties of a glass is the cooling rate with which the sample has been produced. In the past it has been found that structural quantities show a clear dependence on this cooling rate, often described by a logarithmic law [START_REF] Vollmayr | Cooling-rate effects in amorphous silica: A computersimulation study[END_REF][START_REF] Lane | Cooling rate and stress relaxation in silica melts and glasses via microsecond molecular dynamics[END_REF][START_REF] Li | Cooling rate effects in sodium silicate glasses: Bridging the gap between molecular dynamics simulations and experiments[END_REF] and hence one can expect that also the mechanical properties depend on this rate. In Fig. 8 we show that this is indeed the case in that with decreasing rate the failure strain decreases while the failure stress increases. At the same time a decreasing cooling rate leads to an increase of the slope in the elastic regime, i.e. the system becomes stiffer. Although these cooling rate effects are noticeable, we see that once the rate is below 0.3 K/s the effect is minor and hence the results discussed below, obtained with a rate of 0.25 K/s, are only mildly affected by the cooling rate used to produce the glass sample. (Note that the samples used in Fig. 8 were produced using the SHIK potential and by quenching firstly the liquids in the NVT ensemble to a temperature a bit below the glass transition, and then to 300 K in the NPT ensemble [START_REF] Zhang | Fracture of oxide glasses: New insights from atomistic simulations[END_REF]. These samples contained 600,000 atoms.)

Fracture behavior and elastic moduli

We now discuss how the elastic moduli and the fracture behavior depend on the potential, the ensemble, and the way the Coulomb interaction is handled. In Fig. 9a we compare the stressstrain curves for BKS with the one for SHIK (obtained in the N P T ensemble). One sees that show the stress-strain curve and the tangent modulus for silica, respectively. The green symbols are the experimental data from Ref. [START_REF] Krause | Deviations from linearity in the dependence of elongation on force for fibers of simple glass formers and of glass optical lightguides[END_REF]. The green dashed line in panel (b) is the polynomial fit from Ref. [START_REF] Gupta | Intrinsic failure and non-linear elastic behavior of glasses[END_REF] to experimental data and its extrapolation to large strains. Panels (c) and (d): Same quantities as in panels (a) and (b) for the case of NS3.

the BKS glass is stiffer and stronger than the SHIK glass, and also breaks in a more ductile manner.

Also included in the graph is the data for the SHIKc potential with the full Coulomb interaction (blue dashed line). One sees that the results are very close to the SHIK curve (blue full line), i.e. when the Coulomb term is replaced by the expression given in Eq. ( 2), and hence one can conclude that Wolf truncation is a good approximation, a result which from the computational point of view is most advantageous. (Already for 38400 atoms one gains a factor of around 3 in CPU time.)

From the stress-strain curve one can obtain directly the tangent modulus which is defined as

E t = dσ d , (7)
where σ is the stress in the pulling direction and the strain. This quantity is shown in Fig. 9b for the various potentials. In agreement with the curves from panel (a) we find that for all strains the BKS glass has a tangent modulus that is significantly larger than the one predicted by the SHIK potential. Interestingly, however, both potentials predict the same failure strain, i.e. f ≈ 0.17, in excellent agreement with the experimental value of 0.18 [START_REF] Griffioen | Optical fiber mechanical reliability[END_REF]. For the SHIK potential we have also included the data as obtained by using the full Coulomb term and we see that the corresponding curve (dashed blue line) agrees very well with the one in which the Wolf term is used.

In addition, we have included in panels (a) and (b) of Fig. 9 the experimental data from

Refs. [START_REF] Krause | Deviations from linearity in the dependence of elongation on force for fibers of simple glass formers and of glass optical lightguides[END_REF] and [START_REF] Gupta | Intrinsic failure and non-linear elastic behavior of glasses[END_REF]. One sees that these data sets agree very well with the predictions of the SHIK potential which is also able to reproduce accurately the Young's modulus, i.e.

E = lim →0 E N P T t ( ) , (8)
with a value equal to 72 GPa, in very good agreement with the experimental value given by 73 GPa [START_REF] Bansal | Handbook of glass properties[END_REF]. The prediction of the BKS potential is a bit less satisfactory, 86 GPa, in agreement with the findings from Ref. [START_REF] Yuan | Molecular dynamics simulation of amorphous silica under uniaxial tension: From bulk to nanowire[END_REF]. For the failure strength the BKS and SHIK potentials predict 17.6 GPa and 12.8 GPa, respectively, and the latter value is in excellent agreement with the experimental value of 12.6 GPa [START_REF] Griffioen | Optical fiber mechanical reliability[END_REF].

The stress-strain curve of NS3 as predicted by the different potentials is shown in Fig. 9c).

Surprisingly we find that under N P T conditions the HO and Teter potentials show no sign of fracture even if the sample is stretched to 100%. This shows that these potentials have a serious flaw in that they strongly overestimate the ductile behavior of NS3. Also the stressstrain curve from the GS potential is not realistic in that the stiffness in the elastic regime is strongly underestimated, the failure stress is too small, and that the fracture is way too ductile.

Hence we conclude that these three potentials should not be used to study the fracture behavior of NS3. A much more reasonable stress-strain curves is found for the SHIK potential which shows a relatively brittle fracture. This brittleness is, however, less pronounced than the one found in SiO 2 , see panel (a), in agreement with the expectation that the addition of Na will make the glass more ductile. The failure stress for the SHIK potential is around 6 GPa (see Table 2 for exact values), which compares reasonably well with the experimental values that are between 7.5 GPa and 11.7 GPa [START_REF] Lower | Inert failure strain studies of sodium silicate glass fibers[END_REF][START_REF] Kurkjian | Intrinsic strength and the structure of glass[END_REF] (this latter value is only an upper limit, see Ref. [START_REF] Lower | Inert failure strain studies of sodium silicate glass fibers[END_REF]).

The Young's modulus of the NS3 glass, given by Eq. ( 8), is around 55 GPa for the SHIK, HO and Teter data (see Table 2 for exact values), in good agreement with the experimental value of 59.8 GPa [START_REF] Bansal | Handbook of glass properties[END_REF], while the GS potential predicts a value of only 26 GPa, thus way too small.

A more notable difference between the various potentials is found for the strain dependence of the tangent modulus, shown in Fig. 9d: While for the HO and Teter potentials, E t decreases basically in a linear manner, the SHIK and GS potentials show at intermediate strain a plateau before they decrease to zero. (Note that the data for the GS potential has been multiplied by a factor of 2.0 in order to bring it on the same scale as the other curves.) This plateau is also directly visible in the stress-strain curves in that they show at around = 0.07 a marked bend (see Fig. 9c). Elsewhere we will show that this rapid change in the effective stiffness of the sample is related to a change in the plastic behavior of the sample on the atomic scale [START_REF] Zhang | Fracture of oxide glasses: New insights from atomistic simulations[END_REF].

The fact that neither the Teter nor the HO potential predict breaking if strained to 100% (see Fig. 9c), is astonishing since, from a structural point of view, these potentials give reasonable predictions. Thus one might wonder whether this unrealistic behavior is related to the ensemble used during the tensile loading. In Fig. 10 we thus show the stress-strain curve for the two potentials but now in the N V T ensemble, i.e. the sample size orthogonal to the pulling direction is not allowed to change. We see that the so obtained stress-strain curves are completely different from the ones in the N P T ensemble in that the formers indicate that the sample breaks in a rather brittle manner. A similar ensemble effect has also been observed for the case of silica [START_REF] Pedone | Molecular Dynamics Studies of Stress-Strain Behavior of Silica Glass under a Tensile Load[END_REF].

Furthermore we note that in the elastic regime the curve becomes somewhat steeper, i.e. in the N V T ensemble the system is stiffer than in the N P T ensemble, and this change is also seen for the three other potentials, indicating that release of the constraint in the orthogonal directions decreases the tangent modulus, a result that certainly makes sense.

Finally we mention that, for the BKS as well as SHIK potential for silica, the fracture strain decreases significantly when one switches to the N V T ensemble although the brittle nature of the fracture is independent of the ensemble. (Note that, as for NS3, the system is stiffer when working in the N V T ensemble.) From these results we thus conclude that the choice of the ensemble is important, and thus the N V T ensemble should be avoided for this type of measurements.

In the linear elastic regime one can characterize the mechanical behavior of the sample by the elastic constants given by the Young's modulus from Eq. ( 8 e Ref. [START_REF] Lower | Inert failure strain studies of sodium silicate glass fibers[END_REF] (strength at 77 K, and the failure stress is overestimated, see discussion section in Ref. [START_REF] Lower | Inert failure strain studies of sodium silicate glass fibers[END_REF]). f Ref. [START_REF] Kurkjian | Intrinsic strength and the structure of glass[END_REF] (strength at 77 K)

C 11 = lim →0 E N V T t ( ) . (9)
From E and C 11 one can then obtain the bulk modulus B, the shear modulus G, and the Poisson's ratio ν using the following relations [START_REF] Bansal | Handbook of glass properties[END_REF][START_REF] Zhao | In-situ high temperature Raman and Brillouin light scattering studies of sodium silicate glasses[END_REF]:

E = 9B(C 11 -B) (3B + C 11 ) (10)
G = 3(C 11 -B) 4 (11) ν = E 2G -1 . (12)
The resulting values for the different potentials as well as the failure stress and failure strain are reported in Tab. 2. In this table we have also included the experimental values and one recognizes that overall the SHIK potential with the Wolf truncation is likely the most reliable potential. Note that this potential has the merit to be applicable not only to NS3 but, with the same set of parameters, also to SiO 2 and therefore also for compositions between the two systems [START_REF] Sundararaman | New interaction potentials for alkali and alkaline-earth aluminosilicate glasses[END_REF].

Summary and Discussion

The goal of this work has been to investigate how the properties of the glasses, in particular their mechanical behavior, depend on the details of the interaction potential used to simulate the system, and to also to study the influence of the ensemble in which the simulations are done. In agreement with previous studies we have found that the structural quantities show a relatively mild dependence on the interaction potentials, with a notable exception for the GS potential that predicts a phase separation into sodium rich/poor regions. The temperature dependence of the density as well as the distribution of the different silicon species show instead a noticeable dependence on the potential and hence these quantities are more sensitive indicators for evaluating the quality of a potential.

A very strong dependence on the potential is found for the fracture behavior of the sample:

While the Teter and HO potentials predict a plastic deformation at constant pressure even if the strain is 100%, the samples generated by the GS and SHIK potentials fail at a much smaller strain. Therefore the two former potentials should not be used for this type of studies. In addition, we have found that both structural and mechanical properties are not realistic when using the GS potential. If the sample is kept at constant geometry in the directions orthogonal to the pulling directions all potentials predict a brittle fracture.

Our results demonstrate that certain structural quantities show a dependence on the potential, but that the properties of the fracture show a much stronger one. The latter dependence is, however, definitely weaker than the ones found for dynamical quantities (see Ref. [START_REF] Hemmatti | Comparison of pair potential models for the simulation of liquid SiO 2 : Thermodynamic, angular-distribution, and diffusional properties[END_REF], for example). This can be understood from the fact that static properties depend mainly on the position of the particles close to the local minimum of the potential energy, while fracture and dynamic properties depend also on the position and height of the saddle points connecting these minima. Since for many systems there is only very limited experimental information about dynamic quantities, while elastic and inelastic measurements are much more abundant in the literature, the mechanical properties are thus a very attractive way to probe whether or not a potential is reliable.

A further important result of our study is the dependence of the fracture behavior on the system size. As in previous studies we have found that the structural quantities do depend only in a mild manner on the size of the system. In contrast to this, our study suggests that a reliable description of the fracture behavior can only be obtained for systems that have more than 75,000 particles, i.e. a box size of the order 10 nm for the case of silica and L = 16 nm (300,000 atoms) for NS3. For smaller sizes, the elastic regime shows almost no dependence on the system size, while this dependence becomes very strong after the failure point and the system will present a fracture behavior that is much more ductile than a sufficiently large system.

Finally we have also documented the fact that the fracture behavior of the glass samples can show a strong dependence on the nature of the boundary conditions. If the traction is done under constant cross section condition, the samples will break in a brittle manner. However, if one fixes the pressure, as it is the case as in an experimental set-up, some potentials predict a much more ductile fracture or even no fracture at all up to a strain of 100%. This shows that for this type of studies one should use the N P T ensemble and a potential that gives a realistic fracture behavior.

Figure 1 :

 1 Figure 1: Distance dependence of the various interaction potentials [panels a), c), and e)] and their corresponding forces [panels b), d),and f)] for the SiO, SiSi and OO pairs. Also included in the latter panels are the corresponding pair correlation function as predicted by the SHIK potential at 300 K (symbols).

Figure 2 :

 2 Figure2: Mean squared displacement for silicon in the silica (a) and the NS3 (b) samples during the N P T equilibration at high temperatures and at zero pressure. The different curves correspond to the potentials used and the dashed line shows the diffusive behavior. For SiO 2 the temperature is 3600 K. For NS3 the temperature is 3000 K except for the GS curve which is for 2100 K.

Figure 3 :Figure 4 :

 34 Figure 3: Temperature dependence of the density during the quench at zero pressure for silica and NS3, panels (a) and (b), respectively.

Figure 5 :

 5 Figure 5: SiOSi angle distribution for silica (a) and NS3 (b) glasses at 300 K. The vertical green lines indicate the range of most probable SiOSi angle extracted from experimental studies (see Ref.[START_REF] Malfait | 29 Si NMR spectroscopy of silica glass: T 1 relaxation and constraints on the SiOSi bond angle distribution[END_REF]).

Figure 6 :

 6 Figure 6: Distributions of oxygen speciations, panel (a), and Q n species, panel (b), for the NS3 glass at 300 K. The experimental data in panels (a) and (b) is from Ref. [40], and Ref. [41], respectively.

Figure 8 :

 8 Figure 8: Stress-strain curves as predicted by the SHIK potential for glasses at 300 K, produced with different cooling rates. Panels (a) and (b) are for silica and NS3, respectively.

Figure 9 :

 9 Figure 9: Mechanical responses of the glasses under uniaxial tension at 300 K. Panels (a) and (b)show the stress-strain curve and the tangent modulus for silica, respectively. The green symbols are the experimental data from Ref.[START_REF] Krause | Deviations from linearity in the dependence of elongation on force for fibers of simple glass formers and of glass optical lightguides[END_REF]. The green dashed line in panel (b) is the polynomial fit from Ref.[START_REF] Gupta | Intrinsic failure and non-linear elastic behavior of glasses[END_REF] to experimental data and its extrapolation to large strains. Panels (c) and (d): Same quantities as in panels (a) and (b) for the case of NS3.

Figure 10 :

 10 Figure 10: Influence of simulation ensemble on fracture at 300 K. The orange and magenta curves are the same as in Fig. 9c), i.e. N P T ensemble. The violet and blue dashed curves are the data obtained in the N V T ensemble.

Table 1 :

 1 Parameters for the various potentials. The oxygen charges for the HO and SHIK potentials, marked by a star, are given for the NS3 composition (see text for details).

	Parameters	BKS	GS	Teter	HO	SHIK
	A SiO [eV]	1388.7730	9022.8533	1844.7458	1742.1231	1120.5290
	A OO	18003.7572 50306.4259	13702.9050 10631.4994	23107.8476
	A SiSi	0.0000	0.0000	0.0000	865032008.0130 2797.9792
	A NaO		120304.5810 4383.7555	1854.3947	1127566.0
	A SiNa		0.0000	0.0000	81407.6194	495653.0
	A NaNa		0.0000	0.0000	2558.6814	1476.9000
	B SiO [ Å]	0.3623	0.2650	0.3436	0.3513	0.3457
	B OO	0.2052	0.1610	0.1938	0.2085	0.1962
	B SiSi	1.0000	1.0000	1.0000	0.0657	0.2269
	B NaO		0.1700	0.2438	0.2603	0.1450
	B SiNa		1.0000	1.0000	0.1175	0.1847
	B NaNa		1.0000	1.0000	0.1692	0.2935
	C SiO [eV• Å6 ] 175.0000	85.0927	192.5800	212.9333	26.1321
	C OO	133.5381	46.2981	54.6810	69.9590	139.6948
	C SiSi	0.0000	0.0000	0.0000	23.1044	0.0000
	C NaO		0.0000	30.7000	0.0000	40.5620
	C SiNa		0.0000	0.0000	0.0000	0.0000
	C NaNa		0.0000	0.0000	0.0000	0.0000
	q Si [e]	2.4	1.89	2.4	2.4	1.7755
	q Na [e]		0.4725	0.6	0.88	0.6018
	q O [e]	-1.2	-0.945	-1.2	-1.28	-0.9328
	R cut [ Å]	5.5/12.0	11.0/12.0	8.0/12.0	8.0/12.0	8.0/10.0

These

graphs demonstrate that the different potentials and forces depend strongly on the chosen set of parameters which gives thus a first indication that the predicted glass properties will depend on

  • , 147.7 • , 147.1 • and 147.7 • , and 152.3 • for the SHIK,

Table 2 :

 2 , 56 d 69.5 a 37.2 a 24.3 a , 22 d 0.232 a , 0.25 d 20.85 e 11.71 e , 7.5 f 2.431 a Elastic constants (Young's modulus E, longitudinal modulus C 11 , bulk modulus B), Poisson's ratio ν, failure strain f and stress σ f (GPa), and density ρ at 300 K, compared with experimental data.

	) and the longitudinal modulus

a Ref.

[START_REF] Bansal | Handbook of glass properties[END_REF]

b Ref.

[START_REF] Griffioen | Optical fiber mechanical reliability[END_REF]

c Ref.

[START_REF] Smith | DOE contract# DEAC040dpoo789[END_REF]

d Ref.

[START_REF] Januchta | Elasticity, hardness, and fracture toughness of sodium aluminoborosilicate glasses[END_REF]

Acknowledgements Z.Z. acknowledges financial support by China Scholarship Council (NO. 201606050112). This work was granted access to the HPC resources of CINES under the allocation A0050907572 attributed by GENCI (Grand Equipement National de Calcul Intensif).