A reactive neural network framework for water-loaded acidic zeolites

Under operating conditions, the dynamics of water and ions confined within protonic aluminosilicate zeolite micropores are responsible for many of their properties, including hydrothermal stability, acidity and catalytic activity. However, due to high computational cost, operando studies of acidic zeolites are currently rare and limited to specific cases and simplified models. In this work, we have developed a reactive neural network potential (NNP) attempting to cover the entire class of acidic zeolites, including the full range of experimentally relevant water concentrations and Si/Al ratios. This NNP has the potential to dramatically improve sampling, retaining the (meta)GGA DFT level accuracy, with the capacity for discovery of new chemistry, such as collective defect formation mechanisms at the zeolite surface. Furthermore, we exemplify how the NNP can be used as a basis for further extensions/improvements which include data-efficient adoption of higher-level (hybrid) references via Δ-learning and the acceleration of rare event sampling via automatic construction of collective variables. These developments represent a significant step towards accurate simulations of realistic catalysts under operando conditions.


Introduction
Zeolites are a class of microporous aluminosilicates with tremendous structural and chemical diversity, which originates from the myriad stable three-dimensional arrangements of covalently connected silica/alumina tetrahedra.This makes zeolites a versatile material class with applications ranging from thermal energy storage to gas separation and water purification, but predominantly in heterogeneous catalysis.[1,2] The presence of aluminium, and in particular the necessary charge compensation add another layer of complexity to the structural characterisation of these materials, but are crucial to the catalytic function of zeolites.For example, acidic zeolites, in the form of protonated (H-form) aluminosilicates (H-AS) are one of the cornerstones of industrial petrochemical processes.[3] Recently, great experimental and theoretical efforts have been made to go beyond the traditional applications of zeolites, for example in converting sustainable bio-feedstocks into chemicals.[4,5,6] A further critical consideration for both existing and emerging applications is the interaction between H-AS and water.This relationship governs many features of H-AS including i) proton solvation, and thus acidity, [7,8] ii) hydrolytic bond dissociation and defect formation, which controls catalyst durability and activity, [9] iii) water mobility and clustering in zeolite pores [10] and iv) the synthesis of zeolites from precursor gels containing silica fragments, water and cations [11].Owing to the microporous nature of zeolites, this interaction is not adequately viewed as a simple bulk-liquid interface, but rather a collection of complex binding, cluster-ing, exchange and reaction steps between variously sized water clusters and an inhomogeneous surface that is complicated by topology-dependent confinement effects [12].As a result, the proper mechanistic understanding of H-AS water-zeolite interactions is still lacking, limited to either static calculations at ultra high vacuum conditions [13,14] or exploratory dynamical (AIMD) simulations of narrow scope.[15,16,17,18] These investigations demonstrate the importance of capturing dynamics under operating conditions, being able to discover unexpected reaction mechanisms and defective species that hitherto eluded structural identification, but are not sufficiently economical for a global exploration of structural and reactive space.
A state-of-the-art tool for accelerating the reactive sampling in zeolite-water systems, and thus reaching experimentally relevant timescale or realistic levels of model complexity is the class of reactive analytical potentials, for example ReaxFF.[19] However, due to their fixed functional form, these potentials have limited transferability to systems with different chemical composition.[20] Therefore, they frequently require re-parameterization for a specific system for fine-tuning.[21] An emerging alternative to analytical force fields is represented by machine learning potentials (MLPs), which interpolate the potential energy surface (PES) at the level of an ab initio training set.[22,23] Two paradigms dominate the MLP field currently: i) training of MLP that cover large parts of the chemical space with dozens of elements, but with limited coverage of the configuration space, e.g., not considering all relevant chemical reactions with the associated transition states, e.g., OC22.[24] and ii) active learning procedures to accelerate simulations for a specific system [25,26], which capture the details of the PES, including transition states, but have little or no transferability to systems with different chemical composition.Hence, MLPs that are able to simultaneously cover the broad chemical and configurational space needed for a class of materials such as H-AS zeolites, including the complexity of framework and water-framework based reactive transitions, are currently missing.
In this work, we developed reactive global neural network potentials (NNP) for an entire material class, namely, H-AS zeolites.These potentials capture the chemical space from dense silica and alumina polymorphs, through water-containing H-AS zeolites of all experimentally relevant Si/Al ratios, to bulk water and water gas-phase clusters.We observed excellent transferability among unseen framework topologies, with consistently high accuracy with respect to the training reference data.Generalization tests showed hitherto unseen chemical species and processes, including a collective hydrolysis mechanism at the surface of a zeolite nanosheet.
Finally, we show that the learned representations of the NNP baseline models can be used for data-efficient learning of higher (hybrid) DFT level corrections (∆-learning) [27] for specific use cases, in addition to developing machine learned collective variables for the acceleration of rare event sampling.[28]

Database generation and training of the general H-AS NNPs
One of the challenges in training general NNPs for a material class is creating an interpolation grid that captures relevant parts of the configuration and chemical space.The computational procedure employed in this work is summarized in Figure 1a.The bulk of the database is derived from 500 short (10 ps) ab initio ab initio molecular dynamics (AIMD) trajectories (at PBE+D3(BJ) level) [29,30,31] using a set of H-AS zeolite models.This structure set contains 150 zeolites constructed using ten topologies (three existing and seven hypothetical) with varying Si/Al ratios and water loading at three temperatures ranging from 1200 K to 3600 K (see Supplementary Table 1) to sample both low-and high-energy parts of the potential energy surface (PES).In addition, nine AIMD runs were performed for bulk water at three densities (0.9-1.1 g cm −3 ) and at three temperatures (300-900 K).All AIMD trajectories were then subsampled by Farthest Point Sampling (FPS) [32,33] using metric based on the smooth overlap of atomic positions (SOAP) [34] kernel to obtain a collection of de-correlated and structurallydistinct configurations.These structures were used for SCAN+D3(BJ) [35] single-point (SP) calculations to create the bulk of the training database.To systematically sample states close to the equilibrium structures, lattice deformations (see Supplementary methods) were applied to the optimized structures of the H-AS zeolite models mentioned above and these structures were then used for SCAN+D3(BJ) single point (SP) calculations.Finally, to further diversify the database, the same lattice deformations were used for alumina and ice polymorphs as well as water clusters [36] (see Methods section).The resulting database was used to train (an ensemble of six) SchNet-based NNPs achieving average test root mean square errors (RMSE) of 5.3 meV atom −1 and 186 meV Å−1 for energies and forces, respectively (see Supplementary Table 2).These errors are similar to other reactive (rotationally invariant) MLPs [37,38,39].Figure 1b provides a low-dimensional representation of the training database.It shows a t-distributed stochastic neighbor embedding (t-SNE) [40] plot of the averaged SchNet representation vectors (see Supplementary methods) to visualize the structural and chemical diversity of the training database, ranging from water-free alumina and silica systems through various water-loaded H-AS zeolites to bulk water and small clusters.The t-SNE components of the averaged SchNet representations change smoothly with the chemical composition of the structures as well as with their total energy (see Supplementary Figure 1).In addition, all generalization tests (see below) lie within the generated interpolation grid of the database and cover a wide range of application cases for zeolite modeling.Generalization tests are highlighted in red.c Energy error distribution ∆E r (see Eq 1) of the NNPs in comparison with ReaxFF [41].

Generalization tests and exploration of configuration space
To properly test the generalization abilities of the trained NNPs, we employed a series of simulations for systems outside of the training domain, i.e., including: i) MD simulations at ambient conditions that probe the performance of NNPs for close-to-equilibrium structures and ii) high-temperature MD simulations (supplemented by nudged elastic band (NEB) transition path searches) to assess the NNP quality for highly activated, reactive events.The systems considered sample the chemical and structural space of water-loaded H-AS zeolites (see Figure 1b) varying in water and aluminum content, as well as in the zeolite topology (FAU, GIS and MFI zeolite frameworks which are not seen during the training) and dimensionality of the H-AS systems (three-dimensional crystal, zeolite layer or a zeolitic molecular fragment interacting with bulk water).Further details about the performed MD test simulations including discussion of some chemically relevant observations are provided below in the "Sampling equilibrium properties" and "NNP robustness at high temperatures" sections below.Here, we focus on overall performance of our trained NNP in these generalization tests.To compare energies of all computational methods across the chemical H-AS space, we used the energies E r of the hypothetical formation reaction: with α-quartz, corundum (α-Al 2 O 3 ), and a single water molecule in the gas phase as reference structures.Adoption of E r allows for benchmarking methods across broader chemical space, as exemplified by Hautier et al. [42].Alternatively, we also expressed the errors of relative energies ∆E with respect to a reference configuration for each model system, e.g., the initial structure of an MD trajectory (see Supplementary Table 3).Obviously, the force errors as intensive properties are independent of the reference and can be compared directly across the chemical space Table 1 summarizes the RMSEs of energies and forces for all test cases (2700 structures in total) of the NNP with respect to SCAN+D3(BJ) reference and Figure 1c shows the total energy error ∆E r distribution (see Supplementary Figure 2 for ∆E r and ∆E distributions for each system separately).Table 1 and Figure 1c also show the performance of the reactive analytical force field ReaxFF specialized for water-loaded H-AS systems.[41] Table 1.Root mean square errors of reaction energies ∆E r (see Eq 1) and forces for the test cases shown in Fig. 2  The total NNP errors are similar to other state-of-the-art (rotationally invariant) MLPs [37,38,39] for the modeling of reactive events and outperform ReaxFF by more than an order of magnitude in both energy and forces.More importantly, the NNP calculated reaction energies E r are consistent over the entire range of chemical H-AS compositions and configurations.
Only in the case of GIS(T =3000K)+24H 2 O, energy and force errors are about twice as high compared to the other test cases.Such higher errors were also obtained for MLPs when applied to simulations with a large number of reactive events at extreme temperatures.[38,37] To put these NNP errors in context, note that standard GGA-level DFT functionals show, for 27 formation reactions that involve silica, an RMSE of about 28 meV atom −1 with respect to experiment.[42] Therefore, the NNPs safely retain the (meta)GGA-level DFT quality for the description of the water-loaded H-AS systems.By contrast, ReaxFF shows, e.g., a relatively low energy RMSE for the high silica FAU tests (FAU(Si/Al=47)+nH 2 O) but ten times higher energy errors for FAU(Si/Al=1).This difference arises from shifted ∆E r energy error distributions for ReaxFF and, to a far lesser extent, for the NNPs (see Supplementary Figure 2).Such almost constant shifts are removed when taking a structure with the same chemical composition as the reference, that is, comparing relative energies ∆E with DFT.As an example, the relative energy RMSE ∆∆E of GIS(T =3000K)+24H 2 O (NNP: 6 meV atom −1 ; ReaxFF: 83 meV atom −1 ) are about two times lower than ∆E r (see Supplementary Table 3).In addition, the NNPs provide an order of magnitude higher force quality than RexFF across all generalization test cases.This comparison demonstrates the ability of the NNPs for general modeling of structure, properties and chemical reactivity across the class of H-AS materials.

Sampling of equilibrium properties
The dynamic behaviour of zeolite-confined water containing solvated protons is of high interest due to (potential) application of zeolites in water purification [10,43], heat storage [44] or reaction optimization under humid conditions (such as biomass conversion).[7] The ability to realistically model the (water-loaded) H-AS systems close to equilibrium is crucial for understanding of many of their signature properties such as acidity, (water) adsorption and diffusion, or relative stability as a function of topology, water content and aluminum distribution and concentration.The role of water loading and aluminum concentration was probed using equilibrium MD simulations of water-loaded zeolite with FAU (faujasite) topology -an industrially important zeolite topology unseen in the NNP training -under standard conditions (300 K).
We considered model systems with the (theoretically) lowest and highest possible Si/Al ratios in a (primitive) FAU unit cell, namely, a single Brønsted acid site (BAS) with Si/Al=47 and Si/Al=1 according to Löwenstein's rule that prohibits the formation of Al-O-Al pairs.In the case of Si/Al=47 (FAU(Si/Al=47)+nH 2 O), three water loadings n were tested, from single water through a water tetramer to full water loading of FAU with 48 water molecules (approximate water density of 1 g cm −3 ).For Si/Al=1, full water loading with 48 molecules per FAU unit cell was chosen to focus on extensive sampling of BAS protonation and deprotonation events, a key reactive event characterizing these strong solid acids.
In the case of FAU (Si/Al=47) model, the single water molecule (n=1) remains adsorbed at the BAS throughout the 1 ns MD simulation, in line with the very strong interaction between BAS and water molecule amounting to -79 kJ mol −1 calculated here and with adsorption energies reported in the literature (-70 to -80 kJ mol −1 in CHA [45]).We quantified the degree of solvation by calculating the minimum distance of Al-O FW -Si framework oxygens to all hydrogen atoms (see Supplementary Figure 3).The proton was considered solvated if it is closer to a water oxygen than to Al-O FW -Si.Only very few solvated states (less than 3% of the trajectory) were observed during the 1 ns run, in line with previous (shorter) AIMD simulations.[45] However, the water tetramer (n=4) is already able to deprotonate the BAS, but similarly to single water, the tetramer stays close to the framework Al during the 1 ns MD trajectory (on average 3.1 Å between Al and the cluster center-of-mass, see Supplementary Figure 3).At full water loading (n=48), the proton rapidly leaves the BAS and stays solvated with an average distance of 7.3 Å (ranging from 3-10 Å) from the Al-O FW -Si (see Supplementary Figure 3).However, the evaluation of the confined water dynamics herein is hindered by finite-size effects, due to a small primitive cell of FAU chosen primarily for benchmarking purposes.Hence, we refer the interested reader to our preliminary work [46] using the NNPs presented here on water diffusion in FAU using appropriately sized FAU unit cell (cubic cell with edge length of 25 Å) that is prohibitively large for carrying out routine DFT calculations.
A particularly challenging case is FAU with Si/Al=1.The MD trajectory contains several protonation and deprotonation events that cannot be described by analytical, non-reactive force fields.Also the ReaxFF parameterization used in this work shows considerably higher errors compared to FAU(Si/Al=47) (see Table 1 and Supplementary Figure 2).On the other hand, the errors of the NNP for this challenging case increase only mildly and are well below the test errors for the training database (see above).
To check the NNP generality further, we employed the "inverse" models to water-loaded zeolites, namely the fragments of zeolites (silicic acid Si(OH) 4 and aluminium hydroxide Al(OH) 3 ) solvated in bulk-like water (using a simulation box containing 96 water molecules).Such systems are relevant for modeling potential precursors of zeolite synthesis or products of zeolite degradation in hot liquid water (de-silication and de-alumination).[5,9,47] Since both processes take place under rather harsh hydrothermal conditions (temperature above 100 • C and pressure above 10 bars), the test simulations were carried out at 500 K.In both cases, the NNPs accurately reproduce the SCAN+D3(BJ) energies and forces.Similar to the previously discussed system FAU(Si/Al=1)+48H 2 O, the NNP energy errors (4 meV atom −1 ) are mainly connected to the offset of E r (see Supplementary Figure 2).These results suggest that NNPs enable us to run reliable large-scale equilibrium MD simulations across the chemical and configuration space of H-AS and water with the (meta)GGA level of accuracy.

NNP performance for highly activated reactive events
Modeling of chemical reactions at the H-AS-water interface needs a robust interpolation of the relevant transition states.However, MLPs are expected to have only limited capability to reliably describe configurations and energetics in extrapolated or sparsely interpolated regions of the potential energy surface, [22] which often coincide with the high energy transition state configurations.Therefore, we tested the trained NNPs by performing MD simulations at very high temperatures (at 1600 and 3000 K) for an unbiased assessment of the NNP quality and robustness for modelling reactive processes.We chose two systems that were not part of the training dataset: an interface model of siliceous MFI slab in interaction with the bulk-like water and water-loaded GIS with Si/Al=1.In addition, we tested the accuracy of the NNPs using static nudge-elastic band (NEB) calculations, that are used to locate specific transition pathways, for four reactions relevant for H-AS zeolite based catalysts (see below).The first generalization test was a model for the external interface of siliceous MFI with bulk water.This zeolite model resembles 2D MFI nanosheets which were successfully prepared by exfoliation of a multilamellar MFI zeolite.[48] To sample reactive events at the external zeolite surface, we performed an exploratory MD run at 1600 K for 1 ns.No extrapolation was detected using the ensemble of NNPs.As expected, significantly increased temperature leads to increased probability of the highly activated reactive events to take place and we do observe a relevant chemical reaction taking place over the course of 5 ps, in which a silanol defect is created at the external MFI surface (see Figure 2a).In our previous work, [15] a similar reactive event was found to be characterized by a free energy of activation amounting to approx.80 kJ mol −1 (at 450 K).The observed reaction starts at the intersection of the bulk water phase with the MFI main channel (along the crystallographic b-direction).Firstly, a water molecule adsorbs at a surface Si site (Figure 2b).The autoprotolysis of a nearby water molecule leads to the transfer of a proton from the adsorbed water molecule to the formed hydroxide ion creating an additional surface silanol group at the five-fold coordinated Si atom (Figure 2c).The remaining hydronium ion shuttles the excess proton together with the surrounding water molecules to a framework oxygen bound to the five-fold coordinated Si (Figure 2d).This process finally leads to the cleavage of the Si-O(H) bond, creating a silanol defect in axial position to the previously formed surface silanol.Our exploratory MD simulation revealed a feasible reaction mechanism for silanol defect creation at the external zeolite surface involving the autoprotolysis of water, which would be challenging to find by biased dynamics simulations with human designed CVs.These findings are in line with previous experimental studies on the hydrolysis of MFI zeolites in hot liquid water that suggest that zeolite degradation predominantly starts at the external surface in which water autoprotolysis and silanol groups play a crucial role.[49,50] To confirm that the defect creation process observed at the NNP level is reliable, we performed SCAN+D3(BJ) SP calculations using 100 snapshots comprising the reaction steps depicted in Figures 3b-3e.The NNPs proved very accurate for this test case with an energy RMSE of 1.4 meV atom −1 (see Table 1), in contrast to the employed ReaxFF parameterization with a rather broad error distribution (see Figure 2f) and an RMSE of 38 meV atom −1 .
The second particularly challenging generalization test was a GIS zeolite model (Si/Al=1) loaded with waters (24 molecules), which was molten at 3000 K for 2 ns to sample multiple highly activated reactive events taking place simultaneously (see also Supplementary Figure 4).
We obtained a stable MD trajectory of the liquid H-AS state with thousands of bond-breaking events (e.g., Si-O and Al-O bond hydrolysis, aluminol and silanol formation) over the entire simulation time, without detecting extrapolation using the trained ensemble of NNPs.However, this test case shows higher energy RMSE by around a factor of two when compared to the other test cases (see Table 1).Similar trends of increased energy errors have also been observed for other MLPs applied to high temperature MD runs of the liquid state of strongly (covalently) bound materials (see e.g.Refs [37,38]).Even though the NNP accuracy mildly deteriorates at these extremely high temperatures, they proved robust in these simulations of a large variety of highly activated reactive events.
Lastly, we tested the accuracy of trained NNPs on specific elemental reactions in water-loaded H-AS zeolites with well-known transition states: a proton jump with and without water, [51,52,53] and water-assisted bond breaking mechanisms of the Si-O and Al-O(H) bonds [14,54,55] (see Supplementary Figure 5 and Supplementary Table 4 as well as Figure 3).All NEB calculations were performed for FAU; an industrially relevant zeolite topology that was not part of the training dataset.For quantification of the NNP error, we used SCAN+D3(BJ) SP calculations on all generated NEB images.On average, the relative NNP energies only slightly deviate from their DFT reference with an RMSE of about 6 kJ mol −1 .Such small errors for activation barriers can be considered to lie within DFT accuracy.[56] The generalization tests of the herein trained NNPs presented above demonstrate that the trained NNPs are robust and general interpolators for simulations across the chemical and configuration space spanned by the water-loaded H-AS zeolites and that they are able to retain the reference-level (SCAN+D3(BJ)) DFT quality not only for close-to-equilibrium simulations but also for highly activated reactive events.

Extensions of the NNP model
Obtaining a general NNP model that describes water-loaded H-AS zeolitic systems with (meta)GGA DFT quality with several orders of magnitude speedup is clearly beneficial.However, with such a robust baseline NNP model available, it is possible to construct extensions that can improve either the accuracy of the description or the efficiency of sampling of the reactive events of interest.

Improving the baseline model using ∆-learning
To improve the accuracy of the benchmark level, one can employ the well-known ∆-learning concept.[27] In this way, one can train a correction model on top of the baseline model, using a computationally more demanding but more accurate level of theory for a small subset of datapoints that covers the system of interest.We demonstrate the applicability of the concept for two model reactions discussed above -a proton jump and water-assisted scission of Al-O(H)-Si bonds.[14] (see Figures 4c-d) For the higher-level reference, we chose the range-separated hybrid DFT functional ωB97X complemented with the empirical dispersion correction D3(BJ): a functional that shows considerably better performance for water cluster binding energies and reaction barriers [57,56,58] than our baseline reference SCAN+D3(BJ) functional.First, we generated a small ωB97X-D3(BJ) database containing 500 structures taken from the biased (NNP level) MD runs of a H-jump in water-free CHA (between O2 and O3, see Methods section) taken from Ref. [28].Next, we trained a correction (∆NNP) to the atomic energies of the NNP baseline model by using a simple linear regression (see Method section).It turned out that 150 training structures were sufficient to reach (test set) RMSEs of 1.3 meV atom −1 and 69 meV Å−1 for energy and forces, respectively (see Supplementary Figure 6).We then tested whether the trained ∆NNP model is capable to model a proton jump and Al-O(H) bond dissociation mechanism in a zeolite topology different from those included in the training set.Again, we chose FAU with a single Al defect (Si/Al=47) as a test case.Figure 3a and 3b show the results of the NNP level NEB calculations along with the corresponding DFT energies.Figure 3c and 3d depict the structures for both reaction paths and Table 2 also compares the relative energies of the proton jump in CHA (O2-O3) and FAU (O1-O4).Not surprisingly, both the general baseline and ∆NNP model are in very good agreement with their DFT reference in case of the proton jump in CHA with less than 3 kJ/mol deviation.The ∆NNP accuracy of the relative energies only slightly deteriorates (about 5 kJ mol −1 error) when applied to FAU.When comparing the reaction energies ∆E r (see Eq 1), the ∆NNP model shows an almost constant shift to its DFT reference (see Supplementary Table 5).This offset shows that the ∆NNP model is capable of describing local atomic environments (e.g., of BAS) but with less transferability to different compositions and zeolite frameworks as the baseline NNP model.red, H: white).e, f Estimated free energy profiles using ML collective variables.
To check how the ∆NNP model performs for different but related reactions to proton jumps, we repeated the NEB calculations at the ∆NNP level for the Al-O(H) bond dissociation (Al-O2 bond cleavage, see Method section) mechanism suggested by Silaghi et al. [14] (Figure 3b and 3d).The ∆NNP model reproduces the ωB97X-D3(BJ) relative energies ∆E of the reaction path with near-DFT accuracy albeit with somewhat higher energy errors (less than 8 kJ mol −1 ) compared to the CHA proton jump and the NNP baseline model (see Table 3).To calculate ∆NNP energy and force errors, we performed single-point calculations for the biased dynamics runs of the reactions (in FAU) shown in Figure 3 (see Method section).Similar to the NEB calculations, low RMSEs for the relative energies ∆E of less than 2 meV atom −1 were obtained (see Supplementary Table 6).However, the reaction energies ∆E r (see Eq 1) show an almost constant offset (see Supplementary Figure 7), as observed for the baseline NNP and ReaxFF (see Supplementary Figure 2), leading to larger RMSEs of up to 18 meV atom −1 .
This offset shows that the ∆NNP model does not retain the generality of the NNP baseline model across different zeolite topologies.However, the low ∆E errors indicate that the ∆NNP model correctly interpolates the local environments of the BAS in a low water (or water-free) regime.
In addition, the ∆NNP model shows low force errors of around 100 meV Å−1 for both reactions (see Supplementary Table 6) which is comparable to the accuracy of the NNP baseline model (see Table 1).This finding supports the conclusion that the ∆NNP model is capable to model the local atomic environments in FAU even when trained on only 150 structures of a different zeolite (CHA).For comparison, Bocus et al. [53] trained specialized SchNet NNPs on biased dynamics trajectories of proton jumps in CHA using more than 100k training points.
When applying these NNPs to different zeolite topologies, they obtained up to twice as high force errors (around 200-250 meV Å−1 ).These results indicate that the learned representation vectors of the more general NNP baseline model contain information on the H-AS PES, allowing a data-efficient ∆-learning of higher level corrections.Therefore, the ∆NNP exhibits fairly good extrapolation robustness even when applied to similar reaction pathways in other zeolite topologies.Accelerating rare event sampling using baseline model representations In the previous section, we tested the ∆-learned model for accurate (hybrid DFT) modeling using static calculations of known reaction mechanisms.However, prior investigations have shown the unforeseen and highly collective nature of water-involved reaction mechanisms as well as the sizable role of temperature effects.[15,16,18] Both imply the need for a tool with the ability to effectively discover and sample transition pathways.For effective sampling of the activated reactive events, which are therefore rare on the timescales accessible even for the NNP-accelerated simulations, one typically adopts a biasing along a low-dimensional representation of the reactive process, i.e., along the reaction coordinates or collective variables (CVs).
However, good CVs can be difficult to construct in case of unknown, possibly complex, reaction pathways.
Our recent work [28] shows how the end-to-end learned atomic representations of our baseline NNP model can be used to automatically generate robust machine-learned CVs (ML-CVs).In this approach, the structures of the reactant, product, and perhaps also tentative transition states are first represented using the atomic representations of the herein-trained baseline NNP model.
Next, these representations are used as an input for the dimensionality reduction model (variational autoencoder), which generates a low-dimensional (typically one-or two-dimensional) latent space from which the model attempts to reconstruct the input representation vectors as precisely as possible.As a result, the latent low-dimensional space effectively distinguishes products from reactants, i.e., it represents the reactive coordinate, or collective variable.We showed previously that learned ML-CVs coupled with the baseline NNP enable efficient sampling of the free energy surface for a proton jump and Si-O bond hydrolysis in CHA zeolite.[28] Here, we test this procedure using the aforementioned proton jump and Al-O2(H) bond dissociation mechanism in FAU with Si/Al=47 which is outside of the NNP training domain, in contrast to CHA (see Methods section).Figure 3e and 3f show the estimated free energy profiles calculated with ML-CVs using well-tempered metadynamics [59] simulations (see Supplementary Information and Methods).The free energy barrier (approx.110±10 kJ mol −1 at 300 K) of the proton jump in FAU is somewhat higher compared to the static calculations (84 kJ mol −1 , see Table 2).This is in line with previous calculations [52,53] which showed increasing reaction barriers with temperature (up to 20 kJ mol −1 from 0 K to room temperature).In case of the Al-O(H) bond dissociation the activation free energy is about 80 kJ mol −1 , similar to the barrier found by the NEB simulations (see Table 3).Hence, with the baseline NNP model, one can not only accelerate the evaluation of energies (and forces) necessary for sampling the water-loaded H-AS systems but also use it to automatically generate ML-CVs accelerating the sampling of a particular reactive process.

Conclusions
In this work, we developed a neural network potential (NNP) for the entire class of protonic aluminosilicate (H-AS) zeolites, which are one of the cornerstones of existing petrochemical processes, [3] as well as one of the main candidates for emerging applications of zeolites in sustainable chemistry.[2] Our NNP provides a general approximation of the potential energy surface of the H-AS zeolites, including reactive interactions with water, capturing both closeto-equilibrium structures and high-energy bond-breaking scenarios.The ability to cover large portions of both the chemical and configurational space of this material class was demonstrated using multiple generalization tests that ranged from zeolite surfaces varying in water and aluminum content to zeolite fragments solvated in bulk-like water and a high-temperature melt of the aluminosilicate zeolite GIS.These tests confirm the outstanding transferability of the NNP, which are able to maintain consistent accuracy close to the reference meta-GGA DFT level, across the entire H-AS material class, outperforming state-of-the-art analytical reactive force fields for water-loaded H-AS zeolites [41] by at least one order of magnitude.Moreover, in some of these tests we observed hitherto unseen chemical processes and species, which confirms the capability of the NNP for exploration and discovery of novel reactive pathways, in addition to acceleration of configuration space sampling.
Furthermore, we showed that the NNP can be used as a robust baseline model that allows for extensions including: i) data-efficient adoption of higher-level (range-separated hybrid DFT) description via ∆-learning [27] and ii) acceleration of reactive event sampling using automatic construction of collective variables, via end-to-end learned atomic representations [28].Hence, the baseline model with its extensions constitutes a broader ML-based framework within which one can simulate H-AS materials in a comprehensive bias-free fashion with tunable accuracy.
We expect that the tools developed in this work will enable large-scale simulations of H-AS zeolites to tackle long-lasting challenges in the field, ranging from understanding the mechanistic underpinnings of zeolite hydrothermal (in)stability to the determination of the character of active species and defects under operating conditions.

Dataset generation
Covering the chemical and configuration space of H-AS zeolites requires a structurally distinct set of zeolite frameworks with different water loadings and Si/Al ratios.In our previous publication, [37] we used Farthest Point Sampling [32,33] together with the smooth overlap of atomic positions [34] descriptor (SOAP-FPS) to find a subset of siliceous zeolites that optimally covers the structural diversity of existing and more than 300k hypothetical zeolites.From this subset, ten zeolites were selected, three existing (CHA, SOD, MVY) and seven hypothetical zeolite frameworks (see Supplementary methods).These frameworks were used for construction of 150 initial structures combining four water loadings (from 0 to ∼1.1 g cm −3 ) with three Si/Al ratios between ∼1-32 (in protonic form) and water-loaded purely siliceous zeolites (see Supplementary Table 1).We also added a two-dimensional silica bilayer (12 Å vacuum layer) used in ref. [37] with three different water-loadings to the initial structure set.All 153 initial configurations were then optimized under zero pressure conditions.
Next, the entire structure set was equilibrated for 10 ps at 1200, 2400, and 3600 K using ab initio MD (AIMD) simulations to sample reactive events at higher energies.Sampling of the low energy parts of the potential energy surface (PES) used 210 unit cell deformations applied to all optimized structures (see Supplementary methods).Apart from microporous structures and two-dimensional H-AS, we also added the same set of 210 lattice deformations for six dense H-AS polymorphs, namely, four alumina polymorphs α-Al 2 O 3 (corundum), [60] θ-Al 2 O 3 , [61] γ-AlO(OH) (Boehmite), [62] and α-Al(OH) 3 (Gibbsite), [63] as well as two aluminosilicate polymorphs, Si 3 Al 2 O 12 H 3 (H 3 O-Natrolite) [64] and Al 2 Si 2 O 5 (OH) 4 (Dickite).[65] Additionally, we (SOAP-FPS) subsampled AIMD trajectories of zeolite CHA taken from previous publications [15,66] to further extend the structure database.These trajectories are equilibrium MD runs of non-Löwenstein pairs (Al-O-Al) with various water loadings (0, 1, 15 water molecules) [66] and biased AIMD runs of Si-O(H) and Al-O(H) bond cleavage mechanisms.[15] For interpolation of the interactions in pure water, we performed AIMD simulations (10 ps) for bulk water with 64 water molecules at three densities (0.9, 1.0, 1.1 g cm −3 ) and at three temperatures (300, 600, 900 K).In addition, we used single water and water clusters in vacuo taken from the BEGDB database [36]  Finally, the unit cells of two ice polymorphs (Ice II [68] and Ice I h [69]) were deformed in the same way for sampling of low-energy structures of crystalline water.
All AIMD simulations and structure optimizations were performed at the computationally less demanding PBE+D3(BJ) [29] level employing the dispersion correction of Grimme et al. (D3) [30] along with Becke-Johnson (BJ) [31] damping.The AIMD equilibration used the canonical (NVT) ensemble along with a 1 fs time step, the Nosé-Hoover thermostat [70,71] and with hydrogen being replaced by tritium.Structurally diverse configurations were extracted from the MD trajectories using SOAP-FPS (see Supplementary methods).These decorrelated MD structures were used, together with the generated set of lattice deformations, for single-point (SP) calculations at the (metaGGA) SCAN-D3(BJ) level.[35] The resulting SCAN+D3(BJ) reference dataset contained 248 439 structures.
An ensemble of six SchNet [72] NNPs was trained on the final SCAN+D3(BJ) database.Finally, two hundred structures were extracted from the MD run for the NNP and ReaxFF error evaluation.

Generalization tests and reaction path searches
All MD simulations used a time step of 0.5 fs with hydrogen being replaced by deuterium, employing the Nosé-Hoover thermostat.[70,71] The final generalization test set collected from all trajectories contains 2700 configurations for SP calculations at the SCAN+D3(BJ) and ReaxFF [41] level allowing the energy and force error evaluation shown in Table 1 and Figure 1c.
Lastly, we conducted NNP performance tests for modeling of reaction pathways in FAU (primitive unit cell, Si/Al=47) using NNP level climbing image NEB calculations.We chose four reaction pathways: i) a proton transfer with one water molecule and without water (between O1-O4), [51,52,53] and ii) water-assisted bond breaking mechanisms of the Si-O2 and Al-O2(H) bonds [14,54,55] (see Figure 3 and Supplementary Figure 5).In addition, we tested a water-free proton transfer in CHA between O2 and O3.The numbering of the symmetry inequivalent oxygen atoms (see Supplementary Figure 8) is consistent with the labeling of the zeolite frameworks in the IZA database (available under: iza-structure.org/databases).Energies at the SCAN+D3(BJ) were then obtained by SP calculations for all NEB images.

∆-learning
We applied the ∆-learning approach [27] to improve the accuracy of our baseline SCAN+D3(BJ) model to the (hybrid DFT) ωB97X-D3(BJ) level.First, we generated a training set using a subset of an (NNP level) biased dynamics run of a proton jump in CHA between O2 and O3, taken from Ref. [28].These structures were selected by FPS using the euclidean distance of the (baseline) SchNet NNP representation vectors averaged over each MD snapshot.SP calculations were then applied to 500 extracted configurations to obtain energies and forces at the ωB97X-D3(BJ) level.
The ∆NNP correction of the atomic energies ∆E i to the NNP baseline model (SCAN+D3(BJ)) was obtained by linear regression of the SchNet representation vectors x i of each atom i with the (column) weight vector w i and bias b i : ∆E i = w T x i + b i .We tested the ∆NNP model performance using 250 randomly chosen structures from the dataset and convergence tests showed that 150 training points give sufficiently low test set errors (see Supplementary Figure 6).To test the ∆NNP quality, we repeated the NEB calculations described above for the water-free proton jump and Al-O2(H) bond hydrolysis in FAU as well as the proton transfer in CHA (O2-O3) without water.In addition, we performed two hundred SP calculations for structures taken from the biased dynamics runs of the proton jump (O1-O4) and the water-assisted Al-O2(H) bond cleavage to improve the force error statistics of the ∆NNP model.

Biased dynamics using ML collective variables
Collective variables guiding proton transfer (O1-O4) and Al-O2(H) bond dissociation reactions in FAU described in the manuscript were trained using a variational autoencoder build on top of NNP-generated representation vectors.
To train the proton jump reaction, we used 3500 data points simulated from equilibrium MD runs on reactants and the same number on products.The data generated by using intuitively chosen CV and steered dynamics method were available for verification (see Supplementary Figure 9) but were not used during training.The NNP representations were calculated and saved in a cache and an encoder producing a single CV was trained together with the decoder.
We trained on 40 epochs with a learning rate of 10 −4 .The encoder generating the CV was a simple linear layer on top of pre-trained representations.To reduce the complexity of the task, only the 30 representation elements were chosen that maximized the variance between reactant and product configurations.For the biased dynamics, we used well-tempered metadynamics from the PLUMED package.[74,59] The parameters for the simulation are in Supplementary Table 7.The simulation was run for 1 800 000 steps using a 0.5 fs timestep.
Al-O(H) bond dissociation was trained similarly to the proton jump case.We used the same number of data points obtained by running unbiased trajectories in both end states.The encoder was again only linear and was trained for 60 epochs with a learning rate 2.10 −4 .200 representation elements were pre-selected from the representation vectors generated by the NNP.
Parameters for the well-tempered metadynamics method are reported in Supplementary Table 7.We ran the biased dynamics for both test reaction for 1 500 000 timesteps of 0.5 fs.To avoid simulating degrees of freedom that are unimportant for the Al-O(H) bond dissociation reaction we introduced some restrains to the biased dynamics.We required the distance between the free water molecule and the aluminium atom to be at most 2.2 Å to avoid it diffusing away.
We also fixed two hydrogen atoms to their corresponding oxygens to avoid permutations that would complicate the process.In the same fashion we disallowed the formation of the hydrogen bond (Figure 3d-FS) for different O and H pairs.
Otherwise, the setup is standard for the representation-based collective variables and a more thorough description of the autoencoder architecture and workflow can be found in [28].

Fig. 1
Fig. 1 Training and testing of general H-AS NNPs. a Computational workflow for creation of the SCAN+D3(BJ) database and application of the trained NNPs to test their generality.The end-to-end learned representations are used for ∆-learning and construction of ML collective variables.b t-distributed stochastic neighbor embedding (t-SNE) plot of the average representation vectors of all configurations in the training database (color codes shown on the left).

Fig. 2
Fig. 2 Surface defect creation in an 2D-MFI nanosheet.a Snapshot of 2D-MFI from an exploratory 1 ns MD run at 1600 K to sample reactive events (Si: yellow, O: red, H: white).b-e Reaction steps of silanol defect creation at the 2D-MFI-water interface: b water adsorption on a surface Si and water autoprotolysis; c proton transfer from the adsorbed water to the hydroxide ion; d migration of the hydronium ion and adsorption on a framework oxygen; e Si-O bond hydrolysis, creating a silanol defect in axial position to the formed surface silanol.f NNP and ReaxFF energy error distribution of 100 snapshots comprising the reaction steps b-e.

Fig. 3
Fig. 3 Reaction path modeling using ∆-learning and ML collective variables.Reaction paths of a proton jump (a, c, e) and an Al-O(H) bond dissociation (b, d, f) in FAU.a, b Static NNP simulations and corresponding DFT energies for the NNP baseline level (SCAN+D3(BJ)) and the ∆NNP level (ωB97X-D3(BJ)) trained on 150 structures (taken from Ref [28]) of a proton jump in CHA.c, d Atomic structures along the reaction path (Si: yellow, Al: grey, O: The six independent training runs used different, randomly split parts of the DFT dataset with approximately 80% of the datapoints as training set and 10% as validation and test set, respectively.We used the same SchNet hyperparameters (6 Å cutoff, 6 interaction blocks, 128 feature vector elements, 60 Gaussians for distance expansion) and loss function for training of energies and forces (trade-off 0.01) as in our previous publication.[37]Minimization of the loss function used mini-batch gradient descent along with the ADAM optimizer[73] and four structures per batch.If the loss function for the validation set did not decrease in three subsequent epochs, the learning rate was lowered (from 10 −4 to 3 • 10 −6 ) by factor 0.75.
Testing of the NNP accuracy, robustness and generality used a series of MD and NEB calculations of systems that were not included in the training database.To test the NNP performance at close-to-equilibrium (low temperature) conditions, we performed four MD runs (1 ns, 300 K) for zeolite FAU (primitive unit cell, 48 T-sites) at different chemical compositions.Three of the MD runs used a single Al atom (and BAS) per unit cell (Si/Al=47) and three water loadings (1, 4, 48 water molecules).The fourth run was performed with 24 Al per unit cell (Si/Al=1) and 48 water molecules.From every MD trajectory, 500 were selected for subsequent SCAN+D3(BJ) and ReaxFF SP calculations.As an "inverse" test case to three-dimensional zeolites, we chose silicic acid Si(OH) 4 and aluminum hydroxide Al(OH) 3 solvated in bulk water (96 water molecules).Both systems were equilibrated at hydrothermal conditions 500 K for 1 ns.Subsequently, two hundred configurations were selected from both MD runs for accuracy evaluation.To check the NNP performance and robustness for the sampling of reactive events, we first constructed a model of the external MFI-water interface.The starting point was an orthorhombic (96 T-site) MFI unit cell with one silanol nest at T-site T9 for exploratory, high-temperature MD runs to sample reactive events at the internal and external MFI-water interface.Next, an MFI(010) surface model was created by adding a 12 Å vacuum layer and cleaving the Si-O bonds between the T-sites T7, T9, T10, and T12 (lattice plane with lowest number of bridging O) yielding eight silanol groups on both surfaces.After addition of 165 water molecules, the model was equilibrated for 2 ns at 1600 K.One hundred structures were selected from the trajectory that include the surface defect creation shown in Figure2for SP calculations.As an extreme case to test the NNP robustness at very high temperatures, we simulated the liquid state of an H-AS-water model system at 3000 K.The initial configuration was a model of GIS (32 T-site unit cell) with Si/Al = 1 and 24 water molecules which was equilibrated for 2 ns.
edges the support of the Primus Research Program of Charles University (PRIMUS/20/SCI/004) and that of Czech Science Foundation (23-07616S).Computational resources were provided by the e-INFRA CZ project (ID:90140), supported by the Ministry of Education, Youth and Sports of the Czech Republic.
at the NNP and ReaxFF level.

Table 2 .
Relative energies ∆E [kJ mol−1] of the proton jump in CHA and FAU at the (∆)NNP and DFT level.

Table 3 .
Relative energies ∆E [kJ mol −1 ] of the Al-O(H) bond dissociation in FAU at the (∆)NNP and DFT level.
[67]isomersfrom (H 2 O) 2 to (H 2 O) 10 , available under: begdb.org)andfourisomers of (H 2 O) 20 .[67]Allclusters were first optimized (constant volume conditions) with a unit cell ensuring a distance between equivalent periodic images of at least 1 nm.Then the aforementioned 210 lattice deformations were applied to all optimized clusters.