Structural signatures for thermodynamic stability in vitreous silica: Insight from machine learning and molecular dynamics simulations

The structure-thermodynamic stability relationship in vitreous silica is investigated using machine learning and a library of 24 157 inherent structures generated from melt-quenching and replica exchange molecular dynamics simulations. We find that the thermodynamic stability, i.e., enthalpy of the inherent structure $\left(e_{\mathrm{IS}}\right)$, can be accurately predicted by both linear and nonlinear machine learning models from numeric structural descriptors commonly used to characterize disordered structures. We find that short-range features become less indicative of thermodynamic stability below the fragile-to-strong transition. On the other hand, medium-range features, especially those between 2.8 and $\sim 6$ Å, show consistent correlations with $e_{\mathrm{IS}}$ across the liquid and glass regions, and are found to be the most critical to the stability prediction in this study among features from different length scales. Based on the machine learning models, a set of five structural features that are the most predictive of the silica glass stability is identified.

Figure 1

$e_{\mathrm{IS}}$ distributions of the inherent structures from REMD and melt-quenching MD simulations, respectively normalized by the total number of structures from each technique. The number of structures from REMD are more than 10 times of that from the melt-quenching MD. The black peaks from the melt-quenching simulations correspond to the inherent structures collected below $T_{\mathrm{f}}$ that have similar inherent enthalpies. The inset shows the distribution in a broader range with the $e_{\mathrm{IS}}$ of $\alpha$-quartz noted by the dashed line.

Figure 2

Plots of $e_{\mathrm{IS}}$ vs the numbers of coordination defects: (a) ${\mathrm{O}}^{\mathrm{III}}$, (c) ${\mathrm{O}}^{\mathrm{I}}$, (d) ${\mathrm{Si}}^{\mathrm{VI}}$. (b) shows the slope of $e_{\mathrm{IS}}$ with respect to the number of ${\mathrm{O}}^{\mathrm{III}}$ from (a), calculated using moving average. A transition can be observed at $e_{\mathrm{IS}}$ corresponding to a fictive temperature of around 3100 K. The dashed lines and the shaded region in (b) are added to guide the eye.

Figure 3

Distribution of (a) the Si-O-Si bond angles and (b) the torsion angles in silica structures generated by REMD, melt-quenching MD with cooling rates of 0.1 and 100 K/ps. The $e_{\mathrm{IS}}$ of the 100 K/ps and 0.1 K/ps melt-quenched structures are 0.052 and 0.010 eV/atom above the REMD structure, respectively. The plots have been smoothed for clarity. The dashed line in (a) shows the data from x-ray diffraction experiments [33, 34].

Figure 4

Plots of $e_{\mathrm{IS}}$ vs (a) the mean, (b) the variance of the Si-O-Si bond angles, and (c) the full width at half maximum (FWHM) of the main peak in the Si-O-Si bond angle distribution.

Figure 5

(a) Simulated total pair distribution functions (PDF) of the three silica structures generated by REMD, melt-quenching MD with cooling rates of 0.1 and 100 K/ps. (b)–(d) show the partial PDFs. The vertical dashed lines show the peaks positions identified by experiments [39, 40].

Figure 6

Plots of $e_{\mathrm{IS}}$ vs (a) the FWHM of the first peak in the total PDF, (c) the height of the first peak in the total PDF, and (e) the height of the second peak in the total PDF. We smooth the total PDFs using a window length of 0.1 Å before extracting the peak information. (b) and (d) show the derivatives of $e_{\mathrm{IS}}$ with the width and the height of the first peak in PDF, i.e., the slopes in (a) and (c), respectively. A transition can be observed here in (b) and (d), similar to Fig. 2.

Figure 7

(a) Simulated structure factors of the three silica structures generated by REMD, melt-quenching MD with cooling rates of 0.1 and 100 K/ps. (b)–(d) Plots of $e_{\mathrm{IS}}$ vs (b) the height, (c) the FWHM, and (d) the location of the FSDP. The dashed line in (c), which is a guide for the eye, shows a transition in the trend.

Figure 8

Ring size distributions (i.e., the numbers of 3- to 10-member rings) in silica glasses generated by REMD, and melt-quenching MD with cooling rates of 0.1 and 100 K/ps. The extent of the variations in the ring size distribution for glass structures with a given stability can be seen in Fig. 9.

Figure 9

Plots of $e_{\mathrm{IS}}$ vs the numbers of primitive rings: 3,4,5,6,8,10-member rings in (a)–(f), respectively. The plots for 7- and 9-member rings are shown in SM, which are similar to those for 6,8,10-member rings [32].

Figure 10

Plots of $e_{\mathrm{IS}}$ vs density. The structures from REMD and melt-quenching MD are distinguished by color.

Figure 11

Performances of (a) the lasso model, and (b) the artificial neural network model. The $e_{\mathrm{IS}}$ predicted by the models agree well with the $e_{\mathrm{IS}}$ directly calculated using the force field on both the training and testing data.

Figure 12

Features with the largest absolute values of the weights in (a) the lasso model and (b) the ridge model.

Figure 13

Mean-squared errors of the ANN model trained (a) without a feature group and (b) only with a feature group. The number of features in each group is noted in parentheses. Feature Set 1 and Feature Set 2 are the two most important groups of features found by the lasso method and the decision tree model, respectively. The dashed lines show the original MSE with the full set of inputs.

Figure 14

Mean squared errors of the ANN model trained without or only with features grouped based on their length scale. The number of features in each group is noted in parentheses.