# Neural Network-Based Interatomic Potential for the Study of Thermal and Mechanical Properties of Siliceous Zeolites

Luca Brugnoli, Maxime Ducamp, François-Xavier Coudert

► To cite this version:

Luca Brugnoli, Maxime Ducamp, François-Xavier Coudert. Neural Network-Based Interatomic Potential for the Study of Thermal and Mechanical Properties of Siliceous Zeolites. Journal of Physical Chemistry C, 2024, 128 (47), pp.20512-20522. 10.1021/acs.jpcc.4c07365. hal-04780016

HAL Id: hal-04780016

https://hal.science/hal-04780016v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Neural Network-Based Interatomic Potential

for the Study of Thermal and Mechanical Properties of Siliceous Zeolites

Luca Brugnoli,<sup>a,*</sup> Maxime Ducamp,<sup>a</sup> François-Xavier Coudert<sup>a,*</sup>

<sup>a</sup> Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie Paris, 75005 Paris,
France

*Emails: luca.brugnoli1190@gmail.com; fx.coudert@chimieparistech.psl.eu

Abstract

Molecular dynamics simulations of zeolites are commonly employed for the characterization of their
framework dynamics and response to the application of temperature and pressure. While classical
interatomic potentials are commonly used for this task, they offer a description of the interactions in
the system with limited accuracy. Density Functional Theory, meanwhile, is accurate but its high
computational expense limits its scalability for large systems or long dynamics. Recent advances in
machine learning interatomic potentials, trained on computational data obtained at the quantum
chemical level, offer a promising alternative combining high accuracy with computational efficiency.
In this study, we developed an MLIP specifically for pure silica zeolites, trained on data from high-
temperature *ab initio* MD simulations across various zeolitic topologies. This MLIP was then applied
to predict structural properties, thermal expansion, and pressure response of different zeolites,
demonstrating its potential for accurate and generalizable in simulations of topologies beyond its
initial training set.
# 1. Introduction

Zeolites are a family of natural and artificial porous aluminosilicates that rank among the most widely-produced crystalline nanoporous materials in the chemical industry, where they are used in particular for their active role in catalysis, adsorption, and ion exchange. Computational modelling plays a key role in understanding their fundamental properties, linking their physical and chemical behaviour at the macroscopic scale to their structure and dynamics at the atomic scale. With this goal, different levels of chemical description have been used in molecular dynamics (MD) simulations of zeolites, ranging from classical interatomic potentials – with low accuracy but high transferability – to first-principles or *ab initio* schemes based on Density Functional Theory (DFT) – which feature high accuracy along with high computational expense. While DFT calculations are feasible on relatively small zeolitic systems, their high CPU cost makes them unsuitable for applications at large time or length scales, or in high-throughput computational workflows (for example, identifying materials with a specific property by screening a database).

In the specific case of pure silica zeolites, i.e. $SiO_2$ polymorphs, a large number of classical force fields have been used in the literature, with different functional forms or parameterization. Combariza et al. reviewed various existing interatomic potentials¹, and highlighted that accurate atomistic molecular dynamics simulations have to account for a delicate balance of interactions. Yet, at the same time, it is known that a detailed description of the framework dynamics is key in the understanding of the macroscopic behaviour of zeolites: many of their properties arise from the interplay between the motions of the framework with an external constraint, which can be a change in temperature, the application of a mechanical stress, or the adsorption of guest molecules (in the gas or liquid phase). Computational studies have shown that systematic DFT calculations in the linear response regime under the application of pressure,² or in the quasi-harmonic approximation (QHA),³ allow the calculation of mechanical and thermal properties of pure-silica zeolites. However, these
approaches only capture a fraction of the configuration space of zeolitic frameworks, and the
computational cost of the more general *ab initio* MD approach limits its uses.

Machine-learned potentials, or machine learning interatomic potentials (MLIP), have garnered
significant attention in recent years due to their ability to capture complex interatomic interactions
with high accuracy, as well as their computational efficiency. They are data-driven models, typically
deriving the energy and atomic forces of a given molecular configuration through trained neural
networks, rather than relying on a physics-informed and parameterized analytical decomposition of
the interactions at play. They can achieve an accuracy that is close to the reference data that they are
trained on, while retaining sufficient flexibility to describe chemically and structurally complex
environments. Over the past few years, their use has been rapidly growing and they have been applied
to a large and diverse number of systems, ranging from isolated molecules to materials, and from
organic bonds to inorganic coordination.⁴ They have been very recently used, in particular, to describe
silica systems such as amorphous glasses and zeolites. In this work, we have created and published a
data set of trajectories, energies, and forces obtained from high-temperature *ab initio* molecular
dynamics on a diverse set of zeolitic topologies (see **Figure 1**). Trained on this new large-scale data
set, we optimized a machine learning interatomic potential (MLIP) for pure silica zeolites, which we
then used to simulate the structure, the thermal expansion properties and the effects of pressure on a
selection of zeolites, both inside and outside its original training set, in order to probe its
generalizability.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Unit cells and labels of the pure silica zeolite frameworks considered in this work. Green and blue labels are used respectively for the frameworks used for the training of the potential and for an external validation of the potential. FAU, labelled in black, was used as external validation without *ab initio* data reference.</figcaption></figure>

## 2. Computational Methods

We give in this section the key details of the computational procedures followed in the present work. In addition, for full reproducibility of our work, representative input files for all types of calculations performed are made available online at https://github.com/fxcoudert/citable-data

### 2.1. DFT calculations

The *ab initio* data set used for training the machine learning interatomic potential (MLIP) was generated by periodic DFT calculations using the QUICKSTEP/CP2K software.⁵,⁶ This code makes use of a hybrid basis set formed by Gaussian functions and plane waves, where the wave functions
are described by the Gaussian basis with the addition of an auxiliary plane waves basis to describe
the density. The Strongly Constrained and Appropriately Normed (SCAN) meta-GGA functional<sup>7</sup> was
chosen for its good performance across a diverse range of chemical systems and properties of interest,
in particular it reproduces satisfactorily the structural and elastic properties and the relative stability
of different silica polymorphs<sup>8,9</sup>. The inner core of Si and O were described by the Goedecker-Teter-
Hutter (GTH) type of pseudopotential,<sup>10-12</sup> with triple-ζ valence polarized Gaussian basis sets
employed for the valence electrons.

Two types of geometry optimisation were performed on the conventional unit cells of the zeolite
frameworks considered in this study. In one case, starting from the literature geometries (IZA
Structure Commission, Database of Zeolite Structures: https://www.iza-structure.org/databases/), we
obtained four distinct distorted cells from each structure, in which the cell parameters varied by –3%,
–1.5%, 0% and +1.5%. In the other cases, full optimization of the system was achieved by performing
an iterative optimisation cycle, alternating between the optimization of internal coordinates and the
optimization of cell parameters, until both converged. We found this approach more robust in finding
the lowest energy structure compared to the simultaneous optimization of both lattice and internal
coordinates.

In all the latter calculations, the SCF convergence threshold was settled at 10<sup>-9</sup> Hartree, while the
geometry optimization was considered achieved when maximum atomic forces and atomic
displacement were both below 10<sup>-4</sup> Bohr and 10<sup>-4</sup> Hartree/Bohr, respectively. QUICKSTEP makes
use of a multi-grid system to regulate the accuracy of the calculations: 4 grids with a cut-off of 900
Ry on the main grid and a relative cut-off of 100 Ry were enough to optimize most of the structures
and to perform AIMD on them. To fully optimize all the structures, a more accurate multigrid was
required, formed by 5 grids, a main cutoff of 1800 Ry and a relative cutoff of 300 Ry.

**2.2. Ab initio Molecular Dynamics Simulations**
In order to generate a diverse and physically representative set of zeolite configurations, we
performed *ab initio* molecular dynamics (AIMD) simulations, coupling Born-Oppenheimer
dynamics with evaluation of the energy and atomic forces at the DFT level. AIMD simulations were
run on the 4 strained structures of each of the 17 zeolite frameworks considered. QUICKSTEP/CP2K
with the same computational setting adopted in the relaxation of the distorted lattices was employed.
The equations of motion were integrated with a timestep of 0.5 fs.

The unit cell of the systems was fixed and the simulations were carried out in the canonical ensemble (N, V, T) at 1000 K, adopting the CSVR<sup>13</sup> thermostat with a time constant of 1 ps to regulate the temperature. Each system was simulated for a duration of 20 ps. Snapshots of each system were sampled every 20 timesteps, resulting in a total of 1000 configurations per strained structure, which were used to build the data set. We found that during the dynamics, for some strained structures, a few configurations could have very large forces and those outliers would negatively impact the training of the force field. Therefore, configurations with atomic forces higher than 15 eV Å<sup>-1</sup> were excluded from the data set. Among the 64 considered structures, one was found to exhibit dynamically unstable behaviour (the most strained cell of JST), leading to its exclusion from the data set to ensure the integrity and reliability of the training data.

**2.3 Interatomic Potential training**

The Interatomic Potential (IP) was trained with the use of Neural Equivariant Interatomic Potentials (NequIP),¹⁴ an equivariant graph neural network approach for learning interatomic potentials from *ab initio* calculations for molecular dynamics simulations. NequIP represents a significant advance over prior non-equivariant models, by utilizing higher-order tensors instead of relying solely on scalar features. This allows NequIP to capture the complexities of molecular systems more effectively, encoding bond angles, dihedral angles, and other angularly dependent features in a more flexible manner — something of high interest in this work, given the strongly directional nature of the
interactions. Compared to other deep neural network models, this approach is more efficient,
requiring a smaller training set to develop accurate Machine-Learned Interatomic Potentials (MLIPs).

The development of the NequIP model depends on several *hyperparameters*, that determine the accuracy of the potential and the efficiency of the training procedure. In this work, we tuned these hyperparameters to improve the model's performance and capture sufficiently the physics of the systems described. Specifically, before training our final model, we explored the effects of the cutoff radius (*r*<sub>max</sub>), which regulates the range of interatomic interactions considered, and the maximum rotational order (*L*<sub>max</sub>), which limits the rank of the tensor product considered by NequIP. Additionally, we investigated hyperparameters that impact the learning behaviour itself (i.e., the optimization of the ML model): we looked at the influence of the initial learning rate (LR) and the number of interaction blocks in the network on the training speed.

The tuning process involved systematically adjusting the parameters and assessing their effects on
the model's accuracy and computational efficiency. In order to avoid overfitting and limit data
leakage, this optimization of the hyperparameters was conducted on a subset of the training set,
focusing on zeolite frameworks SOD, BEC, and JST. The selected zeolites represent a diverse range
of structures in terms of size and density; they range from small (SOD, 36 atoms) to medium (BEC,
96 atoms) and large (JST, 144 atoms) unit cell zeolites, and in density, span from 1.8 g/cm³ (SOD) to
1.55 g/cm³ (BEC) and to 1.4 g/cm³ (JST). This also allowed us to restrain the computational load
during systematic tests, and efficiently explore a wider range of values for the hyperparameters.

**2.4. Classical molecular dynamics simulations**

We then used our trained interatomic potential (IP) to perform molecular dynamics (MD) simulations
with the LAMMPS code,¹⁵ using the *pair_nequip* pair style provided by the NequIP developers
and available online at https://github.com/mir-group/pair_nequip.
For each zeolite framework, we created a 2 × 2 × 2 supercell based on the crystalline unit cell obtained from the IZA database. The equations of motion were integrated using a timestep of 0.5 fs in all simulations. The following procedure was followed for all classical simulations:

1. Relaxation. Initially, each structure was relaxed through 200 ps of constant-pressure (N, P, T) MD simulation at a temperature of 20 K and pressure of 0 GPa.

2. *Temperature ramping*. The conditions were then switched to constant-volume (N, V, T), and the temperature was gradually increased to the target value over a 10 ps simulation.

3. *Equilibration*. Subsequently, the conditions were switched back to (N, P, T), and 200 ps of dynamics were performed at the target temperature and pressure.

In all (N, P, T) simulations, pressure and temperature were controlled using the Nosé-Hoover thermostat¹⁶,¹⁷ and the Nosé-Hoover barostat¹⁸, with damping time scales of 50 and 500 fs, respectively. In the (N, V, T) simulations, only the Nosé-Hoover thermostat was employed. For simulations at pressures higher than 1 GPa, the use of a pressure ramping approach would increase the stability of the simulation, allowing to reach a higher pressure (from 3 to 10 GPa, depending on the system) almost quasistatically, compared to initializing the system directly at the target pressure (which would fail above the range of 3–4 GPa).

These simulations explored a temperature range between 20 and 900 K for evaluating the thermal properties. Mechanical properties were investigated at 20, 300 and 600 K and in a pressure range of up to 0–10 GPa. Properties of interest were computed as averages from the last 50 ps of each simulation. Furthermore, the elastic constants were computed at 20 K using the ELASTIC_T scripts from LAMMPS, applying a strain of 2% to the simulation box.

The effects of the simulation box size were tested for CHA using unit cells, 2 × 2 × 2, and 3 × 3 × 3 supercells on both thermal expansion and elastic constants, revealing negligible differences.
# Results and Discussion

## 3.1. Force Field Training and Hyperparameters Tuning

Out of the datasets obtained for the 17 frameworks, 11 were used for training the interatomic potential, while the remaining 6 datasets served as a test set to evaluate the performance of the trained potential and its generalizability across unseen zeolitic topologies — an important quality for the exploration of novel materials. To optimize the potential training procedure, a subset of the dataset comprising data from SOD, BEC, and JST was selected. This subset comprised a total of 12,300 configurations, which were subsequently divided into 11,100 configurations for training and 1,200 configurations for validation in the cross-validation scheme. In all the trainings, the loss function was minimized adopting the same scaled weights on forces and energies.

Hyperparameters such as $r_{max}$ (cutoff radius for interatomic interactions) and $L_{max}$ (maximum rotation order considered by NequIP) were tested and optimized. Other hyperparameters used default values as suggested by the NequIP developers. As depicted in **Figures S1** and **S2**, increasing $r_{max}$ improves the accuracy of the potential with a convergence reached with less iterations. Furthermore, higher accuracies and faster convergences are achieved by increasing the number of layer blocks from 3 (**Figure S1**) to 5 (**Figure S2**).

The highest accuracy is obtained with a $r_{max}$ of 6 Å and 5 layer blocks. However, adopting a $r_{max}$ of 5 Å provides a very close accuracy, with an increase of the MAE of the forces from 30.2 to 30.4 meV Å⁻¹ and an increase of the MAE of the energy per atom from 0.130 to 0.164 meV. In addition, this choice saves about the 20% (22 hours) of the total GPU time of the training, which would result in a much more significative GPU time saving on the training on the whole data set, and will result in a more efficient model in production simulations as well.
The results of these tests show that both $r_{max}$ and $L_{max}$ have a significative impact on the final accuracy of the trained model. While a value of $r_{max} = 4$ Å may intuitively seem large enough to describe the atomic environment around Si and O in a zeolitic framework, we find that larger values of 5 and 6 Å lead to a further improvement of both forces and energies, as shown in **Figure S3**. Moreover, by increasing the number of interaction layers in the model from 3 to 5, there is an overall improvement of the performance, while the differences between models obtained with different $r_{max}$ are decreased. Given the similar performances of the models with $r_{max} = 5$ Å and $r_{max} = 6$ Å, and the lower computational cost of the former, we settled on $r_{max} = 5$ Å as a compromise for all the following tests.

Another parameter that affects the accuracy of the model is the value of $L_{max}$ which represents the maximum rotation order considered by NequIP. Typically, good results are obtained with a value of 1, while increasing it to 2 can lead to improved accuracy at the expense of speed and memory requirements. In our evaluation, we assessed the balance between the increase in accuracy achieved with a higher $L_{max}$ and the associated increase in computational cost. Based on the training curves depicted in **Figure S4**, we see that increasing $L_{max}$ to 2 resulted in a modest improvement in the accuracy of the model for forces and energies, with a decrease of approximately 0.9 meV Å⁻¹ and 0.04 meV. However, this improvement came at the cost of a significant increase in simulation time, approximately 90% longer (88 hours vs 171 hours on our hardware). Additionally, it is worth noting that the MAE of the forces from training and validation diverged in the case of the model with $L_{max} = 2$. Considering the limited improvement in accuracy and the substantial increase in computational cost associated with a higher $L_{max}$ value, we decided to keep it fixed at 1.

Lastly, we also examined the effect of the initial learning rate (LR) used during training. From the tuning test conducted we concluded that the most balanced set of parameters to use in the training which included the full data set would be constituted by $r_{max} = 5$ Å, 5 interaction blocks, $L_{max} = 1$, and an initial LR = 0.005.
With the selected hyperparameters, the final training was performed on the full set formed by the
atomic configurations, atomic forces and unit cell energies for 41,000 configurations, split into 37,000
and 4,000 between training and validation, respectively, from the frameworks ANA, BEC, FER, JNT,
JST, LTA, NPT, PUN, RHO, SOD and WEI. The resulting training and validation curves were
reported in **Figure S5**. Each training tasks was performed on a single GPU, the NVIDIA V100 Tensor
Core GPU. The training of the test models required between 4 and 9 days each, while the training of
the final model required 21 days.

All training was originally performed with the main (stable) development branch of NeqiIP, which
does not include the computation of the stress, and the potentials resulting from this branch were
unable to describe the volume variation of the zeolites, and could not be used in constant-pressure
MD simulations. To overcome this shortcoming, the ML potential obtained from the main branch
was updated by performing a single step of a training adopting the stress-branch version of NequIP,
which encoded the evaluation of the virial on the final form of the deployed potential. To check that
this update did not mess up the potential, the errors on the training set were computed with both the
versions, which resulted in negligible differences on all the metrics considered (below the third
decimal point).

The training curves obtained on the larger and more diversified training set shown higher MAE than
those obtained from the initial subset, as expected from increasing the heterogeneity of the dataset.
Nevertheless, the computed errors are well within the range for the MLIPs on siliceous systems.⁸,⁹,¹⁹

## 3.2 Generalizability of the potential

The performance of the model was further assessed on the single dataset for each zeolite-type and on
each of the four considered structure deformations. The whole dataset exhibited a mean absolute error
(MAE) on forces of 32.2 meV Å⁻¹ while for the single frameworks datasets the MAE ranges between
29.8 and 34.4 meV Å⁻¹. Concerning the errors for the energies per particle, the MAE for the training
set is 0.51 meV while for the single framework datasets is of 0.46 meV or below, with the exception
of JNT and WEI with higher MAE values (of 1.38 and 1.51 meV, respectively). As comparison, with
another MLIP recently developed with DeePMD-kit²⁰ on pure siliceous zeolites, the reported MAE
on forces and energies were respectively 39 meV Å⁻¹ and 7.8 meV per particle, on a training set of
about 280,000 configurations from AIMD simulations on 219 different topologies from the IZA
database¹⁹. Another work adopted SchNet²¹ to train a MLIP adopting both quartz polymorphs, well
know zeolites topologies and zeolites extracted from the Deem database²²,²³, reporting MAE values
for energies of 155 meV Å⁻¹ and 3.83 meV per particle.⁹

By inspecting the metrics for the dataset of each deformed framework, we see that the accuracy may
change with the volume of the unit cell, within the same framework type. For most cases, the
deformed structures with the largest strain – lowest volume (labeled “2”) and with the largest volume
(labelled “3”) – exhibit the highest and the lowest MAE on the forces, respectively. A similar trend
can be observed with the MAE of the energies, although with more fluctuations. Within the two
frameworks with the highest MAE, JNT and WEI, the internal differences between the deformed
structures are more striking than the other case.

The performance of the MLIP on the evaluation dataset aligned with that on the training set, showing
good generalization; except for the TON framework, which displays MAE on forces and energy per
particle are respectively 4-5 times and 40-80 times higher than those on the training set, indicating it
as possible outlier. The collected metrics computed on the training and on the evaluations training set
are reported in **Tables S1 and S2**.

The trained MLIP was then assessed on the prediction of the lattice parameters of different zeolitic
frameworks, where crystal shapes were constrained in both static DFT and MD simulations. **Tables
1 and 2** reports the comparison between the results of the two methods. We observe that unit cell
volumes are predicted with a MAE of 1.21%, considering structures from both training and test set,
while errors lower than 0.5% were achieved for 8 out of the 17 frameworks. This displays better accuracy than empirical zeolite force fields used in the field, indicating the advantage of the ab initio data-based force field optimization.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Correlation plots between the relaxed volume computed at the DFT (SCAN) level and with the trained MLIP for the different framework used in the training set (blue dots) and as evaluation set (red dots).</figcaption></figure>

**Figure 2** shows the DFT and the MLIP computed equilibrium volumes for the different frameworks, including both the training and evaluation sets. We note here that the cubic system RHO exhibited the largest deviation from the DFT results, with a volume deviation of 7.46%. This is a special (and interesting) case because DFT optimization revealed two local minima structures with different volumes, while the MD relaxation resulted in a structure with a volume close to the metastable, low-
volume configuration. Despite this complication, the trained MLIP demonstrated reasonably accurate performance, considering it was trained solely on high-energy structures without specific consideration of the energy minima.

<table><caption>**Table 1:** Comparison between the lattice parameters and volume of the zeolite frameworks used in the training set obtained from DFT minimization and from MLIP molecular dynamic relaxation.</caption><thead><tr><th>Framework</th><th>Method</th><th>a (Å)</th><th>b (Å)</th><th>c (Å)</th><th>α (°)</th><th>β (°)</th><th>γ (°)</th><th>V (Å<sup>3</sup>)</th><th>ΔV/V (%)</th></tr></thead><tbody><tr><td rowspan="2">JNT</td><td>MLIP</td><td>8.231</td><td>14.179</td><td>13.495</td><td>90.00</td><td>87.125</td><td>90.00</td><td>1573</td><td rowspan="2">-0.37</td></tr><tr><td>DFT</td><td>8.237</td><td>14.167</td><td>13.602</td><td>90.00</td><td>87.125</td><td>90.00</td><td>1579</td></tr><tr><td rowspan="2">NPT</td><td>MLIP</td><td>13.940</td><td>13.940</td><td>13.940</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2709</td><td rowspan="2">-0.45</td></tr><tr><td>DFT</td><td>13.961</td><td>13.961</td><td>13.961</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2721</td></tr><tr><td rowspan="2">PUN</td><td>MLIP</td><td>14.713</td><td>8.659</td><td>18.972</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2417</td><td rowspan="2">0.44</td></tr><tr><td>DFT</td><td>14.647</td><td>8.621</td><td>19.058</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2406</td></tr><tr><td rowspan="2">WEI</td><td>MLIP</td><td>11.782</td><td>10.260</td><td>9.972</td><td>90.00</td><td>90.00</td><td>90.00</td><td>1205</td><td rowspan="2">-1.79</td></tr><tr><td>DFT</td><td>11.919</td><td>10.346</td><td>9.954</td><td>90.00</td><td>90.00</td><td>90.00</td><td>1227</td></tr><tr><td rowspan="2">ANA</td><td>MLIP</td><td>13.366</td><td>13.366</td><td>13.367</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2388</td><td rowspan="2">-0.06</td></tr><tr><td>DFT</td><td>13.369</td><td>13.369</td><td>13.369</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2389</td></tr><tr><td rowspan="2">LTA</td><td>MLIP</td><td>11.826</td><td>11.826</td><td>11.826</td><td>90.00</td><td>90.00</td><td>90.00</td><td>1654</td><td rowspan="2">-0.21</td></tr><tr><td>DFT</td><td>11.834</td><td>11.834</td><td>11.834</td><td>90.00</td><td>90.00</td><td>90.00</td><td>1657</td></tr><tr><td rowspan="2">JST</td><td>MLIP</td><td>15.072</td><td>15.072</td><td>15.072</td><td>90.00</td><td>90.00</td><td>90.00</td><td>3424</td><td rowspan="2">1.96</td></tr><tr><td>DFT</td><td>14.975</td><td>14.975</td><td>14.975</td><td>90.00</td><td>90.00</td><td>90.00</td><td>3358</td></tr><tr><td rowspan="2">FER</td><td>MLIP</td><td>18.502</td><td>13.915</td><td>7.336</td><td>90.00</td><td>90.00</td><td>90.00</td><td>1888</td><td rowspan="2">-0.43</td></tr><tr><td>DFT</td><td>18.575</td><td>13.926</td><td>7.333</td><td>90.00</td><td>90.00</td><td>90.00</td><td>1897</td></tr><tr><td rowspan="2">BEC</td><td>MLIP</td><td>12.643</td><td>12.643</td><td>12.848</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2054</td><td rowspan="2">-2.45</td></tr><tr><td>DFT</td><td>12.606</td><td>12.606</td><td>13.250</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2106</td></tr><tr><td rowspan="2">SOD</td><td>MLIP</td><td>8.661</td><td>8.661</td><td>8.661</td><td>90.00</td><td>90.00</td><td>90.00</td><td>649.8</td><td rowspan="2">1.52</td></tr><tr><td>DFT</td><td>12.606</td><td>12.606</td><td>13.250</td><td>90.00</td><td>90.00</td><td>90.00</td><td>640</td></tr><tr><td rowspan="2">RHO</td><td>MLIP</td><td>14.524</td><td>14.524</td><td>14.524</td><td>90.00</td><td>90.00</td><td>90.00</td><td>3064</td><td rowspan="2">-7.46</td></tr><tr><td>DFT</td><td>14.905</td><td>14.905</td><td>14.905</td><td>90.00</td><td>90.00</td><td>90.00</td><td>3311</td></tr></tbody></table>
<table><caption>**Table 2:** Comparison between the lattice parameters and volume of the zeolite frameworks used in the evaluation set obtained from DFT minimization and from MLIP molecular dynamic relaxation.</caption><thead><tr><th>Framework</th><th>Method</th><th>a (Å)</th><th>b (Å)</th><th>c (Å)</th><th>α (°)</th><th>β (°)</th><th>γ (°)</th><th>V (Å³)</th><th>ΔV/V (%)</th></tr></thead><tbody><tr><td rowspan="2">MOR</td><td>MLIP</td><td>17.865</td><td>20.095</td><td>7.381</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2650</td><td rowspan="2">-0.31</td></tr><tr><td>DFT</td><td>17.789</td><td>20.173</td><td>7.406</td><td>90.00</td><td>90.00</td><td>90.00</td><td>2658</td></tr><tr><td rowspan="2">MFI</td><td>MLIP</td><td>20.024</td><td>19.674</td><td>13.099</td><td>90.00</td><td>90.00</td><td>90.00</td><td>5204</td><td rowspan="2">-0.45</td></tr><tr><td>DFT</td><td>20.103</td><td>19.596</td><td>13.211</td><td>90.00</td><td>90.00</td><td>90.00</td><td>5228</td></tr><tr><td rowspan="2">GIS</td><td>MLIP</td><td>9.539</td><td>9.539</td><td>9.886</td><td>90.00</td><td>90.00</td><td>90.00</td><td>900</td><td rowspan="2">-0.51</td></tr><tr><td>DFT</td><td>9.570</td><td>9.570</td><td>9.872</td><td>90.00</td><td>90.00</td><td>90.00</td><td>904</td></tr><tr><td rowspan="2">CHA</td><td>MLIP</td><td>9.285</td><td>9.285</td><td>9.285</td><td>93.87</td><td>93.87</td><td>93.87</td><td>795</td><td rowspan="2">0.06</td></tr><tr><td>DFT</td><td>9.286</td><td>9.286</td><td>9.286</td><td>94.06</td><td>94.06</td><td>94.06</td><td>794</td></tr><tr><td rowspan="2">HEU</td><td>MLIP</td><td>17.346</td><td>17.466</td><td>7.326</td><td>90.00</td><td>116.10</td><td>90.00</td><td>1993</td><td rowspan="2">0.82</td></tr><tr><td>DFT</td><td>17.246</td><td>17.550</td><td>7.310</td><td>90.00</td><td>116.68</td><td>90.00</td><td>1977</td></tr><tr><td rowspan="2">TON</td><td>MLIP</td><td>13.536</td><td>17.122</td><td>5.044</td><td>90</td><td>90</td><td>90</td><td>1169</td><td rowspan="2">-1.32</td></tr><tr><td>DFT</td><td>14.012</td><td>17.114</td><td>4.940</td><td>90</td><td>90</td><td>90</td><td>1185</td></tr></tbody></table>

### 3.3 Thermal expansion simulations

Understanding the thermal expansion behaviour of zeolites has important implications for practical applications, in part for process engineering but also because the diffusivity and reactivity of the guest molecules in the structure are strongly influenced by the size of the pores, which can contract or expand with variations in temperature.

Zeolites are known to be one of the few classes of materials that can feature negative thermal expansion (NTE), where the material contracts upon heating, a difference from most materials which expand upon heating. In zeolites, as in some other framework materials, the NTE behaviour has been attributed to the presence of Rigid Unit Modes (RUMs) originating from their corner-sharing tetrahedral structure. When heated, these RUMs permit the 'flexing' or rotation of these tetrahedral units within the zeolitic framework, causing an overall contraction of the material – thus leading to NTE. Essentially, the thermal energy allows the rotation of the corner-sharing tetrahedra, resulting in a smaller overall volume even as temperature rises. Both experimental and in-silico investigations found exceptions to this behaviour, i.e., zeolites with positive thermal expansions (PTE), but since
their rarity and the presence of RUMs in all of these frameworks, is still debate that PTE in these
systems is due to the nano-porous nature of the material instead of the framework itself.

The characterization of the thermal behaviour of zeolites is made experimentally difficult by the
combination of several factors. First, the analysis of thermal expansion requires the presence of single
crystal data, which are difficult to obtain on fully evacuated/activated samples; in some cases in the
literature, organic templates or solvent molecules are still present in the pores, which affects the
thermal behaviour. Secondly, several zeolitic frameworks undergo phase transitions upon heating,
which makes their thermal expansion harder to analyse. Therefore, the computational study of the
thermal behaviour of zeolites is an area of active research. In a previous work in our group, the thermal
expansion coefficient of the pure-silica zeolites were systematically computed through DFT in the
quasi-harmonic approximation (QHA)³, i.e., in a low-temperature regime, finding that all the
investigated frameworks displayed NTE within the temperature range of 10–300 K. The exploitation
of DFT and QHA to investigate the thermal properties of the whole zeolitic database is a very
demanding task from a computational point of view, therefore an alternative cheaper and as accurate
approach is desired. The MLIP was trained on trajectories at different volumes and degrees of strain,
therefore the rotation of the RUMs should be captured by the model.

We investigated the thermal expansion behaviour of pure-silica zeolite frameworks FAU, CHA, LTA
and SOD through a series of NPT simulations ranging from 20 to 1000 K. We note that volume
fluctuations during the MD increases with the temperature, making the results less reliable in the
higher end of the temperature range. The thermal expansion curves and computed thermal expansion
coefficient are reported respectively in **Figure 3** and **Table 3** and we now discuss them individually
against the known experimental data.

**SOD.** Sodalite shows positive thermal expansion (PTE) up to 800 K, then expansion stops and it starts to contract between 900 and 1000 K, in apparent agreement with observations which found the inversion of the expansion behaviour at about 700 K²⁴,²⁵. It must be noted that while experimental
measurements found a PTE for sodalite minerals in this wide range of T,²⁶ the study of Leardini *et al*²⁴ found for pure silica sodalite PTE between 273 and 573 K and NTE above 653 K. These differences were ascribed to the neutral nature of the organic template they used in the synthesis, which once calcinated did not interfere with the framework, allowing the contraction of the structure.

**LTA.** For the pure-silica LTA zeolite (experimentally called ITQ-29), PTE was observed from 20 K to 400 K; then the volume is stationary in the range 400-600 K, then NTE from 600 K to 1000 K. Unfortunately, the only experimental results reported on the thermal expansion of LTA are limited to the low temperature regime, between 100-300 K, where a strong NTE was observed, with a mean thermal expansion coefficient (300 to 100 K) of –22.1 10⁻⁶ K⁻¹.²⁷ At this stage, it is unclear whether the origin of this discrepancy at low temperature comes from an inadequacy of the MLIP (which was trained mostly on high-temperature data) or from the DFT calculations themselves.

**CHA.** For chabazite, NTE is found for the whole T range explored and show the largest (negative) thermal expansion coefficient $\alpha_V$ among the investigated frameworks –13.07 10⁻⁶ K⁻¹ between 300 and 900 K, in qualitative agreement with the experimental value (–28.5 10⁻⁶ K⁻¹).²⁸,²⁹ At difference with the experiments, simulations found a larger expansion rate at low temperature. Moreover, CHA, being a non-cubic framework, has an anisotropic expansion.

**FAU.** For pure-silica faujasite, NTE is observed up to 900 K, but of about an order of magnitude lower than CHA in the range 20-900 K range, in agreement with experiments³⁰,³¹ Above 900 K the behaviour change to PTE.
<table><caption>**Table 3:** Computed volumetric and linear thermal expansion coefficients for different temperature ranges. Experimental values are reported in parentheses for the appropriate temperature range.</caption><thead><tr><th>Zeolite</th><th>Temp Range (K)<br>(Exp. Range)</th><th>α<sub>v</sub> (MK<sup>-1</sup>)</th><th>α<sub>a</sub>(MK<sup>-1</sup>)</th><th>α<sub>b</sub> (MK<sup>-1</sup>)</th><th>α<sub>c</sub> (MK<sup>-1</sup>)</th></tr></thead><tbody><tr><td>CHA</td><td>20-300</td><td>-28.12</td><td>-14.46</td><td>-14.90</td><td>1.21</td></tr><tr><td>CHA</td><td>300-600</td><td>-15.67</td><td>-6.30</td><td>-6.57</td><td>-2.79</td></tr><tr><td>CHA</td><td>600-900</td><td>-10.52</td><td>-5.34</td><td>-14.53</td><td>-4.49</td></tr><tr><td>CHA</td><td>900-1000</td><td>-8.24</td><td>-2.33</td><td>-5.77</td><td>-0.21</td></tr><tr><td>CHA</td><td>20-1000</td><td>-16.80</td><td>-7.59</td><td>-10.81</td><td>-1.45</td></tr><tr><td>CHA</td><td>300-900 (293–873<sup>32</sup>, 308–753<sup>29</sup>)</td><td>-13.07 (-28.5<sup>32</sup>, -21.22<sup>29</sup>)</td><td>-5.29 (-8.24<sup>32</sup>, -9.11<sup>29</sup>)</td><td>-9.79 (-8.24<sup>32</sup>, -9.11<sup>29</sup>)</td><td>-2.89 (-13.3<sup>32</sup>, -3.08<sup>29</sup>)</td></tr><tr><td>SOD</td><td>20-300</td><td>35.87</td><td>11.96</td><td>11.96</td><td>11.96</td></tr><tr><td>SOD</td><td>300-600</td><td>62.83</td><td>20.94</td><td>20.94</td><td>20.94</td></tr><tr><td>SOD</td><td>600-900</td><td>43.57</td><td>14.52</td><td>14.52</td><td>14.52</td></tr><tr><td>SOD</td><td>900-1000 (1033–1123<sup>24</sup>)</td><td>-22.11 (-14.0<sup>24</sup>)</td><td>-7.37 (-4.7<sup>24</sup>)</td><td>-7.37 (-4.7<sup>24</sup>)</td><td>-7.37 (-4.7<sup>24</sup>)</td></tr><tr><td>SOD</td><td>20-900</td><td>39.72</td><td>13.24</td><td>13.24</td><td>13.24</td></tr><tr><td>SOD</td><td>20-1000</td><td>34.85</td><td>11.62</td><td>11.62</td><td>11.62</td></tr><tr><td>LTA</td><td>20-300 (100–300<sup>33</sup>)</td><td>12.34 (-22.1<sup>33</sup>)</td><td>4.11 (-7.37<sup>33</sup>)</td><td>4.11 (-7.37<sup>33</sup>)</td><td>4.11 (-7.37<sup>33</sup>)</td></tr><tr><td>LTA</td><td>300-600</td><td>4.92</td><td>1.64</td><td>1.64</td><td>1.64</td></tr><tr><td>LTA</td><td>600-900</td><td>0.10</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>LTA</td><td>900-1000</td><td>-3.81</td><td>-1.27</td><td>-1.27</td><td>-1.27</td></tr><tr><td>LTA</td><td>20-1000</td><td>3.20</td><td>1.07</td><td>1.07</td><td>1.07</td></tr><tr><td>FAU</td><td>20-300 (25-298<sup>30</sup>)</td><td>-5.40 (-15.8<sup>30</sup>)</td><td>-1.80 (-2.5<sup>30</sup>)</td><td>-1.80 (-2.5<sup>30</sup>)</td><td>-1.80 (-2.5<sup>30</sup>)</td></tr><tr><td>FAU</td><td>20-600 (25-573<sup>30</sup>)</td><td>-4.93 (-12.6<sup>30</sup>)</td><td>-1.64 (-4.2<sup>30</sup>)</td><td>-1.64 (-4.2<sup>30</sup>)</td><td>-1.64 (-4.2<sup>30</sup>)</td></tr><tr><td>FAU</td><td>300-600 (298-573)</td><td>-0.23 (-13.6)</td><td>-0.08 (-4.6)</td><td>-0.08 (-4.6)</td><td>-0.08 (-4.6)</td></tr><tr><td>FAU</td><td>300-900 (298–923<sup>31</sup>)</td><td>-0.76 (-11.8<sup>31</sup>)</td><td>-0.25 (-4.0<sup>31</sup>)</td><td>-0.25 (-4.0<sup>31</sup>)</td><td>-0.25 (-4.0<sup>31</sup>)</td></tr><tr><td>FAU</td><td>900-1100 (923–1123<sup>31</sup>)</td><td>1.76 (30)</td><td>0.59 (10)</td><td>0.59 (10)</td><td>0.59 (10<sup>31</sup>)</td></tr><tr><td>FAU</td><td>20-1000</td><td>-2.33</td><td>-0.78</td><td>-0.78</td><td>-0.78</td></tr></tbody></table>

<figure><img src="image_4.png" /></figure>
**Figure 3:** Plot of the percentage variation of the volume with the temperature, between 20 and 1000K, for the different frameworks considered.

**3.4 Isothermal compression simulations: phase transition and amorphization**

Due to their high porosity, zeolites undergo very important structural changes in response to the application of mechanical pressure. There are several examples of phase transitions, involving tilting of the constituent tetrahedra or changes in the ring structures. Pressure may also induce pore collapse and transition to an amorphous state, which can retain some topological order from the parent crystalline phase and is distinct from common silicate glasses.

We investigated the structural effects of compression on the frameworks SOD, CHA, and FAU through a series of (N, P, T) simulations at increasing values of applied pressure. The influence of temperature was evaluated by performing these series of simulations at 20 K, 300 K and 600 K. The range of pressure investigated shows a variation of the cell volume that can exceed that imposed to the structure for the generation of the training data (i.e., between +1.5% and -3.0%). We also performed calculations of the second-order elastic tensor at 20 K, in order to demonstrate the potential for the MLIP to describe anisotropic mechanical properties, and to investigate the influence of system size by comparing 2 x 2 x 2 supercells and 3 x 3 x 3 supercells. These results are reported and discussed in the Supporting Information.

**SOD.** The mechanical properties of siliceous sodalite were experimentally found to be significantly influenced by the nature of the contained template species, even in cases of supposedly weak interactions.<sup>34-37</sup> For instance the experimental bulk modulus $B_0$ ranges between to 24 GPa<sup>37</sup> and 132 GPa<sup>36</sup>. This variability is justified by the existence of two low-pressure phases with different compressibility, whose relative stability at ambient pressure is influenced by the template species. *Ab initio* investigations at the PW91 DFT level on ideal, template-free silica sodalite found that the zero-pressure phase belongs to the high symmetry $Im\bar{3}m$ space group and undergoes a phase transition at
0.12 GPa to the Im$\overline{3}$m group, with a loss of inversion symmetry.³⁸ Additionally, two other phase
transitions were identified with further compression up to about 4 GPa, preserving the lattice
symmetry I$\overline{4}$3m. The I$\overline{4}$3m phase was found to be less compressible, with a computed bulk modulus
of 113 GPa, closely resembling the experimental value of ~132 GPa. The three high-pressure I$\overline{4}$3m
phases became progressively more compressible, with the bulk modulus decreasing from 67 to 15
GPa. However, the agreement with experimental results is only partial: for example, in the case of
dioxolane silica SOD, a transition from cubic to rhombohedral was observed at pressures higher than
1.2 GPa, with the trigonal space group being either R3m or R3.³⁷ This discrepancy with the *ab initio*
calculations is likely due to the constraints adopted in the simulations, which preserved the cubic
symmetry of the lattice. All the studies indicated the increase in the tilting of the silica tetrahedra as
main mechanism of compression underlying the phase-transition Im$\overline{3}$m/ I$\overline{4}$3m, which is supposed to
be continuous. On all-silica SOD, the pressure-induced amorphization has not been reported in any
experimental work.

Our simulations studies found that the behaviour under pressure of sodalite changes with the
temperature: no phase transitions were observed at 20 K in the whole range of P explored (from –1
to 10 GPa), and the structure of minima is located very close to zero pressure, as show in **Figure 4**.
Similar results were found at 300 K, with the minima structure shifted at about 0.2 to 0.3 GPa. At 600
K, the low-pressure part of the PES curve is difficult to interpret, but a continuous transformation is
observed below 0.5 GPa, with a decrease of volume of about 8% (the experimental one of the
dioxolane silica sodalite in the same P range at 300 K is 3%)³⁷. No further transitions were found at
higher pressure: the isotropic constraints adopted in these simulations prevent the breaking of the
lattice symmetry, thus the observation of the experimentally characterized cubic to rhombohedral
transition. A recent work which used the Consistent Valence Force Field characterized the pressured
induced amorphization of SOD at 300 K, finding it occurring as a discontinuous transition a 6.81
GPa.³⁹
As shown in the **Table 4**, bulk modulus at 20 K and 300 K (obtained through third-order Birch-Murnaghan equation of state, see Methods) are very close, 25.0 and 24.6 GPa, whereas at 600 K it decreases to 10.0 GPa. The computed low temperature values are a good match for the experimental values obtained for the lower symmetry I4\3m phase than the ground state, larger, Im3\3m phase. The negative value of $B_0'$ at 600 K underlines the characterized low-pressure phase transition.

<table><caption>**Table 4:** Bulk modulus and pressure derivative of the bulk modulus obtained by fit of the third order Birch-Murnaghan equation of states, and the pressure at which a transition was observed.</caption><thead><tr><th>Framework</th><th><i>T</i> (K)</th><th><i>B</i><sub>0</sub> (GPa)</th><th><i>B</i><sub>0</sub>'</th><th><i>P</i><sub>t</sub> (GPa)</th></tr></thead><tbody><tr><td rowspan="3">SOD</td><td>20</td><td>25.0</td><td>3.43</td><td>—</td></tr><tr><td>300</td><td>24.6</td><td>6.94</td><td>—</td></tr><tr><td>600</td><td>10.0</td><td>-7.44</td><td>0 - 0.5</td></tr><tr><td rowspan="3">CHA</td><td>20</td><td>20.0</td><td>10.65</td><td>1.4</td></tr><tr><td>300</td><td>47.4</td><td>21.90</td><td>1.0</td></tr><tr><td>600</td><td>52.1</td><td>30.21</td><td>0.8 - 1.0</td></tr><tr><td rowspan="3">FAU</td><td>20</td><td>61.7</td><td>13.78</td><td>0.8 - 3.4</td></tr><tr><td>300</td><td>45.3</td><td>16.40</td><td>1.8 - 3.0</td></tr><tr><td>600</td><td>37.0</td><td>7.54</td><td>2.6 - 3.0</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Pressure/Density (a), Pressure vs. Energy, b) and Volume vs. Energy (c) curves for the SOD framework in the -1 to +2 GPa range.</figcaption></figure>
**FAU.** The experimental characterization under pressure of the purely siliceous faujasite is quite limited. One work investigated the effect of pressure on the IR spectra of the highly siliceous zeolite Y (Si:Al = 100), finding a phase transition from about 1.7 GPa to 2.2 GPa associated with the amorphization of the network<sup>40</sup>. The hypothesized mechanism of amorphization involve primarily the deformation of the double six-membered rings forming the sodalite cages in the network, up to rupture or severe distortion.

Another work characterized the elastic properties of a pure siliceous faujasite, deriving a bulk modulus $B_0$ of 38 GPa, the same of α-quartz, and the start of the amorphization at 2.2 GPa.<sup>41</sup> Our simulations revealed a temperature-dependent relationship for the onset pressure of amorphization and the pressure range within which it occurs, as shown in **Figure 5**. At 20 K, amorphization initiates between 0.8 and 3.2 GPa; beyond this pressure, the structure collapses, resulting in a 10% volume reduction, followed by unstable dynamics upon further pressure increases. At 300 K, amorphization occurs within a narrower pressure range of 1.8-2.6 GPa, preceding structural collapse. At 600 K, the transition appears to take place between 2.4 and 2.6 GPa, just before the structural collapse. The computed $B_0$ values, fitted between -1 and 1.6 GPa, shown an increasing of the compressibility of the zeolite with $T$, from 61.74 to 45.33 to 37.03 GPa in the range of 20 to 600 K, with the latter values in good agreement with the experimental value.
<figure><img src="image_6.png" /><figcaption>**Figure 5:** (a) Pressure/Density, (b) Pressure/Energy, and (c) Volume/Energy curves for the FAU framework in the -1 to +3.5 GPa range.</figcaption></figure>

**CHA.** The experimental characterization on the pure silica chabazite under pressure has found a phase transition at room temperature from the rhombohedral space group R$\overline{3}$m to the triclinic P$\overline{1}$ occurring
at 1.42 GPa.⁴² The triclinic symmetry of the cell is preserved up to 3.47 GPa, then the complete, irreversible amorphization is reached at 5.33 GPa.

The simulation results for chabazite under pressure revealed a phase transition within the 0.8 to 1.4 GPa range (see **Figure 6**). In the case of CHA, temperature influences the processes, albeit with a more complex interpretation due to the anisotropic structure. As depicted in **Figure 7**, at 20 K, compressibility shifts from negative to positive strain, with compression affecting solely the *c*-axis contraction of the hexagonal cell until reaching the phase transition at 1.4 GPa. This transition is characterized by a volume decrease of approximately 15% and a reduction in the *c*-axis, accompanied by a minor expansion of lattice parameters *a* and *b*, indicative of auxetic behaviour. Beyond this pressure, the structure continues to contract predominantly along the *c*-axis, leading to simulation failure at pressures exceeding 6 GPa. Similarly, at higher *T*, volume contracts prior to the phase transition, occurring between 0.6 and 1.0 GPa at 300 K, and between 0.8 and 1.0 GPa at 600 K, predominantly affecting the *c*-axis. Post-transition, the *a* and *b* lattice parameters slightly increase, yet under additional stress, they exhibit more significant contraction than at lower temperatures. The calculated value of *B*₀ indicates a reduction in compressibility with increasing temperature but not experimental data is available on this for comparison.
<figure><img src="image_7.png" /><figcaption>**Figure 6:** (a) Pressure/Density, (b) Pressure/ Energy, and (c) Volume/Energy curves for the CHA framework in the -1 to +6 GPa range.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Ratios *V*₀/*V*, *a*₀/*a*, *b*₀/*b*, *c*₀/*c* plotted against the applied pressure for the CHA framework at different values of temperature.</figcaption></figure>

**3. Conclusions**

In this work, we developed a machine learning interatomic potential (MLIP) for pure silica zeolites, trained on a diverse set of high-temperature *ab initio* molecular dynamics data encompassing various zeolitic topologies. The MLIP was successfully utilized to simulate the structure, thermal expansion properties, and the effects of pressure on selected zeolites, demonstrating high accuracy and
transferability. Notably, the ability of the MLIP to generalize to other zeolitic frameworks was confirmed by its accurate performance on zeolite structures not included in the training set, indicating its potential applicability to a wide range of pure silica zeolite frameworks. The generalizability of our MLIP suggests a promising path for the exploration and development of novel zeolitic materials, with broad implications for the chemical industry and the nanoporous materials modelling community. Indeed, molecular dynamics simulations of zeolites are commonly employed for the characterization of their framework dynamics and response to the application of temperature and pressure, but existing classical interatomic potentials can only describe interactions with limited accuracy, while Density Functional Theory, meanwhile, is accurate but has very high computational cost. The use of MLIP offers an interesting middle-ground choice in the computational toolbox of materials scientists and physical chemists. This will allow, in future work, to revisit the questions of high-throughput screening of zeolites for energetic, mechanical, and thermal stability, as well as the question (still widely open) of experimental feasibility of zeolitic frameworks.

**Acknowledgments.**

This work was supported by the Agence Nationale de la Recherche (project MATAREB, ANR-18-CE29-0009-01), and access to HPC platforms was provided by a GENCI grant (A0150807069).

**Supporting Information.** Tables of error metrics, learning curves, parity plots, thermal and mechanical properties, influence of system size.

**References**

(1) Combariza, A. F.; Gomez, D. A.; Sastre, G. Simulating the Properties of Small Pore Silicazeolites Using Interatomic Potentials. *Chem. Soc. Rev.* **2013**, *42* (1), 114–127. https://doi.org/10.1039/C2CS35243E.

(2) Coudert, F.-X. Systematic Investigation of the Mechanical Properties of Pure Silica Zeolites: Stiffness, Anisotropy, and Negative Linear Compressibility. *Phys. Chem. Chem. Phys.* **2013**, *15* (38), 16012. https://doi.org/10.1039/c3cp51817e.
(3) Ducamp, M.; Coudert, F.-X. Systematic Study of the Thermal Properties of Zeolitic Frameworks. *J. Phys. Chem. C* **2021**, *125* (28), 15647–15658. https://doi.org/10.1021/acs.jpcc.1c03975.

(4) Mortazavi, B.; Silani, M.; Podryabinkin, E. V.; Rabczuk, T.; Zhuang, X.; Shapeev, A. V. First-Principles Multiscale Modeling of Mechanical Properties in Graphene/Borophene Heterostructures Empowered by Machine-Learning Interatomic Potentials. *Advanced Materials* **2021**, *33* (35), 2102807. https://doi.org/10.1002/adma.202102807.

(5) Kühne, T. D. et al. CP2K: An Electronic Structure and Molecular Dynamics Software Package - Quickstep: Efficient and Accurate Electronic Structure Calculations. *The Journal of Chemical Physics* **2020**, *152* (19), 194103. https://doi.org/10.1063/5.0007045.

(6) Lippert, G.; Hutter, J.; Parrinello, M. A Hybrid Gaussian and Plane Wave Density Functional Scheme. *Mol. Phys.* **1997**, *92* (3), 477–487. https://doi.org/10.1080/00268979709482119.

(7) Sun, J.; Ruzsinszky, A.; Perdew, J. P. Strongly Constrained and Appropriately Normed Semilocal Density Functional. *Phys. Rev. Lett.* **2015**, *115* (3), 036402. https://doi.org/10.1103/PhysRevLett.115.036402.

(8) Erhard, L. C.; Rohrer, J.; Albe, K.; Deringer, V. L. A Machine-Learned Interatomic Potential for Silica and Its Relation to Empirical Models. *npj Comput Mater* **2022**, *8* (1), 1–12. https://doi.org/10.1038/s41524-022-00768-w.

(9) Erlebach, A.; Nachtigall, P.; Grajciar, L. Accurate Large-Scale Simulations of Siliceous Zeolites by Neural Network Potentials. *npj Comput Mater* **2022**, *8* (1), 1–12. https://doi.org/10.1038/s41524-022-00865-w.

(10) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. *Phys. Rev. B* **1996**, *54* (3), 1703–1710. https://doi.org/10.1103/PhysRevB.54.1703.

(11) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. *Phys. Rev. B* **1998**, *58* (7), 3641–3662. https://doi.org/10.1103/PhysRevB.58.3641.

(12) Krack, M. Pseudopotentials for H to Kr Optimized for Gradient-Corrected Exchange-Correlation Functionals. *Theor Chem Acc* **2005**, *114* (1–3), 145–152. https://doi.org/10.1007/s00214-005-0655-y.

(13) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. *The Journal of Chemical Physics* **2007**, *126* (1), 014101. https://doi.org/10.1063/1.2408420.

(14) Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E(3)-Equivariant Graph Neural Networks for Data-Efficient and Accurate Interatomic Potentials. *Nat Commun* **2022**, *13* (1), 2453. https://doi.org/10.1038/s41467-022-29939-5.

(15) Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in 't Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a Flexible Simulation Tool for Particle-Based Materials Modeling at the Atomic, Meso, and Continuum Scales. *Computer Physics Communications* **2022**, *271*, 108171. https://doi.org/10.1016/j.cpc.2021.108171.

(16) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. *J. Chem. Phys.* **1984**, *81* (1), 511–519. https://doi.org/10.1063/1.447334.

(17) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. *Phys. Rev. A* **1985**, *31* (3), 1695–1697. https://doi.org/10.1103/PhysRevA.31.1695.

(18) Melchionna, S.; Ciccotti, G.; Lee Holian, B. Hoover NPT Dynamics for Systems Varying in Shape and Size. *Molecular Physics* **1993**, *78* (3), 533–544. https://doi.org/10.1080/00268979300100371.

(19) Sours, T. G.; Kulkarni, A. R. Predicting Structural Properties of Pure Silica Zeolites Using Deep Neural Network Potentials. *J. Phys. Chem. C* **2023**, *127* (3), 1455–1463. https://doi.org/10.1021/acs.jpcc.2c08429.

(20) Wang, H.; Zhang, L.; Han, J.; E, W. DeePMD-Kit: A Deep Learning Package for Many-Body Potential Energy Representation and Molecular Dynamics. *Computer Physics Communications* **2018**, *228*, 178–184. https://doi.org/10.1016/j.cpc.2018.03.016.

(21) Schütt, K. T.; Arbabzadah, F.; Chmiela, S.; Müller, K. R.; Tkatchenko, A. Quantum-Chemical Insights from Deep Tensor Neural Networks. *Nat Commun* **2017**, *8* (1), 13890. https://doi.org/10.1038/ncomms13890.

(22) Pophale, R.; Cheeseman, P. A.; Deem, M. W. A Database of New Zeolite-like Materials. *Phys. Chem. Chem. Phys.* **2011**, *13* (27), 12407–12412. https://doi.org/10.1039/C0CP02255A.

(23) Michael Deem's PCOD and PCOD2 Databases of Zeolitic Structures. https://doi.org/10.5281/zenodo.4030231.
(24) Leardini, L.; Martucci, A.; Cruciani, G. The Unusual Thermal Expansion of Pure Silica Sodalite Probed by in Situ Time-Resolved Synchrotron Powder Diffraction. *Microporous and Mesoporous Materials* **2012**, *151*, 163–171. https://doi.org/10.1016/j.micromeso.2011.10.042.

(25) Zheng, Z.; Guliants, V. V.; Misture, S. Sodalites as Ultramicroporous Frameworks for Hydrogen Separation at Elevated Temperatures: Thermal Stability, Template Removal, and Hydrogen Accessibility. *J Porous Mater* **2009**, *16* (3), 343–347. https://doi.org/10.1007/s10934-008-9206-y.

(26) Taylor, D. The Thermal Expansion of the Sodalite Group of Minerals. *Mineral. mag. j. Mineral. Soc.* **1968**, *36* (282), 761–769. https://doi.org/10.1180/minmag.1968.036.282.02.

(27) Carey, T.; Corma, A.; Rey, F.; Tang, C. C.; Hriljac, J. A.; Anderson, P. A. The Effect of Extra Framework Species on the Intrinsic Negative Thermal Expansion Property of Zeolites with the LTA Topology. *Chem. Commun.* **2012**, *48* (47), 5829. https://doi.org/10.1039/c2cc30582h.

(28) Woodcock, D. A.; Lightfoot, P.; Villaescusa, L. A.; Diaz-Cabañas, M.-J.; Camblor, M. A.; Engberg, D. Negative Thermal Expansion in the Siliceous Zeolites Chabazite and ITQ-4: A Neutron Powder Diffraction Study. *Chem. Mater.* **1999**, *11* (9), 2508–2514. https://doi.org/10.1021/cm991047q.

(29) Martínez-Iñesta, M. M.; Lobo, R. F. Investigation of the Negative Thermal Expansion Mechanism of Zeolite Chabazite Using the Pair Distribution Function Method. *J. Phys. Chem. B* **2005**, *109* (19), 9389–9396. https://doi.org/10.1021/jp044165w.

(30) Attfield, M. P. Strong Negative Thermal Expansion in Siliceous Faujasite. *Chem. Commun.* **1998**, No. 5, 601–602. https://doi.org/10.1039/a707141h.

(31) Leardini, L.; Quartieri, S.; Vezzalini, G.; Arletti, R. Thermal Behaviour of Siliceous Faujasite: Further Structural Interpretation of Negative Thermal Expansion. *Microporous and Mesoporous Materials* **2015**, *202*, 226–233. https://doi.org/10.1016/j.micromeso.2014.10.006.

(32) Lightfoot, P.; Woodcock, D. A.; Maple, M. J.; Villaescusa, L. A.; Wright, P. A. The Widespread Occurrence of Negative Thermalexpansion in Zeolites. *J. Mater. Chem.* **2001**, *11* (1), 212–216. https://doi.org/10.1039/B002950P.

(33) Carey, T.; Tang, C. C.; Hriljac, J. A.; Anderson, P. A. Chemical Control of Thermal Expansion in Cation-Exchanged Zeolite A. *Chem. Mater.* **2014**, *26* (4), 1561–1566. https://doi.org/10.1021/cm403312q.

(34) Braunbarth, C. M.; Behrens, P.; Felsche, J.; van de Goor, G. Phase Transitions and Thermal Behaviour of Silica Sodalites. *Solid State Ionics* **1997**, *101–103*, 1273–1277. https://doi.org/10.1016/S0167-2738(97)00401-3.

(35) Werner, S.; Barth, S.; Jordan, R.; Schulz, H. Single Crystal Study of Sodalite at High Pressure. *Zeitschrift für Kristallographie - Crystalline Materials* **2010**, *211* (3), 158–162. https://doi.org/10.1524/zkri.1996.211.3.158

(36) Compression Mechanism in Trioxane Silica Sodalite, [Si<sub>12</sub>O<sub>24</sub>]·2C<sub>3</sub>H<sub>6</sub>O<sub>3</sub>. *Zeitschrift für Kristallographie - Crystalline Materials* **1994**, *209* (6), 517–523. https://doi.org/10.1524/zkri.1994.209.6.517.

(37) Knorr, K.; Braunbarth, C. M.; Van De Goor, G.; Behrens, P.; Griewatsch, C.; Depmeier, W. High-Pressure Study on Dioxolane Silica Sodalite (C3H6O2)2[Si12O24]—Neutron and X-Ray Powder Diffraction Experiments. *Solid State Communications* **2000**, *113* (9), 503–507. https://doi.org/10.1016/S0038-1098(99)00516-5.

(38) Knorr, K.; Winkler, B.; Milman, V. Compression Mechanism of Cubic Silica Sodalite [Si<sub>12</sub>O<sub>24</sub>]: A First Principles Study of the Im ¯3m to I43m Phase Transition. *Zeitschrift für Kristallographie - Crystalline Materials* **2001**, *216* (9), 495–500. https://doi.org/10.1524/zkri.216.9.495.20347.

(39) Li, N.; Shi, Z.; Zhai, S.; Zhou, N.; Zhang, P.; Arramel, A.; Bennett, T. D.; Yue, Y. Unveiling the Amorphization of Sodalite Topology Zeolitic Imidazolate Frameworks and Zeolites by Pressure and Stress. *Journal of the American Ceramic Society* **2023**, *106* (10), 6117–6129. https://doi.org/10.1111/jace.19228.

(40) Havenga, E. A.; Huang, Y.; Secco, R. A. An Investigation of the Effect of High Pressure on the Structure of Siliceous Zeolite Y. *Materials Research Bulletin* **2003**, *38* (3), 381–387. https://doi.org/10.1016/S0025-5408(02)01061-9.

(41) Colligan, M.; Forster, P. M.; Cheetham, A. K.; Lee, Y.; Vogt, T.; Hriljac, J. A. Synchrotron X-Ray Powder Diffraction and Computational Investigation of Purely Siliceous Zeolite Y under Pressure. *J. Am. Chem. Soc.* **2004**, *126* (38), 12015–12022. https://doi.org/10.1021/ja048685g.

(42) Confalonieri, G.; Quartieri, S.; Vezzalini, G.; Tabacchi, G.; Fois, E.; Daou, T. J.; Arletti, R. Differential Penetration of Ethanol and Water in Si-Chabazite: High Pressure Dehydration of Azeotrope Solution.
*Microporous and Mesoporous Materials* **2019**, *284*, 161–169.
https://doi.org/10.1016/j.micromeso.2019.04.032.
TOC Graphic

<figure><img src="image_9.png" /></figure>