Competitive Role of CH4–CH4 and CH− π Interactions in C6H6–(CH4

Article pubs.acs.org/JPCA

Competitive Role of CH4−CH4 and CH−π Interactions in C6H6−(CH4)n Aggregates: The Transition from Dimer to Cluster Features M. Albertí,*,† A. Aguilar,† J. M. Lucas,† and F. Pirani‡ †

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

‡

ABSTRACT: The intermolecular methane−methane and benzene (Bz)−methane interactions formulated in this paper are suitable to investigate systems of increasing complexity. The proposed CH4−CH4 and Bz−CH4 potential energy functions are indeed applied to study some macroscopic properties of methane and important features of both small Bz−(CH4)n (n > 1−10) clusters and Bz surrounded by several CH4 molecules. Relevant parameters of the interaction, 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 proposed potential model also allows one to isolate the role of the different intermolecular energy contributions. The spatial distribution of the CH4 molecules in the clusters is investigated by means of molecular dynamics simulations under various conditions, even when methane phase transition from liquid to gas is likely to occur. In addition, several properties, such as radial distribution functions, density values, and mean diffusion coefficients, are analyzed in detail.

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. In particular, at low temperature, such interactions can control the formation of the precursor state, by collision of reagents, which selectively evolves toward the final state of products. Therefore, the understanding of nonbonded or weak interactions is one of the major goals of chemistry.1 Particular attention has been paid to molecular aggregates involving aromatic rings for which noncovalent intermolecular interactions control basic phenomena such as the formation of weak hydrogen bonds,2,3 the competitive solvation of ions by different partners,4−6 molecular recognition and selection processes,7−9 and chemical processes occurring in biological systems.4,5,10 Some aggregates, like those containing cations and aromatic rings, are characterized by strong electrostatic components playing a prominent role in the whole interaction.7−9,11 However, the description of cation−π systems in terms of pure electrostatic forces only provides a qualitative picture of the situation, due to the key role played by other interaction components, like exchange or size repulsion, induction, and dispersion attraction,12,13 which scale in a © 2012 American Chemical Society

different way with the intermolecular distance. As a matter of fact, molecules containing π-electrons can interact with ions but also with neutral species such as, for instance, rare gases14−19 and other molecules with permanent multipole20−24 or without permanent dipole and quadrupole moments,25−36 and they can participate in the formation of cation−π, π−π, and XH−π (X = O, S, N, C, ...) intermolecular bonds, whose strength depends on the nature of X. For instance, the OH−π and the NH−π interactions are stronger than the CH−π ones.35 Actually, the strength of CH−π lies between the weakest class of hydrogen bonds and van der Waals bonds.31 However, while the electrostatic component is often dominant in hydrogen bonds, the dispersion forces can provide an important contribution to the CH−π interactions, with the electrostatic effects (which depend on the polarity of the donor group26,35) being, in general, small. Moreover, in some cases, the role of possible charge transfer (CT) effects, emerging at intermediate and short intermolecular distances, being dependent on the orbital overlap, must be properly taken into account, since they Received: March 12, 2012 Revised: May 4, 2012 Published: May 16, 2012 5480

dx.doi.org/10.1021/jp3023698 | J. Phys. Chem. A 2012, 116, 5480−5490

The Journal of Physical Chemistry A

Article

the decomposition of the molecular polarizability in bond polarizabilities,52,53 each one having parallel and perpendicular components related with the dimension and shape (often ellipsoidal) of the electronic charge distribution around the considered bond.52 This model has been recently applied to study several systems containing Bz, and the role of the XH−π interaction has been investigated in Bz−H2O.54 The Bz−H2O interaction energy was formulated by decomposing the polarizability of Bz (see refs 52, 53 for the decomposition details) in 12 bond polarizabilities assigned to each molecular bond of Bz (CC and CH). On the contrary, the polarizability of the small H2O molecule was not decomposed. Then, the nonelectrostatic Bz−H2O interaction was expressed as a sum of H2O-bond contributions (H2O−CC and H2O−CH), each one described by means of an improved Lennard-Jones (ILJ) function (see the next section and refs 55, 56). The same procedure was first applied to describe the weak Bz−rare gas14−16 and Bz−ion interactions (see, for instance, refs 57−63). In our model, the assumption of the additivity of the bond polarizability allows also one to associate effective values of polarizability not only to molecular bonds but also to atoms and groups of atoms within the molecules,48,57,64−68 whose sum must be, in any case, equal to the proper value of the whole molecular polarizability. Due to the fact that many body effects, associated with the formation of chemical bonds, modify substantially the electronic charge distribution, the effective atomic polarizability is different of that of free atom. From the polarizability of the bonds and that associated with atoms,14−16 ions (see, for instance, ref 12), or groups of atoms56,64−68 the basic parameters of the interaction, which are completely transferable57 since having a physical meaning, can be easily determined by exploiting widely used correlation formulas of general validity.69,70 The reliability of the present formulation was also extensively tested for several systems through comparisons of the model predictions, concerning both the total interaction energy and its decomposition in leading components, with ab initio results obtained at various levels of theory.71−76 In particular, for systems involving closed shell ions and benzene molecules, where a strong induction attraction is operative together with a pronounced electrostatic contribution, the theoretical energy decomposition has been obtained according to the Kitaura−Morokuma scheme77−79 and the symmetry-adapted perturbation theory (SAPT)80 for cation−benzene71 and anion−benzene,72 respectively. In the present paper, we are interested in investigating the energy and directionality of the weak Bz−CH4 interaction and in its competition with the CH4−CH4 one. The paper is organized as follows: in section 2, the potential energy function is described and the meaning of its parameters is analyzed. In section 3, the Bz−CH4 interaction is presented and analyzed as a function of the relative orientation of the partners considering two formulations of the potential energy, based on different decompositions of the molecular polarizability. The combined and competitive role of CH4−CH4 and Bz−CH4 interaction is emphasized in section 4 from the MD simulations carried out on (CH4)n and Bz−(CH4)n clusters. Finally, concluding remarks are reported in section 5.

provide important stabilization contributions to the intermolecular bond, and/or a reduction of the strength of the shortrange repulsion between involved partners. In spite of the fact that the extent of electrostatic forces depends on the nature (acidity) of the alkane CH bond interacting with the π cloud, the related energy contribution is always weaker than the dispersion one. When electrostatic components are dominant, the effect of the dispersion forces may be hidden. On the contrary, when the electrostatic interaction is weak or null, an accurate description of the dispersion and of other attractive components appears to be absolutely necessary. In the last years, some experimental37−39 and theoretical25,31−36,40 methodologies have been applied to investigate the CH−π interaction nature. The presence of alkyl (CH bonds) and phenyl (π-electron) groups in biomolecules suggests that the CH−π interaction can be crucial in defining the characteristics of several processes in which biomolecules are involved. Due to the presence of aliphatic groups in the side chains of alanine, valine, leucine, etc. and aromatic groups in phenylalanine, tyrosine, etc., the possible simultaneous participation of several CH bonds in the intermolecular interaction can originate its enhancement.41,42 However, the weakness of individual CH−π interaction contributions makes difficult the assessment of the relative role of the different components, and in this case, the use of small prototype systems, like Bz−CH4, can be very instructive. However, as the potential energy surfaces depend on 3N − 6 variables (see, for instance, refs 43−47), calculations of Bz−CH4 (N = 17) at the ab initio level have been only performed for selected configurations of the system. Tsuzuki et al.31,33,35,36 have widely investigated CH−π interactions proving that systematic CCSD(T) calculations, carried out using very large basis sets, provide accurate intermolecular interaction energies. The authors indicate that an accurate evaluation of the CH−π interaction energy requires a high computationally demanding CCSD(T) level of calculation using a large basis set near saturation. Although the Bz−methane and Bz−acetylene clusters can be investigated at the CCSD(T)/aug-cc-pVTZ level, such computationally demanding calculation is still not practical for larger systems.36 In general, the theoretical study of relatively large and/or complex systems requires the use of potential energy functions associated with small parts of the system that can be conveniently assembled to provide the description of the global interaction (force field). Often, the force field is implemented in molecular dynamics (MD) programs, and in this case, not only an accurate description of the interaction at equilibrium, derived from accurate ab initio calculations and/or high level experimental data, 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.48 Accordingly, to study large systems, for which the computational cost is too high, semiempirical and empirical methods can be very useful. However, a careful inspection of the potential energy functions is required before their use in MD simulations. In the last years we have modeled some nonbonded interactions on the basis of the assumption that polarizability is the key property to scale both dispersion/induction attraction and size repulsion in atomic and molecular systems, being suitable to define both the probability of induced dipole formation and the particle volume or size (see, for instance, refs 49−51). In the adopted potential energy model we also exploit

2. POTENTIAL ENERGY FUNCTION The total potential energy, Vtotal, is described, as usual, as a sum of electrostatic, Vel, and nonelectrostatic, Vnel, components, which must be considered as “effective interaction terms”, since 5481



Article

model, the balancing of a positive (repulsion) and a negative (attraction) term. Such formulation, exploiting a proper modulation with r (the distance between two interaction centers) of the relative strength of the two terms, removes most of the LJ inadequacies, specially in the asymptotic region.55 Relevant parameters of each pair interaction are the binding energy (well depth), ε; the equilibrium distance, r0; and β. The latter, controlling the strength of attraction and repulsion and the steepness of the repulsive wall, is a shape parameter. The VILJ function is then formulated as follows

the role of less important contributions (often providing opposite effects) is included. In our calculations, both Bz and CH4 have been kept rigid (CC and CH distances equal to 1.39 and 1.09 Å) and Vel is calculated easily from the charge distributions of both Bz and CH4 by means of Coulombic potentials. For Bz, the values of the charges and their positions are chosen to reproduce the correct components of the Bz quadrupole moment (see, for instance, ref 12), in which two negative charges of −0.04623 au have been placed on each C atom (above and below the aromatic plane separated by 1.905 Å) and one positive charge of 0.09246 au on each H atom. For CH4, charges of −0.62632 and 0.15658 au have been assigned to C and H atoms, respectively, which are consistent with its octupole moment. Vnel, in general, is more difficult to formulate, because it results from the balancing of both dispersion and induction attraction, dominant at long range, with exchange (size) repulsion, dominant at short range.69,70 In our potential model Vnel is obtained from the consideration that polarizability is the key property to scale simultaneously attraction and repulsion in weakly bound molecular systems, which can be justified by the good correlation of polarizability and relevant interaction properties.69 The interaction between small molecules, such as CH4 or H2O, can be described directly considering each molecule as a unique center with an assigned polarizability coincident with the average value of the molecular polarizability.54 However, when larger molecules are involved, it is convenient to decompose the molecular polarizability in several values assigned to different groups on the molecule, and in this case, several interaction centers with assigned values of the polarizability are considered. In the present study, CH4 has been then considered as an indecomposable group with the interaction center placed on the C atom (here represented by Cm). The choice of a single center for methane coincides with the recent proposal of the AMPF model potential of water56,64−66 that has been successfully used to predict the main properties of water (liquid and small clusters). On the contrary, the Bz polarizability has been decomposed into 12 contributions,52,53 assigned to the twelve bonds of Bz (each one having a parallel and a perpendicular component to the bond axis), whose weighted sum provides the proper value of the molecular polarizability. In this way, each bond is considered as an independent interaction group, whose interaction center has been considered in the middle of the CC and CH bonds (see, for instance, refs 12, 14). Then, the full interaction has been defined as a sum of CH4−CC and CH4−CH pair contributions. However, those systems containing a high number of CH4 molecules have been investigated by further decomposing bond polarizabilities in effective atomic polarizabilities,81,82 whose values are consistent with the literature data.83,84 These effective atomic polarizability values, as has been indicated above, differ from those of the isolated atoms and their sum provides, as for the previous decomposition, the proper average value of the molecular polarizability. In this case, the 12 interaction centers, here represented by Ceff and Heff, are placed on the C and H atoms and, then, the interaction between Bz and any CH4 is defined as a sum of CH4−Ceff and CH4−Heff contributions. Each pair interaction, independently of the decomposition of the molecular polarizability, is described by the ILJ function,55,56 VILJ, which involves, like the traditional LJ

VILJ

⎡ ⎢ m = ε⎢ ⎢ β + 4.0 r ⎢⎣ r0

r

2

⎛ r0 ⎞ β+ 4.0( r0 ) ⎜ ⎟ ⎝r⎠

2

( ) −m ⎤ β + 4.0( ) ⎛ r ⎞ ⎥ ⎥ − ⎝r⎠ ⎥ β + 4.0( ) − m ⎥⎦ r r0

r r0

2

2

m

⎜

0⎟

(1)

Accordingly with the bond decomposition of the Bz polarizability, the Bz−CH4 nonelectrostatic interaction is calculated as a sum of pair interactions of CH4-bond type (six Cm−CC and six Cm−CH), each one expressed by means of eq 1. The perpendicular and parallel components of the polarizability assigned to the CC and CH bonds,12 together with the molecular polarizability of CH4 equal to 2.56 Å3, assigned to the Cm center, allow the determination of the perpendicular and parallel components of the well depth (ε⊥ and ε∥) and of the intermolecular equilibrium distance (r0,⊥ and r0,∥) associated to each CH4-bond pair by using the abovementioned correlation formulas (see, for instance, refs 12, 15 and references therein). The ε and r0 values to be used in eq 1 are obtained as simple trigonometric combination of parallel and perpendicular components, whose coefficients depend on the angle that a given bond forms with r.⃗ 15 The values of ε⊥, ε∥, r0,⊥, and r0,∥ are given in Table 1, while the choice of the β value will be illustrated below. Table 1. Perpendicular and Parallel Components of the Well Depth (ε⊥, ε∥) and of the Equilibrium Distances (r0⊥, r0∥) for the Cm−CC and Cm−CH Interaction Pairsa CH4 bond

ε⊥/kcal mol−1

ε∥/kcal mol−1

r0⊥/Å

r0∥/Å

Cm−CC Cm−CH

0.1039 0.1206

0.1374 0.1019

4.006 3.792

4.292 3.997

a Cm is used to describe the interaction center placed on the C atom of CH4.

The average CC and CH bond polarizabilities can be further decomposed in atomic effective polarizabilities taking into account that a given atom on the molecule can participate in more than one bond. In this case, the Bz−CH4 nonelectrostatic interaction is calculated as a sum of 12 effective pair interactions (six Cm−Ceff and six Cm−Heff) also represented by the ILJ function but with ε and r0 derived directly from the polarizability assigned to the Ceff and Heff centers. The corresponding potential parameters are given in Table 2. Vnel, for the CH4−CH4 system, is also described by means of the ILJ function, and the related relevant parameters, εCH4−CH4 = 0.335 kcal mol−1 and r0,CH4−CH4= 4.06 Å, have been obtained 5482



Article

Table 2. Well Depth (ε) and the Equilibrium Distances (r0) for the Cm−Ceff and Cm−Heff Interaction Pairsa εCm−Ceff/kcal mol−1

εCm−Heff/ kcal mol−1

r0,Cm−Ceff/Å

r0,Cm−Heff/Å

0.1547

0.0811

4.093

3.743

a

Cm is used to describe the interaction center placed on the C atom of CH4, and Ceff and Heff are the effective interaction centers on the atoms of the Bz molecule.

using only the molecular polarizability of methane indicate above. In order to describe properly the dependence of the longrange attraction on r, for the present neutral−neutral systems the m parameter has been taken equal to 6. The VILJ (see eq 1) formulation suggests an enlargement of the negative region of the potential energy where the attraction dominates the repulsion (essentially the well located at intermediate values of r) by lowering β, while a reduction of the strength of the repulsion is expected in the region of the potential where the repulsion dominates the attraction (small r values). Moreover, for weak noncovalent intermolecular potentials, β assumes values practically independent of the relative orientation of the involved partners (being the interaction of the same nature) and for systems of the same type its value falls in a limited range of variation.55 For typical neutral−neutral van der Waals interactions (size dispersion plus dispersion attraction) β = 9 has been found to be an appropriate value.55 For such systems, the lowering of β with respect to 9 can indirectly include the effect of less important additional attractive components, such as CT in the perturbative limit, and its variation, when large molecules are involved, can account also for changes in the interaction reference centers.65 Therefore, its value can partially depend on the type of environment56 and on the framework of the centers adopted to represent the interaction (see below). 2.1. Influence of the β Parameter. Bearing in mind the inclusion (in an effective way) of possible orientationally averaged CT effects, a lowering of β to 7.0 (see eq 1), together with the parameters given in Table 1, has been used as a first test of the potential model predictive capability. Calculations for approaches from CH4 to Bz [with one of the CH bonds of CH4 placed along the C6 rotational axis of Bz and pointing toward its center of mass (c.m.)] have been performed. This particular approach has been considered because the experimental evidence of such equilibrium structures.28,31,85 In Figure 1 are represented Vtotal (top panel), Vnel, and Vel (medium panel) as a function of the distance from the C atom of CH4 to the Bz c.m., R. Moreover, as Vnel can be separated in two main contributions, the long-range “effective” attraction (Vattr) and the short-range “effective” repulsion (Vrep), which arise from the negative and positive terms in eq 1, respectively, Vattr and Vrep, are also reported in the same figure (lower panel). The different energy contributions at the equilibrium distance, Req = 3.653 Å, can be observed with the help of the vertical dashed line, going from the top to the lower panel of the figure. The red points indicate the energy values of the different contributions at the equilibrium configuration. The predicted value of Vtotal (−1.43 kcal mol−1) is very similar to that calculated by Tsuzuki et al.35 (−1.45 kcal mol−1), which corresponds to the CCSD(T) interaction energy at the basis set limit. In this study the authors adopt the SAPT80 scheme to obtain detailed information on the intermolecular interaction. The corresponding results allow one to compare several energy

Figure 1. Atom-bond formulation of the interaction with β = 7.0: Energy contributions, for approaches of CH4 along the C6 rotational axis of Bz (with an H atom of CH4 pointing to the Bz c.m.), as a function of R (the distance from the C atom of CH4 to Bz c.m.). Top panel: Total potential energy, Vtotal vs R. Middle panel: Nonelectrostatic, Vnel, and electrostatic, Vel, interaction energies vs R. Bottom panel: Vnel repulsive and attractive energy components, Vrep and Vattr, vs R. At the equilibrium distance (red points) Req = 3.653 Å, Vtotal = −1.43, Vnel = −1.19, Vel = −0.24, Vattr = −2.74, and Vrep = 1.54 (all energy contributions in kcal mol−1).

components. In particular, the energy predictions of −1.19 and −0.24 kcal mol−1 for Vnel and Vel (medium panel of Figure 1), respectively, are also in agreement with the values of −1.20 and −0.25 kcal mol−1 obtained by applying the SAPT methodology.35 However, the individual strength of Vattr and Vrep components indicates that predictions are overestimated by about 20% and that the predicted attraction is excessive. Such discrepancy could be partially or totally artificial, since it can arise from the difficulty to independently evaluate Vattr and Vrep, especially at intermediate values of r, because of their incomplete separability due to non-negligible damping effects of attraction attributed to overlap effects. Nevertheless, in order to assess the flexibility of the model, changes in Vattr, Vrep, and Vtotal have been evaluated for four values of β. In particular, Vattr and Vrep are represented as a function of R in Figure 2 and results at equilibrium are given in Table 3. Figure 2 shows that upon increasing β the function becomes globally more repulsive. As stressed above, such effect is more evident at short distances, i.e., when the repulsion dominates the attraction. Moreover, from Table 3 and accordingly with the expectations, it appears that at intermediate R values, the effect of increasing β is more evident for Vattr and Vrep than for Vtotal. The results given in Table 3 and their comparison with ab initio data35 suggest that β values between 8.5 and 9 should give good energy predictions.

3. BZ−CH4 INTERMOLECULAR INTERACTION The value of β = 8.75 is used to calculate the Bz−CH4 interaction, and the definitive results are compared again with 5483



Article

Figure 2. Atom-bond formulation of the interaction: Vrep (continuous line) and Vattr (dotted line) predictions, as a function of R, using different values of the β parameter.

Figure 3. Atom-bond formulation of the interaction with β = 8.75: Orientation dependence of the total interaction energy, Vtotal, and electrostatic energy, Vel, in the Bz−CH4 cluster. The ϕ angle equal to zero corresponds to an hydrogen atom of CH4 pointing toward the center of mass of Bz and with the C atom placed on the C6 rotational axis of Bz. The distance from the C atom of CH4 to the Bz c.m is the equilibrium value of 3.667 Å.

Table 3. Vattr, Vrep, and Vtotal Interaction Energies (in kcal mol−1) Calculated at the Equilibrium Distance from the C Atom of CH4 to the c.m. of Bz (Req, in Å) and for Different Values of the β Parameter Vattr

Vrep

Vtotal

Req

β

−2.74 −2.45 −2.39 −2.25

1.54 1.18 1.17 1.08

−1.43 −1.42 −1.42 −1.41

3.653 3.661 3.665 3.669

7.0 8.0 8.5 9.0

ab initio calculations. Results, for the same approach from CH4 to Bz as in the previous section, indicate that not only Vtotal (−1.42 kcal mol−1) and its electrostatic (−0.24 kcal mol−1) and nonelectrostatic (−1.18 kcal mol−1) contributions are in agreement with ab initio data, but also Vattr and Vrep, equal to −2.31 kcal mol−1 and 1.13 kcal mol−1, respectively.35 These results also suggests indirectly that in the present case CT plays a negligible role. When an aromatic molecule interacts with a saturated hydrocarbon, dispersion is the major source of attraction, which can be associated with the high value of the polarizability of involved partners.86 However, although the electrostatic contribution is small in comparison with the dispersion one, it plays an important role in determining the geometry of the clusters, because electrostatics exhibits a pronounced dependence on the geometry of approach and then is most responsible for the CH−π interaction directionality. The main advantage of having an analytical potential energy function is the possibility to evaluate, in a straightforward way, stable and nonstable system configurations, as well as to investigate the orientational dependence of the interaction energy. At first, the rotation of CH4 above the aromatic plane has been considered, with the C atom placed always on the C6 rotational axis of Bz (at the same value of the intermolecular distance as that in ref 35) by varying the angle ϕ that one of the CH bonds of methane forms with the C6 rotational axis of Bz (see Figure 3). In particular, seven values going from ϕ = 0° (the CH bond along the rotational axis of Bz) to ϕ = 135°, at regular intervals of 15°, have been analyzed (the same interval that in Figure 5 of ref 35). The agreement between semiempirical predictions and ab initio calculations is very good.35 Approaches from CH4 to Bz, in which the C of CH4 lies far from the C6 axis, have been also investigated. To this purpose, several values of θ, the angle that the r ⃗ vector forms with the C6 rotational axis of Bz (see Figure 4) going from 0° to 90°, have been considered, and results are given in Table 4. It has been found that when θ increases, the distance from the C atom of CH4 to the Bz c.m., associated to

Figure 4. CH4 approaching to Bz with a CH bond forming the θ angle with the z axis. For θ = 90° the CH4 molecule is aligned along a CH bond of Bz.

Table 4. θ Angle (see Figure 4), Distance from the C atom of CH4 to the c.m. of Bz (R, in Å), and Nonelectrostatic (Vnel) and Electrostatic (Vel) Interaction Energiesa θ/deg

R

Vnel

Vel

Vattr

Vrep

Vtotal

0 15 30 45 60 75 90

3.677 3.791 4.097 4.540 5.073 5.530 5.717

−1.18 −1.07 −0.88 −0.70 −0.55 −0.44 −0.40

−0.24 −0.27 −0.27 −0.18 −0.03 0.07 0.10

−2.31 −2.13 −1.85 −1.47 −1.03 −0.73 −0.62

1.13 0.97 0.97 0.77 0.49 0.29 0.22

−1.42 −1.34 −1.15 −0.88 −0.58 −0.37 −0.30

a

The attractive (Vattr) and repulsive (Vrep) components of Vnel (Vnel = Vattr + Vrep) and Vtotal (Vtotal = Vnel + Vel) are also given. All energies are given in kcal mol−1.

the minimum energy structure, also increases. Accordingly, the potential energy becomes more repulsive going from θ = 0° to θ = 90°. As can be seen from the third and fourth columns of Table 4, the increase of θ causes a decrease (less negative values) of both, Vnel and Vel, components. However, for small angles θ ≤ 30°, Vel remains almost unaffected. 5484



Article

3.1. CH4-Bond vs CH4-Effective Atom Representations. The main motivation of using the molecule-effective atom representation of the potential is the reduction of the simulation time, which is desirable to investigate large systems. This reformulation of the potential energy interaction was also used successfully to study the Bz−H2O system.54 Due to the similarities of Bz−H2O and Bz−CH4 it can be expected that such decomposition is also applicable to Bz−CH4. However, the new formulation requires an additional inspection of the energy contributions. A value of Vtotal equal to −1.45 kcal mol−1, in a quite good agreement with the more accurate CH4bond formulation of the interaction, is obtained when the value of β = 8.75 is used together with the parameters of Table 2. In spite of this, by analyzing again the attractive and repulsive contributions it has been observed that the new description, involving new reference centers, enhances the attractive character of the potential energy function, originating an overestimation of the dispersion energy. As before, thanks to the flexibility of the function, this overestimation can be compensated by slightly increasing the value of β. As a matter of fact, the proposed molecule-effective atom representation with β = 9.4 (which is close to the upper limit of the β variation range55) gives results of Vtotal = −1.44 kcal mol−1, Vattr = −2.31 kcal mol−1, Vrep = 1.08 kcal mol−1, and Vel = −0.22 kcal mol−1, in a good agreement with ab initio results.35 The equilibrium distance from the C atom of CH4 to the Bz c.m. is equal to 3.818 Å, which represents an increase of less than 4% in comparison with the value obtained (3.667 Å) using the CH4bond formulation of the interaction. Both values fall in the range of distances calculated for the Bz−CH4 cluster.30,87 Other geometries have been also analyzed with satisfactory results, with binding energies in agreement with the previous formulation, thus justifying the use of the CH4-effective atom components to investigate larger systems.

expansion of the system, simulations assuming a microcanonical ensemble (NVE) have been performed to equilibrate the system at T = 40 K. Additional equilibrations, using as initial configuration that of the last step of the previous run, have been carried out at increasing values of T. In order to compare our predictions with results available in the literature,91 additional NVE trajectories for the (CH4)n system at a density of 0.347 g cm−3 and T = 126.8 K have been also performed. Under these conditions the system exhibits a liquid behavior in good agreement with available data.91 The radial distribution functions (RDFs) calculated at these conditions are shown in Figure 5.

4. MOLECULAR DYNAMICS SIMULATIONS OF THE (CH4)n AND Bz−(CH4)n SYSTEMS Some properties of the (CH4)n and the Bz−(CH4)n systems have been characterized by means of MD simulations performed using the DL-POLY MD program.88 Both, Bz and CH4, as indicated in section 2, have been considered as rigid bodies. A time step of 1 fs has been considered in all simulations, integrating the trajectories up to a final time of 15 ns. The system has been equilibrated along an additional 1 ns, and during the equilibration period, which has been excluded from the statistical analysis at the end of any trajectory, the velocities of all atoms are rescaled to match the input temperature. The electrostatic energy has been calculated by means of the Ewald sum, with a precision of 10−5. It has been observed that the time step chosen is large enough to keep the fluctuations of Etotal (the sum of kinetic and potential energy available for the system at a given temperature) well below 10−5 kcal mol−1. 4.1. (CH4)n System. The predictions for the CH4−CH4 interaction energy using β = 8.5 and parameters given in section 2 are on the right side of experimental determinations.89 Moreover, the reliability of the CH4−CH4 interaction has been further investigated by taking into account an ensemble of 448 molecules of CH4 placed into a cubic box of 30.4553 Å size. The volume of the cubic box and the number of molecules are consistent with the density of methane in liquid phase equal to 0.42262 g cm−3 at 1.013 bar of pressure and 112 K of temperature (boiling point).90 In order to avoid an initial

Figure 5. Atom-effective atom formulation of the interaction with β = 9.4: C−C (top), C−H (middle), and H−H (bottom) RDFs for methane at liquidlike conditions (d = 0.347 g cm−3 and T = 126.8 K).

However, most MD simulations have been performed at constant values of the number of molecules, N; pressure, p; and temperature, T (NpT ensemble). The effect of increasing T has been analyzed at constant values of N and p. Bearing in mind that the experimental values of melting, boiling, and critical temperatures for methane fall in the 90−190 K range,90 temperatures within this range have been considered with particular attention. This procedure allows one to estimate the temperature at which a change of phase is likely to occur.48 After equilibration, for temperature values below 135 K, the predicted density at p = 1.013 bar is higher than 0.42262 g cm−3. In Figure 6, the behavior of both density (top panel) and mean diffusion coefficient (lower panel) values are represented for the temperature range investigated. A gradual decrease of the density is observed when increasing T, and the value of 0.42262 g cm−3 is attained at T ≈ 135 K. In spite of the fact that at the boiling point (T = 112 K) the predicted density is about 9% higher than the experimental value,90 the description of the liquid phase, at least at the investigated conditions, can be considered quite realistic (the potential parameters have been derived for CH4−CH4 gas-phase interactions and only 448 5485



Article

Figure 6. Atom-effective atom formulation of the interaction with β = 9.4: Density (top panel) and mean diffusion coefficients (bottom panel) of methane as a function of the temperature at p = 1.013 bar. Red dashed lines show the linear dependence of both density and mean diffusion coefficients on T (followed only at T values lower than 135 K) and blue dashed lines indicate the nonlinear dependence at higher temperatures.

Figure 7. Atom-effective atom formulation of the interaction with β = 9.4: Density (top panel) and mean diffusion coefficients (bottom panel) of methane as a function of the temperature at p = 45.96 bar.

As was observed for Bz−H2O,54 the two formulations give very similar results also for the Bz−(CH4)n (n > 1−10) aggregates. At a defined value of total energy, Etotal, trajectories have been run by increasing the number of methane molecules. The initial configuration for a given Bz−(CH4)n+1 cluster has been constructed by adding a new CH4 molecule to the Bz− (CH4)n configuration obtained in the last step of the run. Then, trajectories for a specific cluster have been run over increasing values of Etotal, by starting a new equilibration period from the final configuration, velocities, and forces of the previous run. The mean value of both kinetic (and temperature, T) and potential energy contributions, Ekin and Ecfg, respectively, are obtained from the corresponding values monitored over the MD trajectory. This means that Ecfg is defined as the average of the total potential energy over all the accessible configurations of the system for a given value of Etotal (or T). At each mean temperature, the minimum value of Ecfg, the corresponding distance between the C atom of methane and the Bz c.m., and the temperature associated with the step have been monitored. In spite of the fact that the characteristics of the global minima should be established by means of global optimization studies (see, for instance, refs 92, 93), in order to obtain an estimation of the equilibrium energy, the mean values of several energy contributions have been calculated at different values of total energy (or T) and then extrapolated at 0 K. The mean values of the potential energy (Ecfg) and of the corresponding nonelectrostatic (Enel) and electrostatic (Eel) contributions (extrapolated at 0 K) are given in Table 5. Results, as expected, indicate that the interaction is mainly dominated by nonelectrostatic forces. The most stable configuration for the Bz−(CH4)2 system is the one in which the two methane molecules are placed in opposite sides of the aromatic plane (1|1). However, the configuration energy in which the two methane molecules are placed on the same side of the aromatic ring (2|0) only differs by 0.3 kcal/mol. In the (2|0) configuration, one of the methane molecules is displaced from the C6 rotational axis of Bz. At low temperatures, only the initial (1|1) configuration is observed along the trajectory, with the two methane molecules rotating around the CH bond pointing toward the Bz c.m. However, when the number of the CH 4 molecules increases,

molecules are taken into account). The figure shows that at T ≈ 135 K, both the density an the mean diffusion coefficient modifies their dependence on T, which for T > 135 K does not follow the linear dependence observed at lower temperatures. Moreover, it has been observed that for temperatures below 180 K, the proposed temperature value at the beginning of the trajectory is easily reached along the equilibration, while for the higher ones further attempts to increase T cause its decrease, and then, the system suddenly expands (this originates a high decrease of the density, which takes values of the same order of magnitude as methane gas). The previous results of density and diffusion coefficients seem to point out a liquid- to gas-phase transition within the range of investigated temperatures. In order to further analyze the behavior of the ensemble of methane molecules, several MD simulations at a pressure of 45.96 bar, which is the critical point pressure,90 have been performed at increasing values of T, starting from T = 110 K. Results show a density decrease and a diffusion coefficient increase with T, as shown in Figure 7. As can be observed, in general, an increase of 5 K on T causes a small and gradual density decrease, but when T is increased from 190 to 195 K, the variation of the density values becomes more pronounced. The sudden changes in both diffusion coefficients and density could indicate T and p conditions close to those of the critical point (T = 190.45 K and p = 45.96 bar).90 4.2. Bz−(CH4)n Aggregates (n > 1−10). In order to analyze possible competitive and cooperative effects on CH4− CH4 and CH−π interaction, an increasing number n of CH4 molecules have been considered. Bearing in mind the maximum value of n used in our MD simulations, the molecule-effective atom decomposition of the interaction has been used and the results have been carefully tested in some favorable cases. In particular, to ensure that this formulation of the interaction allows one to obtain results similar to those derived from the more accurate atom-bond one, the smaller Bz−(CH4)n (n > 1− 10) clusters have been analyzed, considering a NVE ensemble of particles, using both formulations of the interaction (atombond parameters of Table 1 with β = 8.75 and moleculeeffective atom formulation parameters of Table 2 with β = 9.4). 5486



Article

trajectory. Results of Figure 8 have been obtained from MD simulations at 75 K. On the contrary, the isosurfaces shown on the right-hand side of Figure 8, derived from MD simulations at 30 K for the Bz−(CH4)7 system, clearly indicate a high probability density of the CH4 molecules in the Bz plane. In general, for n > 3, frequent interconversions, with the CH4 molecules moving from one side of the aromatic plane to the other one, are observed. By further increasing the number of methane molecules, most of them cannot be placed in the most favorable position with respect to Bz. This means that the interaction of Bz with such molecules is very weak, being on the same order of magnitude (or even lower) as the CH4−CH4 interaction, favoring the mobility of the molecules and, accordingly, the isomerization processes between different (m| n) structures. 4.3. Bz Molecule Surrounded by Methane. The systems containing a high number of methane molecules have been studied by considering a NpT ensemble of 512 molecules of CH4 surrounding one Bz molecule, and as usual, MD calculations have been performed at p = 1.013 bar and for increasing values of T within the 40−200 K range. As can be expected, the increase of T and, accordingly, of Ekin allows one to explore more repulsive zones of the potential energy surface, promoting an enhancement of Ecfg (less negative values), which is due to variations of both Enel and Eel. However, some differences in the dependence of Enel and Eel with T are observed: while the most important variation of Eel is observed at the lowest temperatures investigated, that of Enel exhibits a different behavior. In Figure 9, Eel is represented as a function

Table 5. Mean Potential Energy (Ecfg) and the Corresponding Nonelectrostatic (Enel) and Elecrostatic Contributions (Eel) for the the Bz−(CH4)n (n > 1−10) Aggregates Obtained Extrapolating Etotal at T = 0 Ka

a

n

Ecfg

Enel

Eel

2 3 4 5 6 7 8 9 10

−2.79 −4.20 −5.34 −6.93 −8.56 −9.85 −12.15 −13.80 −15.28

−2.35 −3.57 −4.45 −5.85 −7.33 −8.40 −10.48 −11.94 −13.31

−0.44 −0.63 −0.89 −1.08 −1.23 −1.45 −1.67 −1.86 −1.97

All energies are given in kcal mol−1.

interconversion between different structures becomes frequent, even at low temperatures. Due to the increase (more negative) of the CH4−CH4 interaction contribution with the number of molecules, for Bz−(CH4)3, the energy difference between the most stable (2|1) and (3|0) configurations is very low (≈ 0.2 kcal/mol). However, the relatively large value of the distance between two methane molecules placed on opposite sides of the aromatic ring makes their interaction difficult and, accordingly, the (2|1) to (3|0) isomerization is unlikely. In the (2|1) configuration, the two molecules placed on the same side of the Bz plane occupy similar positions as in the (2|0) one, while in the (3|0) configuration the three molecules tend to be placed in such a way as to form an equilateral triangle. The changes originated by an increase of the number of CH4 molecules have been also studied by considering the threedimensional (3D) probability density, evaluated by introducing our simulation data (properly shifted to the inertial reference frame) into the VMD visualization program.94 The Volmap tool, included in the mentioned package, allows one to visualize two-dimensional (2D) slices from a volumetric data set and isosurfaces of (mass weighted) constant probability density for the methane molecules. As an example, on the left-hand side of Figure 8 a two-dimensional map of the probability density of CH4 molecules around Bz, for the Bz−(CH4)3 system, is shown. As can be seen in the figure, the isosurfaces placed on opposite sides of the aromatic plane are not connected, so that the probability density for CH4 in the Bz plane is null, indicating that no isomerizations takes place along the

Figure 9. Atom-effective atom formulation of the interaction with β = 9.4: Eel as a function of T for the system formed by benzene surrounded by 512 methane molecules.

Figure 8. Atom-effective atom formulation of the interaction with β = 9.4: Two-dimensional maps of the probability density of the CH4 molecules around benzene: the case of 3 (T = 74 K) and 7 (T = 30 K) CH4 molecules in panels a and b, respectively. 5487



Article

Figure 10. Atom-effective atom formulation of the interaction with β = 9.4: Two-dimensional maps of the probability density of the CH4 molecules around benzene at T = 160 K, with x, y, and z axis on the right, middle, and left, respectively. The rotational axis of Bz is the z axis.

of T showing that such an energetic term tends to a near constant value (estimated to be ≈ −3.5 kcal mol−1). Bearing in mind that at 177 K, for the ensemble of 448 molecules of CH4, a value of Eel = −2.69 kcal mol−1 has been found, the plateau observed in Figure 9 should correspond to temperatures close to a transition from liquid to gas phase. Moreover, MD simulations reveal some particularities on the values of the mean diffusion coefficients of benzene and methane as a function of T. At temperatures below the melting point of methane, the Bz mean diffusion coefficient is slightly larger than that of methane but on the same order of magnitude, while in the 100−140 K temperature range the coefficients are nearly identical, indicating that at temperatures lower than 140 K, the benzene is transported by methane along its spreading. By further increasing T, the methane diffusion coefficient rises and that of benzene diminishes and the corresponding mean diffusion coefficients become very different at the higher T investigated, consistent with their different mass. These facts could indicate that the benzene and methane movements seem to be independent. In particular, at 180 K methane behaves as in the gas phase, with a mean diffusion coefficient of 22.9 × 10−9 m2/s, in good agreement with results shown in section 4.1. It is interesting to note that the dependence of both diffusion coefficients and Eel with T seems to confirm once more a methane transition from liquid to gas phase. The distribution of CH4 around Bz has been also investigated. It has been found that, accordingly with the characteristics of the Bz−CH4 interaction, methane molecules tend to occupy axial positions (close to the rotation axis) more than Bz plane positions, as can be appreciated in Figure 10, where the three x, y, and z (axial) two-dimensional (2D) slices of volumetric data are represented on the left, middle, and right sides, respectively.

of the total energy to ab initio results indicates that the representation of the Vattr and Vrep energy components can be improved by varying, within a limited range, the value of the β parameter. To this scope, ab initio calculations are of fundamental importance and allow one to investigate the potentiality of the model and to improve its predictions. In particular, the possibility of comparing different energy contributions with reliable results is of fundamental importance to improve and to extend the formulation of our method to more complex systems. The present formulation of the potential has been used to analyze the dependence of the intermolecular interaction energy on the relative orientation of CH4 and Bz molecules. Results seem to indicate that CT effects play a minor role in Bz−CH4 interaction. The agreement between the prediction of the electrostatic energy and the ab initio values confirms the validity of its representation in terms of punctual charges distributed on the molecular frame, the charge distribution on the Bz molecule being consistent with its quadrupole moment and that on CH4 with its octupole moment. Bearing in mind the importance of extending the present formulation of the potential to study large systems, a careful comparison of the CH4 bond and the CH4 effective atom descriptions has been performed, and the simpler formulation has been used to carry out MD simulations of the Bz−(CH4)n clusters, which have allowed investigation of the competitive role of isomerization and dissociation processes when increasing both the number of molecules surrounding Bz and the temperature of the aggregate. Isomerization processes have been mainly observed by considering small Bz−(CH4)n (n = 1−10) aggregates. In spite of the simplicity of the model, the study of the methane system at increasing values of T seems to indicate that the potential model, even using an ensemble of only about 500 molecules, is able to describe also some macroscopic characteristics of methane in the condensed phase.

■

5. CONCLUSIONS The potential model based on the decomposition of molecular polarizabilities has been exploited to investigate the role of CH−π interaction in the Bz−CH4 cluster. The sum of 12 molecule bond components (six CH4−CC and six CH4−CH), each one represented by the ILJ function, has been used to describe the nonelectrostatic component, Vnel, of the total interaction energy. The relevant parameters of the ILJ functions, having a physical meaning, have been derived from the molecular polarizability of the methane molecule and from the polarizability associated with the Bz bonds. It has been found that the present formulation gives binding energies in agreement with high quality ab initio results. However, the comparison of the values associated with various components

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS M.A., A.A., and J.M.L. acknowledge financial support from the Ministerio de Educación y Ciencia (Spain, Project CTQ201061709) and to the Comissionat per a Universitats i Recerca del DIUE (Generalitat de Catalunya, Project 2009-SGR 17). Thanks are also due to the Centre de Serveis Científics i Acadèmics de Catalunya CESCA and Fundació Catalana per a 5488



Article

(34) Tsuzuki, S.; Honda, K.; Fujii, A.; Uchimaru, T.; Mikami, M. Phys. Chem. Chem. Phys. 2008, 10, 2860−2865. (35) Tsuzuki, S.; Fujii, A. Phys. Chem. Chem. Phys. 2008, 10, 2584− 2594. (36) Fujii, A.; Shibasaki, K.; Kazama, T.; Itaya, R.; Mikami, N.; Tsuzuki, S. Phys. Chem. Chem. Phys. 2008, 10, 2836−2843. (37) Gord, J. R.; Garrett, A. W.; Bandy, R. E.; Zwier, T. S. Chem. Phys. Lett. 1990, 171, 443−450. (38) Gotch, A. J.; Pribble, R. N.; Ensminger, F. A.; Zwier, T. S. Laser Chem. 1994, 13, 187−205. (39) Plevin, M. J.; Bryce, D. L.; Boisbouvier, J. it Nature Chemistry 2010, 2, 466−471. (40) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Fujii, A. J. Phys. Chem. A 2006, 110, 10163−10168. (41) Umezawa, Y.; Nishio, M. Bioorg. Med. Chem. 1998, 6, 493−504. (42) Umezawa, Y.; Nishio, M. Bioorg. Med. Chem. 1998, 6, 2507− 2515. (43) Albertí, M.; Solé, A.; Aguilar, A. J. Chem. Soc., Faraday Trans. 1991, 87, 37−44. (44) Albertí, M.; Pietro, M.; Aguilar, A. J. Chem. Soc., Faraday Trans. 1992, 88, 1615−1619. (45) Albertí, M.; Sayós, R.; Solé, A.; Aguilar, A. J. Chem. Soc., Faraday Trans. 1991, 87, 1057−1068. (46) Bargueño, P.; Jambrina, P. G.; Alvariño, J. M.; Hernández, M. L.; Aoiz, F. J.; Menéndez, M.; Verdasco, E.; González-Lezana, T. J. Phys. Chem. A 2009, 113, 14237−14250. (47) Martínez, T.; Hernández, M. L.; Alvariño, J. M.; Laganà, A.; Aoiz, F. J.; Menéndez, M.; Verdasco, E. Phys. Chem. Chem. Phys. 2000, 2, 589−597. (48) Albertí, M.; Pirani, F. J. Phys. Chem. A 2011, 115, 6394−6404. (49) Gough, K. M. J. Chem. Phys. 1989, 91, 2424−2432. (50) Israelachvili, J. N. In Intermolecular and Surface Forces; Academic Press: London, 1991; Chapter 5. (51) Ghanty, T. K.; Ghosh, S. K. J. Phys. Chem. 1996, 100, 17429− 17433. (52) Denbigh, K. G. Trans. Faraday. Soc. 1940, 36, 936−948. (53) Smith, R. P.; Mortensen, E. J. J. Chem. Phys. 1960, 32, 502−507. (54) Albertí, M.; Faginas Lago, N.; Pirani, F. Chem. Phys. 2012, 399, 232−239. (55) Pirani, P.; Brizi, S.; Roncaratti, L. F.; Casavecchia, P.; Cappelletti, D.; Vecchiocattivi, F. Phys. Chem. Chem. Phys. 2008, 10, 5489−5503. (56) Faginas Lago, N.; Huarte Larrañaga, F.; Albertí, M. Eur. Phys. J. D 2009, 55, 75−85. (57) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. J. Phys. Chem. A 2010, 114, 11964−11970. (58) Albertí, M.; Aguilar, A.; Lucas, J. M.; Laganà, A.; Pirani, F. J. Phys. Chem. A 2007, 111, 1780−1787. (59) Huarte-Larrañaga, F.; Aguilar, A.; Lucas, J. M.; Albertí, M. J. Phys. Chem. A 2007, 111, 8072−8079. (60) Albertí, M.; Castro, A.; Laganà, A.; Moix., M; Pirani, F.; Cappelletti, D. Eur. Phys. J. D 2006, 38, 185−191. (61) Huarte-Larrañaga, F.; Aguilar, A.; Lucas, J. M.; Albertí, M. Theor. Chem. Acc. 2011, 128, 757−767. (62) Albertí, M.; Huarte-Larrañaga, F.; Aguilar, A.; Lucas, J. M.; Pirani, F. Phys. Chem. Chem. Phys. 2011, 13, 8251−8258. (63) Albertí, M.; Aguilar, A.; Lucas, J. M.; Cappelletti, D.; Laganà, A.; Pirani, F. Chem. Phys. 2006, 328, 221−228. (64) Albertí, M.; Aguilar, A.; Bartolomei, M.; Cappelletti, D.; Laganà, A.; Lucas, J. M.; Pirani, F. Lect. Notes Comput. Sci. 2008, 5072, 1026− 1035. (65) Albertí, M.; Aguilar, A.; Cappelletti, D.; Laganà, A.; Pirani, F. Int. J. Mass. Spectr. 2009, 280, 50−56. (66) Albertí, M.; Aguilar, A.; Bartolomei, M.; Cappelletti, D.; Laganà, A.; Lucas, J. M.; Pirani, F. Phys. Script. 2008, 78, 058108(1)− 058108(7). (67) Paolantoni, M.; Faginas Lago, N.; Albertí, M.; Laganà, A. J. Phys. Chem. A 2009, 113, 15100−15105.

la Recerca for the allocated supercomputing time. M.A. also thanks the Ministerio de Educación y Ciencia for the mobility project PR2010-0243. F.P. acknowledges financial support from the Italian Ministry of University and Research (MIUR) for PRIN Contracts.

■

REFERENCES

(1) Chandra Dey, R.; Seal, P.; Chakrabarti, S. J. Phys. Chem. A 2009, 113, 10113−10118. (2) Müller-Dethelfs, K.; Hobza, P. Chem. Rev. 2000, 100, 143−167. (3) Alonso, J. L.; Antolínez, S.; Bianco, S.; Lesarri, A.; L’opez, J. C.; Caminati, W. J. Am. Chem. Soc. 2004, 126, 3244−3249. (4) Cabarcos, O. M.; Weinheimer, J. C.; Lisy, J. M. J. Chem. Phys. 1998, 108, 5151−5154. (5) Cabarcos, O. M.; Weinheimer, J. C.; Lisy, J. M. J. Chem. Phys. 1999, 110, 8429−8435. (6) Morais-Cabral, J. H.; Zhou, Y.; MacKinnon, R. Nature 2001, 414, 37−42. (7) Kumpf, R. A.; Dougherty, D. A. Science 1993, 261, 1708−1710. (8) Dougherty, D. A. Science 1996, 271, 163−168. (9) Ma, J. C.; Dougherty, D. A. Chem. Rev. 1997, 97, 1303−1324. (10) Jalbout, A. F.; Adamovicz, L. J. Chem. Phys. 2002, 116, 9672− 9676. (11) Tsuzuki, S.; Yoshida, M.; Uchimaru, T.; Mikami, M. J. Phys. Chem. A 2001, 105, 769−773. (12) Albertí, M.; Castro, A.; Laganà, A.; Moix, M.; Pirani, F.; Cappelletti, D.; Liuti, G. J. Phys. Chem. A 2005, 109, 2906−2911. (13) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. Theor. Chem. Acc. 2009, 123, 21−27. (14) Albertí, M.; Castro, A.; Laganà, A.; Pirani, F.; Porrini, M.; Cappelletti, D. Chem. Phys. Lett. 2004, 392, 514−520. (15) Pirani, F.; Albertí, M.; Castro, A.; Moix Teixidor, M.; Cappelletti, M. Chem. Phys. Lett. 2004, 394, 37−44. (16) Albertí, M. J. Phys. Chem. A 2010, 114, 2266−2274. (17) Siglow, K.; Neuhauser, R.; Jürgen Neusser, H. J. Chem. Phys. 1999, 110, 5589−5599. (18) Hobza, P.; Bludsky, O.; Selzle, H. L.; Schlag, E. W. Chem. Phys. Lett. 1996, 250, 402−408. (19) Brupbacher, Th.; Makarewicz, J.; Bauder, A. J. Chem. Phys. 1994, 101, 9736−9746. (20) Duan, G.; Smith, V. H., Jr.; Weaver, D. F. Mol. Phys. 2001, 99, 1689−1699. (21) Tauer, T. P.; Derrick, M. E.; Sherrill, C. D. J. Phys. Chem. A 2005, 109, 191−196. (22) Hermida-Ramón, J. M.; Cabaleiro-Lago, E. M.; RodríguezOtero, J. J. Chem. Phys. 2005, 122, 204315(1)−204315(5). (23) Cabaleiro-Lago, E. M.; Rodríguez-Otero, J.; Peña-Gallego, A. J. Chem. Phys. 2008, 129, 084305(1)−084305(8). (24) Grabowski, S. J.; Leszczynski, J. Chem. Phys. 2009, 355, 169− 176. (25) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2000, 122, 3746−3753. (26) Reimann, B.; Buchhold, K.; Vaupel, S.; Brutschy, B.; Havlas, Z.; Špirko, V.; Hobza, P. J. Phys. Chem. A 2001, 105, 5560−5566. (27) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Phys. Chem. A 2002, 106, 4423−4428. (28) Fujii, A.; Morita, S.; Miyazaki, M.; Ebata, T.; Mikami, N. J. Phys. Chem. A 2004, 108, 2652−2658. (29) Munshi, P.; Row, T. N. G. J. Phys. Chem. A 2005, 109, 659−672. (30) Ringer, A. L.; Figgs, M. S.; Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2006, 110, 10822−10828. (31) Shibasaki, K.; Fujii, A.; Mikami, N.; Tsuzuki, S. J. Phys. Chem. A 2006, 110, 4397−4404. (32) Tekin, A.; Jansen, G. Phys. Chem. Chem. Phys. 2007, 9, 1680− 1687. (33) Shibasaki, K.; Fujii, A.; Mikami, N.; Tsuzuki, S. J. Phys. Chem. A 2007, 111, 753−758. 5489



Article

(68) Albertí, M.; Faginas Lago, N.; Laganà, A.; Pirani, F. Phys. Chem. Chem. Phys. 2011, 13, 8422−8432. (69) Cambi, R.; Cappelletti, D.; Liuti, G.; Pirani, F. J. Chem. Phys. 1991, 95, 1852−1861. (70) Pirani, F.; Maciel, G. S.; Cappelletti, D.; Aquilanti, V. Int. Rev. Phys. Chem. 2006, 25, 165−199. (71) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Cappelletti, D.; Coletti, C.; Re, N. J. Phys. Chem. A 2006, 110, 9002−9010. (72) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Coletti, C.; Re, N. J. Phys. Chem. A 2009, 113, 14606−14614. (73) Li, S.; Cooper, V. R.; Thonhauser, T.; Puzder, A.; Langreth, D. C. J. Phys. Chem. A 2008, 112, 9031−9036. (74) Ma, J.; Alfè, Michaelides, A.; Wang, E. J. Chem. Phys. 2009, 130, 154303(1)−154303(6). (75) Feller, D. J. Phys. Chem. A 1999, 103, 7558−7561. (76) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2000, 122, 11450−11458. (77) Kitaura, K.; Morokuma, K. Int. J. Quantum Chem. 1976, 10, 325−340. (78) Kitaura, K.; Morokuma, K. Chemical Applications of Electrostatic Potentials; P., Politzer, D. G., Truhlar, Eds.; Plenum Press: New York, 1981. (79) Chen, W.; Gordon, M. S. J. Phys. Chem. 1996, 100, 14316− 14328. (80) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. Rev. 1994, 94, 1887−1930. (81) Lillestolen, J. C.; Wheatley, R. J. J. Phys. Chem. A 2007, 111, 11141−11146. (82) Wheatley, R. J.; Lillestolen, T. C. Mol. Phys. 2008, 106, 1545− 1556. (83) Ewig, C. S.; Waldman, M.; Maple, J. R. J. Phys. Chem. A 2002, 106, 326−334. (84) Gavezotti, A. J. Phys. Chem. B 2003, 107, 2344−2353. (85) Ramos, C.; Winter, P. R.; Sterans, J. A.; Zwier, T. S. J. Phys. Chem. A 2003, 107, 10280−10287. (86) Tsuzuki, S.; Uchimaru, T.; Tanabe, K. J. Mol. Struct. THEOCHEM 1994, 307, 107−118. (87) Sherrill, C. D.; Takatani, T.; Hohenstein, E. G. J. Phys. Chem. A 2009, 113, 10146−10159. (88) http://www.cse.scitech.ac.uk/software/DL-POLY/ (89) Reid, B. P.; O’Loughlin, J.; Sparks, R. K. J. Chem. Phys. 1985, 83, 5656−5662. (90) http://encyclopedia.airliquide.com/encyclopedia.asp (91) Hernández de la Peña, L.; van Zon, R.; Schofield, J. J. Chem. Phys. 2007, 126, 074106(1)−074106(12). (92) Schulz, F.; Hartke, B. ChemPhysChem 2002, 3, 98−106. (93) Llanio-Trujillo, J. L.; Marques, J. M. C.; Pereira, F. B. J. Phys. Chem. A 2011, 115, 2130−2138. (94) Humphrey, W.; Dalke, A.; Schulten, K. VMD-Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38.

5490


Competitive Role of CH4–CH4 and CH− π Interactions in C6H6–(CH4

Recommend Documents