6156
J. Phys. Chem. B 2007, 111, 6156-6160
On the Performance of DFT and Interatomic Potentials in Predicting the Energetics of (Three-Membered Ring-Containing) Siliceous Materials Martijn A. Zwijnenburg,*,† Furio Cora´ ,†,‡ and Robert G. Bell†,‡ DaVy Faraday Research Laboratory, The Royal Institution of Great Britain, 21 Albemarle Street, London W1S 4BS, United Kingdom, and UniVersity College London, Department of Chemistry, Christopher Ingold Laboratories, 20 Gordon Street, London WC1H 0AJ, United Kingdom ReceiVed: February 7, 2007; In Final Form: March 30, 2007
We compare the enthalpies of transition for a range of SiO2 phases, including siliceous zeolites and dense polymorphs, calculated using three different interatomic potentials (Sanders-Leslie-Catlow (SLC), SastreGale (SG), van Beest-Kramers-van Santen (BKS)), and from B3LYP periodic DFT calculations, with the experimentally measured values. It is found that the calculated results show a linear correlation with the measured values but that they often either considerably underestimate or overestimate enthalpy differences compared to experiment. Care should thus be taken when comparing experimental and calculated results. A linear rescaling of the calculated enthalpies to put the data on the same energy scale is proposed. Furthermore, it is found that when comparing enthalpies of transitions for materials containing three-membered rings, for which there is no experimental data available, the values, rescaled to the experimental energy scale, are very similar for both DFT and interatomic potentials (except for the BKS potential). The latter result suggests that the energetics of three-membered ring containing materials is well described using both approaches. Finally, we discuss the transition enthalpies of four three-membered ring containing siliceous materials and demonstrate that three-membered ring containing materials are not necessarily energetically disadvantageous but do become so progressively with increasing number of three-membered rings.
Introduction In recent years there has been a great deal of computational work aimed at characterizing the properties (e.g., energetic stability, pore-volume, flexibility) of hypothetical zeolite frameworks, mostly in their siliceous forms.1-5 The main objectives have been to assess the likelihood of their successful (hydrothermal) synthesis and their application in novel processes. There are some particular difficulties associated with this research: the number of unique frameworks to optimize is very large (10 000+),4,6-10 making interatomic potential (IP) calculations the only tractable approach. However, a large part of the set of hypothetical frameworks comprises those which contain threemembered rings (3-MR), for which the reliability of IP methods is currently untested, in contrast to the case of materials with only larger rings (g4). The latter fact is due to the lack of experimental thermochemical data for siliceous versions of frameworks containing three-membered rings (which in turn is due to the fact that such frameworks typically have only been prepared with Al, Be, or Zn doping). Similar problems occur when attempting to rationalize the relationship between topology, geometry and energetics, and this has led us to exclude 3-MR containing frameworks in some of our previous work.11 Finally, a recent comparison of IP and DFT calculations on silica nanoclusters rich in small (three- and two-membered) rings,12 suggests that IPs may indeed struggle with describing such small rings. High-temperature lead borate solution calorimetry was performed by Navrotsky and co-workers on 13 all-silica zeolites * Corresponding author. † The Royal Institution of Great Britain. ‡ University College London.
(corresponding to the AFI, AST, CFI, CHA, FAU, FER, IFR, ITE, ISV, MEL, MTW, MWW, and STT framework types) and 3 dense silica polymorphs (Quartz, Cristobalite, and Tridymite) to obtain enthalpies of transition at 298 K.13 The latter quantity is the difference in enthalpy between a given siliceous material and R-quartz (normalized per mole of SiO2) at 298K and is indicative of the material’s thermodynamic stability at room temperature, as entropy differences between siliceous materials and R-quartz are known to be small and to span a narrow range (3.2-4.2 J (K-1 mol-1).14 The entropic contribution to the Gibbs free energy (the product of room temperature and the entropy differences) is thus at least an order of magnitude smaller than its enthalpic counterpart. The set of calculations which most closely mirrors the experiment would be a minimization of the Gibbs free energy of a siliceous material, to obtain the equilibrium structure at 298 K, followed by a single point enthalpy calculation, including the vibrational contribution to the internal energy, to obtain the material’s enthalpy at 298 K. Gibbs free energy minimizations, however, are computationally exceptionally expensive at the DFT level of theory. In the case of empirical interatomic potentials, the reliability of Gibbs free energy calculations rests on how the parameters have been derived. Most often the empirical potentials have been fitted to reproduce with static calculations (0 K, or more strictly athermal as zero-point motion is ignored) a set of experimental observables corresponding to high T measurements15-16 (for silica potentials typically those measured for R-quartz at room temperature).17-19 In such cases, adding thermal contributions to the enthalpy will result in double counting terms. A reliable Gibbs free energy calculation with IPs therefore requires the potential to be fitted to “0 K” structural data (atomic positions, elastic constants, dielectric constants etc)
10.1021/jp071060v CCC: $37.00 © 2007 American Chemical Society Published on Web 05/11/2007
Energetics of Siliceous Materials and to describe well the experimental phonon spectra, especially the low-energy low-frequency (rigid unit) modes (neither of which is generally the case (either due to lack of lowtemperature data or by design)). (Historically it was often seen as advantageous to have a potential which would reproduce room-temperature structures by means of static enthalpy minimizations.16) Also because, as noted above, experiments suggest that entropic effects are small, we thus chose instead to minimize the enthalpy of a structure at atmospheric pressure (or in the case of the DFT calculations the internal energy) and use its minimum value to calculate the theoretical enthalpy of transition. (The internal energy is a rather good approximation of the enthalpy at ambient pressure, as the PV term is very small compared to the internal energy differences considered (e.g., for Sanders-Leslie-Catlow optimized quartz 0.0022 kJ/mol, SLC optimized CHA 0.0039 kJ/mol and SLC optimized FAU 0.0045 kJ/mol compared to internal energy differences of 1-20 kJ/mol).) The reference temperature for this theoretical enthalpy of transition depends on the calculation method. DFT calculations are athermal (no temperature effects and zero-point motion included), so the reference temperature is approximately zero, while interatomic potentials have implicit temperature effects when fitted to non “0 K” data,15-16 in which case the reference temperature should be roughly equal to the temperature at which the experimental data were obtained. Similar calculations were previously performed by Henson20 (IP), Civalleri21 (DFT), Catti22 (DFT), and Astala et al.23 (DFT) for small subsets of the experimentally thermodynamically characterized silica polymorphs (respectively 8, 3, 3, and 3). In this communication we compare the enthalpies of transition predicted from the three most commonly used silica interatomic potentialssSanders-Leslie-Catlow (SLC), Sastre-Gale (SG), van Beest-Kramer-van Santen (BKS)sand from B3LYP periodic DFT calculations with the experimental values, and discuss the usage of linear rescaling to map the calculated enthalpies onto the experimental energy scale. Following on from this, we compare the rescaled enthalpies predicted by the interatomic potentials with those from DFT calculations for 3-MR containing materials. The DFT method is based on a number of assumptions; however the nature of these suggests that the DFT calculated rescaled enthalpies do not suffer from any major spurious dependence on the ring-size, as can occur in interatomic potentials which are fitted to a limited group of structures, often not more than one. The DFT data can thus be used as reference values in the absence of experimental data. Finally, we provide estimates of the enthalpies of transition for the siliceous realizations of three archetypal 3-MR containing zeolite framework types: MEI, OSO and RWY and one hypothetical framework: dt1_35,24 the only one of the nine frameworks corresponding to a uninodal simple tiling that has not yet been synthesized. Computational Methodology In the atomistic calculations, the interaction between silicon and oxygen ions was described respectively by the SLC,17-18 SG25 and BKS19 potentials. All initial structures of the 15 siliceous materials for which enthalpy of transition values are available and the three-ring containing frameworks studied were taken from the Atlas of Zeolite Framework Types26 and subsequently minimized without symmetry constraints using a constant pressure optimization algorithm as implemented in the program GULP27 (i.e., both atomic positions and cell parameters were optimized). A rational function optimizer (RFO)28 was employed to guarantee that true minima were found.
J. Phys. Chem. B, Vol. 111, No. 22, 2007 6157 Eleven siliceous materials (AFI, AST, CHA, FAU, FER, ISV, MEI, OSO, RWY, dt1_35, and Cristobalite) in addition to R-quartz were also optimized in periodic Density Functional calculations, employing the CRYSTAL03 code29 and the B3LYP hybrid exchange functional.30 The basis sets used are composed of localized atomic orbitals, expressed in a series of Gaussian type functions, and are of double valence plus polarization quality for each atomic species.31 All elements are treated at the all-electron level. We used a Monkhorst-Pack shrinking factor of 6 for R-quartz, reduced in other polymorphs according to the respective unit-cell sizes, and truncation thresholds of (7 7 7 7 14) for the Coulomb and exchange series, while SCF convergence thresholds were set to 10-7 Hartree for both eigenvalues and total energies. These tolerances ensure high numerical accuracy in the SCF calculations. The results of the geometry optimizations have been checked against the root-mean-square (rms) and absolute value of the largest component for both gradients and nuclear displacements. Optimization was considered complete when the four conditions were simultaneously satisfied for both fractional coordinates and unit cell parameters, using the values of 1.2 × 10-3 for the maximum gradient and 1.8 × 10-3 for the maximum displacement (all in au). All enthalpies quoted for the siliceous materials are, as discussed above, relative enthalpies compared to R-quartz (comparable to the experimentally measured enthalpies of transition), normalized to the number of T-atoms per unit cell. To rescale calculated enthalpies to the experimental energy scale we employed linear fits obtained through standard least-squares linear regression of the experimentally determined transition enthalpies to the calculated values. Results and Discussion Figure 1A-D shows the correlation between the experimental and calculated enthalpies of transition for the three interatomic potentials considered (15 structures) and the B3LYP periodic DFT calculations (7 structures). Just as previously observed by Henson et al.20 for the SLC potential and a more limited set of structures, we observe a good agreement between the experimental and calculated energetic ordering of the polymorphs, where the SLC and BKS calculated enthalpies are higher than experiment (considerably higher in the case of BKS) and the SG and B3LYP values lower. The correlation between calculated and experimental values is shown to be well described by a linear model (r2 SLC, 0.893; r2 SG, 0.866; r2 BKS, 0.889; r2 B3LYP, 0.964), where the relatively better fit for B3LYP is partly due to the smaller number of structures considered. When we rescale the calculated enthalpies to the experimental energy scale using these linear fits, the root-mean-square errors (RMSE) of the different methods (SLC, 1.2; SG, 1.4; BKS, 1.2; B3LYP, 0.8 kJ/mol) are small and in the same order of magnitude as the experimental 95% confidence interval for the measured enthalpies ((1-2 kJ/mol).13 For all the computational methodologies considered, the predicted energetic ordering of the polymorphs is slightly different from the experimental one, most notably for the SG potential which finds cristobalite and trydimite to lie lower in enthalpy than quartz (by 0.9 and 0.2 kJ/mol, respectively). Parts of these small differences (although probably not the relative misplacements of cristobalite, tridymite and quartz for the SG potential) are, as for the RMSE, likely to stem from the uncertainty in the experimental values, especially for such pairs of polymorphs where the experimental 95% confidence interval of the enthalpy considerably overlaps. Even though the non-free energy calculations are a rudimentary
6158 J. Phys. Chem. B, Vol. 111, No. 22, 2007 Zwijnenburg et al.
Figure 1. The experimental measured enthalpies of transition vs their calculated counterparts (in kJ/mol SiO2) for the three interatomic potentials considered (A, SLC; B, SG; C, BKS) and the B3LYP periodic DFT calculations (D). The dashed line in each figure is the line on which all points would lie if both methods agreed perfectly (i.e., the line for Hexpsc ) Hcalcsc).
Energetics of Siliceous Materials
J. Phys. Chem. B, Vol. 111, No. 22, 2007 6159
Figure 2. The rescaled IP enthalpies for the three-membered ring containing materials vs the rescaled B3LYP enthalpies (all in kJ/mol SiO2).
TABLE 1: Coefficients for the Linear Scaling Equations (Hexpsc ) a + bHcalcsc) for the Three Interatomic Potentials Considered and the B3LYP Periodic DFT Calculations SLC GS BKS B3LYP
a
b
r2
0.674 0.865 0.216 1.527
0.626 3.121 -0.342 -0.435
0.893 0.866 0.889 0.964
representation of the experiment, all the different methods thus do a good job of predicting the energetic ordering of the siliceous materials, although the absolute values of the calculated enthalpies of transition are either too small or too large. Caution is thus advisable when comparing experimental enthalpies of transition with calculated enthalpy values for hypothetical polymorphs, and rescaling to a common scale (e.g., the experimental energy scale or another reference scale) to ensure a fair comparison of values is almost always critical. The parameters of the linear rescaling equations which were fitted as part of this work are given in Table 1 but it is important to note that these will slightly differ for different choices of computational parameters (e.g., cut-offs, accuracies, grid-sizes, etc.). Figure 2 shows how the rescaled IP enthalpies for the 3-MR containing materials compare with the rescaled B3LYP enthalpies. For all three IPs the match is best at low enthalpies of transition (i.e., closest to the range of enthalpy to which the scaling equations were fitted) and becomes less good for larger enthalpies of transition (i.e., for less thermodynamically stable materials). Even so, the rescaled SLC and SG values follow the rescaled B3LYP values closely over the whole range considered, while only BKS gives substantially lower rescaled enthalpy of transitions for the less thermodynamically stable and 3-MR rich OSO and RWY materials. These results strongly suggest that both the SLC and SG IPs describe 3-MR containing materials as accurately as they do larger ring materials, without any spurious dependence on the ring-size. In contrast, the BKS IP, which performs comparably to the other IPs for the synthesized (more stable) materials, appears considerably to over-stabilize materials rich in three-membered rings. The BKS model should thus be used only cautiously for such systems. This includes systems where three-membered rings might be
TABLE 2: Rescaled Estimates of the Enthalpy of Transition of Three-Membered Ring Materials Hexpsc (kJ/mol SiO2) MEI dt1_35 OSO RWY
12.5-13.9 19.5-22.7 34.7-42.3 70.7-83.0
expected to arise during the course of the simulation, such as dynamics of amorphous and reconstructed silica surfaces. Finally, Table 2 gives predictions of the experimental enthalpy of transition of siliceous three-membered ring materials (as the range of rescaled enthalpies given by SLC, SG and B3LYP). Siliceous MEI is predicted to have an enthalpy of transition comparable to that of siliceous FAU, while the three other threemembered ring containing materials lie higher in enthalpy than the currently synthesized materials. Furthermore, the data suggests that there is a correlation between the number of threemembered rings (and the relative distance between them) in a material and the enthalpies of transition, as the latter increases in the order MEI (isolated three-membered rings) < dt1_35 (isolated three-membered rings but higher density than in MEI) < OSO (extended network of corner-sharing three-membered rings) < RWY (extended network in which each threemembered ring is part of a supertetrahedron). (When describing zeolite frameworks as so-called simple tilings the 3 memberedring concentration increases from 1/20 for MEI, via 2/7 for dt1_35 to 24/31 for RWY (the OSO framework can only be described as a more complicated non-simple tiling).) The dependence of the enthalpy of transition on the local environment of the three-membered rings in the material is in good agreement with recent work by Bromley and co-workers who observe a similar dependence for the energetic strain of small (three- and two-membered) rings in silica.32 Three-membered ring containing siliceous materials are thus not, as suggested previously,33 necessarily energetically much more unfavorable than those with only 4-membered and larger rings, but do become so progressively with increasing number and connectivity of three-membered rings. Conclusions Through a comparison of experimental thermochemical data and periodic DFT and interatomic potential calculations we
6160 J. Phys. Chem. B, Vol. 111, No. 22, 2007 demonstrate that the Sanders-Leslie-Catlow and Sastre-Gale interatomic potentials describe the energetics of siliceous threemembered ring containing materials as accurately as for those containing only larger rings and can thus be reliably used to assess the synthetic viability of hypothetical zeolite frameworks. For these two potentials, in contrast to the van Beest-Kramervan Santen potential, we do not observe any spurious effect due to the ring size. Furthermore, we observe that a rescaling of the calculated enthalpy data to the experimental energy scale is critical for a fair comparison of experimental and calculated enthalpies of transition. Finally, we estimated the transition enthalpies of four three-membered ring containing siliceous materials and demonstrate that three-membered ring containing materials are not necessarily energetically disadvantageous but do become so progressively with increasing number of threemembered rings. Acknowledgment. We like to thank Drs. D.W. Lewis and G. Sastre for organising the workshop that stimulated us to write this paper and we kindly acknowledge Drs. S.T. Bromley, G. Sastre, B. Slater, A.A. Sokol, S.M. Woodley and Prof. J.D. Gale for stimulating discussion. M.A.Z. acknowledges the European Commission for a Marie Curie Intra-European Fellowship (MEIF-CT-2005-010326). References and Notes (1) Foster, M. D.; Delgado Friedrichs, O.; Bell, R. G.; Almeida Paz, F. A.; Klinowski, J. Angew. Chem., Int. Ed. 2003, 42, 3896. (2) Foster, M. D.; Simperler, A.; Bell, R. G.; Delgado Friedrichs, O.; Almeida Paz, F. A.; Klinowski, J. Nat. Mater. 2004, 3, 234. (3) Foster, M. D.; Delgado, Friedrichs, O.; Bell, R. G.; Almeida Paz, F. A.; Klinowski, J. J. Am. Chem. Soc. 2004, 126, 9769. (4) Treacy, M. M. J.; Rivin, I.; Balkovsky, E.; Randall, K. H.; Foster, M. D. Microporous Mesoporous Mater. 2004, 74, 121. (5) Simperler, A.; Foster, M. D.; Delgado Friedrichs, O.; Bell, R. G.; Almeida Paz, F. A.; Klinowski, J. Acta Crystallogr. 2005, B61, 263. (6) Treacy, M. M. J.; Randall, K. H.; Rao, S.; Perry, J. A.; Chadi, D. J. Z. Kristallogr. 1997, 212, 768. (7) Delgado Friedrichs, O.; Dress, A. W. M.; Huson, D. H.; Klinowski, J.; Mackay, A. L. Nature 1999, 400, 644.
Zwijnenburg et al. (8) Boisen, M. B.; Gibbs, G. V.; O’ Keeffe, M.; Bartelmehs, K. L. Microporous Mesoporous Mater. 1999, 29, 219. (9) Woodley, S. M.; Catlow, C. R. A.; Batlle, P. D.; Gale, J. D. Chem. Commun. 2004, 22. (10) Earl, D. J.; Deem, M. W. Indus. Eng. Chem. Res. 2006, 45, 5449. (11) Zwijnenburg, M. A.; Bromley, S. T.; Foster, M. D.; Bell, R. G.; Jansen, J. C.; Maschmeyer, Th. Chem. Mater. 2004, 16, 3809. (12) Bromley, S. T. Phys. Status Solidi A 2006, 203, 1319. (13) Piccione, P. M.; Laberty, C.; Yang, S. Y.; Camblor, M. A.; Navrotsky, A.; Davis, M. E. J. Phys. Chem. B 2000, 104, 10001. (14) Piccione, P. M.; Woodfield, B. F.; Boerio-Goates, J.; Navrotsky, A.; Davis, M. E. J. Phys. Chem. B 2001, 105, 6025. (15) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J. Phys. Chem. 1996, 100, 7352. (16) Williams, D. E. J. Comp. Chem. 2001, 22, 1154. (17) Sanders, M. J.; Leslie, M. J.; Catlow, C. R. A. J. Chem. Soc. Chem. Com. 1984, 1271. (18) Schro¨der, K. P.; Sauer, J.; Leslie, M. J.; Catlow, C. R. A.; Thomas, J. M. Chem. Phys. Lett. 1992, 188, 320. (19) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. ReV. Lett. 1990, 64, 1955. (20) Henson, N. J.; Cheetham, A. K.; Gale, J. D. Chem. Mater. 1994, 6, 1647. (21) Civalleri, B.; Zicovich-Wilson, C. M.; Ugliengo, P.; Saunders, V. R.; Dovesi, R. Chem. Phys. Lett. 1998, 292, 394. (22) Catti, M.; Civalleri, B.; Ugliengo, P. J. Phys. Chem. B 2000, 104, 7259. (23) Astala, R.; Auerbach, S. M.; Monson, P. A. J. Phys. Chem. B 2004, 108, 9208. (24) Delgado-Friedrichs, O.; Huson, D. H. Discrete Comput. Geom. 2000, 24, 179. (25) Sastre, G.; Gale, J. D. Chem. Mater. 2003, 15, 1788. (26) Baerlocher, C.; Meier, W. M.; Olson, D. H. Atlas of Zeolite Framework Types; Elsevier: Amsterdam, 2001 (updates available at http:// www.iza-structure.org/). (27) Gale, J. D.; Rohl, A. L. Mol. Simul. 2003, 29, 291. (28) Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. J. Phys. Chem. 1985, 89, 52. (29) Saunders, V. R.; Dovesi, R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Harrison, N. M.; Doll, K.; Civalleri, B.; Bush, I. J.; D’Arco, Ph.; Llunell, M. CRYSTAL03 User’s Manual; Universita’ di Torino: Torino 2003. (30) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (31) Zicovich-Wilson, C. M.; Dovesi, R. J. Mol. Catal. A: Chem. 1997, 119, 449. (32) Bromley, S. T.; Moreira, I. P. R.; Illas, F.; Wojdel, J. C. Phys. ReV. B 2006, 73, 134202. (33) Earley, C. W. J. Phys. Chem. 1994, 98, 8693.
