The accuracy challenge of the DFT-based molecular assignment of 13C MAS NMR characterization of surface intermediates in zeolite catalysis

The influence of the model and method choice on the DFT predicted 13 C NMR chemical shifts of zeolite surface methoxide species has been systematically analyzed. Twelve 13 C NMR chemical shift calculation protocols on full periodic and hybrid periodic–cluster DFT calculations with varied structural relaxation procedures are examined. The primary assessment of the accuracy of the computational protocols has been carried out for the Si–O(CH 3 )–Al surface methoxide species in ZSM-5 zeolite with well-defined experimental NMR parameters (chemical shift, d ( 13 C) value) as a reference. Different configurations of these surface intermediates and their location inside the ZSM-5 pores are considered explicitly. The predicted d value deviates by up to (cid:2) 0.8 ppm from the experimental value of 59 ppm due to the varied confinement of the methoxide species at different zeolite sites (model accuracy). The choice of the exchange–correlation functional (method accuracy) introduces (cid:2) 1.5 ppm uncertainty in the computed chemical shifts. The accuracy of the predicted 13 C NMR chemical shifts for the computational assignment of spectral characteristics of zeolite intermediates has been further analyzed by considering the potential intermediate species formed upon methane activation by Cu/ZSM-5 zeolite. The presence of Cu species in the vicinity of surface methoxide increases the prediction uncertainty to (cid:2) 2.5 ppm. The full geometry relaxation of the local environment of an active site at an appropriate level of theory is critical to ensure a good agreement between the experimental and computed NMR data. Chemical shifts ( d ) calculated via full geometry relaxation of a cluster model of a relevant portion of the zeolite lattice site are in the best agreement with the experimental values. Our analysis indicates that the full geometry optimization of a cluster model at the PBE0-D3/6-311G(d,p) level of theory followed by GIAO/PBE0-D3/aug-cc-pVDZ calculations is the most suitable approach for the calculation of 13 C chemical shifts of zeolite surface intermediates.


Introduction
DFT calculations for the prediction of 13 C NMR chemical shifts for different surface intermediates could provide crucial support to mechanistic NMR studies in zeolite catalysis. DFT calculations have been widely employed to compute the chemical shifts of various organic molecules, 1 solids, 2 and organometallic compounds. 3 NMR chemical shift calculation has found its own niche in zeolite catalysis, namely for the calculation of 27 Al NMR to support the assignment of the experimental spectra. 4 This technique is being actively used in combination with 27 Al 3Q MQMAS NMR spectroscopy for the investigation of Al atom distribution throughout the zeolite framework. 4,5 However, this has not yet been extensively applied for the mechanistic studies of hydrocarbon conversion on zeolites, raising, therefore, concerns regarding the accuracy and spectral resolution of the DFT-based assignments. There are a number of issues that could limit the accuracy of the calculated chemical shifts. First, the choice of the method of calculations, in particular the choice of exchange-correlation functional and the basis set for the electronic structure calculations, is important. 6,7 Secondly, the selection of a model to represent the structure and local zeolite environment of the examined species is also important. 8 The screening of both the calculation method limitations and chosen model effect on the accuracy of calculations could determine the most reliable approach for the prediction of 13 C NMR chemical shifts based on DFT analysis of the species adsorbed on zeolites.
High resolution solid-state 13 C MAS NMR is a very efficient tool for studying the adsorption and transformations of hydrocarbon species in zeolites. [9][10][11][12][13] In particular, this technique can provide information on the pathways of methane chemical activation and conversion in the pores of metal-containing zeolite catalysts. [14][15][16][17][18][19][20][21][22][23][24] The detection and correct identification of the methane activation products, surface bound intermediates, based on observed 13 C NMR chemical shifts could give valuable insights into the reaction mechanism. The combination of the experimental techniques and computational chemistry methods is a powerful approach for the accurate description of the surface products.
Recent studies have successfully employed the 13 C MAS NMR technique to study methane transformation on coppercontaining zeolites. [25][26][27][28][29] The reported spectroscopic data have shown the presence of few different signals with similar 13 C chemical shifts in the range of 50-65 ppm. 30,31 The unambiguous assignment of such signals to specific molecular configurations based solely on the experimental results is very challenging. It is thus not surprising that different assignments have been proposed for the same signals by different research groups. [25][26][27][28][29] In particular, the origin of the signal at 62-63 ppm is highly disputed. 26,28,29,[32][33][34] Initially, this signal was assigned to a [Cu(m-OCH 3 )Cu] structure. 26 Other studies suggested that this signal should be attributed to methanol adsorbed to a dual Cu + ([Cu(m-CH 3 OH)Cu]) 27 or a mono-Cu + site. 28 Herein, we present a study on the advantages and limitations of common DFT strategies for rationalizing the 13 C NMR spectra of zeolite adsorbed reaction intermediates. ZSM-5 zeolite is one of the most widely studied and applied zeolite materials in catalysis and has been selected as the main object of this computational investigation. Both model and methodological inaccuracies have been examined by calculating 13 C NMR chemical shifts of the Si-O(CH 3 )-Al methoxides, chosen as the reference with a reliably determined chemical shift of 59 ppm, placed at the different locations within the zeolite framework. The determination of the prediction ability of the examined calculation protocols for 13 C NMR chemical shifts of various proposed intermediates of the methane activation on copper-containing zeolites has been carried out.

Models
Synthetic MFI-type zeolites (ZSM-5) can exist in the form of either monoclinic (P2 1 /n) or orthorhombic (Pnma, Pn2 1 a, or P2 1 2 1 2 1 ) phase. A reversible phase transformation is possible under various conditions. [35][36][37] For simplicity, the orthorhombic MFI unit cell with lattice parameters of a = 20.241 Å, b = 20.015 Å and c = 13.439 Å as optimized using DFT with an all-silica MFI periodic model was used for the periodic DFT calculations. 38 To compensate for the positive charge of the extra-framework species, one or two framework Si atoms in MFI unit cell were substituted by Al. Si-O(CH 3 )-Al methoxides were placed along the 12 different locations in the zeolite framework (Fig. S1, ESI †). In the case of structures containing BAS, H + was bonded to the Si(T1)-(O À )-Al(T2) oxygen atom. Cu + cations in all models were coordinated to the two oxygen atoms connected to the Al atom at the for the second one. The resulting Si/Al ratios in the lattice of the periodic zeolite models were 95 (1 Al models) and 47 (2 Al models). After the optimization of the fully periodic zeolite structures, cluster models representing the local zeolite environment at the zeolite site were directly constructed from the periodic models.
''Two-ring'' cluster models were constructed for the zeolite frameworks where methoxide located near the channel intersection (Fig. 1a). ''Three-ring'' models were constructed in the cases where methoxide located in the middle zeolite ''ring'' between the other two zeolite ''rings'' (Fig. 1b). For the models where methoxide is located in the cage between the two intersecting zeolite channels, the two-channel clusters were constructed (Fig. 1c). Dangling O atoms were substituted by H atoms at the B1.6 Å distance. 39 To avoid the presence of Al-H bonds in the cluster models, the zeolite framework was extended by one or two Si atoms near the aluminum atoms.

Methods
Periodic DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP 5.3.5). 40 The exchangecorrelation (XC) energy term was described using the generalized gradient approximation PBE functional. 41 To take into account the van der Waals interactions, semi-empirical D3 dispersion correction by Grimme 42 with Becke-Johnson damping 43 was applied. A plane-wave basis set with a cutoff energy of 450 eV in combination with the projected augmented wave (PAW) method was used. 44 Brillouin zone-sampling was restricted to the G point.
The spin-unrestricted calculations of 13 C NMR chemical shifts were carried out using the cluster models directly constructed from the zeolite framework. All finite model cluster calculations were carried out using the Gaussian 16 program (Revision B.01). 45 Two different computational strategies were implemented for the chemical shift calculations of the cluster models directly constructed from optimized periodic models, namely, (i) singlepoint electronic structure (SP) and chemical shift calculation directly on the geometry from periodic DFT or (ii) full optimization (OPT) of the cluster. The terminal framework O atoms in the clusters were replaced by H atoms at the respective positions. 39 Six DFT functionals were used for the cluster calculations: the generalized gradient approximation PBE functional was chosen in order to be consistent with the method used in periodic DFT calculations, hybrid exchange-correlation functionals PBE0 46 and HSE06, 47-50 long-range corrected wB97xD functional 51 and hybrid mPW1PW91 52,53 and B3PW91 54 functionals. D3BJ 42,43 dispersion corrections were implemented to PBE, PBE0 and B3PW91 functionals. In total, 12 different calculation protocols have been examined. 6-311G(d,p) 55,56 basis was used for the cluster geometry optimization. During the cluster geometry optimization all dangling H atoms were fixed in their positions at a Si-H distance of B1.6 Å. 39 The gauge-independent atomic orbital (GIAO) 57 method along with aug-cc-pVDZ basis set was implemented for the NMR chemical shift calculation of the zeolite cluster 13 C atoms. 58 The choice of the functionals and basis set is based on the analysis of the 13 C chemical shift prediction ability of the various techniques carried out in the study. 6

Results and discussion
Zeolite-bound Si-O(CH 3 )-Al methoxide species are characterized using a 13 C chemical shift at 56-59 ppm, depending on the zeolite framework type. 16,22,26,[32][33][34]59,60 Methoxide groups formed within the MFI framework give a signal at about 59 ppm on the surface of H-ZSM-5. 34 Therefore, such species can be used as the benchmark to investigate various possible uncertainties of the multiple DFT calculation protocols. In this study, we explicitly examine how the local environment affects the predicted chemical shifts by placing the Al atom into different T-sites throughout the framework of ZSM-5 zeolite (MFI framework type). There are twelve distinct T-sites in the MFI framework, and ten of them are examined here to analyze the influence of the model, i.e., the position of the Al-site and therefore the location of Si-O(CH 3 )-Al species, on the calculated chemical shift. Note, that there is more than one possible position for the methyl group for each T-site occupied by Al, i.e., the methyl group can be attached to two Si-O-Al-O-Si bridges related to the same Al-site. Therefore, the studied models were denoted further as T(Al)-site-T(Si)-sites''. T1-T5, T7-T11, T9-T10, and T12-T12 methoxides are shown in Fig. 2 as examples. All examined optimized structures of the methoxides are presented in Fig. S1 (ESI †). The relative stabilities of these intermediates computed using periodic PBE calculations and chemical shift values are reported in Tables S2 and S3 (ESI †), respectively.
Firstly, the relative stability of the methoxide intermediates at different zeolite sites was examined and the results are summarized in Fig. 3. The T9-T6 configuration provides the highest stability for the surface methoxide, with stabilities of the alternative structures being within approximately 30 kJ mol À1 . An exceptionally low stability is predicted for the T7-T7 configuration (DE = 82 kJ mol À1 ). However, this structure was still taken into account and its chemical shifts were calculated to see the correlation between the relative energy and the chemical shift.
The calculation of the 13 C NMR chemical shifts of methoxide intermediates by means of a linear response as implemented in VASP 5.3.5 using the full periodic zeolite models (see Section S2, ESI †) gave unsatisfactory results. The respective computed d( 13 C) values were 91.6 ppm (T1-T2 methoxide), 89.6 ppm (T7-T8 This journal is © the Owner Societies 2020 Phys. Chem. Chem. Phys., 2020, 22, 24004--24013 | 24007 methoxide), 91.1 ppm (T9-T9 methoxide) and 89.9 ppm (T9-T9 methoxide), which are all much higher than the experimental value of 59 ppm. An improved accuracy of the predicted chemical shifts for the periodic system could in principle be achieved using a higher level of theory and more stringent computational parameters. The associated increase in the demand for computational resources rendered such a direct periodic modelling approach impractical for this case. Therefore, in this work a cluster modelling approach was used to accurately calculate the NMR spectral parameters of intrazeolite intermediates.
In order to calculate the methoxide chemical shift values, two main approaches were implemented. The single point (denoted as SP) approach was implemented by using the clusters directly constructed from the periodic model without further geometry relaxation. The main point of the SP calculations is to preserve the geometry features derived from periodic calculations which are on the higher ''model'' level. Another protocol is the cluster geometry optimization (denoted as OPT). This approach enables carrying out structural relaxation on a higher level of theory (hybrid GGA) compared to the periodic calculations (GGA), which, in turn, yields more accurate local geometry in the context of the current computational method. However, during the relaxation of the cluster model, essential structural features due to the extended zeolite structure could be lost resulting in artefacts that could not be found in the periodic models. 8 To demonstrate possible inaccuracies arising from the choice of calculation method, 6 different functionals were implemented. Since the same basis set (aug-cc-pVDZ) was used for all calculations, the protocols were named ''FUNCTIONAL(SP)'' or ''FUNCTIONAL(OPT)'' for the single point and optimized models, respectively.
The results of chemical shift calculations for the surface methoxide species are summarized in Fig. 4. One can see that for all methods the 13 C NMR chemical shift of the most unstable T7-T7 methoxide is predicted to be much higher than those of other configurations. For the rest of the species, there is no correlation observed between the relative stability and the value of the chemical shift. On comparing the results obtained from PBE-D3(SP) and PBE0-D3(SP) calculation protocols (Fig. 4), methoxide 13 C NMR chemical shifts calculated using the PBE-D3 method are slightly higher than those obtained using the PBE0-D3 method. The mean chemical shift values are equal, 62.1 AE 0.8 ppm and 61.5 AE 0.7 ppm for the PBE-D3 and PBE0-D3 methods, respectively. Note, these two values are  For PBE-D3(OPT), the deviation of the chemical shift values increases. The notable increase of the T12-T12 methoxide chemical shift is related to a significant change of the methoxide orientation during the PBE-D3/6-311G(d,p) geometry optimization (Fig. S4, ESI †).
The values obtained using the PBE0-D3(OPT) protocol are much closer to the experimental ones with the mean value of 59.0 AE 0.5 ppm. Similar trends are observed for the data obtained using the wB97xD, HSE06 (Fig. 4b), mPW1PW91 and B3PW91-D3 (Fig. 4c) functionals. The results obtained from SP procedures with different exchange-correlation functionals are very close to each other. The XC-functional related uncertainty for the methoxide species derived from the method choice is AE0.3 ppm.
After cluster geometry optimization, all chemical shift values are decreased, with the HSE06(OPT) and mPW1PW91(OPT) results being closer to the experimental value. The configuration of the methoxide species inside the zeolite pores has a profound impact on the predicted chemical shift. In particular, we note a strong increase of the predicted chemical shift value for the T12-T12 methoxide when wB97xD and B3PW91-D3 functionals are used. The mean values of the predicted chemical shifts for different computational methods and model configurations are summarized in Table 1.
The observed differences between the experimental chemical shifts and the values predicted for different methoxide configurations are attributed to the fact that the experimental data can be viewed as a statistical average over all different possible configurations inside the zeolite pores, whereas the calculations provide This journal is © the Owner Societies 2020 Phys. Chem. Chem. Phys., 2020, 22, 24004--24013 | 24009 the deterministic prediction for one specific configuration within the ensemble. In this respect, the computed values using the hybrid GGA functionals with OPT procedure appear to fall very well within the experimental ranges of 13 C NMR chemical shifts for Si-O(CH 3 )-Al in ZSM-5 zeolite (Fig. 5). 31 On the other hand, the GGA-level PBE-D3(OPT) procedure gives far less accurate results. Fig. 6 presents the XC-related uncertainty of the mean 13 C chemical shift values for each calculated methoxide (Fig. 6).
All mean values and deviations are reported in Table S4 (ESI †). It is clearly seen that the obtained chemical shifts are not significantly influenced by the chosen calculation DFT functional in the case of SP calculations. The value deviations are in the range of 0.2-0.5 ppm. Concerning mean 13 C chemical shifts calculated using cluster OPT procedures, it is clearly seen that XC-related uncertainties are larger compared to single point calculations. The results obtained demonstrate that the deviation of the chemical shifts calculated using the cluster optimization model is in the range of 0.8-1.5 ppm except in the case of T12-T12 methoxide (AE3.7 ppm) because of a conformational change due to the cluster relaxation.
The role of structural relaxation is very important for the 13 C NMR chemical shift calculations. 61 It was previously observed that the level of theory employed for the geometry optimization could give rise to deviations in the calculated 13 C chemical shift values for transition metal complexes by up to 15 ppm. 7 In our case, it follows from the reported results that cluster geometry optimization on a higher level of theory (h-GGA) compared to periodic calculations (GGA) is the key to obtaining accurate 13 C NMR chemical shift values. However, cluster geometry optimization could lead to a local conformation change that can significantly change the predicted chemical shifts. The current analysis identifies the PBE0-D3(OPT) technique as providing the best results in terms of the agreement between theory and experiment (59 ppm experimental value vs. 59 AE 0.5 ppm calculated one).
Next, we analyzed the performance of different DFT methods for the prediction of 13  All structures containing [Cu-O-Cu] 2+ sites are believed to be possible products of methane activation using the Cu/ZSM-5 catalyst (Fig. S3, ESI †). The calculated chemical shifts are presented in Fig. 7.
The general trends observed for the model case of the surface methoxides (Fig. 4) are well reproduced for the various hydrocarbon species on Cu/ZSM-5 (Fig. 7). The chemical shifts obtained using the SP scheme depend only slightly on the choice of the XC functional. The exception is SP-PBE, where uniformly higher d values are predicted. Cluster geometry relaxation leads to a significant decrease in the resulting 13 C NMR chemical shifts. SP calculations predict chemical     shifts of 52-53 ppm for methanol adsorbed on BAS, which is in good agreement with the experimental values of 53-54 ppm. 26,28,29,[32][33][34] After the cluster relaxation in OPT scheme, predicted shifts decrease the values to 49-51 ppm, which are close to the chemical shift of ''free'' methanol. 62 For BAS adsorbed DME, SP calculations predict a difference of about 4 ppm between the chemical shifts of the two carbon atoms caused by the differences in O-C interatomic distances in the optimized periodic zeolite models. The mean chemical shift values for each C atom are 68-69 ppm and 64-65 ppm. After the cluster geometry optimization, a conformational change of dimethyl ether takes place (Fig. S5, ESI †). The chemical shift values of two different atoms obtained using PBE-D3, HSE06 and mPW1PW91 techniques are equalized to 58-59 ppm after geometry optimization. However, calculations using PBE0-D3, wB97xD and B3PW91-D3 functionals resulted in slight differences between d values for C-1 (58-59 ppm) and C-2 (60-61 ppm) atoms.
For methanol adsorbed on the Cu + cation, the SP calculations give d values of 63-64 ppm which are closer to the experimental assignment of the signal at 62 ppm 27,28 than the chemical shifts of 58-59 ppm obtained using the OPT procedure. Interestingly, in the case of DME adsorbed on a Cu + cation, no significant difference between the chemical shifts of C-1 and C-2 is observed in spite of the pronounced conformational change during the cluster geometry optimization (Fig. S6, ESI †). The chemical shifts predicted using SP   2+ ) is characterized by the calculated chemical shifts in the range of 70-81 ppm, which is much larger than the previous assignment of the signal at 63 ppm to such sites made on the basis of DFT calculations with a simplistic 10 T-atom cluster model. 27 The calculated 13  The influence of the XC-related uncertainty can be demonstrated by examining the mean 13 C chemical shift values and their deviations for each calculated species (Fig. 8). All mean values are reported in Table S6 (ESI †). For the di-Cu-containing structures and DME on single Cu + site, the XC-related uncertainties are in the range of AE1.5 and AE2.4 ppm (SP), AE2.2 and AE3.3 ppm (OPT). This is probably due to the higher sensitivity of the results for the models containing two heavy atoms to the choice of XC functional.
Indeed, the XC-related uncertainty in the predicted chemical shifts is quite small for DME/BAS, despite the significant conformational change happening during the geometry optimization. Our analysis suggests that the OPT procedure is the most reliable computational protocol in providing computed chemical shifts in line with the experimental data.
The presented results show that the accurate calculation of the 13 C NMR chemical shifts of intrazeolite intermediates requires the use of the h-GGA level of theory. This is in line with the previous computational studies on adsorption and reactivity of zeolite-based catalysts. Although general trends in the energetics of hydrocarbon activation by H-forms 63,64 and TM-modified zeolites 65,66 can be predicted well using the GGA and meta-GGA methods, the quantitative analysis of such systems requires more accurate h-GGA calculations. 63,67 From a set of XC functionals considered here, similar performance is shown for almost all methods with the exception of the GGA PBE-D3 functional. The optimal performance in predicting chemical shifts for C 1 intermediates in zeolite pores is observed for the PBE0-D3(OPT) protocol.

Conclusion
Model and XC functional-related uncertainties of twelve 13 C NMR chemical shift calculation protocols were analyzed. The influence of the model inaccuracy and local zeolite confinement was demonstrated using the 13 C NMR chemical shift DFT calculations of Si-O(CH 3 )-Al methoxide benchmark species located in 12 different positions throughout the ZSM-5 zeolite framework. The analysis was also extended to the possible methane activation products on the surface Cu/ZSM-5 zeolite.
It is established that the effect of local confinement on the calculated chemical shifts influences chemical shifts values up to AE0.9 ppm. The XC-related uncertainty is AE0.3 ppm in the case of single point calculations. However, the cluster geometry optimization increases the value range up to AE1.5 ppm. This effect significantly increases (up to AE3.3 ppm) for the case of structures containing copper atoms. The results obtained indicate that cluster geometry optimization on the hybrid GGA level of theory is the crucial step to obtain local geometry that provides the basis for accurate prediction of the chemical shift values of Si-O(CH 3 )-Al species. Out of the examined calculation protocols, the pPBE-D3BJ/450eV//PBE0-D3BJ/6-311G(d,p)//GIAO-PBE0-D3BJ/aug-cc-pVDZ technique is the most appropriate one for the calculation of 13 C NMR chemical shifts of the hydrocarbon species formed on the zeolite surface.

Conflicts of interest
There are no conflicts to declare.