Accurate and transferable neural network potentials for reactive simulations of siliceous zeolites

The study addresses whether SchNet-based neural network potentials (NNPs) can accurately and transferably model the potential energy surface (PES) of silica across the broad structural diversity of zeolites, including reactive transformations. While SchNet is known to model energies and forces well, its robustness for materials science tasks such as diffusion, phase stability and transitions, phonon properties, and especially reactive phase transformations in zeolites has been unclear. Previous silica machine-learning potentials (MLPs) did not capture the extensive diversity of zeolite frameworks and their reactive transformations. The central aim is to train reactive SchNet NNPs for silica, spanning low-density zeolites to high-pressure polymorphs, and to use an ensemble for active learning to iteratively refine coverage of relevant configurations, enabling large-scale simulations with ab initio accuracy. The work also reoptimizes large zeolite databases to inform future machine learning on structure–stability–property relationships.

Prior studies highlighted strong performance of dispersion-corrected PBE and SCAN DFT for silica; PBE+D3 and PBE+MBD align well with experimental zeolite structures and phase-transition enthalpies among GGAs, with no prior SCAN benchmarks for siliceous zeolites. Here, SCAN+D3 shows slightly better performance than PBE+D3 for structures (e.g., density MAD 0.2 vs 0.3 Si nm−3) and good energetics (MAD ~2–3 kJ mol−1 versus experiment). For the Stone–Wales defect in a silica bilayer, PBE+D3 slightly underestimates and SCAN+D3 slightly overestimates barriers relative to B3LYP (up to −0.67 eV and +0.59 eV, respectively). Prior silica NNPs included limited polymorphs/surfaces and smaller cells (<144 atoms), with higher force errors (~200 meV Å−1), limiting generality. Other MLP frameworks (e.g., moment tensor, Gaussian approximation potentials) have reported meV/atom-level energy and ~0.1–0.3 eV Å−1 force errors for various materials, but comprehensive, reactive modeling across zeolites remained lacking. Analytical force fields (SLC, ReaxFF) and tight-binding DFT (GFN0-XTB) are widely used but show significantly higher errors for silica PES modeling compared with modern NNPs.

Computational workflow: (i) selection of a structurally diverse subset of hypothetical zeolites from the Deem database (>330k frameworks) using farthest point sampling (FPS) with SOAP kernel-based similarity; (ii) sampling of low- and high-energy PES regions via unit cell deformations and molecular dynamics (MD) using selected zeolites, dense silica polymorphs, surface models, and amorphous silica; (iii) DFT single-point calculations on FPS-sparsified configurations; (iv) initial SchNet NNP ensemble training; (v) iterative active learning using query-by-committee to detect extrapolation during structure optimizations and high-T/P MD, followed by DFT labeling and retraining, until no extrapolation detected; (vi) final NNP-level applications: reoptimization of Deem and IZA databases, glass melting, zeolite amorphization (including frameworks not in training), and prediction of equilibrium and vibrational properties. Dataset generation: Initial DFT dataset (PBE+D3 single points) included 10 hypothetical zeolites selected by FPS (plus a detected hypothetical silica bilayer and an α-quartz (001) surface model), five existing zeolites (CHA, SOD, IRR, MVY, MTF), and six silica polymorphs (α-quartz, α-cristobalite, tridymite, moganite, coesite, stishovite). All were optimized at zero pressure (PBE+D3), followed by 210 lattice deformations each. Additional near-equilibrium sampling used 10 ps ReaxFF MD at 600 and 1200 K per system, with 200 diverse structures selected per trajectory via FPS. High-energy sampling employed ReaxFF MD for melting/simulated annealing of β-cristobalite (2×2×2 supercell): density scaled from 2.3 to 2.2 g cm−3, optimized, equilibrated 100 ps at 6000 K, then cooled to 3000 K in three steps with 100 ps at each, plus an extra 100 ps at 3000 K. FPS selected 1000 diverse configurations. Ten 3000 K snapshots were quenched (PBE+D3) at constant volume; lowest-energy amorphous structure optimized at zero pressure; 210 lattice deformations applied. For zeolite amorphization (ZA) mimicking thermal collapse, 8 hypothetical zeolites were equilibrated at 1200 K for 100 ps (ReaxFF), then volume contracted in 10 steps to reach 2.2 g cm−3, with 100 ps equilibration after each step; FPS selected 1000 diverse structures per zeolite. PBE+D3 single points provided energies and forces for initial NNP training (ensemble of six). Active learning and refinement: Extrapolation detection via committee: simulations with one leading NNP, with single-point predictions evaluated by the remaining five NNPs. If deviations exceeded 10 meV atom−1 (energy) or 750 meV Å−1 (forces) from the ensemble average, corresponding PBE+D3 labels were added, and the ensemble retrained. To boost equilibrium coverage, Deem (>331k structures) and IZA (235 frameworks) databases were reoptimized at zero pressure. β-cristobalite was equilibrated at 4800 K for 1 ns to augment liquid silica sampling. ZA dataset extended using the same protocol for LTA and SOD. Final PBE+D3 dataset comprised ~33k structures with up to 400 atoms per cell. SCAN+D3 single points were computed on the final database to train NNPscan. Training used energies and forces; test RMSE ~4.7 meV atom−1 (energy) and 147 meV Å−1 (forces). Reference methods and models: DFT with PBE/SCAN plus dispersion corrections (D3 and MBD tested in benchmarks; SCAN+D3 selected as primary reference). SchNet NNPs trained at PBE+D3 (NNPpbe) and SCAN+D3 (NNPscan) levels; ensembles of six (eNNPpbe, eNNPscan). Comparative methods included SLC, ReaxFF (Fogarty et al.), and GFN0-XTB. Evaluation datasets and protocols: Independent test set (not used for training) of 1460 configurations from NNPscan simulations: (i) randomly sampled near-equilibrium structures from optimized zeolite databases; (ii) silica bilayer Stone–Wales defect pathway structures; (iii) high-energy snapshots from glass melting and ZA. Data available at Zenodo (doi:10.5281/zenodo.5827897). NNP extrapolation indicators included ensemble energy spread; deviations increased significantly once spread exceeded ~8 meV atom−1. Applications: Final NNPs used for reoptimization of Deem and IZA databases, for simulating glass melting and zeolite amorphization in FAU (not in training), and for predicting equilibrium structures and vibrational properties for comparison with DFT and experiment.

- DFT benchmarking: SCAN+D3 and PBE+D3 yield relative energies in very good agreement with experimental phase-transition enthalpies (MAD ~2–3 kJ mol−1). SCAN+D3 slightly outperforms PBE+D3 for structural metrics (e.g., density MAD 0.2 vs 0.3 Si nm−3). For the Stone–Wales defect in a silica bilayer, PBE+D3 barriers are slightly underestimated (up to −0.67 eV) and SCAN+D3 slightly overestimated (up to +0.59 eV) relative to B3LYP. - NNP accuracy vs other PES approximations (with respect to SCAN+D3): • eNNPscan (ensemble): EQ energies RMSE 3.95 meV atom−1 (MAE 2.83); EQ forces RMSE 0.067 eV Å−1 (MAE 0.046). All-data energies RMSE 5.49 meV atom−1 (MAE 3.83); forces RMSE 0.265 eV Å−1 (MAE 0.155). • sNNPscan (single): EQ energies RMSE 4.20 meV atom−1 (MAE 2.99); EQ forces RMSE 0.070 eV Å−1 (MAE 0.048). All-data energies RMSE 5.44 meV atom−1 (MAE 3.90); forces RMSE 0.303 eV Å−1 (MAE 0.175). • SLC: EQ energies RMSE 111 meV atom−1 (MAE 88.0); EQ forces RMSE 3.431 eV Å−1 (MAE 2.612). All-data energies RMSE 313 meV atom−1 (MAE 207); forces RMSE 4.166 eV Å−1 (MAE 3.206). • ReaxFF: EQ energies RMSE 78.4 meV atom−1 (MAE 56.4); EQ forces RMSE 2.996 eV Å−1 (MAE 1.266). All-data energies RMSE 136 meV atom−1 (MAE 88.8); forces RMSE 8.533 eV Å−1 (MAE 2.789). • GFN0-XTB: EQ energies RMSE 106 meV atom−1 (MAE 57.5); EQ forces RMSE 0.787 eV Å−1 (MAE 0.302). All-data energies RMSE 201 meV atom−1 (MAE 127); forces RMSE 3.391 eV Å−1 (MAE 0.735). These results demonstrate roughly an order-of-magnitude improvement in energy and force accuracy of NNPscan over analytical force fields and tight-binding DFT. - Database reoptimization and synthesizeability: Reoptimization of the Deem and IZA databases with NNPs provided high-accuracy energetics and revealed over 20,000 additional hypothetical zeolites within the thermodynamically accessible range for synthesis, offering richer data for downstream ML studies. - Reactive transformations: NNPs accurately reproduce high-temperature bond-breaking events (e.g., glass melting at ~4800 K) with transition-state motifs (fivefold Si) consistent with DFT. Zeolite amorphization simulations (LTA, SOD in training; FAU as an out-of-training test) showed bond-breaking in LTA near densities 2.1–2.2 g cm−3 and no cleavage in FAU below ~2.2 g cm−3, matching DFT and indicating strong transferability. - Extrapolation behavior: At densities above ~2.2 g cm−3 (up to ~10% beyond silica glass density), NNP energy errors rose (up to ~40 meV atom−1), with increased ensemble prediction spread (>~8 meV atom−1) signaling extrapolation. Even then, energies maintained qualitative agreement with SCAN+D3 (Pearson r ≈ 0.99 for FAU, 0.66 for LTA), and errors remained substantially smaller than other PES models (e.g., ReaxFF RMSE ~136 meV atom−1).

The study demonstrates that SchNet NNPs trained with active learning on a diverse silica dataset can interpolate the PES across low- and high-energy regions relevant to zeolite structures, equilibria, and reactive transformations, achieving DFT-level accuracy. The NNPs notably surpass analytical force fields and tight-binding DFT by about an order of magnitude in both energy and force accuracy, directly addressing the need for accurate, transferable, and reactive modeling tools in zeolite research. Robustness under high-temperature, bond-breaking scenarios (glass melting) and during structural collapse simulations indicates reliable handling of reactive events. The reoptimization of extensive zeolite databases and identification of many thermodynamically accessible hypothetical frameworks underscores the utility of these NNPs for guiding materials discovery and for studying structure–stability–property relationships at scales impractical for direct DFT. The controlled extrapolation detection with committee variance provides a principled pathway for continual refinement, ensuring the model remains reliable as it explores new regions of configuration space.

Reactive SchNet NNPs trained on a broad, actively learned silica dataset deliver DFT-accurate energies and forces across zeolite structures, enabling predictive modeling of equilibrium properties, thermodynamic stabilities, and reactive/non-reactive phase transitions. Compared with SLC, ReaxFF, and GFN0-XTB, the NNPs provide an order-of-magnitude improvement in accuracy. They facilitated high-fidelity reoptimization of the Deem and IZA databases and uncovered over 20,000 additional hypothetical zeolites within the accessible thermodynamic range, offering rich data for future ML studies of structure–stability–property relationships. Active learning and ensemble-based extrapolation detection ensure robustness and efficient dataset expansion. Future work includes extending the NNPs to heteroatom-containing zeolites (e.g., Al) and guest species (e.g., water) to enable realistic operando modeling under synthesis and operating conditions with ab initio accuracy.

- The NNPs inherit the accuracy limitations of the reference DFT methods; they cannot outperform SCAN+D3/PBE+D3. - Zeolite amorphization simulations mimic thermal collapse and required artificial compression; purely siliceous LTA and FAU have no reported thermal collapse experimentally. - Extrapolation occurs at densities above ~2.2 g cm−3 (up to ~10% beyond glass density), with increased errors, though still smaller than other PES approximations. - Training focused on pure silica; transferability to heteroatom-containing frameworks and guest-loaded systems requires further dataset extension and retraining. - Detailed publication metadata (e.g., exact title, full author list, affiliations, and publication date) are not provided in the text excerpt.