Alkali-Ion Microsolvation with Benzene Molecules - The Journal of

Article pubs.acs.org/JPCA

Alkali-Ion Microsolvation with Benzene Molecules J. M. C. Marques,*,† J. L. Llanio-Trujillo,† M. Albertí,‡ A. Aguilar,‡ and F. Pirani§ †

Departamento de Química, Universidade de Coimbra, 3004-535 Coimbra, Portugal IQTCUB, Departament de Química Física, Universitat de Barcelona, 08028-Barcelona, Spain § Dipartimento di Chimica, Università di Perugia, 06123-Perugia, Italy ‡

ABSTRACT: The target of this investigation is to characterize by a recently developed methodology, the main features of the first solvation shells of alkaline ions in nonpolar environments due to aromatic rings, which is of crucial relevance to understand the selectivity of several biochemical phenomena. We employ an evolutionary algorithm to obtain putative global minima of clusters formed with alkali-ions (M+) solvated with n benzene (Bz) molecules, i.e., M+−(Bz)n. The global intermolecular interaction has been decomposed in Bz−Bz and in M+−Bz contributions, using a potential model based on different decompositions of the molecular polarizability of benzene. Specifically, we have studied the microsolvation of Na+, K+, and Cs+ with benzene molecules. Microsolvation clusters up to n = 21 benzene molecules are involved in this work and the achieved global minimum structures are reported and discussed in detail. We observe that the number of benzene molecules allocated in the first solvation shell increases with the size of the cation, showing three molecules for Na+ and four for both K+ and Cs+. The structure of this solvation shell keeps approximately unchanged as more benzene molecules are added to the cluster, which is independent of the ion. Particularly stable structures, so-called “magic numbers”, arise for various nuclearities of the three alkali-ions. Strong “magic numbers” appear at n = 2, 3, and 4 for Na+, K+, and Cs+, respectively. In addition, another set of weaker “magic numbers” (three per alkali-ion) are reported for larger nuclearities.

1. INTRODUCTION In the last decades, the study of clusters has attracted the attention of many researchers from several areas, which is due to their peculiar properties in comparison with bulk matter. Specifically, clusters may be regarded as intermediate entities that connect the single atom or molecule with the corresponding bulk phase. The knowledge of the arrangement of atoms or molecules that minimize the interaction energy in the cluster is a key-step to understand chemical phenomena, such as microsolvation and crystal growth. Special attention has been paid to molecular aggregates involving aromatic rings for which noncovalent intermolecular interactions control basic phenomena typically governed by the combination of various components (see for instance ref 1). It has been recognized2−4 that interactions between aromatic rings contribute, among other, to the base-pair stacking in DNA, to protein structure or to the formation of supramolecular architectures. Particularly important structural motifs arising in the interaction between two benzene molecules are the parallel-displaced, T-shaped edge-to-face, and eclipsed face-to-face geometries. The former two were shown to be the most stable structures5 by employing the CCSD(T) ab initio method at the complete basis-set limit approach. The latest experimental and computational works concerning the understanding and application of noncovalent interactions involving aromatic rings have been recently reviewed by Diederich and collaborators,6 which update a previous review on the same topic by the same group.7 Also in the same review,6 the importance of cation−π interactions for binding biologically relevant molecules (e.g., protein−DNA © 2012 American Chemical Society

complexes) has been highlighted. In turn, previous theoretical work has established8−11 that the strength of the cation−π interaction arises from two major contributions, i.e., electrostatic attraction as well as induction. However, the critical balance of these and other components, as dispersion attraction and size/exchange repulsion, is crucial to assess the selectivity on controlling static and dynamic phenomena. The quantification, at the ab initio level, of the cation−π interactions when other interaction contributions are present, for instance in solvated systems, is a difficult task and the construction of suitable semiempirical potential energy surface (PES) models to describe the main features of ionic systems, properly reproducing all energetic details of the dissociative channels, is therefore an important issue in the dynamics study of these systems. In these cases, the PES model should ensure the accurate assessment of the strength and relative role played by the various components of the interaction. Considering that bound clusters are mainly dominated by pairwise additive twobody forces,12 the approximation most often used in the construction of potential energy functions (to be used later in dynamic simulations) is to describe the noncovalent interactions as a pairwise summation of atom−atom model potentials. However, the need for including contributions other than the pure atom−atom ones has been repeatedly discussed in the literature.13−16 Received: March 5, 2012 Revised: April 16, 2012 Published: April 19, 2012 4947

dx.doi.org/10.1021/jp302136u | J. Phys. Chem. A 2012, 116, 4947−4956

The Journal of Physical Chemistry A

Article

molecules. As it is known, in the most stable configuration of the M+−Bz aggregates19 the cation is placed along the C6 rotational axis of Bz. In principle, due to the charge distribution of Bz, cations tend to tie the Bz rings on a pile. However, π−π interactions should distort the pile. The investigation of the microsolvation effects of Bz molecules surrounding alkali cations is interesting to analyze the competition between cation−π and π−π interactions. The plan of the paper is as follows. In section 2, we present the potential model used for the title systems, while the evolutionary algorithm (EA) applied for the global optimization of the cluster structures is briefly described in section 3. The results of the global optimization are presented and discussed in section 4, where a comparison of the microsolvation clusters for the three alkali ions is given. Finally, the main conclusions are gathered in section 5.

In recent times, some of us have collaborated in the development of an alternative representation of the PES, based on the partition of the molecular polarizability, which is assumed the key property to scale both attraction and repulsion in molecular systems (see for instance refs 17 and 18). The reliability of the model to describe cation−π interactions has been proved by comparing the model predictions for the M+− Bz systems with high quality ab initio results.19 Moreover, it has been used in molecular dynamics (MD) simulations of M+−Bz aggregates solvated by Ar atoms (see for instance refs 20 and 21). In particular, these studies have analyzed the isomerization processes and the disposition of the Ar atoms around the M+− Bz aggregates. However, as the system becomes larger, it is difficult to establish the characteristics of the global minimum. Due to the existence of several minima having similar energies, in MD simulations the isomerizations are frequent, even at low temperatures, and the global minima must be characterized by means of global optimization studies. Although many methods have become available in the last decades22−35 for the global optimization of atomic and molecular clusters, this continues to be, in general, a very hard task that cannot be performed in a routine way. In the last five years or so, two of the present authors have been actively involved in the development of completely unbiased evolutionary algorithms (EAs) that seek the global minimum structure of atomic31,36−39 and molecular40 clusters; also, the performance of another unbiased optimization algorithm, the so-called big-bang method, has been recently analyzed in detail by one of these authors.41,42 Our EAs have been successfully applied to discover the global minima of several cluster systems, such as argon43 (up to Ar78), short-ranged Morse37 (up to 80 atoms), mixed argon− krypton38 and binary Lennard-Jones38,39 (up to 55 atoms), and cyclohexane solvated by Ar atoms.44 In the application to binary Lennard-Jones clusters, the EA was able to discover a new global minimum (not yet reported in the Cambridge Cluster Database45) for the cluster with 38 atoms whose size ratio between unlike atoms is 1.05.38 In turn, our recently proposed EA38 for the global optimization of molecular clusters has been shown to be reliable for discovering the minimumenergy structure of clusters of water (up to 20 molecules), benzene (up to 30 molecules), and cation benzene (up to 20 molecules). Microhydration of alkali ions is a relevant subject that has been addressed by employing global optimization methods. Hartke and collaborators46−49 have applied their own EA to conduct a thorough global optimization investigation on the structural trends of the Na+(H2O)n, K+(H2O)n, and Cs+(H2O)n microsolvation clusters. Different global minimum geometry structures were observed depending on the cation. In addition, Wales and co-workers50 have obtained the global minimum structures of Li+(H2O)n (n ≤ 20) and discovered that such clusters tend to have the ion located on the surface with a preferential coordination number of four water molecules in tetrahedral positions. They have also shown that, in similar clusters where Li+ is substituted by Ca2+, a different structural feature arises: Ca2+ occupies a position in the interior of the cluster and it is coordinated by eight water molecules in a square antiprismatic geometry. In turn, the influence of the ion size has been also observed for the Ar−M+ aggregates and for the Ar-solvated M+−Bz systems.21 In the present study, we are interested in employing our recently proposed EA40 on a properly formulated interaction potential to study the microsolvation features of alkali cations when solvated by benzene

2. POTENTIAL ENERGY SURFACE The potential model has been described before,1 and here we only outline its main characteristics and its extension to the present systems. The total potential energy, V, is decomposed in electrostatic, Vel, and nonelectrostatic, Vnel, contributions. To calculate Vel, we have considered the same Bz charge distribution used in our precedent studies (see for instance ref 1) which reproduce the quadrupole moment of Bz. This ensures an accurate description of both, the cation quadrupole and the quadrupole−quadrupole interactions at long-range. Vnel is described as a sum of pair-potential contributions (of two or three-body type) involving interaction centers, placed on different molecules, that are separated by a distance r. The pair contributions are all described by means of an improved Lennard-Jones (ILJ) function,51,52 which removes most of inadequacies in the asymptotic regions presented by the traditional and extensively used LJ model. In particular, the correct description of the long-range attraction, overestimated by LJ, appears to be important especially when several partners are involved in the formation of aggregates. 2.1. Ion−Benzene Interaction. The nonelectrostatic contribution (V nel) for the ion-benzene interaction is decomposed in a sum of twelve ion-bond terms (six M+−CC and six M+−CH). Such formulation exploits the representation of the molecular polarizability, affecting both dispersion and induction attractions, as sum of bond polarizability components53,54 6

Vnel =

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

+

i

i

(1)

i=1

Each term in eq 1 is described by means of the ILJ function, by assuming that the dispersion/induction centers are placed on the center of the bonds and on the cation. The ILJ function is expressed as

VILJ

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

2

2

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

r r0

4948

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

2

2

m

⎜

0⎟

(2)



Article

where ε and r0 are calculated from the corresponding parallel (ε∥ and r0,∥) and perpendicular (ε⊥ and r0,⊥) components,1 derived from two limiting approaches from the cation to the bond center and using the wheigthed sum, 2

2

ε = ε cos α + ε⊥ sin α

(3)

r0 = r0, cos 2 α + r0, ⊥ sin 2 α

(4)

parameters (see Table 2) are directly derived from the effective values of the atomic polarizabilities. According with the Bz Table 2. Potential Parameters Used to Define the Benzene− Benzene Nonelectrostatic Interaction Component

where α is the angle of approach. The first term in eq 2 describes the size repulsion while the second one is for the combination of dispersion and induction attraction. Due to the additional β parameter and to the dependence on r of both coefficients and the exponent in the first term, the ILJ function becomes more flexible with respect to the traditional LJ model. Moreover, the relevant well depth (ε) and equilibrium distance (r0) parameters, directly related to polarizability and charge of involved partners,1,51 assume a transferable character. The potential parameters are given in Table 1. Vel is calculated by applying the Coulomb law between

interaction

r0/Å

ε/meV

m

β

C−C C−H H−H

4.073 3.733 3.099

3.340 1.720 1.610

6 6 6

9.0 8.5 9.0

charge distribution given above, the electrostatic component of the Bz−Bz interacion is expressed as follows: 12

Vel =

12

∑∑ i=1 j=1 12

+

qCqC 4πε0r(Ch)i (Ch)j 6

∑∑ j=1 i=1

Table 1. Potential Parameters Used to Define the Ion− Benzene Nonelectrostatic Interaction Component interaction

r0,⊥/Å

r0,∥/Å

ε⊥/meV

ε∥/meV

m

β

Na+−CC Na+−CH K+−CC K+−CH Cs+−CC Cs+−CH

2.848 2.601 3.266 3.044 3.638 3.445

3.149 2.808 3.547 3.240 3.894 3.632

33.01 62.15 22.95 39.97 18.20 29.42

102.2 62.73 75.77 42.70 64.11 32.57

4 4 4 4 4 4

8.5 8.5 8.5 8.5 10.0 10.0

12

∑ i=1

6

qCqM+ 4πε0r(Ch)i M +

+

∑ i=1

(5)

where r(Ch)iM+ is the distance from M+ to any charge associated to the C atoms and rHiM+ is the distance from M+ to the H atoms. 2.2. Benzene−Benzene Interaction. As it was made previously to investigate large systems,55,56 the Bz−Bz interaction has been constructed by further decomposing the CC and CH bond polarizability in effective atomic components (whose values are different, often minor, from those of the corresponding isolated gas phase atoms), in order to account for the many body effects arising from the formation of chemical bonds. However, their sum must be consistent with the value of the molecular polarizability. Each atom with an assigned effective polarizability is, then, considered a center from which a partial contribution to the full intermolecular interaction can be defined. Accordingly, in this case Vnel is defined as follows: 6

Vnel =

6

∑ ∑ (VC − C i

i=1 j=1

j

∑∑ i=1 j=1 6

+

6

∑∑ i=1 j=1

qCqH 4πε0r(Ch)i Hj qHqH 4πε0rHiHj

(7)

3. EVOLUTIONARY ALGORITHM FOR GLOBAL OPTIMIZATION The EA proposed by two of us40 is employed in this work to search for the global minima of alkali-ion microsolvation clusters. The motivation for applying our EA is 2-fold. First, we want to test the performance of the EA for different kinds of molecular aggregates, and hence, those arising in the microsolvation of alkali-ions with benzene constitute a natural choice after the recent application to benchmark systems.40 Second, EAs are considered to be state-of-the-art algorithms for global optimization, and in particular, our algorithm has already shown to give good results for other molecular clusters.40 Moreover, the clusters studied in this work are difficult to optimize, especially when the number of benzene molecules surrounding the alkali-ion is greater than 15. Just to give an example, the EA localized 4 minima structures of K+−(Bz)16 whose difference in energy to the putative global minimum is less than 1 kJ mol−1. Nonetheless, the structural differences to the lowest-energy minimum are non-negligible. Indeed, by employing a recently developed methodology for obtaining the best superposition between pairs of structures (see ref 57 and references therein), we calculated the values 0.75, 1.14, 1.21, and 1.84 Å for the root-mean-square deviation (rmsd) arising in the comparison of each of the four local minima with the global minimum. Although the previous applications40 of this EA have used only site−site interaction potentials (while the present energy functions show a more complicated angular dependence1), it is straightforward to put on the same framework to discover the lowest-energy arrangement of benzene molecules around an alkali ion. Since the details of the EA were already presented in the original paper,40 we just describe here the main features of the method as applied to the title systems. The EA search begins with a random generation of a set of cluster structures (designated as the population) that are tentative solutions for the global optimization problem. Then, an iterative process is carried on, so that in each iteration (socalled generation) the whole population is replaced by a new set of solutions (designated as offspring) that arise from the

qHqM+ 4πε0rHiM+

4πε0r(Ch)j Hi

+

6

where r(Ch)i(Ch)j, r(Ch)iHj (or r(Ch)jHi) and rHiHj represent the distance between two charges of different molecules associated to C atoms, to C and H atoms, and to H atoms, respectively.

the charge of the cation (1 au) and the charge distribution on Bz, which consists in 12 negative charges qC = −0.04623 au placed on the C atoms, above and below the aromatic ring and separated by 1.905 Å and 6 positive charges of qH = +0.09246 au on each H atom. Accordingly, Vel is expressed as follows: Vel =

qCqH

12

+ VCi − Hj + VHi − Cj + VHi − Hj) (6)

and each term is described by means of the ILJ function (eq 2); i and j label the two interacting Bz molecules. The potential 4949



Article

Figure 1. Global minimum structures of the Na+−(Bz)n clusters. Plots were produced with VMD program.67

Figure 2. Global minimum structures of the K+−(Bz)n clusters. Plots were produced with VMD program.67

application of genetic operators. Specifically, the generation of offspring is achieved by performing the following steps: (i) Application of tournament selection, a probabilistic operator widely adopted by EAs, to choose the current best solutions from the population (i.e., the parents of the next population). (ii) Application of the genetic operators to the parent solutions selected in (i). First, crossover is applied to all pairs of parents. Afterward, mutation is applied to each gene (i.e., coordinate of the system) of the resulting solution. Both crossover and mutation operators are applied with prespecified probabilities. The present EA employs the simulated binary crossover (SBX)58 and sigma mutation36 that are state-of-the-art operators for numeric optimization. (iii) Local optimization and evaluation of the offspring. Before evaluating the quality of the offspring created in

(ii), each new solution is relaxed to the nearest local optimum by applying the quasi-Newton L-BFGS method.59,60 To avoid a deterioration of the best solution from one generation to the other, the EA applies an elitism operator whenever all offspring solutions have worse quality than the parent cluster structures. In practice, the elitism operator substitutes the worst-quality offspring by the best-quality parent. The iterative process ends when a predetermined number of function evaluations (or, equivalently, a given number of generations) is reached. Since the EA optimization is a stochastic, essentially unbiased, approach, it is important to run several times the generational process by departing from randomly generated populations in order to guarantee that the method is statistically meaningful. In general, we have employed 30 runs for each cluster size. Finally, the putative 4950



Article

Figure 3. Global minimum structures of the Cs+−(Bz)n clusters. Plots were produced with VMD program.67

distance on the Na+−Bz dimer (∼2.25 Å), indicating that Na+ is placed at the intersection point of the rotational axes of the two Bz molecules (thus forming a nonstacked structure). For Na+−(Bz)3, the third Bz molecule is placed at a slightly larger distance from Na+ (∼3.1 Å) and also the distance of the other two is slightly enlarged. By increasing the ion size, the interaction energy of the M+−Bz fragment decreases19 and, accordingly, the distance from M+ to the c.m. of Bz increases. This allows a reorientation of the Bz molecules which can form structures in which the cation is placed along the rotational axes of the three Bz molecules. As a matter of fact, while Na+ can not accommodate three Bz molecules at the same distances, the K+−(Bz)3 cluster forms an equilibrium structure having similar distances from K+ to the c.m. of the three aromatic rings [see Figure 4b]. As regarding Cs+−(Bz)3 in Figure 4c, it presents the same equilibrium distances from Cs+ to the c.m. of the three Bz molecules as for the Cs+−(Bz) and Cs+−(Bz)2 clusters. Moreover, the addition of a fourth Bz molecule forms a Cs+−(Bz)4 aggregate for which the four distances from Cs+ to the c.m. of any Bz are very similar. The Cs+−(Bz)4 structure at the equilibrium is very compact, and no fluctuations in the distances from Cs+ to the Bz c.m. are observed when the number of Bz molecules increases. This means that the core structure is built stepwise: for n = 1, the alkali ion is placed along the C6 rotation axis, above (or below) the benzene ring; for n = 2, the additional Bz molecule goes to the opposite side, in a nonstacked position (although it tends to a more stacked structure as the size of the ion diminishes); the addition of a third Bz molecule leads to a structure with a triangular shape (i.e., with Bz molecules occupying three of the faces of a tetrahedron where the ion is in the center). Although a similar pattern is observed for the three alkali ions, the Na+−(Bz)3 cluster is more asymmetric since it shows, as it has been indicated above, two similar dc.m. distances and a third one slightly bigger. This asymmetric core motif of the sodium-ion clusters essentially remains unchanged as n increases, but for both n = 9 and 13 global minima, the two coinciding dc.m. distances of the core structure split to accommodate a third Bz molecule in a closer position to the ion than for other nuclearities. The core motif for potassium and cesium ion

global minimum of the cluster is, then, assigned to the best solution obtained during the total number of runs employed. For the EA, we have employed the set of parameters previously used for benzene clusters:40 population size, 100; tournament size, 5; crossover rate, 0.75; mutation rate, 0.1; σ, 0.1. Conversely, we found that the 50000 evaluations per run used for benzene clusters can be reduced in the present work, since the best candidate structure to be the global minimum of microsolvated alkali-ion clusters is reached very early in each search. Such difference may be due to the strong ion-benzene interaction (not present in the pure benzene clusters) which basically “freezes” the orientational degrees of freedom of Bz molecules that are close to the alkali-cation (especially those in the first solvation shell); probably, this effect leads to a decrease in the number of the low-energy local minima of the M+-(Bz)n clusters in comparison with those for (Bz)n, which reduces the searching effort of the EA. Indeed, the global minimum is always obtained within 5000 (10000) evaluations for clusters with n ≤ 15 (n > 15).

4. RESULTS AND DISCUSSION In the equilibrium configuration of the smaller aggregates, the Bz molecules tend to be placed in such a way to allocate the cation along their C6 rotational axis. This can be seen clearly for the M+−(Bz)2 aggregates. However, when the number of Bz molecules increases, the competition between M+−π and π−π interactions originates a distortion from the preferred geometry of the M+−Bz fragment. Besides the observation of the structures in Figures 1−3, the elucidation of the growing pattern of the ion clusters benefits from the representation of the distances that separate the ion from the center of mass (c.m.) of each benzene molecule (hereafter designated as dc.m.); this is shown in Figure 4 for sodium-ion [panel a], potassiumion [panel b], and cesium-ion [panel c] clusters. We observe from this figure that the growing pattern of the ion clusters is, essentially, based on the addition of solvent to a core structure composed by a central ion surrounded by a small set of Bz molecules. In Figure 4a, it can be observed that, for Na+−(Bz)2 (at equilibrium), the distance from Na+ to the c.m. of the two Bz molecules is the same and coincides with the equilibrium 4951



Article

same happens for the analogous Na+ cluster ion.62 Then, it appears that the first solvation shell of this mixed cluster is formed by three (two) benzene and one (two) water molecules in the case of K+ (Na+). Such a difference in the first solvation shell was attributed62 to a smaller Na+ ionic radius (by about 30%) in comparison to the K+ one. This experimental achievement is compatible with the present theoretical result that shows, in the first solvation shell, one fewer benzene molecule for Na+ than for K+ (see above). Moreover, Lisy and collaborators mention in ref 61 an unpublished work carried out by themselves on K+−(Bz)n (n = 1−6) to suggest that the first solvation shell of these clusters comprises four benzene molecules, which is in perfect agreement with our result (cf. Figure 2). It is interesting to note in Figures 1−3 that the Bz microsolvation of Na+ tends to present more stacked Bz molecules than the corresponding ones for K+ and Cs+ ions. In general, the Na+−(Bz)n clusters are the most compact ones, although they show out layer distances for Na+−(Bz)15 and Na+−(Bz)18. Also, among the molecules placed outside the core structure, the Na+−(Bz)n clusters presents a more scattered pattern for the dc.m. values, which indicates a less definite (i.e., more diffuse) second solvation-shell than for K+ and Cs+. For Na+, the second solvation-shell extends from ∼4.5 up to ∼8.5 Å, whereas for both K+ and Cs+, the corresponding range is about 6−8 Å. This may be attributed to the strongest Na+−Bz interaction, which leads to a subtle balance between a more effective attraction of Bz molecules beyond the first solvationshell and the organization of the solvent to maximize the increasingly important Bz−Bz contribution as n increases. It is worth noting from Figure 5 that most of the global minimum structures show the alkali-ion displaced from the

Figure 4. Scatter plot of the distances between the center of the benzene rings and alkali-ion as a function of the cluster size: (a) sodium ion, (b) potassium ion, and (c) cesium ion. See text.

clusters is achieved by adding the fourth Bz molecule to the base face of the n = 3 tetrahedron structure. We notice that such structure presents approximately the same value for all dc.m. distances, that are slightly increased in comparison to the n = 3 cluster (especially in the case of K+). However, the dc.m. distances of the core Bz molecules split into different values for most of the potassium ion clusters with n > 5; an exception occurs only for n = 14, where very similar dc.m. values are obtained for the four molecules forming the core. Conversely, the cesium ion clusters show essentially the same dc.m. value for the core molecules, which is independent of n. Experimental insight about alkali-ions microsolvated with benzene can be extracted from the work of Lisy and collaborators61,62 on the gas-phase cluster ions of the form M+−(Bz)n(H2O)m, with M = Na or K. These authors have employed infrared and mass spectrometric techniques to show that the K+−Bz interaction is strong enough to lead to a partial dehydration of the ion when a third benzene molecule is added to K+−(Bz)2(H2O)2,61,62 while there is no evidence that the

Figure 5. Distance of the ion to the center of the cluster as a function of n. See text.

center of the cluster. Indeed, the ion position only coincides with the geometrical center of the cluster for a small set of nuclearities: n = 3 and 4 for K+, and n = 4 and 8 for Cs+; in the case of Na+, the position of the ion never coincides with the center of the cluster, though it becomes very close (i.e., dion−cluster < 0.5 Å) for n = 3 and 21. It is also apparent from Figure 5 that, though showing a significant oscillation with n, bigger ions tend to be more displaced from the center of the cluster, especially for n > 10. This may be related with the fact that the M+-Bz interaction diminishes as the size of the ion increases, which leads to a structure where the number of Bz− 4952



Article

in the right scale of the experimental determinations.64 As it can be seen in the figure, by increasing n, the M+−π contribution at first strongly increases (more negative values) but, from a certain value of n, its increasing rate significantly diminishes. Obviously, the π−π interaction is enhanced (more negative) when the number of Bz molecules increases. Accordingly, by increasing n, the π−π contribution is expected to become more important than the M+−π one. As a matter of fact, the M+−π contribution dominates the increase in the energy of the cluster with n only for small nuclearities (i.e., essentially for the formation of the first solvation shell), and the π−π interaction becomes the major contribution for n > 18 (n > 19) in the case of Cs+ (Na+ and K+). In Figure 7, we show the second energy difference, i.e.

Bz close-contacts tend to be maximized (i.e., in the case of K+ and Cs+ clusters). However, there is not any unambiguous trend for the comparison among the three alkali-ion clusters. A similar statement has been made by Beyer63 about the influence of monovalent ions on establishing interior vs peripheral structures of hydrated metal-ion clusters. The cluster energy contributions, arising from ion-benzene (M+−π) and benzene−benzene interactions (including CH− CH, CH−π, and π−π components) globally indicated, for simplicity, as π−π, are represented in Figure 6 as a function of the cluster size (n). Similar trends for both interactions have been observed for the clusters of the three alkali ions, and reported values for the M+−(Bz)n contributions for small n are

Δ2 E = En + 1 + En − 1 − 2En

(8)

Figure 7. Second energy difference for the microsolvated alkali-ion− benzene clusters. The calculated values are represented by solid dots (for sodium), asterisks (for potassium) and full squares (for cesium). For clearness, the results are shown in two panels and the points are linked by three different lines (see inserted key); the zero is indicated by the dotted line.

which is an approximation of the curvature of the energy, En, as a function of the cluster size; hence, such a quantity is related to the stability of each n cluster in comparison to its n − 1 and n + 1 neighbor clusters, and it is usually compared with experimental mass spectral intensities. Thus, maxima of Δ2E correspond to the so-called “magic number” clusters. It is apparent in Figure 7 that alkali-ion clusters present very stable structures for n = 2, 3, and 4 respectively for Na+, K+, and Cs+. These particularly stable structures are due to the ability to surround the alkali-ion with a specific number of solvent molecules at a distance close to the ion-benzene equilibrium geometry (which minimizes the strain energy of the cluster); such equilibrium values are 2.25, 2.65, and 3.06 Å for Na+−Bz, K+−Bz, and Cs+−Bz interactions, respectively. Indeed, we observe in Figure 4 that the maximum number of benzene molecules close to the ion, which does not significantly distort the corresponding equilibrium geometry, is 2, 3, and 4 respectively for Na+, K+, and Cs+. Clearly, this pattern follows the increasing ion size. In addition, other magic numbers appear: n = 11, 14, and 18 for Na+ clusters; n = 7, 13, and 17 for K+ clusters; n = 10, 13, and 19 for Cs+ clusters. It is apparent from Figure 4 that the n = 13 (for K+ and Cs+) or the n = 14 (for Na+) magic numbers are associated with structures that are more compact than the corresponding neighbor-size clusters. Such effect is especially relevant for Na+ and Cs+ that occupy more central positions in the cluster than in the corresponding neighbor-size clusters (cf. Figure 5). In contrast, the n = 13 magic-number structure for K+ is more (less) compact than the

Figure 6. Contributions for the global minima energy of alkali-ion clusters: (a) Na+−(Bz)n, (b) K+−(Bz)n, and (c) Cs+−(Bz)n. The energy contribution from the ion−benzene interactions is represented by the dotted line, whereas energy from the benzene−benzene pair potentials is shown with a dash line. Also shown by the solid line is the energy of the global minima. See text. 4953



Article

corresponding n = 14 (n = 12) structure; also in this case, we observe in Figure 5 that the K+ ion is less (more) displaced from the center of the cluster than the n = 14 (n = 12) one. In agreement with these observations, we show in Figure 8 that

where the ion tends to occupy a place on the surface of the cluster;47,66 however, K+ is an intermediate situation whose smaller microhydration clusters are similar to the Na+−(H2O)n structures, whereas bigger ones approach more closely the global minimum geometries of the Cs+−(H2O)n aggregates.46 A far-reaching orientational influence of Na+ on the benzene solvent is also apparent in Figure 1. This “structure-breaking” effect on solvent leads the sodium ion to be close to as many benzene molecules as possible, despite eventual less preferred Bz−Bz orientations, which generally favors a central position of Na+ in the cluster (cf. Figure 5). Conversely, more favorable Bz−Bz orientations have a larger contribution in the case of K+ and Cs+ and, hence, the global minimum structure tends to sacrifice an ion−Bz proximity (especially for those benzene molecules beyond the first solvation shell) in order to gain some energy due to Bz−Bz interactions. Because of this, the K+ and Cs+ ions are located off-center in the cluster structures (mainly for n > 10), as displayed in Figure 5. Another interesting outcome due to the above-mentioned balance between ion−solvent and solvent−solvent interactions is apparent for the first solvation shell of the alkali ions in both water and benzene. Indeed, the number of water molecules directly linked to the ion has shown46 to oscillate as n increases (and the shell closure of K+ and Cs+ only occurs at n = 16 and 18, respectively). In contrast, the M+−Bz clusters keep the number of benzene molecules in the first solvation shell as n increases. Since the size of benzene is much greater than water, the shell closure for M+−Bz clusters occurs for a smaller number of solvent molecules than in the corresponding M+− (H2O)n aggregates.

Figure 8. Ion−benzene (dash line) and benzene−benzene (solid line) contributions for the second energy difference represented in Figure 7: (a) sodium-ion clusters, (b) potassium-ion clusters, and (c) cesium-ion clusters. See text.

5. CONCLUSIONS We have performed a global geometry optimization of clusters resulting from the alkali-ion microsolvation with benzene molecules. Specifically, three alkali-ion systems (i.e., Na+, K+, and Cs+) solvated with benzene up to 21 molecules have been studied. The ion−benzene and benzene−benzene interactions were represented by physically meaningful model potentials that, in particular, ensure an accurate description of the longrange part. In such a model, the interactions are decomposed into electrostatic and nonelectrostatic contributions, with the latter being formulated in terms of a sum of ILJ functions. The equilibrium configuration of the smaller aggregates allows to see that the Bz molecules tend to be placed in such a way to allocate the cation close their C6 rotational axis. This can be seen clearly for the M+−(Bz)2 aggregates. When the number of Bz molecules increases, not all the Bz molecules can adopt the configuration of minimum energy of the M+−Bz fragment. This can be interpreted on the basis of the different interaction strength of the M+−Bz aggregates, which diminishes when the size of the cation raises. Indeed, as in the microsolvation of alkali-ions with water molecules,50,65 “structure-breaking” effect on benzene is larger for the smaller cations (e.g., Na+). It is apparent in the present work that the number of Bz molecules allocated in the first solvation shell increases with the size of the cation (three Bz molecules for Na+ and four for both K+ and Cs+). This is compatible with previous experimental work involving microsolvated clusters of Na+ and K+.61,62 In turn, various global minimum structures arise for the three alkali ions, which present special stability in relation to the corresponding neighbor size clusters (so-called “magic numbers”). Strong “magic numbers” of the M+−(Bz)n clusters appear at n = 2, 3, and 4 for M = Na, K, and Cs, respectively.

the appearance of such magic numbers for Na+ and Cs+ are due to the ion−Bz interaction, whereas the Bz−Bz interaction presents a negative contribution. In addition, the ion-Bz interaction is also responsible for the magic numbers at n = 2 for Na+−(Bz)n, n = 3 and 7 for K+−(Bz)n, and n = 4 for Cs+− (Bz)n; in the case of K+−(Bz)n clusters, the Bz−Bz interaction has also a small contribution for appearing the n = 3 magic number. Conversely, the remaining above-mentioned magic numbers (i.e., n = 11 and 18 for Na+, n = 17 for K+, and n = 10 and 19 for Cs+) are due to the Bz−Bz interaction, whereas ion−Bz interaction presents, now, a negative contribution. Finally, we will comment on similarities and differences between the present results and the achievements pointed out for alkali ions solvated with water molecules. While water molecules tend to form hydrogen-bond networks which compete with the most favorable ion−water orientation, the ion−benzene interaction is clearly much stronger than the benzene−benzene one. Such competition in sodium−water clusters leads to global minima where Na+ is placed in the center (with most of the solvent molecules oriented toward it) or, in contrast, to structures with most of the water molecules belonging to hydrogen-bond networks and the ion in an offcenter position;49,65 a similar pattern is also observed for Li+− (H2O)n clusters.50 Due to a lesser “structure-breaking” effect on solvent, the larger-size K+ and Cs+ ions can be surrounded by several water molecules that are able to orient in such a way which allows for an easy formation of hydrogen-bond networks 4954



Article

This set of “magic numbers” has been attributed to the strong ion−benzene interaction. Additional weaker “magic numbers” (three per alkali ion) are observed for larger values of n; among each set of weak “magic numbers”, two (one) are (is) due to the Bz−Bz (M+−Bz) interaction. Although these weak “magic numbers” may appear in the mass spectrum of M+−Bz aggregates, the present theoretical result has to be looked at with caution. In fact, the relative stability of these “magic numbers” does not exceed more than a few kJ/mol, whereas the clusters were studied at T = 0 K without any consideration of the corresponding zero-point energies. Due to the fluxional character of these microsolvation clusters, several low-energy local minima are expected to compete with the global minimum as the temperature is increased. In a future work, it would be interesting to apply molecular dynamics simulations to investigate the real importance of those weak “magic number” structures.

■

(11) Soteras, I.; Orozco, M.; Luque, F. J. Phys. Chem. Chem. Phys. 2008, 10, 2616. (12) Leutwyler, S.; Jortner, J. J. Phys. Chem. 1987, 91, 5558−5568. (13) Schmidt, M.; Calvé, J. L.; Mons, M. J. Chem. Phys. 1993, 98, 6102−6120. (14) Ondrechen, M. J.; Berkovitch-Yellin, Z.; Jortner, J. J. Am. Chem. Soc. 1981, 103, 6586−6592. (15) Mons, M.; Courty, A.; Schmidt, M.; Calvé, J. L.; Piuzzi, F.; Dimicoli, I. J. Chem. Phys. 1997, 106, 1676−1686. (16) Easter, D. C.; Bailey, L.; Mellot, J.; Tirres, M.; Weiss, T. J. Chem. Phys. 1998, 108, 6135−6143. (17) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F. Theor. Chem. Acc. 2009, 123, 21. (18) Albertí, M.; Faginas-Lago, N.; Laganà, A.; Pirani, F. Phys. Chem. Chem. Phys. 2011, 13, 8422. (19) Albertí, M.; Aguilar, A.; Lucas, J. M.; Pirani, F.; Cappelletti, D.; Coletti, C.; Re, N. J. Phys. Chem. A 2006, 110, 9002−9010. (20) Albertí, M.; Aguilar, A.; Lucas, J. M.; Cappelletti, D.; Laganà, A.; Pirani, F. Chem. Phys. 2006, 328, 221−228. (21) Huarte-Larrañaga, F.; Aguilar, A.; Lucas, J. M.; Albertí, M. J. Phys. Chem. A 2007, 111, 8072−8079. (22) Deaven, D. M.; Ho, K. M. Phys. Rev. Lett. 1995, 75, 288−291. (23) Hartke, B. J. Phys. Chem. 1993, 97, 9973−9976. (24) Gregurick, S. K.; Alexander, M. H.; Hartke, B. J. Chem. Phys. 1996, 104, 2684. (25) Niesse, J. A.; Mayne, H. R. J. Comput. Chem. 1997, 18, 1233− 1244. (26) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111. (27) Roberts, C.; Johnston, R. L.; Wilson, N. T. Theor. Chem. Acc. 2000, 104, 123. (28) Leary, R. H. J. Global Optim. 2000, 18, 367. (29) Locatelli, M.; Schoen, F. Comput. Optim. Appl. 2003, 26, 173. (30) Shao, X.; Cheng, L.; Cai, W. J. Comput. Chem. 2004, 25, 1693. (31) Pereira, F. B.; Marques, J. M. C.; Leitão, T. Tavares, J. In Proceedings of the 2006 IEEE Congress on Evolutionary Computation; Vancouver, CEC, 2006; Vols. 1−6, pp 2270−2277. (32) Takeuchi, H. J. Chem. Inf. Model. 2006, 46, 2066. (33) Grosso, A.; Locatelli, M.; Schoen, F. Math. Program. Ser. A 2007, 110, 373. (34) Rossi, G.; Ferrando, R. J. Phys.: Condens. Matter 2009, 21, 084208. (35) Dieterich, J. M.; Hartke, B. Mol. Phys. 2010, 108, 279−291. (36) Pereira, F. B.; Marques, J. M. C.; Leitão, T. Tavares, J. In Advances in Metaheuristics for Hard Optimization, Springer Natural Computing Series; Siarry, P., Michalewicz, Z., Eds.; Springer: Berlin, 2008; pp 223−250. (37) Pereira, F. B.; Marques, J. M. C. Evol. Intel. 2009, 2, 121−140. (38) Marques, J. M. C.; Pereira, F. B. Chem. Phys. Lett. 2010, 485, 211−216. (39) Pereira, F. B. Marques, J. M. C. In Proceedings of the 2010 IEEE Congress on Evolutionary Computation; Barcelona, CEC, 2010; pp 1−7. (40) Llanio-Trujillo, J. L.; Marques, J. M. C.; Pereira, F. B. J. Phys. Chem. A 2011, 115, 2130−2138. (41) Marques, J. M. C.; Pais, A. A. C. C.; Abreu, P. E. J. Comput. Chem. 2010, 31, 1495. (42) Marques, J. M. C.; Pais, A. A. C. C.; Abreu, P. E. J. Comput. Chem. 2012, 33, 442−452. (43) Marques, J. M. C.; Pereira, F. B.; Leitão, T. J. Phys. Chem. A 2008, 112, 6079. (44) Abreu, P. E.; Marques, J. M. C.; Pereira, F. B. Comput. Theor. Chem. 2011, 975, 83−91. (45) The Cambridge Cluster Database, http://www-wales.ch.cam.ac. uk/ccd.html, accessed in March, 2012. (46) Schulz, F.; Hartke, B. ChemPhysChem 2002, 3, 98. (47) Hartke, B. Angew. Chem., Int. Ed. 2002, 41, 1468. (48) Schulz, F.; Hartke, B. Phys. Chem. Chem. Phys. 2003, 5, 5021. (49) Schulz, F.; Hartke, B. Theor. Chem. Acc. 2005, 114, 357. (50) González, B. S.; Hernández-Rojas, J.; Wales, D. J. Chem. Phys. Lett. 2005, 412, 23−28.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by Fundaçaõ para a Ciência e Tecnologia (FCT), Portugal, under Project PTDC/QUI/ 69422/2006, which is financed by Programa Operacional Factores de Competitividade (COMPETE) of QREN and FEDER programs. We are grateful for the provision of supercomputer time on the Milipeia cluster which is hosted at Laboratório de Computaçaõ Avançada, Universidade de Coimbra. M.A. and A.A. acknowledge financial support from the Ministerio de Educación y Ciencia (Spain, Project CTQ2010-61709/BQU) and to the Comissionat per a Universitats i Recerca del DIUE (Generalitat de Catalunya, Project 2009-SGR 17). Also thanks are due to the Centre de Serveis Cientı ́fics i Acadèmics de Catalunya CESCA 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.

■

REFERENCES

(1) Albertí, M.; Castro, A.; Laganà, A.; Moix, M.; Pirani, F.; Cappelletti, D.; Liuti, G. J. Phys. Chem. A 2005, 109, 2906−2911. (2) Fasan, R.; Dias, R. L. A.; Moehle, K.; Zerbe, O.; Obrecht, D.; Mittl, P. R. E.; Grütter, M. G.; Robinson, J. A. ChemBioChem 2006, 7, 515−526. (3) Malinovskii, V. L.; Samain, F.; Häner, R. Angew. Chem., Int. Ed. 2007, 46, 4464−4467. (4) Hänni, K. D.; Leigh, D. A. Chem. Soc. Rev. 2010, 39, 1240−1251. (5) Lee, E. C.; Kim, D.; Jurecka, P.; Tarakeshwar, P.; Hobza, P.; Kim, K. S. J. Phys. Chem. A 2007, 111, 3446−3457. (6) Salonen, L. M.; Ellermann, M.; Diederich, F. Angew. Chem., Int. Ed. 2011, 50, 4808−4842. (7) Meyer, E. A.; Castellano, R. K.; Diederich, F. Angew. Chem., Int. Ed. 2003, 42, 1210−1250. (8) Dougherty, D. A. Science 1996, 271, 163. (9) Tsuzuki, S.; Yoshida, M.; Uchimaru, T.; Mikami, M. J. Phys. Chem. A 2001, 105, 769. (10) Tsuzuki, S.; Mikami, M.; Yamada, S. J. Am. Chem. Soc. 2007, 129, 8656. 4955



Article

(51) Faginas-Lago, N.; Huarte-Larrañaga, F.; Albertí, M. Eur. Phys. J. D 2009, 55, 75−85. (52) Pirani, F.; Brizi, S.; Roncaratti, L. F.; Casavecchia, P.; Cappelletti, D.; Vecchiocattivi, F. Phys. Chem. Chem. Phys. 2008, 10, 5489. (53) Denbigh, K. G. Trans. Faraday Soc. 1940, 36, 936−948. (54) Smith, R. P.; Mortensen, E. M. J. Chem. Phys. 1960, 32, 502− 507. (55) Albertí, M.; Pirani, F. J. Phys. Chem. A 2011, 115, 6394−6404. (56) Albertí, M.; Faginas-Lago, N.; Pirani, F. J. Phys. Chem. A 2011, 115, 10871−10879. (57) Marques, J. M. C.; Llanio-Trujillo, J. L.; Abreu, P. E.; Pereira, F. B. J. Chem. Inf. Model. 2010, 50, 2129−2140. (58) Deb, K.; Agrawal, R. B. Complex Systems 1995, 9, 115−148. (59) Nocedal, J. Math. Comp. 1980, 35, 773−782. (60) Liu, D.; Nocedal, J. Math. Program. B 1989, 45, 503−528. (61) Cabarcos, O. M.; Weinheimer, C. J.; Lisy, J. M. J. Chem. Phys. 1998, 108, 5151−5154. (62) Cabarcos, O. M.; Weinheimer, C. J.; Lisy, J. M. J. Chem. Phys. 1999, 110, 8429−8435. (63) Beyer, M. K. Mass. Spectrom. Rev. 2007, 26, 517. (64) Ma, J. C.; Dougherty, D. A. Chem. Rev. 1997, 97, 1303−1324. (65) Hartke, B.; Charvat, A.; Reich, M.; Abel, B. J. Chem. Phys. 2002, 116, 3588. (66) Hartke, B.; Schulz, F. ChemPhysChem 2002, 3, 98−106. (67) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38.

4956


Alkali-Ion Microsolvation with Benzene Molecules - The Journal of

Recommend Documents