Benzene–Hydrogen Bond (C6H6–HX

Feb 13, 2014 - Benzene−Hydrogen Bond (C6H6−HX) Interactions: The Influence of ... Moreover, the strength of the C6H6−HX interaction has been ana...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/JPCA

Benzene−Hydrogen Bond (C6H6−HX) Interactions: The Influence of the X Nature on their Strength and Anisotropy M. Albertí,*,† A. Aguilar,† F. Huarte-Larrañaga,† J. M. Lucas,† and F. Pirani‡ †

IQTCUB, Departament de Química Física, Universitat de Barcelona, Barcelona, Spain Dipartimento di Chimica, Università di Perugia, Perugia, Italy



ABSTRACT: The intermolecular potential energy of the C 6 H 6 −SH 2 and C 6 H 6 −NH 3 dimers is formulated as combination of independent electrostatic and nonelectrostatic contributions. The relevant parameters of the nonelectrostatic terms, derived from molecular polarizability components, have been proved to be useful to describe in a consistent way both size repulsion and dispersion attraction forces. The representation adopted for the electrostatic contribution asymptotically reproduces the dipole quadrupole interaction. To test the validity of the proposed potential formulation, the features of the most stable configurations of the systems predicted have been compared with the available ab initio and experimental data. Moreover, the strength of the C6H6−HX interaction has been analyzed comparing the obtained results with the corresponding ones for the C6H6−H2O and C6H6−CH4 systems, investigated previously with the same methodology. Information on the relative orientation dependence of the partners, arising from the anisotropy of the intermolecular interaction, evaluated at different intermolecular distances, has been also obtained. Such information is crucial to evaluate sterodynamics effects in bimolecular collisions. in proteins.8−12 Studies on XH−π systems reveal that molecules containing π-electrons can participate easily in the formation of XH−π bonds, whose strength depends on the nature of X. In particular, the OH−π and the NH−π interactions are stronger than the CH−π ones.1 Actually, the binding in the CH−π system lies between the weakest class of hydrogen and van der Waals bonds.13 This feature can be attributed to the low polar character of the CH bond in several hydrocarbon molecules. When electrostatic effects are negligible, dispersion forces provide an important contribution to the total interaction potential energy. In the case of weak and moderate hydrogen bond interactions, the assessment of the relative role of the different interaction components is not trivial. In the last years the interaction energies were often calculated using the symmetry adapted perturbation theory (SAPT)14−16 and SAPT(DFT)17,18 methodologies, which allow a proper partition of the long-range potential for specific configurations into basic components evaluated at various levels of perturbation theory. Therefore, the detailed computational investigation of small prototype systems, can be very instructive and crucial for the modeling of systems at increasing complexity. Usually, the weak interactions as those involved in the present systems need to be calculated using very large basis sets (even for small systems, the full potential energy surface depends on 3n − 6 variables, with n being the number of

1. INTRODUCTION Noncovalent interactions, dominant at large and intermediate intermolecular distances, affect both the static properties of weakly interacting partners and the stereodynamics of several phenomena. Such interactions are important in many fields of chemistry, physics, and biology since they determine structures and reorientation dynamics, as well as basic properties of liquids, molecular crystals, and biological systems.1−4 Crucial aspects to be investigated regard the nature of the intermolecular bond, the dynamics of collisions under a variety of conditions, the formation of weakly bound aggregates, as well as the growth and the structure of solvation shells. Among all intermolecular interactions, the case of hydrogen bonds has a paradigmatic importance. They are ubiquitous in chemistry, biology, and material science5,6 and can be grouped into three categories (weak, moderate, and strong)7 and different interaction components, as electrostatics, polarization and charge transfer, combine to define their binding energies. The relative contributions of such components vary widely among the different categories of hydrogen bonds. Although great attention has been focused on moderate hydrogen bonds with binding energies between 4 and 15 kcal mol−1,8 another type of hydrogen bond, involving aromatic π-systems, has awakened interest, mainly due to their abundance in biological systems. In particular, many experimental and theoretical studies reporting details on the XH−π (X = C, O, N, S) systems can be found in the literature. CH−π interactions are able to stabilize adducts between aromatic residues (see for instance ref 5.) and statistical analyses of protein databanks have shown that OH−π, NH−π, and SH−π interactions are frequently present © 2014 American Chemical Society

Received: November 6, 2013 Revised: February 13, 2014 Published: February 13, 2014 1651

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

component (Vnel).34 Such formulation, exploiting the balancing of positive and negative contributions, properly includes the combination of size repulsion, dispersion and induction attraction, mixed terms, and damping effects. This characteristic has been properly tested in ref 35, where a comparison with the SAPT energy decomposition was performed. Therefore, the so parametrized Vnel can be added to the electrostatic (Vel) component, which can be represented under different ways.36 In the present study, we have used a point charge distribution to maintain a compact representation of the interaction and to make a proper comparison with other systems previously analyzed in the same way. In the present study, we are interested in the semiempirical formulation of the C6H6−SH2 and C6H6−NH3 dimers. In order to test the accuracy of the interaction model for these systems, equilibrium geometries and energies are compared with available ab initio data. The ultimate purpose, however, concerns a complete investigation of the XH-C6H6 anisotropic interaction as a function of the nature of X and the evaluation of the energy barriers controlling the stereodynamics of collisions under a variety of conditions. The results obtained for C6H6−SH2 and C6H6−NH3 are also compared with those predicted also by our model and previously published for the C6H6−H2O37 and C6H6−CH438 systems. The paper is structured as follows: in section 2, we outline the formulation of the semiempirical PES, in section 3, the model predictions for C6H6−SH2 and C6H6−NH3 are presented. In section 4, the XH−C6H6 interaction and its anisotropy are analyzed, and concluding remarks are given in section 5.

atoms) and ab initio calculations can only be performed for selected configurations. For instance, Tsuzuki et al.1,13,19,20 investigated the C6H6−CH4 interaction energy and equilibrium geometry using different basis sets at selected geometries. These detailed studies demonstrated that a highly computationally demanding CCSD(T) level of calculation (using a large basis set near saturation) is required in order to obtain a reliable description of the system. Large basis sets have been also exploited to investigate the origin of the attraction and the directionality of the NH−π interaction.21 Selected geometries have also been considered to compare the performance of several ab initio methods in describing the C6H6−H2O dimer binding energy curve.22 Driven by the high computational cost of CCSD(T) level of calculation, required to obtain a good description of the dispersion contribution, Deligkaris and Rodriguez23 have proposed the correction of DFT interaction energies by introducing an empirical dispersion term valid for a defined range of intermolecular distances. Sherrill et al.24 use large, correlation-consistent basis sets to properly describe the potential energy curves for C6H6−SH2, C6H6−CH4, and C6H6−C6H6 (the later accurately investigated at SAPT level).25 In that work, the performance of several methods is tested in the calculation of nonbonded interactions, comparing their computational cost. High quality ab initio calculations are indeed of fundamental importance, providing valuable information on strength and nature of weak intermolecular bonds. However, if such an effort is required to characterize energy and structure of the most stable configuration, much more difficulties are met in the characterization of the less stable configurations. On the other hand, a good description of the interaction cannot be limited only to the geometry and energetic of some stable configurations. A reliable formulation of the whole potential energy surface (PES) is absolutely necessary in molecular dynamics (MD) simulations. In particular, reorientation effects promoted by the interaction anisotropy could play a crucial role on the dynamical evolution of the approaching and removing partners. Often, in MD simulations of a complex system, the involved PES is constructed by considering different groups of the interaction components, which are conveniently assembled to provide the description of the global interaction (force field). In this case, not only an accurate description of the system at equilibrium, derived from accurate ab initio calculations and/or high level experimental data,26 is needed but also the adoption of potential energy functions able to properly describe the energetic behavior of the system in the whole configuration space.27 In the last years, some of us have proposed a modification of the Lennard-Jones potential energy, the Improved LennardJones (ILJ) function,26,28 useful to describe the combination of some effective interaction components in neutral and ionic systems (see for instance refs 29−33). The basic parameters involved can be anticipated by correlation formulas given in terms of fundamental physical properties of the interacting partners, such as the electronic polarizability α (which can be properly decomposed for molecules in bond and/or effective atomic components) and the charge distribution over the molecular frame, determining total molecular charge and permanent multipole. The ILJ function and the related correlation formulas allow the description of systems at increasing complexity. In particular, the formulation adopted and the involved potential parameters allow to represent in a suitable and compact formulation the nonelectrostatic

2. POTENTIAL ENERGY SURFACE The total intermolecular potential energy (Vtotal) for the C6H6− SH2 and C6H6−NH3 systems has been modeled by assuming the separability in electrostatic, Vel, and nonelectrostatic, Vnel, partial contributions (Vtotal = Vel + Vnel). The electrostatic component, Vel, has been calculated by applying the Coulomb law to any pair of point charges placed on different molecules. The charge distribution on the SH2 and NH3 molecules, with point charges placed on the atoms, are compatible with the geometry and the dipole moment of the molecules calculated at MP2FC/6-311G* level.39 In particular, for the SH2 molecule a charge of −0.324 au has been placed on the S atom and charges of 0.162 au have been placed on the H atoms, while for NH3, charges of −0.996 au and 0.332 au have been placed on the N and H atoms, respectively. On the other hand, as in our previous studies on ion-C6H6 systems (see for instance ref 40), a total of 18 point charges have been distributed on the C6H6 molecule frame (six placed on the H atoms and the remaining 12 at a fixed distances from C atoms on both sides of the aromatic ring). Such distribution has been chosen from the consideration that, asymptotically, Vel must correspond to the ion quadrupole interaction.40 This leads to a charge of +0.09245 au on each H atom and to two negative charges of −0.04623 au separated by 1.905 Å on each C atom, placed on opposite sides of the benzene plane. Following the same procedure previously adopted to study the C6H6−H2O37 and C6H6−CH438 systems, the non electrostatic component of the interaction energy, Vnel, has been decomposed into 12 molecule-bond interaction terms. The intermolecular potential representation as a sum of terms, each one related to the pair interaction between specific centers can be justified at theoretical level considering distributed molecular response properties36 and the decomposition of the polar1652

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

izability in local contributions.41−43 To this end, the C6H6 molecular polarizability (αC6H6 = 10.32 Å3)44 has been decomposed in bond polarizabilities (αCC and αCH),41,45 each one having parallel and perpendicular components related with the dimension and the shape (often ellipsoidal) of the electronic charge distribution around the considered bond.41 On the contrary, the polarizability of both SH2 and NH3 (αSH2= 3.98 Å3 and αNH3= 2.16 Å3),46 due to their small value in comparison with that of C6H6, has not been decomposed and only the corresponding average value, centered on the S and N atoms, has been considered. Accordingly, the C6H6-M (where M is here SH2 or NH3) interactions are expressed as a sum of molecule-bond (M-CC and M-CH) terms. Thus, the overall potential energy, Vtotal, for a C6H6-M system is expressed as 6

Vtotal = Vel +

Table 1. They have been obtained using the molecular and bond polarizability components. The procedure adopted is described in detail in ref 51 (see also ref 40). Table 1. Perpendicular and Parallel Components of the Well Depth (ε⊥, ε∥) and of the Equilibrium Distances (r0⊥, r0∥) for the Different Molecule-Bond Pairsa

a

i

i=1

i

i=1

(1)

The VM−CC and VM−CH terms are represented by means of the ILJ function26,47 as

(2)

where r corresponds to the distance from the S or N atom to the center of a bond and γ is the angle between the r vector and the bond axis. The m parameter is taken equal to 6, the typical value for neutral−neutral interactions. The well depth, ε(γ), and the equilibrium distance, r0, for each interaction pair are modulated from the corresponding perpendicular (ε⊥, r0⊥) and parallel (ε∥, r0∥) values ε(γ ) = ε⊥ sin 2 γ + ε cos2 γ

(3)

r0 = r0 ⊥ sin 2 γ + r0 cos2 γ

(4)

which are derived by exploiting the values of αSH2 or αNH3 with the parallel and perpendicular components of αCC and αCH.48 The angular dependence allows us to take into account the variation of the well depth and the equilibrium distance for different approaches of the molecule (SH2 or NH3), considered as a “pseudoatom”, to a particular bond (ε(γ) and r0(γ)). The first term in eq 2 (positive) represents the size-repulsion contribution associated with each molecule-bond pair, while the second one (negative) provides the effective dispersion plus induction attraction. The n(r,γ) exponent, defining simultaneously the falloff of the molecule-bond repulsion and the strength of the attraction is expressed as ⎛ r ⎞2 n(r , γ ) = β + 4.0⎜ ⎟ ⎝ r0(γ ) ⎠

ε∥/meV

r0⊥/Å

r0∥/Å

4.981 5.481 4.142 4.906

6.868 4.720 5.379 4.110

4.139 3.947 3.954 3.731

4.404 4.147 4.249 3.938

1 meV = 0.02306 kcal mol−1

3. RESULTS In this section, the potential energy model described in the previous section is used to calculate the interaction energy profiles as the molecule (SH2 or NH3) and the benzene approach under selected configurations. The energy minima along these interaction curves are identified and characterized (geometries) comparing these data with ab initio and experimental data available in the literature. Subsequently, this static study is complemented by performing MD simulations of both C6H6−SH2 and C6H6−NH3 systems, by allowing SH2 and NH3 to explore the whole phase space around the benzene molecule. The molecules have been kept rigid throughout this study. The geometry adopted for C6H6 has been the same used in our previous studies (see for instance ref 38), with D6h symmetry and the CC and CH distances equal to 1.390 and 1.090 Å, respectively. On the other hand, for SH2 and NH3 we used geometries obtained from ab initio optimizations employing a MP2FC/6-311G* level of theory, which report for dihydrogen sulfide a SH distance of 1.341 Å and a HSH angle of 93.437°.39 At the same level of accuracy, the ammonia molecule exhibits a NH distance of 1.0105 Å and a HNH angle of 107.274°.39 In order to compare more directly our static results with the similar ones found in the literature, the potential energy values are given in kcal mol−1. 3.1. Energy and Geometry Predictions for C6H6−SH2. In order to obtain a general idea of the shape of the PES and the possible equilibrium geometries we have designed seven possible approaches between the SH2 and C6H6 molecules. In all the cases, the relative orientation of both fragments is kept fixed and only the distance between the S atom and the benzene center of mass (c.m.) is varied. The charge distribution of C6H6 favors out of plane approaches, especially with SH2 situated along the C6 axis of C6H6 (axial approaches), with the two H atoms placed at shorter distances from the aromatic ring than S. The potential energy curves for the selected configurations of the system are shown in Figure 1 as a function of the distance from the S atom to the C6H6 c.m. Panels a and b correspond to perpendicular approaches of SH2 to the aromatic plane, while panels c and d describe the on plane ones (equatorial approaches). More in detail, in panel a,

⎡ ⎛ r0(γ ) ⎞n(r , γ ) m VILJ(r , γ ) = ε(γ )⎢ ⎜ ⎟ ⎢⎣ n(r , γ ) − m ⎝ r ⎠ m⎤ n(r , γ ) ⎛ r0(γ ) ⎞ ⎥ − ⎜ ⎟ n(r , γ ) − m ⎝ r ⎠ ⎥⎦

ε⊥/meV

In the present study, in order to investigate the capability of the model to describe the C6H6−SH2 and the C6H6−NH3 intermolecular interactions through the molecular polarizability, the same value of the β parameter used to study C6H6−CH4, equal to 8.75, has been considered.

6

∑ VM − (CC) + ∑ VM − (CH)

M-bond SH2−CC SH2−CH NH3−CC NH3−CH

(5)

where β is an adjustable parameter related to the hardness of the interacting partners, which adds flexibility to the potential energy function in comparison to the Lennard-Jones one.26 In particular, the decrease of β can indirectly include stabilization effects due to a charge transfer in the perturbative limit.27,35,49,50 All the parameters of the ILJ function describing the C6H6−SH2 and the C6H6−NH3 interactions are given in 1653

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

Figure 1. Potential energy curves for the selected approaches of SH2 to the benzene ring indicated in the figure as a function of the distance from S to the C6H6 c.m.

Table 2. Intermolecular Interaction Energy (Vtotal) at Equilibrium and the Corresponding Distance from S Atom to the Center of C6H6 (RS−C6H6) for Approaches of SH2 with a HS Bond on the Rotational Axis of C6H6 (Cs Symmetry) and for Approaches in Which the C2 Axis of SH2 and the C6 One of C6H6 Are Coincident (C2v Symmetry) symmetry group

Vtotal/kcal mol−1

RS−C6H6/Å

methodology

Cs C2v C2v C2v C2v C2v C2v

−2.71 −2.72 −2.64 −2.74 −2.81 −2.83 −2.85

3.693 3.675 3.80 3.80 − 3.80 3.80

present present CCSD(T)/aug-cc-pVTZ52 CCSD(T)/aug-cc-pVQZ52 CBS CCSD(T)52 CBS CCSD(T)24 CCSD(T)/aug-cc-pVQZ//CCSD(T)/CBS53

from the literature are also included for comparison. As can be seen, the predictions of the model are in quite good agreement with available ab initio data. The overall similarity of the plots in parts a and b of Figure 1 suggests that the minima should be connected through low rotational barriers. In order to check this feature, we have explored the rotation of SH2 in the yz plane being C6H6 confined in the xy plane maintaining the axial configuration and the Bz-SH2 intermolecular distance close to the equilibrium geometry (3.680 Å). The rotation angle ϕ has been taken as the one formed by one SH bond and the C6 rotational axis of C6H6. The PES values are represented in Figure 2 as a function of the rotation showing a minimum energy of −2.804 kcal mol−1 at a rotational angle of 17°. Identical results are obtained when SH2 rotates in the xz plane. As expected, small variations on the potential energy are observed when ϕ varies from 0 to 90°. These results are indicative of the flatness of the PES for axial approaches.

the interaction energy is illustrated for two approaches in which the C2v symmetry axis of SH2 coincides with the C6 symmetry axis of C6H6. In one of them, the H atoms of SH2 are aligned with a HCCH axis of C6H6, while in the second one the H atoms of SH2 point toward the center of two CC parallel bonds. Panel b in Figure 1 contains the potential energy curve corresponding to an axial approach with one of the SH bonds of the SH2 molecule pointing along the C6 symmetry axis of C6H6. As regarding equatorial approaches, panel c contains the values of the potential energy for two coplanar approaches of SH2 to C6H6 in which the S atom of the SH2 molecule points toward the center of a CC bond. On the contrary, in panel d, the S atom point toward an H atom of C6H6. As it can be seen, perpendicular approaches are energetically favored with respect to coplanar ones. Moreover, the minimum potential energy values associated with the different axial approaches, with one or both H atoms pointing to the aromatic ring, are very similar indicating that, for these approaches, the PES is very flat. In particular, observing the two curves in Figure 1a, it is evident that the potential energy values are nearly independent of the relative position of the H atoms of SH2 in the bidentate conformation, being the differences in energy of about 0.01 kcal mol−1. This result is in agreement with the previous work of Tauer et al.,52 which found a negligible barrier for the twisting of SH2 above C6H6. This interaction feature allows the SH2 molecule to rotate freely around its C2v axis. Clearly, minimum energy equilibrium geometries can be identified from the curves shown in Figure 1, parts a and b, one with C2v symmetry (a) and another one with Cs symmetry. The total interaction energy, Vtotal and the internuclear distance R at equilibrium for these adducts are given in Table 2, where some results available

Figure 2. Potential energy represented as a function of the rotation of SH2 in the yz plane at a fixed value of RS−C6H6 = 3.680 Å. For ϕ = 0°, one SH bond is aligned along the C6 rotational axis of C6H6. 1654

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

Figure 3. Potential energy curves for the selected approaches of NH3 to benzene indicated in the figure as a function of the distance from the N atom to the C6H6 c.m.

Table 3. Interaction Energy (Vtotal) and Distance from N Atom to the Center of C6H6 (RN−Bz) for Axial Approachesa complex

Vtotal/kcal mol−1

RN−C6H6/Å

methodology

tridentate tridentate monodentate monodentate monodentate monodentate monodentate monodentate

−2.17 −1.70 −2.32 −2.37 −2.22 −2.37 −2.43 −2.47

3.506 3.6 3.505 3.6 3.6 − 3.60 3.60

present MP2/cc-pVQZ21 present MP2/cc-pVQZ21 ECCSD(T)limit (cc-pVDZ)21 CCSD(T)/AVTZ55 CCSD(T)/aug-cc-pVTZ53 CCSD(T)/aug-cc-pVQZ53

a Tridentate complex refers to a confirmation with both molecules symmetry axis aligned and monodentate to the configuration with one NH aligned with the benzene C6 axis (see text).

3.2. Energy and Geometry Predictions for C6H6−NH3. In order to characterize the C6H6−NH3 PES, we have proceeded in an analogous way to what has just been explained in the case of SH2, building potential energy curves for intermolecular approaches at fixed relative configurations. The potential energy curves for some selected configurations of the C6H6−NH3 complex are shown in Figure 3. Panels a and b in Figure 3 represent potential energy curves corresponding to two axial approaches of NH3 perpendicular to the benzene plane, while the on plane (equatorial) ones are shown in panels c and d. Similarly to the case of dihydrogen sulfide, axial approaches present a larger interaction energy and more stable minima than equatorial approaches. However, these differences are less noticeable in the case of ammonia. The representations in panel a correspond to approaches of NH3 in which the C3 symmetry axis of the molecule coincides with the C6 one of C6H6 and, therefore, the three hydrogen atoms are equidistant to the benzene plane (tridentate configuration), while in panel b) is described an approach of NH3 in which a NH bond coincides with the C6 symmetry axis of C6H6 (monodentate configuration). On the other hand, in panel c are represented the values of the potential energy for two equatorial approaches of NH3 to C6H6, with the N atom of NH3 pointing toward an H atom of C6H6, while in panel d, the N atom of NH3 point toward the center of a CC bond of the aromatic molecule. As it can be seen in Figure 3a, the potential energy curves for the two axial configurations with NH3 differently oriented above the benzene plane are almost identical, reflecting again the flatness of the PES with respect to rotation of the ammonia molecule around its C3v axis. These results are in agreement with the high-resolution optical and

microwave spectra of C6H6−NH3 in the gas phase (the appearance of two peaks suggests that the ammonia can freely rotate in C6H6−NH3).54 Moreover, they also agree with the high-level ab initio calculations carried out by Tsuzuki et al.21 The results of the intermolecular interaction energy Vtotal and the intermolecular distance R corresponding to the equilibrium structures for the tridentate and monodentate complexes represented in panels a and b of Figure 3, respectively, are listed in Table 3, where some results, available from the literature, are included for comparison. Our potential energy model predicts that the monodentate complex is slightly more stable than the bidentate (not shown) and tridentate ones. However, the model predicts smaller differences of the binding energy for the different geometries of the adducts than those obtained from ab initio methods. In particular, our potential energy predictions at the equilibrium geometry for the mono and tridentate complexes differ only for 0.25 kcal mol−1, whereas this energy difference is increased to 0.67 kcal mol−1 by MP2 calculations using a large basis set.21 Furthermore, the C6H6−NH3 PES has been investigated more in detail exploring the rotation of NH3 around C6H6. The rotation angle ϕ has been taken as the one formed by one NH bond and the C6 rotational axis of C6H6. In Figure 4 the interaction energy is represented as a function of ϕ at a fixed distance from N to the C6H6 c.m. of 3.505 Å. Results indicate that the most stable geometry of the C6H6− NH3 is the one in which one of the NH bonds (the closest to the aromatic plane) forms an angle, ϕ, of about 18° with the C6 rotational axis of C6H6, being he minimum of the potential energy equal to −2.491 kcal mol−1. 1655

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

Figure 4. Potential energy represented as a function of the rotation of NH3 at a fixed value of the distance from N to the center of mass of C6H6 (see text). ϕ = 0° corresponds to a monodentate configuration.

3.3. Molecular Dynamics Simulations. The static results of the previous sections indicate that the most stable structures correspond to axial configurations in which neither symmetry axis of the molecules nor the bond are completely aligned with the C6 rotational axis of C6H6. The preferred relative position of both SH2 and NH3 respect to C6H6 is further investigated by performing molecular dynamics (MD) simulations of the C6H6−SH2 and C6H6−NH3 dimers in gas phase using the DL−POLY program.56 Bearing in mind that at low temperatures, if the initial configuration of the system is not far from equilibrium, the distribution of distances, angles and configuration energies, Ecfg, probed along the trajectory provides information on near-equilibrium structures, MD trajectories have been run in the 10−50 K temperature range, analyzing the corresponding statistics. Simulations have been carried out at constant total energy, using a time step of 0.001 ps and running trajectories for 10 ns (after being equilibrated for 0.5 ns). The statistics outcome of a MD trajectory of the SH2−C6H6 dimer at a temperature value of 10 K is shown in Figure 5. The figure focuses on the intermolecular (S−C6H6) distance distribution (top panel), the angular distribution (middle panel) with α being the angle between the intermolecular (S−C6H6) axis and the C2v symmetry axis, and the average (potential) configuration energy (bottom panel). Clear maxima can be identified in all three distributions shown in the figure at about of 3.75 Å (top panel), 30° (middle panel) and −2.8 kcal mol−1 (lower panel). The C6H6−NH3 system has been studied analogously, running MD trajectories and extracting statistical averages. The results of this analysis are shown in Figure 6. The distribution of the distances from N to the benzene c.m., RN−C6H6, evaluated along the trajectory (top panel of Figure 6) shows a narrow peak at about 3.55 Å. The middle panel of Figure 6 shows two well-defined peaks in the distribution of the distances from the three H atoms of NH3 to the center of mass of C6H6, RH−C6H6. One of them is placed at shorter distances than RN−C6H6, while the other is located at larger distances. These observations indicate that one of the H atoms of NH3 tends to point toward the aromatic plane, while the other two tend to be positioned equidistantly at larger distances from C6H6. The interaction energy associated with these configurations, as can be derived from the energy distribution shown in the lower panel of Figure 6, is of about −2.6 kcal mol−1. A better estimation of the total potential energy at equilibrium, V(total)eq, can be obtained from a linear extrapolation to T = 0 K of the mean values of Ecfg ( Ecfg ), evaluated along trajectories run at several (low) temperatures (see, for

Figure 5. Relevant distributions for the C6H6−SH2 system derived from MD simulations at T = 10 K. Top panel: distances from the S atom of SH2 to the center of mass of C6H6; middle panel: angles formed by the dimer intermolecular axis and the C2 axis of SH2 and lower panel: configuration energy values (for the definition of distances and angles see the top panel of Figure 8).

Figure 6. Relevant distributions for the C6H6−NH3 system derived from MD simulations at T = 10 K. Top panel: distances from the N atom of NH3 to the center of mass of C6H6. Middle panel: distances from the H atoms of NH3 to the center of mass of C6H6. Lower panel: configuration energy values (for the definition of distances and angles see the lower panel of Figure 8). 1656

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

theoretical works addressed to characterize computationally the SH−π interaction, can be found in the literature,24,52,53,61−67 we have only found one experimental (rotational spectroscopy) report providing information on the geometry of the C6H6− SH2 system. Therefore, the model predictions of the binding energy can be only compared with theoretical results. For the C6H6−NH3 dimer, the predicted values of the RN−C6H6 distance and of the angle α (formed by the vector going from the C6H6 c.m. to N and the C3 axis of NH3) amount to 3.523 Å and 50.2°, respectively, while the value of the angle (not shown in Figure 8) formed by the C6 axis of C6H6 and the

instance, ref 31). Under these conditions, without phase changes of the system, Ecfg varies linearly with T. In Figure 7, Ecfg is represented as a function of temperature for the C6H6− SH2 and C6H6−NH3 systems in the top and lower panels, respectively.

Figure 7. Mean configuration energy, Ecfg , of C6H6−SH2 (top panel) and C6H6−NH3 (lower panel) as a function of temperature for very low T values. The labeled arrow indicates the potential energy value extrapolated at T = 0 K.

The process of running consecutive MD simulations at decreasing lower temperatures can be regarded as a simulated annealing minimization.57 However, given that this procedure can run into problems whenever the thermal energy is of order of the isomerization energy barriers, we judged better to extrapolate to 0 K values rather than running MD trajectories below 5 K. The extrapolated value of Ecfg at T = 0 K, represented by V(total)eq, is equal to −2.881 and −2.618 kcal mol−1 for the C6H6−SH2 and the C6H6−NH3 systems, respectively. These values are very similar to the lowest potential energy values sampled by the MD trajectory and equal to −2.843 and −2.606 kcal mol−1, respectively. Accordingly, it can be expected that the configurations in the MD trajectory with an interaction potential energy equal to −2.843 kcal mol−1 for C6H6−SH2 and equal to −2.606 kcal mol−1 for C6H6−NH3 should be similar to the equilibrium ones. Such configurations have been extracted from the simulations and are represented in Figure 8. For the C6H6−SH2 dimer, the predicted values of the RS−C6H6 distance and of the angle α equal to 3.692 Å and 29.2°. Similar values of RS−C6H6 and α are obtained at MP2/aug-cc-pVDZ level (3.625 Å and 27°), with an interaction energy of −3.05 kcal mol−1.58 Tauer et al.,52 applying the coupled-cluster theory through perturbative triple substitutions, CCSD(T),59 with the aug-cc-pVDZ basis set, estimate values of RS−C6H6= 3.8 Å and α = 30°, with an interaction energy of −2.74 kcal mol−1 (complete basis set extrapolation yield a CCSD(T) interaction energy of −2.81 kcal mol−1). The geometrical parameters coming out of the MD simulations are also in agreement with the data reported from the measured rotational spectrum (RS−C6H6 = 3.818 Å and α = 28.5°).60 Unfortunately, while some

Figure 8. Equilibrium geometry of the C6H6−SH2 (top panel) and of C6H6−NH3 (lower panel) systems.

C3 one of NH3 is equal to 58.8°, which is in good agreement with that of about 60° obtained from a high-resolution optical and microwave spectra of benzene−ammonia.54 Moreover, the difference between the predicted value of the RN−C6H6 distance (3.523 Å) and that derived from the experiment (equal to 3.590 Å) is less than a 2%. The fact that the tilted structure with an H atom pointing toward the aromatic plane is the most stable one can be explained pointing out that the bonding between ammonia and benzene is mostly electrostatic in origin. As has been explained before,68 C6H6 and NH3 exhibit a charge distribution compatible with quadrupole and dipole moments, respectively. The asymptotic dipole−quadrupole interaction favors approaches of NH3 with the C3 axis pointing straight toward the C6H6, while the quadrupole−quadrupole interaction favors approaches of NH3 with the C3 axis parallel to the aromatic plane. The observed tilt is a compromise between these two electrostatic terms. 1657

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

4. STRENGTH AND ANISOTROPY OF THE XH−C6H6 INTERACTION The present results for SH2 and NH3 can be used to complement previous studies on CH4 and H2O and analyze the strength of the XH-C6H6 interaction. The results for the three first clusters X (C, N, O) provide information on the evolution of the interaction nature for some first-row hydrides, while by comparing H2O−C6H6 and SH2−C6H6, we analyze the interaction between benzene and two hydrides of the same column of the periodic table. One one hand, the first-row hydrides are characterized by an increase of the dipole moment, μ going from CH4, without permanent dipole moment, to NH3 and H2O with experimental values of μ equal to 1.47 and 1.85 D, respectively.69 On the contrary, the molecular polarizability, α decreases in going from CH4 to H2O (αCH4 > αNH3 > αH2O). On the other hand, the polarizability of SH2 is higher than that of CH4 and much higher than that of H2O, while the SH2 dipole moment is approximately a half of that of H2O. Accordingly, differences can be expected in both the XH-C6H6 electrostatic and non electrostatic contributions, depending on the nature of the central atom X. In previous studies, using the same potential model described in section 2, some of us predicted interaction energies for the equilibrium configurations of H2O−C6H637 and CH4−C6H638 equal to −2.975 and −1.420 kcal mol−1, respectively. Accordingly, the relative order of the binding energies is H2O−C6H6 > SH2−C6H6 > NH3−C6H6 > CH4−C6H6. These results are in agreement with the ab initio data of Cheney et al. on some hydrogen-bond complexes involving benzene as an H acceptor.70 The mentioned study concludes that, although the magnitude of the binding energy varies using different basis sets and including electronic correlation, the qualitative trends for the series of systems remains essentially unchanged. This trend on the binding energies emphasized above is also consistent with the electronegativity order of the proton-donating atoms (C < N < O). This suggests that the increase of the strength of the XH−π interaction from X = C to X = O is mainly due to electrostatic effects (see also ref 21). However, when the electrostatic component in the XH-C6H6 system is weak or null, the interaction is dominated by dispersion attraction, although this effect may be hidden when electrostatic components are dominant. For instance, the CH4 molecule, without permanent dipole and quadrupole moments (having only an octupole moment), at large and intermediate distances basically interacts with C6H6 through dispersion forces.38 As a matter of fact, experimental77 and theoretical1,71−76 evidence on the dominance of the dispersion on the CH−π interaction can be found in the literature. The NH3−C6H6 and the H2O−C6H6 systems, due to the permanent dipole moments of NH3 and H2O, are bound by both dispersion (weaker than for CH4− C6H6) and electrostatic forces. This originates an extra stability in the NH3−C6H6 and the H2O−C6H6 clusters with respect to that of CH4−C6H6. The CH4−C6H6, NH3−C6H6 and H2O− C6H6 interaction energies are compared in Figure 9: for a intermolecular distance scan in the axial configuration. The obtained results are very similar to those of Tsuzuki et al.,21 which argue that the different size of the van der Waals radii of the proton-donating atoms (C, 1.75 Å; N, 1.55 Å; O, 1.40 Å), which scales approximately as the cubic square of the molecular polarizability, would be one of the main causes of the different intermolecular distances at the minimum of the potential energy. Moreover, confirming the prediction of our

Figure 9. CH4−C6H6, NH3−C6H6, and H2O−C6H6 intermolecular interaction potential energy curves varying the intermolecular distance maintaining an axial configuration.

model, the authors of ref 21 indicate that the electrostatic contribution is mainly responsible for the magnitude of the total interaction. On the other hand, SH2, with a permanent dipole moment of 0.97 D and a molecular polarizability higher than that of NH3 and H2O, forms SH2−C6H6 complexes stabilized again by the combination of both dispersion and electrostatic forces. However, while the electrostatic energy contribution accounts for about a 55% of the total interaction in the SH2−C6H6 dimer, it contributes a 76% in the H2O−C6H6 case. However, due to the larger polarizability, dispersion forces provide a higher contribution in SH2−C6H6 than in H2O−C6H6. The H2O−C6H6 and SH2−C6H6 interaction energies, for the cut of the PES corresponding to the most favorable approach, are compared in Figure 10, where noticeable differences can be seen in the repulsive zone of the PES.

Figure 10. H2O−C6H6 and SH2−C6H6 intermolecular interaction potential energy curves varying the intermolecular distance maintaining an axial configuration.

The equilibrium distance values, RX−C6H6, predicted by our model equal to 3.391 and 3.692 Å for H2O−C6H6 and SH2− C6H6, respectively, are in agreement with the distances reported by ab initio calculations, equal to 3.35 and 3.80 Å for H2O− C6H6 and SH2−C6H6, respectively.53 For the CH4−C6H6, NH3−C6H6, H2O−C6H6, and SH2− C6H6 series, the relevant interaction energy values are given in Table 4, where the predicted data are compared with the converged results of the CCSD(T) interaction energies toward Table 4. Interaction Energies at Equilibrium

1658

CH4−C6H6

NH3−C6H6

SH2−C6H6

H2O−C6H6

methodology

−1.420 −1.39 −1.41

−2.618 −2.47 −2.50

−2.881 −2.79 −2.85

−2.975 −3.09 −3.16

present aug-cc-pVQZ53 CBS53

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

the complete basis set limit obtained by Crittenden et al.53 The agreement, more than acceptable, could be even improved by a fine-tuning of the value of the non transferable β parameter (here maintained fixed), but this is outside the purpose of the present work. The strength of the XH−π interaction in the first row hydrides benzene compounds follows the electronegativity order of the proton-donating atoms (C < N < O). This indicates that the increase of the strength of the XH−π interaction going from X = C to X = O is mainly due to electrostatic forces.21 However, aspects other than the electronegativity, such as the molecular structure and polarizability, need to be involved. For instance, the electronegativity of S and C is very similar but the C6H6−SH2 interaction energy is different than the C6H6−CH4 one.24 Actually, the strength of the C6H6−SH2 and the C6H6−H2O interactions is very similar, with total interaction energies equal to 2.85 and 3.16 kcal mol−1 respectively,53 while lower values are reported for C6H6−CH4. From the previous results, it seems that when X is N, O or S, the interaction is dominated by an electrostatic contribution, while when X is C, the interaction has a more dispersive character. The small relevance of the electrostatic contribution in the interaction between CH4 (highly symmetric) and C6H6 is maintained even when polar groups are introduced in the benzene core, giving support to the small relevance of electrostatic effects in alkyl−benzene interactions.78 One of the main advantages of the proposed formulation of the interaction energy is the possibility to explore the PES for energetically less favorable approaches of the partners, analyzing anisotropic effects. In particular, the axial approaches in which SH2, NH3, H2O, and CH4 are rotated (in a plane perpendicular to the aromatic ring) about 180° with respect to the geometry of the most stable structures have been investigated. Results indicate that the molecular rotation generates a general increase of the intermolecular distance and an inverse order of the adducts stability (see Figure 11).

Figure 12. Rotational energy barriers for SH2−C6H6 (top panel) and H2O−C6H6 (lower panel) calculated at 5 Å, 7 Å and at the equilibrium intermolecular distance.

and energy, specially when the barrier high becomes comparable with the rotational energy of the partners.

5. CONCLUDING REMARKS A potential model based on the decomposition of molecular polarizabilities has been exploited to formulate the PES for the NH3−C6H6 and SH2−C6H6 systems and to investigate the role of the XH−π interaction along the CH4−C6H6, NH3−C6H6, H2O−C6H6 and SH2−C6H6 series of adducts. The nonelectrostatic component of the total interaction for each system has been constructed by summing 12 molecule-bond (CC and CH) components expressed by means of the ILJ function. The relevant parameters of the ILJ functions, having a physical meaning, have been derived from the average molecular polarizability of NH3 and SH2 and from the polarizability associated with the CC and CH bonds of C6H6. On another hand, the electrostatic component has been represented by applying the Coulomb law to 18 point charges placed in the C6H6 frame and 4 and 3 point charges compatible with the dipole moments of the NH3 and SH2 molecules, respectively. It has been found that the present formulation gives binding energies associated with specific configurations of the systems in agreement with high quality ab initio results. It has also been found that the lowest binding energy is that of CH4−C6H6, arising almost exclusively by dispersion forces. By increasing the number of electrons on the hydrides, the number of XH bonds decrease but the interaction becomes stronger. Accordingly, the NH−π interaction is weaker and stronger than the OH−π and the CH−π ones, respectively. On the other hand, the comparison of the OH−π and the SH−π interactions reveals the influence of both, the nonelectrostatic (depending on the molecular polarizability) and the electrostatic (depending on the dipole and quadrupole moments) interaction components. Reorientation effects induced by the anisotropy of intermolecular noncovalent forces in these prototype systems have been also investigated. Therefore, the obtained PES’s can be used to define the stereodynamics of molecular collisions under a variety of conditions by performing molecular dynamics simulations out molecular dynamic simulations.

Figure 11. H2O−C6H6, SH2−C6H6, NH3−C6H6, and CH4−C6H6 intermolecular interaction potentials for the represented molecular approaches.

The previous results suggest lower orientational effects for C6H6−CH4 in comparison with those present in the remaining systems. In particular, these effects have been investigated in detail at long-range for the two systems having similar symmetry: H2O−C6H6 and SH2−C6H6. The corresponding results are presented in Figure 12, where the value of ϕ = 0 corresponds to a XH bond placed along the y > 0 axis. Figure 12 shows some barriers on the PES’s and their dependence on the intermolecular distance. These anisotropy effects in the long-range region of the interaction can induce orientational effects in the colliding partners, depending on the temperature 1659

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

Potential of Carbon Dioxide Dimer from Symmetry-Adapted Perturbation theory. J. Chem. Phys. 1999, 110, 3785−3803. (17) Misquitta, A. J.; Szalewicz, K. Intermolecular Forces from Asymptotically Corrected Density Functional Description of Monomers. Chem. Phys. Lett. 2002, 357, 301−306. (18) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Dispersion Energy from Density-Functional Theory Description of Monomers. Phys. Rev. Lett. 2003, 91, 033201(1)−033201(4). (19) Shibasaki, K.; Fujii, A.; Mikami, N.; Tsuzuki, S. Magnitude and Nature of Interactions in Benzene-X (X = Ethylene and Acetylene) in the Gas Phase: Significantly Different CH/π Interaction of Acetylene as Compared with those of Ethylene and Methane. J. Phys. Chem. A 2007, 111, 753−758. (20) Fujii, A.; Shibasaki, K.; Kazama, T.; Itaya, R.; Mikami, N.; Tsuzuki, S. Experimental and Theoretical Determination of the Accurate Interaction Energies in Benzene-Halomethane: the Unique Nature of the Activated CH/π Interaction of Haloalkanes. Phys. Chem. Chem. Phys. 2008, 10, 2836−2843. (21) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. Origin of the Attraction and Directionality of the NH/π Interaction: Comparison with OH/π and VH/π Interactions. J. Am. Chem. Soc. 2000, 122, 11450−11458. (22) Ma, J.; Alfè, D.; Michaelides, A.; Wang, E. The Water-Benzene Interaction: Insight from Electronic Structure theories. J. Chem. Phys. 2009, 130, 154303(1)−154303(6). (23) Deligkaris, C.; Rodriguez, J. H. Correction to DFT Interaction Energies by an Empirical Dispersion Term Valid for a Range of Intermolecular Distances. Phys. Chem. Chem. Phys. 2012, 14, 3414− 3424. (24) Sherrill, C. D.; Takatani, T.; Hohenstein, E. G. An Assessment of Theoretical Methods for Nonbonded Interactions: Comparison to Complete Basis Set Limit Coupled-Cluster Potential Energy Curves for the Benzene Dimer, the Methane Dimer, Benzene-Methane, and Benzene-H2S. J. Phys. Chem. A 2009, 113, 10146−10159. (25) Podeszwa, R.; Bukowski, R. Szalewicz. Potential Energy Surface for the Benzene Dimer and Perturbational Analysis of π-π Interactions. J. Phys. Chem. A 2006, 110, 10345−10354. (26) Pirani, P.; Brizi, S.; Roncaratti, L. F.; Casavecchia, P.; Cappelletti, D.; Vecchiocattivi, F. Beyond the Lennard-Jones Model: a Simple and Accurate Potential Function Probed by High Resolution Scattering Data Useful for Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2008, 10, 5489−5503. (27) Albertí, M.; Pirani, F. Features of Ar Solvation Shells in Neutral and Ionic Clustering: The Competitive Role of Two-Body and ManyBody Interactions. J. Phys. Chem. A 2011, 115, 6394−6404. (28) Faginas Lago, N.; Huarte-Larrañaga, F.; Albertí, M. On the Suitability of the ILJ Function to Match Different Formulations of the Electrostatic Potential for Water-Water Interactions. Eur. Phys. J. D 2009, 55, 75−85. (29) Albertí, M.; Costantini, A.; Laganà, A.; Pirani, F. Are Micelles Needed to Form Methane Hydrates in Sodium Dodecyl Sulfate Solutions? J. Phys. Chem. B 2012, 116, 4220−4227. (30) Albertí, M. Rare Gas-Benzene-Rare Gas Interactions: Structural Properties and Dynamic Behavior. J. Phys. Chem. A 2010, 114, 2266− 2274. (31) Albertí, M.; Faginas Lago, N.; Laganà, A.; Pirani, F. A Portable Intermolecular Potential for Molecular Dynamics Studies of NMANMA and NMA-H2O aggregates. Phys. Chem. Chem. Phys. 2011, 13, 8422−8432. (32) Albertí, M.; Pacifici, L.; Laganà, A.; Aguilar, A. A Molecular Dynamics Study for the Isomerization of Ar Solvated (benzene)2-K+ heteroclusters. Chem. Phys. 2006, 327, 105−111. (33) Albertí, M.; Aguilar, A.; Pirani, F. Cation-π-Anion Interaction in Alkali Ion-Benzene-Halogen Ion Clusters. J. Phys. Chem. A 2009, 113, 14741−14748. (34) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. A Generalized Formulation of Ion-π Electron Interactions: Role of the Nonelectrostatic Component and Probe of the Potential Parameter Transferability. J. Phys. Chem. A 2010, 114, 11964−11970.

The extension of this formulation to aggregates formed by C6H6 and two or more molecules containing XH bonds is not only possible but is also useful to investigate both the stability and other basic features of the solvation shells of clusters formed by hydrogenated molecules around benzene.



AUTHOR INFORMATION

Corresponding Author

*(M.A.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS M.A., A.A. and F.H.-L. acknowledge financial support from the Ministerio de Educación y Ciencia (Spain, Project CTQ201016709) and the Generalitat de Catalunya (2009SGR-17). Also thanks are due to the Center de Supercomputació de Catalunya CESCA-C4 and Fundació Catalana per a la Recerca for the allocated supercomputing time. F.P. acknowledges financial support from the Italian Ministry of University and Research (MIUR) for PRIN Contracts.



REFERENCES

(1) Tsuzuki, S.; Fujii, A. Nature and Physical Origin of CH/π Interaction: Significant Difference from Conventional Hydrogen Bonds. Phys. Chem. Chem. Phys. 2008, 10, 2584−2594. (2) Buckingham, A. D.; Fowler, P.; Hutson, J. M. Theoretical Studies of van der Waals Molecules and Intermolecular Forces. Chem. Rev. 1988, 88, 963−988. (3) Chalasinski, G.; Szczesniak, M. M. State of the ab Initio Theory of Intermolecular Interactions. Chem. Rev. 2000, 100, 4227−4252. (4) Meyer, E. A.; Castellano, R. K.; Diederich, F. An Efficient Algorithm for the Density-Functional Theory Treatment of Dispersion Interactions. Angew. Chem., Int. Ed. 2003, 42, 1210−1250. (5) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond in Structural Chemistry and Biology; Oxford University Press: New York, 2001. (6) Grabowsky, S. G. Hydrogen Bonding: New Insights; Springer Verlag: New York, 2006. (7) Geffrey, G. A. Introduction to Hydrogen Bonding; Oxford University Press: New York, 1997. (8) Saggu, M.; Levinson, N. M.; Boxer, S. G. Experimental Quantification of Electrostatics in X−H···π Hydrogen Bonds. J. Am. Chem. Soc. 2012, 134, 18986−18997. (9) Burley, S. K.; Petsko, G. A. Aromatic-Aromatic Interaction: A Mechanism of Protein Structure stabilization. Science 1985, 229, 23− 28. (10) Steiner, T.; Koellner, G. Hydrogen Bonds with π-Acceptors in Proteins: Frequencies and Role in Stabilizing Local 3D Structures. J. Mol. Biol. 2001, 305, 535−557. (11) Brandl, M.; Weiss, M. S.; Jabs, A.; Suhnel, J.; Hilgenfeld, R. CH···π Interactions in Proteins. J. Mol. Biol. 2001, 307, 357−377. (12) Ringer, A. L.; Senenko, A.; Sherrill, D. Models of S/π Interactions in Protein Structures: Comparison of the H2S-Benzene Complex with PDB Data. Protein Sci. 2009, 16, 2216−2223. (13) Shibasaki, K.; Fujii, A.; Mikami, N.; Tsuzuki, S. Magnitude of the CH/π Interaction in the Gas Phase: Experimental and Theoretical Determination of the Accurate Interaction Energy in Benzenemethane. J. Phys. Chem. A 2006, 110, 4397−4404. (14) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887−1930. (15) Mas, E. M.; Szalewicz, K.; Bukowski, R.; Jeziorski, B. Pair Potential for Water from Symmetry-Adapted Perturbation Theory. J. Chem. Phys. 1997, 107, 4207−4218. (16) Bukowski, R.; Sadlej, J.; Jeziorski, B.; Jankowsky, P.; Szalewicz, K.; Kucharski, S. A.; Williams, H. L.; Rice, B. M. Intermolecular 1660

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

(35) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Coletti, C.; Re, N. Atom-Bond Pairwise Additive Representation for Halide-Benzene Potential Energy Surfaces: an Ab Initio Validation Study. J. Phys. Chem. A 2009, 113, 14606−14614. (36) Stone, A. J.; Misquitta, A. J. Atom-Atom Potentials from Ab Initio Calculations. International Reviews in Physical Chemistry 2007, 26, 193−222. (37) Albertí, M.; Faginas Lago, N.; Pirani, F. Benzene Water Interaction: From Gaseous Dimers to Solvated Aggregates. Chem. Phys. 2012, 399, 232−239. (38) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. Competitive Role of CH4-CH4 and CH-π Interactions in C6H6-(CH4)n Aggregates: The Transition from Dimer to Cluster Features. J. Chem. Phys. A 2012, 116, 5480−5490. (39) http://cccbdb.nist.gov/MP2FC/6-311G* (40) Albertí, M.; Castro, A.; Laganà, A.; Moix, M.; Pirani, F.; Cappelletti, D.; Liuti, G. A Molecular Dynamics Investigation of RareGas Solvated Cation-Benzene Clusters Using a New Model Potential. J. Phys. Chem. A 2005, 109, 2906−2911. (41) Denbigh, K. G. The Polarizabilities of Bonds-I. Trans. Faraday. Soc. 1940, 36, 936−948. (42) Zitto, M. E.; Caputo, M. C.; Ferraro, M. B. Resolution of Molecular Polarizabilities of CH3-X and CH3- CH2-X Derivatives into Atomic Terms. |it J. Chem. Phys. 2001, 114, 4053−4057. (43) Lillestolen, T. C.; Wheatley, R. J. First-Principles Calculation of Local Atomic Polarizabilities. J. Phys. Chem. A 2007, 111, 11141− 11146. (44) Hirschfelder, J. O.; Curtis, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids: Whiley: New York, 1967. (45) Smith, R. P.; Mortensen, E. J. Bond and Molecular Polarizability Tensors. I. Mathematical Treatment of Bond Tensor Additivity. J. Chem. Phys. 1960, 32, 502−507. (46) Olney, T. N.; Cann, N. M.; Cooper, G.; Brion, C. E. Absolute Scale Determination for Photoabsorption Spectra and the Calculation of Molecular Properties using Dipole Sum-Rules. Chem. Phys. 1997, 59−98. (47) Pirani, F.; Albertí, M.; Castro, A.; Moix, M.; Cappelletti, D. Atom-Bond Pairwise Additive Representation for Intermolecular Potential Energy Surfaces. Chem. Phys. Lett. 2004, 394, 37−44. (48) Pirani, F.; Cappelletti, D.; Liuti, G. Range Strength and Anisotropy of Intermolecular Forces in Atom-Molecule Systems: An Atom Bond Pairwise Additive Approach. Chem. Phys. Lett. 2001, 350, 286−296. (49) Cappelletti, D.; Ronca, E.; Belpassi, L.; Tarantelli, F.; Pirani, F. Revealing Charge Transfer in Gas-Phase Water Chemistry. Acc. Chem. Res. 2012, 45, 1571−1580. (50) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Cappelletti, D.; Coletti, C.; Re, N. Atom-Bond Pairwise Additive Representation for Cation-Benzene Potential Energy Surfaces: an Ab Initio Validation Study. J. Phys. Chem. A 2006, 110, 9002−9010. (51) Pirani, F.; Cappelletti, D.; Liutu, G. Range, Strength and Anisotropy of Intermolecular Forces in Atom-Molecule Systems: An Atom-Bond Pairwise Additive Approach. Chem. Phys. Lett. 2001, 350, 286−296. (52) Tauer, T. P.; Derrick, M. E.; Sherrill, C. D. Estimates of the ab Initio Limit for Sulfur-π Interactions: The H2S-Benzene Dimer. J. Phys. Chem. A 2005, 109, 191−196. (53) Crittenden, D. L. A Systematic CCSD(T) Study of Long-Range and Noncovalent Interaction between Benzene and a Series of Firstand Second-Row Hydrides and Rare Gas Atoms. J. Phys. Chem. A 2009, 113, 1663−1669. (54) Rodham, D. A.; Suzuki, S.; Suenram, R. D.; Lovas, F. J.; Dasgupta, S.; Godard, W. A., III; Blake, G. A. Hydrogen Bonding in the Benzene-Ammonia Dimer. Nature 1993, 362, 735−737. (55) Bloom, J. W. G.; Raju, R. K.; Wheeler, S. E. Physical Nature of Substituent Effects in XH/π Interactions. J. Chem. Theory Comput. 2012, 8, 3167−3174. (56) http://www.ccp5.ac.uk/DL_POLY_CLASSIC/.

(57) Kirkpatrick, S.; Gelatt, C. D., Jr.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671−680. (58) Wang, Y.; Paulus, B. A Comparative Electron Correlation Treatment in H2S-Benzene Dimer with DFT and Wavefunction-Based ab Initio Methods. Chem. Phys. Lett. 2007, 441, 187−193. (59) Raghavachari, K.; Trucks, K.; Pople, J. A.; Head-Gordon, M. A Fifth-Order Perturbation Comparison of Electron Theories. Chem. Phys. Lett. 1989, 157, 479−483. (60) Arunan, A.; Emilsson, T.; Gutowsky, H. S.; de Oliveira, G.; Dykstra, C. E. Rotational spectrum of the weakly bonded C6H6-H2S dimer and comparisons to C6H6-H2O dimer. J. Chem. Phys. 2002, 117, 9766−9776. (61) Cabaleiro-Lago, E. M.; Carrazana-Garcia, J. A.; RodríguezOtero, J. Study of the Interaction between Water and Hydrogen Sulfide with Polycyclic Aromatic Hydrocarbons. J. Chem. Phys. 2009, 130, 234307(1)−234307(11). (62) Hermida-Ramon, J. M.; Cabaleiro-Lago, E. M.; RodríguezOtero, J. Theoretical Characterization of Structures and Energies of Benzene-(H2S)n and (H2S)n (n = 1−4) clusters. J. Chem. Phys. 2005, 122, 204315(1)−204315(5). (63) Cabaleiro-Lago, E. M.; Rodríguez-Otero, J.; Peña-Gallego, A. Computational Study on the Characteristics of the Interaction in Naphthalene...(H2X)n=1,2 (X = O, S) Clusters. J. Phys. Chem. A 2008, 112, 6344−6350. (64) Cabaleiro-Lago, E. M.; Rodríguez-Otero, J.; Peña-Gallego, A. Characteristics of the Interaction of Azulene with Water and Hydrogen Sulfide: A Computational Study. J. Chem. Phys. 2008, 129, 084305(1)−084305(8). (65) Duan, G. L.; Smith, V. H., Jr.; Weaver, D. F. Characterization of Aromatic-Thiol π-Type Hydrogen Bonding and PhenylalanineCysteine Side Chain Interactions through ab Initio Calculations and Protein Database Analyses. Mol. Phys. 2001, 99, 1689−1699. (66) Tarakeshavar, P.; Soon Choi, H.; Kim, K. S. Olefinic vs Aromatic π-H Interaction: A Theoretical Investigation of the Nature of Interaction of First-row Hydrides with Ethene and Benzene. J. Am. Chem. Soc. 2001, 123, 3323−3331. (67) Mons, M.; Dimicoli, I.; Tardivel, B.; Piuzzi, F.; Brenner, V.; Millié, P. Energetics of a Model NH-π Interaction: the Gas Phase Benzene-NH3 Complex. Phys. Chem. Chem. Phys. 2002, |it 4, 571−576. (68) Bauschlicher, C. W., Jr.; Ricca, A. Binding of NH3 To Graphite and To a (9,0) Carbon Nanotube. Phys. Rev. 2004, 70, 115409(1)− 115409(6). (69) Weast, R. Handbook of Chemistry and Physics; The Chemical Rubber Co.: Cleveland, OH, 1967. (70) Cheney, B. V.; Schulz, M. W.; Cheney, J.; Richards, G. Hydrogen-Bonded Complexes Involving Benzene as an H Acceptor. J. Am. Chem. Soc. 1988, 110, 4195−4198. (71) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. High-Level ab Initio Calculations of Interaction Energies of C2H4-CH4 and C2H6-CH4 Dimers: A Model Study of CH/π Interaction. J. Phys. Chem. A 1999, 103, 8265−8271. (72) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. The Magnitude of the CH/π Interaction between Benzene and Some Model Hydrocarbons. J. Am. Chem. Soc. 2000, 122, 3746−3753. (73) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. The Interaction of Benzene with Chloro- and Fluoromethanes: Effects of Halogenation on CH/π Interaction. J. Phys. Chem. A 2002, 106, 4423−4428. (74) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Fujii, A. Magnitude and Directionality of the Interaction Energy of the Aliphatic CH/π Interaction: Significant Difference from Hydrogen Bond. J. Phys. Chem. A 2006, 110, 10163−10168. (75) Morita, S.; Fujii, A.; Mikami, N.; Tsuzuki, S. Origin of the Attraction in Aliphatic C-H/π Interactions:? Infrared Spectroscopic and Theoretical Characterization of Gas-Phase Clusters of Aromatics with Methane. J. Phys. Chem. A 2006, 110, 10583−10590. (76) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Fujii, A. CH/π Interactions in Methane Clusters with Polycyclic Aromatic Hydrocarbons. Phys. Chem. Chem. Phys. 2008, 10, 2860−2865. 1661

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662

The Journal of Physical Chemistry A

Article

(77) Fujii, A.; Hayashi, H.; Park, J. W.; Kazama, T.; Mikami, N.; Tsuzuki, S. Experimental and Theoretical Determination of the Accurate CH/π Interaction Energies in Benzene-Alkane Clusters Correlation Between Interaction Energy and Polarizability. Phys. Chem. Chem. Phys. 2011, 13, 14131−14141. (78) Ribas, J.; Cubero, J.; Luque, F. J.; Orozco, M. Theoretical Study of Alkyl-π and Aryl-π Interactions. Reconciling Theory and Experiment. J. Org. Chem. 2002, 67, 7057−7065.

1662

dx.doi.org/10.1021/jp410917x | J. Phys. Chem. A 2014, 118, 1651−1662