ARTICLE pubs.acs.org/JPCA
Features of Ar Solvation Shells in Neutral and Ionic Clustering: The Competitive Role of Two-Body and Many-Body Interactions Margarita Albertí*,† and Fernando Pirani‡ † ‡
IQTCUB, Departament de Química Física, Universitat de Barcelona, Barcelona, Spain Dipartimento di Chimica, Universita di Perugia, Perugia, Italy ABSTRACT: The semiempirical methodology, introduced to describe noncovalent intermolecular interactions in atom/ionmolecule systems, is here extended to investigate a prototype cluster, formed by benzene (Bz) and closed-shell ions (Naþ and/or Cl), surrounded by neutral species (Ar), forming solvation shells. The involved multidimensional potential energy surface (PES) is assumed to depend on a critical balancing of some effective interaction components. In particular, for the Ar solvated BzNaþCl system, the nonelectrostatic component of the total interaction has been formulated as a combination of two-, three-, and four-body contributions, each one represented by a proper function, with the fourbody and part of the three-body terms arising from nonadditive induction effects. The proposed formulation, in which the induction is included both implicitly and explicitly, ensures the accurate description of all dissociation channels, leading to simpler clusters and/or pure solvent. Some properties of the solvent, represented by an ensemble of 500 Ar atoms, have been analyzed by performing molecular dynamics simulations at several temperatures. The obtained results have been found to be consistent with experimental observations. In order to investigate propensities, similarities, and differences in the competing clusters, the Ar solvation shells of Bz, BzNaþ, BzCl and BzNaþCl have been characterized. The inspection of the solvation shell of Bz allows one to distinguish between groups of Ar atoms occupying positions on and out of the plane defined by the aromatic ring. Regarding the solvation shells of BzNaþ and BzCl, it has been observed that they are strongly affected by the most stable structures of the unsolvated systems. However, BzNaþ shows more compact solvation shells than BzCl. Finally, important asymmetries, basically promoted by the additional many-body induction effects on the solvent atoms, have been observed in the solvation shells of BzNaþCl.
1. INTRODUCTION The study at the molecular level of the influence of intermolecular forces on static and dynamics properties of clusters is a topic of interdisciplinary relevance. In particular, the study of molecular aggregates, often stabilized by weak noncovalent intermolecular forces, which, despite their weakness, affect the structural properties of both simple and complex systems and the dynamics of several elementary processes, is of great interest in chemistry, physics, and biology.120 Noncovalent intermolecular interactions control basic phenomena, such as molecular recognition and selection processes,2123 the formation of hydrogen bonds,1,24 the competitive solvation of ions by different partners,2426 and material design rationalization.2729 However, the characterization of the noncovalent interaction role, when competing forces are present, is difficult,30 and in most cases, a coherent and accurate formulation of the involved force field is nowadays unavailable, at least at the level of accuracy and detail required in molecular dynamics (MD) simulations, in which the system explores equilibrium and nonequilibrium regions of the phase space. As it has been indicated before,31 the use of an ab initio potential energy surface (PES) often implies the interpolation of calculated points using spline functions or fitting them to an adequate function, and several studies, such as a global minimum search of atomic and moleculer clusters31,32 or the r 2011 American Chemical Society
analysis of solvation processes, are too computer demanding when the most accurate ab initio potentials are used. Therefore, accurate semiempirical methods, adopting suitable analytical functions to properly describe the potential energy in the whole configuration space, are of fundamental importance. Moreover, in order to ensure an affordable computational cost, the analytical functions should be as simple as possible while ensuring the accurate assessment of the strength and relative role played by the various components of the interaction. However, this accurate assessment is a difficult task,33 and experiments and/or high-level ab initio calculations are crucial to ensure quality and reliability of the adopted analytical representations. Moreover, the combination of semiempirical methods with experimental and ab initio data allows one to perform critical comparisons and to obtain information useful to improve the semiempirical models. In the last years, we have developed a semiempirical method based on the decomposition of the molecular polarizability,34 whose reliability has been probed by means of high-level ab initio calculations3538 and high-resolution scattering data.39 The model has been adopted to study several systems, mainly involving aromatic Received: March 31, 2011 Revised: May 11, 2011 Published: May 31, 2011 6394
dx.doi.org/10.1021/jp202995s | J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
V ¼ VBzNaþ Cl þ
n1 n
n
∑ ∑ VAr Ar þ i∑¼ 1 VAr ðBzNa Cl Þ i¼1 j > i i
j
i
þ
ð1Þ
0
VBzNaþ Cl ¼ VBzNaþ þ VBzCl þ VNaþ Cl þ Vind
ð2Þ
where both VBzNaþ and VBzCl combine independent nonelectrostatic (Vnel) and electrostatic (Vel) interaction components.3335,4042,54 Vnel represents effective potential interactions that implicitly include the induction effects operative in the BzNaþ and BzCl dimers. On the contrary, for convenience, VNaþCl in eq 2 is not further decomposed. Accordingly, VNaþCl is also an effective term but including, in this case, nonelectrostatic, electrostatic, and induction contributions.46 This means that V0ind in eq 2 does not represent the total induction component but only those additional effects, explicitly considered, arising from the simultaneous presence of the two ions.46 This approach, in which an important part of induction is incorporated in the effective VBzNaþ, VBzCl, and VNaþCl potential functions, ensures the correct partition of the interaction energy in the binary compounds when the ternary one dissociates (crucial in MD simulations). The Vnel component, associated with either VBzNaþ or VBzCl, is described by applying the concept of molecular bond polarizability additivity,55 for which the polarizability of the Bz molecule can be decomposed in 12 ellipsoids of polarizability associated with the Bz chemical bonds. Then, for VBzNaþ (or VBzCl), Vnel is described as a sum of 12 individual ionbond contributions (see, for instance, refs 35, 37, and 38), each one represented by an improved Lennard-Jones (ILJ) function,34,39 whose form is " m r0 ðγÞ nðr, γÞ VILJ ðr, γÞ ¼ εðγÞ nðr, γÞ m r m nðr, γÞ r0 ðγÞ ð3Þ nðr, γÞ m r with the n(r,γ) exponent expressed as 2 r nðr, γÞ ¼ β þ 4:0 r0 ðγÞ
ð4Þ
The relevant transferable potential parameters,56 ε and r0, describing the well depth and its location for any ionbond pair are modulated by means of a simple trigonometric weighted sum of the respective parallel and perpendicular values.33 The weights depend on the angle of approach, γ, of the ion to a given bond (see panel (a) of Figure 1). For each, an ionbond pair has been considered with m = 4. It must be indicated that VNaþCl has been also represented by eq 3, after removing the angular dependence (ε and r0 no longer depend on γ) and assuming m = 1. As stressed above, V0ind arises from the nonadditive behavior of induction, and its expression has been isolated by subtracting the induction energies of the isolated BzNaþ and BzCl binary compounds from the total induction. Then, each bond contribution to Vind0 is given again as a combination of its parallel, V0||,ind, and perpendicular, V0^,ind, components, defined as46 q1 q2 0 V , ind ¼ 2 2 cos γ1 cos γ2 R ð5Þ r1 r2 6395
)
2. POTENTIAL ENERGY SURFACE The total intermolecular potential energy, V, of Bz NaþCl solvated by n Ar atoms can be decomposed as
where the first term, corresponding to the intermolecular interaction of the isolated ternary aggregate, includes Bzion (VBzNaþ and VBzCl) and ionion (VNaþCl) interaction contributions and additional four-body induction effects (V0ind).46 VBzNaþCl is then expressed as
)
molecules. In particular, it has been applied to investigate benzenealkaline ion (BzMþ) and benzenehalide ion (BzX) systems,33,4044 even including two Bz molecules.45 The previous results encouraged the extension of the methodology to more complex systems, such as the BzMþX aggregates, for which the inclusion of additional induction effects due to the simultaneous presence of the two ions was carefully taken into account.46 The interaction energy, predicted by the model at selected geometries, was found to be in good agreement with results and trends of ab initio calculations.47 Moreover, the induction energy was introduced in such a way to predict the correct value of the potential energy in all dissociation limits, thus opening the possibility of studying the ternary compound by means of MD simulations, even at temperatures where the dissociation probability is high. The inclusion of induction effects appears to be necessary to investigate several processes, such as the ion selectivity occurring in biological channels, which often depends on three-body interactions involving ions of different sign,4852 for instance, in the ammonia channel where the interaction between the cation and negatively charged amino acid residues plays a key role.47 Research activity on molecular aggregates involving aromatic and cyclic molecules, considered prototypical systems to investigate important processes, some of them indicated above, has significantly grown.1123,53 In particular, BzMþX aggregates are good candidates to evaluate the mentioned competing effects between interaction forces of different nature. Two of the related fragments, BzMþ and BzX, involve the individual and selective roles of cationπ36 and anionπ37 interactions, which asymptotically are mainly governed by ionmolecular quadrupole interactions, while the third fragment, MþX, is dominated by ionion interactions. Moreover, a further assessment can be favored when the investigation involves a rare gas (Rg)-solvated environment, in which the RgRg interaction exclusively depends on van der Waals forces. In the present investigation, MD simulations of the BzNaþCl aggregate and of all dissociation products in a solvation environment represented by an ensemble of about 500 Ar atoms previously thermalized are carried out. Due to the high number of Ar atoms, the solvation shells of Bz, BzNaþ, BzCl, and BzNaþCl have been calculated by representing the BzAr interaction as a combination of effective atomatom contributions rather than of independent atombond pair potentials.34 Moreover, in the solvated ternary aggregate, further three-body (non additive) induction effects, arising because of the simultaneous presence of argon atoms and two ions, have been explicitly considered by including a new term on the potential energy representation. The paper is structured as follows; in section 2, we outline the formulation of the semiempirical PES, and in section 3, the results, including the investigation of some properties of pure solvent and of solvated Bz, BzNaþ, BzCl, and BzNaþCl aggregates, are analyzed. Concluding remarks are given in section 4.
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Table 4. Well Depth (ε), Equilibrium Distance (r0), and β and m Parameters Defining the Two-Body ArAr, ArNaþ, ArCl, ArC, and ArH Interactions pairs ArAr ArNaþ ArCl
Figure 1. (a) Ionbond potential variables used to describe Vnel by means of eq 3. (b) Variables used in the description of nonadditive induction effects in eqs 5 and 6.
Table 1. Perpendicular and Parallel Components of the Well Depth (ε^, ε||) and of the Equilibrium Distances (r0^, r0||) for the Different IonBond Pairs and β and m Parameters ε^/meV
ε||/ meV
Na CC
33.01
102.20
2.848
3.149
8.5
4
NaþCH
62.15
62.73
2.601
2.808
8.5
4
ClCC
16.37
59.64
3.832
4.073
7.0
4
ClCH
25.48
28.60
3.655
3.839
7.0
4
ionbond þ
r0^/ Å
r0||/ Å
β
m
ArC ArH
ε/meV
r0/Å
β
m
12.370
3.757
9.0
6
2.730
8.9
4
3.710
7.8
4
3.982 3.524
9.0 9.0
6 6
167.0 64.87 5.668 3.611
neutralneutral interactions compared to those of the ionneutral ones. However, the new atomatom parameters have been selected to reproduce the main characteristics provided by atombond formulation of the BzAr interaction. Moreover, to explicitly account for nonadditive induction effects, caused by the simultaneous presence of a cation and anion in the proximity of Ar atoms, the corrective three-body term V00ind has been included in the formulation. This procedure, followed to investigate the unsolvated NaþBzCl ionic trimer,46 is also applied here to obtain again an accurate description of all dissociation limits of the solvated system. Accordingly, the interaction of one Ar atom with NaþBzCl is expressed as 00
Table 2. Well Depth (ε), Equilibrium Distance (r0), and β and m Parameters Defining the NaþCl Interactions MþX
ε/meV
r0/Å
β
m
NaþCl
5740
2.32
8.0
1
Table 3. Parallel, r||, and Perpendicular, r^, Components of the CC and CH Bond Polarizabilities33
a
R||/Å3 a
R^/Å3 a
CC
1.78
0.38
CH
0.67
0.49
1 Å3 = 6.748 au.
0
V^, ind ¼
q1 q2 sin γ1 sin γ2 cos ωR^ r1 2 r2 2
ð6Þ
where q1 and q2 are the charges of the two ions and R and R^ the parallel and perpendicular CC and CH bond polarizabilities, respectively. The variables r1, r2, γ1, and γ2 are represented on the panel (b) of Figure 1, and ω is the dihedral angle that the plane defined by the cation and a given bond forms with the plane defined by the anion and the same bond. All parameters involved in the previous equations are given in Tables 13, and they have been obtained and tested, as discussed in refs 35, 38, and 46. In the present paper, the interactions in which Ar is involved are described, when it is possible, by means of pair potentials. Here, the usual atombond formulation of the BzAr interaction, exploiting ArCC and ArCH pair potentials,34,54 has been substituted by effective atomatom (ArC and ArH) terms. This choice, motivated by the large number of Ar atoms considered, can be justified by taking into account the less direct role of the bond polarizability tensors on the basic features of
VArðBzNaþ Cl Þ ¼ VArC þ VArH þ VArNaþ þ VArCl þ Vind
ð7Þ All of the involved interactions associated with the ArAr, ArNaþ, ArCl, ArC, and ArH pairs have also been described by means of the ILJ function (eq 3) after removing the angular dependence of ε and r0 and considering m = 6 and 4 for the neutralneutral and neutralion pairs, respectively. The used parameters, reported in Table 4, have been obtained following the guidelines reported before.5761 The reliability of the adopted ionAr pair potentials was discussed before by a direct comparison of the model predictions with results of theoretical calculations and/or experiments.39,60 Moreover, the adopted ArAr potential properly reproduces detailed scattering and spectroscopic observables.39 The V00ind induction energy contribution is attractive for insertions of Ar atoms between the ionic ClNaþ pair and is repulsive for approaches of Ar atoms by the side of Naþ and Cl. V00ind depends on the ion charges, q1 and q2, on the Ar polarizability, RAr, on the distances of the Ar atoms to the cation and the anion, R1 and R2, respectively, and on the angle, θ, formed by R1 and R2 vectors, and it is represented by means of the following equation 00
Vind ¼
q1 q2 RAr cosðθÞ R1 2 R2 2
ð8Þ
In the present study, the value of RAr = 1.64 Å3 (11.07 au) has been considered. According with the previous considerations, the PES is represented by a sum of two-body (ArAr, ArNaþ, ArCl, ArC, ArH, NaþCl), three-body (NaþCC, NaþCH, ClCC, ClCH, ArNaþCl), and four-body (NaþCCCl, NaþCHCl) nonelectrostatic components plus the electrostatic contributions originating from the interaction of each ion with the charge distribution on the benzene molecular frame (see, for instance, refs 33, 35, and 38). 6396
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Configuration potential energy of the ensemble of 500 Ar atoms evaluated in the temperature range from 75 to 145 K. Green and red lines are used to highlight the behavior changes.
Figure 2. The Ar density (top panel) and self-diffusion coefficient (lower panel) as a function of the temperature, T. Green and red lines are used to highlight the behavior changes.
3. RESULTS AND DISCUSSION The present investigation focuses on some basic features arising from the Ar solvation of Bz monomer, of BzNaþ and BzCl dimers, and of BzNaþCl trimer, for which the intermolecular potential described in the previous section has been used. MD simulations have been performed considering a NpT ensemble of particles placed inside of a cubic box bearing periodic boundary conditions. All simulations were performed using the DL_POLY program,62 in which the analytical PES and its derivatives were implemented. A time step of 0.1 fs for a total simulation time of 2 ns was used. To handle the geometry constraints and to integrate the equations of motion, the SHAKE method and the Verlet algorithm were used, respectively. The total intermolecular potential energy was calculated considering a cutoff radius of the half-length of the simulation box, and the electrostatic energy, when present, was evaluated using the Ewald sum.63 At first, an ensemble of 500 Ar atoms, placed inside of the cubic box of 27.2 Å of size, was thermalized. In order to investigate the solvation of the Bz, BzNaþ, BzCl, and BzNaþCl systems, 15 Ar atoms placed on the center of the box were removed and replaced by the corresponding unsolvated system. This means that Bz, BzNaþ, BzCl, and BzNaþCl are solvated by a total of 485 Ar atoms. Further thermalizations of the different solvated systems were also performed. Each simulation provides the values of the energy contributions and of all atom positions along the trajectory, from which the mean values of the energy contributions (represented by Ei) and of the spatial distributions can be derived. Characteristics of the Solvent. MD calculations of an ensemble of 500 Ar atoms have been performed at an initial temperature (T) of 5 K. Then, by using as an initial configuration the one obtained in the last step of the previous simulation, runs at increasing T have been carried out. Bearing in mind that the experimental values of melting, boiling, and critical temperatures
fall in the 80155 K range,64 temperatures within this range have been considered with particular attention. MD results show that the predicted density values of Ar at T e 145 K are of the same order of magnitude as the experimental density of the liquid (d = 1.394 g cm3 at T = 87.28 K and p = 1 atm). Moreover, a linear decrease of the density when T increases (T e 120 K), followed by a sharper decrease of such a property at higher temperatures (T > 120 K), has been observed. Interestingly, the change in the slope has not only been observed for the density values but also for the mean diffusion coefficient. As can be seen in Figure 2, the behavior of both the density (top panel) and mean diffusion coefficient (lower panel) change at about the same temperature. Specifically, at T ≈ 120 K, a sharp decrease of the density is observed (see the top panel of Figure 2), and at the same temperature, the mean diffusion coefficient modifies its dependence on T, which for T > 120 K does not follow the linear dependence on T observed at lower temperatures. Moreover, MD simulations indicate that at the 80145 K range, the proposed values of T at the beginning of the simulation are easily achieved at the end of the equilibration period. On the contrary, when temperatures higher than 145 K are suggested, the system is not able to reach the proposed temperature, and at the end of the simulation, T is always lower than the proposed value. In particular, at 1 atm of pressure, the system cannot reach the temperature of 155 K but only values of about 134 K. At the same time, a suddenly expansion of the system with a noticeable decrease of the simulated density (assuming values of the same order of magnitude as Ar gas, d = 1.92 103 g cm3)64 is observed. Other interesting solvent features come from the dependence on T of the configuration (potential) energy, Ecfg. For T e 145 K, MD simulations show that Ecfg increases with T (see Figure 3). However, the dependence of Ecfg on T looks different in the 75 e T e120 and 120 e T e 145 K temperature ranges. As can be seen in Figure 3, from 75 to 120 K, Ecfg varies linearly with T but with a different slope than from 120 to 145 K. Moreover, important changes in the values of Ecfg are observed by increasing T. For instance, at temperatures of about 145 K, Ecfg ≈ 23 eV is obtained, but further attempts to increase T cause a decrease of the temperature (T ≈ 140 K), accompanied by a noticeable change of Ecfg (Ecfg ≈ 0.3 eV). From this, variation of Ecfg is possible to estimate the order of magnitude of the enthalpy of vaporization, ΔHvap ≈ 6 kJ/mol, which is in reasonable agreement with the experimental value of 6.4 kJ/mol.64 6397
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Figure 4. The Ar density as a function of the temperature, T, at p = 48.98 bar. The red point is not calculated (see the text).
Figure 5. RDFs for (top) ArAr and (bottom) BzAr.
The previous results of density, diffusion coefficient, and Ecfg seem to point out a phase transition within the range of investigated temperatures. In order to further analyze the behavior of the solvent, several MD simulations at a pressure of 48.98 bar, which is the critical point pressure,64 have been performed. Therefore, additional MD simulations at increasing values of T have been carried out starting from T = 150 K. It has been observed that the density decreases by increasing T, as shown in Figure 4. In fact, an increase of about 10 brings about a small gradual density decrease, but when the T changes from 161 to 162 K, a pronounced decrease of density, passing from 0.85 to 0.29 g cm3, is observed. Taking into account that the critical temperature and pressure have been experimentally found at 151 K and 48.98 bar with an Ar density equal to 0.5377 g cm3, results of Figure 4 and particularly the red point, could indicate T and p conditions and density values close to those of the critical point. The previous results reflect that the adopted potential model is adequate to describe important characteristics of the solvent.
Figure 6. Isosurface (1.81, 1.3, and 1.0 isovalues in the top, middle, and lower panels of the figure, respectively) plots of the Ar probability density in the solvated BzNaþ aggregate, derived from MD results obtained at T = 275 K.
Solvation of Benzene. Due to the high number of Ar atoms considered, as indicated above, the new and simpler atomatom formulation of BzAr interaction has been adopted. Therefore, the first step of the present study has been focused on the comparison of atomatom and atombond34,54 predictions for the BzAr dimer, concerning both the energy and the geometry at equilibrium. Only small differences of about 2 and 1% on the potential energy and distances at equilibrium, respectively, have been found. This result indicates that the decomposition of the molecular polarizability of Bz in effective atomic polarizabilities can be exploited to describe the basic features of the BzAr 6398
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Figure 8. BzAr, ArAr, and ArNaþ RDFs in the Ar-solvated BzNaþ aggregate.
Figure 7. Isosurface (5.18, 3.25, and 1.81 isovalues in the top, middle, and lower panels of the figure, respectively) plots of the Ar probability density in the solvated BzCl aggregate, derived from MD results obtained at T = 275 K.
interaction. However, it must be indicated that such decomposition has been adopted in light of the previous and more accurate atombond partition of the Bz polarizability33 and using data from the literature.65,66 As has been indicated above, MD simulations have been performed by solvating Bz by an ensemble of 485 thermalized Ar atoms. The first simulation of solvated Bz was performed at T = 100 K. Then, following the same procedure as that used in the previous section, MD simulations at increasing values of T were carried out. The radial distribution functions (RDFs) for ArAr and Arcenter of mass (c.m.) of Bz, predicted at two values of T, are shown in the top and lower panels of Figure 5, respectively. As can be expected, RDFs become less structured with increasing T.
Moreover, Figure 5 shows more defined RDF peaks for ArAr than those for BzAr. The shoulder observed in the first peak (see the lower panel of the figure) can be explained by taking into account the RDF characteristics observed for some small BzArn clusters. For instance, the RDF of BzAr16,40 obtained at conditions close to the equilibrium, show three peaks that can be associated with three different groups of Ar atoms. The first peak at ∼3.6 Å has been assigned to the two Ar atoms lying in opposite positions along the C6 symmetry axis of Bz. The second, peaking close to 4.7 Å, corresponds to the Ar atoms placed above and below the benzene plane. Finally, the third peak at ∼5.4 Å has been associated with Ar atoms occupying in-plane positions. The first composed structure, observed in the lower panel of Figure 5, can then be associated with those Ar atoms defining the first solvation shell (in-plane and out-of-plane positions), while the second peak corresponds to a second solvation shell. Therefore, it is interesting to note that the atomatom formulation of the BzAr interactions used here reproduces the main characteristics of the RDF functions obtained by means of the ionbond formalism.40 Solvation of BzNaþ and BzCl. For the reasons stressed above, the accurate description of BzNaþ and BzCl requires the use of the ionbond formalism. This approach was already exploited in our previous studies of the solvation of BzNaþ and BzCl by Ar atoms performed by considering only a few 6399
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Table 5. Configuration Energy Including Additional Induction Effects Caused by the Simultaneous Presence of a Cation and Anion in the Proximity of Ar Atoms, Ecfgind, and Configuration Energy without Including the Additional Induction Effects (see text), Ecfg, for Some Small BzNaþClArn Aggregates n 1
2
4
6
8
Ecfgind/mev
6599
6745
6868
7099
8245
Ecfg/mev
6752
7025
7413
7625
7925
Figure 9. BzAr, ArAr, and ArCl RDFs in the Ar-solvated BzCl aggregate.
number of solvents (e30),40,41,43,44,67,68 obtaining information about some microsolvation effects. In fact, the study of small BzNaþArn and BzClArn clusters showed different propensities of the Ar atoms solvating BzNaþ and BzCl dimers.40,41,67,68 In the present investigation, the number of solvent atoms has been increased up to 485, and the spatial disposition of several Ar atoms around the Bzion dimers has been investigated. Because of the different equilibrium structure of the BzNaþ aggregate, with the cation placed along the C6 symmetry axis with respect to the BzCl aggregate and the anion occupying in-plane positions,35,37,38 the Ar solvation shell around BzCl tends to occupy a higher volume than that around BzNaþ. For instance, at T = 75 K, the maximum distance from Bz c.m. to Ar atoms is 5 Å larger for BzCl than that for BzNaþ. Accordingly, the Ar atoms are distributed in a more compacted way around BzNaþ than around BzCl, as shown in Figures 6 and 7, where the three-dimensional (3D) probability densities of Ar atoms at 275 K are represented for the solvated BzNaþ and BzCl, respectively. The 3D probability densities are obtained by introducing our simulation data (properly shifted to the inertial reference frame (see ref 67 for details)) into the Visual Molecular Dynamics (VMD) program. The Volmap tool included in the mentioned package has allowed us to visualize isosurfaces of constant Ar probability density.69 Some RDFs of Ar-solvated BzNaþ and BzCl systems are shown in Figures 8 and 9, respectively. In particular, the BzAr, ArAr, and Arion RDFs are represented in the top, middle, and lower panels of the figures, respectively. In agreement with
Figure 10. BzNaþClAr4 equilibrium configuration calculated with (lower panel) and without (top panel) the inclusion of V00ind (see the text).
the 3D density probabilities shown in Figures 6 and 7, larger distances are exhibited by the solvated BzCl (less compacted solvation shell). Moreover, as expected when comparing Figures 8 and 9, while ArAr RDFs (middle panels) appear to be very similar, the largest differences are observed for the BzAr RDFs (top panels of the figures). Despite this, independently of the solvated ionic dimer (BzNaþ or BzCl), the position of the first peak for the BzAr RDFs is approximately the same, 3.7 Å for solvated BzNaþ and 3.6 Å for solvated BzCl. These values are very close to the corresponding equilibrium distance of the BzAr dimer.54 Moreover, the different heights of the peaks can be understood, bearing in mind that while Naþ and Ar compete for positions close to the C6 symmetry axis of Bz, Cl and Ar do not. As a matter of fact, the first peak corresponds to 6400
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Figure 12. (top) ArNaþ, (middle) ArAr, and (lower) ArCl RDFs of the BzNaþCl aggregate solvated by Ar atoms.
Figure 11. Kinetic energy, Ekin (top panel), and configurational energy, Ecfg (lower panel), represented as a function of total energy, Etotal, and temperature, T.
those Ar atoms of the first solvation shell placed below and above the aromatic ring and close to the C6 rotational axis of Bz. Therefore, as indicated before, Ar atoms compete for preferred positions with Naþ and not with Cl. The Arion RDFs represented in the lower panels of Figures 8 and 9 reflect again the more compacted solvation shell of BzNaþ in comparison with that of BzCl. On one hand, the larger distances observed in BzCl in comparison with those in BzNaþ and, on another hand, the second peak observed in the upper panels of Figures 8 and 9, placed at 4.4 and 4.7 Å, respectively, suggest more compacted solvation shells for BzNaþ than those for BzCl. The analysis of all RDFs shown in Figures 8 and 9 allows identification of the same peaks observed in the previous studies of smaller BzNaþArn40 and BzClArn67 clusters. Despite the position of the RDF peaks, reflecting a preferred distance of the Ar atoms from the ions and from the c.m. of Bz, which seems to indicate that the main features of the microsolvation processes remain unchanged when the number of Ar atoms increase, other approaches such as, for instance, global minimum optimization methods32 could help to elucidate information about such an issue. Solvation of BzNaþCl. In order to investigate the role played by the additional induction effects (see eq 8), some small BzNaþClArn aggregates have been also investigated. The
simplest BzNaþClAr aggregate has been considered, and the effects due to the simoultaneous presence of Ar, Naþ, and Cl have been analyzed. It has been found that, independently of the inclusion of the additional induction contribution, the aggregate mantains the same equilibrium configuration. Ar, Naþ, and Cl are placed on the same side of the aromatic plane, showing only small differences in both the ArNaþ (rArNaþ) and ArCl (rArCl) interatomic distances. When the induction energy, V00ind (represented in eq 8) is considered, a shortening of about 2% and an enlongation of about 4% are observed for rArCl and rArNaþ, respectively. These small changes in the equilibrium geometry bring about variations in the potential energy, with a V00ind contribution at equilibrium of about 150 meV. The remaining energy contributions to Etotal are not modified when additional induction effects are included. By increasing the number of Ar atoms, it has been observed that the mean contribution of any Ar atom to V00ind diminishes. In Table 5, the configuration energy values for BzNaþClArn (n = 1, 2, 4, 6, 8) with and without including V00ind, represented by Ecfgind and Ecfg, respectively, are given. On another hand, important effects on the equilibrium structures due to the role of V00ind have been observed. In fact, without the inclusion of the additional induction effects, the equilibrium structure of the small aggregates is more symmetric. The loss of symmetry, resulting from the inclusion of V00ind, can be apreciated in Figure 10, where the equilibrium structures of BzNaþClAr4 are represented. As can be seen, in both cases, the cation is placed along the C6 rotational axis of the aromatic compound. Moreover, when V00ind is not included, one of the Ar atoms is also placed close to the rotational axis of Bz, while 6401
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
Figure 14. Isosurface (0.5 isovalue) of the Ar probability density in the solvated BzNaþCl aggregate, obtained without including induction effects and derived from MD simulations pertormed at T = 275 K. The orange and red symbols represent Naþ and Cl, respectively.
Figure 13. Isosurface (1.03, 0.90, and 0.74 isovalues in the top, middle and lower panels of the figure, respectively) plots of the Ar probability density in the solvated BzNaþCl aggregate, derived from MD simulations pertormed at T = 275 K. The orange and red symbols represent Naþ and Cl, respectively.
the remaining atoms and the anion form a square around the cation (see the top panel of Figure 10). Conversely, when
including V00ind due to the attractive contribution for insertions of Ar atoms between the ClNaþ pair, the Ar atoms tend to point to the center of the ionic pair (see the lower panel of Figure 10). The inclusion of V00ind leads then to a destabilization of the small aggregates, which tend to diminish by increasing the number of Ar atoms. Following the same procedure used for simpler systems, MD simulations of the Ar-solvated BzNaþCl aggregate have been performed at increasing values of T. The mean values of Ecfg and the kinetic energy (Ekin) are represented in Figure 11 as a function of T (and of Etotal). As can be observed in Figure 11, increases of Etotal because of the stability of the BzNaþCl trimer when both ions are on the same side of the aromatic ring promote mainly an enhancement of Ekin. It has been also observed that the interaction energy of ClNaþ, BzNaþ, and BzCl remains almost unaffected by increasing T. Thus, despite the low-energy contribution of BzCl, due to the position out of the plane of the anion, the ionic trimer BzNaþCl is highly stabilized. In fact, the cation placed very close to the rotational axis of Bz interacts strongly with both Bz and Cl. MD simulations do not evidence changes in the cation position when T increases, while the anion rotates around the C6 axis of Bz even at very low values of T. This dynamics is not surprising considering that six equivalent equilibrium structures of the BzNaþCl aggregate, connected by small barriers, exist. Accordingly, the change in Ecfg can be attributed to the variation in the interactions of Ar atoms, bring about lessstructured RDFs when T increases (see Figure 12). An analysis of the Ar 3D density probabilities (Figure 13) shows that in the presence of both ions, Ar atoms tend to occupy preferably positions close to Cl. This behavior, probably promoted by additional induction effects, emerges also in the corresponding solvation shells, which do not show any symmetry around Bz, and looks like an incomplete solvation shell of the BzNaþCl aggregate. The solvation shells close to NaþCl are very compact, and to better visualize their characteristics in Figure 13, the Ar isoprobability density plots are represented in a transparent fashion. Three different isovalues are plotted in the figure, whose results have been derived from MD simulations performed at T = 275 K. In order to better analyze the influence of solventinduction effects, MD simulations at 275 K have been also performed without the inclusion of V00ind (see eq 8). While no noticeable differences in the energy have been observed, suggesting that induction does not provide the main contribution to V, it has been found that V00ind determines the high asymmetry observed in the 3D density probability of Ar atoms, as can be seen by 6402
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A
ARTICLE
comparing the results shown in Figure 13 with those given in Figure 14, where the solvation shell has been derived from MD simulations performed without the inclusion of such a three-body component.
’ AUTHOR INFORMATION
4. CONCLUSIONS The semiempirical methodology, providing the noncovalent intermolecular potential V as a sum of nonelectrostatic and electrostatic contributions, has been generalized to describe more complex systems. In this paper, the bigger system analyzed is that formed by two ions (Naþ and Cl) interacting with Bz, solvated by Ar atoms. The idea that V can be formulated as a combination of a few leading components has been extended to incorporate many-body effects, mainly originating from the presence of more than one ion in the cluster. More specifically, following some considerations concerning the decomposition of molecular polarizability into effective atomic or bond tensor components, V has been defined as sum of two-, three-, and fourbody components, each one having a physical meaning. In particular, the four-body and part of the three-body components have been ascribed to nonadditive induction effects. It has been emphasized that the proposed formulation ensures the correct representation of all dissociation channels, where the potential energy depends on the specific interaction of the involved partners. Therefore, such formulation appears to be suitable to carry out MD simulations. A study of some characteristics of the pure solvent has been performed, considering the involved PES given as a sum of twobody interactions, each one described by means of the ILJ function, depending only on the interatomic distance. MD results, obtained at increasing values of T, seem to indicate a transition from the liquid to gas phase. Despite the limitations arising from the use of an ensemble of only 500 Ar atoms, the potential model is able to describe some macroscopic characteristics of Ar in the condensed phase. RDFs of solvated Bz allow one to distinguish between those Ar atoms placed around the symmetry axis above and below the aromatic ring and those occupying preferably in-plane positions. Moreover, the solvation of both BzCl and BzNaþ shows the same basic characteristics observed when analyzing the microsolvation of smaller BzClArn and BzNaþArn aggregates. However, it has also been observed that the solvation shell of BzNaþ is more compact than that of BzCl. Finally, the investigation of the Ar solvation of the BzNaþCl aggregate has evidenced a noticeable asymmetry of the solvation shell. The observed tendency of the solvated Ar BzCl system, in which the Ar atoms tend to occupy positions close to Cl, is enhanced in the case of BzNaþCl. This propensity has been attributed to three-body induction effects appearing because of the simultaneous presence of Naþ, Cl, and Ar. The extension of the potential model to include nonadditive induction effects caused by the simultaneous presence of a cation and anion in the proximity of Ar atoms has been performed to assess the reliability of the PES to describe the unsolvated systems by comparing the predicted results for BzNaþ,35,36 BzCl,37,38 and BzNaþCl46 by means of accurate ab initio calculations. Further progress of our study will be achieved by substituting Ar atoms by water molecules, whose polarizability value (1.47 Å3) is very close to that of Ar (1.64 Å3). Results of investigations on small systems, N-methylacetamidewater70 and benzenewater,71 appear to be very promising.
’ ACKNOWLEDGMENT M. Albertí acknowledges financial support from the Ministerio de Educacion y Ciencia (Spain, Project CTQ2010-16709, and Mobility Program, Project PR2010-0243). Also, thanks are due to the Centre de Supercomputacio de Catalunya CESCA-C4 and Fundacio 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.
Corresponding Author
*E-mail:
[email protected].
’ REFERENCES (1) Alonso, J. L.; Antolínez, S.; Bianco., S; Lesarri, A.; Lopez, J. C.; Caminati, W. J. Am. Chem. Soc. 2004, 126, 3244. (2) M€uller-Dethelfs, K.; Hobza, P. Chem. Rev. 2000, 100, 143. (3) Aquilanti, V.; Cornicchi, E.; Moix Teixidor, M.; Saendig, N.; Pirani, F.; Cappelletti, D. Angew. Chem., Int. Ed. 2005, 44, 2336. (4) Dykstra, C. E.; Lisy, J. M. J. Mol. Struct.: THEOCHEM 2000, 500, 375. (5) Ondrechen, M. J.; Berkovitch-Yellin, Z.; Jortner, J. J. Am. Chem. Soc. 1981, 103, 6586. (6) Fried, L. E.; Mukamel, S. J. Chem. Phys. 1991, 96, 116. (7) Schmidt, M.; Le Calve, J.; Mons, M. J. Chem. Phys. 1993, 98, 6102 and references therein. (8) Mons, M.; Courty, A.; Schmidt, M.; Le Calve, J.; Piuzzi, F.; Dimicoli, I. J. Chem. Phys. 1997, 106, 1676. (9) Easter, D. C.; Bailey, L.; Mellot, J.; Tirres, M.; Weiss, T. J. Chem. Phys. 1998, 108, 6135. (10) Schmidt, M.; Mons, M; Le Calve, J. Chem. Phys. Lett. 1991, 177, 371. (11) Brupbacher, Th.; Makarewicz, J.; Bauderm, A. J. Chem. Phys. 1994, 101, 9736. (12) Hobza, P.; Bludsky, O.; Selzle, H. L.; Schlag, E. W. Chem. Phys. Lett. 1996, 250, 402. (13) Lenzer, T.; Luther, K. J. Chem. Phys. 1996, 105, 10944. (14) Bernshtein, V.; Oref, I. J. Chem. Phys. 2000, 112, 686. (15) Vacek, J.; Konvicka, K.; Hobza, P. Chem. Phys. Lett. 1994, 220, 85. (16) Vacek, J.; Hobza, P. J. Phys. Chem. 1994, 98, 11034. (17) Dullweber, A.; Hodges, M. P.; Wales, D. J. J. Chem. Phys. 1998, 106, 1530. (18) Riganelli, A.; Memelli, M.; Lagana, A. Lect. Notes Comput. Sci. 2002, 2331, 926. (19) Zoppi, A.; Becucci, M.; Pietraperzia, G.; Castellucci, E.; Riganelli, A.; Albertí, M.; Memelli, M.; Lagana, A Clustering properties of rare gas atoms on aromatic molecules, 16th International Symposium on Plasma Chemistry, Taormina, Italy, June 2003; p 22. (20) Koch, H.; Fernandez, B.; Makariewicz, J. J. Chem. Phys. 1999, 111, 198. (21) Kumpf, R. A.; Dougherty, D. A. Science 1993, 271, 163. (22) Dougherty, D. A. Science 1996, 261, 1708. (23) Ma, J. C.; Dougherty, D. A. Chem. Rev. 1997, 97, 1303. (24) Cabarcos, O. M.; Weinheimer, C. J.; Lisy, J. M. J. Chem. Phys. 1998, 108, 5151. (25) Cabarcos, O. M.; Weinheimer, C. J.; Lisy, J. M. J. Chem. Phys. 1999, 110, 8429. (26) Morais-Cabral, J. H.; Zhou, Y.; MacKinnon, R. Nature 2001, 414, 37. (27) Lehn, J. M. Supramolecular Chemistry; VCH: Weinheim, Germany, 1995. 6403
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
The Journal of Physical Chemistry A (28) Meyer, E. A.; Castellano, R. K.; Diederich, F. Angw. Chem., Int. Ed. 2003, 42, 1210. (29) Lopinski, G. P.; Wayner, D. D. M.; Wolkow, R. A. Nature 2000, 406, 48. (30) MacGillivray, L. R.; Atwood, J. L. J. Chem. Soc., Chem. Commun. 1997, 477. (31) Marques, J. M. C.; Pereira, F. B.; Leit~ao, T. J. Phys. Chem. A 2008, 112, 6079. (32) Llanio-Trujillo, J. L.; Marques, J. M. C.; Pereira, F. B. J. Phys. Chem. A 2011, 115, 2130. (33) Albertí, M.; Castro, A.; Lagana, A.; Moix, M.; Pirani, F.; Cappelletti, D.; Liuti, G. J. Phys. Chem. A 2005, 109, 2906. (34) Pirani, F.; Albertí, M.; Castro, A.; Moix, M.; Cappelletti, D. Chem. Phys. Lett. 2004, 394, 37. (35) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Cappelletti, D.; Coletti, C.; Re, N. J. Phys. Chem. A 2006, 110, 9002. (36) Coletti, C.; Re, N. J. Phys. Chem. A 2006, 110, 6563. (37) Coletti, C.; Re, N. J. Phys. Chem. A 2009, 113, 1578. (38) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Coletti, C.; Re, N. J. Phys. Chem. A 2009, 113, 14606. (39) Pirani, P.; Brizi, S.; Roncaratti, L. F.; Casavecchia, P.; Cappelletti, D.; Vecchiocattivi, F. Phys. Chem. Chem. Phys. 2008, 10, 5489. (40) Albertí, M.; Aguilar, A.; Lucas, J. M.; Lagana, A.; Pirani, F. J. Phys. Chem. A 2007, 111, 1780. (41) Huarte-Larraga~na, F.; Aguilar, A.; Lucas, J. M.; Albertí, M. J. Phys. Chem. A 2007, 111, 8072. (42) Albertí, M.; Aguilar, A.; Lucas, J. M.; Cappelletti, D.; Lagana, A.; Pirani, F. Chem. Phys. 2006, 328, 221. (43) Albertí, M.; Castro, A.; Lagana, A.; Moix., M; Pirani, F.; Cappelletti, D. Eur. Pys. J. D 2006, 38, 185. (44) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. Theor. Chem. Acc. 2009, 123, 21. (45) Albertí, M.; Pacifici, L.; Lagana, A.; Aguilar, A. Chem. Phys. 2006, 327, 105. (46) Albertí, M.; Aguilar, A.; Pirani, F. J. Phys. Chem. A 2009, 113, 14741. (47) Kim, D.; Lee, E. C. X.; Kim, K. S.; Tarakeshwar, P. J. Phys. Chem. A 2007, 111, 7980. (48) Izatt, R. M.; Pawlak, K.; Bradshaw, J. S.; Bruening, R. L. Chem. Rev. 1991, 91, 1721. (49) Izatt, R. M.; Pawlak, K.; Bradshaw, J. S.; Bruening, R. L. Chem. Rev. 1995, 95, 2529. (50) Islam, M. S.; Pethrick, R. A.; Pugh, D.; Wilsom, M. J. J. Chem. Soc., Faraday Trans. 1998, 94, 39. (51) Nathanson, G. M.; Davidovits, P.; Worshop, D. R.; Kolb, C. E. J. Phys. Chem. 1996, 100, 13007. (52) Laskin, A.; Gaspar, D. J.; Wang, W.; Hunt, S. W.; Cowin, J. P.; Colson, S. D.; Finlayson-Pitts, B. J. Sience 2003, 301, 340. (53) Roncero, O.; Villarreal, P.; Delgado-Barrio, G.; Gonzalez Platas, J.; Breton, J. J. Chem. Phys. 1998, 109, 9288. (54) Albertí, M.; Castro, A.; Lagana, A.; Pirani, F.; Porrini, M.; Cappelletti, D. Chem. Phys. Lett. 2004, 392, 514. (55) Pirani, F.; Cappelletti, D.; Liuti, G. Chem. Phys. Lett. 2001, 350, 286. (56) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. J. Phys. Chem. A 2010, 114, 11964. (57) Cambi, R.; Cappelletti, D.; Liuti, G.; Pirani, F. J. Chem. Phys. 1991, 95, 1852. (58) Cappelletti, D.; Liuti, G.; Pirani, F. Chem. Phys. Lett. 1991, 183, 297. (59) Aquilanti, V.; Cappelletti, D.; Pirani, F. Chem. Phys. 1996, 209, 299. (60) Pirani, F.; Maciel, G. S.; Cappelletti, D.; Aquilanti, V. Int. Rev. Phys. Chem. 2006, 25, 165. (61) Capitelli, M.; Cappelletti, D.; Colonna, G.; Gorse, C.; Laricchiuta, A.; Liuti, G.; Longo, S.; Pirani, F. Chem. Phys. 2007, 338, 62. (62) The DL_POLY Molecular Simulation Package. http://www.cse. clrc.ac.uk/ccg/software/DL_POLY/index.shtml.
ARTICLE
(63) Allen, M. P.; Tildesley, D. J.; Computer Simulation of Liquids; Clarendon Press: Oxford, 1989. (64) Material Safety Data Sheet Liquid Argon. http://www.uigi.com/ MSDS_liquid_Ar.html. (65) Ewig, C. S.; Waldman, M.; Maple, J. R. J. Phys. Chem. A 2002, 106, 326. (66) Gavezotti, A. J. Phys. Chem. B 2002, 107, 2344. (67) Huarte-Larraga~ na, F.; Aguilar, A.; Lucas, J. M.; Albertí, M. Theor. Chem. Acc. 2011, 128, 757. (68) Albertí, M.; Huarte-Larraga~ na, F.; Aguilar, A.; Lucas, J. M.; Pirani, F. Phys. Chem. Chem. Phys. 2011, 13, 8251. (69) Humphrey, W.; Dalke, A; Schulten, K. VMD — Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33. (70) Albertí, M.; Faginas Lago, N.; Lagana, A.; Pirani, F. Phys. Chem. Chem. Phys. 2011, 13, 8422. (71) Albertí, M.; Faginas Lago, N.; Pirani, F. 2011 submitted.
6404
dx.doi.org/10.1021/jp202995s |J. Phys. Chem. A 2011, 115, 6394–6404
