A density functional theory study of the structure of pure-silica and aluminium-substituted MFI nanosheets

The layered MFI zeolite allows a straightforward hierarchization of the pore system which accelerates mass transfer and increases its lifetime as a catalyst. Here, we present a theoretical study of the structural features of the pure-silica and aluminium-substituted MFI nanosheets. We have analysed the effects of aluminium substitution on the vibrational properties of silanols as well as the features of protons as counter-ions. The formation of the two-dimensional system did not lead to appreciable distortions within the framework. Moreover, the effects on the structure due to the aluminium dopants were the same in both the bulk and the slab. The principal differences were related to the silanol groups that form hydrogen-bonds with neighbouring aluminium-substituted silanols, whereas intra-framework hydro-gen-bonds increase the stability of aluminium-substituted silanols toward dehydration. Thus, we have complemented previous experimental and theoretical studies, showing the lamellar MFI zeolite to be a very stable material of high crystallinity regardless of its very thin structure. & 2016 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license


Introduction
Zeolites are microporous aluminosilicate crystals widely recognised for their extraordinary characteristics as solid acids [1][2][3], size-selective molecular sieves [4,5], and ion-exchange matrices [6,7].Accordingly, the practical applications of these materials are numerous, ranging from water treatment [8] and membranes for gas permeation [9] to fuel and solar cell materials [10,11 ] .Specifically, zeolites find their largest relevance as acid-base catalysts, accounting for more than 40% of the catalysed industrial processes [12], with the environmental and petrochemical sectors among the main users [13].
The confinement of molecules inside the zeolite's microporous network affects the catalytic rate [14] and also determines the reaction selectivity by controlling the movement of reagents, intermediaries or products throughout the pore system [15].However, although the mass-transfer constrain forms part of the exceptional features of zeolites as catalysts, it can also have an adverse effect, when the products get trapped inside the pores, thereby provoking undesired secondary reactions [16].Several schemes have been developed to increase the mobility of molecules toward and from the active sites by combining micropore (o2 nm) and mesopore (2-50 nm) frameworks.This hierarchization of the zeolite pore structure allows an improvement in the catalyst efficiency, owing to faster diffusion to, and better accessibility of, the active sites [17].
The conventional three-dimensional (3D) zeolite structure may be confined into a two-dimensional (2D) one by reducing the size of the crystal along a specific direction to just a few nanometres.These layered materials are obtained through different methods, including direct synthesis or modification of layered zeolite precursors [18,19].The interlayer separation can be tailored within the mesoscale which, together with the micropore network, allows the hierarchization of the system.Eleven zeolite frameworks are currently known to have 2D versions [19].Among them, the MFI type [20] has been the subject of extensive experimental research, mainly by Ryoo and co-workers who reported the first method to synthetise MFI nanosheets.Their method employs diquaternary ammonium-type surfactants as structure-directing agents (SDA) with hydrophobic and hydrophilic functionalities [16].Later works have focused on the control of the thickness and mesopore dimension, together with the function of each segment of the surfactant [21][22][23].
MFI nanosheets show remarkable thermal, hydrothermal and mechanical stability [16,23,24].Furthermore, this material retains the acid strength displayed by the bulk crystal, which, in conjunction with the highly enhanced external surface, allows the treatment of larger molecules than normally fit inside the pores [23,25].Although this structure is very thin ( $ 2 nm), it conserves the crystallinity within the a-c plane and hence the pore system.Therefore, the nanosheet continues to show selectivity in the Beckmann rearrangement [26] and the epoxidation of large molecules [25], as well as the disproportionation and alkylation of toluene [27] and the selective oxidation of benzene [28].
Layered MFI has an effective lifetime as a catalyst.The nanosheets combine small pore lengths along the straight channels (b-axis) with an enhanced availability of active sites at the external surface (see Fig. 1).As a result, not only is the residence time reduced, but the poisoning effects of secondary reactions taking place inside the pores are also diminished [16,26,28].Accordingly, MFI nanosheets are an exceptional material which maintains the bulk-type features, but with additional improvements, such as increased mass-transfer and accessibility for bulky compounds.
Although the experimental work on MFI nanosheets is extensive and has shown the interesting potentials of this type of structure, computational tools can provide a fundamental understanding of their properties at the atomic level, thereby helping to speed up the route to MFI applications.For instance, Varoon et al. reported the optimisation of the layered structure of MWW and MFI zeolites for different thicknesses using Car-Parinello molecular dynamics [29], although they investigated the pure-silica structure under relaxation restrictions.Park et al. have carried out interatomic potential calculations to understand the effect of surfactants in the final arrangement of the nanosheets [22].Molecular dynamics studies have examined the capacity of MFI to act as a membrane for seawater desalination [30] and argon adsorption in layered zeolites [31].
The properties of layered MFI remain of interest, both from a fundamental point of view and for new applications.In this study, we have carried out a systematic examination of a number of MFI nanosheets, including the substitution of silicon atoms in each independent T-site by aluminium.Furthermore, a complete vibrational analysis was performed of both silanols and bridging oxygens binding a proton as counter-ion.Finally, the main features in the electronic structure before and after aluminium inclusion were analysed by means of the projected density of states (PDOS) scheme.

Computational methods
We have carried out the present work using two techniques: (i) atomistic simulations based on the Born model of ionic solids [32] and (ii) Density Functional Theory (DFT) calculations [33,34].The first method was used to study possible terminations of the MFI zeolite surface perpendicular to the [010] direction (the normal vector to the MFI nanosheet plane is collinear with [010] [16]).For this purpose we employed the METADISE code [35], considering the surface dipoles in accordance with the approach by Tasker [36,37].The model definition is outlined in Section 1S of the Supplementary Information (ESI).By using these interatomic potentials the energy of the system comprised two main contributions: one arising from the Coulombic interactions computed using the Parry technique [38,39], and the other from the shortrange repulsions and Van der Waals attractions defined by Buckingham potentials [40].In the parameterisation of the interatomic interactions we have used the potential model for SiO 2 derived by Sanders et al. [41], together with the compatible parameters for the hydroxide ion by Baram and Parker [42].Sander and coworkers have shown the suitability and accuracy of the Born model (initially designed for the study of ionic solids) to simulate covalent materials in the form of SiO 2 polymorphs (α-Quartz, α-Cristobalite, Coesite and α-Tridymite) [41].The main modifications consisted of the addition of harmonic bond-bending terms to describe the rigidity of the SiO 2 tetrahedra, and the inclusion of the shell model to effectively describe the polarizability of the O atoms within the structure [41].Later works have shown the validity and transferability of the model when studying different zeolite frameworks [43][44][45].
For the DFT calculations we have employed the Vienna Ab-initio Simulation Package (VASP) [46][47][48][49] to carry out full structural relaxations of the pure-silica and Al-substituted MFI within the slab model.The maximum kinetic energy to describe the valence electrons was set to 500 eV, while their nodal features and interactions with the internal electrons of the atoms were considered through the projector-augmented-wave method (PAW) [50,51].The Generalised Gradient Approximation (GGA) within the Perdew, Burke and Ernzerhof (PBE) functional [52] accounted for the exchange and correlation effects of the electronic system under the Born-Oppenheimer approximation [53].The Grimme method [54] was applied with the DFT methodology (DFT þD2) to include the effect of the long-range interactions [55][56][57][58].The Monkhorst-Pack [59] scheme was used to generate the grid for the numerical integration over the Brillouin zone; we have only considered the Gamma point due to the large size of the MFI unit cell.The threshold was set at 10 À 5 eV for the electronic iterations and any movement of the ions was stopped when the forces were smaller than 0.03 eV/Å.The threshold for the electronic iterations was further decreased to 10 À 7 eV in order to improve the accuracy of the forces when performing the vibrational analysis.The electronic convergence was improved via smearing of the electronic states by a fictitious temperature of k B T¼ 0.1 eV, although the final total energies were reported at 0 K [60].The O-H frequencies were calculated using the finite difference technique, where positive and negative displacements along each Cartesian coordinate were performed for each atom; these displacements (70.015Å) are sufficiently small to keep the system within the harmonic limit.The Hessian matrix (second derivative of the energy with respect to the atomic displacements) was then calculated and its diagonalization allowed us to obtain the vibrational frequencies.
The slabs were constructed by truncating the bulk system along the [010] direction as shown in Fig. 1.Two values for the thickness were used: (i) a single unit cell, or equivalent, consisting of two pentasil layers (312 atoms), and (ii) a second slab represented by three pentasil layers (456 atoms).The dangling Si-O bonds were saturated with H, thus generating silanol groups.The area of the slab surface was 272.5 Å 2 .In order to minimise any interactions between periodic images, we set the thickness of the vacuum layer that separates successive slabs at 20 Å.In addition, the correction for the dispersion forces was considered up to a maximum distance of 15 Å (this cut-off was used for both bulk and slab calculations for consistency).
The visualisation of the zeolite structure was generated with the code Visualisation for Electronic and Structural Analysis (VESTA 3) [61].

Bulk structure
The pure-silica MFI zeolite has a monoclinic structure (P2 1 /n) which undergoes a phase transition toward an orthorhombic symmetry (Pnma) at temperatures above the range 317-325 K; these temperatures decrease with Al content [62].The orthorhombic phase was employed throughout the present work, mainly justified by the synthesis conditions and the use of this material as a catalyst above 320 K.
The starting atomic positions and cell parameters were extracted from the structural database of the International Zeolite Association (IZA) [63] and further optimised whilst keeping the orthorhombic symmetry; Fig. 2 shows the MFI unit cell.For the GGA calculations under periodic boundary conditions, the energy cut-off for the plane wave basis set (500 eV) was still too small to avoid the effects of Pulay stresses if we had chosen the simultaneous optimisation of the volume cell and atomic positions.However, this problem was overcome by performing fixed-volume relaxations of a set of lattices within the range V 0 70.1V 0 , where V 0 was the volume of the cell from the IZA database.The Pulay stress is almost isotropic [64], and hence cancels when the correlation of the lattice energy versus the volume is fitted to an equation of state; in the present case the Birch-Murnaghan equation [65].An additional advantage of this method is that the Bulk modulus is obtained directly from the fitting as an adjustable parameter.The optimisation of the lattice using the Interatomic Potential (IP) technique was carried out by optimising simultaneously the volume and the atomic positions.The final volume of the relaxed IP lattice was then rescaled ( 710% of the volume) to create the set of fixed-volume relaxations to obtain the Bulk modulus by fitting to the Birch-Murnaghan equation.Table 1 summarizes the results of the optimisations using the experimental work of Quartieri et al. [66] for comparison.The distribution of distances and angles for each method are shown in Fig. 1S.
The IP method accurately reproduced the cell parameters of the MFI framework.Although an underestimation of the three vectors was noted, the largest discrepancy (along b) remained below 1%.However, we observed a large discrepancy in the calculated bulk modulus with respect to the experimental reference (see Table 1).Pure PBE overestimated almost isotropically the experimental cell vectors by between 1.3 and 1.7% while the bulk modulus was 33% larger in comparison with the experimental value.Inclusion of the long-range correction was necessary to improve the DFT outcome [67]; the largest mismatch was reduced to less than 1% for PBEþD2, with a and b longer than the experimental values and c slightly shorter.In addition, the bulk modulus was only marginally different from the reported one by less than 1%.As such, the PBEþD2 combination appears suitable for the description of further distortions within the cell derived from the replacement of Si by Al atoms.
The IP Si-O distances spanned a range of 1.587 to 1.606 Å with an average value of 1.597 Å, with the main distribution peaks staying around 1.60 Å (see Fig. 1S).In the case of PBE and PBEþ D2 practically the same minimum and maximum values were found, extending from 1.615 to 1.631 Å with average values of 1.624 and 1.625 Å respectively.The intra-tetrahedral angle distribution for IP ranged from 106.4 to 112.4 °which was larger than PBE (107.0-111.4 °) and PBE þD2 (107.5-111.8 °).However the three average angles were 109.5°,i.e. essentially the ideal tetrahedral angle.Finally, the inter-tetrahedral angles ranged over 26°for IP and 35°F ig. 2. MFI zeolite unit cell.The oxygen atoms were deleted for a clearer view.The framework is created by linking the silicon atoms.The section in red represents the twelve non-equivalent T-sites with their numeration.for PBE and PBEþD2 whose mean numbers were 152, 154 and 149°respectively.
Overall, the Si-O distances differed most when IP was compared with PBE and PBEþD2.The better performance of PBE þD2 over PBE to predict the Si-O character may be associated with the smaller inter-tetrahedral angles.The use of these techniques in the rest of this work was split in two: (i) IP was used to study the relative stability of the hydrated terminations of the (010) surface (see next section) and (ii) the DFTþ D2 was employed to perform all calculations concerned with Al substitution.

Nanosheet structure
The MFI nanosheet was first reported as a layered structure composed of three pentasil sheets [16].Further work showed the possibility of reducing its thickness to two pentasils layers leading to a single row of micropores along the [001] direction [23].W e have carried out a full geometry optimisation with both thicknesses and compared their structural features with the bulk system.
Previous experimental [68] and computational [29,31] studies of layered MFI have not considered any other surface termination than the full pentasil layer.However, we have also investigated other terminations to identify surfaces with maximum numbers of Si-O-Si bridges and minimum numbers of hydroxy groups per terminal Si.As a result, 203 different structures were obtained, each of them carrying a single hydroxy group per Si atom and ranging in the number of hydroxy groups per unit cell from eight to twenty depending on the termination.These surfaces were optimised using IP with five cells representing the bulk and a single cell for the surface (see ESI for more information).We have plotted the energies against the number of hydroxy groups per unit of area, revealing the full pentasil layer as the most likely hydrated surface for the MFI zeolite, perpendicular to the [010] direction (see Section 3S, ESI).Hence, the full pentasil layer with eight hydroxy groups per unit cell at its surface was chosen as the slab model for the remaining calculations employing DFT.
The orientation of the H atoms in the surface silanol groups are expected to change constantly under experimental conditions.However, if this effect is considered, the number of possible configurations to be examined is too large.We have therefore considered two different options: (i) to orientate the H atoms randomly, or (ii) to orientate the H atoms in the same direction toward their nearest silanol O atoms as shown in Fig. 3.The second alternative was chosen and as a result two H configurations were investigated (see Fig. 3).
Non-terminal atoms did not suffer any appreciable distortion when the optimised slabs were compared with the bulk (see Fig. 2S).The distances and angles essentially spanned the same ranges as the bulk-type structure: 1.61 and 1.63 Å for the Si-O distances and 107-112°and 134-170°for the intra-and intertetrahedral angles, respectively.Consequently, the average values remained unchanged, with a 0.001 Å variation for the Si-O distances, and constant 109.5°and150°values for the intra-and inter-tetrahedral angles, respectively.
Two methods have been used experimentally to obtain the nanosheet thickness: (i) AFM showing a thickness of 3.4 70.3 nm [29], which is very similar to the 3.0 nm value calculated here if the hydroxy O atoms are used as limits for the three-pentasil slab, and (ii) transmission electron microscopy (TEM) images (2.5 nm) [68] which are also close to the 2.7 nm obtained here considering the Si heights for the same slab.As to the two-pentasil nanosheets, a maximum value of 1.7 nm, with an average of 1.5 nm, is reported by TEM studies [23] which corresponds to our 1.7 nm if the Si atoms are assumed as slab limits.Overall the nanosheets do not suffer any expansion or contraction in the b direction as a result of the crystal truncation.
Finally, the relative energies of the two hydroxy orientations depicted in Fig. 3 showed a slight preference for configuration B independent of the slab thickness; the two-pentasil configuration B was favoured by 14 kJ/mol, decreasing in the three-pentasil to 5 kJ/mol relative to configuration A. It is perhaps worth mentioning that the distances between the silanol H atoms and their nearest silanol hydroxy groups were smaller in configuration B.
The above results confirm a highly stable two-dimensional framework with practically no modification compared to the bulk crystallinity [29].The minor structural variations in the bond distances and angles are restricted to the atoms related to the silanol groups, which will be analysed in subsequent sections.The overall trends are also shown not to be dependent on the thickness of the sheets.

Bulk structure
The Si substitution by Al atoms within the 3D MFI framework has been studied before using both cluster methods [69][70][71][72] and those based on periodic boundary conditions [67,[73][74][75].W e performed a similar analysis in order to compare these results with those derived from the Al-substituted 2D structure.To this end, we systematically studied the substitution of one Si by Al per orthorhombic unit cell in all non-equivalent T-sites (12 configurations), allowing full relaxation of the atomic positions.To determine the energetic stability of the doped systems, we have mainly focused on the location of the proton near each individual Al-centred tetrahedron, instead of using an Al-siting probability determination, because the Al occupation of a particular T-site is not random but depends highly on the synthesis conditions and Si/Al ratio [76].The protons are added by the calcination of the ammonium-exchanged zeolite with the Al atoms already in the framework.
The position of the proton as counter-ion was also varied among the four O atoms that form the Al tetrahedron.As such, a total of 48 structures were fully optimised using the PBE þ D2 method.A summary of the final structural parameters is compiled in Fig. 4. We observed that the first neighbouring layer around the Al atom captures the main distortions within the structure.
The biggest distortions in the Si-O distances are observed for those tetrahedra sharing corners with Al.They ranged between 1.586 and 1.715 Å, which represents a 0.029 Å decrease of the lower limit and a 0.084 Å increase of the upper limit when compared with the pure-silica bulk (see Fig. 4).The contraction of the Si-O distances was mainly due to Si binding the non-protonated O atom that is directly bonded to Al, ranging from 1.586 to 1.610 Å with 1.599 Å as the mean.The expansion of the Si-O distances was related to Si atoms binding the protonated O, varying between 1.692 and 1.715 Å with an average of 1.704 Å.Further expansions and contractions were found, although to a lesser extent, for the Si-O distances that connect the nearest Si atoms to the Al atom with the rest of the framework (see Fig. 4).Those atoms beyond the first Si-coordination sphere of the Al atom conserved features very close to that of the pure-silica structure.Nevertheless, expansion of some Si-O distances was found for Si atoms at more than one O from the Al, due to H-bond formation between the proton as counterion and a framework O atom.This additional elongation of the Si-O bond will be analysed in detail when discussing the silanol groups.The Al-O bonds showed similar trends, but escalated according to the covalent radius ratio Al/Si: the Al-O bond varied between 1.690 and 1.730 Å, while the Al-O(H) stayed within 1.884-1.951Å.
As to the intra-tetrahedral angles, the main changes were also limited to those Si atoms at less than 4.0 Å from the Al.The (Si-)O-Si-O(-Al) intra-tetrahedral angles varied from 105.2°to 114.3°,a n average 1.1°above the mean of the pure-silica bulk.The inverse trend was noted for (Si-)O-Si-OH(-Al) angles, at an average 3.3°s maller than that corresponding to the pure-silica structure.Meanwhile, the Al atoms showed the same tendency: intra-tetrahedral angles formed by the protonated O were always smaller than those without proton, differing their average angles by 13.3°( see Fig. 4).Thus, either the Si or the Al tended to get closer to a plane formed by the three O atoms of the tetrahedron opposite to the OH group, thereby decreasing the corresponding intra-tetrahedral angles; the inverse trend occurred with the non-protonated O atoms.
The main feature in the inter-tetrahedral angles was the narrower distribution for the OH groups (see Fig. 4), whose angles varied between 129°and 151°, i.e. 11°below the mean of the puresilica framework.This is a clear indication that OH tends to have inter-tetrahedral angles closer to 120°which corresponds to the ideal value of a tri-coordinated species.

Nanosheet structure
The differences encountered between both H configurations (configurations A and B in Fig. 3) for the pure-silica nanosheets were negligible.The same applies for the slab thickness.Therefore, to avoid duplication of the calculations, we chose the two-pentasil slab with H configuration B (see Fig. whose pure-silica version was 14 kJ/mol more stable than A, to examine the effect of a single Al substitution on the nanosheet framework. The number of non-equivalent T-sites for the two-pentasil slab increased to 24.For each framework T-site (20 sites), the proton may be bonded to four different O atoms.In the case of substitutions taking place in silanol positions (4 sites, see Fig. 3) the number of possible positions for the counter-ion reduces to three (one of the O atoms forms part of the silanol group); a total of 92 structures.
Practically none of the features described in the previous subsection changed when the slab framework was compared with the bulk-type structure.The average Si-OH(-Al) distance increased to 1.703 Å, i.e. the opposite effect to that observed for Si-O(-Al) distances; in this case the average value dropped to 1.600 Å.The (Si-)O-Si-OH(-Al) and (Si-)O-Al-OH(-Si) intra-tetrahedral angles tended to decrease (106.3°forSi and 102.2°forAl on average).The inverse trend occurred when the O atom was not protonated but served as a bridge between a Si and an Al (110.7°forSi and 115.7°f or Al in average).As in the bulk, the inter-tetrahedral angle was prone to decrease towards the ideal value (120°) for protonated bridging O atoms (see Fig. The main distortions appeared, however, when the silanols and the Al-substituted silanols were analysed.

Silanols and Al-substituted silanols
Four sets of silanol groups were derived from the Al substitution: (i) without Al as a nearest neighbour, (ii) silanols with Al as a nearest neighbour connected through the non-protonated O atom, (iii) silanols with Al as a nearest neighbour connected through the protonated O atom and (iv) the Al atom replacing the Si at the silanol position.Fig. 5 shows the structural schemes that represent each type of silanol, while in Table 2 we have compiled the values of O-H vibration frequency and O-H bond distance corresponding to the silanol hydroxyl groups.
The O-H bond distances of the silanol groups within the puresilica two-pentasil slab with H configuration B (see Fig. varied within a very narrow interval, from 0.969 to 0.970 Å (see Table 2), Fig. 5. Structural schemes representing the different types of silanol groups at the external surface of the zeolite: silanol of the pure-silica slab; type (i), silanol whose silicon atom does not have the aluminium as nearest neighbour; type (ii), silanol whose silicon atom is directly connected to the aluminium through the non-protonated oxygen; type (iii), silanol whose silicon atom is directly connected to the aluminium through the protonated oxygen; and type (iv), aluminium-substituted silanol.

Table 2 Summary of the O-H vibrational frequencies [ν(O-H)] and O-H bond distances (O-
H dist.) of the different types of silanol groups within the pure-silica and aluminium-substituted slabs.Three values are given for each entry: minimum, average (bold numbers within parenthesis) and maximum values.The silanol types are represented in Fig. 5. and the stretching frequency of these O-H groups remained in the range of 3792-3815 cm À 1 .The average value was 3805 cm À 1 , although this is 60 cm À 1 above the experimental value [25,26].The overestimation of the O-H frequencies is inherent to the harmonic approximation used here to perform the vibrational analysis.Previous reports have shown that the frequencies decrease and approach the experimental references if the anharmonicity is taken into consideration [77][78][79].However, the suitability of the finite-difference method has been reported before in the case of mordenite surfaces using the PW91 functional [80].Furthermore, the comparison between isolated and non-isolated hydroxy groups is of interest when studying their accessibility and protondonor capacity; the shift provoked by the H-bond interaction (500 cm À 1 ) [81,82] is well above the 60 cm À 1 error.

ν(O-H) (cm
Once the Al substitution has taken place, some of the nonsubstituted silanol O-H distances suffered a slight elongation up to 0.973 Å.The related O-H frequencies decreased to a lower value of 3735 cm À 1 .There were no appreciable differences in the average values between those silanols with Al as nearest neighbour and those at more than two bridging O atoms away from the Al (see Table 2).However, there were two exceptions to the average as a consequence of H-bond formation between the non-substituted silanol H atom and the Al-substituted silanol O atom as shown in Fig. 6a and b.In both cases, the silanol O-H frequencies dropped below 3300 cm À 1 and the O-H distances increased above 0.990 Å.The first resulted from the vicinal interaction of silanols when the Al occupies the T9 position (Fig. 6a), and the second derived from the interaction between the Al-substituted silanol T10 and the silanol T7 two O atoms farther away (Fig. 6b).However, no H-bond formation between the Al-substituted silanol H and the non-substituted silanol O atoms was observed.As shown in Fig. 6c and d,

Table 3
Classification of the protons according to the framework oxygen binding them within the BULK and SLAB.The classification ranges from Inaccessible protons to Accessible protons.Also, accessible protons may be further divided into those with and without intra-framework H-bonds.Additionally, the oxygen-hydrogen bond distance (O-H dist.), its stretching frequency [ν(O-H)] and the hydrogen distance to the nearest framework oxygen atom are provided a The oxygens corresponding to the 12 T-sites closer to the surface termination are listed in the table.The most notable variations were obtained for this set of oxygen atoms.
b Oxygens binding protons that form intra-framework H-bond.The channel, sinusoidal (sin) or straight (str), where the proton is located is set within parenthesis.The column labelled as OH•••O states within parenthesis with which oxygen the hydrogen is forming the H-bond.
c Oxygens binding protons without intra-framework H-bond interactions.
the (Al)O-H•••OH(Si) distance was always larger than 2.5 Å.Meanwhile, the O-H stretching frequencies for the Al-substituted silanols increased with respect to the conventional silanol average by at least 20 cm À 1 .This highlights a weaker Al-OH surf bond which may evolve into Lewis acid formation by dehydration.
Regarding the silanol-related structural parameters (see Fig. 5S), the Si-OH surf distances were the properties that differed most compared to the bulk values.In the case of the pure-silica slab, an elongation of the Si-OH surf bond is noted with an average of 1.637 Å versus 1.625 Å of the bulk.The Al inclusion increased this value further, if the silanols were bridged through the nonprotonated O to the Al atom.The opposite tendency occurred with the protonated bridging O, reducing the Si-OH surf distance to an average of 1.628 Å.
Intra-framework H-bonds also caused the elongation of the Si-O distances for Si located more than two O atoms away from the Al.For instance, in Fig. 6e, if the Al is substituting the silanol T12-O24 and the proton is on O12, it forms a 1.755 Å H-bond with O9.Consequently, the T9-O9 and T10-O9 distances increase to 1.654 and 1.671 Å, respectively.This effect is repeated for every intraframework H-bond formed in the bulk and slab; O atoms with the correct environment to establish this kind of interaction are compiled in Table 3.
Overall, the changes in the zeolite structure after the Al substitution have a local character that does not extend beyond the first sphere of tetrahedra around the Al centre, which applies identically to both the 2D and 3D systems.Also, the crystalline structure formed by corner-sharing SiO 4 tetrahedra is highly stable, making the 2D segment indistinguishable from the equivalent portion within the 3D framework (see Figs. 2S, 4 and 5S).Therefore, differences in chemical behaviour will only be derived from the morphology modification, because the increase in the external area versus the total volume will allow larger molecules to have access to the active sites.In addition, Lewis centres will easily form at the external surface once the dehydration of Al-substituted silanols has taken place; in the 2D system the role of the Lewis acidity will increase relative to that of the Brønsted acids.

Lewis centre formation
The dehydration process of Al-substituted silanols, with the subsequent formation of Lewis centres, was considered by analysing Al-silanol T9(Al)-O25.This position has the additional feature of an intra-framework H-bond between O9 and O12, when O9 (connected to T9) binds the proton; the O9-H•••O25 distance was 1.694 Å after optimisation.Also, this silanol is in a vicinal position with respect to silanol T10-O26 (see Fig. 6c).We have studied three different paths for the dehydration reaction, each accounting for a different framework O atom that simultaneously binds the proton and the Al atom (for T9 site, O8, O9 and O18).We followed the sequence O8, O18 and O9, as shown in Fig. 7.
When O8 is protonated, the O8-H bond is almost parallel to the surface, and thus far from O25 (see Fig. 7a).The first stage of the dehydration is featured by the vertical tilt of the proton with an activation energy of merely 2 kJ/mol, which points to a constant flipping of the O8-H bond at normal temperatures.The energy barrier for the proton transfer toward O25 is 13 kJ/mol.Thereafter, the lattice energy decreases to a stable configuration where a water molecule is adsorbed on the Lewis centre.This final arrangement is 25 kJ/mol more stable than the initial Al-silanol.Bučko and co-workers, studying Al-silanols at the (001) surface of mordenite, found an energy barrier of 10 kJ/mol for the proton transfer, and a further system stabilisation of 45 kJ/mol following the dehydration.These results allowed the authors to explain why it is less likely to find Brønsted acids at the external surface of zeolites: the less stable external acid sites are destroyed even at modest temperatures, releasing water during the process [85].
In the second dehydration path (proton attached to O18, see Fig. 7b), the O18-H bond is almost perpendicular to the surface, placing the proton close to O25.From here, we have found a single energy barrier of 11 kJ/mol to be overcome to form water. Again, the configuration with water adsorbed on the Lewis centre is 28 kJ/mol more stable than the Brønsted acid.
The third case involved the deprotonation of O9, whose proton forms a strong H-bond (1.694 Å) with the neighbour O12 (Fig. 7c).We have used the H configuration A (see Fig.  The formation of a 2MR with a vicinal silanol group, once the Lewis centre has been created, is considered (d).The Nudge Elastic Band Method (NEB) was used to optimise the structures along the reaction path [83,84].The aluminium atom is represented by a black ball, oxygens by dark grey balls, hydrogens by light pink balls, and silicons are located at the interception of the light grey sticks.
Al-silanol T9(Al)-O25 and silanol T10-O26 when O9 is binding the proton (see Fig. 6a).Although the orientation of the O9-H bond is parallel to the surface, with the proton far from O25, we did not find an intermediate configuration which puts the proton in a more favourable position to be transferred, as it was the case for the path in Fig. 7a.Also, the energy barrier increases from around 10 kJ/mol in the previous two configurations to 20 kJ/mol (see Fig. 7c).We related this outcome to the occurrence of the H-bond O9-H•••O12, which makes the proton less prone to transfer, while stabilizing the Al-silanol.However, this effect does not invert the stability relation between the dehydrated arrangement and the Alsilanol, with the latter continuing to be less stable by 12 kJ/mol.
We have also studied the formation of two-membered T-O rings (2MR), which occurs when the tri-coordinated Al atom (Lewis centre) is in a vicinal position to a silanol group, hereafter Lewis/silanol configuration, when the silanol hydroxy group is no longer shared by the Al and the Si.Bučko et al. showed that the 2MR further stabilizes the Lewis centre by 9 kJ/mol, with an energy barrier from the Lewis/silanol configuration toward the 2MR of 26 kJ/mol.Also, these authors found that the most favourable position for the silanol H atom is to remain attached to the silanol O atom even after the 2MR formation [85].In the present work, we used the slab with H configuration B (see Fig. 3) turning the O-H bond of the vicinal silanol T10-O26 to favour the Al-OH(-Si) interaction.Our results show practically the same energy barrier (25 kJ/mol in Fig. 7d) as in Ref. [85].However, although the final 2MR continues to be more stable, the energy difference with respect to the initial Lewis/silanol configuration (3 kJ/mol) is less pronounced than in Ref. [85].
Overall, the dehydration products of Al-silanols are energetically favoured, although the presence of intra-framework H-bonds may affect the kinetics of the process.Vicinal silanol groups further promote the formation of Lewis centres since the formation of 2MR stabilizes the freshly formed Lewis acid.In the particular case of the pentasil layer, which forms the MFI nanosheets, half of the silanols are grouped as vicinal silanols (see Fig. 3).

Protonated framework oxygen atoms
We studied the probability (P i ) of the proton binding to each of the O atoms of the Al-centred tetrahedron using Boltzmann statistics: = ∑ () where N is the degeneracy of each structure (in our present case equal to 1), k is the Boltzmann constant, and T is the equilibrium temperature, chosen as 550 °C, which corresponds to the temperature used in a catalytic fast pyrolysis process [86,87].The index i of the partition sum goes from 0 to 3 in the case of Al-substituted silanols, or from 0 to 4 in any other framework substitution.We assumed that the thermal effects were identical within each set of three or four structures, where only the O atom binding the proton changes, with lesser modifications of the rest of the structure during the optimisation.This allowed us to use the calculated lattice energies at 0 K, relative to the most stable configuration inside the set, instead of the free energies at 823 K (550 °C); the entropic effects, zero-point energy contributions and variation in the internal energies cancel out under the constantvolume condition of our calculations.Boltzmann distributions allowed us to analyse the tendency of the proton to bind to each of the four O atoms around a single Al-centred tetrahedron for the 3D and 2D systems; Fig. 6S shows the individual distributions at a temperature of 550 °C.These results indicate that Al substitutions occurring in the interior of the nanosheets are less prone to change the most probable O to bind the proton compared to the 3D system; just T3 and T6 varied.Meanwhile, the inclination of the proton to bind a different O atom in each Al-substituted T-site was more appreciable in sites exposed at the external surface: six out of twelve distributions changed the relative order of the four O atoms from the bulk to the slab (see Fig. 6S).
At the same time, the acid O-H stretching frequencies varied between 2953 and 3730 cm À 1 (see Table 3).This broad range included the contributions of protons forming intra-framework H-bonds (i.e.H-bond O12-H•••O9 in Fig. 6e) whose frequencies ranged from 2953 to 3364 cm À 1 and isolated protons (without intra-framework H-bond interactions) between 3608 and 3730 cm À 1 .Fig. 7S shows a complete comparison of the 3D and 2D frequencies, where the main variations are derived from disruption/formation of intra-framework H-bonds (protons binding O7 and O10).
In addition, protons immersed in intra-framework H-bonds can be divided into accessible or inaccessible protons depending on whether they are able or not of interacting with molecules loaded in the main channels through any hindrance of the surrounding framework atoms.Meanwhile, all isolated protons are accessible from the pore system.Table 3 lists a compilation of the above classification relating the type of proton with the O atom binding it.
First, inaccessible protons are bonded to O atoms O4, O9, O10, O12, O13 and O16 where H-bond formation can be established between the pairs (O9, O12) and (O13, O16) with O-H distances ranging from 1.000 to 1.014 Å and the O-H stretching frequency from 2953 to 3212 cm À 1 (see Table 3).In the slab, the protons bound to O9 or O12 become exposed on the external surface and, due to the H-bond interaction O9(12)-H•••O12(9), the O-H bond is approximately parallel to this surface.Also, a new H-bond is established between O10 and O4 decreasing the O-H frequency from 3509 to 3364 cm À 1 .
Accessible protons forming intra-framework H-bonds are bonded to O atoms O3, O6, O7, O14 and O19.H-bond interactions are formed for the pairs (O3, O6) in the sinusoidal channel, and (O7, O7) and (O14, O19) in the straight channel.The O-H distance varied between 0.976 and 1.004 Å and the O-H stretching frequencies ranged from 3147 to 3366 cm À 1 .As a consequence of the formation of the 2D system, the H bound to O7 was no longer involved in an H-bond interaction and instead remained in the entry (or mouth) of the straight channel (its frequency increased from 3366 to 3708 cm À 1 ).Meanwhile, the H involved in the H-bond O3(6)-H•••O6(3) along the sinusoidal channel became approximately parallel to the external surface in the 2D structure.
Finally, isolated protons (distances to their nearest neighbouring O atom remained above 2.41 Å, see Table 3) are bonded to O atoms O1, O2, O5, O8, O11, O15, O17, O18, O20, O21, O22, O23, O24, O25 and O26.For these H atoms the O-H distance changed between 0.975 and 0.982 Å and the stretching frequency from 3608 to 3730 cm À 1 .Once the 2D structure is formed, more than half of the O atoms available for isolated protons are at the mouth of the straight pores (see Table 3).

Projected density of states
We have used the PDOS analysis to compare the differences in the electronic structure between the 3D and 2D systems.The 2p states of O atoms O24(-T12) and O12(-T8) were selected to measure the effects of the Al substitution in T sites exposed at the external surface of the slab.O24 is part of a silanol group and O12 remains part of the slab framework.Fig. 8 shows the zeolite segment where both O sites are located.
Fig. 9a shows the O(2p) PDOS for silanol O24 (bound to T12) within the pure-silica two-pentasil slab.Its most intense peak is located close to the Fermi energy (set at 0.0 eV), which overlaps with theoccupiedbandedgeofthetotalDOS(greybackground).Also,the O24(2p) PDOS spreads over a wide range resembling each of the total DOS features.The Al substitution at site T12 (Fig. 9b) shows the appearance of two narrow and intense peaks within 1.0 eV from the Fermi energy.In addition, the O(2p) PDOS for O24(-Al,T12) are much more localised (see Fig. 9b), pointing to a weaker Al-OH surf bond for the Al-substituted T12 silanol, which keeps the 2p states mainly confined in the O atom.This explains the much higher tendency of the Al-substituted silanol O atoms to form H-bonds with vicinal silanols (see Fig. 6 and related discussion).
Meanwhile, the analysis of framework O atoms shows the opposite behaviour when one of its Si is replaced by Al.In the present case, O12 (bound to T8 and T12) was chosen to represent this feature.As for the pure-silica silanol O24, O12 has an intense peak just under the Fermi energy when there is no Al substitution (Fig. 9c).However, once the Al is placed in T8 (Fig. 9d) there is a depression in the PDOS close to the uppermost occupied edge as a result of the three-fold coordination (see Fig. 8b).The O12 atom when protonated forms intra-framework H-bonds with O9 (see Fig. 6f), although the results shown in Fig. 9c and d are applicable to O atoms without any possibility of H-bond interactions.

Conclusions
The catalytic performance of MFI depends on its composition and structure which may be as thin as a single unit cell length along the [010] direction.We have carried out a detailed analysis of the pristine and Al-doped two-and three-dimensional systems, using both IP and DFT calculations, where the Van der Waals correction was necessary for a precise reproduction of the experimental cell parameters and the bulk modulus of the zeolite.
We have modelled two different slab thicknesses, and as shown before by other experimental and theoretical reports, the puresilica MFI nanosheet structure does not change appreciably with respect to the 3D system.Once we have replaced the different T-sites by an Al atom, the same distortions were observed in the nanosheet framework as in the bulk structure.
The silanol groups are also affected by the Al substitution.When Al atoms were located at silanol positions, their tetrahedron distorted in such a way that allowed strong H-bond formations with vicinal silanol groups with a red shift of approximately 500 cm À 1 of the O-H vibrational stretching frequency.The H-bond was always established between the silanol H atom and the Alsubstituted silanol O atom; we never observed the inverse interaction.The O-H bond for Al-substituted silanol was slightly stronger than its Si counterpart, as shown by the vibrational analysis, which is a consequence of the weaker Al-OH interaction.Intra-framework H-bonds have an important role in the kinetics of the dehydration process of Al-substituted silanols.Without these intra-framework H-bonds the energy barrier is around 10 kJ/mol, but when they are present it may increase to 20 kJ/mol.In all the studied configurations, the Al-substituted silanols were less stable than the system with one water molecule adsorbed on the Lewis centre, making the dehydration of Brønsted acids a natural process at the external surface of the zeolite.The Lewis acid may be further stabilized by the formation of two-membered T-O rings with vicinal silanol groups.
The accessibility of protons as counter-ions for molecules loaded in the pore system can be classified according to the framework O atom binding them, (i) inaccessible protons, (ii) accessible protons forming H-bonds with vicinal framework O, and (iii) accessible protons without any kind of H-bond interaction.When the 3D structure is truncated, half of the O atoms binding accessible protons were located at the pore mouth; the rest were forming part of the silanol groups or intra-framework H-bonds.The only appreciable variations in the O-H stretching frequencies for protons as counter-ions, when the 3D-to-2D modification occurs, were associated with disruption/formation of intra-framework H-bonds.
Finally, the O(2p) PDOS analysis for framework O atoms revealed that the peak intensities close to the Fermi energy decrease when protonated due to the Al substitution.However, a different trend is obtained for the silanol hydroxy O when a silanol Si is substituted by Al.In this case, the density of states close to the Fermi energy shows narrower and more intense peaks, together with a shift to higher energies.

Fig. 1 .
Fig. 1.(a) Lateral view of the three-pentasil slab (the framework is created by linking the silicon atoms, oxygen atoms are not shown).The pentasil layer at the centre of the slab is contained between dotted lines.The twenty four silicon atoms per unit cell at the external surface and those of the straight channel are highlighted with darker bonds.The sinusoidal channels run along the interception of two pentasil layers.(b) Top view of the three pentasil slab.The sinusoidal channel which connects neighbouring straight channels is highlighted with darker bonds.

Fig. 3 .
Fig. 3. Pure-silica slabs formed by two (a) and three (b) pentasil layers; the black lines denote the supercell box.The silanol hydrogen atoms are oriented towards their nearest silanol groups.Following this scheme two configurations for the hydrogen atoms are obtained, labelled along the text as configuration A (c) and B (d).The silicon atoms are located at the interception of the grey sticks and the framework oxygens were eliminated for a better view.The silanol oxygens are represented with dark grey balls and the hydrogens with light pink balls.Each silanol is labelled according to the numbers of the T-site and oxygen atom that carries the hydrogen.

Fig. 4 .
Fig. 4. Scattering of the T-O distances (a), O-T-O angles (b) and T-O-T angles (c) for the framework T-sites and oxygen atoms in the Al-substituted bulk (horizontal lines in black) and two-pentasil slab (horizontal lines in blue).The horizontal lines establish the scattering ranges for the distances and angles; the average values are marked with small vertical lines.The structural scheme, besides each set of lines, highlights in red the distance or the angle to which the range of values is referring.

Fig. 6 .
Fig. 6.Local structure of different aluminium-substituted sites highlighting the features of specific hydrogen-bonds.The silicon atoms are located at the interception of the light grey sticks, and the non-important oxygens are not shown for a better view.Silanol oxygens and framework oxygens involved in the hydrogenbonds are represented by dark grey balls, the hydrogens by light pink balls and the aluminium by a black ball.
for this specific configuration to avoid the formation of the H-bond between the

Fig. 7 .
Fig. 7. Energy profile of the dehydration reaction of aluminium-substituted silanol T9-O25, when the proton is binding the oxygen atoms O8 (a), O18 (b) and O9 (c).The formation of a 2MR with a vicinal silanol group, once the Lewis centre has been created, is considered (d).The Nudge Elastic Band Method (NEB) was used to optimise the structures along the reaction path[83,84].The aluminium atom is represented by a black ball, oxygens by dark grey balls, hydrogens by light pink balls, and silicons are located at the interception of the light grey sticks.

Fig. 9 .
Fig. 9. Projected density of states (PDOS) of the O(2p) states (black) and total DOS (grey, reduced 100 times with respect to O(2p) projections).(a) Silanol oxygen O24 in pure-silica slab.(b) Silanol oxygen O24 when the aluminium substitution takes place at the T12 site of the slab (see Fig.8a).(c) Framework oxygen O12 at the external surface of the pure-silica slab.(d) Protonated oxygen O12 when the aluminium substitution takes place at the T8 site of the slab (see Fig.8b).The Fermi energy (0.0 eV) is highlighted by a dotted line.

Fig. 8 .
Fig. 8. Local structure of the aluminium-substituted zeolite framework.(a) Aluminium substitution at the T12 site, the silanol oxygen atom O24 is labelled for reference within the text.(b) Aluminium substitution at the T8 site, the protonated oxygen atom O12 is labelled for reference within the text.Silanol oxygens and framework oxygens are represented by dark grey balls, the hydrogens by light pink balls and the aluminium by a black ball.The silicons are located at the interception of the light grey sticks.

Table 1
Cell vectors and bulk modulus (K) of the orthorhombic MFI zeolite (Pnma) calculated with Interatomic Potentials (IP), pure PBE and GGA together with the longrange dispersion correction (PBE þD2).