Accurate and efficient potential energy surface (PES) models are crucial for large-scale simulations of materials, especially for complex systems like zeolites, where the structural diversity poses significant challenges for traditional methods. While machine learning (ML) approaches, such as neural network potentials (NNPs), offer promising solutions for achieving high accuracy and transferability, the robustness of these models for reactive phase transformations of zeolites remains largely unexplored. This research addresses this gap by training reactive SchNet NNPs for silica, encompassing the wide range of structures present in zeolites, from low-density zeolites to high-pressure polymorphs. The goal is to create a model capable of accurate and transferable simulations of siliceous zeolites, including their reactive phase transformations, with ab initio accuracy. The improved accuracy of these NNPs is compared against state-of-the-art analytical force fields and tight-binding DFT methods. The successful training of such an NNP provides a valuable computational tool for the design and synthesis of zeolites.

Previous studies have demonstrated the accuracy of SchNet NNPs in modeling energy and forces. However, their transferability and accuracy for material science questions, such as diffusion, phase stability, phonon properties, and reactive phase transitions of zeolites, remain largely unknown. Several studies have employed Multilayer Perceptrons (MLPs) to model the PES of silica using polymorphs, surface models, and amorphous configurations. However, none of these studies incorporated the structural diversity of zeolites and their reactive phase transformations in their training data. This work aims to overcome these limitations by training NNPs on a significantly more diverse dataset that includes the vast structural diversity of zeolites, offering a far more general model of the silica PES.

The methodology involves a multi-stage computational workflow (Figure 1). It begins with selecting a structurally diverse subset of zeolites from the Deem database using farthest point sampling (FPS) with SOAP as the similarity metric. This subset, along with dense silica polymorphs, surface models, and amorphous silica structures, undergoes unit cell deformations and molecular dynamics (MD) simulations at both low and high temperatures. DFT single-point calculations are performed on an FPS-sparsified set of configurations. An initial ensemble of NNPs is trained and iteratively refined using active learning. This involves structure optimizations of the Deem and IZA databases, extrapolation detection using a query-by-committee approach, and MD simulations of high-temperature and high-pressure phase transitions (melting of β-cristobalite, zeolite amorphization of LTA and SOD). The process continues until no extrapolation is detected. Two reference DFT methods, PBE + D3 and SCAN + D3, are employed for training. The resulting NNPs (NNPpbe and NNPscan) are then used for reoptimization of the Deem and IZA databases, simulations of glass melting and zeolite amorphization (including FAU, not in the training set), and predictions of equilibrium structures and vibrational properties. The accuracy of the NNPs is compared with experimental data and other PES approximations (SLC potential, ReaxFF force field, GFN0-XTB). The dataset generation includes selecting zeolites from the Deem database via FPS, performing unit cell deformations and MD simulations, and extracting diverse structures using FPS. High-energy configurations are sampled using MD simulations for melting and annealing of β-cristobalite and zeolite amorphization. Single-point calculations at the PBE + D3 level are performed for initial NNP training. NNP refinement involves extrapolation detection and iterative retraining until no extrapolation is detected. The Deem and IZA databases are also optimized, and additional simulations are conducted to enhance data diversity. Finally, NNPscan (trained on SCAN+D3 data) is chosen as the reference due to its superior performance over PBE+D3 for equilibrium properties and reactive transformations.

The trained NNPs exhibit remarkable accuracy, with RMSEs below 4.2 meV/atom for energies and 0.070 eV/Å for forces in equilibrium structures, showing good generalization capabilities. This accuracy is an order of magnitude better than analytical potentials (SLC, ReaxFF) and GFN0-XTB. The NNPs successfully model reactive phase transitions, including glass melting at 4800 K and zeolite amorphization (LTA and SOD). Even for FAU, not included in the training set, the NNPs produce qualitatively correct results, demonstrating transferability. Although extrapolation is observed at densities above 2.2 g/cm³, the energy errors remain significantly lower than those of other methods. The reoptimization of the Deem database using NNPs revealed >20,000 additional hypothetical zeolites within the thermodynamically accessible range. The NNPs demonstrate good agreement with DFT and experimental results for equilibrium structures, phonon properties, and phase transformations under extreme conditions. Table 1 summarizes the RMSE and MAE of energies and forces for NNPscan, eNNPscan, SLC, ReaxFF, and GFN0-XTB compared to SCAN + D3. Figure 2 illustrates the energy error distributions. Figure 5 shows the energetics of Stone-Wales defect formation calculated using DFT (B3LYP) and NNPscan, highlighting the accuracy of the NNP in modeling bond-breaking events. The density range of the simulated zeolite collapse was 1.3 to 2.4 g cm⁻³. At densities below 2.2 g cm⁻³, MD simulations showed bond-breaking events in LTA (2.1–2.2 g cm⁻³) and no bond cleavage in FAU. For both zeolites, the NNP energies and forces showed no extrapolation and agreed well with DFT results at densities <2.2 g cm⁻³. Only further compression to artificially high densities up to 10% beyond silica glass density resulted in NNP extrapolation. However, the difference between NNP and DFT energies was even in the extrapolation region at least three times lower (<40 meV atom⁻¹) than the RMSEs of the other PES approximations (e.g., 136 meV atom⁻¹ for ReaxFF).

The results demonstrate that SchNet NNPs achieve accuracy comparable to other state-of-the-art MLPs, not only for near-equilibrium structures but also for high-energy bond-breaking scenarios. The significantly improved accuracy compared to analytical force fields and tight-binding DFT underscores the effectiveness of this approach. The inclusion of the large structural diversity of zeolites in the training data contributes to the high transferability and generalizability of the model. The observed extrapolation at higher densities provides a clear path for future model refinement through active learning. The reoptimization of the Deem and IZA databases using these highly accurate NNPs is a significant contribution, providing a refined dataset for future ML studies on structure-stability-property correlations. This research contributes significantly to the in silico design and discovery of zeolites. The future extension of this model to include heteroatoms (like Al) or guest molecules (like water) promises to enable more realistic modeling of zeolites under synthesis and operating conditions.

This study successfully trained highly accurate and transferable SchNet NNPs for modeling siliceous zeolites, achieving an order of magnitude improvement over other PES approximations. The NNPs accurately predict equilibrium properties, vibrational properties, and reactive phase transitions. The reoptimized Deem and IZA databases provide valuable resources for future research. Future work should focus on extending the model to include heteroatoms and guest molecules for a more comprehensive representation of zeolites under realistic conditions.

The study acknowledges that the NNPs cannot outperform the reference DFT method used for training data generation. The zeolite amorphization simulations, while demonstrating bond-breaking events, are simplified models and may not fully capture the complexity of real-world thermal zeolite collapse. Extrapolation is observed at higher densities, though the errors remain lower than those of other methods, indicating areas for future model improvement and dataset expansion.