We have proposed a combination of density functional theory calculations and interatomic potential-based simulations to study the structural, electronic, and mechanical properties of pure-silica zeolite Linde Type A (LTA), as well as two titanium-doped compositions. The energetics of the titanium distribution within the zeolite framework suggest that the inclusion of a second titanium atom with configurations Ti(Si)0Ti, Ti(Si)1Ti, and Ti(Si)2Ti is more energetically favorable than the mono-substitution. Infra-red spectra have been simulated for the pure-silica LTA, the single titanium substitution, and the configurations Ti(Si)0Ti and Ti(Si)2Ti, comparing against experimental benchmarks where available. The energetics of the direct dissociation of water on these Lewis acid sites indicate that this process is only favored when two titanium atoms form a two-membered ring (2MR) sharing two hydroxy groups, Ti(OH)2Ti, which suggests that the presence of water may tune the distribution of titanium atoms within the framework of zeolite LTA. The electronic analysis indicates charge transfer from H2O to the Lewis acid site and hybridization of their electronic states.

Zeolite Linde Type A (LTA) was originally reported as an aluminosilicate material with formula Na12[(AlO2)12(SiO2)12]27H2O.1,2 The periodic building unit (PerBU) of LTA comprises three-dimensionally interconnected β-cages [Fig. 1(a)] that generate an α-cage at the center of a cubic unit cell with Pm3¯m symmetry [Fig. 1(b)]. An alternative PerBU to the β-cage consists of double 4-membered rings [D4R, Fig. 1(c)].3 The pore system of LTA is characterized by eight-membered ring (8MR) apertures with a diameter smaller than 7 Å connecting α-cages along the three dimensions [Fig. 1(e)]. The application of LTA as a molecular sieve was identified early on following its discovery by Breck and collaborators, who reported how the ion exchange of Na+ by other cations with charges 1+ and 2+ was able to tailor the adsorption capacity and selectivity for a wide array of small molecules.1 The relatively large size of the α-cages in LTA, combined with small apertures formed by 8MR, leads to attractive applications such as the formation and confinement of metal clusters, which are protected by the zeolite framework against poisoning, thereby allowing their selectivity towards the conversion of small molecules.4 

FIG. 1.

Representation of the (a) β-cage, (b) α-cage, and (c) double 4-membered ring (D4R). The cubic unit cell is represented by black lines; framework O atoms are deleted for clarity and Si atoms are represented by orange sticks. (d) Six-membered ring (6MR) and (e) eight-membered ring (8MR) together with relevant O–O distances obtained after geometry optimization (this work); framework O atoms are represented in red and Si atoms by orange sticks.

FIG. 1.

Representation of the (a) β-cage, (b) α-cage, and (c) double 4-membered ring (D4R). The cubic unit cell is represented by black lines; framework O atoms are deleted for clarity and Si atoms are represented by orange sticks. (d) Six-membered ring (6MR) and (e) eight-membered ring (8MR) together with relevant O–O distances obtained after geometry optimization (this work); framework O atoms are represented in red and Si atoms by orange sticks.

Close modal

Although zeolite LTA is mainly recognized by its high Al content (known as zeolite A), driving its applications as molecular sieve and sorbent,5 pure-silica LTA has been synthesized directly using self-assembled organic structure-directing agents (SDAs).6 Recently, high silica zeolite LTA, with Si/Al ratios of 11 and 16, has shown exceptional durability in the selective catalytic reduction of NOx with NH3 (NH3-SCR) under hydrothermal conditions at temperatures as high as 750 °C.7 Furthermore, aluminosilicate LTA with Si/Al ratios between 12 and 42 has shown promising results related to the selectivity of C5 and C6 products in the methanol-to-olefins (MTO) reaction.8 Boal et al. have also proposed the synthesis of titanosilicate LTA, suggesting that this material may be suitable for the oxidation of small molecules, taking advantage of the 8MR pore openings, the relatively large cavities, and the presence of intra-framework sites with Lewis acidity.8 Titanium-substituted zeolites have been used previously as catalysts for the epoxidation and oxidation of different organic molecules, such as phenol, ethylene, and bulkier olefins.9–12 

A number of theoretical reports have analyzed Ti substitution in different zeolites, such as silicalite, using cluster models and periodic boundary conditions (PBC).13–17 However, a systematic study of Ti substitution in zeolite LTA is still missing. The recent synthesis of high silica, Ti-substituted zeolite LTA8 provides the rationale for theoretical research to elucidate the local structure, electronic properties, and distribution of Ti(iv) within the framework of LTA, and this fundamental information may be of use in the design of future improved catalysts.

In the present work, we have examined single and double substitutions of Si by Ti in the T-sites of zeolite LTA. The combination of density functional theory (DFT) calculations with atomistic simulations allowed us to explore in detail the electronic structure of these systems whilst validating the use of low Si/Ti ratios against more realistic, lower Ti concentrations. We also considered the interaction of water with the Lewis acid site and its possible dissociation over the active site. In addition, the infrared (IR) spectra of relevant structures were simulated, which allowed us to propose a theoretical benchmark that may help us to interpret the experimental analysis of the zeolite structure.

We have performed the structural and electronic analysis of pure-silica and Ti-substituted zeolite LTA using DFT calculations as implemented in the Vienna Ab initio Simulation Package (VASP).18–21 In addition, we have included atomistic simulations that employ the Born model of ionic solids,22 to investigate a more realistic Si/Ti ratio to complement the DFT results.

The DFT calculations were carried out under periodic boundary conditions, employing a basis set of plane waves to explicitly treat the valence electrons, whilst the interaction with the inner core of the atom and nodal features were described by the projector-augmented-wave (PAW)23,24 method. The accuracy of the basis set was guaranteed by a kinetic energy cutoff of 550 eV. We have used the generalized gradient approximation (GGA), as proposed by Perdew, Burke, and Ernzerhof (PBE), to account for the exchange and correlation contributions to the electronic energy.25 The effect of long-range dispersion forces was considered using Grimme’s atomic pairwise method that includes three-body terms to avoid the overbinding of only two-body methods.26–28 Furthermore, we have chosen the Becke-Johnson damping function to avoid the double counting of the correlation effects at short distances.29 The Brillouin zone was sampled with a k-point mesh of 3 × 3 × 3, generated with the Monkhorst–Pack scheme.30 We have used the linear tetrahedron method, with Blöchl’s correction, to account for the occupation of the electronic bands and the integration over the k-point mesh during the structural optimization.31 This method was employed because the systems analyzed here are insulators, with all bands fully occupied or fully empty, thus avoiding the method’s drawback of not being variational with respect to partial occupancies. The thresholds for the electronic and ionic iterations were set at 10−6 eV and 0.01 eV/Å, respectively, ensuring proper convergence.

The GGA vibrational frequencies of the zeolite framework were obtained by the diagonalization of the Hessian matrix computed by density functional perturbation theory.32 The signal intensity of each normal mode was estimated using the following equation:

IIRα=iNτiνiα2,
(1)

where τi is the Born effective charge tensor of the ith atom, calculated by density functional perturbation theory and νiα is the atomic displacement in the αth eigenvector divided by the square root of its mass. The IR spectra were simulated by widening each signal with a Gaussian function.

The interatomic potential (IP) calculations were performed using the General Utility Lattice Program (GULP),33,34 where the energy of the system comprises Coulombic contributions, calculated according to the Ewald summation method,35 and short range repulsions and dispersion forces described by Buckingham potentials.36 For the potential model of the SiO2 units that form the zeolite framework, we have used the parametrization proposed by Sanders et al.,37 which was initially developed for α-quartz but has shown satisfactory transferability to zeolites.38–41 The presence of Ti within the zeolite framework was described by the potentials reported by Jentys and Catlow.14 The IP vibrational frequencies were computed after the diagonalization of the Hessian matrix, which was obtained by the finite difference method that keeps the system within the harmonic approximation. The signal intensity of each normal mode was evaluated by the following equation:

IIRα=iNqiνiα2,
(2)

where qi is the charge of the ith atom and νiα is the atomic displacement in the αth eigenvector divided by the square root of its mass.

Within the GGA, we optimized the LTA unit cell by using a set of fixed-volume calculations, where only the atomic positions were allowed to relax. The equilibrium volume was then evaluated by fitting the Birch-Murnaghan equation of state to the correlation between the lattice energy versus the lattice volume.42 This procedure also gave the isotropic bulk modulus as a fitted parameter. The optimization of the cell and the atomic positions with the IP method took place at constant pressure, thus permitting the simultaneous relaxation of the cell volume and the atomic coordinates in a single calculation. The unit cell of zeolite LTA has 24 T-sites, all of them equivalent by an operation of symmetry within the Pm3¯m space group.

Two different cells were used for the calculations. The first one, the 1 × 1 × 1 unit cell containing 24 tetrahedral sites (or equivalently: 24 TO2 units), was used for the GGA study of pure-silica and Ti-substituted zeolite LTA. For this cell size, the single Ti substitution gave a Si/Ti ratio of 23 whilst the double Ti substitution decreased it to 11. These values are at least 4 to 5 times smaller than those reported experimentally.8 However, a bigger cell would have been computationally extremely demanding within the electronic structure calculations. Therefore, we have complemented the GGA study with IP simulations, which allow the examination of larger systems using much less computational resource. Both the GGA and the IP methods have shown to be suitably accurate to analyze the thermodynamic and physical properties of zeolites.41,43–45

The second type of cell was formed by the 2 × 2 × 2 expansion of the 24-TO2 unit cell. As such, the number of TO2 units increased to 192, with a Si/Ti ratio of 191 and 95 for single and double substitutions, respectively, much closer to the experimental value of 100.8 All non-equivalent double substitutions of Ti in the zeolite framework were generated using the Site Occupancy Disorder (SOD) program, obtaining 12 structures for the 1 × 1 × 1 cell and 91 structures for the 2 × 2 × 2 supercell.46 In order to simplify the calculations, the presence of the SDA and its effect on the Ti substitution were not considered.

The structural visualizations presented in this work were generated by the Visualization for Electronic and Structural Analysis (VESTA 3).47 

The GGA elastic constants were obtained following the protocol of deformations proposed by Ravindran et al., where the volume is held constant during the orthorhombic and monoclinic shear distortions.48 Owing to the cubic symmetry of zeolite LTA, the elastic constants that we needed to evaluate reduced to C11, C12, and C44 only. The elastic constants were obtained by fitting a second order polynomial equation to the correlation between the lattice energy and the applied strain,48 

ED1=E0+V0τ1δ+V0C112δ2,
(3)
ED7=E0+V0(τ1τ2)δ+V0(C11C12)δ2,
(4)
ED4=E0+2V0τ4δ+2V0C44δ2,
(5)

where EDn is the set of lattice energies obtained after applying the deformation Dn suggested by Ravindran and collaborators,48 E0 is the energy at the equilibrium volume V0, τi is an element in the stress tensor, and δ is the applied strain (between −3% and 3% to minimize the contribution of higher order terms in the Taylor expansion).

The bulk modulus (B) and shear modulus (G) were obtained following the Voigt approximation, where the external strain is equated to the uniform strain in a polycrystalline ensemble.49 The values of B and G correspond to the resistance of the material to fracture and plastic deformation, respectively. The Voigt equations use the elastic constants for the calculation of B and G,

B=C11+2C123,
(6)
G=C11C12+3C445.
(7)

In addition, Young’s modulus (Y) and Poisson’s ratio (σ), which give information regarding the rigidity of the crystal, and the shear anisotropic factor (A), which measures the level of anisotropy of the material, were calculated from the following equations:

Y=9BG3B+G,
(8)
σ=3B2G6B+2G,
(9)
A=2C44C11C12.
(10)

In the case of the IP method, the elastic parameters were directly provided by the GULP following Voigt’s approximation.

After structural relaxation of pure-silica LTA within the GGA, we obtained a cubic unit cell parameter of 11.950 Å, which was only 0.8% longer than the reported value by Boal et al. (see Table I).8 The IP method decreased the difference even further, finishing with a value after optimization of 11.847 Å, only slightly shorter than the experimental parameter. The small discrepancy between the cell parameters evaluated by the GGA and IP was because of longer Si–O distances for the GGA compared to the IP method, by an average of 0.022 Å. A similar outcome has been reported for the optimization of the unit cell of zeolite MFI, where the GGA average value was 0.112 Å longer than the IP prediction.41 

TABLE I.

Optimized cell parameter a (Å), average and standard deviation of the SiO bond distance (Å), and OSiO and SiOSi angles (°). Elastic constants C11, C12, and C44 (GPa), bulk modulus B (GPa), shear modulus G (GPa), B/G ratio, Young’s modulus Y (GPa), Poisson’s ratio σ, and shear anisotropic factor A, obtained with the GGA and the IP method.

GGAIPExpt.
a 11.950 11.847 11.857a 
SiO 1.621 ± 0.007 1.599 ± 0.007 1.602 ± 0.006a 
OSiO 109.5 ± 0.9 109.5 ± 0.7 109.5 ± 0.4a 
SiOSi 148 ± 4 150 ± 4 150 ± 4a 
C11 51 102 … 
C12 27 60 … 
C44 19 25 … 
BBMb 42 … 22.1c 
Bd 35 74 22.1c 
G 16 23 … 
B/G 2.2 3.2 … 
Y 42 58 … 
σ 0.30 0.37 … 
A 1.6 1.2 … 
GGAIPExpt.
a 11.950 11.847 11.857a 
SiO 1.621 ± 0.007 1.599 ± 0.007 1.602 ± 0.006a 
OSiO 109.5 ± 0.9 109.5 ± 0.7 109.5 ± 0.4a 
SiOSi 148 ± 4 150 ± 4 150 ± 4a 
C11 51 102 … 
C12 27 60 … 
C44 19 25 … 
BBMb 42 … 22.1c 
Bd 35 74 22.1c 
G 16 23 … 
B/G 2.2 3.2 … 
Y 42 58 … 
σ 0.30 0.37 … 
A 1.6 1.2 … 
a

From the work of Boal et al.8 

b

Value obtained by fitting to the Birch-Murnaghan equation of state.42 

c

Value for zeolite A, formula: [Na12[(AlO2)12(SiO2)12]27H2O], from the work of Arletti et al.50 

d

Value obtained using the Voigt’s approximation.49 

We calculated a bulk modulus B of 42 GPa by fitting the GGA-calculated energies to the Birch-Murnaghan equation,42 decreasing to 35 GPa when employing the GGA-calculated elastic constants within the Voigt equation (see Table I).49 The IP method led to a value for B of 74 GPa. These results followed the same trend of zeolite MFI, where a B of 18 GPa was evaluated by GGA, rising to 54 GPa with IP,41 with the GGA value being much closer to the experimental 18.2 GPa.51 The validation of our results is more difficult in the case of pure-silica zeolite LTA because most of the experimental measurements provided in the literature come from zeolite-A, which has the same framework of LTA but with a Si/Al ratio of 1 and loaded with 6–12 extra-framework cations per unit cell to compensate for the negative charge of the framework. For instance, Arletti et al. have reported a value of 22.1 GPa for zeolite-A,50 whilst Niwa et al. have measured values between 25.3 and 34.4 GPa depending on the pressure-transmitting medium.52 Our calculated B is larger than the experimentally measured values, but this may be attributed to the modifications in bonding strengths when the Si/Al ratio changes from 1 to infinity, which makes the structure more resistant to deformations; an experimental benchmark using pure-silica LTA is necessary to validate our results.

The values of the calculated shear moduli G, presented in Table I, of 16 and 23 GPa for GGA and IP, respectively, mean that zeolite LTA has a relatively low resistance to plastic deformation.53,54 Huang and Havenga have actually reported the reversibility of the pressure-induced amorphization of zeolite-A in the range of 15–25 kbar, although they showed that the presence of extra-framework cations is essential for recovering the crystalline structure after the deformation.55G calculated in this study considered a pure-silica LTA without extra-framework cations, which should therefore make its G value larger than the one for zeolite-A.

At the same time, the ratio B/G measures the level of ductility (if high) or brittleness (if low) of a material, with the turning point between ductility and brittleness at approximately 1.75.53 We obtained values of 2.2 and 3.2 for B/G, using the GGA and IP method, respectively. Both types of calculations depict zeolite LTA as ductile, although the larger value of IP should be due to the overestimation of B.

The Poisson ratio σ is a measure of the ionic or covalent character of the bonds within the material. Values smaller than 0.1 are attributed to covalent bonding, whilst a ratio above 0.25 suggests ionic interactions.56 The GGA and IP led to similar values of σ of 0.30 and 0.37, respectively, with both methods giving values above 0.25. This result also validates the use of the Born model of ionic solids to describe the largely ionic zeolite structure.22 

We have also calculated the shear anisotropic factor A that evaluates the level of anisotropy of the crystal for values smaller or larger than unity. In the present study, GGA had an A of 1.6, reducing to 1.2 in the case of IP. Accordingly, pure-silica LTA should show a significant degree of anisotropy according to the GGA. However, a previous report has shown that A may be overestimated by the GGA methodology, where a theoretical value of 1.34 is predicted for magnetite (Fe3O4) against an experimental measurement of 1.13.57 

We have studied the stability and distribution of single and double substitutions of Si by Ti using Boltzmann statistics, a method that has been used in the past to analyze single and double substitutions of aluminum in zeolite MFI.58 The 1 × 1 × 1 unit cell was used for the GGA calculations, whilst an expanded 2 × 2 × 2 supercell was used with the IP method in order to complement the GGA results. Figure 2 shows the Boltzmann distributions at 125 °C, which is the synthesis temperature used by Boal et al. to obtain the titanosilicate.8 

FIG. 2.

Boltzmann distributions at 125 °C for double Ti-substituted zeolite LTA as a function of the separation between the two Ti atoms within the framework. Segments of the optimized structures are shown for the peaks with significant probability. The 0_b configuration, although with a contribution of approximately 0.0001 to the distribution, is also shown due to its later-discussed relevance in the dissociative adsorption of water. (a) The GGA using a 24-TO2 unit cell, (b) the IP method using a 192-TO2 supercell and considering only double substitutions within the same unit cell out of eight available, and (c) the IP method using a 192-TO2 supercell and considering the full set of double substitution along the supercell. The number within the label of each configuration specifies the number of Si atoms that bridge the two Ti atoms; the letter is only to organize the configurations within each group.

FIG. 2.

Boltzmann distributions at 125 °C for double Ti-substituted zeolite LTA as a function of the separation between the two Ti atoms within the framework. Segments of the optimized structures are shown for the peaks with significant probability. The 0_b configuration, although with a contribution of approximately 0.0001 to the distribution, is also shown due to its later-discussed relevance in the dissociative adsorption of water. (a) The GGA using a 24-TO2 unit cell, (b) the IP method using a 192-TO2 supercell and considering only double substitutions within the same unit cell out of eight available, and (c) the IP method using a 192-TO2 supercell and considering the full set of double substitution along the supercell. The number within the label of each configuration specifies the number of Si atoms that bridge the two Ti atoms; the letter is only to organize the configurations within each group.

Close modal

The GGA calculations identified five double Ti-substituted structures with significant probabilities out of the twelve non-equivalent substitutions in the 1 × 1 × 1 unit cell. Owing to its small size, only the arrangements TiTi, TiSiTi, and TiSiSiTi were obtained for this cell, where configuration 2_a in Fig. 2(a) had the highest prevalence (the number in the label specifies the number of bridging Si atoms placed between the two Ti atoms). In 2_a, the two Ti atoms were evenly placed at opposite sides of the six-membered ring (6MR) of zeolite LTA [see Fig. 2(a)], separated by a TiTi distance of 6.2 Å. This ordering is similar to the one reported by Bare et al., who used EXAFS to study the substitution of Sn in β-zeolite and highlighted the non-random substitution of Sn in the zeolite framework.59 We did not consider the presence of the SDA in this work; therefore, the stability of the 6MR substitution may be related not only to the employed SDA but also to the capacity of this configuration to minimize the framework deformations. Other important structures placed both Ti atoms in the same D4R, labeled as 0_a and 1_a in Fig. 2(a), with TiTi distances of 3.2 and 4.6 Å, respectively, and in two different D4Rs, labeled as 1_b and 2_b, with TiTi distances of 5.3 and 7.2 Å, respectively.

We generated two different Boltzmann distributions for the IP method. The first one, shown in Fig. 2(b), only considered the substitutions within the same unit cell, despite using a 2 × 2 × 2 supercell, to allow a direct comparison with the GGA results. The second distribution used the full spectrum of double substitutions (91 in total) and is shown in Fig. 2(c). A close agreement between the IP and GGA profiles can be noted from the comparison of Figs. 2(a) and 2(b), regardless of the difference in the level of theory and cell size. Again, 2_a was by far the most stable structure, although the intensity of 0_a decreased to almost zero and an additional configuration labeled as 2_c arose [see Fig. 2(b)]. The full spectrum of double substitution shown in Fig. 2(c) confirmed once more that 2_a is the most stable arrangement when two Ti are placed in the zeolite framework. In Fig. 2(c), 2_a was compared against structures with TiTi separations of more than 15 Å, which can be considered as two fully isolated species.

In the present work, we are not primarily concerned about the absolute stability of each type of substitution, but instead we have compared the relative stability of the double substitution against the single substitution by using the following equation:

Δn=ETinLTA+n1ELTAnETi1LTA,
(11)

where ETinLTA is the energy of the n Ti-substituted zeolite, ELTA is the energy of the pure-silica zeolite, and ETi1LTA is the energy of the single Ti-substituted zeolite. The derivation of Eq. (11) from the absolute stabilities is shown in the supplementary material. This equation can be interpreted as the energetics of a process that moves n Ti atoms from n single substituted unit cells to a sole unit cell, leaving the other (n − 1) as pure-silica frameworks. By using Eq. (11) with the GGA energies and n = 2, the calculated values of Δ2 were −11, −9, −8, −18, and −7 kJ/mol for configurations 0_a, 1_a, 1_b, 2_a, and 2_b, respectively. The other seven substitutions yielded positive values. This means that the double Ti substitutions represented in Fig. 2(a) are more favorable than the single substitution per 1 × 1 × 1 unit cell. The IP method showed equivalent outcomes. These results are different from a similar analysis by Yang et al., who studied Ti substitution in β-zeolite by DFT and found that the single substitution is more stable than the double substitution of Ti in the 6MR.60 Nevertheless, this apparent discrepancy may be related to the topology of the 6MR: the ring is planar in LTA whilst it is bent in β-zeolite, which translates into different capacities to support framework distortions. We considered pure-silica LTA, single Ti-substituted LTA, and double Ti-substituted LTA with configurations 0_a and 2_a, as fair representations of the structural diversity; hence, we have only used these structures to continue the analysis in this work.

The methodology developed from Eq. (11) can be used to determine the limit of Ti aggregation in zeolite LTA by adding increasing amounts of Ti atoms in the structure and evaluating the resulting stability of the n-substituted framework against the (n − 1) substitution. However, this analysis is beyond the scope of the present work and will be discussed in future publications.

The SiO bond length and the OSiO and SiOSi angles of Ti-substituted zeolite LTA did not manifest significant modifications when compared with the pure-silica structure (see Table S1 of supplementary material). The GGA- and IP-calculated lengths of the TiO bonds were longer than the SiO bonds by approximately 0.2 Å. The GGA-calculated TiOTi angle of the 0_a structure was 14° to 24° smaller than the average SiOSi and SiOTi angles. This reduction was due to the simultaneous increase in length of the two TiO bonds, which forced the bridging O atom away from the imaginary line that connects both Ti atoms, in consequence reducing the TiOTi angle. In contrast, the IP method did not show marked differences among the SiOSi, SiOTi, and TiOTi angles (see Table S1 of the supplementary material). At the same time, upon Ti substitution, the difference between the lowest unoccupied band (LUB) and the Fermi energy decreased (see Table S2 of the supplementary material). In addition, the electronic charge density associated with the LUB was projected out of the D4Rs into the cavities following the Ti substitution, as shown in Fig. S1 of the supplementary material. These results point to the Lewis acid character that is induced in the structure by the Ti atoms.

1. Vibrational analysis

Infrared (IR) spectroscopy is a powerful experimental technique that provides structural information specific to each material. The theoretical simulation of the IR spectrum enables the correlation of variations in the spectrum profile with explicit modifications in the material structure without undesired interferences, which then may complement and help us to understand the experimental measurements.

We have simulated the theoretical IR spectrum of pure-silica LTA and validated it against the equivalent experimental report by Huang and Jiang.61 Figure 3 shows the simulated IR spectra using the GGA and the IP method, with the vibrational frequencies of the main bands listed in Table II. The theoretical IR spectra of pure-silica LTA present four main bands for both the IP method and the GGA, resembling the experimental spectrum within the limits of the theoretical approximations.61 The most intense signal, located at 928 cm−1 for IP and 1073 cm−1 for GGA, corresponds to the asymmetric SiO stretching, which has an experimental value of 1096 cm−1,61 i.e., very close to the GGA value. The other three bands are located within the range 650–400 cm−1, reproducing the experimental signals with deviations of −31 to +6 cm−1 for both the IP method and the GGA. Within this range, the experimental bands at 621 and 477 cm−1 are related to the D4R, whilst the signal at 407 cm−1 is associated with the 6MR.61 A close inspection of the calculated vibrations in Table II for pure-silica LTA indicates a regular underestimation of the frequencies compared with the experiment, a trend that has been noted previously by Huang and Jiang.61 

FIG. 3.

Simulated infrared spectra of relevant structures (from bottom to top): pure-silica LTA (pure-Si), LTA with a single Ti substitution (1Ti), configuration 0_a, and configuration 2_a. The labeling of each spectrum contains within parenthesis the level of theory used in the simulation, i.e., the GGA or the IP method.

FIG. 3.

Simulated infrared spectra of relevant structures (from bottom to top): pure-silica LTA (pure-Si), LTA with a single Ti substitution (1Ti), configuration 0_a, and configuration 2_a. The labeling of each spectrum contains within parenthesis the level of theory used in the simulation, i.e., the GGA or the IP method.

Close modal
TABLE II.

Vibrational frequencies (cm−1) of the simulated IR spectra in Fig. 3.

IPaGGAExpt.61 GGA
Pure-Si 1Ti 0_a 2_a 
     1168 w 
     1143 w 
928 s 1073 s 1096 s 1067 s 1063 s 1062 s 
   978 w 972 m 973 m 
     905 w 
    755 w  
627 w 597 w 621 w 593 w 591 w 585 w 
464 w 449 w 477 m 450 w 452 w 446 w 
384 w 376 w 407 m 375 w 375 w 374 w 
IPaGGAExpt.61 GGA
Pure-Si 1Ti 0_a 2_a 
     1168 w 
     1143 w 
928 s 1073 s 1096 s 1067 s 1063 s 1062 s 
   978 w 972 m 973 m 
     905 w 
    755 w  
627 w 597 w 621 w 593 w 591 w 585 w 
464 w 449 w 477 m 450 w 452 w 446 w 
384 w 376 w 407 m 375 w 375 w 374 w 
a

Weak (w), medium (m), and strong (s) intensity.

The GGA was better than the IP method in reproducing the experimental spectrum of pure-silica LTA; therefore, we have only analyzed single and double Ti substitutions by GGA-simulated IR spectra. In addition, the IP method manifested a lack of sensibility towards the Ti substitution, showing practically no variations when one or two Ti atoms were included within the zeolite framework.

The single Ti substitution causes a red-shift in the band of the asymmetric SiO stretching by 6 cm−1, which coincides with the softening of this band observed in the experiment when Si atoms are replaced by Al, shifting by almost 90 cm−1 from pure-silica LTA to zeolite A.61 The stretching for TiO appears at 978 cm−1 with lesser shoulders at 951 and 920 cm−1. This signal is shifted with respect to the asymmetric SiO stretch by 89 cm−1, corresponding to the weaker and longer bond length of TiO when compared to SiO. The rest of the spectrum does not show significant variations after the single substitution.

The double substitution represented by structure 0_a in Fig. 2 enhances the intensity of the asymmetric TiO stretching as expected and produces an additional weak band at 755 cm−1. The signal at 755 cm−1 is related to the stretching of the TiO bond associated with the O atom bridging the two Ti atoms. As described above, the angle TiOTi in 0_a decreases when compared with the pure-silica value, moving the bridging O atom further away from the axis linking the two Ti atoms, which, combined with longer TiO bonds, makes the TiO(Ti) stretching 217 cm−1 weaker than TiO(Si). The band at 755 cm−1 disappears again in the spectrum of structure 2_a, where the two Ti atoms are separated by two Si atoms, in agreement with our interpretation.

The spectrum of 2_a presents weak peaks above the asymmetric Si–O band at 1168 and 1143 cm−1. These signals are related to SiO(Si) stretching modes associated with the D4R where the Ti atoms are substituted. Equivalent signals tend to appear in the spectrum of the single Ti-substituted LTA, but their intensities are very weak, almost negligible (see Fig. 3). However, the spectrum of structure 0_a lacks those bands, which may help us to differentiate structures Ti(Si)0Ti from those with a configuration Ti(Si)n>0Ti. Additionally, narrower TiOSi angles cause softening of the TiO(Si) stretching by 68 cm−1, producing a very weak band at 905 cm−1 in the spectrum of 2_a. For instance, while the TiOSi angle is approximately 165° for the band at 973 cm−1, it is reduced to 138° for the signal at 905 cm−1.

The Lewis acidity can be probed by studying the adsorption strength of small molecules at the active sites.13,60,62 We have chosen water as a probe molecule owing to its presence during the synthesis of the zeolite,8 which may also influence the distribution of Ti atoms along the framework positions. Pure-silica LTA, single Ti-substituted LTA, and configurations 0_a and 2_a were used to analyze the adsorption of water. Table III lists the highest adsorption energies of a water molecule in the α-cage and the β-cage (see Fig. 4), including the Van der Waals (VdW) contribution to these energies, for pure-silica LTA and single and double Ti substitutions.

TABLE III.

Adsorption energy (Eads.) of a water molecule in pure-silica LTA (pure-Si), single Ti-substituted LTA (1Ti), configuration 0_a, and configuration 2_a (kJ/mol); contribution of dispersion corrections (VdW) to the adsorption energy (%).

α-cageβ-cage
Eads.VdWEads.VdW
pure-Si −21 66 −32 68 
1Ti −43 41 −52 46 
0_a −43 42 −48 48 
2_a −37 46 −51 46 
α-cageβ-cage
Eads.VdWEads.VdW
pure-Si −21 66 −32 68 
1Ti −43 41 −52 46 
0_a −43 42 −48 48 
2_a −37 46 −51 46 
FIG. 4.

Representation of the adsorption of water in (a) the α-cage and (b) the β-cage, taking as an example the optimized structures of single Ti-substituted LTA. The framework O atoms are deleted for clarity; the water O atom is represented in red, H atoms in white, Ti atoms in light blue, and the silicon atoms by orange sticks.

FIG. 4.

Representation of the adsorption of water in (a) the α-cage and (b) the β-cage, taking as an example the optimized structures of single Ti-substituted LTA. The framework O atoms are deleted for clarity; the water O atom is represented in red, H atoms in white, Ti atoms in light blue, and the silicon atoms by orange sticks.

Close modal

The adsorption of water in pure-silica LTA released energies of 21 and 32 kJ/mol for loadings in the α-cage and the β-cage, respectively, where the VdW contributions accounted for 66% and 68%, respectively. The water molecule did not have a preferred adsorption site, as expected from its interaction with a completely siliceous framework.60 However, the reduced space of the β-cage, compared to the α-cage, produced stronger VdW and H-bond interactions between the water molecule and the surrounding framework atoms.

As expected, the Ti substitution increased the interaction strength of water with the zeolite. The exothermic energies for adsorption in the α-cage were estimated at 37–43 kJ/mol for single and double Ti substitutions, increasing to 48–52 kJ/mol for loadings in the β-cage (see Table III). The VdW contribution to the adsorption energy decreased from 66% to 68% for pure-silica LTA to 41%–48% for the Ti-substituted structures, which is anticipated for water-Ti interactions with more of a chemisorption character. The Ti–Owater distances ranged from 2.363 to 2.424 Å, with the β-cage producing smaller distances than the α-cage by 0.003–0.058 Å.

We have also analyzed the possibility of water dissociation on the Lewis site by breaking one of the OH bonds of water and transferring the proton to the TiO4 tetrahedron. However, in most cases, either the proton returned to the OH group binding the Ti atom, re-forming a water molecule, or the structure was more than 60 kJ/mol less stable than the non-dissociated system. Water dissociation was more likely when a two-membered ring (2MR) was formed, where both OH groups bound simultaneously the Ti atom and its nearest-neighbor Si atom from a different D4R, as shown in Fig. 5(b). We did not observe the formation of 2MRs when the Ti and the Si atoms belonged to the same D4R. We attribute this outcome to the inability of a single D4R to accommodate the distortions provoked by the formation of the 2MR. After observing the importance of the 2MR formation, we also optimized the pure-silica 2MR [Fig. 5(a)] and the double Ti-substituted 2MR [Fig. 5(c)] in order to analyze the dependence of the 2MR stability on the extent of the Ti substitution (for relative energies of the process, see Table IV).

FIG. 5.

Close-ups of the optimized structures of two-membered rings (2MR) with different degrees of Ti substitution. (a) 2MR for pure-silica zeolite, (b) 2MR for single Ti-substituted zeolite, and (c) 2MR for double Ti-substituted zeolite. (d) Location of the 2MR within the framework using the optimized structure of pure-silica zeolite as an example, the close-up specifies the labels used to distinguish the different atoms of the 2MR. The framework O atoms are deleted for clarity; the O atoms of the 2MR are represented in red, H atoms in white, Ti atoms in light blue, and the Si atoms by orange sticks.

FIG. 5.

Close-ups of the optimized structures of two-membered rings (2MR) with different degrees of Ti substitution. (a) 2MR for pure-silica zeolite, (b) 2MR for single Ti-substituted zeolite, and (c) 2MR for double Ti-substituted zeolite. (d) Location of the 2MR within the framework using the optimized structure of pure-silica zeolite as an example, the close-up specifies the labels used to distinguish the different atoms of the 2MR. The framework O atoms are deleted for clarity; the O atoms of the 2MR are represented in red, H atoms in white, Ti atoms in light blue, and the Si atoms by orange sticks.

Close modal
TABLE IV.

T–O distance (Å) (T = Si and Ti), OO distance (Å), O–H distance (Å), and relative energy against non-dissociated water (kJ/mol) of different 2MRs.

2Si-2MRSi,Ti-2MR2Ti-2MR
SiO1a 1.861, 1.866 1.859 … 
SiO2 1.799, 1.795 1.736  
TiO1 … 2.010 1.991, 2.036 
TiO2 … 2.043 2.007, 1.950 
O1O2 2.147 2.254 2.422 
H1O1 0.971 0.971 0.971 
H2O2 0.973 0.971 0.970 
Rel. en.b +108 (−31) +50 (−52) +15 (−71) 
2Si-2MRSi,Ti-2MR2Ti-2MR
SiO1a 1.861, 1.866 1.859 … 
SiO2 1.799, 1.795 1.736  
TiO1 … 2.010 1.991, 2.036 
TiO2 … 2.043 2.007, 1.950 
O1O2 2.147 2.254 2.422 
H1O1 0.971 0.971 0.971 
H2O2 0.973 0.971 0.970 
Rel. en.b +108 (−31) +50 (−52) +15 (−71) 
a

When two values are given for the same entry of TO distance, the first number refers to the distance T1O and the second refers to T2O, following the labeling in Fig. 5.

b

The adsorption energy of non-dissociated water is indicated within parenthesis.

The pure-silica 2MR (labeled as 2Si-2MR) has been experimentally detected by Hunger et al. using magic-angle spinning nuclear magnetic resonance (MAS-NMR).63 The authors observed that when tetrapropylammonium (TPA), acting as the SDA, was used to synthesize zeolite ZSM-5, internal silanols were formed as defects in a concentration of up to 8%. Sokol et al., employing DFT under periodic boundary conditions, found a relatively low energy requirement for the hydrolysis of the Si-O-Si linkages of sodalite with the subsequent formation of 2Si-2MR.64 They suggested that the 2MR may undergo dehydrogenation, creating triplet peroxide species, which transforms into a singlet peroxy bridge with an OO distance of approximately 1.58 Å. The authors further proposed that the migration of these peroxy bridges could be the reason for the increase in the Brønsted and Lewis acidity, which makes the redox chemistry of these materials even more complex.64 

In this study, the dissociation of adsorbed water that led to the formation of 2Si-2MR increased the energy of the system by +108 kJ/mol (see Table IV). The SiO bond lengths of 2Si-2MR remained between 1.795 and 1.866 Å, with an O1O2 distance of 2.147 Å, in agreement with previous calculations.64 

The substitution of one Si by Ti within the 2MR (labeled as Si,Ti-2MR) increased the stability of the structure, which was apparent from the reduction of the energy difference relative to non-dissociated water. In this case, the Si,Ti-2MR was +50 kJ/mol less stable than the system with non-dissociated water, representing a marked decrease in energy difference when compared with the value of +108 kJ/mol found for 2Si-2MR (see Table IV). Nevertheless, Si,Ti-2MR was still hardly stable, with a mere relative energy of −2 kJ/mol against the standard of water in the gas phase plus the bare zeolite without adsorbate. In this regard, it has been reported that Si,Ti-2MR is the least stable configuration when compared with models that reflect different levels of hydration and inversion of the tetrahedra SiO4 and TiO4 such as tripodal, tetrapodal, and bipodal configurations.65 

When we considered the 2MR formed by two Ti atoms (labeled as 2Ti-2MR), the stability of 2Ti-2MR increased significantly. As shown in Table IV, 2Ti-2MR was only 15 kJ/mol less stable than the system with non-dissociated water, in contrast to the values of +108 and +50 kJ/mol for 2Si-2MR and Si,Ti-2MR, respectively. In addition, the adsorption energy of dissociated water forming the 2Ti-2MR structure was −56 kJ/mol, while the adsorption of non-dissociated water, which served as the precursor for the formation of 2Ti-2MR, released 71 kJ/mol. These adsorption energies were calculated as relative energies against water in the gas phase and the corresponding non-hydrated Ti configuration in 2Ti-2MR (labeled as 0_b in Fig. 2). Configuration 0_b was 19 kJ/mol less stable than 0_a; however, upon water adsorption in the β-cage, 0_b was 4 kJ/mol more stable than 0_a. This difference increases if the 2Ti-2MR formation is considered: 2Ti-2MR(0_b) was 32 kJ/mol more stable than 2Ti-2MR(0_a), which is due to the higher structural flexibility provided by locating the two Ti atoms in different D4Rs (0_b case) than in the same D4R (0_a case). Therefore, we reason that the presence of water during the synthesis process can stabilize 0_b, increasing its probability of occurrence against configurations that seem more stable in non-hydrated conditions, such as 0_a and 2_a. As such, the presence of configuration 0_b in zeolite LTA could have important consequences for the performance of this material as a catalyst for redox reactions: two Ti atoms, at close distance from each other and with the capacity to interact simultaneously with the same molecule, may adsorb and stabilize more effectively peroxide species than single Ti-substituted frameworks.17 

1. Electronic structure analysis

We used the single Ti-substituted LTA and configuration 0_b, with water adsorbed in their respective β-cages [see Figs. 6(a) and 6(b)], as representative models to examine the electronic structure of hydrated Ti-LTA. The system water/0_b represented in Fig. 6(b) was the precursor for the formation of 2Ti-2MR(0_b), with an adsorption energy of −71 kJ/mol.

FIG. 6.

Adsorption of water in the β-cages of (a) single Ti-substituted LTA and (b) configuration 0_b. The framework O atoms are deleted for clarity (except for the bridging O atom of 0_b). Electronic charge density difference (Δρ) isosurfaces (0.0015 bohr−3) of the systems (c) water/1Ti and (d) water/0_b (positive density in yellow and negative density in cyan); framework O atoms directly binding the Ti atoms are shown in red; the Bader charge of the water molecule is indicated below each image. Charge density isosurfaces (0.0015 bohr−3) corresponding to the LUB of the systems (e) water/1Ti and (f) water/0_b. The O atoms are represented in red, H atoms in white, Ti atoms in light blue, and Si atoms by orange sticks.

FIG. 6.

Adsorption of water in the β-cages of (a) single Ti-substituted LTA and (b) configuration 0_b. The framework O atoms are deleted for clarity (except for the bridging O atom of 0_b). Electronic charge density difference (Δρ) isosurfaces (0.0015 bohr−3) of the systems (c) water/1Ti and (d) water/0_b (positive density in yellow and negative density in cyan); framework O atoms directly binding the Ti atoms are shown in red; the Bader charge of the water molecule is indicated below each image. Charge density isosurfaces (0.0015 bohr−3) corresponding to the LUB of the systems (e) water/1Ti and (f) water/0_b. The O atoms are represented in red, H atoms in white, Ti atoms in light blue, and Si atoms by orange sticks.

Close modal

The electronic charge density difference (Δρ) is a simple but effective method to observe the movement of electronic charge within a system constituted by several individual components, such as water and Ti-LTA for the water/1Ti and water/0_b systems. Equation (12) is used to obtain (Δρ),

Δρ=ρ(water TiLTA)ρ(water)ρ(TiLTA),
(12)

where ρ(water TiLTA) is the electronic charge density of the system water/Ti-LTA, and ρ(water) and ρ(TiLTA) are the densities of water and Ti-LTA, respectively, calculated separately but retaining the same box size and calculation conditions as for water/Ti-LTA. Figures 6(c) and 6(d) show Δρ(water 1Ti) and Δρ(water 0_b), respectively, where we clearly observe the movement of charge density from the OH bonds of water towards the H2OTi region. The determination of Bader atomic charges confirmed this charge transfer.66–68 The water molecule adsorbed in 1Ti yielded a charge of +0.030 e, whilst the value for water in 0_b was +0.042 e. The increase in the positive electronic charge of water from 1Ti to 0_b is expected; water interacts with both Ti atoms simultaneously in 0_b, and the LUB of 0_b efficiently surrounds the water molecule [see Fig. 6(f)]. The shape of the positive section of Δρ(water 0_b) shows that the electronic charge moves evenly towards both Ti atoms.

Figure 7 shows the projected density of states (PDOS) onto the O atom of water in the gas phase and water in the systems water/1Ti and water/0_b. Also shown are the PDOS onto the Ti atoms of these systems. We observed that the main indications of hybridization between the electronic states of water and the Lewis acid site are provided by the analysis of the water molecule’s orbitals with symmetry 1b1, which is the highest occupied molecular orbital (HOMO), 3a1 (HOMO-1) and 1b2 (HOMO-2) [see Fig. 7(a)].

FIG. 7.

Projected density of states (PDOS) onto the valence orbitals of [(a)-(c)] O atoms and (d) Ti atoms. (a) Oxygen atom of a water molecule in the gas phase; the charge density isosurfaces (0.0100 bohr−3) of the molecular orbitals of water with symmetry 1b2, 3a1, and 1b1 are provided for each peak. (b) Oxygen atom of the water molecule adsorbed in the β-cage of single Ti-substituted LTA. (c) Oxygen atom of the water molecule adsorbed in the β-cage of 0_b. (d) Titanium atoms of the single Ti-substituted LTA (black line) and 0_b structure (red line) when water is adsorbed in the β-cage (both Ti atoms in 0_b show equivalent profiles). The vertical red dashed line indicates the Fermi energy (shifted to 0.0 eV). The charge density isosurfaces of the most intense peaks within the 1b1 regions (indicated by vertical red arrows) are shown at the right-hand side: (*) 179th band of water/1Ti (0.0008 bohr−3) and (**) 118th band of water/0_b (0.0015 bohr−3).

FIG. 7.

Projected density of states (PDOS) onto the valence orbitals of [(a)-(c)] O atoms and (d) Ti atoms. (a) Oxygen atom of a water molecule in the gas phase; the charge density isosurfaces (0.0100 bohr−3) of the molecular orbitals of water with symmetry 1b2, 3a1, and 1b1 are provided for each peak. (b) Oxygen atom of the water molecule adsorbed in the β-cage of single Ti-substituted LTA. (c) Oxygen atom of the water molecule adsorbed in the β-cage of 0_b. (d) Titanium atoms of the single Ti-substituted LTA (black line) and 0_b structure (red line) when water is adsorbed in the β-cage (both Ti atoms in 0_b show equivalent profiles). The vertical red dashed line indicates the Fermi energy (shifted to 0.0 eV). The charge density isosurfaces of the most intense peaks within the 1b1 regions (indicated by vertical red arrows) are shown at the right-hand side: (*) 179th band of water/1Ti (0.0008 bohr−3) and (**) 118th band of water/0_b (0.0015 bohr−3).

Close modal

The molecular orbital 1b1 corresponds to the non-bonding lone pair of water [see charge density isosurfaces in Fig. 7(a)], and it shows the largest evidence of hybridization with the states of the Lewis acid. The peak that represents 1b1 in the gas phase scatters over energy regions that span from −2.7 to 0.0 eV for water/1Ti and from −3.0 to 0.0 eV for water/0_b (with the Fermi energy shifted to 0.0 eV). The 3a1 orbital also scatters but covering a smaller range of energy. These regions match the profiles of PDOS(Ti), which is confirmation of the states’ hybridization. For instance, the visualization of the electronic charge density corresponding to the most intense peaks within the 1b1 region of systems water/1Ti and water/0_b (indicated with red arrows in Fig. 7) clearly shows continuous isosurfaces that go from the water molecule to the Ti atoms.

At the same time, the feature of 1b2 as a single and intense peak is retained even after the adsorption of water on the Lewis site, in a section of energy where there are no peaks associated with PDOS(Ti). This means that 1b2 remains practically intact, without mixing with the electronic states of the Lewis acid site.

2. Vibrational analysis

We have simulated the IR spectra of the 2MRs with the GGA methodology and compared them with the profiles belonging to the equivalent structures without adsorbates, as shown in Fig. 8. The main feature related to the formation of the 2MRs, which may serve to identify this structure in any experiment, is the splitting of one of the bands associated with the D4R. For instance, the signal at 598 cm−1 in the pure-silica spectrum evolves into two weak bands at 604 and 569 cm−1 after the 2Si-2MR formation (indicated by small vertical arrows in the spectra). The same tendency is followed by the single Ti-substituted structure, where two weak bands at 604 and 571 cm−1 are observed for Si,Ti-2MR. The spectrum of 2Ti-2MR shows a similar splitting but it tends to be slightly larger than in the previous two cases. The separation between the two weak bands is 35 and 33 cm−1 for 2Si-2MR and Si,Ti-2MR, respectively, increasing to 54 cm−1 for 2Ti-2MR.

FIG. 8.

Simulated IR spectra of the different 2MRs and their comparison with their equivalent non-hydrated structures. (From bottom to top) Pure-silica LTA (pure-Si), and 2Si-2MR; single Ti-substituted LTA (1Ti) and Si,Ti-2MR; configuration 0_b and 2Ti-2MR. The splitting of one of the bands associated with D4R is marked with vertical arrows.

FIG. 8.

Simulated IR spectra of the different 2MRs and their comparison with their equivalent non-hydrated structures. (From bottom to top) Pure-silica LTA (pure-Si), and 2Si-2MR; single Ti-substituted LTA (1Ti) and Si,Ti-2MR; configuration 0_b and 2Ti-2MR. The splitting of one of the bands associated with D4R is marked with vertical arrows.

Close modal

We have performed density functional theory and interatomic potential-based calculations to study the pure-silica zeolite LTA and its Ti-substituted forms. The calculated mechanical properties of the pure-silica LTA predict that this material is more rigid than zeolite-A, which has the same framework type but with a Si/Al ratio of 1. The IP method tends to overestimate the value of the elastic moduli compared with the GGA, a trend that agrees with previous calculations on zeolite MFI.

The energetics of the Ti-substituted LTA suggests that the addition of a second Ti atom within the framework, at close distance from the first Ti, is a favorable process. In fact, both the GGA and the IP method predict that the global minimum corresponds to a double Ti substitution in a 6MR, with the two Ti atoms symmetrically opposing each other with a configuration Ti(Si)2Ti.

The interaction of a water molecule with Ti has a calculated adsorption energy between −37 and −52 kJ/mol, where the adsorption in the β-cage is stronger than that in the α-cage, by at least 4 kJ/mol. The direct dissociation of water is only favored when it leads to the formation of a 2MR between T-sites of different D4Rs. The stability is especially high for the 2MR Ti(OH)2Ti, with an energy that is comparable to non-dissociated water adsorbed at the Ti sites. The electronic analysis shows a net charge transfer from water to the Lewis acid site and a strong hybridization between their electronic states.

The GGA-simulated infrared spectrum of pure-silica LTA shows good agreement with experimental measurements, and it is capable of reproducing the four main bands within an error of −30 to −20 cm−1. The Ti substitution causes the occurrence of new weak signals between 1200 and 700 cm−1, depending on the configuration of the Ti atoms. The main feature related to the formation of the 2MRs is the splitting of the weak band at approximately 600 cm−1, which is associated with the D4R building unit, into two weak bands that are separated by 33 to 54 cm−1.

See supplementary material for the tabulated structural parameters of zeolite LTA after the Ti incorporation. In addition, the energies of the LUB for pure-silica LTA, the single Ti substitution and the structures 0_a and 2_a, and the isosurfaces of the electronic charge density associated with the LUB are shown. Furthermore, the derivation of Eq. (11) from the expression to evaluate the absolute stability of the Ti-substitution is presented.

This work was performed using the computational facilities of the Centre for High Performance Computing in Cape Town (CHPC). Via our membership of the UK’s HEC Materials Chemistry Consortium, which is funded by EPSRC [Grant No. EP/L000202], this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). This work has also used the computational facilities of the Advanced Research Computing at Cardiff (ARCCA) Division, Cardiff University, and HPC Wales. The authors acknowledge the UK Economic and Social Research Council (ESRC) [Grant No. ES/N013867/1] and South Africa National Research Foundation for funding under the UK-SA Newton Fund. P.E.N. thanks the South African Research Chair Initiative of the Department of Science and Technology and the National Research Foundation. C.E.H.-T. acknowledges UCL and EPSRC [Grant No. EP/K009567/2] for a scholarship, and N.H.d.L. thanks the Royal Society for an Industry Fellowship. C.E.H.-T. also thanks Katlego Phoshoko and Kenneth Kgatwane for the invaluable assistance and Professor C. R. A. Catlow for a very fruitful discussion. All data created during this research are openly available from Cardiff University Research Portal at http://doi.org/10.17035/d.2016.0011512285.

1.
D. W.
Breck
,
W. G.
Eversole
,
R. M.
Milton
,
T. B.
Reed
, and
T. L.
Thomas
,
J. Am. Chem. Soc.
78
,
5963
(
1956
).
2.
T. B.
Reed
and
D. W.
Breck
,
J. Am. Chem. Soc.
78
,
5972
(
1956
).
3.
E. L.
Uzunova
and
G. S.
Nikolov
,
J. Phys. Chem. A
104
,
5302
(
2000
).
4.
Z.
Wu
,
S.
Goel
,
M.
Choi
, and
E.
Iglesia
,
J. Catal.
311
,
458
(
2014
).
5.
O.
Cheung
and
N.
Hedin
,
RSC Adv.
4
,
14480
(
2014
).
6.
A.
Corma
,
F.
Rey
,
J.
Rius
,
M. J.
Sabater
, and
S.
Valencia
,
Nature
431
,
287
(
2004
).
7.
D.
Jo
,
T.
Ryu
,
G. T.
Park
,
P. S.
Kim
,
C. H.
Kim
,
I. S.
Nam
, and
S. B.
Hong
,
ACS Catal.
6
,
2443
(
2016
).
8.
B. W.
Boal
,
J. E.
Schmidt
,
M. A.
Deimund
,
M. W.
Deem
,
L. M.
Henling
,
S. K.
Brand
,
S. I.
Zones
, and
M. E.
Davis
,
Chem. Mater.
27
,
7774
(
2015
).
9.
S.
Inagaki
,
Y.
Tsuboi
,
M.
Sasaki
,
K.
Mamiya
,
S.
Park
, and
Y.
Kubota
,
Green Chem.
18
,
735
(
2016
).
10.
J.
Kim
,
J.
Chun
, and
R.
Ryoo
,
Chem. Commun.
51
,
13102
(
2015
).
11.
B.
Tang
,
W.
Dai
,
X.
Sun
,
N.
Guan
,
L.
Li
, and
M.
Hunger
,
Green Chem.
16
,
2281
(
2014
).
12.
X.
Lu
,
W. J.
Zhou
,
H.
Wu
,
A.
Liebens
, and
P.
Wu
,
Appl. Catal., A
515
,
51
(
2016
).
13.
G.
Yang
,
L.
Zhou
, and
X.
Han
,
J. Mol. Catal. A: Chem.
363-364
,
371
(
2012
).
14.
A.
Jentys
and
C. R. A.
Catlow
,
Catal. Lett.
22
,
251
(
1993
).
15.
C. M.
Barker
,
D.
Gleeson
,
N.
Kaltsoyannis
,
C. R. A.
Catlow
,
G.
Sankar
, and
J. M.
Thomas
,
Phys. Chem. Chem. Phys.
4
,
1228
(
2002
).
16.
P. E.
Sinclair
,
G.
Sankar
,
C. R. A.
Catlow
,
J. M.
Thomas
, and
T.
Maschmeyer
,
J. Phys. Chem. B
101
,
4232
(
1997
).
17.
J.
To
,
A. A.
Sokol
,
S. A.
French
, and
C. R. A.
Catlow
,
J. Phys. Chem. C
112
,
7173
(
2008
).
18.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
19.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
49
,
14251
(
1994
).
20.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
21.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
22.
K.
Born
and
M.
Huang
,
Dynamical Theory of Crystal Lattices
(
Oxford University Press
,
Oxford
,
1954
).
23.
P. E.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
).
24.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
25.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
26.
S.
Grimme
,
J. Comput. Chem.
27
,
1787
(
2006
).
27.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
28.
A.
Ambrosetti
,
D.
Alfè
,
R. A.
DiStasio
, Jr.
, and
A.
Tkatchenko
,
J. Phys. Chem. Lett.
5
,
849
(
2014
).
29.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
,
J. Comput. Chem.
32
,
1456
(
2011
).
30.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
(
1976
).
31.
P. E.
Blöchl
,
O.
Jepsen
, and
O. K.
Andersen
,
Phys. Rev. B
49
,
16223
(
1994
).
32.
S.
Baroni
,
S.
de Gironcoli
,
A.
Dal Corso
, and
P.
Giannozzi
,
Rev. Mod. Phys.
73
,
515
(
2001
); e-print arXiv:0012092v1 [arXiv:cond-mat].
33.
J. D.
Gale
and
A. L.
Rohl
,
Mol. Simul.
29
,
291
(
2003
).
34.
J. D.
Gale
,
Z. Kristallogr. - Cryst. Mater.
220
,
552
(
2005
).
35.
36.
R. A.
Buckingham
,
Proc. R. Soc. A
168
,
264
(
1938
).
37.
M. J.
Sanders
,
M.
Leslie
, and
C. R. A.
Catlow
,
J. Chem. Soc., Chem. Commun.
Issue 19,
1271
(
1984
).
38.
F. M.
Higgins
,
N. H.
de Leeuw
, and
S. C.
Parker
,
J. Mater. Chem.
12
,
124
(
2002
).
39.
L.
Tang
,
L.
Shi
,
C.
Bonneau
,
J.
Sun
,
H.
Yue
,
A.
Ojuva
,
B.-L.
Lee
,
M.
Kritikos
,
R. G.
Bell
,
Z.
Bacsik
,
J.
Mink
, and
X.
Zou
,
Nat. Mater.
7
,
381
(
2008
).
40.
N. J.
Henson
,
A. K.
Cheetham
, and
J. D.
Gale
,
Chem. Mater.
6
,
1647
(
1994
).
41.
C. E.
Hernandez-Tamargo
,
A.
Roldan
, and
N. H.
de Leeuw
,
J. Solid State Chem.
237
,
192
(
2016
).
43.
R.
Grau-Crespo
,
A. G.
Peralta
,
A.
Rabdel Ruiz-Salvador
,
A.
Gómez
, and
R.
López-Cordero
,
Phys. Chem. Chem. Phys.
2
,
5716
(
2000
).
44.
R.
Grau-Crespo
,
E.
Acuay
, and
A. R.
Ruiz-Salvador
,
Chem. Commun.
Issue 21,
2544
(
2002
).
45.
C. E.
Hernandez-Tamargo
,
A.
Roldan
, and
N. H.
de Leeuw
,
J. Phys. Chem. C
120
,
19097
(
2016
).
46.
R.
Grau-Crespo
,
S.
Hamad
,
C. R. A.
Catlow
, and
N. H.
de Leeuw
,
J. Phys.: Condens. Matter
19
,
256201
(
2007
).
47.
K.
Momma
and
F.
Izumi
,
J. Appl. Crystallogr.
44
,
1272
(
2011
).
48.
P.
Ravindran
,
L.
Fast
,
P. A.
Korzhavyi
,
B.
Johansson
,
J.
Wills
, and
O.
Eriksson
,
J. Appl. Phys.
84
,
4891
(
1998
).
49.
W.
Voigt
,
Lehrbuch der Kristallphysik
(
B. G.
Teubner
,
Leipzig, Berlin
,
1928
).
50.
R.
Arletti
,
O.
Ferro
,
S.
Quartieri
,
A.
Sani
,
G.
Tabacchi
, and
G.
Vezzalini
,
Am. Mineral.
88
,
1416
(
2003
).
51.
S.
Quartieri
,
R.
Arletti
,
G.
Vezzalini
,
F.
Di Renzo
, and
V.
Dmitriev
,
J. Solid State Chem.
191
,
201
(
2012
).
52.
K.
Niwa
,
T.
Tanaka
,
M.
Hasegawa
,
T.
Okada
,
T.
Yagi
, and
T.
Kikegawa
,
Microporous Mesoporous Mater.
182
,
191
(
2013
).
53.
54.
C.
Sanchez-Valle
,
S. V.
Sinogeikin
,
Z. A. D.
Lethbridge
,
R. I.
Walton
,
C. W.
Smith
,
K. E.
Evans
, and
J. D.
Bass
,
J. Appl. Phys.
98
,
053508
(
2005
).
55.
Y.
Huang
and
E. A.
Havenga
,
Chem. Phys. Lett.
345
,
65
(
2001
).
56.
I. R.
Shein
and
A. L.
Ivanovskii
,
Scr. Mater.
59
,
1099
(
2008
).
57.
A.
Roldan
,
D.
Santos-Carballal
, and
N. H.
de Leeuw
,
J. Chem. Phys.
138
,
204712
(
2013
).
58.
A. R.
Ruiz-Salvador
,
R.
Grau-Crespo
,
A. E.
Gray
, and
D. W.
Lewis
,
J. Solid State Chem.
198
,
330
(
2013
).
59.
S. R.
Bare
,
S. D.
Kelly
,
W.
Sinkler
,
J. J.
Low
,
F. S.
Modica
,
S.
Valencia
,
A.
Corma
, and
L. T.
Nemeth
,
J. Am. Chem. Soc.
127
,
12924
(
2005
).
60.
G.
Yang
,
E. A.
Pidko
, and
E. J. M.
Hensen
,
J. Phys. Chem. C
117
,
3976
(
2013
).
61.
Y.
Huang
and
Z.
Jiang
,
Microporous Mater.
12
,
341
(
1997
).
62.
Y.
Qiao
,
Z.
Fan
,
Y.
Jiang
,
N.
Li
,
H.
Dong
,
N.
He
, and
D.
Zhou
,
Chin. J. Catal.
36
,
1733
(
2015
).
63.
M.
Hunger
,
J.
Kärger
,
H.
Pfeifer
,
J.
Caro
,
B.
Zibrowius
,
M.
Bülow
, and
R.
Mostowicz
,
J. Chem. Soc., Faraday Trans. 1
83
,
3459
(
1987
).
64.
A. A.
Sokol
,
C. R. A.
Catlow
,
J. M.
Garces
, and
A.
Kuperman
,
J. Phys. Chem. B
102
,
10647
(
1998
).
65.
J.
To
,
A. A.
Sokol
,
S. A.
French
, and
C. R. A.
Catlow
,
J. Phys. Chem. C
111
,
14720
(
2007
).
66.
G.
Henkelman
,
A.
Arnaldsson
, and
H.
Jónsson
,
Comput. Mater. Sci.
36
,
354
(
2006
).
67.
E.
Sanville
,
S. D.
Kenny
,
R.
Smith
, and
G.
Henkelman
,
J. Comput. Chem.
28
,
899
(
2007
).
68.
W.
Tang
,
E.
Sanville
, and
G.
Henkelman
,
J. Phys.: Condens. Matter
21
,
084204
(
2009
).

Supplementary Material