Quantifying the internal stress in over-constrained glasses by molecular dynamics simulations

Topological constraint theory classifies network glasses into three categories, viz., flexible, isostatic, and stres- sed–rigid, where stressed–rigid glasses have more topological constraints than atomic degrees of freedom. Such over-constrained glasses are expected to exhibit some internal stress due to the competition among the re- dundant constraints. However, the nature and magnitude of this internal stress remain poorly characterized. Here, based on molecular dynamics simulations of a stressed–rigid sodium silicate glass, we present a new technique allowing us to directly compute the internal stress present within a glass network. We show that the internal stress comprises two main contributions: (i) a residual entropic stress that depends on the cooling rate and (ii) an intrinsic topological stress resulting from the over-constrained nature of the glass. Overall, these results provide a microscopic picture for the structural instability of over-constrained glasses.

connected to each other by some radial and angular constraints (the chemical bonds) [6]. In Phillips' approach, molecular networks comprise two kinds of constraints, viz., the radial bond-stretching (BS) constraints that maintain the inter-atomic distances fixed around their average values and the angular bond-bending (BB) constraints that prevent the inter-atomic angles from significantly deviating from their average values. As per Maxwell's stability criterion, glasses can be classified into three categories based on their connectivity: (i) flexible, if they comprise fewer independent constraints than atomic degrees of freedom (i.e., 3 per atom), (ii) stressed-rigid, if they exhibit more independent constraints than atomic degrees of freedom, and (iii) isostatic, if all the atomic degrees of freedom are exactly balanced by the inter-atomic constraints. Phillips predicted that isostatic glasses exhibit optimal glass-forming ability [6].
Based on these seminal works, Phillips and Gupta are the founding fathers of topological constraint theory (TCT) [2,7]. Since then, TCT has enabled the development of various analytical models relating the connectivity of glass networks to their engineering properties, including glass-forming ability [2,5,8], glass transition temperature [9,10], hardness [11][12][13], fracture toughness [14], and dissolution rate [15,16]. These models rely on the level of simplification offered by TCT, which reduces the complex, disordered atomic structure of glasses into simpler topological networks comprised of nodes [2,7].
Although they adopt different viewpoints, the approaches from Gupta and Phillips are mathematically equivalent and rely on the same fundamental idea, namely, that there exists an optimal degree of atomic connectivity for which the propensity for vitrification is maximum. This can be understood as follows. On the one hand, due to the deficit of constraints, flexible glasses feature some internal floppy modes, that is, some internal low-energy degrees of freedom-which may facilitate the structural reorganization that is needed for crystallization [17,18]. On the other hand, stressed-rigid glasses are over-constrained, so that some constraints are mutually redundant and may compete against each other-just like the angles of a triangle with three fixed edges cannot be changed to arbitrary values. This competition among constraints in over-constrained glasses has been expected to result in the formation of some internal stress, which arises from some internal "frustration" within the atomic network [5,6] (see Fig. 1). This stress is expected to play a key role to explain the propensity for stressed-rigid glasses to crystallize. For instance, (Na 2 O) x (SiO 2 ) 1-x glasses have been noted to be stressed-rigid for 0 < x < 20% [19,20]. This also corresponds to the range of compositions wherein sodium silicate glasses exhibit low glassforming ability [8,21]. This likely arises from the fact that the internal stress resulting from mutually-dependent constraints provides a thermodynamic driving force for crystallization. In turn, the formation of a crystalline structure wherein all the constraints are compatible with each other (due to the periodic structure of crystals) can alleviate internal stress, that is, by converting mutually incompatible independent constraints into compatible dependent constraints.
However, the nature and magnitude of this internal stress remain poorly characterized. Although it was noted that the internal stress may manifest itself in the Raman scattering patterns of glasses [22], it remains challenging to experimentally probe its magnitude. Alternatively, atomistic simulations offer a direct access to the atomic positions and inter-atomic forces and, hence, should allow one to compute the internal stress within a given atomic network. However, although the manifestation of internal stress has been studied in some "toy models" of glasses (e.g., triangular 2D lattices or spring networks [23][24][25][26][27]), no general method has been proposed thus far to isolate and quantify the internal stress in more realistic glasses.
Here, we present a new computational technique based on molecular dynamics (MD) simulations aiming to describe and quantify the internal stress in glasses and illustrate its application to an over-constrained sodium silicate glass. We show that, although the atomic network of the glass is macroscopically at zero pressure, it exhibits some internal stress, which manifests itself as some inter-atomic bonds being under compression, whereas others are under tension. Specifically, we show that the internal stress comprises two distinct contributions, namely, a residual stress of entropic origin arising from the fast cooling rate and an intrinsic topological stress resulting from the over-constrained nature of the atomic network. These results establish a microscopic foundation behind the structural instability of over-constrained glasses. By directly tracking the magnitude of stress within the network, this technique paves the way toward an accurate tracking of rigidity transitions (i.e., transition from unstressed-to-stressed networks) in structural glasses as a function of composition, temperature, or pressure [9,28,29].

Preparation of the glasses
To establish our conclusions, we perform here some MD simulations of a sodium silicate glass (Na 2 O) x (SiO 2 ) 100-x with x = 10%. This specific composition is here chosen as it belongs to the stressed-rigid domain of binary sodium silicate glasses (i.e., 0 < x < 20%) [19,20]. The initial liquid configuration consists of 3000 atoms, among which 200 are Na atoms, 900 are Si atoms, and 1900 are O atoms. These 3000 atoms are first randomly placed in a cubic simulation box while ensuring the absence of any unrealistic overlap. The system is then melted at 3000 K for 100 ps in the NPT ensemble at zero pressure to ensure the loss of the memory of the initial configuration. To assess the effect of the thermal history of the glass on its internal stress, we subsequently quench the liquid from 3000 to 0 K with five select cooling rates, namely, 1000, 100, 10, 1, and 0.1 K/ps in the NPT ensemble at zero pressure. All simulations are conducted with the LAMMPS package [30]. We adopt the well-established force-field parametrized by Teter [31], as it has shown to offer a realistic description of the structure and properties of sodium silicate glasses [32][33][34][35][36][37][38][39][40]. The short-range interaction cutoff is chosen as 8.0 Å and Coulomb interactions are evaluated by the Ewald summation method with a cutoff of 12 Å [32]. A timestep of 1 fs and the Noosé-Hoover thermostat/barostat is used for all simulations [41,42]. Fig. 2 shows a slice snapshot of the sodium silicate glass simulated herein.

Simulations of isolated atomic clusters
To isolate the internal stress acting in the atomic network of a glass, our approach consists in comparing the state of Si atoms present the network to those belonging to isolated atomic clusters (i.e., where no internal stress resulting from redundant constraints can be formed). To this end, we simulate a series of isolated atomic clusters. Each of the clusters considered herein consists of a distinct Q n unit, that is, an SiO 4 tetrahedral unit connected to n other Si atoms, i.e., comprising n bridging oxygen (BO) atoms and 4 -n non-bridging oxygen (NBO) atoms (see Fig. 3). The initial configuration of each cluster is first created by manually placing the atoms at their designated position, that is, based on the average values of the SieO and NaeO bond lengths, as well as those of the OeSieO, SieBOeSi, and SieNBOeNa inter-atomic angles [34]. The atoms are placed at the center of a large cubic box with fixed volume (with a side of 50 Å) and with no periodic boundary conditions. The clusters are then relaxed at 300 K for 100 ps and subsequently cooled from 300 K to 0 K in the NVT ensemble with a cooling rate of 1 K/ps while fixing to zero the linear and angular momentum of the cluster to avoid any drift of the system. The force-field and other simulation parameters are here similar to those used to simulate the bulk glasses (see Section 2.1).

Computation of the internal stress
To determine the magnitude of the internal stress acting within the network, we adopt here the "stress per atom" framework formulated by Thompson et al. [43]. Note that, strictly speaking, stress is ill-defined for individual atoms and can only be meaningfully calculated for an ensemble of atoms. Nevertheless, we define here the stress σ i of each atom i as being its contribution to the total virial of the system: where V i , m i , v i , and r i are the volume, mass, velocity, and position of the atom i, respectively, and F i is the resultant of the force applied on the atom i by all the other atoms in the system. Here, we define the volume V i of each atom based on its Voronoi volume. Several previous studies have utilized this approach to quantify the extent of local instabilities within the atomic network [44][45][46][47][48]. Note that, although the decomposition of the virial in terms of individual atomic contributions does not have a clear physical interpretation, it provides a convenient means to capture the local state of stress (i.e., compression or tension) experienced by each atom. Here, this local stress does not result in any macroscopic stress as the glass is maintained at zero pressure, that is, local compressive and tensile stress mutually compensate each other. By convention, a positive stress refers here to a state of tension, whereas a negative stress denotes a state of compression.
In the following, we calculate the internal stress within the atomic network as the difference between the stress experienced by the atoms in the bulk glass (referred hereafter as the "network stress") and that experienced in the isolated clusters (referred hereafter as the "reference stress").

Reference stress in isolated clusters
We first focus on the computation of the reference stress, that is, the local stress experienced by the atoms in the simulated isolated Q n clusters. Fig. 4(a) first shows the Voronoi volume of the central Si atom in each isolated cluster. We observe that the volume of Si atoms gradually increases from Q 1 to Q 4 unit. This is in agreement with the fact that Si-BO bonds are slightly more elongated than Si-NBO bonds, as BO atoms are equally attracted by the two neighboring Si atoms they connect to [38]. In contrast, Si atoms in Q 0 units exhibit the highest volume, which may arise from the fact that their Voronoi volume is not well defined due to the lack of surrounding atoms.
Based on these Voronoi volumes, the reference stress per Si atom can be calculated based on Note that, here, we focus on the Si atoms as the Voronoi volume of the O atoms is not well defined in the isolated clusters. Fig. 4(b) shows the reference stress per Si atom for each type of Q n unit. We observe that Si atoms are systematically under tension (which, by convention, corresponds here to a positive stress). The origin of this state of tension is depicted in Fig. 5 and can be understood as follows. Si and O atoms are mutually attracted to each other. However, the O atoms surround Si atoms mutually repulse each other due to mutual Coulombic repulsion (which results in a tetrahedral environment around the Si atoms). These two effects compete against each other since the attraction between Si and O brings O closer to each other, whereas their mutual repulsive brings them apart from each other. Overall, at equilibrium, the balance between these two forces result in the formation of an equilibrium structure wherein (i) the O atoms experience a state of compression as they are forced to partially overlap with each other whereas, (ii) in turn, the central Si atom undergoes a state of tension as the mutual repulsion among O atoms tend to stretch the central Si. This illustrates the fact that, even though the overall atomic structure is at equilibrium, some bonds may actually experience a local stress of tension or compression.
Further, we note that the reference stress per Si atom decreases with the number of surrounding BO atoms, reaches a minimum for Q 3 units, and slightly increases in Q 4 units. This supports the idea that Q 3 units represent an optimal state of connectivity for Si atoms. This is in agreement with the fact that sodium silicate glasses systematically exhibit an excess of Q 3 units with respect to what would be expected from a random model [34,35,49]. This also agrees with the fact that the fraction of Q 3 units in sodium silicate tends to increase upon relaxation, which suggests that this Q 3 structure energetically favored [33]. Finally, this echoes the fact that a silicate glass entirely made of Q 3 units would be topologically equivalent to a B 2 O 3 (or As 2 Se 3 ) glass, which is characterized by an isostatic state of rigidity [5]. It is worthwhile to note that the stress values presented in Fig. 4(b) should not be compared to typical values of glass strength. Indeed, the specific stress values are somewhat arbitrary (as they depend on which convention is used to define the volume of each atom). Also, inter-atomic forces are concentrated on very small surfaces and, hence, should indeed correspond to very large stress values.

Network stress in bulk glass networks
We now focus on the state of stress experienced by the atoms in the bulk network. Fig. 6(a) shows the Voronoi volume distribution of each atomic species in the sodium silicate glass (for a cooling rate of 1 K/ps). As expected, we find that Si atoms exhibit a volume that is fairly similar to those obtained in the isolated clusters (~7 Å 3 ). In turn, Na atoms show the largest average volume, which is in agreement with the fact that these atoms are typically found to be located in some empty pockets within the silicate network [37]. Finally, we observe that the volume of the NBO atoms is slightly larger than that of BO, which arises from the fact that the OeNa bond is larger than that of OeSi. Besides the average value of these distributions, their broadness is also of interest. Overall, the distribution of the volume of the Si atoms is the sharpest, which is in agreement with the fact that these atoms exhibit a well-defined tetrahedral environment. In contrast, BO and NBO atoms exhibit a significantly broader distribution, which echoes the fact that the inter-tetrahedra SieOeSi angle exhibits a greater variability than that of the intra-tetrahedron OeSieO angle [33]. Finally, we note that Na atoms exhibit the broadest distribution, which suggests that their local environment exhibits a great variability.
The distribution of the network stress of each atomic species is shown in Fig. 6(b). Interestingly, we find that Si and O atoms are systematically under tension and compression, respectively. This agrees with the interpretation presented in Fig. 5. As expected, the average network stress undergone by Si atoms exhibits an order of magnitude that is similar to those belonging to isolated clusters (although the network and reference stress differ from each other, see below). In turn,   Fig. 6. (a) Voronoi volume and (b) network stress distribution of each atomic species in a sodium silicate glasses quenched with a cooling rate of 1 K/ps. Positive and negative stress values denote a state of tension and compression, respectively. The data are fitted with Gaussian distributions. Note that the normalizations of select distributions are adjusted to improve the overall readability of the plot.
Na atoms exhibit an average network stress that is close to zero, which is in agreement with the fact that these atoms act as network modifiers, occupy the empty space left within the silicate network, and do not create any strong bond with the surrounding atoms [50].

Contributions to the internal stress
In the following, we focus on the characterization of the internal stress acting within the network, that is, the difference between the network stress in the bulk glass (see Fig. 6(b)) and reference stress in the isolated clusters (see Fig. 4(b)). We focus here on the Si atoms as the computation of their reference stress is more straight-forward than for the other atomic species. Nevertheless, note that, since the overall glass network is at zero pressure, it is equivalent to focus on the atoms under tension (i.e., Si) or on the atoms under compression (i.e., O atoms). To this end, we compute the network stress of each Si atom individually and determine which Q n structure it belongs to, that is, to infer its reference stress.
As shown in Fig. 7, we observe that the distribution of the Q n units hardly depends on the cooling rate. This contrasts with the trend observed in a (Na 2 O) 30 (SiO 2 ) 70 glass (simulated with the same method and forcefield), wherein the fraction of Q 3 units was found to increase (at the expense of Q 2 and Q 4 units) upon decreasing cooling rate [33]. This can be explained by the fact that the (Na 2 O) 30 (SiO 2 ) 70 glass belongs to the flexible domain. Therefore, its structure exhibits some degree of flexibility and can change upon relaxation toward lower fictive temperatures. In contrast, the sodium silicate considered herein belongs to the stressed-rigid domain and, hence, exhibits an atomic network that is fairly locked and unable to significantly reorganize upon decreasing cooling rate. Overall, the fact that Q n distribution only weakly depends on the cooling rate indicates that, in the present case, the average reference stress remains fairly constant with decreasing cooling rate (see Fig. 8(a)).
Finally, we compute the average internal stress based on the difference between the network and reference stress experienced by the Si atoms (see Fig. 8(a)). We find that, in contrast with the case of the reference stress, the average network stress per Si decreases upon decreasing cooling rate. As shown in Fig. 8(a), the dependence of the network stress σ on the cooling rate γ can be fairly well captured by a power law, in agreement with the prediction from mode-coupling theory [51]: where A and δ are some fitting parameters. Nevertheless, the extrapolation toward lower cooling rates values reveals that the network stress does not converge toward the reference stress and, rather, eventually plateaus σ = σ 0 at low cooling rate. This allows us to discriminate two contributions to the internal stress (see Fig. 8(b)), viz., (i) a residual entropic stress that depends on the cooling rate (σ res = (Aγ) 1/δ ) and (ii) an intrinsic topological stress resulting from the over-constrained nature of the glass σ int = σ 0 . The first contribution to the internal stress (residual entropic stress) results from the structural defects forming upon the fast quenching of the glass (e.g., small silicate rings [52]). This arises from the fact that, when the cooling rate is high, the glass "freezes" in some high-energy domains of the energy landscape, that is, the glass structure is virtually similar to that of a high-temperature supercooled liquid (i.e., high fictive temperature) [53]. In turn, as the cooling rate decreases, the system can reach a lower energy state before freezing, thereby achieving a lower fictive temperature and, hence, a lower residual entropic stress. This residual stress is here found to contribute up to about 3 GPa (at the highest cooling rate) to the total internal stress. However, when extrapolated toward infinitely small cooling rate, the internal stress converges toward a constant value, which allows us to isolate the intrinsic topological stress resulting from the over-constrained nature of the present glass (σ 0 ). In contrast to the residual stress, this topological stress arises from the fact that the glass simulated herein comprises more constraints than degrees of freedom. As a consequence, all the constraints cannot be satisfied at the same time-so that some constraints are experiencing a state of tension, while others are under compression (wherein tension and compression mutually compensate each other so that the overall glass is at zero pressure). This intrinsic topological stress is here found to be around 5 GPa per Si atom.   resulting average internal stress (difference between the network and reference stress) per Si atom as a function of the cooling rate. The data are fitted by a power law (Eq. 2) and extrapolated toward lower cooling rate values.

Conclusion
Overall, our MD-based stress computation method allows us for the first time to capture the magnitudes of the stress acting inside isolated atomic clusters (that is, acting within the SiO 4 polytopes) and, based on this knowledge, to isolate the internal stress within the glass network (that is, the stress acting in between the rigid polytopes of the network rather than inside the polytopes themselves). Although a fraction of the internal stress results from the fast cooling of the glass, the main contribution to the internal stress has a topological origin and arises from the redundant constraints found in over-constrained stressed-rigid networks. Overall, these results offer a microscopic interpretation to the notion of "internal stress" commonly used in topological constraint theory. By directly computing the magnitude of such internal stress, our new method should facilitate the accurate tracking of flexible-to-rigid transitions (i.e., transition from unstressed-to-stressed networks) in glasses as a function of composition, temperature, or pressure-without relying on any unproven guessed regarding the number of constraints created by each atom.

Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: