Acid strength of zeolitic Brønsted sites—Dependence on dielectric properties

Catalysis Today xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Catalysis Today journal homepage: www.elsevier.com/locate/cattod

Acid strength of zeolitic Brønsted sites—Dependence on dielectric properties ⁎

Marcin Rybicki , Joachim Sauer Institut für Chemie, Humboldt-Universität zu Berlin, Unter den Linden 6, 10099 Berlin, Germany

A R T I C LE I N FO

A B S T R A C T

Keywords: Brønsted acid Deprotonation energy Dielectric constant DFT calculations Acidity descriptors

The dependence of deprotonation energies of zeolitic Brønsted acid sites on the dielectric constants of these materials is investigated using quantum chemical calculations. The deprotonation energy decreases as the reciprocal value of average static dielectric constant of the zeolite framework increases with a correlation coefficient R2 equal to 0.74 (252 data points for 19 zeolite frameworks). We decompose the deprotonation energy of these nanoporous materials into an intrinsic deprotonation energy and a proton “solvation” energy within zeolites which is the origin of the observed 1/ε dependence. The intrinsic deprotonation energy is shown to be a superior acidity descriptor. In agreement with experiments, it indicates that two-dimensional zeolites (nanosheets) are not more acidic than three-dimensional ones (bulk materials). We use our large data set to examine the dependence of the intrinsic deprotonation energy on structure parameters (AleOeSi bond angle, OeH bond distance, confinement coefficient) and on the structure relaxation energy but for none of them a correlation has been found.

1. Introduction Zeolites are ordered nanoporous aluminosilicates which play a crucial role in the petrochemical industry as catalyst for alkane cracking [1]. As pure silica analogues, they are also of big interest in the electronics industry, since the rapid miniaturisation of microprocessors requires development of insulators with very low dielectric constants [2]. The dielectric constant of most zeolites, with their porous structure, is lower than that of traditional dense silicon dioxide, and their mechanical strength and chemical stability are very high [2–4]. Thus these materials represent a very promising alternative as substrates of integrated circuits. The catalytic activity of protonated zeolites originates from AleO (H)eSi groups which represent Brønsted acid sites. Quantification of their acidity is very important for the development of better catalysts. The acid strength of zeolitic Brønsted sites can be characterised in different ways [5]. Most closely related to catalysis are test reactions that use a representative feed molecule such as n-hexane cracking [6]. Since proton transfer is the initial activation step of Brønsted catalysis, the heat of ammonia adsorption measured either by temperature programmed desorption [7] or calorimetry [8] is often used as an acidity descriptor. If interactions between the reacting partners before and after the reaction are neglected, or cancel each other, the deprotonation energy can be used as a reactivity parameter that does not refer to any particular reaction partner [9]. For molecules in the gas phase, mass spectrometry yields proton ⁎

transfer equilibrium constants that are converted into gas phase acidity and gas phase basicity scales, e.g., ref. [10]. For surface sites, such experiments are not possible and a linear relation has been suggested [11,12] between the OH frequency shift on complex formation with a weak base molecule and the energy of deprotonation (Bellamy–Hallem–Williams plot) [13]. However, the OH stretching frequency shift on adsorption of CO [14,15] is an acidity measure that refers to an “early” stage of the reaction, whereas ammonia adsorption and deprotonation refer to a “late” stage of the catalytic reaction when the proton has been already transferred. Hence, results derived using different base molecules and different reference acids have been found to scatter over a range of 90 kJ/mol for a given zeolite (H-ZSM-5) [16]. The deprotonation energy of a zeolite (ΔEDP ) is the energy of the process, Z-OH → Z-O− + H+

(1) −

where Z-OH represents the whole zeolite, Z-O is the deprotonated zeolite and H+ represents the proton in vacuum. The energy of the proton in vacuum is zero, thus the deprotonation energy is an energy difference between the deprotonated and the protonated zeolite:

ΔEDP = EZO−−EZOH

(2)

Quantum chemical calculations of deprotonation energies can be performed for both molecules and solids on equal footing [9], and they have been extensively used to get insight into the factors that determine the acidity differences between different sites in a given zeolite or

Corresponding author. E-mail address: [email protected] (M. Rybicki).

https://doi.org/10.1016/j.cattod.2018.04.031 Received 19 February 2018; Received in revised form 12 April 2018; Accepted 15 April 2018 0920-5861/ © 2018 Elsevier B.V. All rights reserved.

Please cite this article as: Rybicki, M., Catalysis Today (2018), https://doi.org/10.1016/j.cattod.2018.04.031

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

and, therefore, may depend on the dielectric constant of the zeolite. Since the dielectric constant of zeolites is proportional to the T-site density [34], which in turn is proportional to the average pore volume of the zeolite, one can also expect that the deprotonation energy is proportional to the inverse of the average pore volume. The dependence of the deprotonation energy of Brønsted sites in various zeolitic frameworks on their T-site density was recently investigated using periodic density functional theory (DFT) by Jones and Iglesia [26] and a linear dependence of deprotonation energy on the T-site density has been observed. However, the authors discussed this behaviour as a technical issue of the particular periodic DFT implementation. Another global parameter that has been employed to characterizes zeolites, the average Sanderson’s electronegativity [28], was also shown to depend linearly on the refractive index [35], which is the square root of the dielectric constant of the material. The aim of this work is to investigate the expected 1/ ε dependence of the deprotonation energy of zeolites. We will perform calculations on different zeolites with a large variety of dielectric constants. We will show that the dielectric constant belongs to the global factors that determine to some extent the deprotonation energies of Brønsted sites. Having a large set of structural data, we also check if deprotonation energies correlate with simple structure parameters like AleOeSi bond angles or OH bond distances. Previously it has been concluded that there is not a single structure parameter that can explain the variation of deprotonation energies [16]. Given the recent interest in descriptors for computer-assisted materials design [36] we reconsider this case for our large data set.

Fig. 1. The stepwise model of the deprotonation process within a zeolite. For definitions, see text.

between zeolite frameworks with different structures or composition [16–26]. However, it was not possible to correlate these differences with a single parameter, e.g. the SieOeAl bond angle. With respect to different Si/Al ratios it was found that local effects, the number of AlO4 tetrahedra in the second coordination sphere of the particular site, is important for the Brønsted acidity [8,19,21,27], while early attempts failed to identify a global parameter that determines the acidity, such as Sanderson’s electronegativity [28,29]. When comparing thin-film models [30,31] of two-dimensional zeolites with three-dimensional bulk zeolites (H-CHA), we found that the significantly lower deprotonation energies of two-dimensional zeolites are related to a significantly lower effective dielectric constant [32]. The stepwise model of the deprotonation process (see Fig. 1) may explain this relation. In this model, the acidic proton is initially shifted to infinite distance from the parental Brønsted site, but it still stays inside the zeolite pore. In this first step the OeH bond of the Brønsted site has to be broken, intr which is done at the cost of ΔEDP , which we call the intrinsic deprotonation energy. In the second step the proton inside the pore, which may be treated as an ion “solvated” inside the dielectric medium, has to be transferred to the vacuum. This “solvation” needs to be overcome, thus the energy of −ΔEsolv,H+ has to be provided in the second step of the deprotonation process. The total deprotonation energy, therefore, may be subdivided into two terms: intr ΔEDP = ΔEDP −ΔEsolv,H+

2. Models Since we are interested in the dependence of the deprotonation energy on the static dielectric constant of zeolites, we prepared structures of various framework densities. The latter correlates both with the static and the high-frequency dielectric constant of zeolites [34]. The selected silicate framework structures were taken from the IZA (International Zeolite Association) database, and the database of hypothetical zeolites [37]. Two two-dimensional structures were added: the 2dH bilayer of SiO4 tetrahedra [38], and the simplest model of a MFI nanosheets [39], constructed from an MFI unit cell (2dMFI) by cutting the topmost pentasil layer (1PL) and saturating the dangling bonds with protons (2dMFI-1PL). The selected model systems are summarized in Table 1. The static dielectric constants of these pure silica models range from 1.43 (2dMFI-1PL) to 3.71 (MVY). The corresponding Brønsted site models were prepared by exchanging selected silicon atoms with aluminium, and adding a H atom to one of the oxygen atoms of the AlO4 tetrahedron. This way, we prepared 252 different Brønsted site models. The QM cluster models for use in QM:MM calculations were prepared as in our previous work [32], namely by selecting only the tetrahedral atoms which were within a defined cut-out distance from aluminium. If inserting an additional T-atom results in forming a closed ring of T-sites, this new atom was also included in the cluster. A cut-out radius of 600 pm already gives deprotonation energies converged with respect to the cluster size [32]. Migues et al. [40] who used a similar approach also found fast convergence with the cluster size. The QM clusters (Fig. 2) were always terminated with OH groups with a fixed OH distance of 96.1 pm.

(3)

The “proton solvation”, ΔE solv,H should depend on the dielectric constant of the material, which may be easily explained with the classical Born model [33]. In this approach an ion of radius rion and charge e∙zion is transferred from the vacuum (dielectric constant ε0 ) to the medium with dielectric constant ε . A simple electrostatic analysis of this process results in the following equation: +,

Born ΔGsolv =−

2 2 NA z ion e 1 1 ⎛1− ⎞ = C ⎛1− ⎞ 8πε0 rion ⎝ ε ⎠ ⎝ ε⎠

(4)

The solvation Gibbs free energy of the Born model is proportional to 1/ ε with the constant, C , which depends on the charge and the size of the ion. In the Born model, the constant C is the same for the H+ ions in all zeolites. Therefore, the energy of the proton solvation should depend only on the dielectric constant of the surrounding medium. If the simple Born model can qualitatively describe the “solvation” of the proton within zeolites, the total deprotonation energy should follow the equation:

1 intr ΔEDP = ΔEDP −C ⎛1− ⎞ ⎝ ε⎠

3. Methods Zeolites are crystalline materials, and our QM:MM calculations take the periodicity of the system into account. Applying periodic boundary conditions to the deprotonation process implies that a negative charge is created in every unit cell. This causes divergence of the electrostatic energy for the deprotonated system, ZO−. This is also an unrealistic situation because in the catalytic process we do not observe

(5) 2

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

3.1. Embedded cluster QM:MM calculations

Table 1 Investigated zeolite frameworks with their T-site densities, ρT (10−3 Å−3), calculated average static dielectric constants, εav (unitless), Al/Si ratio, number of distinguishable T-sites, NT , and number of possible Brønsted sites, NB . Model 2dMFI-1PL 2dH RWY JSR OBW IRR LTA PCODb FAU CHA BEA LEV SOD MOR MFI FER RWR AEN MVY

a

ρT

εav

Si/Al

NT

NB

6.4 7.1 7.4 11.4 11.8 11.2 13.8 11.5 12.8 14.7 15.3 15.6 16.4 17.2 17.7 18.3 18.7 19.7 21.0

1.43 1.55 1.62 2.14 2.16 2.19 2.51 2.52 2.56 2.85 2.85 2.97 3.06 3.24 3.26 3.40 3.64 3.66 3.71

95 63 47 95 75 51 191 191 191 95 255 215 95 95 191 71 127 191 95

8 1 1 4 5 4 1 7 1 1 9 2 1 4 12 4 2 3 2

32 2 3 15 13 13 3 26 4 4 32 6 1 14 48 13 5 12 6

To calculate deprotonation energies of zeolites, we used a mechanical embedding scheme [16,45], which divides the whole system (S) into the inner (I) and the outer region (O). The dangling bonds of the inner part are saturated with link atoms, and the inner part together with the link atoms forms the cluster (C). The total energy E (S) is obtained as follows [16,45]:

E (S) = EQM (C) + EMM (S)−EMM (C)

To the quantum mechanical energy of the cluster, EQM (C) , the results of the periodic description of the lattice by the force field, EMM (S) , are added. The third contribution, EMM (C) , eliminates approximately the double counting of the contributions coming from atoms in the inner region and artificial contributions from link atoms. The forces acting on the atoms are calculated accordingly [45,46]. The link atoms are not moved according to the force acting on them. They are instead fixed on the bond which they terminate [45,46]. If we apply the mechanical embedding scheme, Eq. (8), to calculate the zeolite deprotonation energy and introduce the notation CH, SH, C− and S−, for the cluster C and the total system S of the protonated and deprotonated form, respectively, we obtain:

a We investigate only surface Brønsted sites. Because the AleOH “silanol” termination is not stable, the number of reasonable T-sites is reduced to 8 with respect to the bulk MFI. b Structure number 8,299,368 from database of hypothetical zeolites [37].

ΔEDP = ΔEDP(QM) + ΔEDP(MM)

E−− εav

ΔEDP(MM) =

1 ⋅Tr (ε) 3

EMM (S−)−EMM (C−)−[EMM (SH)−EMM (CH)]

(9b) (9c)

If the cluster is large enough, the structural distortion upon deprotonation decays within the cluster region and the ΔEDP(MM) term consists only of electrostatic interactions, i.e. long-range contributions ΔELR [16]. Hence,

ΔEDP = ΔEDP(QM) + ΔELR

(10)

Effectively, the deprotonation energy of the whole system is the sum of the deprotonation energy of a cluster calculated on the QM level and the long-range term representing the electrostatic interaction between the inner and the outer part of the system. This long range term is described at the MM level. Thus, the choice of an appropriate force field is crucial for the success of QM:MM calculations. 3.2. Details of QM:MM calculations We used the MonaLisa package [47], with DFT for the QM part and the DFT-parametrized polarizable shell-model potential of Sierka and Sauer [43] for the MM part. Only in the case of 2dMFI-1PS system, we introduced an additional type of surface oxygen and hydrogen centres, with parameters optimized for the systems containing surface silanol groups. In Table S1 of the Supporting information, we summarize force field parameters involved in our calculations. All MM calculations were performed using the GULP program [44]. The DFT calculations of energies and forces were performed within the TURBOMOLE package [48] with the B3LYP functional, corrected for dispersion interactions using the D2 parametrization of Grimme [49]. For all ab initio calculations the TZVP basis set named “def2” in the TURBOMOLE library have been applied [50,51]. The B3LYP functional yields deprotonation energies in a very good agreement with the results of CCSD(T) calculations [32]. Deprotonation energies obtained with this combination of QM and MM methods converge to a constant value with increasing cluster size [32], which indicates that the polarizable shell-model force field describes the long-range electrostatic interactions as accurately as the B3LYP functional. The unit cell parameters of all pure-silica structures were optimized using the shell-model potential. These parameters were used to construct the corresponding zeolite models, and were kept fixed during the

(6)

where ΔEDP(PBC) is deprotonation energy calculated within periodic boundary conditions, and εav is an average static isotropic dielectric constant of the system, which is calculated as 1/3 of the trace of the static dielectric tensor of the system, ε:

εav =

(9a)

with:

ΔEDP(QM) = EQM (C−)−EQM (CH)

simultaneous release of a proton in every unit cell. The divergence of the electrostatic energy can be easily eliminated by adding a uniform background charge to the deprotonated unit cell. The artificial interaction between the charged defects can be reduced using the supercell method or approximately eliminated using the method proposed by Leslie and Gillan (LG method) [42]. We have recently shown that for zeolites both of these methods give virtually the same results [32]. Since the LG method of aperiodic corrections is computationally less demanding, it is used in this study to calculate deprotonation energies. Only in the case of two-dimensional zeolites (H-2dH and H-2dMFI-1PS) the supercell method was applied. The LG method relies on the macroscopic approximation. It assumes that the difference between the energy per unit cell for the periodic defect and for the isolated defect is equal to the interaction energy of a periodic array of point charges, E−−, immersed in a dielectric medium with the dielectric constant ε of the system without defect. Therefore, the LG-corrected deprotonation energy can be expressed as: LG ΔEDP = ΔEDP(PBC)−

(8)

(7)

All dielectric constants of the 3D zeolites were calculated using the polarisable shell-model potential [43] within the GULP [44] program at the structures optimized with this force field. This approach was proven to give reliable dielectric constants in good agreement with ab initio results [34]. In the two-dimensional systems we determined the effective dielectric constant from the extrapolation procedure described in ref. [32]. In the following sections we will refer only to (LG- or supercell-) corrected deprotonation energies, and from now on we will call them simply deprotonation energies, ΔEDP . 3

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

Fig. 2. Selected zeolite models: QM:MM models – left column, QM clusters – right column. Grey sticks – part of the system described with shell-model potential; colour balls and sticks - clusters selected for high-level B3LYP+D2 calculations. Color code: green – Brønsted hydrogen, white – hydrogen, red – oxygen, yellow – silicon, grey – aluminium. Link atoms indicated with white sticks. Pictures prepared using the Jmol program [41]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

structure optimisations of protonated, and deprotonated forms of zeolites. For all structure optimisations we used the following convergence criteria: maximum energy change of 1 × 10−5 eV, maximum gradient component of 5 × 10−3 eV/Å, and maximum displacement component of 1 × 10−3 Å.

C = −320.5 ± 10.7 kJ/mol intr ΔEDP,fit

(13b)

may be treated as The fitted intrinsic deprotonation energy, an average intrinsic deprotonation energy of the set of investigated zeolites, whereas C (1−1/ ε ) is an average solvation energy of the proton for a given framework structure. Since the dielectric constant of twodimensional materials was calculated using the extrapolation procedure [32], which is loaded with bigger uncertainties than the values determined for 3D zeolites using analytical expression, we fitted also the points from Fig. 3 excluding these of two-dimensional zeolites. Accordingly, the obtained fitted curve follows the equation:

4.1. Deprotonation energies and dielectric properties The deprotonation energies of the most stable proton position for a given AlO4 site in a given zeolite framework are plotted as a function of the inverse dielectric constant in Fig. 3 (the same plot for all investigated Brønsted site models is presented in Fig. S1 of Supporting information). The figure shows that deprotonation energies are indeed linearly dependent on 1/ εav , with the following least square fit (broken line, Fig. 3 top):

ΔEDP(3D) = 1333.2 ± 7.8−(332.7 ± 21.6) ∙1/εav

(14)

with correlation coefficient R2 equal to 0.51 and the average intrinsic deprotonation energy of 1000 ± 29 kJ/mol. The slope of this line is, within the fitting uncertainty, virtually the same as the slope for the full set of structures (Eq. (13)). However, the correlation coefficient of this curve is smaller, which originates from: (i) removing eight data points of 2dMFI-1PS which fit very well to the linear dependence from Eq. (13) and (ii) removing the data point of the 2dH zeolite which breaks the counterbalance with the strongly outlying data point of the 3D zeolite RWY.

(11)

2

and a correlation coefficient R equal to 0.74. Comparison of this equation with Eq. (5), which one may convert into the form: intr ΔEDP = ΔEDP,fit −C + C⋅1/ εav

= 1346.6 ± 4.5−C = 1026.1 ± 15.2 kJ/mol intr ΔEDP,fit ,

4. Results and discussion

ΔEDP = 1346.6 ± 4.5−(320.5 ± 10.7) ∙1/εav

(13a)

(12)

results in the following deprotonation parameters: 4

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

Fig. 3. Total (top) and intrinsic (bottom) deprotonation energies of the most stable OH Brønsted sites for a given Al position in a given zeolite framework as a function of the inverse of their static average dielectric constant. Black: IRR (circles), MOR (squares), CHA (triangles) and FAU (diamonds); red: MFI (circles), BEA (squares), OBW (triangles) and 2dMFI-1PS (diamonds); green: RWY (circles), JSR (squares), FER (triangles) and LTA (diamonds); blue: AEN (circles), 2dH (squares), SOD (triangles) and PCOD (diamonds); brown: MVY (circles), LEV (squares) and RWR (diamonds). Broken line: least square fit of the data. Numerical data included in Table S2 of Supporting information. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

The strong scatter of calculated points around the fitting lines reflects the large variations of the intrinsic deprotonation energies of Brønsted sites in different zeolites. The intrinsic deprotonation energy obtained by subtracting, from the calculated deprotonation energy, the proton solvation energy (Eq. (13b)) carries the uncertainty of the fitted C parameter which is ± 10.7 kJ/mol. The detailed error analysis (see Section S3 in Supporting information) indicates that, using this methodology, the uncertainty of the intrinsic deprotonation energy for all investigated 3D zeolites is almost constant, between 8 and 9 kJ/mol, whereas in the case of 2D zeolites it is higher, ± 17 and ± 19 kJ/mol for H-2dH and H2dMFI-1PL, respectively. The intrinsic deprotonation energies calculated from Eq. (5) are presented in Fig. 3 as a function of the zeolite inverse dielectric constant. This plot indicates that most of the zeolites with industrial applications, i.e. MFI, FAU, CHA, MOR and BEA, feature very similar intrinsic deprotonation energies, within a narrow range of about 1030 ± 15 kJ/mol. This similarity of acid strength has already been noticed in the literature [26,52]. Nevertheless, there are also zeolites with considerably higher (e.g. RWR, SOD, LTA, RWY) and lower (e. g. MVY, PCOD and 2dH) intrinsic deprotonation energies. Even within the same zeolite, the acidities of Brønsted sites at different crystallographic

aluminium positions may scatter significantly, like in IRR, where deprotonation energies of T3 and T2 sites differ by 60 kJ/mol. The differences between the intrinsic deprotonation energies of various Brønsted sites do not exceed 80 kJ/mol in most cases. Therefore, we do not expect dramatic differences in intrinsic energy barriers for the reactions involving the transfer of zeolitic Brønsted proton. This is supported by DFT (B3LYP) studies by Zheng et al. [53]. They observed a very good linear correlation between the energy barriers of monomolecular cracking of propane and the deprotonation energies of the Brønsted site of different 8T cluster models. According to this theoretical dependence, the differences of intrinsic deprotonation energies observed by us, up to 80 kJ/mol, would correspond to rather small differences of the corresponding propane cracking energy barriers, up to 10 kJ/mol. This has also been observed in experiments, where the intrinsic enthalpy barriers of the monomolecular n-hexane cracking reaction are constant within 10 kJ/mol for FAU, MOR and MFI zeolites (226 ± 12, 235 ± 14 and 236 ± 12 kJ/mol respectively) [54]. Worth mentioning is the large difference of the deprotonation energies of SOD and 2dH zeolitic frameworks, 260 kJ/mol. Even after eliminating the 1/ εav dependence, the difference between the resulting intrinsic deprotonation energies is large, 158 kJ/mol. This is surprising, since the local structure of the Brønsted site of both zeolites is very 5

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

similar (they have very similar vertex symbols [55] of the T-site). We did not find a simple reason for this difference, but the analysis provided in the next subsection suggests that the difference results from the interplay of structure relaxation on deprotonated ion and electrostatic/ electronic effects. Comparison of deprotonation energies of two- and three-dimensional MFI zeolites suggests that 2D zeolites are more acidic than the 3D structures. The deprotonation energy of the Al4-O(H)-Si7 site in MFI is 110 kJ/mol higher than that of corresponding site in the 2dMFI-1PS framework. However, the difference between the intrinsic deprotonation energies of MFI and 2dMFI-1PS is much smaller, about −10 kJ/ mol, which indicates that the acidity of these two structures is virtually the same. This is in line with the catalytic performance for cycloalkane cracking of MFI nanosheets and bulk zeolites [56], which suggest similar acid strength of 2D and 3D structures. This similar acidity of 2D and 3D materials was explained by a compensation effect between binding energy of protonated base molecules on the deprotonated zeolite surface and the total deprotonation energy by Sauer [57]. This shows that the intrinsic deprotonation energy, which to some extent reflects this compensation, may be a better descriptor of zeolite acidity than the total deprotonation energy. We note, however, that for systems with equal dielectric constants, especially for the Brønsted sites of the same zeolite, intrinsic and total deprotonation energies provide the same information, with only a constant shift of solvation energy. Therefore, the comparison of acidity of different Brønsted sites of the same zeolite may be done based on the total deprotonation energy. 4.2. Deprotonation energies and structure parameters To identify the factors affecting the acidity of zeolites we analysed the dependence of deprotonation energies of all Brønsted site models on selected structural parameters. Fig. 4 shows the dependence of the total and intrinsic deprotonation energies of Brønsted sites on the AleOeSi angle. Whereas there is no correlation of the total deprotonation energy on the AleOeSi angle, the smaller scatter of the intrinsic deprotonation energies results in a very weak trend (R2 = 0.12). The same analysis performed separately for each zeolitic framework (Fig. S2 in Supporting information) indicates, that in almost all studied systems the total deprotonation energy decreases with increase of the AleOeSi angle, which was already observed in an earlier work [58]. However, with exception of LEV, JSR, IRR and MVR structures, the correlations are very weak (R2 < 0.7), but better than that observed for the whole set of Brønsted sites. The fitted curve slopes are different for each zeolite framework, and for most cases, are between −0.7 and −2.0 kJ mol−1 degree−1. This suggests that the AleOeSi angle is indeed one of the factors that affects the intrinsic acidity of the Brønsted sites, but, as already noted in the earlier studies [16] not the only one. The plot of the total deprotonation energies as a function of the OeH distance, performed for each zeolite framework separately (Fig. S3 in Supporting information), indicates that there is no simple functional dependence, which was already reported in the literature [16,19,26]. Since, according to Badger’s rule [59], the OeH stretching frequency is proportional to the OeH distance, our results imply that the OeH frequency is not a good measure of the acidity. The relation between the total deprotonation energies and the confinement coefficients, ηO , has also been analysed (Fig. S4 in Supporting information). The latter is defined as [26]:

ηO =

∑ ri−6 i

Fig. 4. Total (top) and intrinsic (bottom) deprotonation energies of the all Brønsted site models as a function of AleOeSi angle. Colour code as in Fig. 3.

(15)

Fig. 5. Dependence of the OeH bond distance of all investigated Brønsted sites on the confinement coefficient. Red points indicate Brønsted sites involved in hydrogen bonding. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

where ri is a distance between the proton and zeolite O atoms (excluding the O atom of the OeH group). The functional dependence is the same as for the attractive part of the Lennard–Jones potentials. There is no simple functional dependence of deprotonation energy on this coefficient. Notably, there is a very good linear correlation between 6

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

an electronic contribution that includes also electrostatic effects. Both of these effects are apparently very sensible to the ring topology of the framework around the Brønsted sites.

OeH distance and the proton confinement coefficient (see Fig. 5), which was first noticed by Jones and Iglesia [26]. We confirm this dependence here for a much larger set of zeolites. Deprotonation energies seem to be unaffected by the hydrogen bonding of the proton with the framework oxygen atoms (see red points in Fig. S2 in Supporting information), defined on the basis of the OHO distance (dOH − O < 2.4 Å ). However, structures involved in hydrogen bonds are characterized by longer OeH distances and constraint coefficients (red points in Fig. 5) and, therefore, higher stretching OeH frequencies. This all shows that OeH bond distances and stretching frequencies are not good descriptors of the acidity, but are sensible to the shape of the framework wall around the Brønsted site. To gain further insight into the factors affecting the deprotonation energies we decompose it into the “vertical” deprotonation energy and rel the relaxation energy, ΔEDP , defined as: rel ΔEDP = E (ZO−//ZO−)−E (ZO−//ZOH)

5. Conclusions Our ab initio studies confirm that deprotonation energies of zeolites correlate linearly with the reciprocal value of the average dielectric constant. This allows the subdivision of the total deprotonation energy into two energetic terms: (i) intrinsic deprotonation energy and (ii) “solvation energy” of the Brønsted proton (see Fig. 1). The latter strongly depends on the dielectric constant of the medium, namely it is proportional to its inverse value. Therefore, in zeolites with small dielectric constants, like 2D networks or large pore zeolites, the “solvation” energies of the Brønsted site protons are larger (smaller in absolute values) than in materials with a high network density (large dielectric constant). Because the “solvation energy” of the proton enters with a minus sign into deprotonation energy (solvation has to be overcome), a lower dielectric constant results in a lower deprotonation energy. This suggests that zeolites with a low T-site density should be more acidic than those with a high T-density. However we find that the intrinsic deprotonation energy is a better descriptor of the acidic properties of zeolites, than the total deprotonation energy. Since the intrinsic deprotonation energies of most of the popular zeolites, like MFI, FAU, CHA, MOR or BEA, are very similar, they are expected to have a very similar acidity. Two-dimensional zeolites are not, in general, more acidic than the bulk materials. However, the particular two-dimensional zeolite, 2dH, seems to be the most acidic out of the materials studied here. We have not been able to identify a single parameter such as the SieOeAl angle that would determine the intrinsic deprotonation energy. We conclude that the values of the deprotonation energy result from a delicate interplay of structural, electrostatic and electronic effects. So far, it was not possible to find one, or even a few parameters, which may describe the acidity differences observed between these complex materials. It will be interesting to see if machine learning can help [61]. For the time being quantum chemical calculations of deprotonation energies cannot be avoided but these days this is possible even for a large number of structures.

(16)

where E (ZO−//ZO−) is the energy of the optimized deprotonated zeolite, whereas E (ZO−//ZOH) denotes the energy of deprotonated Brønsted site at the optimized structure of protonated zeolite. The “vertical” deprotonation energy may be calculated as: v rel ΔEDP = ΔEDP−ΔEDP

(17)

For most of the investigated zeolite frameworks the relaxation energy correlates well with the AleOeSi angle (Figs. S5 and S6 respectively in Supporting information), whereas the vertical deprotonation energy is rather independent of it (with exceptions of LEV and JSR). The total deprotonation energies of Brønsted sites do not correlate well with the relaxation energy (Fig. S7 in Supporting information), while they depend linearly on the “vertical” deprotonation energy (Fig. S8 in Supporting information). This indicates that in many cases deprotonation energy differences between Brønsted sites are dominated by the differences of the “vertical” values, while the relaxation effects are very similar. This “vertical” deprotonation energy strongly depends on the local structure of the zeolitic framework, but also on the electronic structure of both the parent zeolite and the non-relaxed deprotonated structure. In some cases the “vertical” deprotonation energies of different Brønsted sites are very similar, while the relaxation energies greatly differ. Such an example is mordenite (MOR) for which the total deprotonation energy correlates very well with the relaxation energy (R2 = 0.93, see Fig. S7 in Supporting information). The structure relaxation of the charged defects created upon deprotonation strongly depends on the ring structure of the T-site. As it was shown before [60], the Brønsted sites which participate in six-membered rings are able to relax their structures more easily than those participating in smaller or bigger rings, therefore they should be more acidic. To gain further insight we look at different components of the deprotonation energies of SOD and 2dH, which have similar local structures (OeH bond distance, AleOeSi bond angle and the sequence of the silicon rings in which these T-sites are located) but very different intrinsic deprotonation energies. The relaxation energy of SOD, −69 kJ/ mol, is the highest one (lowest in the absolute value) out of the all Brønsted sites investigated here, whereas the relaxation energy of 2dH is much lower, −126 kJ/mol, but still higher than the values obtained for most of the structures investigated here. This is surprising, since with both of them acidic hydrogen is located on the shared edge between two six-membered rings, which should make the relaxation of these deprotonated structures similar. The vertical ‘intrinsic’ deprotonation energies (‘intrinsic’ deprotonation energy minus relaxation energy) for SOD and 2dH, 1146 and 1052 kJ/mol, respectively, differ also strongly. The difference between intrinsic deprotonation energies of both zeolites is 151 kJ/mol, out of which 57 kJ/mol originate from the anion relaxation (relaxation energy) whereas the other 94 kJ/mol, originate from the intrinsic ‘vertical’ deprotonation energy. The latter is

Acknowledgments This work has been supported by German Science Foundation (DFG) within a Reinhart Koselleck grant to J.S. The Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN) is acknowledged for a computational grant (bec00160). M.R. would like to thank Stephen Leach for discussions and language help. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.cattod.2018.04.031. References [1] W. Vermeiren, J.-P. Gilson, Top. Catal. 52 (2009) 1131–1161. [2] Z. Li, M.C. Johnson, M. Sun, E.T. Ryan, D.J. Earl, W. Maichen, J.I. Martin, S. Li, C.M. Lew, J. Wang, M.W. Deem, M.E. Davis, Y. Yan, Angew. Chem. 118 (2006) 6477–6480. [3] H.K. Hunt, C.M. Lew, M. Sun, Y. Yan, M.E. Davis, Microporous Mesoporous Mater. 130 (2010) 49–55. [4] M. Sun, W. Maichen, R. Pophale, Y. Liu, R. Cai, C.M. Lew, H. Hunt, M.W. Deem, M.E. Davis, Y. Yan, Microporous Mesoporous Mater. 123 (2009) 10–14. [5] J. Sauer, J. Mol. Catal. 54 (1989) 312–323. [6] W.O. Haag, R.M. Lago, P.B. Weisz, Nature (London) 309 (1984) 589–591. [7] M. Niwa, N. Katada, Chem. Rec. 13 (2013) 432–455. [8] W.E. Farneth, R.J. Gorte, Chem. Rev. 95 (1995) 615–635.

7

Catalysis Today xxx (xxxx) xxx–xxx

M. Rybicki, J. Sauer

12 (2013) 191. [37] M.M.J. Treacy, I. Rivin, E. Balkovsky, K.H. Randall, M.D. Foster, Microporous Mesoporous Mater. 74 (2004) 121–132. [38] J.A. Boscoboinik, X. Yu, B. Yang, F.D. Fischer, R. Włodarczyk, M. Sierka, S. Shaikhutdinov, J. Sauer, H.-J. Freund, Angew. Chem. Int. Ed. 51 (2012) 6005–6008. [39] M. Choi, K. Na, J. Kim, Y. Sakamoto, O. Terasaki, R. Ryoo, Nature 461 (2009) 246–249. [40] A.N. Migues, A. Muskat, S.M. Auerbach, W. Sherman, S. Vaitheeswaran, ACS Catal. 5 (2015) 2859–2865. [41] Jmol: an open-source Java viewer for chemical structures in 3D. http://www.jmol. org/. [42] M. Leslie, M.J. Gillan, J. Phys. C 18 (1985) 973. [43] M. Sierka, J. Sauer, Faraday Discuss. 106 (1997) 41–62. [44] J.D. Gale, J. Chem. Soc. Faraday Trans. 93 (1997) 629–637. [45] U. Eichler, C. Kölmel, J. Sauer, J. Comput. Chem. 18 (1997) 463–477. [46] M. Sierka, J. Sauer, J. Chem. Phys. 112 (2000) 6983–6996. [47] F. Bischoff, M. Alessio, M. Rybicki, M. John, J. Sauer, Multi-Level Energy Landscapes: The MonaLisa Program, Humboldt-University, Berlin, 2017https:// www.chemie.hu-berlin.de/de/forschung/quantenchemie/monalisa. [48] TURBOMOLE V6.4; A Development of the University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, (2012) http://www.turbomole.com. [49] S. Grimme, J. Comput. Chem. 27 (2006) 1787–1799. [50] F. Weigend, M. Häser, H. Patzelt, R. Ahlrichs, Chem. Phys. Lett. 294 (1998) 143–152. [51] F. Weigend, R. Ahlrichs, Phys. Chem. Chem. Phys. 7 (2005) 3297–3305. [52] F. Haase, J. Sauer, Microporous Mesoporous Mater. 35–36 (2000) 379–385. [53] A. Zheng, B. Han, B. Li, S.-B. Liu, F. Deng, Chem. Commun. 48 (2012) 6936–6938. [54] S.M. Babitz, B.A. Williams, J.T. Miller, R.Q. Snurr, W.O. Haag, H.H. Kung, Appl. Catal. A 179 (1999) 71–86. [55] M. O’Keeffe, S.T. Hyde, Zeolites 19 (1997) 370–374. [56] Y. Seo, K. Cho, Y. Jung, R. Ryoo, ACS Catal. 3 (2013) 713–720. [57] J. Sauer, Faraday Discuss. 188 (2016) 227–234. [58] A. Redondo, P.J. Hay, J. Phys. Chem. 97 (1993) 11754–11761. [59] J. Cioslowski, G. Liu, R.A. Mosquera Castro, Chem. Phys. Lett. 331 (2000) 497–501. [60] M. Sierka, J. Sauer, J. Phys. Chem. B 105 (2001) 1603–1613. [61] B.R. Goldsmith, M. Boley, J. Vreeken, M. Scheffler, L.M. Ghiringhelli, New J. Phys. 19 (2017) 013031.

[9] W.J. Mortier, J. Sauer, J.A. Lercher, H. Noller, J. Phys. Chem. 88 (1984) 905–912. [10] J.E. Bartmess, J.A. Scott, R.T. McIver Jr., J. Am. Chem. Soc. 101 (1979) 6046. [11] E.A. Paukshtis, R.I. Soltanov, E.N. Yurchenko, React. Kinet. Catal. Lett. 16 (1981) 93–96. [12] A.P. Evgenii, E.N. Yurchenko, Russ. Chem. Rev. 52 (1983) 242. [13] L.J. Bellamy, H.E. Hallam, R.E. Williams, J. Chem. Soc. Faraday Trans. 54 (1958) 1120. [14] A. Zecchina, C.O. Arean, Chem. Soc. Rev. 25 (1996) 187–197. [15] K.I. Hadjiivanov, G.N. Vayssilov, Adv. Catal. 47 (2002) 307–511. [16] U. Eichler, M. Brändle, J. Sauer, J. Phys. Chem. B 101 (1997) 10035–10050. [17] M.J. Janik, J. Macht, E. Iglesia, M. Neurock, J. Phys. Chem. C 113 (2009) 1872–1885. [18] H.V. Brand, L.A. Curtiss, L.E. Iton, J. Phys. Chem. 97 (1993) 12773–12782. [19] L. Grajciar, C.O. Arean, A. Pulido, P. Nachtigall, Phys. Chem. Chem. Phys. 12 (2010) 1497–1506. [20] A.G. Pel’menshchikov, E.A. Paukshtis, V.G. Stepanov, V.I. Pavlov, E.N. Yurchenko, K.G. Ione, G.M. Zhidomirov, S. Beran, J. Phys. Chem. 93 (1989) 6725–6730. [21] M. Sierka, U. Eichler, J. Datka, J. Sauer, J. Phys. Chem. B 102 (1998) 6397–6404. [22] C. Lo, B.L. Trout, J. Catal. 227 (2004) 77–89. [23] J. Sauer, Acidic catalysis by zeolites and the active site concept, in: W. Reschetilowski, W. Hönle (Eds.), On Catalysis, VWB Verlag für Wissenschaft und Bildung, Leipzig, 2010, pp. 136–161. [24] A. Redondo, P.J. Hay, J. Phys. Chem. 97 (1993) 11754. [25] M. Brändle, J. Sauer, J. Am. Chem. Soc. 120 (1998) 1556–1570. [26] A.J. Jones, E. Iglesia, ACS Catal. 5 (2015) 5741–5755. [27] B. Gil, E. Broclawik, J. Datka, J. Klinowski, J. Phys. Chem. 98 (1994) 930–933. [28] W.J. Mortier, J. Catal. 55 (1978) 138–145. [29] J.A. Rabo, Catal. Rev. 23 (1981) 293–313. [30] J.A. Boscoboinik, X. Yu, B. Yang, F.D. Fischer, R. Włodarczyk, M. Sierka, S. Shaikhutdinov, J. Sauer, H.-J. Freund, Angew. Chem. 124 (2012) 6107–6111. [31] J.A. Boscoboinik, X. Yu, E. Emmez, B. Yang, S. Shaikhutdinov, F.D. Fischer, J. Sauer, H.-J. Freund, J. Phys. Chem. C 117 (2013) 13547–13556. [32] M. Rybicki, J. Sauer, Phys. Chem. Chem. Phys. 17 (2015) 27873–27882. [33] P.J. Atkins, J. De Paula, J. Keeler, Physical Chemistry, Oxford Univ. Press, Oxford, 2017. [34] R. Poloni, J. Kim, J. Mater. Chem. C 2 (2014) 2298–2300. [35] K. Van Genechten, W. Mortier, P. Geerlings, J. Chem. Soc. Chem. Commun. (1986) 1278–1279. [36] S. Curtarolo, G.L.W. Hart, M.B. Nardelli, N. Mingo, S. Sanvito, O. Levy, Nat. Mater.

8