An NMR Crystallographic Investigation of the Relationships between the Crystal Structure and 29Si Isotropic Chemical Shift in Silica Zeolites.

NMR crystallography has recently been applied to great effect for silica zeolites. Here we investigate whether it is possible to extend the structural information available from routine NMR spectra via a simple structure-spectrum relationship. Unlike previous empirically derived relationships that have compared experimental crystal structures for (often disordered) silicates with experimental NMR spectra, where the structure may not be an accurate representation of the material studied experimentally, we use NMR parameters calculated by density functional theory (DFT) for both model Si(OSi(OH)3)4 clusters and also extended zeolitic SiO2 frameworks, for which the input structure corresponding to the NMR parameters is known exactly. We arrive at a structure-spectrum relationship dependent on the mean Si-O bond length, mean Si-O-Si bond angle, and the standard deviations of both parameters, which can predict to within 1.3 ppm the 29Si isotropic magnetic shielding that should be obtained from a DFT calculation. While this semiempirical relationship will never supersede DFT where this is possible, it does open up the possibility of a rapid estimation of the outcome of a DFT calculation where the actual calculation would be prohibitively costly or otherwise challenging. We also investigate the structural optimization of SiO2 zeolites using DFT, demonstrating that the mean Si-O bond lengths all tend to 1.62 Å and the distortion index tends to <2.0°, suggesting that these metrics may be suitable for rapid validation of whether a given crystal structure represents a realistic local geometry around Si, or merely a bulk average with contributions from several different local geometries.

Introduction

29Si NMR spectroscopy has long been a key tool in the structural characterization of silicate-based zeolites, owing to its moderate natural abundance and receptivity, spin quantum number I = 1/2, and, most importantly, its sensitivity to small changes in the local structure.[1−3] It is, for example, well-known that the ranges of chemical shifts observed for Si(OH)4–(OT) (i.e., Q silicate species, where T = Si) are distinct for different values of n and, in aluminosilicate zeolites, the chemical shift for a given Q Si species differs by ∼7 ppm per next-nearest neighbor Al atom substituted on the T site.[3,4] In one of the most important recent examples of the power of solid-state NMR spectroscopy to provide structural information on silicate zeolites, Brouwer et al. demonstrated that it is possible to solve such structures using only a unit cell determined from crystallographic measurements and the build-up curves from 29Si double-quantum NMR experiments to provide distance restraints.[5] However, this sort of approach is extremely time-consuming, owing to the requirement to record a series of experiments where 29Si double-quantum coherences must be excited between spin pairs at natural abundance (i.e., only 0.22% of all Si pairs). It would, therefore, be desirable to have some means of extracting information from the simple one-dimensional 29Si MAS NMR spectra of zeolites and relating this in some way to their structure.

In the past, this goal has led to many proposed links between the 29Si isotropic chemical shift, δiso, and a variety of structural parameters including the mean Si–O–Si bond angle (⟨θSiOSi⟩, in deg) and the mean Si–O bond length (⟨rSiO⟩, in Å) in a range of silicate minerals, zeolites, and glasses.[3] Examples include the relationship based on a set of 20 silicates with δiso = 875⟨rSiO⟩ – 1509, albeit with substantial scatter[6] and a set of four silica polymorphs and a silicalite precursor that exhibit a very different relationship of δiso = 325.8⟨rSiO⟩ – 633.[7] Other relationships proposed between 29Si δiso and ⟨rSiO⟩ are typically closer to the former than the latter, with δiso = 1447⟨rSiO⟩ – 2432 for sodium and potassium feldspars,[8] δiso = 1218⟨rSiO⟩ – 2058 for several silicates and quartz,[9] δiso = 1372⟨rSiO⟩ – 2312 for albite, natrolite, and two silica polymorphs,[10] δiso = 1187⟨rSiO⟩ – 2014 for a selection of silicates,[11] and δiso = 1126⟨rSiO⟩ – 1909 for Mg2SiO4,[12] leading to a range of descriptions that follow the same general trend; i.e., δiso moves downfield as ⟨rSiO⟩ increases. Hochgräfe et al. used this trend to great effect in the assignment of 29Si resonances in three siliceous zeolites.[13] However, these relationships typically exhibit significant scatter, as shown in Figure a, which suggests that the dependence on a single parameter may be an oversimplification. This is not surprising, given the range of materials from which data points have been taken, the uncertainty associated with the structural parameters and the fact that these relationships aim to describe the magnetic shielding interaction by a single structural parameter.

Figure 1

Plots of published relationships between 29Si chemical shift and (a) mean Si–O bond length[6−12] and (b) mean Si–O–Si bond angles.[7,8,14−19] The lines represent lines of best fit to the experimental data (where data are available). Experimental data are not shown for ref (11), which is a reanalysis of existing literature data and determined the relationship indicated by the broken gray line in part a, and ref (18), where the data points are all included in the analysis of ref (16) but the relationship discussed in the main text is shown as the broken gray line in part b. Experimental data are not available for ref (10) (green line in part a), which is a conference abstract and the numerical values do not appear to have been published elsewhere since. (c) Schematic representation of a general Si–O–T motif, showing the distances and angles used by Sherriff et al. to calculate the contribution to δiso from the dipole moment of the O–T bond.[26]

Similarly, given that it is known that the electronegativity of the Si–O bond relates to the Si–O–Si bond angle,[3] many relationships (as shown in Figure b) have been reported between 29Si δiso and the mean Si–O–Si bond angle, ⟨θSiOSi⟩. These relationships include δiso = −0.603⟨θSiOSi⟩ – 20.8 for four silica polymorphs and a silicalite precursor,[7] δiso = −1.17⟨θSiOSi⟩ + 68.6 for sodium and potassium feldspars,[8] δiso = −0.619⟨θSiOSi⟩ – 18.7 for 13 silica polymorphs and zeolites,[14] δiso = −0.533⟨θSiOSi⟩ – 10.7 for Si(OAl)4 in nine zeolites,[15] δiso = −0.579⟨θSiOSi⟩ – 25.3 for six zeolites,[16] δiso = −0.563⟨θSiOSi⟩ – 9.62 for three sodium disilicate polymorphs,[17] δiso = −0.686⟨θSiOSi⟩ – 8.29 for Si(OSi)4 in three zeolites,[18] δiso = −0.609⟨θSiOSi⟩ – 20.6 for silicalite-1,[19] δiso = −0.79⟨θSiOSi⟩ + 18.18 for 13 leucites and related compounds,[20] and δiso = −0.62⟨θSiOSi⟩ – 1.09 for 33 sodalites with different cage contents.[21] Müller et al. reported a gradient of −0.57 ppm per degree for three dense phases of SiO2 when combined with data for the isostructural AlPO4 phases,[22] although for just the SiO2 phases (using numerical data from Smith and Blackwell[7]), the relationship is δiso = −0.622⟨θSiOSi⟩ – 17.8. As in the case of the relationship between 29Si δiso and ⟨rSiO⟩, although there is a general trend for a decrease in δiso with increasing ⟨θSiOSi⟩, there is a large variation in the gradients and y-intercepts for the assumed linear correlations. It is clear from Figure b that the data discussed above generally fall into several sets of near-parallel lines but with significant scatter for each grouping.

Relationships between 29Si δiso and several other geometric and geometric–electronic parameters have also been investigated, including the mean Si–T distance (T = Si, Al, Ge, etc.),[7,20] the mean O–Si–O angle, ⟨θOSiO⟩,[7−9] sec(⟨θSiOSi⟩), and cos(⟨θSiOSi⟩)/[1 – cos(⟨θSiOSi⟩)],[7,14,15,17,23] and the bond strengths or electronegativities of the adjacent T cations (later modified to account for variation in the Si–O–X angles).[6,8,11,12,14,19,24,25] Sherriff et al. proposed a more complicated but, in principle, universally applicable relationship, where the major contribution to the 29Si shielding was assumed to be from the magnetic susceptibility of the bond between O and the next-nearest neighbor T atom.[26] Their relationshipwhereand the angles and distances, θ, r, R, and D are as shown in Figure c, and r0 is the length of the bond of unit valence (tabulated by Brown and Altermatt,[27] except for Si and Al, for which the respective values of 1.64 and 1.62 Å were redetermined by Sherriff et al.[26]), provided a reasonable prediction of δiso, reported for a range of Si-containing motifs in minerals (giving a root mean squared deviation of 0.66 ppm over 60 silicates).

A major source of possible error in all of the relationships discussed above is that they often seek to compare experimental crystallographic data with experimental NMR spectra. While this is, of course, the ultimate aim of these relationships: to be able to determine structural parameters from an NMR spectrum (or to predict an NMR spectrum from an experimental structure), the techniques are sensitive to structure on very different length scales. As an example, the “Si–O” bonds reported for an aluminosilicate typically (unless the Al is well ordered) represent the weighted mean Si–O and Al–O bond lengths (typically ∼1.6 and 1.7 Å, respectively), whereas the 29Si NMR spectrum will be sensitive to only the Si–O bond lengths. With some of the relationships mentioned above reporting a variation in chemical shift of ∼1000 ppm per Å, a 1 pm error in bond length can have an effect similar to substitution of a neighboring Si for Al (cf. ∼10 ppm per pm and ∼7 ppm per Al). While this may be a somewhat extreme example (or may, in fact, suggest that the substitution of a single Al leads to an increase in ⟨rSiO⟩ of ∼1 pm), other smaller errors relating to temperature effects are also relevant, with magic angle spinning (MAS) NMR spectra typically recorded at just above room temperature whereas crystal structures are typically obtained at lower temperature where thermal motion is reduced. Therefore, in order to determine whether there is any true worth in attempting to relate simple geometric parameters such as mean bond lengths and angles to the experimental NMR spectrum, in this work, we compare the NMR parameters calculated for exactly known model systems (small clusters and extended zeolite-like solids), where the experimental errors in both structure and chemical shift referencing are removed.

The use of empirical structure–spectrum relationships has largely been superseded by the use of quantum-chemical calculations, most notably using density functional theory (DFT). Periodic planewave DFT approaches have made highly accurate calculations of NMR parameters of extended periodic solids almost a routine accompaniment to solid-state NMR spectroscopy.[28−31] At their most basic, the calculations can confirm an assignment, but the ease with which a structural model can be manipulated, perhaps to investigate cation or anion substitution or motion, means that these calculations provide extremely detailed insight into a range of challenging systems that exhibit complex spectra arising from nonperiodic features.[32,33] However, there remain systems where it is too costly to apply DFT calculations of NMR parameters, most notably in molecular dynamics calculations, where the simulation of just a few ps of motion can lead to a “trajectory” comprising many thousands of structural snapshots. Applying DFT calculations to all of these would rapidly lead to computational costs on the order of CPU decades, which is unfortunate, since it is the dynamic processes occurring within many materials, including zeolites, that are of most interest to their applications and NMR should be ideally placed to study these, owing to its sensitivity both to local structure and motion spanning ∼12 orders of magnitude.[33,34] Therefore, in order to provide a bridge between structures and materials where NMR spectra are likely to be of most interest and DFT calculations would prove too costly, we attempt here to determine whether there are any underlying structure–spectrum relationships that can be used to predict NMR parameters with near DFT-level accuracy, without invoking costly computation.

Building on our earlier work on calcined aluminophosphates (AlPOs),[35] in this work, we consider the effect of various local structural parameters on the 29Si δiso for a series of simple model clusters and zeolitic SiO2 frameworks. We show that, by considering multiple geometrical parameters simultaneously, a more robust relationship between spectra and structural parameters can be obtained. Ultimately, we hope that the relationship we have determined will find application in more disordered materials, or in molecular dynamics simulations, where it may not be feasible to calculate NMR parameters using relatively costly DFT methods.

Computational Details

DFT Calculations

Model cluster DFT calculations were carried out using Gaussian 03 (revision D.01)[36] using the continuous set of gauge transformations (CSGT) method to calculate the NMR parameters. The B3LYP hybrid GGA functional was used, with the 6-311+G(d,p) basis set employed for H and O and the aug-pcS-2 basis set (which has been optimized to provide accurate nuclear magnetic shielding parameters)[37] for Si. Prior to the calculation of the NMR parameters, the structures of the clusters were optimized to an energy minimum, with the parameters specified in the text constrained to their stated values. Calculations were carried out using either a local cluster comprising four Intel Core i7–930 quad-core processors with 6 GB memory per core or the EaStCHEM Research Computing Facility comprising a 198-node (2376-core) Intel Westmere cluster with 2 GB memory per core and QDR Infiniband interconnects.

Periodic DFT calculations were performed using version 16.11 of the planewave CASTEP code,[38] which employs the GIPAW algorithm[39] to reconstruct the all-electron wave function in the presence of a magnetic field. The generalized gradient approximation (GGA) PBE[40] functional was employed, and core–valence interactions were described by ultrasoft pseudopotentials,[41] which were generated on the fly. Wave functions were expanded as planewaves with a kinetic energy smaller than a cutoff energy of 60 Ry (816 eV). Integrals over the first Brillouin zone were performed using a Monkhorst–Pack grid with a k-point spacing of 0.04 2π Å–1. Where optimization of the structure to an energy minimum was carried out, this used the same cutoff energy and k-point spacing as above, and with all atomic coordinates and unit cell parameters allowed to vary. Calculations were performed using the EaStCHEM Research Computing Facility, comprising a 54-node (1728-core) Intel Broadwell cluster with 4 GB memory per core and FDR Infiniband interconnects at the University of St Andrews.

Calculations generate the absolute shielding tensor, σ, in the crystal frame. From the principal components of the symmetric part of σ, it is possible to generate the isotropic shielding, σiso = (1/3) Tr{σ}. The isotropic chemical shift is given (assuming σref ≪ 1) by δiso = −(σiso – σref)/m, where σref is a reference shielding, here (for the CASTEP calculations) 289.13 ppm for 29Si, and m is a scaling factor, ideally 1 but, here, 1.3652. The values for σref and m were determined by comparing experimental and calculated chemical shifts for MFI- and FER-type SiO2.[42,43]

Linear Regression

Multivariate linear regression was carried out using the MATLAB[44] routines described in the Supporting Information (S1). All values generated by MATLAB are truncated to five significant figures.

Results and Discussion

Model Cluster Calculations

Using an approach shown earlier to be successful for AlPOs,[35] the influence of ⟨θSiOSi⟩ and ⟨rSiO⟩ on the calculated 29Si σiso was investigated using several series of model Si(OSi(OH)3)4 clusters, shown in Figure a. These clusters allow systematic (and independent) variation of ⟨θSiOSi⟩ and ⟨rSiO⟩ for the central Si, without considering the longer-range effects of an extended zeolitic framework. For investigations into the effect of ⟨θSiOSi⟩, the central SiO4 tetrahedron was fixed with the ideal Si–O length of 1.62 Å and O−Si−O angles of 109.47°, while ⟨θSiOSi⟩ was varied according to Table . When only ⟨θSiOSi⟩ is varied (series 1), there is a strong linear correlation (R2 = 0.972) between σiso and ⟨θSiOSi⟩, with a gradient of 1.04 ppm per degree, which is remarkably similar to that found previously for 31P in AlPOs (1.05 ppm per degree variation in ⟨θPOAl⟩).[35] However, as also observed earlier, there is some deviation from this straight line as the angle approaches 180°. The relationship(where the stated coefficients give σiso in ppm) provides an improved correlation coefficient (R2 = 0.9988) and, crucially, the deviation from the straight line is now less dependent on the angle. The term cos(⟨θSiOSi⟩)/[cos(⟨θSiOSi⟩) – 1][23] gave a poorer value of R2 (0.9844) and was not considered further. Parts b and c of Figure show plots of 29Si σiso against cos(⟨θSiOSi⟩) and the standard deviation of θSiOSi, σ(θSiOSi), for series 1–8. In series 2–6, ⟨θSiOSi⟩ was kept constant at 140° while σ(θSiOSi) was varied as indicated in Table . As can be seen from the inset in Figure b, a difference is observed of up to −6.8 ppm (series 3, n = 5) in σiso relative to the corresponding point of series 1 (n = 4), in which ⟨θSiOSi⟩ = 140° and σ(θSiOSi) = 0. This is similar to our earlier observation for AlPOs that the individual bond angles contribute to 31P σiso, rather than simply the mean bond angle. Series 7, where both ⟨θSiOSi⟩ and σ(θSiOSi) were varied systematically (see Table ), provides further evidence that the individual θSiOSi, rather than just ⟨θSiOSi⟩, are of importance. As can be seen from Figure b, there is a strong linear relationship between σiso and cos(⟨θSiOSi⟩), withalthough there is significant deviation from linearity toward lower cos(⟨θSiOSi⟩) (higher σ(θSiOSi)). To further investigate the contributions of ⟨θSiOSi⟩ and σ(θSiOSi), in series 8, the bond angles were all randomly generated (see the Supporting Information (S2) for values). From Figure b, it can be seen that series 1 and 8 have a very similar relationship between σiso and cos(⟨θSiOSi⟩), with series 8 described byThis similarity to eq suggests that cos(⟨θSiOSi⟩) is a reasonably good predictor of σiso, although, clearly, the variation in individual bond angles leads to some scatter in the shielding for a given mean bond angle (R2 for series 8 is 0.9738, and the mean absolute error (MAE) in σiso calculated by DFT and from eq is 1.01 ppm). Using multivariate linear regression (see the Supporting Information (S1) for more details), the contributions of both cos(⟨θSiOSi⟩) and σ(θSiOSi) to σiso can be determined, withwhich increases R2 to 0.9966 and reduces the MAE to 0.38 ppm for series 8.

Figure 2

(a) Example of a Si(OSi(OH)3)4 cluster used to investigate the dependence of 29Si σiso on the systematic variation of the structural parameters, ⟨θSiOSi⟩ and ⟨rSiO⟩. Atoms are colored blue (Si), red (O), and gray (H). Plots of 29Si σiso calculated for Si(OSi(OH)3)4 clusters against (b) cos(⟨θSiOSi⟩) and (c) σ(θSiOSi). For details of the bond angles used in the model clusters, see Table . The inset in part b shows only values for series 2–6, and series 1 (n = 5), with ⟨θSiOSi⟩ = 140°.

Table 1

Relationships Describing the Systematic Variation of Si–O–Si Bond Angles (θSiOSi(i)), in the Series of Model Si(OSi(OH)3)4 Clusters Studied Here (See Figure a for an Example)a

series	θSiOSi(i) (deg)	N
1	θSiOSi(1) = θSiOSi(2) = θSiOSi(3) = θSiOSi(4) = ⟨θSiOSi⟩ = 115 + 5n	12
2	θSiOSi(1) = θSiOSi(2) = 140, θSiOSi(3) = 140 + 5n, θSiOSi(4) = 140 – 5n	7
3	θSiOSi(1) = θSiOSi(2) = 140 + 5n, θSiOSi(3) = θSiOSi(4) = 140 – 5n	6
4	θSiOSi(1) = 105 + 5n, θSiOSi(2) = θSiOSi(3) = θSiOSi(4) = 140 + (140 – θSiOSi(1))/3	15
5	θSiOSi(1) = 150, θSiOSi(2) = 105 + 5n, θSiOSi(3) = θSiOSi(4) = (410 – θSiOSi(2))/2	9
6	θSiOSi(1) = 120, θSiOSi(2) = 175 – 5n, θSiOSi(3) = θSiOSi(4) = (440 – θSiOSi(2))/2	13
7	θSiOSi(1) = 105 + 5n, θSiOSi(2) = θSiOSi(3) = θSiOSi(4) = 140	15
8	all angles randomly generated,b 107.06 ≤ θSiOSi(i) ≤ 174.96	20

The angles are expressed for the nth member of the series, and the number of clusters in the series, N, is given. For the central SiO4 tetrahedron, the Si–O bonds were fixed at 1.62 Å and the O–Si–O angles at 109.47°.

For a full list of the randomly generated angles, see the Supporting Information (S2).

As discussed above, many attempts have also been made to link σiso with the mean Si–O bond length, ⟨rSiO⟩.[3,6−12] This was investigated using a second set of model clusters, where all O–Si–O and Si–O–Si bond angles were constrained to 109.47 and 140°, respectively, and ⟨rSiO⟩ was varied systematically as given in Table . It can be seen from Figure that, when only ⟨rSiO⟩ is allowed to vary and all other structural parameters are kept constant (series 9), σiso and ⟨rSiO⟩ are related by the quadratic functionwith R2 = 0.9995. This is similar to our previous finding for 31P in calcined AlPOs.[35] In series 10 and 11, the value of ⟨rSiO⟩ was fixed at 1.62 Å, while the standard deviation in the Si–O bond lengths, σ(rSiO), was systematically varied (see Table ). This resulted in differences of up to −1.4 ppm (series 11, n = 6) in σiso relative to the corresponding point of series 9 (n = 0), in which σ(rSiO) = 0. As above for the Si–O–Si bond angles, this suggests that the 29Si σiso is sensitive to the individual Si–O bond lengths, rather than just their average value. In series 12, both ⟨rSiO⟩ and σ(rSiO) were varied systematically (see Table ) and, as can be seen from Figure a, the σiso values for this series are in reasonably good agreement with eq , although it must be noted that the range of σ(rSiO) for series 12 is relatively small compared to those for series 10 and 11, where larger deviations from eq are observed. In series 13, all Si–O bond lengths were randomly generated between 1.45 and 1.85 (see the Supporting Information (S2) for details), giving a maximum σ(rSiO) of 0.14 Å. From Figure a, it can be seen that the data from series 13 describe a very rough parabola, with the best-fit quadratic functionwith a correlation coefficient of R2 = 0.84. Using multivariate linear regression, it is possible to account for variation in both ⟨rSiO⟩ and σ(rSiO), withwhich is very close to eq (in the limit of σ(rSiO) = 0) and improves R2 to 0.98 for series 13.

Table 2

Relationships Describing the Systematic Variation of Si–O Bond Lengths (rSiO(i)), in the Series of Model Si(OSi(OH)3)4 Clusters Studied Here (See Figure a for an Example)a

series	rSiO(i) (Å)	N
9	rSiO(1) = rSiO(2) = rSiO(3) = rSiO(4) = ⟨rSiO⟩ = 1.61 + 0.01n	15
10	rSiO(1) = rSiO(2) = 1.62, rSiO(3) = 1.61 + 0.01n, rSiO(4) = 1.63 – 0.01n	7
11	rSiO(1) = rSiO(2) = 1.61 + 0.01n, rSiO(3) = rSiO(4) = 1.63 – 0.01n	7
12	rSiO(i) = 1.56 + 0.01n, rSiO(2) = rSiO(3) = rSiO(4) = 1.62	10
13	all lengths randomly generated,b 1.45 ≤ rSiO(i) ≤ 1.85	40

The lengths are expressed for the nth member of the series, and the number of clusters in the series, N, is given. For the central SiO4 tetrahedron, the O–Si–O angles were fixed at 109.47° and all Si–O–Si bond angles were fixed at 140°.

For a full list of the randomly generated lengths, see the Supporting Information (S2).

Figure 3

Plots of 29Si σiso calculated for Si(OSi(OH)3)4 clusters against (a) ⟨rSiO⟩ and (b) σ(rSiO). For details of the bond lengths used in the model clusters, see Table . The inset in part a shows only values for series 10 and 11, and series 9 (n = 1), with ⟨rSiO⟩ = 1.62 Å.

Model SiO2 Frameworks

From the model cluster calculations, it can be seen that both the mean Si–O–Si bond angles and Si–O bond lengths, as well as the standard deviations in their values, influence the 29Si σiso, which goes some way to explaining why many of the relationships between a single structural parameter and 29Si chemical shift in the literature disagree to some extent and are not generally transferrable. To investigate whether these findings are relevant in the extended periodic structures of zeolites, where variation in all of these parameters may occur simultaneously and independently, calculations were carried out on a series of model zeolitic SiO2 polymorphs using the periodic planewave code, CASTEP.[38] There is, of course, a large and well-documented effect on the 29Si δiso as the number of next-nearest neighbor Si species is changed, either as a function of condensation (e.g., Q2 Si(OSi)2(OH)2 vs Q4Si(OSi)4 species) or as a function of cation substitution (e.g., Q4Si(OSi)4 vs Q4Si(OSi)3(OAl) species),[3,4] and thus, to avoid complications arising from this, structures taken from the literature (with international zeolite association framework topology codes[45] of EDI, ITG, JBW, MTT, SFE, THO, and VET—see the Supporting Information (S3) for further details) were converted to idealized models where all framework T atoms were 100% occupied by Si. This also provided a charge-neutral framework that allowed for removal of the extraframework cations and H2O within the pores, leading to a set of 7 microporous SiO2 structures containing 49 crystallographically unique Si atoms. The structure of the dense phase α-quartz, containing one unique Si site, was also included. NMR parameters were calculated for these structures before and after optimization to an energy minimum, leading to the consideration of 100 unique Si atoms. As discussed below, the structures showed significantly greater variation in the Si–O bond lengths and O–Si–O bond angles prior to optimization, so all structures were considered here in order to ensure that the study was as widely applicable as possible to the various types of structures that may be encountered in real materials of interest.

Parts a and b of Figure plot the calculated 29Si δiso for the set of 16 structures (i.e., prior to and post optimization) against cos(⟨θSiOSi⟩) and ⟨rSiO⟩, respectively, and it can be seen that there is a strong linear correlation with cos(⟨θSiOSi⟩)but a less apparent correlation with ⟨rSiO⟩ for the “real” data. While the dependence on cos(⟨θSiOSi⟩) is similar to that in eqs and 5 (note the change in sign arises from changing from σiso to δiso), there is still some scatter (R2 = 0.89) and the MAE is 1.23 ppm, which is insufficient to provide a generally useful link between an NMR spectrum and a given structure. As described in the Supporting Information (S1), multivariate linear regression was used to generate a relationship dependent on multiple structural parameters, givingwhere, as noted in the Supporting Information (S1), the ⟨rSiO⟩2 term was discarded, as this is effectively collinear with ⟨rSiO⟩ over the relevant range of Si–O bond lengths (see below). Equation can be seen to contain coefficients whose magnitudes (accounting again for the change in sign from σiso to δiso) are very similar to those found in eqs , 6, and 9 (for cos(⟨θSiOSi⟩), σ(θSiOSi), and σ(rSiO), respectively), with the degree of similarity especially surprising given that the earlier equations relate to calculations for very simple model systems carried out using a completely different code and level of theory. Figure c shows a plot of 29Si σiso calculated using CASTEP against that from eq . It can be seen that there is excellent agreement, with R2 now increased to 0.945 and the MAE reduced to 0.97 ppm. The MAE is now affected mainly by the unoptimized MTT structure, which contains unusually short ⟨rSiO⟩ (1.568–1.594 Å) and the “Al” sites of unoptimized JBW and THO, which have unusually long ⟨rSiO⟩ (1.675–1.749 Å), leading to a discrepancy between CASTEP and eq of up to 4.9 ppm for very short bonds and 4.4 ppm for very long bonds. It could, therefore, be suggested that a quadratic dependence on ⟨rSiO⟩ may be required to make the structure–spectrum relationship more generally useful. However, as discussed below, it is unlikely that such extremes of Si–O bond lengths would be observed in real SiO2 zeolites.

Figure 4

Plots of δiso calculated by CASTEP against (a) ⟨rSiO⟩, (b) cos(⟨θSiOSi⟩), and (c) δiso predicted from the structure by eq for the series of 16 model zeolitic SiO2 frameworks discussed in the text.

To gain some insight into the errors in the coefficients in eq , the structure–spectrum relationship was recalculated for each of the 12870 unique combinations of 8 structures selected from the 16 considered here (see the Supporting Information (S1) for details). Histograms showing the distribution of coefficients determined this way are shown in Figure , and it can be seen that, when only 8 structures are considered, there is significant uncertainty in many of these values, depending on the structure set chosen. However, as shown in the Supporting Information (S4), the distributions of coefficients are essentially independent of one another, with the exception of the coefficient for ⟨rSiO⟩, which is strongly correlated with the intercept (R2 = 0.9929), so that it remains unclear to what extent this structural parameter actually influences δiso. As discussed above, this may result from the approximation that 29Si δiso depends only linearly on ⟨rSiO⟩, leading to sets including one or more structures with unusually long or short Si–O bonds contributing to spurious values of the coefficient. We do not, however, observe any coefficients for ⟨rSiO⟩ approaching the ∼1000 ppm Å–1 mentioned in the Introduction. This significant variation goes some way to explaining the distribution of relationships within the literature, where, depending on the subset of zeolites chosen, it would be possible to obtain very disparate structure–spectrum relationships. It is particularly worthy of note that 17% of the relationships determined found no dependence on ⟨rSiO⟩, whereas the coefficients for cos(⟨θSiOSi⟩) were always nonzero and very similar, suggesting that the contribution from ⟨θSiOSi⟩ to δiso is more universally applicable.

Figure 5

Histograms showing the distribution of values for the coefficients in eq . The values were determined by repeating the parametrization of eq for each of the 12 870 possible combinations of 8 of the 16 model SiO2 frameworks as described in the Supporting Information (S1).

It is worth comparing the results of eq to the relationship of Sherriff et al.,[26] that was also reported to be universal and parametrized using (experimental) data for a wide range of structure types. As discussed in the Supporting Information (S5), for the test set of SiO2 frameworks discussed above, there is significant scatter from the ideal 1:–1 correspondence expected. However, the Sherriff model performs remarkably well for the optimized structures (R2 = 0.96) and very poorly for the unoptimized structures (R2 = 0.50), indicating it may suffer from some overparameterization and be less applicable to more unusually distorted frameworks (see below) than our own model, which was parametrized using a set of structures that included some with more extreme distortions.

Structural Changes upon Optimization

As discussed above, there is some indication that at the extremes of ⟨rSiO⟩ there may be a quadratic relationship between this term and σiso. However, upon optimization of the eight structures considered here, it was observed that the individual Si–O bonds all fall within the range from 1.603 to 1.643 Å (see Figure a), indicating that the very long and short bonds observed for the unoptimized structures above arise from the fact that these were derived from experimental structures for (alumino)silicates containing guest cations and water molecules within the pores. The optimum value of rSiO observed here is in good agreement with the mean value of 1.597(26) Å reported by Wragg et al.[46] in a study of 35 experimental zeolite structures (although, since these structures were not optimized, a range from 1.54 to 1.67 Å was observed for individual bond lengths), with the slight increase observed in the DFT calculations possibly arising from thermal motion of the O atoms[47] (since the structures in the DFT calculations were effectively at 0 K, whereas the experimental structures were obtained at finite temperature). In pure calcined SiO2 polymorphs, then, such a variation in bond lengths is much less likely. However, while the bond lengths of the SiO4 tetrahedra tended toward all being equal, the O–Si–O bond angles did not necessarily optimize to closer to the ideal tetrahedral angle, θ0 = 109.471° but, rather, the distortion indextended to fall within the range of DI ≤ 2.0°, as shown in Figure b (the point on the dotted gray line indicating DI = 2.0° is Si3 of the VET structure, for which DI changed from 3.312 to 1.991° on optimization). When DI was very small in the initial structure, optimization often led to an increase but never above the threshold of 2°. There is no optimum value of θSiOSi, as shown in Figure c, since this parameter is strongly dictated by the framework topology.[46] These observations suggest that, at least for pure silicates, the values of ⟨rSiO⟩ and DI might be used as an indicator for an unrealistic structure. However, the situation becomes more complicated when considering, for example, an aluminosilicate with fractional occupancy of Si sites by Al, which has longer bonds to O and may also be higher coordinate, leading to a superposition of several contributions to the final “SiO4” tetrahedron in the crystal structure, and a wider distribution of ⟨rSiO⟩ and DI might be expected for (disordered) substituted frameworks. Such experimental crystal structures will not, of course, represent accurate descriptions of the true local geometry, even if they are correct for the long-range average structures.

Figure 6

Plots of (a) ⟨rSiO⟩, (b) distortion index, DI, and (c) ⟨θSiOSi⟩ for the set of model silicate framework structures discussed in the main text before (red) and after (blue) optimization using CASTEP. In part b, the dotted gray line indicates the threshold of DI = 2.0°. In all parts, the x axis serves only to separate the distinct Si species.

Applications to Siliceous Zeolites

There are many examples of pure SiO2 zeolites in the literature, where the combination of detailed 29Si homonuclear correlation NMR spectroscopy, high-quality crystallographic data, and, in some cases, DFT calculations has been used to provide a full spectral assignment.[2] Here, we provide two examples to demonstrate the ability of eq to help provide both spectral assignment based on the crystal structure and structural validation based on the NMR spectrum.

The structure of the monoclinic form of MFI-type SiO2 ZSM-5 was determined by van Konningsfeld et al.[48] and contains 24 crystallographically distinct Si species with ⟨rSiO⟩ and ⟨θSiOSi⟩ covering the relatively narrow ranges of 1.589–1.601 Å and 147.11–158.83°, respectively. From the 29Si NMR spectrum of the material, Fyfe et al.[42] were able to resolve and assign 16 resonances or groups of resonances (within a shift window of only ∼7 ppm) based on homonuclear 29Si double-quantum correlation spectra. Figure plots the 29Si chemical shifts predicted from eq (using the structure of van Konnigsveld et al.) against the corresponding experimental values. There is good agreement in the order of the shifts, although the predicted values have an overall spread of ∼8 ppm and an offset of ∼1.2 ppm. Figure also shows the 29Si chemical shifts calculated by DFT (again using the structure of van Konnigsveld et al. without optimization), and the agreement between calculation and experiment is very good. This example demonstrates that, when a high-quality crystal structure is available, eq can be used to provide at least an initial assignment, even when the structure contains many distinct Si sites.

Figure 7

Plots of 29Si δiso predicted by eq (red points) and calculated by CASTEP (blue points), against the experimental values[42] for SiO2-MFI. The gray line indicates the ideal 1:1 correspondence.

Morris et al. determined the structure of siliceous ferrierite (FER topology) and recorded high-resolution one- and two-dimensional 29Si NMR spectra of the material.[43] Five resonances were observed, corresponding to the five crystallographic Si sites, and these could be partially assigned using double-quantum correlation spectroscopy. The final two sites, Si4 and Si5, were assigned on the basis of a correlation between δiso and cos(⟨θSiOSi⟩)/cos(⟨θSiOSi⟩) – 1. The filled circles in Figure show the 29Si δiso predicted by eq for the experimental structure of SiO2 ferrierite. The experimental points (shown by crosses in Figure ) cover a smaller shift range than predicted, and agreement with calculation is poor. On closer inspection, the experimentally determined structure is likely to be unrealistic, with DI > 2.0 for four of the five Si sites. Morris et al. also optimized the structure using a force field method, leading to DI < 2.0 for four of the five Si sites. Despite this optimization, the MAE in the shifts predicted by eq (not shown) actually increases from 1.19 ppm for the experimental structure to 1.27 ppm after optimization. The open circles in Figure represent δiso calculated by CASTEP for the experimental structure, and it can be seen that eq predicts these well, even though agreement with the experimental shifts is poor. This confirms that the structures are likely to be unrealistic, rather than that eq cannot predict the values obtained by DFT. Upon optimization of the structure using CASTEP, the DI is reduced to below 1.5° for all five Si sites. From this optimized structure, CASTEP calculates values of δiso in excellent agreement with experiment (open squares in Figure ) and eq predicts very similar values (filled squares in Figure ), with a MAE of just 0.82 ppm (cf. the 0.92 ppm for the CASTEP values), although the order of the shifts for Si1 and Si5 is reversed. This example demonstrates that, even where a structure is an unrealistic representation of the material, eq is able to rapidly predict the outcome of the DFT calculation and, therefore, any large discrepancies between the experimental and predicted δiso most likely indicate that the structure must be improved.

Figure 8

Plots of 29Si δiso for the five Si sites in SiO2–FER. Experimental values (from Morris et al.[43]) are shown by crosses, values calculated by eq are shown by filled shapes, and empty shapes show values calculated by CASTEP. The values are calculated from the experimental structure (expt., circles) and the DFT-optimized structure (opt., squares).

Re-Examining the Literature Data

Using eq , it is possible to predict δiso from the crystallographic structures of the tectosilicates for which spectral data[7,8,14−16,19] was shown in Figure . Note that data for other classes of silicates were not considered here, as these contain Si with lower degrees of condensation. Where possible, the experimental crystallographic structures referenced in the original spectroscopic studies[48−72] were used here (see the Supporting Information (S6) for further details). Figure a shows a plot of the reported experimental 29Si δiso against that predicted by eq for all 31 tectosilicates (78 Si sites) discussed above. The points all lie reasonably close to the ideal line of 1:1 correspondence, although there is significant scatter, with a MAE of 4.7 ppm and a maximum deviation of 22.1 ppm. However, the greatest deviations are for the data reported by Newsam[15] (highlighted in red in the figure), which is unsurprising, since the experimental data were reported for Si(OAl)4 resonances, whereas eq inherently assumes Si(OSi)4 species. Figure b shows that, when these points are not included in the plot, much better agreement is now obtained, with a MAE of 2.5 ppm and a maximum deviation of 11.3 ppm. As demonstrated above, at least for pure silicate zeolites, structures or sites with DI > 2.0° are unlikely to represent an energetic minimum and can, therefore, be considered poor descriptions of the true local structure (for whatever practical reason). In the present data set, there are 15 SiO4 tetrahedra with DI > 2.0° (highlighted in blue in Figure a and b) and, when these are also removed from the plot, as shown in Figure c, the MAE drops to 1.5 ppm and the maximum deviation is now 6.2 ppm. While this MAE may not appear to be particularly low (certainly not within the <1 ppm accuracy required for interpreting some 29Si spectra of zeolites, as in the examples above), it is actually surprisingly small, given that the experimental structures include those determined for aluminosilicates where (in several cases) the Al sites in the framework were not located and the cations and water molecules in the pores (where present) were not considered in the chemical shift prediction. Furthermore, the structures were not optimized to an energy minimum (as would be carried out when the DFT-based prediction of accurate NMR parameters would be required[13,73]) and there are several cases of small (∼1 ppm) discrepancies between 29Si δiso values reported for the same Si site in the same material by different authors. Given the number of accumulated experimental errors present in the data set, it is, in fact, more remarkable that such a simple structure–spectrum relationship as eq can predict the experimental results so closely.

Figure 9

Plots of 29Si δiso predicted by eq against the experimental values for the tectosilicates shown in Figure .[7,8,14−16] (a) Plot including all reported data points, with the red points corresponding to Si(OAl)4 sites reported by Newsam.[15] (b) The same plot as part a but with the red points omitted. The blue points correspond to structures with a distortion index greater than 2.0°. (c) The same plot as part b but with the blue points omitted. For all parts, structural parameters were taken from the literature references cited in the original spectroscopic works, where possible.[48−72] Further details are given in the Supporting Information (S6). The gray lines indicate the ideal 1:1 correspondence.

The prediction of 29Si δiso from experimental crystal structures that may contain disordered framework substitution, extraframework cations, or water suggests that one such application of the work considered here might be in understanding the NMR spectra of real zeolites, where the disorder is too great to allow the application of meaningful DFT calculations. In such cases, the spectral resonances are generally broadened by this disorder and predicting isotropic shifts to <1 ppm accuracy is probably not required.

Conclusions

The use of empirical structure–spectrum relationships between 29Si δiso and the local bonding geometry around Si in zeolites is an area that has received intense interest from the 1970s until the beginning of the 21st century, when computing methods and hardware became sufficiently powerful to predict accurate NMR spectra from extended periodic crystal structures. However, there remain many structures and questions that DFT calculations are (currently at least) ill suited to handle—for example, where low amounts of Al occupy the tetrahedral sites in a zeolite, a series of large “supercell” calculations may be required to accurately model the distribution of Al ions within the material. In addition, atoms and molecular species such as Brønsted acidic H, water, and disordered (or dynamic) SDA cations may not be located (or, indeed, locatable) by diffraction experiments. In the most interesting case of modeling catalytic processes occurring within zeolites using molecular dynamics, there is a need to be able to provide a link between the thousands of structures generated per MD trajectory and experimental measurements including in situ NMR spectroscopy, which can (at least in principle) provide a rich variety of information on chemical species present, their concentrations, and any dynamics that may be present. In all of these cases, the need to be able to calculate NMR parameters to DFT-level accuracy is clear, but it is also evident that such calculations would be very time-consuming and not necessarily possible on routinely available computing hardware. In light of this, we re-examined the early empirical work that compared experimental NMR parameters with experimental structures (complete with experimental errors in both sets of data), using DFT calculations and more detailed statistical analysis to determine whether it is, indeed, possible to relate the local bonding geometry to the NMR spectrum in a simple way, or whether the disparate relationships reported in the literature were merely the result of chance fluctuations in the structures of the relatively small numbers of zeolites studied in any one case.

DFT calculations were first carried out on small Si(OSi(OH)3)4 clusters to model the immediate bonding environment around Si in a SiO2 zeolite. These clusters allowed ready systematic manipulation of the bonding geometry and revealed that both the mean Si–O bond length and the mean Si–O–Si bond angle have a strong influence on 29Si δiso. It was also clear from these calculations that the standard deviations of Si–O bond lengths and Si–O–Si bond angles influence δiso but to a lesser degree compared to the mean values. This approach was then applied to more realistic model microporous SiO2 frameworks, to investigate whether there were any additional longer-range effects arising from the extended periodic structure. We demonstrated that the relationships determined for the model clusters could be applied almost directly to the periodic frameworks, although the quadratic relationship between the mean Si–O bond length and δiso observed for the cluster compounds was revised to a simple linear relationship as, over the relevant range of bond lengths (i.e., 1.55–1.75 Å), x and x2 are essentially collinear. The final structure–spectrum relationship allowed the prediction to within ∼1 ppm of δiso calculated from DFT-level calculations with only knowledge of the Si–O bond lengths and Si–O–Si bond angles. The relationship was tested first on MFI- and FER-type SiO2 frameworks, and was able to match the order of the experimentally determined spectral assignment for many of the 24 Si sites in the MFI framework, allowing at least a preliminary assignment. Agreement with experiment was poorer for the FER-type material but improved to within 0.82 ppm upon structural optimization with DFT. These results demonstrate that our structure–spectrum relationship can accurately predict the DFT-calculated NMR parameters and, where significant disagreement is observed with the experimental spectrum, this may indicate that the crystal structure requires optimization. The relationship was also tested on published experimental crystal structures and NMR spectroscopic data for a range of tectosilicates and was able to predict the experimental δiso for Si(OSi)4 species to within 1.5 ppm (on average). However, the error was larger for Si(OAl)4 species owing to the known relationship between 29Si δiso and next-nearest neighbor Al/Si substitution. Our relationship was parametrized for SiO2 zeolites only and will require modification to take into account other cation substitutions.

To determine whether a given crystal structure represents a realistic energy minimum, the geometries of the set of model tectosilicates were optimized. This showed that the mean Si–O bond lengths converge to ∼1.62 Å and the distortion index is always below 2.0° for optimized structures. There was, however, no optimum value for the Si–O–Si angles. In other words, the SiO4 tetrahedron will be as close to ideal as possible (to within some tolerance dictated by crystal symmetry and framework topology), whereas the geometry of the connections between the tetrahedra (Si–O–Si linkages) is dictated by the long-range topology of the framework. By removing unrealistic SiO4 tetrahedra from the set of experimental structures and chemical shifts, the accuracy of the predictions was improved to 1.5 ppm, which is remarkable given the number of potential experimental errors in the structures and NMR data, and also the simplicity of the structural model used for the predictions.

This approach will never supersede DFT calculations, where such are possible. However, we envisage that the ability to predict shifts with close to DFT-level accuracy for systems where DFT is impractical or impossible will be a great advantage in providing a stronger link between experimental NMR spectroscopic measurements and structural and mechanistic models for a wide variety of experimental and computational studies.