Thermodynamic rules for zeolite formation from machine learning based global optimization.

While the [TO4] tetrahedron packing rule leads to millions of likely zeolite structures, there are currently only 252 types of zeolite frameworks reported after decades of synthetic efforts. The subtle synthetic conditions, e.g. the structure-directing agents, pH and the feed ratio, were often blamed for the limited zeolite types due to the complex kinetics. Here by developing machine learning global optimization techniques, we are now able to establish the global potential energy surface of a typical zeolite system, Si x Al y P z O2H y-z with 12 T atoms (T: Si, Al and P) that is the general formula shared by CHA, ATS, ATO and ATV zeolite frameworks. After analyzing more than 106 minima data, we identify thermodynamic rules on energetics and local bonding patterns for stable zeolites. These rules provide general guidelines to classify zeolite types and correlate them with synthesis conditions. The machine learning based atomistic simulation thus paves a new way towards rational design and synthesis of stable zeolite frameworks with desirable compositions. This journal is © The Royal Society of Chemistry.

Introduction

Zeolites, a class of crystalline microporous materials with well-defined channels and cages packed with corner-sharing [TO4] tetrahedral units (T: Si, Al, P etc.), are widely used in chemical industries, e.g. as catalysts for converting petrochemicals.[1-7] While millions of hypothetical zeolite structures have been proposed by theory[8-11] and a number of synthetic approaches were innovated for zeolite synthesis, only 252 distinct zeolite structures (IZC-SC database)[12] have been successfully synthesized in the past century. The trial-and-error method prevails in zeolite synthesis via the testing of many conditions, including the feed Si : Al : P ratios, pH and structure-directing agents (SDAs), besides the crystallization temperature and the time conventional to chemical synthesis.[13-16] The physical laws governing the structures of zeolites remain largely a mystery to chemists.

Direct crystallization from an alkaline solution mixture containing silicon, aluminum and alkali metals is the most traditional way to synthesize zeolites, which only leads to aluminosilicate zeolites with a low Si : Al ratio (e.g. <10).[17] The kinetics was suggested to control the liquid-to-solid condensation due to the metastable nature of the product, where the strong alkali (OH−) acts as the key mineralizer to dissolve Si and Al ions.[4,18] However, the replacement of inorganic alkalis by tetramethylammonium found by Barrer et al. in 1961 [19] breaks the Si : Al ratio limitation to obtain even pure SiO2 zeolites (e.g. silicalite-1), which can also be produced using a recently reported solvent-free calcination synthetic route.[5,20] The aluminophosphate zeolite (AlPO) that entered the zeolite family in the 1980s[21] utilizes boehmite, phosphoric acid and organic amine as reagents. All these new findings question the essentiality of strong alkalis and the kinetics-controlled mechanism. Recent decades have witnessed an increasing variety of zeolites owing to the introduction of different SDAs. Some rules of thumb gleaned from synthesis are: (i) the SDA can define the pore size and thus is critical to the zeolite type and shape;[2,3,22] (ii) pH environment affects zeolite formation that aluminosilicate zeolites prefer alkaline condition but phosphate-containing zeolites like weak acidic or near neutral condition;[17,23] and (iii) the compositions of feed Si : Al : P ratios are limited to few choices, such as the silicoaluminophosphate (SAPO) zeolite that has nAl : nSi+P ∼ 1 : 1 and nP > nSi (nT: the number of element T atoms in the crystal). Many underlying questions on zeolites thus arise and three of them must rank top:

(Q1) Is zeolite formation thermodynamically controlled?

(Q2) Why are the Si : Al : P ratios not freely tunable in synthesis?

(Q3) How does pH control the zeolite formation?

To date, few theories are available to answer these questions. The most accepted rule is perhaps the exclusion of the Al–O–Al (known as Löwenstein's rule), P–O–P and Si–O–P patterns in zeolites, which were summarized from a reported zeolite structure.[17] Here by developing a machine learning based global optimization technique, we are able to, for the first time, establish the global potential energy surface (PES) of zeolites as represented by SiAlPO2H using a 12 T system as an example. This leads to quantitative solutions to the key puzzles associated with zeolite structures, compositions and energetics.

Methods

We have utilized a new technique, as implemented in our LASP code,[24] to solve the complex PES problem as encountered in zeolites. A five-element Si–Al–P–O–H global neural network (G-NN) potential is constructed via self-learning global PES data,[25,26] which can fast and accurately evaluate the PES and facilitate exploring the huge dimensionality of the zeolite global PES, in both the element and the configuration space, where the current quantum mechanics and empirical force field calculations fail either in speed or in accuracy. The global PES exploration is achieved via an enhanced stochastic surface walking (E-SSW) method designed for zeolite structures to cope with variable-size micropores in the framework. We note that traditional PES methods, such as molecular dynamics and evolution algorithms,[27] are generally frustrated in the global search for zeolite structures.

E-SSW method

The E-SSW method extends SSW global optimization by implementing explicitly virtual rigid bodies to enhance the PES search towards structures with open space. The SSW algorithm[28] is an unbiased global optimization method that can explore both minima and saddle points on the PES. SSW implements an automated climbing mechanism to manipulate a structural configuration moving smoothly from a local minimum to a high-energy configuration along one random mode direction. The method was initially developed for aperiodic systems,[29] such as molecules and clusters, and has been extended to periodic crystals.[30]

In the E-SSW method, the real PES is transformed after the addition of a tunable external potential U(ri), where i is the index of the rigid body, as shown in eqn (1). We utilize a power-type repulsive potential to create a repulsive sphere with a variable radius, see ESI Fig. S1.†[30]where Ereal is the energy from G-NN potential, rij is the distance between the rigid-body i and the real atom j in the system, and rs controls the repulsive range, varying from 0 to 6 Å in simulation.

E-SSW-NN simulation

Our approach to resolve complex zeolite (Si–Al–P–O–H elements) structures is based on the E-SSW-NN method where the E-SSW-NN combines E-SSW global PES exploration with fast G-NN potential calculations. The G-NN potential[31,32] is generated by iterative self-learning of the plane wave DFT global PES dataset generated from E-SSW search. The E-SSW-NN simulation to explore the PES can be divided into three steps: global dataset generation based on DFT calculations using selected structures from E-SSW simulation, G-NN potential fitting and E-SSW global optimization using G-NN potential. These steps are iteratively performed until the G-NN potential is transferable and robust enough to describe the global PES. The procedure is briefly summarized below.

At first, a global dataset is built iteratively during the iterative self-learning of G-NN potential. The initial data come from the density functional theory (DFT) based E-SSW simulation and the data in the subsequent cycles are from E-SSW-NN PES exploration. In order to maximally cover the likely configurations from all elements, extensive SSW simulations have been carried out for as many as possible structures, compositions and supercell sizes. Overall, these SSW simulations generate more than 107 structures on the PES. The final global dataset that is computed by high accuracy DFT calculation contains 27 135 structures, and is detailed in ESI Table S1.†

Then, the G-NN potential is generated using the method as introduced in our previous work.[31,33] To pursue a high accuracy for the PES, we have adopted a large set of power-type structure descriptors, which contains 343 descriptors for every element, including 129 2-body, 208 3-body, and 6 4-body descriptors, and compatibly, the network involves two-hidden layers (343-50-50-1 net), equivalent to 99 005 network parameters in total. The min–max scaling is utilized to normalization the training data sets. Hyperbolic tangent activation functions are used for the hidden layers, while a linear transformation is applied to the output layer of all networks. The limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method is used to minimize the loss function to match DFT energy, force and stress. The final energy and force criteria of the root mean square errors are around 2.73 meV per atom and 0.103 eV Å−1 respectively. The benchmark between G-NN and DFT results can be found in ESI Tables S2–S4,† which is accurate enough for searching for stable structure candidates.

Finally, E-SSW-NN simulation is performed over a wide range of compositions and structures, both for the global dataset generation and for the final production of the ternary phase diagram. For the E-SSW-NN simulation, each composition is simulated in the unit cells of 36–42 and explored to cover more than 20 000 minima on the PES by E-SSW. Thus, a large variety of structures ranging from crystalline to amorphous structures have been obtained from E-SSW-NN simulation.

DFT calculation

All DFT calculations are performed by using the plane wave VASP code, where electron–ion interaction is represented by the projector augmented wave pseudopotential.[34,35] The exchange–correlation functional utilized is the BEEF-vdw functional[36] because of its accurate estimation of adsorption energies[37] and its explicit inclusion of van der Waals interactions. The kinetic energy cutoff is set as 450 eV. The first Brillouin zone k-point sampling utilizes the Monkhorst–Pack scheme with an automated mesh determined by 25 times the reciprocal lattice vectors. The energy and force criteria for the convergence of the electron density and structure optimization are set at 10−6 eV and 0.02 eV Å−1, respectively. To determine the Gibbs formation energy (Gf) of SiAlPO2H and SiAlPO2Na, ab initio thermodynamics analyses have been performed with respect to quartz–SiO2, quartz–AlPO, α-Al2O3, NaOH solution and H2O at 200 °C and 15.53 bar (saturated vapor). The free energy of H2O is computed by the standard thermodynamics approach by incorporating zero-point energy and entropy contribution to the total energy (see the ESI†).

Results and discussion

Using the G-NN based E-SSW search, we have explored the global PES for every likely zeolite composition in the ternary (Si–Al–P) phase diagram from Al0.5P0.5O2 to SiO2 and to Si0.5Al0.5O2H0.5. For each composition, five different radii were examined, namely rs = 0, 3, 4, 5 and 6 Å, and each simulation with a fixed rs was set to explore 20 000 minima structures starting from known configurations (e.g. quartz and CHA-type zeolite). The minima were confirmed via full structure relaxation after the removal of the external potential and the low energy structures identified were finally verified by using DFT calculations.

Fig. 1a illustrates a representative global PES of AlPO (Al6P6O24), where the framework density (FD, the number of tetrahedral T atoms per 1000 Å3), the key structural feature of the zeolite, is plotted against the total energy of minima. Similar plots for other compositions can be found in ESI Fig. S2.† As denoted in the figure, the global PES can be divided roughly into three regions according to the FD value, i.e. densely packed, caged and layered regions from right to left. The densely packed region, as represented by the quartz structure (P3121) with FD = 24.6 which is the global minimum (GM), has FD values above 20. The zeolite belongs to the caged region with FD values of 12–20. Importantly, four known zeolites, i.e. with the topology of CHA, ATS, ATO and ATV, appear at the bottom of this region. They can be visualized as the three-dimensional (3D) crystalline assembly of the basic structural units, so-called d6r, ats, lau and afi, respectively (Fig. 1b). The energies of these four zeolites are 0.07–0.10 eV per TO2 (Al0.5P0.5O2) formula unit (f.u., i.e. per TO2) relative to the quartz GM, confirming the metastable nature of the caged structures. The layered region represents the most open structure with the FD values of 2–12, where a special 2D layer AlPO phase, denoted as the DL phase, is located at the bottom of this region, 0.17 eV per f.u. above the GM. The DL phase can be visualized as a 2D assembly of the d6r structural unit (Fig. 1b), similar to the CHA zeolite.

Fig. 1

Global PES exploration using E-SSW based on G-NN. (a) Global PES contour plot of Al6P6O24 minima. The x axis is the framework density (FD), the y axis is the total energy of minima with respect to the GM (quartz phase) and the color indicates the density of states (DOS). (b) Key low energy crystalline structures (black dots in (a)) together with their basic building blocks. Pink and blue balls represent the skeleton Al and P atoms, respectively. The yellow sphere indicates the zeolite void. (c) The variation of the GM structures identified from E-SSW global search in the presence of the rigid body at different rs values.

The global PES defines clearly the energy criterion for zeolite formation. Being the crystalline form of the caged structures, the energy of the zeolite has to be lower than that of the DL phase, but higher than that of the densely packed quartz phase, e.g. in between 0 and 0.17 eV per f.u., as indicated by the red lines in Fig. 1a. The universality of the energy criterion is supported by analyzing the zeolite bank. Among 252 known zeolite structures, there are 29 as-synthesized AlPO zeolites and their energies are in between 0.05 and 0.14 eV per f.u. (ESI Table S5†). The structure classification from the global PES and the energy criterion of the zeolite are not limited to AlPO. For SiO2 zeolites, similarly, the presence of the DL phase (0.18 eV per f.u. above quartz) dictates the upper energy bound: 60 as-synthesized SiO2 zeolites are indeed in the energy window from 0.04 to 0.18 eV per f.u. (ESI Table S6†). This finding suggests that the type of stable zeolite structure is strongly restrained by the presence of the corresponding DL phase, which provides the key cause for the fact that the types of known zeolites are far fewer than the theoretical prediction based on the TO4 packing rule.

Since the zeolite is located only in a small region of the global PES, being metastable in nature, it would be interesting to determine how they are formed under experimental conditions. Fig. 1c (data in ESI Table S7†) provides the clue, which illustrates the identified GMs of AlPO in the presence of an external rigid body (a rigid body per Al6P6O24). It shows that the increase of rs rapidly decreases the FD value: too large or too small rs fails to identify the zeolite, not surprisingly, leading to either a layered or a densely packed structure. The zeolite only turns out to be the GM under the suitable rs being applied, i.e. in between 3.5 and 5.5 Å. The four known zeolites, i.e. ATV-, ATO-, ATS- and CHA-types, do emerge as the GM at rs = 3.5, 4, 4.5 and 5–5.5 Å, respectively. Their cage size matches well with the rs value of the rigid body (see Fig. 1b). Indeed, in the experiment these zeolites are synthesized by using selected-size molecular SDAs: for example, the synthesis of ATO-type AlPO utilizes dipropylamine (8 Å diameter, see ESI Fig. S3†).[38]Fig. 1c thus supports that the zeolite may well be the thermodynamically favored product under the synthetic conditions. In particular, the choice of SDAs with suitable sizes can be the key way to condense [TO4] towards the zeolite, instead of the quartz or DL phase.

Now for each sensible rs, we can further examine the composition effect of the Si : Al : P ratio on the zeolite stability. By using the formation free energy (Gf) of the obtained GM structure, Fig. 2a illustrates the ternary phase diagram at rs = 5 Å for different SAPO compositions. Gf is the free energy of the zeolite (without the rigid body) relative to the free energies of quartz–SiO2, quartz–AlPO, α-Al2O3 and H2O at 200 °C and 15.53 bar (saturated vapor), corresponding to the neutral pH in normally used hydrothermal synthesis (in ESI Fig. S4† we discuss the effect of hydrothermal conditions on zeolite thermodynamics). Not surprisingly, due to the presence of the rigid body (rs = 5 Å) in our global PES search, all GMs in Fig. 2a are in the CHA-type framework, but they differ significantly in Gf: it has the minima nearby two vertexes, Al0.5P0.5O2 and SiO2 (∼0.11 eV per f.u.), but yields the maximum (∼0.30 eV per f.u.) at the left-bottom corner (Si0.5Al0.5O2H0.5). Most phases in the map are higher than their corresponding layer structure (0.17–0.18 eV per f.u., dotted lines), only the phases nearby two vertices (Al0.5P0.5O2 and SiO2), i.e. nAl : nSi+P ∼ 1 : 1 and nP > nSi (SAPO) and nSi : nAl > 5 : 1, survive in thermodynamics. As indicated by the brackets in the figure, these Si : Al : P ratios correspond to two CHA-type SAPO-34 and SSZ-13 materials, which are normally synthesized with the aid of SDAs (N,N,N-trimethyladamantammonium) in near neutral pH (pH < 8).[39,40] We note that SSZ-13 can be synthesized not only under neutral but also under alkaline conditions, and under both conditions the SDAs have to be supplied, suggesting the critical role of SDAs.[40-42] Overall, even in the presence of suitable rigid bodies, the Si : Al : P composition is still critical to the stability of zeolites, causing a non-freely tunable Si : Al : P ratio in experiments (ESI Table S8†).

Fig. 2

Thermodynamics of zeolite formation with different Si : Al : P ratios at acidic or neutral pH. (a) Gibbs formation free energy (Gf) contour plot in the ternary phase diagram using the GM identified by E-SSW with rs = 5 Å for each composition (black point). The Si : Al : P ratio is indicated for each composition; (b) correlation between Gf and the linear fitting of Gf (eqn (2)) for all GM data in (a).

Fig. 2b correlates the stability (Gf) with the structure using the GM data in Fig. 2a. By approximating Gf as a function of the proportions of TO4 (monomer, PT) and their linkages (T–O–T′, PTT′), see ESI Table S9,† we can obtain eqn (2) by linear fitting with R2 = 0.96. As Gf is positive in nature, it is no wonder that most terms, including monomers terms, Si–O–P, Al–O–Al, and Si–O–Al terms, have positive energy contributions. But it is important to reveal that Si–O–Si and Al–O–P terms yield negative contributions, suggesting that they are the major driving forces to stabilize the zeolite. The empirical rule in zeolite chemistry, namely, no Si–O–P and Al–O–Al patterns, is clearly manifested by their large positive prefactors, 0.35 and 0.16. In addition, our results identify the positive prefactor for the Si–O–Al term and thus explain the difficulty to incorporate Si element in AlPO and the special Si : Al : P ratio of nSi < nP in SAPO synthesized under neutral pH conditions. This finding on the bonding patterns is general and well confirmed in other zeolite frameworks (ESI Fig. S5†).

We would like to point out that low Si : Al ratio zeolites, while not favorable under acidic and neutral conditions, are known to form under strongly alkali conditions in the experiment (ESI Table S10†).[17] This is not surprising from thermodynamics as the replacement of H by an alkali metal, i.e. Na, changes the energy reference from H2O to NaOH solution. Fig. 3 illustrates the thermodynamic ternary phase diagram for different SAPO compositions under alkaline conditions (in the presence of Na ions for charge balance: the system is charge neutral) and all these calculations were performed by DFT. The structures of each composition are determined by Metropolis Monte Carlo sampling in the same CHA-type framework where the skeleton positions of T (T: Si, Al and P) and the likely positions of Na atoms (at the centers of 6- and 8-membered rings, see ESI Fig. S6† for details) are utilized as the pool of sites for random selection. The final Gibbs formation free energy Gf turns out to be completely different from that under acidic and neutral environments. A high Gf occurs at the compositions without Na, e.g. Si0.5Al0.25P0.25O2 (0.25 eV per f.u.) and the increase of the Na content can significantly stabilize the zeolite, leading to low Gf appearing at the left-bottom corner with a low Si : Al ratio. In particular, the Gf of Si7/12Al5/12O2Na5/12 is −0.01 eV per f.u., which is the most stable composition. As indicated by the brackets in the figure, these low Si : Al ratio zeolites happen to be chabazite materials, which were synthesized in strongly alkaline NaOH solution (pH > 13). Moreover, further considering the presence of water molecules to coordinate with Na ions in the zeolite framework, the zeolite containing Na can further be stabilized to reach even larger exothermic Gf (see, e.g. a hydrated SiAlO4Na zeolite in ESI Fig. S7†). This finding also proves that the zeolite formed under alkaline conditions is again the thermodynamically preferred product, which, unlike that under neural/acidic conditions, occurs even without the introduction of external SDA molecules. The Si–O–Al bonding pattern is thus only favored under alkaline conditions.

Fig. 3

Gibbs formation free energy (Gf) contour plot with different Si : Al : P ratios under alkaline conditions. The Si : Al : P ratio is indicated for each composition. All structures are in the CHA framework.

Both the energy bound and the pH-dependent bonding pattern rule indicate that thermodynamics dictates largely zeolite formation. This knowledge from the global PES leads to the classification of synthetic conditions for zeolites into three types:

• AlPO (and SAPO): SDAs + water + Al2O3 + H3PO4 + (SiO2);

• Pure silica (and aluminosilicates with high Si : Al ratios): SDAs + water + SiO2 + (Al2O3);

• Aluminosilicates with low Si : Al ratios: NaOH + water + SiO2 + Al2O3.

The first two types of zeolites rely on appropriate SDAs to enforce thermodynamics towards microporous structures, while a strong alkali leads to the third type zeolite, although SDAs may not be required.

With this knowledge, we have examined 252 known zeolite frameworks from the IZC-SC database[12] and we can quickly select 63 and 233 of the 252 frameworks that satisfy the energy criterion (<0.17 (0.18) eV per f.u.) for AlPO and SiO2 zeolites, respectively (as listed in ESI Tables S5 and S6†). The same approach is applicable to millions of conceived zeolite topologies in the DEEM PCOD database,[43] which shows that 14 900, only ∼4.5%, hypothetical zeolite structures satisfy the energy criterion. In the experiment, we note that there are 29 AlPO and 60 SiO2 zeolites synthesized, where only five (AST, ATS, AFI, CHA, and ANA) have both AlPO and SiO2 forms. Obviously, from our work there is ample room to synthesize new zeolites, particularly for AlPO (SAPO) and pure silica (high Si : Al ratio aluminosilicates).

Conclusions

By developing machine-learning based atomic simulation techniques, this work explores the global PES of zeolites which reveals the roles of SDAs, Si : Al : P ratios and pH values on the zeolite stability. We find that (i) zeolites, while being metastable in energy, become the thermodynamically stable products under synthetic conditions; (ii) the energy of zeolites is limited by an upper bound determined by the 2D layered phase, which is 0.17–0.18 eV per f.u.; and (iii) the preferred bonding patterns are pH sensitive: the Si–O–P, Al–O–Al and Si–O–Al bonding patterns are not favored under acidic and neutral pH environments, but the Si–O–Al pattern becomes desirable under basic conditions. With these rules, we answer the long-standing questions on zeolite formation, in particular on the Si : Al : P ratio, the choice of SDA molecules and the pH conditions in synthesis. Besides these thermodynamic rules, the advent of machine learning based atomic simulation enables the fast and reliable analyses of unknown zeolite structures and thus opens new avenues towards the rational design of zeolites.

Conflicts of interest

There are no conflicts to declare.