Modelling atomic and nanoscale structure in the silicon–oxygen system through active machine learning

Domain: earth_science
Published: 2024-03-02
Editorial outcome: minor_revision

Abstract
Silicon–oxygen compounds are among the most important ones in the natural sciences, occurring as building blocks in minerals and being used in semiconductors and catalysis. Beyond the well-known silicon dioxide, there are phases with different stoichiometric composition and nanostructured composites. One of the key challenges in understanding the Si–O system is therefore to accurately account for its nanoscale heterogeneity beyond the length scale of individual atoms. Here we show that a unified computational description of the full Si–O system is indeed possible, based on atomistic machine learning coupled to an active-learning workflow. We showcase applications to very-high-pressure silica, to surfaces and aerogels, and to the structure of amorphous silicon monoxide. In a wider context, our work illustrates how structural complexity in functional materials beyond the atomic and few-nanometre length scales can be captured with active machine learning.

Introduction
Elemental silicon and its oxide, silica (SiO2), are widely studied building blocks of the world around us1: from minerals in geology to silicon-based computing architectures; thin-film solar cells in which amorphous silicon is the active material2; or zeolite catalysts based on the SiO2 parent composition3. Some of these materials have a single phase and are precisely defined on the atomic scale, whereas others show longer-ranging, hierarchical structures and varying degrees of disorder. For example, silica aerogels contain pores with sizes of 5–100 nm, leading to very low thermal conductivity and making aerogels promising candidates for thermal insulation4. Under pressure, SiO2 shows amorphous–amorphous transitions to structures exceeding sixfold coordination5, crystallisation from the amorphous phase under shock compression6, and conversely the formation of complex disordered phases from crystalline SiO2 7. Beyond fundamental studies, there is much technological importance in silicon–oxygen phases with nanoscale structure—the interface between Si and SiO2 is essential in silicon metal-oxide semiconductors, and defects at this interface have been investigated for decades8–10.

A material in the binary silicon–oxygen system which is in fact dominated by such interfaces is the so-called silicon monoxide (SiO). The structure of SiO was controversially discussed for long11,12; today, it is known as an nanoscopic mixture of amorphous Si and SiO2 13–15. Initial applications of SiO have been in protective layers for mirrors16 or dielectrics for thin-film capacitors17; more recently, the same material has emerged as a promising anode material for lithium-ion batteries18,19. However, to be able to fully exploit SiO in next-generation energy-storage solutions, it would be valuable to understand the features of the nanoscopic structure on an atomistic level.

To develop atomic-scale models of complex materials such as SiO, molecular-dynamics (MD) computer simulations have become a central research tool. While there are now plenty of interatomic potentials for silicon20–22 and silica23–25, the number of potentials for the mixed (i.e., full binary) system is limited due to its chemical complexity26–29. Alongside established, empirically fitted potentials based on physical models, alternatives based on large datasets and machine learning (ML) have emerged in recent years. These models have been fitted for silicon30 as well as silica31 and also for the more complex silica–water system32. ML potentials promise the accuracy of first-principles methods such as density-functional theory (DFT) for a small fraction of the cost. ML potentials are now firmly established in the field of computational materials science and their application to homogeneous phases has been well documented.

In the present work, we describe a unified computational model for the Si–O system that we have obtained with the help of an active-learning scheme for local environments. We extract representative atomic environments from large-scale simulations and embed them in a melt-quenched amorphous matrix, allowing us to sample representative environments for the fitting of accurate ML potentials. Our final model shows high accuracy across a wide configurational space including high-pressure silica, silica surfaces, and mixtures of silica and silicon. We showcase the usefulness of the method by creating fully atomistically resolved, 10-nm-scale structure models of amorphous and partially crystalline SiO.

Results
Active learning for nanoscale structure
We have developed a comprehensive dataset of atomistic structures and quantum-mechanical reference data for the binary Si–O system, as well as an interatomic potential fitted to that database in the atomic cluster expansion (ACE) framework33–35. We initialised the protocol with two existing datasets for silicon (Bartók et al. 30) and silica (Erhard et al. 31) respectively, and we then gradually explored the relevant configurational space using the active-learning workflow illustrated in Fig. 1. Quantum-mechanical reference ("training”) data for energies and forces were obtained with the strongly constrained and appropriately normed (SCAN)36 exchange–correlation functional for DFT, which shows good performance for elemental silicon37 and the various silica polymorphs31.

Our active-learning workflow follows three main tracks: high-pressure bulk silica, silica surfaces, and non-stoichiometric SiOx systems (Fig. 1a). The individual tracks are kept separate during initial training, i.e., they do not share their newly generated training data; however, in the end, all structures are merged into one comprehensive database.

The single subtracks are further divided into stages. In the first stage, we added initial structures, e.g., for crystalline high-pressure polymorphs or surface models. In the next stage, we fitted moment tensor potential (MTP) models38 to the database and used these MTPs to explore configurational space in MD and to identify new structures by active learning39. Energies and forces for new structures were computed with DFT and added to the database. This process was iterated until the extrapolation threshold (Supplementary Note 1B) was not exceeded during the MD trajectories anymore.

The third stage, highlighted in red in Fig. 1a, is the most important part of our workflow, and is based on large-scale simulations in each track. We used 2–4 MTPs trained on the same database to estimate a per-atom committee error, as is commonly done for neural-network potentials40. For atoms with high uncertainty (Supplementary Note 1C), we extracted the environments into smaller, “DFT-sized” cells by an approach that we call amorphous matrix embedding (Fig. 1b). After identifying an atom with high uncertainty, we cut out a cube containing the corresponding environment of the atom. This cube has a size which is feasible for performing DFT computations; it is generally chosen larger than twice the cut-off of the potential. After extracting the cube, the atoms within the cut-off radius of the atom with high uncertainty are kept fixed. The remaining structure is melted in an ML-MD simulation to create an amorphous matrix and smooth boundaries. Details of the procedure can be found in Supplementary Note 1C. We note very recent, related approaches to isolating fragments for active learning based on minimising the uncertainty for boundary atoms41,42.

The final database was obtained by merging the data of all tracks together, including some additional samples such as clusters and vacancies. This database contains 11,428 structures with a total of ≈1.3 million atoms (Supplementary Table 1). For validation, we held out 5% of these structures from training, selected at random.

Performance
The final potential is a complex non-linear ACE model, obtained by summation of one linear and seven non-linear ACE terms (Methods). This approach allows a more flexible description than just a linear or Finnis–Sinclair-like embedding, at only moderately higher computational expense. The resulting potential has a test-set root mean square error (RMSE) of 16.7 meV atom−1 for energies and 306 meV Å−1 for forces. These errors are averaged over the full dataset, however, and so they are not in themselves sufficient to characterise the quality of the potential. For example, they refer to a highly heterogeneous set of structures, with target energy values spanning more than 8 eV atom−1, and a range of forces of 40 eV Å−1 covered by the database. Furthermore, the numerical accuracy of the potential in certain parts of configurational space (e.g., crystalline polymorphs) is far more important than in others (e.g., liquid and amorphous structures).

In Table 1, we therefore show the performance of our model on different separate test sets. The complex non-linear ACE is compared to our previous silica GAP model described in ref. 31 ("SiO2-GAP-22” in the following), and also to simpler ACE models fitted to the new database using linear and Finnis–Sinclair-like embeddings, respectively. Indeed, the complex non-linear ACE potential is the only one among the three which achieves comparable errors to SiO2-GAP-22 for amorphous and crystalline silica structures. In contrast, for amorphous elemental silicon, mixed-stoichiometry as well as high-pressure phases the complex non-linear ACE is significantly more accurate than SiO2-GAP-22, since these structures are not part of the GAP database. This table therefore indicates the main challenge–and its solution–in the present model compared to the previous GAP: both are highly accurate for crystalline (≈1 meV atom−1) and bulk amorphous (≈5 meV atom−1) SiO2, but our ACE model caters to a much wider range of scenarios outside of the 1:2 stoichiometric composition.

Figure 2 shows the phase diagram of SiO2 calculated by thermodynamic integration43,44 using the ACE potential compared to a CALPHAD phase diagram from the literature45. The ACE and CALPHAD predictions agree well throughout, and for the boundary between quartz and coesite we observe almost quantitative agreement. In contrast, the cristobalite and tridymite phases seem to be over-stabilised. At 0 GPa, the melting point is notably overestimated (about 2400 K, compared to ≈2000 K experimentally45); moreover, the phase stability regions of both phases are more extended than in the reference. To illustrate the sensitivity of the analysis to small errors in predicted energies, we added a fictitious energy penalty of 5 meV atom−1 for cristobalite and tridymite (Supplementary Fig. 1a); in this case, the transition lines agree much better with the CALPHAD reference than before. Further numerical tests showed that the tridymite–cristobalite transition line, in particular, is strongly affected by small shifts in energy (Supplementary Fig. 1b–f). We thus conclude that the quantitative deviation seen in Fig. 2 is due to the inaccuracy of the underlying exchange–correlation functional, rather than indicating a shortcoming of the ACE approach. This is an example of the more general problem that any issues with the ground-truth method (e.g., numerical instabilities) will translate into the ML model.

High-pressure structural transitions of SiO2
Figure 3 characterises high-pressure properties of silica. In Fig. 3a, we show energy–volume curves of α-quartz, coesite, stishovite, α-PbO2-type, and pyrite-type silica as predicted by our ACE model and compared with DFT data, with which they agree well. In addition, we tested the behaviour of the model for rosiaite-type silica, which was recently observed in experiment46 and predicted theoretically47 for direct compression of α-quartz. In contrast to the structures mentioned before, this particular polymorph is not part of the training database. Nevertheless, the ACE model reproduces DFT data for this structure similarly well as for the other polymorphs.

Figure 3b shows an enthalpy–pressure diagram at 0 K. For lower pressures, there is a transition from α-quartz to coesite between 2.5 and 3.0 GPa, consistent with the predicted phase diagram (Fig. 2), followed by a transition to stishovite at 5.5–6.0 GPa. At higher pressures of ≈110 GPa, we observe the transition from stishovite to α-PbO2-type silica. Experimentally, rather than stishovite (rutile type), the structurally closely related CaCl2 (distorted rutile) type polymorph of silica is stable. The transition from CaCl2- to α-PbO2-type silica was observed at 120 GPa and 2400 K48. Given that our enthalpy data correspond to a temperature of 0 K, both values agree well with each other. For the transition of α-PbO2- to pyrite-type silica, our ACE model predicts a pressure of ≈246 GPa, in good agreement with the experimentally determined transition pressure of ≈260 GPa at 1800 K49. Finally, rosiaite-type silica46 is correctly identified as metastable over the pressure ranges studied.

Figure 3c shows the pressure evolution of the average coordination number (CN) of silicon atoms in amorphous silica, extracted from an MD simulation at room temperature and under isostatic pressure. The ACE results agree well with experiment up to about 50 GPa5,50, and with ab initio MD51 results over the whole pressure range. The good agreement with experiment is particularly pronounced for the data from ref. 50. Above 50 GPa, our model underestimates the average CN: at 175 GPa the experimental estimate is about 7; the ACE simulation predicts it to be 6. Importantly, this does not mean that there are no 7-fold coordinated environments, but there remain some 5-fold coordinated atoms as well, lowering the average (Fig. 3d). A possible reason for the good agreement with the ab initio result, but the deviation from experiment, might be the limited timescales in our simulations, which hinder a complete transition into higher-coordinated environments. Moreover, we note that computed X-ray-Raman spectra of the ab initio structures from ref. 51 are in good agreement with experiment indicating a lower CN. Other MD simulations also showed slightly lower CNs than the experimental values52. Figure 3e shows three different 7-fold coordinated environments extracted from the simulations.

A CN of 7 in amorphous silica might be surprising, since silicon is sixfold-coordinated in all crystalline silica polymorphs that are stable in this pressure range. However, the pyrite-type phase, which becomes thermodynamically stable at ≈240–260 GPa, contains silicon atoms with a 6+2-fold environment. A recent study found certain, but limited, similarities between these 7-fold environments in glassy and pyrite-type silica52.

In Supplementary Fig. 2, we show two additional structural fingerprints which have been commonly analysed in experiment: the position of the first sharp diffraction peak and the Si–O bond length. For both cases, our simulations show good agreement with experiment.

SiO2 surfaces and aerogels
Figure 4 tests the ability of the potential to accurately predict surface energies. We begin with validation for different α-quartz surfaces: we created surface slab models, relaxed them with the ACE potential, and evaluated the energetics, and therefore the surface energy per area, with DFT single-point computations (Fig. 4a). The ACE results agree well with DFT, especially considering that the training database does not contain all the surface terminations shown. The ACE model is also able to predict the stability of the reconstructed α-quartz (001) surface (Supplementary Fig. 3a, b). Whilst these surface energies can be computed with DFT, realistic amorphous surface energies are much more difficult to obtain, due to the required system sizes. Therefore, Fig. 4b validates the potential on 125 small-scale surface structures of amorphous SiO2, each containing 192 atoms. The surface models are created based on bulk structures from ref. 31; the latter had been generated in melt–quench simulations with different interatomic potentials and therefore span a range of energies. Regardless of the starting structure, the ACE model captures the surface energy for all slab models very well: the total RMSE is about 0.01 eV/Å2, and only a slight underestimation compared to DFT is seen. Moreover, there are no clear outliers although the various surface energies indicate a large diversity of the surface structures.

The amorphous surfaces shown are already very complex, but in reality they are often not flat as here. They have curvature, for example when occurring inside pores, and such complex structures can no longer be directly validated with DFT. Figure 4c therefore shows how well atomic environments in various porous amorphous structures are covered by the dataset. These structures were prepared by straining amorphous structures at elevated temperatures to the desired density. To validate the performance of the potential on this model, we show the linear extrapolation grade according to the max vol selection41,53. An extrapolation grade above 1 corresponds to atomic environments that have not been covered in the training database. This does not mean that the potential is no longer reliable, as there is a certain range of more or less reliable extrapolation, but as the extrapolation grade increases, non-physical behaviour and failure of the potential becomes more and more likely39. For all porous structures, regardless of density, we find that the maximum extrapolation value is less than 1. Thus, we observe no pronounced extrapolation in any of the cases considered, indicating an accurate description of the potential for a variety of curved surfaces. In Supplementary Fig. 3e, f, we demonstrate the effectiveness of the ACE potential in reconstructing artifical amorphous surfaces. To this end, we cut out a spherical cavity from an amorphous bulk structure. Upon heating, the amount of incorrectly arranged atoms near the surface (<3 Å) decreases over time, indicating that the surface atoms are undergoing significant reconstructions.

Elemental silicon
Whilst our ACE model is designed for the binary Si–O system, we show in Table 2 the performance for diamond-like elemental silicon compared to both DFT and experiment. The bulk modulus is very well reproduced, whereas the vacancy formation energy is underestimated by about 30%. The experimental surface energies are well recovered by the ACE model, but this may be partly due to serendipity, because the SCAN ground-truth data show poorer quality (Table 2). Moreover, the potential captures the reconstruction of the diamond-type silicon (100) surface (Supplementary Fig. 3c, d). We also computed the linear thermal expansion coefficient of diamond-type silicon in the quasi-harmonic approximation, finding almost perfect agreement between the ACE prediction and experiment; in particular, the unusual negative expansion coefficient below 130 K is reproduced (Supplementary Fig. 4). Finally, melt–quench simulations were performed to generate a-Si structures (Supplementary Fig. 5). The agreement with the experimental structure factor is as good as that for a GAP-18-generated structure from ref. 54. In addition, we are able to achieve lower quenching rates with the ACE than with the GAP, and for quenching rates as low as 1010 K s−1, we observed crystallisation.

Compared to Si-GAP-1830, we observe higher errors with respect to the referenced data (Table 2). This is not a principal shortcoming of ACE compared to GAP: it was already shown that it is possible to fit an ACE potential with similar numerical accuracy as Si-GAP-18 to the same training database34. Instead, the lower accuracy might be caused by the extension of the database to a second element, compared to the Si-GAP-18 one, and its strong focus on the SiO2 part of the configurational space. We assume that this causes, in turn, a less accurate description of the configurational space of the elemental species. Larger training databases might help to overcome this issue in the future. Due to the lower accuracy in reproducing the SCAN data, our potential has some shortcomings for higher-pressure structures: the bc8 phase is erroneously predicted to be stable at elevated pressure (Supplementary Fig. 6), and upon compressing a-Si we do not observe the eventual crystallisation that is described by Si-GAP-18 (Supplementary Fig. 7)55. We emphasise that very-high-pressure silicon phases were not the scope of the present work – instead, we focus on the accurate description of ambient-pressure silicon as a constituent part of mixed binary phases and nanostructures.

SiO and mixed silicon–silica systems
Whilst the results so far have served to demonstrate the usefulness of the approach – both in terms of development of datasets and the fitting within the ACE framework – we are now able to study an actual application problem. To this end, Fig. 5a shows structural models of SiO. Experimentally, amorphous SiO is obtained by deposition of SiO from the gas phase56. In contrast, we created our models by melt–quench simulations. SiO phases are known to be metastable with respect to Si and SiO2. For example, a recent DFT-based crystal-structure prediction study explored possible ordered phases of homogeneous SiO, and found that these are metastable compared to a mixture of crystalline Si and SiO2 57. We verified that our ACE potential similarly reproduces the metastability of the ambient-pressure phases (Supplementary Fig. 9), and that it accurately predicts relevant Si–SiO2 interface energies (Supplementary Fig. 10). In good agreement with these results, our melt-quenched structures show a clear segregation between a-Si-like (blue) and a-SiO2-like regions (red). With decreasing quench rate, the number of silicon grains decreases while their size increases. Figure 5b shows the structure factor, S(q), for the structure quenched with 5×1012 K s−1 (data for the other structures are shown in Supplementary Fig. 11), which agrees best with the experimental data from ref. 15. Figure 5c shows the ratio between the volume of the silicon grains divided by the interface area. In the approximation of spherical particles, the grain diameter is d = 6·Vgrains/Ainterface. From this we can estimate average grain diameters between 24 and 54 Å for the tested quench rates. These grain diameters agree very well with transmission electron microscopy measurements, which indicated diameters of 30–40 Å14.

Figure 5d shows the excess energies of the structures referenced to to α-quartz and diamond-type silicon. The SiO structures were relaxed by optimisation of the cell size as well as the atomic positions at 0 K. As experimental reference, we show the standard enthalpy of formation of SiO58. The structures generated using quench rates of 5×1012 and 2×1012 K s−1 have energies comparable to experiment. Indeed, we can even create structures that are energetically more favourable than in experiment, noting again that our procedure to produce the structures deviates significantly from the experimental one.

But is this really an improvement compared to existing, empirically fitted interatomic potentials? We tested the Munetoh potential59 and a charge-optimised many-body (COMB) potential29 for the same procedure to generate structural models of SiO. The Munetoh potential yielded a homogeneous structure without observable segregation into silicon and SiO, and the resulting structure factor (Supplementary Fig. 11c) deviates strongly from experiment. For the COMB potential, we observed pore formation at elevated temperatures, finally resulting in a strongly increased simulation-cell size. Therefore, we only equilibrated our best-matching structure at room temperature and analysed the change in structure factor (Supplementary Fig. 11f): again, we observed a strong deviation from experiment, indicating that the structure is very different from the ACE model prediction.

Crystallisation of silicon in amorphous SiO
Silicon in silicon monoxide is experimentally known to crystallise above 850 °C60. Figure 6 illustrates simulations of such crystallisation processes. The SiO structures shown in Fig. 5 were heated to 1400 K, causing the silicon-rich regions to melt while the silica matrix remained solid. The structures were then quenched to 1200 K within 20 ns (Fig. 6a). Through this cooling process, we noticed crystallisation in the silicon-rich regions of the SiO structure. Details are shown in Supplementary Fig. 12, indicating how crystallisation starts from two seeds, appearing shortly after each other and propagating throughout the structure. Before thermal treatment, all structures show nearly no sign of crystallinity, whereas afterwards, the structures with larger silicon grains do (Fig. 6c–f). This also affects the structure factor: those systems with no or only a small amount of crystallinity (quench rates: 1×1013 and 5×1012 K s−1) show structure factors that are still similar to the experimental structure factor of SiO (Fig. 6g), whereas those with larger amounts of crystalline silicon (1×1012 and 2×1012 K s−1) show distinct S(q) peaks that indicate crystallinity (Fig. 6h). This comparison allows us to exclude the occurrence of large amounts of crystalline silicon in silicon monoxide samples, given the experimental S(q) (cf. Fig. 6g and ref. 15).

Figure 6i shows that all structures are energetically more favourable after the thermal treatment. However, even though one might expect that the stronger crystallised structures gain more energy, we observe no direct connection between the energy gain and the level of crystallinity after heat treatment. The reason is likely that for the fast-quenched structure (1×1013 K s−1), the silicon-rich regions might be not perfectly arranged. Thermal treatment thus lowers the interface energy and the internal energy of the amorphous silicon region. For the crystallised structures there is, on one hand, an energy gain due to crystallisation but, on the other hand, an energy loss due to the higher interface energy between crystalline silicon and the silica matrix compared to that of a-Si and the silica matrix (cf. Supplementary Fig. 13). Based on these interface energies, we constructed a simple physically-based model (Methods) to calculate the energy gain of a spherical silicon inclusion in an amorphous silica matrix by crystallisation. This energy gain is shown in Fig. 6j for interface energies taken from four pairs of manually constructed interfaces. The difference between the models is rather small, indicating that the minimum radius of a spherical inclusion has to be around 9 Å to be energetically favourable in a crystalline state. We note that the manually constructed models do not have the same interface orientations as in the SiO structures and are also not perfectly relaxed, causing an overestimation of the interface energies – as seen, for example, in the difference between the interface energies of the SiO structures and of the manually constructed a-Si–a-SiO2 interface in Supplementary Fig. 13d. However, since for our model the difference between the interface energy between a-Si–a-SiO2 and c-Si–a-SiO2 is the only relevant quantity, we assume that these effects partially cancel.

Methods
Machine-learning potential fitting
We used two frameworks for fitting ML potential models. While constructing the referenced database, we used Moment Tensor Potentials38 with active learning53 as implemented in the MLIP package39. For the final potential fit, we used the nonlinear Atomic Cluster Expansion (ACE)33,34 as implemented in PACEMAKER35. For ACE, we tested a range of combinations of embeddings, and found the following to be suitable:
Ei = ϕ + ϕ + ϕf i, (1)
with ϕ i being atomic properties, which are expanded by the ACE basis functions (for details see ref. 33). The exponents of the embeddings include fractional exponents and higher integer powers of f i 2 1=8, 1=4, 3=8, 3=4, 7=8, 2. We found that especially fractions between 0 and 1 improved the behaviour of the potential. This approach goes beyond the previously suggested linear embedding (only the first term) and Finnis–Sinclair (the first two terms) type embedding34, and is referred to as “complex” embedding in Table 1.

For the expansion of the atomic properties ϕ i we used 600 basis functions with 5700 parameters. As radial basis we employed Bessel functions. A κ value of 0.01, which gives the ratio between force and energy weights value, was used during fitting. For optimisation we used the BFGS algorithm for 2000 steps.

DFT computations
All DFT computations were performed using VASP62,63 and the projector augmented-wave method64,65. We used the SCAN functional36 with an energy cutoff of 900 eV and a k-spacing of 0.23 Å−1. Surface were performed with dipole corrections. We note that these convergence parameters are optimised for silica; however, we found them to be also well converged for mixed phases and for pure silicon structures. Only for very-high-pressure silicon allotropes, a higher k-spacing would provide a relevant advantage; however, since these are not in the scope of the present work, we neglect these inaccuracies.

MD workflows
Simulation protocols were implemented using the atomic simulation environment (ASE)66 and the OVITO Python interface67. While optimisation and small-cell MD were partially performed with ASE, large-scale MD and statics simulations were carried out using LAMMPS68. The time step was 1 fs. For NVT simulations, we used a Nosé–Hoover thermostat with temperature damping constant of 100 fs; for NPT simulations, we added a Nosé–Hoover barostat with a pressure damping constant of 1000 fs.

For quench simulations we used the same protocol as in ref. 31. This protocol starts with a randomisation part at 6000 K for 10 ps (NVT). Afterwards the temperature is immediately reduced to 4000 K and held there for 100 ps (NPT, zero external pressure). From there the melt is quenched with different quench rates to 300 K (NPT, zero external pressure). At this temperature the structure is equilibrated for another 10 ps. In case of ‘hybrid’ simulations the quenches have been performed using the CHIK potential25 and afterward the structure has been equilibrated for another 20 ps with the ACE potential.

As input for the compression simulations we generated amorphous structures using this quenching protocol with a quench rate of 1013 K s−1 and only the ACE potential. The compression was performed under isostatic conditions. In each step, the pressure was initially increased by 1 GPa within 2.5 ps of simulation time, followed by equlibration over 2.5 ps at the new pressure. This procedure was iteratively repeated. Coordination numbers were determined after equilibration.

The aerogel structures were created by a similar protocol as in ref. 31. An initial structure was randomised at 6000 K for 10 ps, instantly cooled to 4000 K and kept there for 100 ps. From this temperature, the liquid was cooled to 300 K with a quench rate of 1013 K s−1. During the equilibration at 4000 K and up to half of the quenching process, the cells were additionally extended to the desired density.

The mixed structures were created using the same protocol as in ref. 31 for producing amorphous structures. The volume of the silicon grains was determined within OVITO by creating bonds between silicon and oxygen atoms (cut-off: 2 Å) and deleting all atoms which have more than one such bond. This deletes the whole silica matrix. The interface area and volume is then determined by using the Construct Surface Mesh modifier (Gaussian density method, resolution: 50, radius scaling: 100%, isovalue: 0.1) on the remaining atoms.

Phase diagram calculations
Thermodynamic integration was carried out as implemented in 43,44. We used the same simulation cells as in 50,000 atoms, as well as 80,000 atoms for the thermodynamic integration to calculate the temperature dependence of the free energy. We checked the influence of the system size on the thermodynamic integration and found it to be negligible. More details on the phase diagram calculations can be found in Supplementary Note 2.

Structure factors
Faber–Ziman structure factors were obtained by summation of the Fourier transformations of the partial radial distribution functions calculated with OVITO. The corresponding partial structure factors were weighted by atomic form factors taken from ref. 69. For the high-pressure structures, we used a cut-off radius of 20 Å for the radial distribution function, and analysed a single snapshot (without time averaging). For the SiO structures, we used a cut-off of 80 Å and an average over 10 snapshots.

Surface energies
In these calculations we consider only stoichiometric slabs. The surface energies, γ, were calculated as
γ = Eslab − N·Eref / A, (2)
where A is the total surface area (at the top and bottom of the slab combined), N is the number of particles in the slab, Eref is the bulk reference energy and Eslab is the potential energy of the slab. The slab energies for Table 2, Fig. 4a (relaxed), and Supplementary Fig. 3 have been calculated for DFT- and ACE-relaxed structures. For Fig. 4b, we used the ACE-relaxed structures also to determine the DFT single-point surface energy. The reference energy for the α-quartz surface is the energy of the optimised α-quartz unit cell per atom, and the reference energy for the diamond surface is obtained for the optimised diamond-type unit cell. The reference energy of the amorphous sample is the bulk energy of the same relaxed amorphous structure without surfaces.

Enthalpy and structural analysis at high pressure
The enthalpy H is given by
H(p) = E(V) + p(V)·V, (3)
where E is the internal energy, p is the pressure, and V is the volume. The volume dependence of the energy was determined by a Birch–Murnaghan fit to the energy–volume curve of each polymorph. p(V) was given by the corresponding derivative. The energy-volume curves were calculated by varying the volume by ±20% for α-quartz and coesite, by ±25% for stishovite and ±30% for all other phases. The corresponding structures were structurally optimised, allowing changes of the positions as well of the box shape, however keeping the volume fixed. For the analysis of the compression MD simulations, coordination numbers were determined by integrating over the first peak of the partial Si–O radial distribution function. The Si–O bond distances are given by the first peak position of the partial Si–O radial distribution function.

Interface energy model for the SiO crystallisation
We build an interface energy based model for a spherical inclusion of silicon with radius r within an amorphous SiO2 matrix. The energy difference between the crystallised silicon and the amorphous silicon is given by,
ΔE = 4πr2(γc-Si-a-SiO2 − γa-Si-a-SiO2) + 4πr3/3(Ea-Si − Ec-Si) / Vatom : (4)
where γc-Si-a-SiO2 and γa-Si-a-SiO2 are the interface energies between crystalline silicon and amorphous silica, and amorphous silicon and amorphous silica, respectively, Ea-Si is the energy of the corresponding amorphous silicon structure (cf. Supplementary Table IV), Ec-Si is the energy of crystalline silicon, and Vatom is the volume per atom. Details about the interface energies we used can be found in Supplementary Note 2.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors report the state of the art of machine learning technique applied to a series of Si–O system. Although underlying technique of the paper is already reported in ref 31, I positively recommend the ms. for publication in Nature Communications with minor revisions because the introduction of large-scale active learning is quite important and the application to a wide range of Si–O system including amorphous silicon monoxide (please see ref. 15; Hirata, A. et al. Atomic-scale disproportionation in amorphous silicon monoxide. Nat. Commun. 7, 11591 (2016)) is very impressive. I believe that the contents are suitable for a wide range of broad readers of Nature Communications because Si–O system is very important as the authors mentioned in the introduction part of the ms.\n1. I recommend the authors to refer the following article because this article reported the landmark study on the structure of amorphous silicon monoxide, although the content is controversial.\nGreaves, G. N. EXAFS and the structure of glass. J. Non-Cryst. Solids 71, 203–217 (1985).\n2. I strongly recommend the authors to calculate neutron structure factor for amorphous silica because neurons are sensitive to oxygen atoms. This is helpful supplemental data for a benchmark. Numerical data is available at the following web site.\nhttp://www.alexhannon.co.uk/DBindex.htm\n3. The authors may be interested in the following article.\nOnodera, Y. et al., Understanding diffraction patterns of glassy, liquid and amorphous materials via persistent homology analyses. J. Ceram. Soc. Jpn. 127, 853–863 (2019).", "reviewer_2": "The authors present the development and application of a machine-learned interatomic potential for Si-O systems. The results are scientifically sound and the paper is written in a compelling and thorough way that makes it easy to follow and enjoyable to read. While I am not 100% sure that the novelty of the work compared to previous works justifies the publication in such a high-impact journal, the importance of the material system combined with the timeliness of the approach does. For this reason, I recommend publishing this paper.\nComments:\n- To claim that the potential is suitable to study surfaces, I think it is critical to show that surface reconstructions and relaxations are properly captured. Fig. 4a only shows DFT single-point results for the unrelaxed and relaxed ACE structures. How far off are the relaxed DFT values? Moreover, does the potential capture reconstructions such as six-membered rings for the a-quartz (001) surface [https://doi.org/10.1039/B701176H] or the 2x1 reconstruction for (001) Si surfaces?\n- Some simple tests in the spirit of Fig. 4a for SiO2/Si interface energies would strongly support that the segregation of SiO2 and Si is properly captured by the potential.\n- A stunning result obtained with the GAP-18 for silicon is that it is able to describe the polyamorphic transition from LDA-Si to HDA-Si – even though no high-pressure amorphous structures were included in the fit data base [https://doi.org/10.1038/s41586-020-03072-z]. Moreover, the GAP-18 potential correctly describes E(V)-curves for the high-pressure st12 and bc8 phases [https://doi.org/10.1103/PhysRevX.8.041048], in contrast to the ACE potential. Since the GAP-18 data base is part of the data base for the present potential, but the major shortcoming of the present potential is that high-pressure phases of pure Silicon are poorly described, I wonder where this loss of transferability comes from. Is this a consequence of the GAP vs ACE framework or does the other fit data impede the transferability? Does the current potential capture the polyamorphic transition? In my opinion, a more through discussion of this aspect might help other researchers in the quest of designing truly general-purpose potentials.\n- How do the results of Fig. 5 depend on the system size?\n- Ref. 3 in the SM “${\\mathrm{SiO}} {2}$” -> SiO2\n- The “construct surface mesh” modifier results in OVITO should depend on the probe sphere radius. For the sake of reproducibility, I recommend adding this parameter to the methods section.", "reviewer_3": "The authors present an ACE Si-O model trained using active learning on DFT(SCAN) reference data. One of the prominent contributions is that the reported ML potential is two orders of magnitude faster than GAP model at the expense of a small reduction in accuracy. The article presents a thorough analysis of the materials’ properties using such model. In general, it provides good results, making it of interest for the materials’ science users. Nevertheless, this work only focuses on the analysis of the models’ characterization and not on a specific problem that, for example, reveals a new understanding on the dynamics of Si-O material.\nSome further comments are presented below:\n-When citing the supplementary material, please state the specific section.\n-Table I is a good overview of the accuracy, but what is more informative to assess accuracy are the forces, this because MD simulations only used forces to sample the phase space. Hence, I would recommend to also include force errors in this table.\n-Fig. 2 Phase diagram: Here, we can see that qualitatively, the ACE model can recover the phase diagram, nevertheless, it shows a severe deficiency predicting the transition liquid <-> solid of ~400K in the 0-1GPa region. This is a common issue in ML potentials trained on undersampled liquid phases. Meaning, the ML mostly learn the solid phases. The same happens with water models, where the simulations tend to generate icy phases instead of the proper liquid state.\n-The article presents an interesting manual-fine tuning of the model via manipulating the different penalties on the loss function for the different phase states. This corrects some of the boundaries in the phase diagram. This could be a valid non-machine-learned way to compensate undersampled phases. Although, a full analysis on the rest of physical properties is required to assess if the model continues to be physically valid.\n-The authors stated: “We thus conclude that the quantitative deviation seen in Fig. 2 is due to the inaccuracy of the underlying exchange–correlation functional, rather than indicating a shortcoming of the ACE approach.” In general, a more elaborated note should be added to discuss the accuracy of the reference DFT data, this because it is well known that, even though the SCAN functional provides accurate energies, its gradients (i.e. forces) are problematic and highly depend on the stability of the implementation. Thereby, this could be a potential source of predictive error beyond the ML model.\n- Table II presents nice summary on the performance of the model. As expected for most of ML potentials, mechanical properties (also Fig. S3) are accurately described by the ML models (GAP and ACE) compared to DFT reference data. Nevertheless, we can see that ACE model consistently underestimate energy related predictions, contrasting to the high accuracy of the GAP model. Most likely, this is related to limitation on the learning capacity of the chosen model. As highlighted by the authors, here, the only comparison that should be considered is ACE to DFT(SCAN), while the apparent agreement of ACE with experiment is due to serendipity. \n-The authors stated “How- ever, for properties such as the bulk modulus and surface energies of diamond, the agreement with experiment is better than that of GAP-18.” This should be removed given that is a misleading comment. The apparent agreement of ACE results with experiment is due to serendipity. There is no fundamental reason why ACE predictions are comparable with experiment. IF ACE was generating predictive simulations, it should recover SCAN predictions."}, "author_response": "Response to Reviewers’ Comments – Manuscript NCOMMS-23-30101\nWethankallthreereviewersfortheirhelpfulandconstructivefeedback. Inthefollowing,\nwequotethereviewers’reportsinfull,andourresponsesareinterspersedinblue. Action\ntaken is indicated in red.\nReviewer #1 (Remarks to the Author):\nThe authors report the state of the art of machine learning technique applied to a series\nof Si–O system. Although underlying technique of the paper is already reported in ref\n31,Ipositivelyrecommendthems.forpublicationinNatureCommunicationswithminor\nrevisions because the introduction of large-scale active learning is quite important and\nthe application to a wide range of Si–O system including amorphous silicon monoxide\n(please see ref. 15; Hirata, A. et al. Atomic-scale disproportionation in amorphous silicon\nmonoxide. Nat. Commun. 7,11591(2016))isveryimpressive. Ibelievethatthecontents\nare suitable for a wide range of broad readers of Nature Communications because Si–O\nsystem is very important as the authors mentioned in the introduction part of the ms.\nResponse: Wethankthereviewerfortheirpositiverecommendation,andforemphasising\nthe importance of the Si–O system studied in this work.\n1. I recommend the authors to refer the following article because this article reported the\nlandmark study on the structure of amorphous silicon monoxide, although the content is\ncontroversial. Greaves, G. N. EXAFS and the structure of glass. J. Non-Cryst. Solids 71,\n203–217 (1985).\nResponse: We thank the reviewer for pointing out this reference.\nAction taken: We now cite the mentioned reference in the Introduction (new ref. 13, p.\n1).\n2. I strongly recommend the authors to calculate neutron structure factor for amorphous\nsilica because neurons are sensitive to oxygen atoms. This is helpful supplemental data\nfor a benchmark. Numerical data is available at the following web site.\nhttp://www.alexhannon.co.uk/DBindex.htm\nResponse: This is a good idea, and indeed provides an additional way of validating the\npotential.\nAction taken: We added Supplementary Fig. 14 (p. S32), which contains neutron and\nx-ray structure factors for amorphous structures generated by the ACE model, as well as\nby the hybrid protocol described in Ref. 31. We added the following comment to the\nDiscussion section:\n“ThepotentialalsounderestimatestheFSDPheightofamorphousSiO (Supplemen-\n2\ntary Fig. 14), as already observed for our earlier GAP model.31 Further research is\nrequired to identify the origin of this underestimation.” (p. 9)\n1\n3. Theauthorsmaybeinterestedinthefollowingarticle. Onodera,Y.etal.,Understanding\ndiffraction patterns of glassy, liquid and amorphous materials via persistent homology\nanalyses. J. Ceram. Soc. Jpn. 127, 853–863 (2019).\nResponse: We thank the reviewer for sharing this article. The methodology is highly\ninterestingandwehopethatwecanusethisinfuturewithourstructuremodelstoachieve\na better understanding of the structure of amorphous silica.\nThe authors present the development and application of a machine-learned interatomic\npotential for Si-O systems. The results are scientifically sound and the paper is written in\nacompellingandthoroughwaythatmakesiteasytofollowandenjoyabletoread. While\nI am not 100% sure that the novelty of the work compared to previous works justifies the\npublicationinsuchahigh-impactjournal,theimportanceofthematerialsystemcombined\nwith the timeliness of the approach does. For this reason, I recommend publishing this\npaper.\nResponse: Thank you very much.\nComments:\n- To claim that the potential is suitable to study surfaces, I think it is critical to show that\nsurface reconstructions and relaxations are properly captured. Fig. 4a only shows DFT\nsingle-point results for the unrelaxed and relaxed ACE structures. How far off are the\nrelaxed DFT values? Moreover, does the potential capture reconstructions such as six-\nmembered rings for the a-quartz (001) surface [https://doi.org/10.1039/B701176H] or\nthe 2x1 reconstruction for (001) Si surfaces?\nResponse: Surfacerelaxationsareindeedanimportantpoint,whichwerealisewedidnot\ncoverinsufficientdetailintheinitialsubmission. Regardingthefirstpoint,inFig. 4a,the\nrelaxation of the α-quartz surfaces by DFT leads to only minor differences in the surface\nenergy,viz.0.14eV/A˚2 for(001)andatmost0.04eV/A˚2 forallothersurfacesconsidered.\nWedidnotobservesignificantstructuralchanges(reconstructionsetc.) duringDFTrelax-\nation. Regarding the second point, in both cases the surface reconstructions mentioned\nby the reviewer are correctly predicted as being energetically more favourable than the\ncleaved surfaces.\nFurthermore,giventhatastrongfocusofourpotentialisonamorphoussilicasurfaces,we\nalsolookedatsurfacereconstructionsofthese. Forthispurpose,wecutoutaspherefrom\na bulk amorphous structure, which introduces a highly unfavourable surface with many\nwronglycoordinatedatoms. WehaveequilibratedthissurfaceandshowinSupplementary\nFig. 3 (p. S22) that the number of wrongly coordinated atoms is strongly decreasing over\nthesimulationtime. Thisalsoindicatesapromisingbehaviourforsurfacereconstructions\nof arbitrary amorphous surfaces.\nAction taken: WeaddedSupplementaryFig.3(p.S22),whichshowsthereconstructions\nas well as the surface energies for the (001) α-quartz and the (100) silicon surface, and\nwe updated Fig. 4a in the main text (p. 6). Moreover, we added three comments to the\nresults part:\n“The ACE model is also able to predict the stability of the reconstructed α-quartz\n(001) surface (Supplementary Fig. 3a–b).” (p. 5)\nand\n“In Supplementary Fig. 3e–f, we demonstrate the effectiveness of the ACE potential\n3\nin reconstructing artifical amorphous surfaces. To this end, we cut out a spherical\ncavity from an amorphous bulk structure. Upon heating, the amount of incorrectly\narranged atoms near the surface (< 3 A˚) decreases over time, indicating that the\nsurface atoms are undergoing significant reconstructions.” (p. 6)\nas well as\n“Moreover, the potential captures the reconstruction of the diamond-type silicon\n(100) surface (Supplementary Fig. 3c–d).” (p. 6)\n- Some simple tests in the spirit of Fig. 4a for SiO2/Si interface energies would strongly\nsupport that the segregation of SiO2 and Si is properly captured by the potential.\nResponse: This is a good point. We constructed interfaces between amorphous and\ncrystallinesilicaandsiliconbymergingmodelsofthesestructurestogether. Wecalculated\ninterfaceenergieswithDFTandtheACEpotential. Theinterfaceenergiespredictedbythe\nACE potential agree well with DFT, especially given the wide range of interface energies\ncovered in our tests.\nAction taken: We added Supplementary Fig. 10 (p. S28), which shows the interface\nenergies for a wide range of silica–silicon interfaces. We show these results both for\nrelaxed and unrelaxed interfaces, and for DFT as well as for ACE.\n- A stunning result obtained with the GAP-18 for silicon is that it is able to describe the\npolyamorphictransitionfromLDA-SitoHDA-Si–eventhoughnohigh-pressureamorphous\nstructureswereincludedinthefitdatabase[https://doi.org/10.1038/s41586-020-03072-\nz]. Moreover, the GAP-18 potential correctly describes E(V)-curves for the high-pressure\nst12 and bc8 phases [https://doi.org/10.1103/PhysRevX.8.041048], in contrast to the\nACEpotential. SincetheGAP-18databaseispartofthedatabaseforthepresentpotential,\nbut the major shortcoming of the present potential is that high-pressure phases of pure\nSilicon are poorly described, I wonder where this loss of transferability comes from. Is\nthis a consequence of the GAP vs ACE framework or does the other fit data impede the\ntransferability? Does the current potential capture the polyamorphic transition? In my\nopinion, a more through discussion of this aspect might help other researchers in the\nquest of designing truly general-purpose potentials.\nResponse: Regardingthefirstquestionaboutwherethelossoftransferabilitycomesfrom,\nwe believe that it is caused by the additional fitting data. The ACE framework is unlikely\nto be the reason, as it has already been successfully used to fit the GAP-18 database\n[Lysogorskiy et al., npj Comput. Mater., 7, 97 (2021)] with similar accuracy as the GAP.\nThe lower accuracy in the silicon part is rather due to the extension to two elements and\nthe focus on SiO in the overall database.\nfrom low density to high density amorphous silicon. However, it does not show the\ncrystallisationaround20GPaobservedforGAP-18. Wewouldliketopointoutthatthisis\nnotthepurposeofthepotential–therearegoodpotentialsavailableforpuresilicon,and\nfor a study of only the elemental system those are likely to be the preferred option at this\npoint. Ofcourse,inthefuture,itwouldbeveryinterestingtoextendthecurrentpotential\n4\nto describe high-pressure silicon, arriving at the “truly general-purpose” type of potential\npointed out by the reviewer.\nAction taken: We revised the discussion regarding the accuracy for silicon and extended\nit (in line with the referee’s comments):\n“ThisisnotaprincipalshortcomingofACEcomparedtoGAP:itwasalreadyshown\nthat it is possible to fit an ACE potential with similar numerical accuracy as GAP-\n18 to the same training database.35 Instead, the lower accuracy might be caused\nby the extension of the database to a second element, compared to the GAP-18 one,\nand its strong focus on the SiO part of the configurational space. We assume that\n2\nthis causes, in turn, a less accurate description of the configurational space of the\nelemental species. Larger training databases might help to overcome this issue in\nfuture. ” (p. 7)\nWeaddedSupplementaryFigure7(p.S25)containinganillustrationofthepolyamorphic\ntransition with the ACE and analysis regarding coordination numbers as well as atomic\nvolume as dependent on the pressure. We also added a short sentence regarding this in\nthe main part:\n“[...],anduponcompressinga-Siwedonotobservetheeventualcrystallisationthat\nis described by Si-GAP-18 (Supplementary Fig. 7).62” (p. 7)\n- How do the results of Fig. 5 depend on the system size?\nResponse: The results of Fig. 5 are nearly independent of the system size. Only for the\nratio of the interface area between Si and SiO we are able to observe small differences\n2\nwith different system sizes. However, for all other properties, differences are negligible.\nAction taken: We added Supplementary Fig. 8 (p. S26) to the manuscript, which shows\nthestructurefactors,interface-to-volumeratiosandenergiesofformationforsystemsizes\nof 43,904, 221,184, 432,000, and 877,952 atoms, respectively.\n- Ref. 3 in the SM \"${\\mathrm{SiO}} {2}$\" − > SiO\nAction taken: We corrected the formatting on p. S38.\n- The “construct surface mesh” modifier results in OVITO should depend on the probe\nsphere radius. For the sake of reproducibility, I recommend adding this parameter to the\nmethods section.\nResponse: We added the specifications for the ConstructSurfaceMesh modifier in the\nMethodssection. Wechangedtheprocedureslightlycomparedtotheinitialsubmissionto\nobtain an improvement in the calculation of the interface area. This improvement is only\nminor and the overall conclusions are similar.\nAction taken: We extended the Methods section (p. 10) and updated Fig. 5c (p. 7).\nThe authors present an ACE Si-O model trained using active learning on DFT(SCAN)\nreference data. One of the prominent contributions is that the reported ML potential\nis two orders of magnitude faster than GAP model at the expense of a small reduction in\naccuracy. The article presents a thorough analysis of the materials’ properties using such\nmodel. In general, it provides good results, making it of interest for the materials’ science\nusers. Nevertheless,thisworkonlyfocusesontheanalysisofthemodels’characterization\nand not on a specific problem that, for example, reveals a new understanding on the\ndynamics of Si-O material.\nResponse: Wethankthereviewerfortheirpositivecomments. Wealsoappreciatethefact\nthattheyarechallengingustoprovideclearerinsightsintoaspecificapplicationproblem,\nwhich we have now addressed during revision.\nSpecifically,wehavenowcarriedoutnewanalysesofthedynamicalcrystallisationprocess\nin SiO structures. We achieved this by heating the SiO structures above the melting point\nof silicon and then quenched it rapidly by 200 K within 20 ns. During the quench, the\nsilicon-rich regions partially crystallised. This process is characterised in the newly added\nFig. 6. Forexample,webuiltasimplethermodynamicmodeltounderstandtheconditions\nunder which the amorphous Si crystallises into the SiO structures.\nAction taken: We added a new figure (Fig. 6, p. 9) along with a paragraph providing\ninformation about it. Additionally, Supplementary Figures 12 (p. S30) and 13 (p. S31)\nhavebeenadded,offeringdetailsabouttheprocessofcrystallisationandinterfaceenergies.\nSupplementary Note 2 contains a description of our thermodynamic model for the Si\ncrystallisation.\nSome further comments are presented below:\n-When citing the supplementary material, please state the specific section.\nAction taken: We added the relevant specific references.\n-TableIisagoodoverviewoftheaccuracy,butwhatismoreinformativetoassessaccuracy\nare the forces, this because MD simulations only used forces to sample the phase space.\nHence, I would recommend to also include force errors in this table.\nResponse: We agree that the forces are crucial for MD simulations. Whilst we had\npreviously shown them in the Supplementary Material (Supplementary Tab. 4), we now\ninclude them in the main part following the referee’s suggestion.\nAction taken: We extended Table 1 by including force errors (p. 4).\n-Fig. 2 Phase diagram: Here, we can see that qualitatively, the ACE model can recover\nthe phase diagram, nevertheless, it shows a severe deficiency predicting the transition\nliquid< − >solidof∼400Kinthe0-1GParegion. ThisisacommonissueinMLpotentials\ntrained on undersampled liquid phases. Meaning, the ML mostly learn the solid phases.\nThe same happens with water models, where the simulations tend to generate icy phases\n6\ninstead of the proper liquid state.\nResponse: We believe that the melting-point underestimation is likely caused by the\nunderlying DFT exchange-correlation functional. We note that we calculated the free\nenergy of the liquid and crystalline phases separately using thermodynamic integration,\ni.e., without a direct simulation of the transition of the two phases, thereby ruling out\ndynamic effects. However, undersampling may indeed result in an imprecise description\nof the liquid phase.\nTo estimate the error of the ACE model for the liquid phase, one can compute energy\nerrors for snapshots from MD simulations which are part of the ACE-MD test set for\namorphous structures (cf. Supplementary Table 4). The energy RMSE for the 3,000 K\nstructures is approximately 3.7 meV/atom. Assuming that the energy of the liquid phase\nisoverestimatedbyACE,correctingitby3.7meV/atomwouldproducethephasediagram\ndepicted in Supplementary Fig. 19 (p. S36). However, the melting point of 2,260 K\nthus obtained still deviates from experiment by over 260 K – which, in our view, can\nbe attributed to the underlying exchange-correlation functional data.\nActiontaken: WeaddedSupplementaryFig. 19(p.S36)aswellasthefollowingparagraph\nto Supplementary Notes 2:\n“Another point is that the melting point is significantly overestimated by the ML\npotential at low pressures. This can, in principle, be due to two reasons, namely the\nexchange–correlation functional and an undersampling of the liquid phase space.\nThe latter would likely correspond to a high error in the prediction of the energy\nof the liquid phase. To assess this, we analyse the 3000 K ACE-MD test set (see\nSupplementary Table III) for the liquid phase, which has an energy RMSE error of\n3.7meV/atom. Whenweassumethattheenergyoftheliquidphaseisoverestimated\nbythisvaluebythepotential,andcorrectthephasediagramaccordingly,weobtain\nthe phase diagram shown in Supplementary Fig. 19. In this phase diagram, the\nmelting point is still overestimated by 260 K, which thus indicates that (at least\npartially) the exchange-correlation functional is likely the reason for the deviation\nof the predicted melting point from experiment.” (p. S12)\n-The article presents an interesting manual-fine tuning of the model via manipulating the\ndifferent penalties on the loss function for the different phase states. This corrects some\nof the boundaries in the phase diagram. This could be a valid non-machine-learned way\nto compensate undersampled phases. Although, a full analysis on the rest of physical\nproperties is required to assess if the model continues to be physically valid.\nResponse: This is a very interesting point indeed. We would like to point out, for full\nclarity, that our Supplementary Fig. 1 analyses the uncertainty of the phase diagram by\nshifting the free energies of certain phases to higher values as a subsequent correction.\nThe interatomic potential is not directly modified in this study. Nevertheless, in future\nwork, it would be interesting to fine-tune the potential itself to see whether that might\nimprove agreement with experimental results.\nActiontaken: WerevisedthecaptionofSupplementaryFig.1(p.S20)tomakeclearthat\nthe potential model itself remains unchanged, and that we added the penalty to the free\nenergies predicted using the model.\n-The authors stated: “We thus conclude that the quantitative deviation seen in Fig. 2\nis due to the inaccuracy of the underlying exchange–correlation functional, rather than\nindicatingashortcomingoftheACEapproach.”Ingeneral,amoreelaboratednoteshould\nbe added to discuss the accuracy of the reference DFT data, this because it is well known\nthat, even though the SCAN functional provides accurate energies, its gradients (i.e.\nforces)areproblematicandhighlydependonthestabilityoftheimplementation. Thereby,\nthis could be a potential source of predictive error beyond the ML model.\nResponse: We thank the reviewer for pointing this out. Indeed, even though the SCAN\nfunctional provides accurate energies and significantly outperformed other functionals in\nour tests for silica (and therefore was our DFT approach of choice), there might still be\nremaining errors for the forces.\nAction taken: We added the sentence\n”This is an example of the more general problem that any issues with the ground-\ntruth method (e.g., numerical instabilities) will translate into the ML model.”\non p. 3 to extend the discussion regarding possible problems with underlying data.\n- Table II presents nice summary on the performance of the model. As expected for most\nof ML potentials, mechanical properties (also Fig. S3) are accurately described by the\nML models (GAP and ACE) compared to DFT reference data. Nevertheless, we can see\nthat ACE model consistently underestimate energy related predictions, contrasting to the\nhigh accuracy of the GAP model. Most likely, this is related to limitation on the learning\ncapacity of the chosen model. As highlighted by the authors, here, the only comparison\nthat should be considered is ACE to DFT(SCAN), while the apparent agreement of ACE\nwith experiment is due to serendipity. The authors stated “How- ever, for properties such\nas the bulk modulus and surface energies of diamond, the agreement with experiment is\nbetter than that of GAP-18.” This should be removed given that is a misleading comment.\nTheapparentagreementofACEresultswithexperimentisduetoserendipity. Thereisno\nfundamental reason why ACE predictions are comparable with experiment. IF ACE was\ngenerating predictive simulations, it should recover SCAN predictions.\nResponse: We agree that it seems like the learning capacity of the model reached a\nlimit, given the broad part of configurational space which is covered by the database.\nWe acknowledge that the comment regarding the agreement with experiment may have\nbeeninadvertentlymisleading,andwerevisedthemanuscriptaccordinglytoensureafair\ncomparison to GAP-18.\nActiontaken: Weextendedthediscussionofshortcomingsofthepotentialforsiliconand\npossible reasons for this (p. 7). We removed the statement mentioned by the referee.\nIn conclusion, we thank all three referees again for their constructive comments and for\ntheir helpful suggestions. We are convinced that the manuscript has improved further\nbased on their feedback.\n8"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have answered all my questions, and made the corresponding changes. Now I recommend the revised manuscript for publication in Nature Communications.", "reviewer_2": "The authors addressed all the points I raised very thoroughly. I recommend the publication in its present form. Excellent work!"}, "author_response": ""}]

Decision letter