Strength of 2D glasses explored by machine-learning force fields

The strengths of crystals and glasses have garnered significant attention in statistical physics for their non-equilibrium nature^1–3 and are highly sensitive to the local arrangement of atoms.^4,5 In contrast to edge cleavage in crystals, failure of glasses often involves nucleation of voids in regions not necessarily at the crack tip. The amorphous network can trap the cracks, adding kinetic contributions to the energy cost of cracking, and the voids ahead of the crack tip could guide its advancement.^6,7 Recent studies show that nanostructured glasses in the form of nanocrystals^8–10 and paracrystals¹¹ can markedly enhance their fracture resistance. Introducing nanocrystalline domains¹² increases fracture toughness by $80%$⁠, from $0.6$ to 1.1 MPa $m 1 / 2$⁠, while paracrystalline glasses exhibit a threefold greater enhancement.¹¹ It is, thus, interesting to explore the failure process of nanostructured glasses by considering their atomic-level structures, beyond the continuum framework.¹³ 

Direct characterization of glass structures remains challenging.^14–18 Structural description of glasses is limited to the short-range order (SRO)¹⁰ and mediate-range order (MRO) parameters.¹¹ Notably, two-dimensional (2D) silica has recently been discussed in the literature since its successful synthesis. The crystalline and amorphous structures can be directly visualized by electron microscopes (EMs).^19–21 Exploring 2D silica, thus, allows direct characterization of their failure behaviors at the atomic level from both theoretical exploration and experimental observations, which remain largely unexplored.

Atomistic simulations are powerful tools to reveal atomic-scale material kinetics.^22,23 The accuracy and efficiency of their prediction are limited by the fidelity of interatomic interaction models from first-principles calculations at the electronic-structure level to empirical force fields. Empirical models, such as Sundarararaman–Huang–Ispas–Kob (SHIK),^24,25 Du,²⁶ Bertani–Menziani–Pedone (BMP),²⁷ were developed for oxide glasses including silica and successfully applied to problems such as material discovery and mechanistic analysis.^28,29 However, the capability of modeling fracture is under question due to the presence of strong lattice distortion and dangling bonds at the crack tip and cleaved edges,^15,30–32 factors that typically do not account for force field parametrization.

Recent advances in developing machine-learning force fields (MLFFs)³³ allow us to revisit the non-equilibrium processes of material failure¹³ by using chemically accurate force fields trained in the Deep Potential-Smooth Edition (DeepPot-SE),³⁴ Neural Equivariant Interatomic Potential (NequIP),³⁵ and Neuroevolution Potential (NEP) frameworks.³⁶ The MLFF approach has recently been applied to fracture problems by expanding the training dataset to include intermediate structures during fracture, which are usually not well modeled by empirical force fields.^13,37 In this work, we develop a neural network force field for fracture (NN- $F 3$⁠) tailored for 2D silica across various nanostructures, including crystals and nanocrystalline (NCG), paracrystalline (Para), and continuous random network (CRN) glasses. We investigate their mechanical responses at the atomic level to gain insights into the kinetics of material failure.

To obtain the quantitative relations between energies, forces on atoms, and the atomic-level structures, spin-polarized DFT calculations are performed using the Spanish Initiative for Electronic Simulations with Thousands of Atoms (SIESTA) package³⁸ with numerical atomic orbitals (NAOs) at the double- $ζ$-plus-polarization (DZP) level.³⁹ Perdew–Burke–Ernzerhof (PBE) parameterization of the generalized gradient approximation (GGA) is used for the exchange-correlation functional.⁴⁰ Troullier–Martins norm-conserving pseudopotentials are chosen for the ion–electron interactions.⁴¹ 

The basis sets and pseudopotentials used in our first-principles calculations are based on Ref. 39. These settings are validated by their consistency in structural parameters of 2D crystalline silica, such as lattice constants, in-plane Si–O–Si bond angles, and bilayer thicknesses, obtained from other codes and settings,^42–44 and experimental results⁴⁵ (see Table S1 in the supplementary material).

The cut-off energy for electron wave functions is 500 Ry. The $k$-space is sampled by a Gamma-centered, $4×3×1$ Monkhorst-Pack grid for the $144$-atom model. For structures with open edges, sampling at the same $k$-point density is used. These settings assure a convergence threshold of 1 meV/atom.

To develop the NN- $F 3$ model for 2D silica, we first extend our approach from crystalline to glassy silica. Building upon our previous work on 2D crystals,³⁷ we first pre-sample the 3D space of basal-plane strain states using a low-accuracy empirical force field (Tersoff⁴⁶). The strain sweeping is conducted at a rate of $1× 10 − 5 ps − 1$ under controlled conditions of 10 K using a Nosé–Hoover thermostat with a damping constant of 1 ps. Structures generated from this pre-sampling process undergo subsequent DFT calculations. We employ the end-to-end, symmetry-preserving DeepPot-SE framework³⁴ and DeePMD-kit⁴⁷ for training the model. Despite these efforts, the NN-F $3$ model trained on this initial dataset shows some inaccuracies. Notably, stress predictions near peak strain and Poisson ratios at high strain deviate from the reference DFT calculations,³⁷ highlighting the constraints inherent in sampling atomic-level structures generated by the Tersoff force field.

To enhance predictions of atomic forces in highly distorted structures at crack tips and cleaved, undercoordinated edges, we expand the initial dataset. Implementing an active learning strategy (Training-Exploration-Labeling^13,48), we iteratively explore the most relevant structures based on predefined criteria related to deviations in atomic forces. Initially, four NN- $F 3$ models are trained using different seeds for random number generation during parameter initialization.⁴⁸ MD simulations using these models are conducted in the Atomic Simulation Environment (ASE),⁴⁹ generating trajectories and computing atomic forces under identical loading conditions as the pre-sampling process. The Query by Committee⁵⁰ algorithm is employed to screen atomic-level structures, selecting those with maximum standard deviations (SD) of atomic forces exceeding 0.1 eV/Å among the four NN- $F 3$ models. Structures containing crack tips and open edges were identified by atoms with coordination numbers below 2 for O atoms. Additionally, structures with an SD of atomic forces greater than 0.2 eV/Å are selected. These screened structures undergo DFT calculations for labeling, enriching the dataset. Consequently, atomic-level structures containing crack tips are incorporated into the dataset following the MD exploration process.

For 2D glassy silica, we employ an active-learning strategy based on the NN- $F 3$ model previously trained for 2D crystalline silica. Initially, we generate all 144-atom inequivalent structures, each with non-isomorphic bonding networks, through bond rotation.⁵¹ This process results in 257 structures distributed across 34 distinct ring distributions. From these, we select $45$ structures for active learning across the full strain space. Structures exhibiting a maximum SD of atomic forces exceeding 0.2 eV/Å (from the four NN- $F 3$ models) are chosen for subsequent DFT calculations. We conduct 30 iterations of active learning to achieve convergence, resulting in a reduction of the maximum SD of atomic forces to below 0.2 eV/Å. The final dataset comprises 264 420 data frames, represented in a sketch map illustrating the diversity of atomic-level structures in the training dataset, encompassing lattices under varied strain states, cleaved edges, and structures containing crack tips [Fig. S1(a) in the supplementary material]. While these 144-atom structures are too small to distinguish NCG, Para, and CRN structures, they effectively capture the essential structural characteristics of both crystalline and glassy silica.

In training the DeepPot-SE framework, we set the sizes of the embedding and fitting networks to $(25,50,100)$ and $(240,240,240)$⁠, respectively, with a cut-off radius of 8.0 Å and a smoothing parameter of 1.0 Å. The batch size is $8$⁠, and the hyperparameters pref_e and pref_f, which determine the weights of energy and force losses in the total loss function, are set to $1.0$ and $10.0$⁠, respectively. We employ Adaptive Moment Estimation (Adam) optimization over $3× 10 6$ batch steps, starting with a learning rate of $0.001$ that exponentially decays to $1.0× 10 − 8$ by the end of training. 50% of our dataset is used for training, with the remaining 5% reserved for validation.

To simulate the fracture of 2D silica, we use the large-scale atomic/molecular massively parallel simulator (LAMMPS).⁵² The structures (NCG, Para, and CRN) are simulated with a cubic box of $200 3$ Å, which contains $∼$15 000 atoms, $∼$4 times larger than previous studies.⁴ The fracture tests were conducted using an athermal quasistatic (AQS) protocol. In AQS, the structure is deformed at a strain rate of $5× 10 − 5$ per step, followed by damped dynamics with a viscous rate of 1  $ps − 1$ to relax the structure under a force-on-atom threshold of $10 − 4$ eV/Å.

The performance of NN-F $3$ is validated for 2D silica crystals and glasses through the predicted energies, forces, and stress–strain relations. The root mean square error (RMSE) of the energy per atom, the interatomic forces, and the in-plane stress are below $0.79$ meV/atom, $37.6$ meV/Å, and 24 N/m, respectively, for the validation dataset [Figs. S2(a)–S2(c) in the supplementary material]. The RMSE of the phonon spectrum measured from the DFT results is below 0.2 meV [Fig. S2(d) in the supplementary material]. The uniaxial stress–strain relations of 2D silica crystals and glasses with structures not in the 45 structures of 2D silica exhibit excellent agreement with the reference DFT calculation results (Figs. S2 and S3 in the supplementary material). These results validate the predictions of equilibrium properties. To verify the accuracy of NN-F $3$ in describing the non-equilibrium fracture process, we study the fracture of crystalline silica using both NN-F $3$ and DFT-based Born–Oppenheimer MD simulations. Runs of 0.3 and 0.4 ps are conducted for the zigzag and armchair crack edges, respectively, revealing almost identical structural evolution in the NN-F $3$ and DFT results (Fig. S4 in the supplementary material).

A critical issue in the study of non-crystalline materials is whether the structural models accurately reflect those observed in experiments, particularly concerning their non-equilibrium properties. The complexities in the chemistry and atomic-level structures of glasses present a challenge in verifying whether structures generated by methods like bond swapping⁵³ or simulated annealing accurately capture the essential structural features. Theoretical models were proposed with a CRN of Si–O tetrahedra in silica and by asserting that alkali or alkaline earth metals are glass modifiers that disrupt the network.⁵⁴ However, a direct comparison with experimental evidence at the atomic level is difficult. Consequently, the agreement is often assessed by the physical properties of the generated structures such as the density, modulus, and radial distribution functions (RDFs).¹⁰ 

2D amorphous materials, such as glassy silica^20,21 and monolayer amorphous carbon (MAC),^55,56 were observed under an EM, which allow for direct determination of the atomic positions. Experimental evidence shows that compared to MAC featuring a sub-crystalline 2D structure and significant out-of-plane displacements,⁵⁷ amorphous silica exhibits a CRN structure and lacks out-of-plane fluctuation for its symmetric three-atom-layer structure.⁴ 

2D glasses in a CRN structure can be generated by consecutive Stone–Wales (SW) rotations, where the topological disorder is controlled and measured through the ring statistics.⁵¹ A Monte Carlo dual-switch procedure is deployed to the network of Si atoms,^4,58 the position of which is adjusted by minimizing a spring-like potential for ring-to-ring interaction in the dual space after each Monte Carlo move. We use a fictitious temperature of $10 − 4$ and $10 4$ Monte Carlo steps and a value of $α=1/3$ for the Aboav–Weaire law,⁵⁹ which aligns with experimental evidence and suggests that small rings tend to sit around larger ones.⁶⁰ Structures of NCG and Para are constructed by preserving crystallites with sizes of 2.8 and 1.4 nm during structural transformation, respectively. Computer-generated NCG, Para, and CRN structures [Figs. 1(a)–1(c)] exhibit identical ring distributions [Fig. 1(d)]. The glassy structures are rich in pentagons and heptagons, with rare quadrilaterals, which agrees well with experimental characterization using scanning tunneling microscopy (STM).¹⁹ Their RDFs [Figs. 1(e) and 1(f)] and angle distribution functions [ADFs, Figs. 1(h) and 1(i)] are also similar and align well with the experimental data. Minor discrepancies in ADF compared to experimental data can be attributed to two factors: the inaccuracy in determining atomic positions and the presence of structural defects within the materials. Pinpointing exact atomic locations, particularly for oxygen atoms, remains challenging, often resulting in uncertainties observed in STM images.¹⁹ Furthermore, defects are commonly overlooked in silica model construction, yet they significantly influence residual stress distribution and O–Si–O bond angles. The RDFs and ADFs of Para, CRN, and NCG structures show minor differences (Fig. S5 in the supplementary material). The RDF comparisons between CRN and Para structures align with previous findings.¹¹ Para structures notably feature a narrower angle distribution (⁠ $SD= 5.73 °$⁠), compared to CRN (⁠ $SD= 6.32 °$⁠) and NCG (⁠ $SD= 6.40 °$⁠) structures, despite similar mean angles (⁠ $108.90 °$ for Para, $108.91 °$ for NCG, and $108.99 °$ for CRN, respectively). These results suggest a higher degree of structural order in the Para structures, possibly indicative of medium-range order (MRO).

Virial coefficients (⁠ $V θ= d θ T V⋅ d θ$⁠) for equilibrium, undeformed samples are computed along tensile directions (⁠ $d θ$⁠), following the definition of normal stresses [Fig. 2(a)]. The virial tensor, which relates to the stress tensor as $σ=− V/a$ (⁠ $a$ is the average area of atoms), characterizes residual stress induced by material heterogeneity. The spatial distribution of atomic sites exhibiting large $| V θ|$ values reveals the presence of “force chains” within the amorphous bonding network of the equilibrium or undeformed state.

There is a strong correlation between nucleation events and the value of $V θ$ [Fig. 2(b)], enabling the prediction of 2D glass strength from their equilibrium structures. Voids typically nucleate at atomic sites with high residual stress. Figure 2(c) shows the cumulative probability density function (CDF) of $V θ$ at the sites of void nucleation, ranging from the smallest (the most stretched) to the largest (the most compressed) values. Unlike a uniform distribution, the shape of CDF indicates that a few highly stretched states of atomic stress in equilibrium determine the strength of 2D glasses under tension.

On the other hand, the Si–O–Si bond lengths (distances between the silicon atoms) measuring the local strain also exhibit a strong correlation with the strain to void nucleation [Fig. 2(d)]. The lengths of critical bonds that break in the damped dynamics follow a narrow distribution with a mean of 3.6 Å and an SD of 0.07 Å, respectively [Fig. 2(e)]. The reduction in the values of virial/bond length descriptors with the size of models indicates the statistical nature [Figs. 2(b)–2(d)]. Heptagonal and octagonal rings are more prone to breakage than hexagonal rings, leading to more frequent void nucleation around larger rings (Fig. S6 in the supplementary material). These arguments made to CRN models apply to NCG and Para as well.

The elastic responses of materials are lost after void nucleation. Catastrophic fracture may be triggered by further increasing the loading amplitude. The strength of glasses is much reduced compared to the crystals, from 26.2 to 10 N/m for CRN. The stiffness decreases from 154.5 to 79.9 N/m. For the glassy phases, although the SROs and MROs (e.g., RDFs, ADFs, ring statistics) are close, the strength of Para and NCG structures increases slightly by $5%$ and $4%$ from CRN, respectively [Fig. 3(a)]. Significantly, the locations of fracture nucleation differ from those of void nucleation [Figs. 3(c)–3(k)]. The CDF of the minimum distances between the sites of void nucleation and the crack suggests a weak correlation (Fig. S7 in the supplementary material). This finding contrasts with crystals, where nucleation typically leads directly to fracture, indicating a more pronounced disorder-trapping effect over lattice trapping, which operates on a larger length and time scale.

The edges cleaved by fracture display distinct characteristics across the three representative glassy structures. In NCG, cracks typically propagate through large crystalline grains ahead of the crack tip and cleave crystalline facets, as the energy cost of crack deflection is relatively high. In contrast, the smaller and more distorted crystalline grains in Para hinder straight crack propagation, resulting in a more tortuous path. Conversely, the highly disordered CRN structures facilitate a relatively smooth fracture path. However, the contrast in the strengths [Fig. 3(a)] and edge roughness [Fig. 3(b)] of NCG, Para, and CRN structures is not significant in our simulations, suggesting a weak nanostructuring effect. The reportedly strong toughening effects observed in experiments may, thus, be attributed to the chemical heterogeneity in the crystalline and glassy phases or an enhanced 3D structuring effect compared to 2D, which awaits further exploration.^11,12

There is a signature difference in the fracture kinetics and patterns in crystals and glasses. In crystals, an advancing crack leaves atomistically smooth crystalline facets behind, although lattice kinks joining these facets in chiral edges can lead to roughening at a length scale larger than the lattice constants. In contrast, in all three amorphous structures, the cracks become ensnared within the disordered phase, where the cleaved edges show roughness at the length scale of a single Si–O–Si bond [Fig. 4(b)].

Crack propagation is facilitated by void nucleation and coalescence in the amorphous phases. In our AQS simulations, the number of voids increases first by their nucleation and then declines by aggregation [Fig. 4(c)]. Propagating cracks are predominantly characterized by the relatively large voids within the samples [Figs. 4(d) and 4(e)]. Following the Griffith theory, we find that the critical crack size is $c cr= 4 γ π E ε 2=15.08$ Å, where the values of $γ$ and $E$ are the edge energy densities and Young’s modulus. The value of $E$ is extracted from the stress–strain relation of the CRN glass [Fig. 3(a)], and $γ$ is the edge energy density calculated for 2D crystalline silica that is the lower bound of the value for glasses. This critical value well separates the propagating and non-propagating voids, validating the linear elastic fracture mechanics theory at the Ångstrom scale. It should be noted that using $γ$ instead of fracture toughness in the estimation excludes the non-equilibrium process during fracture,¹³ thus undervaluing $c cr$⁠. The stress field ahead of the crack tip, averaged over several samples, follows the $r − 1 / 2$ singularity. However, significant structural heterogeneity in each sample dominates the $K$ field and may ultimately determine the fracture strength⁶⁷ [Fig. 4(f)]. Voids can, thus, nucleate at the weak spots that are not necessarily at the crack tip, and the advancing of cracks can be guided by voids nucleated in front of the crack tip and along its path.

Dynamic fracture simulations are conducted to explore the additional inertial effects. A large sample of $100×50 nm 2$ is stretched at strain rates ranging from $10 − 5$ to $10 − 3 ps − 1$ at 0.1 K. The Nosé–Hoover thermostat with a damping constant of 1 ps is used. At a low strain rate of $10 − 5 ps − 1$⁠, fracture proceeds by edge cleavage and lattice trapping and edge cleavage, which are guided by voids nucleated ahead of the crack tip, defining the crack paths [Fig. 4(g), Movie S1 in the supplementary material]. Fragmentation intensifies by void nucleation across the sample at higher strain rates (Fig. S8, Movie S2 in the supplementary material). These phenomena have been supported by numerous simulations^68–70 and experimental studies^71–73 of glassy structures. The kinetic and dynamical effects elucidated here through atomistic simulations provide further insights into the mechanisms of failure.

We discussed the strength and fracture behaviors of 2D silica by developing and using a neural network force field for fracture (NN-F $3$⁠) to simulate the failure process, which features chemical accuracy and the capability to model problems of relatively large size, such as fracture. We find that the strength of 2D glasses is determined by the heterogeneous residual stress and can be reasonably measured by the local virials or bond lengths in the equilibrium state. Beyond this strength, voids are nucleated in the amorphous network but do not always trigger fracture propagation. Only those beyond a critical size defined by the Griffith theory can advance and proceed with continuous coalescence and bridging events. The kinetics of fracture propagation is renormalized by a disorder-trapping effect in addition to lattice trapping and is often guided by voids nucleated ahead of the crack tip. These findings enrich our understanding of material failure at the atomic level, where the effects of disorder and structural heterogeneity are highlighted with the chemical complexity precluded by choosing 2D silica in the study. Further experimental studies, potentially conducted in situ, are expected to add more insights into the problem.^15,16

The supplementary material includes the sketch map of the dataset used to train the neural network force field for fracture (NN-F 3), technical details and validation of NN-F 3, and additional simulation results supporting the discussion.

This study was supported by the National Natural Science Foundation of China (NNSFC) through Grant Nos. 12425201 and 52090032. The computation was performed on the Explorer 1000 cluster system of the Tsinghua National Laboratory for Information Science and Technology.

The authors have no conflicts to disclose.

Pengjie Shi: Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (lead); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Zhiping Xu: Conceptualization (lead); Data curation (supporting); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (supporting); Project administration (lead); Resources (lead); Software (supporting); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Figshare at http://doi.org/10.6084/m9.figshare.25688664, Ref. 74.