Article pubs.acs.org/crystal
Hydrogen-Bond Analysis: Statistical and Computational versus Experimental Position Refinement Matteo Lusi,* Dawie de Villiers, and Catharine Esterhuysen* Department of Chemistry and Polymer Science, University of Stellenbosch, Stellenbosch 7600, South Africa S Supporting Information *
ABSTRACT: The relative accuracy with which O−H···O hydrogen-bonding parameters can be determined using the recently published polynomial neutron-normalization method is compared to that achieved using density functional optimizations at the M062X/6-31+G(d,p), ωB97X-D/aug-cc-pVDZ, and B3LYP/6-31+G(d,p) levels of theory. Calculations were repeated at the Hartree−Fock level utilizing the 6-31+G(d,p) basis set. Pairwise comparisons of the results obtained from these methods with the values measured in neutron diffraction experiments show that the computational and statistical methods are comparable under the conditions applied.
■
INTRODUCTION Hydrogen bonding is arguably the most exploited supramolecular interaction in chemistry and material science,1 while its presence is responsible for a number of chemical, biochemical, and physical phenomena such as the reactivity of protonated species,2 the structure of DNA,3 and the physical properties of water.4 More recently, hydrogen bonding has been utilized in material and pharmaceutical sciences with the aim of producing solids with specifically designed structures and properties.5,6 At the same time, deeper understanding of this interaction has enabled physical chemists to predict the periodic packing of crystalline materials with a high degree of success.7 Progress in these directions has been facilitated by the goldmine of information that is the Cambridge Structural Database (CSD),8 whose exponential growth can be explained by the widespread use of powerful in-house X-ray diffractometers and software that allow effortless structure solution. One obstacle in this development is the poor accuracy with which the positions of hydrogen atoms can be determined by X-ray diffraction (XD). In fact, this technique is not reliable for determining the positions of light elements and, due to the polarization effect of the H-bond, systematically yields an artificial shortening of the covalent bonding distance of hydrogen.9 The most promising procedure to tackle this limitation is deemed to be the quantum mechanical (QM) refinement of the hydrogen atoms’ coordinates. Unfortunately these calculations may be computationally expensive and require specialized knowledge that it is not always in the toolkit of the experimental crystallographer. Another approach is to use empirical adjustments that are based on the mean values observed in crystal structures determined by neutron diffraction (ND), which is sensitive to the positions of atomic nuclei and is not affected by bond polarization.10 This correction is commonly known as neutron normalization © 2014 American Chemical Society
(NN) and can be directly applied while retrieving crystal structures from the database. We recently surveyed the CSD in order to estimate how accurately NN corrects the hydrogen coordinates of H-bonded systems obtained using XD.11 On that occasion it was highlighted that commonly used NN values do not take bond polarization into account and therefore are not ideal to correct for the systematic error obtained with XD. An alternative polynomial normalization method (PN) was proposed that shows better correlation with experimental ND data.11 In this paper the PN method is compared to density functional theory (DFT) calculations, such that the H-bond descriptors (Figure 1) are simultaneously verified against experimental ND data.
Figure 1. Schematic representation of typical H-bond geometry descriptors.
■
RESULTS AND DISCUSSION In our previous work a list of 253 H-bonding interactions was compiled from the 87 structures retrieved from the CSD (version 5.31 February 2010), which were determined with error R < 0.10 at room temperature using ND and have at least Received: March 26, 2014 Revised: May 30, 2014 Published: June 6, 2014 3480
dx.doi.org/10.1021/cg500422p | Cryst. Growth Des. 2014, 14, 3480−3484
Crystal Growth & Design
Article
Figure 2. (a) Plot and linear fitting of values calculated for DH (polynomial normalization: green; M06-2X/6-31+G(d,p): red; ωB97XD/aug-ccpVTZ: yellow; B3LYP/6-31+G(d,p): pink; HF/6-31+G(d,p): blue) versus those observed by neutron diffraction; (b) plot of the error for each method as a function of the DA distance.
Figure 3. (a) Plot and linear fitting of values calculated for HA (polynomial normalization: green; M06-2X/6-31+G(d,p): red; ωB97XD/aug-ccpVTZ: yellow; B3LYP/6-31+G(d,p): pink; HF/6-31+G(d,p): blue) versus those observed by neutron diffraction; (b) plot of the error for each method as a function of the DA distance.
Information for computational details). H-bond descriptors were retrieved manually utilizing the X-Seed graphical interface13 (see Supporting Information for full results). The donor-hydrogen (DH) and hydrogen-acceptor (HA) distances calculated using the three DFT methods show almost identical linear relationships when compared to the values obtained by ND (Figures 2 and 3), while HF produces larger discrepancies. Furthermore, if one considers the Root Mean Square deviations of all the calculated distances from the ND values (Table 1) it is clear that the deviations are, in fact, quite small for the DH values obtained with most methods (the greatest deviations are observed for the HF calculations). In all cases the combination of intercept and angular coefficient indicate a sublinear relationship. This systematic deviation from normality can be verified by the Wilcoxon test (Table 2). The relatively high values of z confirm that there is a median difference between calculated and measured DH values, as opposed to a z-value of zero which would be expected for a normal distribution. The ranked test shows that B3LYP/631+G(d,p) is the best performing of the density functional
one O−H···O interaction. For the present work, 39 (≈ 15%) of those interactions have been selected by applying a randomizing function. For each selected H-bond a cluster of the interacting molecules was built, and the position of the H-bonding hydrogen was optimized at the M06-2X/6-31+G(d,p), ωB97XD/aug-cc-pVTZ, and B3LYP/6-31+G(d,p) levels of theory. The first two methods and basis sets were chosen as they have previously proven successful in describing hydrogen bonding12 and also accurately describe dispersion in noncovalent interactions. The B3LYP method, on the other hand, is a general-purpose method known to yield good results for hydrogen bonding, despite giving a poor description of dispersion effects, since it describes electron correlation at a local level only. In order to identify the role played by dispersion in the modeling of the O−H···O hydrogen bonding parameters, the calculations were repeated at the Hartree−Fock (HF) level utilizing the 6-31+G(d,p) basis set, since this method ignores electron correlation and is thus unable to describe dispersion effects. All the other atomic coordinates were kept fixed at the experimental positions (see Supporting 3481
dx.doi.org/10.1021/cg500422p | Cryst. Growth Des. 2014, 14, 3480−3484
Crystal Growth & Design
Article
sublinear relationship but yields better DH distances than the computational methods, as indicated by the coefficient of determination of 0.84. The Wilcoxon test for PN indicates a zvalue just above one standard deviation, only slightly worse than that obtained for B3LYP. Interestingly, the calculated values of the HA distances obtained using both the computational and the polynomial methods (Figure 3a) appear at first glance to fit the experimental values much better. Once again the slope of the interpolating line is lower than unity, but in this case the coefficients of determination vary between 0.83 for the values calculated at the HF/6-31+G(d,p) level and 0.99 for the M062X/6-31+G(d,p) values. The PN method exhibits a coefficient of 0.93 and an angular coefficient of 1.01. Nevertheless, the Root Mean Square deviations (Table 1) are higher for the HA distances than for the DH distances, indicating that the coefficients of determination are poor indicators of the accuracy of the method. Furthermore, the z-values obtained for the Wilcoxon test (Table 2) also suggest a poorer agreement between calculated and measured values for HA than for DH. The DFT methods exhibit similar levels of accuracy, reproducing the HA values within three standard deviations, with the B3LYP/6-31+G(d,p) level of theory again the best performing of all the methods. PN yields poorer values, showing significant differences between measured and calculated HA distances, although it is still more accurate than HF/ 6-31+G(d,p). The DHA angles are calculated equally well by all the computational methods [Figure 4; coefficients of determination between 0.79 with HF/6-31+G(d,p) and 0.91 with M06-2X/631+G(d,p) and z-values of −2.10 for B3LYP/6-31+G(d,p) to −2.77 for ωB97XD/aug-cc-pVTZ]; however they are extremely poorly described by the polynomial normalization method with a coefficient of determination of 0.30 and a z-value of −4.85. The Root Mean Square deviations confirm the poor fit, with an rms deviation of 14.64° for the DHA angle calculated using the PN method. The poor performance of the PN method can be attributed to the fact that the cosine law used to calculate the DHA angles from the DA, DH, and HA distances is extremely sensitive to small errors.
Table 1. Root Mean Square Deviations of DH (Å), H···A (Å), and D−H···A (deg) Values PN M06-2X/6-31+G(d,p) ωB97XD/aug-cc-pVTZ B3LYP/6-31+G(d,p) HF/6-31+G(d,p)
DH (Å)
HA (Å)
DH···A (deg)
0.02 0.03 0.04 0.03 0.06
0.07 0.04 0.06 0.04 0.07
14.64 3.64 3.97 2.33 4.46
Table 2. Wilcoxon Test: z-Values Calculated for DH (Å), H···A (Å), and D−H···A (deg) Values PN M06-2X/6-31+G(d,p) ωB97XD/aug-cc-pVTZ B3LYP/6-31+G(d,p) HF/6-31+G(d,p)
DH
HA
DH···A
−1.15 −2.12 −2.60 −1.11 −4.99
4.21 3.06 2.82 1.68 4.41
−4.85 −2.46 −2.77 −2.10 −2.70
theory (DFT) methods, able to reproduce calculated DH values to within one standard deviation of the measured value, while the other two DFT methods are slightly less accurate. The very large z-value (−4.99) calculated for HF indicates its poor performance, in agreement with the results shown in Figure 2a where it is clearly the method that exhibits the least linearity. This is further confirmed by plotting the errors between the calculated and measured values over the measured DA distance (Figure 2b). The source of the systematic deviation is clearly at small DA values, i.e., for strong Hbonds, where the DH distance is severely underestimated. This trend is similar, but not as extreme, for the other computational methods. The PN method exhibits a much smaller systematic error in the calculated DH distances, as evidenced by the even scatter of points shown in Figure 2b. The linear regression shown in Figure 2a confirms this, with a slope closer to 1.0 (0.76) and intercept nearer to 0.0 (0.23) than those obtained with the computational methods (the regression values of 1.0 and 0.0 would correspond to an exact agreement between calculated and measured methods). The PN method thus also shows a
Figure 4. (a) Plot and linear fitting of values calculated for DH···A (polynomial normalization: green; M06-2X/6-31+G(d,p): red; ωB97XD/aug-ccpVTZ: yellow; B3LYP/6-31+G(d,p): pink; HF/6-31+G(d,p): blue) versus those observed by neutron diffraction; (b) plot of the error for each method as a function of the DA distance. 3482
dx.doi.org/10.1021/cg500422p | Cryst. Growth Des. 2014, 14, 3480−3484
Crystal Growth & Design
Article
bond lengths and overestimation of HA distances, particularly for strong H-bonding interactions (i.e., at short DA distances). It is expected that using larger molecular clusters and exploring other levels of theory and basis sets would lead to more accurate H-bond descriptors, but again at greater computational expense and requiring a higher level of expertise. Furthermore, since qualitative analysis of H-bonding features in crystallography requires accurate hydrogen atom positions it is desirable to find an easy, fast, and reliable method for applying corrections to the hydrogen atom positions obtained from Xray diffraction data. A valid alternative to expensive QM calculations is represented by empirical normalization based on neutron data; however, the current implementation of this method in most software does not take bond polarization into account. Here we have shown that our previously published11 simple, fast polynomial correction yields DH and HA distances that are better than those obtained using the currently used empirical normalization method and comparable to (and in many cases also superior to) those calculated using QM methods. In cases where angular geometry is important, the new polynomial method is insufficient, and computational methods must be utilized to correctly find the H atom positions. It is clear, however, that the HF/6-31+G(d,p) level of theory poorly reproduces the measured hydrogen-bonding parameters and should be avoided. Instead, the results given here show that, surprisingly, the commonly used B3LYP method yields more accurate H-bonding parameters than those obtained with more modern DFT methods that have an improved description of the dispersion interaction. Since B3LYP is often more robust, exhibiting fewer convergence problems than M06-2X and ωB97XD, this method is to be recommended to experimental crystallographers who are novices in computational chemistry for the determination of corrected hydrogen-atom positions. On the other hand, in studies where only the DH and HA distances are important, or in the qualitative analysis of large data sets where DFT optimization would be too expensive, the new polynomial normalization method should be seen as a viable alternative, since it is extremely quick and easy to use and yields excellent values for these parameters. In addition, it holds the advantage of being largely insensitive to temperature, so it would be applicable for corrections to experimental data measured at all temperatures.
Returning to the data shown in Figure 2b, it is clear that all the computational methods tend to overestimate the short DH distances (corresponding to large DA interactions) and underestimate the longer DH values (short DA interactions). The performance of the computational methods is particularly poor for the longer DH distances involving very strong hydrogen bonds, with the largest deviations from the neutron data obtained for ETHDPH01 (underestimated by 10.6%, 11.5%, and 10.0% at the M06-2X/6-31+G(d,p), ωB97XD/augcc-pVTZ, and B3LYP/6-31+G(d,p) levels of theory, respectively). At the HF/6-31+G(d,p) level, the DH distance is underestimated by 13.5%, thus confirming that, although strong OH···O hydrogen bonds are primarily electrostatic in origin, dispersion nevertheless also plays an important role. For the shorter DH distances (i.e., O−H in weaker hydrogen bonds), the deviations are smaller, with the largest overestimation being for ABINOS01 (4.9%, 4.9%, and 5.2%, respectively), which also corresponds to the shortest OH bond investigated: 0.920 Å. The overestimation is smaller at the HF level (2.5%), suggesting that the DFT methods may be overemphasizing the effect of electron correlation on the OH bond at short distances. It should be pointed out that for both ETHDPH and ABINOS the hydrogen bonds are known to be sensitive to temperature.14 By default, the computational methods used here calculate the structures at absolute zero, so this temperature dependency is not sensibly modeled. On the other hand, we have shown previously11 that, from a statistical perspective, temperature has only a minor effect on the Hbonding parameters calculated with the PN method, and since it successfully yields accurate DH distances for ETHDPH and ABINOS this lends credence to the conclusion that the poor modeling of the ETHDPH and ABINOS H-bond parameters is due to a temperature effect. The maximum positive and negative deviations for PN are −6.1% for ASPARM03 and 4.6% for MGLUCP11, which are substantially less than those obtained with the computational methods. Similarly, the analysis of the HA error (Figure 3b) shows that the greatest deviations are also at short DA distances, but in the opposite direction to the DH distances (Figure 2b), since the calculated values are overestimated. This is to be expected, since as the HA distance shortens the DH distance will correspondingly lengthen for a fixed DA distance. The DH and HA distances show this inverse correlation (see Supporting Information); however, there is no relationship between the DH and HA distances and the DHA angle. The direction of the errors can also be seen from the Wilcoxon test z-values, which all indicate that the calculated DH distances are in general too low (negative deviations), whereas the HA values are too high (positive deviations).
■
ASSOCIATED CONTENT
S Supporting Information *
■
Computational details, a table containing all the H-bonded interactions, and figures showing the relationships between DH, HA, and DHA. This material is available free of charge via the Internet at http://pubs.acs.org.
CONCLUSIONS In summary, hydrogen bonding is of great importance for physical, chemical, and biological systems and requires a careful description that only neutron diffraction data can give. When these data are not available, structures determined from X-ray diffraction data must be corrected. Quantum mechanical refinement is the method of choice for this task; however, the quality of the results produced is proportionate to the sophistication of the model and ultimately to computational expertise and time. The methods utilized in this study are among the most popular used for studying hydrogen bonding, yet it is clear that the results obtained here do not address the problem at a satisfactory level due to underestimation of DH
■
AUTHOR INFORMATION
Corresponding Authors
*Tel: +27-21-808-3345. E-mail:
[email protected]. *Current address: Department of Chemical and Environmental Science, University of Limerick, Limerick, Ireland. E-mail:
[email protected]. Notes
The authors declare no competing financial interest. 3483
dx.doi.org/10.1021/cg500422p | Cryst. Growth Des. 2014, 14, 3480−3484
Crystal Growth & Design
■
Article
ACKNOWLEDGMENTS The authors thank the National Research Foundation (South Africa) and the Claude Leon Foundation for financial support of this work. The Centre for High Performance Computing of South Africa is thanked for computational time.
■
REFERENCES
(1) Gilli, G.; Gilli, P. The Nature of the Hydrogen Bond: Outline of a Comprehensive Hydrogen Bond Theory; Oxford University Press: New York, 2009. (2) Latimer, W. M.; Rodebush, W. H. J. Am. Chem. Soc. 1920, 42, 1419. (3) Watson, J. D.; Crick, F. H. C. Nature 1953, 171, 737. (4) Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381. (5) Aakeröy, C. B.; Seddon, K. R. Chem. Soc. Rev. 1993, 22, 397. (6) Burrows, A. In Supramolecular Assembly via Hydrogen Bonds I; Mingos, D. M., Ed.; Springer: Berlin Heidelberg, 2004; Vol. 108, p 55. (7) Engineering of Crystalline Materials Properties; Novoa, J. J.; Braga, D.; Addadi, L., Eds.; Springer: Dordrecht, 2008. (8) Allen, F. Acta Crystallogr. Sect. B 2002, 58, 380. (9) Allen, F. Acta Crystallogr. Sect. B 1986, 42, 515. (10) Allen, F. H.; Bruno, I. J. Acta Crystallogr. Sect. B 2010, 66, 380. (11) Lusi, M.; Barbour, L. J. Cryst. Growth Des. 2011, 11, 5515. (12) See, for example: Thanthiriwatte, K. S.; Hohenstein, E. G.; Burns, L. A.; Sherrill, C. D. J. Chem. Theory Comput. 2011, 7, 88−96. (13) (a) Barbour, L. J. J. Supramol. Chem. 2001, 1, 189−191. (b) Atwood, J. L.; Barbour, L. J. Cryst. Growth Des. 2003, 3, 3−8. (14) (a) Ananyev, I. V.; Barzilovich, P.Yu.; Lyssenko, K. A. Mendeleev Commun. 2012, 22, 242−244. (b) Jeffrey, G. A.; McMullan, R. K.; Takagi, S. Acta Crystallogr. 1980, B36, 373−377.
3484
dx.doi.org/10.1021/cg500422p | Cryst. Growth Des. 2014, 14, 3480−3484
