Article pubs.acs.org/JPCA
On the Anisotropy of van der Waals Atomic Radii of O, S, Se, F, Cl, Br, and I Hamed Eramian, Yong-Hui Tian,† Zach Fox,‡ Habtamu Z. Beneberu,§ and Miklos Kertesz* Department of Chemistry, Georgetown University, Washington, D.C. 20057-1227, United States S Supporting Information *
ABSTRACT: The Cambridge Structural Database (CSD) was used to obtain flattening factors to describe the overall anisotropy of nonbonding van der Waals (vdW) contacts between several main group elements. The method for obtaining the flattening factors is based on a novel minimization process. Results show that the vdW contact distances are significantly dependent on the environment and the orientations of the surrounding covalently bonded atoms: head-on vdW contacts are generally shorter than sideways contacts in overall agreement with earlier results by Nyburg and Faerman (Acta Crystallogr., Sect. B: Struct. Sci. 1985, 41, 274−279). With the exception of Se, we find flattening factors that are somewhat smaller than those found earlier. High-level ab initio quantum chemical calculations using Ar and Ne as a probe also confirm the flattening effect and its dependency on the environment. A dozen popular long-range corrected and dispersion supplemented density functionals are compared with the CCSD(T) data. While several of them perform quite poorly, four DFTD methods, especially B3LYP-GD3BJ, provided vdW flattening similar to those found by the CCSD(T) theory and experiment. others.9,10 Other definitions of vdW radii have been proposed and used in the literature indicating anisotropy.11−15 Batsanov pointed out that the changes of nonbonded contacts are often incorrectly attributed to hydrogen bonds.14 Badenhoop and Weinhold16 obtained anisotropic radii for a number of elements based on Hartree−Fock (HF) wave functions and their analysis of the steric repulsions with a testing He atom. They obtained in some cases, such as Cl, very good agreement with the anisotropy parameters of Nyburg and Faerman. In other cases, the agreement was poor: for F and O their anisotropy was much larger than the Nyburg−Faerman values; for S, it was much smaller. Similar trends were seen for selected molecules by Ikuta et al. based on the analysis of HF densities.17,18 Batsanov’s anisotropy values are significantly exaggerated for F, N, and O and to a lesser degree for Cl as compared to those of Nyburg and Faerman or the results of this work. A clear advantage of the probe methodology from the theoretical point of view is that it is not biased with respect to any models of halogen bonding19,20 or other models of secondary bonding.21 While clearly such effects can be significant, we focus here on the underlying anisotropy of the bare vdW interaction. According to Nyburg and Faerman, several elements common in organics (S, Se, Cl, Br, I) have a flattened atomic shape characterized by a smaller and a larger vdW diameter (Da and Db for each element) leading to angle-dependent vdW contact distances. A key parameter in characterizing the anisotropy is the flattening factor g of an ellipsoid defined by the lengths of the minor and major axis Da and Db (Da < Db, the
1. INTRODUCTION One of the most successful concepts of structural chemistry is that of the van der Waals radius.1 It posits that the shortest nonbonded contact distances DvdW(I,J) between atoms involved in the contact, I and J, can be well estimated by the sum of the van der Waals radii of the two elements, RvdW(I) and RvdW(J): DvdW (I, J) = R vdW(I) + R vdW(J)
(1)
Pauling’s original values for RvdW(I) have been extensively revised by Bondi.2,3 Bondi’s values are used throughout chemistry, with the exception of the subsequently revised value for hydrogen.4 Rowland−Taylor validated the rest of Bondi’s values in their extensive review of the crystallographic data of the Cambridge Structural Database (CSD)5,6 available in 1996. Most recently, Mantina et al. undertook a theoretical analysis of the van der Waals radii of all main group elements using high-level ab initio quantum chemistry including electron correlation to well describe the dispersion interactions.7 They used a probe atom methodology and found the Bondi− Rowland−Taylor values to be consistent. They used these values to scale the vdW radii for the rest of the main group elements and determined the vdW radii of 16 additional elements. It is clear that the Bondi−Rowland−Taylor values stood the test of time and will remain the pillars of structural analysis on nonbonded interactions. From the beginning, and especially in the book by Bondi, considerations appear3 regarding the basic assumption whether a spherical shape of vdW interactions is an appropriate approximation. Nyburg and Faerman have addressed this issue quantitatively in their analysis of the data available in the CSD at the time.8 van der Waals interactions are intrinsically somewhat anisotropic, more so for some elements than © 2013 American Chemical Society
Received: August 3, 2013 Revised: November 26, 2013 Published: November 27, 2013 14184
dx.doi.org/10.1021/jp4077728 | J. Phys. Chem. A 2013, 117, 14184−14190
The Journal of Physical Chemistry A
Article
We are using for the D versus μ correlation, angles that are close imposing the restriction
shorter contact corresponding to the head-on approach, with the approach angle μ=0): g = (1 − Da /D b) × 100%
|μ1 − μ2 | < 15°
(2)
Note, that g = 0 corresponds to an isotropic vdW shape. This work has two goals. First, we revisit the statistical analysis of the CSD for obtaining updated information on the anisotropy of vdW interactions for selected elements. We define an algorithm that leads to optimized flattening factors for each element making only one assumption: that the vdW shape is an ellipsoid of rotation. This is not the most general form, but there are not sufficient data in the CSD to obtain vdW distances in three directions. Our second goal supplements the CSD analysis with highlevel CCSD(T)22 quantum mechanical calculations using probe atom models and providing insights into the intrinsic nature of the van der Waals anisotropy. As a byproduct, we also obtain information on some popular versions of long-range corrected and dispersion supplemented density functionals.23−26
(3)
ensuring that we are analyzing data that do not mix two substantially different approach angles. The angles were averaged for the statistical analysis. The resulting data sets are sufficiently large for S, F, Br, Cl, I, and O. The data sets are more limited in size for Se and hence the resulting anisotropy value has a larger error. Note that we still have substantially more data points compared to those of Nyburg and Faerman given the growth of the CSD. For example, Nyburg and Faerman had approximately 32 data points for sulfur (D < 4.2 Å) while we have 396 data points for sulfur (D < 4.2 Å).
3. OPTIMIZATION OF THE FLATTENING FACTOR g FOR VAN DER WAALS INTERACTIONS The (D, μ) data points can be presented in a polar coordinate system, illustrated for the S···S contacts in Figure 2.
2. DATABASE SEARCH FOR VAN DER WAALS CONTACT ANISOTROPY INFORMATION CSD searches were set up to generate anisotropy information. The goal was to obtain a scatter plot of D(μ) expressing the distribution of contact distances as a function of the approach angle μ. Figure 1 illustrates the definition of the two independent approach angles μ1 and μ2 and the nonbonded contact distance D. μ is the smaller of μ1 and μ2.
Figure 1. (a) Defined parameters for X = sulfur, selenium, or oxygen. (b) Defined parameters for X = iodine, bromine, chlorine, or fluorine. μ1 and μ2 are the angles being monitored (μ2 is the larger angle). D is the nonbonded distance between the two X atoms.
Figure 2. S···S contact distances, D (measured from the origin) versus approach angle μ (measured relative to the vertical axis) as defined in Figure 1a obtained from the CSD as described in the text. The Bondi vdW sphere corresponds to a circle indicated by a continuous line in the graph. Axes are in angstroms.
The criteria for these searches were based on several factors that aimed at getting the best information on anisotropy. A coordination number of 1 was chosen for the X atoms in Figure 1 to obtain the effects of anisotropy in the most direct form. Also, the following filters were used in the CSD searches to maximize quality of data:27 “3D coordinates determined”, “R factor ≤0.05”, “not disordered”, “no errors”, “not polymeric”, “no ions”, “no powder structures”, and for O, F, Br, Cl, and I “only organics”. Two key issues have to be addressed regarding μ. One is that the data tend to taper off when μ approaches 90° because of steric crowding. Hence, we limited the search to 0° < μ < 90°. More serious is the issue of having two independent approach angles, μ1 and μ2, one for each side of the nonbonded interaction. This muddles the D versus μ correlation because a given contact D could belong to either or both μ1 and μ2. Nyburg and Faerman recognized this problem, but they found that μ1 and μ2 correlated somewhat with one another. They chose to use the smaller value only in their determination of Da and Db. Given the vast number of data available for at least some of our searches today, we choose a more restrictive route.
The innermost data points are related to the vdW contacts. Determining an ellipsoid that best describes these innermost points on the D versus angle plot is not trivial because this is not an ordinary fitting problem. Rather, we are interested in finding an ellipsoid that has the minimum number of points inside it while the two axes of the ellipsoid, Da and Db, are estimated. This optimization is meaningful only under a constraint. We perform an optimization while keeping the average D value constant. This constant is first chosen to be the same as the Bondi value, to be followed by a stability testing with respect to this value. The finding of the best Da and Db values is aided by a transformation whereby the data pairs of X = cos2(θ ), Y = [1/D(θ )]2
(4)
are plotted under these conditions: θ = (μ1 + μ2 )/2, |μ1 − μ2 | < 15° 14185
(5)
dx.doi.org/10.1021/jp4077728 | J. Phys. Chem. A 2013, 117, 14184−14190
The Journal of Physical Chemistry A
Article
In this transformed coordinate system, an ellipsoid becomes a straight line, and the condition that the average D value be unchanged is satisfied by requiring that the line intersect with the Bondi horizontal line at cos2(θ) = 0.5. Figure 3 illustrates
Figure 4. Number of data points above the line given by eq 6 as a function of the slope m based on the data presented in Figure 3. This graph is used to find the optimum slope m by choosing the lowest point on the graph. In this case D is constant and equal to D(Bondi, S···S).
Figure 3. Same data as in Figure 2 for S···S contacts replotted using transformed coordinates X and Y (see eq 4). The Bondi vdW sphere corresponds to a horizontal straight line. The green dashed line has the minimum number of data points above it while keeping the average D value constant (see text). θ is the approach angle as defined in eq 5.
we investigate the dependency of the optimized g value as a function of the average vdW contact D. Figure 5 shows the
the same data as given in Figure 2 together with the Bondi circle, which becomes a horizontal straight line (Bondi line). The data points above the Bondi line are not evenly distributed and tend to aggregate on the right (smaller μ values), expressing the anisotropy. Similar transformed data are shown in Figures S1−S7 of the Supporting Information for the other elements discussed in this paper. We optimize the slope of the line by minimizing the number of points above it while keeping the average constant and close to the Bondi value. The line is defined as Y = m(X − 0.5) + D−2
(6) Figure 5. Optimized g values as a function of D (in angstroms) in eq 6. The vertical line indicates the Bondi value, D(Bondi, S···S).
where D corresponds to the van der Waals contact distance (Bondi value, D(Bondi, S···S) = 3.60 Å in this case). We will return to the choice of D at the end of this section. Note that each m and D pair has a corresponding g value because Da(−2)
(−2)
=D
+ m /2, and
Db(−2)
(−2)
=D
− m/2.
optimized g values as a function of D. (The slope versus D diagram is similar, but g has more direct meaning.) Because g = 0 corresponds to an isotropic vdW shape, it is not surprising that as the contact D value is increased, any trace of anisotropy is lost and the statistics of the approach angle shows no preference; thus, the g value approaches 0. Similar optimized g values as a function of D are given in Figures S8−S14 of the Supporting Information for the other elements discussed in this paper. The variation of the anisotropy parameter is significant. However, most relevant are the values near the Bondi value of D, around 3.60 Å in this case. The number of data diminishes dramatically below 3.4 Å. The respective g values are provided here only for the discussion of the robustness of the approach. For longer distances beyond 3.7 Å the anisotropy disappears because at longer distances contacts cease to be sufficiently close and are not relevant to anisotropy. So while no unique value of g is found, a value of g = 15% with a spread of ±5% appears to describe the vast amount of data quite well for sulfur. Similar analysis has been performed for the models described in
(7)
Using D = D(Bondi, S···S), the resulting optimized line has a positive slope, corresponding to Da = 3.30 Å and Db = 3.99 Å, providing a flattening factor of g = 17.3%. This is a very significant deviation from a spherical shape. An isotropic interaction would have the minimum at a slope of zero. Figure 4 shows the number of points above the fitting line as a function of the slope m with a fixed D value. This is the key element of our optimization process yielding the sought after anisotropy by first obtaining the optimized slope in the transformed scatter plot where the m value at the minimum is chosen. The minimum is well-defined, but there is a range of m values that differ by only a few data points. From this minimum we determined the optimal slope m which in turn determined the optimal g values. While the process outlined is objective and clearly expresses strong anisotropy, we need to consider its robustness. How stable are the resulting optimal g values? Next 14186
dx.doi.org/10.1021/jp4077728 | J. Phys. Chem. A 2013, 117, 14184−14190
The Journal of Physical Chemistry A
Article
Figure 6. Models and geometry parameters used in the quantum mechanical calculations. Panels (a) and (b) define parameters for X = sulfur, selenium, or oxygen. Panels (c) and (d) define parameters for X = bromine, chlorine, or fluorine. μ is the angle being monitored. D is the nonbonded distance between the X and probe Ar (or Ne) atoms.
section Database Search for van der Waals Contact Anisotropy Information for the other elements and will be described and analyzed in Results and Discussion.
D3 dispersion correction, both with “zero damping” (GD(Zero)) and with BJ damping. We use the Gaussian 09 nomenclature: GD3 for the former and GD3BJ for the latter.31
4. QUANTUM CHEMICAL MODELING High-level ab initio quantum chemistry calculations and a variety of density functional theory (DFT) calculations were carried out on simple chemical models for obtaining anisotropy information on van der Waals interactions between an Ar (or Ne) atom serving as a probe and a small chemical species representing the vdW atom and its environment. To obtain a value for the vdW radius of the element of interest, the vdW radius of Ar (1.88 Å) was subtracted from the nonbonded equilibrium distance D. Similar results were obtained by using a Ne probe assuming a vdW radius of 1.54 Å, and the Ne probe data are presented in the Supporting Information in Tables S1 and S2. Ab initio calculations have been used successfully before to determine vdW parameters.6,28−30 The environments affect the vdW anisotropy but at the same time should indicate any intrinsic presence or absence of such anisotropy. Figure 6 shows the chemical models and the definition of the nonbonded distance D and the approach angle μ in each case. We have used the Gaussian 0931 software for the geometry optimizations for the individual molecules followed by rigid scans varying D while keeping the approach angle μ fixed.31 The following quantum mechanical methods have been used. Coupled cluster calculations were performed with the coupled clusters method using single and double excitations with triple excitations via perturbations (CCSD(T)).22,32,33 The aug-ccpVTZ and aug-cc-pVQZ basis sets were used.34,35 While these basis sets obviously do not allow us to obtain results at the complete basis set limit, they are sufficiently large to obtain the qualitative description of the vdW anisotropy that we are seeking here. A few DFT methods were also tested because of their popularity in solving problems involving vdW interactions. The B97D, ωB97, ωB97X, ωB97XD, B3LYP, BLYP, TPSSTPSS, PBEPBE, and M062X density functional methods were used with the aug-cc-pVTZ basis set in combination with an ultrafine grid as defined by the Gaussian 09 program. Some of these density functional methods were used with Grimme’s
5. RESULTS AND DISCUSSION The flattening parameters obtained by the slope optimization process outlined in Optimization of the Flattening Factor g for van der Waals Interactions are shown in Table 1. These Table 1. Observed Flattening Factors g (in %) for Select Group 16 and 17 Elements Based on the Analysis of CSD Structures element
g (this work)
g (Nyburg-Faerman)
F Cl Br I O S Se
a
6 11 16 17 0 21 21
