Effect of Anions on Static Orientational Correlations ... - ACS Publications

Jan 19, 2008 - for all three systems at temperatures 90 K above the corresponding melting ...... all three systems go to a plateau of β ) 1 after 1.5...
0 downloads 0 Views 361KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 1743-1751

1743

Effect of Anions on Static Orientational Correlations, Hydrogen Bonds, and Dynamics in Ionic Liquids: A Simulational Study Baofu Qiao,†,‡ Christian Krekeler,‡ Robert Berger,† Luigi Delle Site,‡ and Christian Holm*,†,‡ Frankfurt Institute for AdVanced Studies (FIAS), Johann Wolfgang Goethe-UniVersita¨t, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany, and Max Planck Institute for Polymer Research, Ackermannweg 10, D-55128 Mainz, Germany ReceiVed: July 26, 2007; In Final Form: NoVember 14, 2007

Three different ionic liquids are investigated via atomistic molecular dynamics simulations using the force field of Lopes and Pa´dua (J. Phys. Chem. B 2006, 110, 19586). In particular, the 1-ethyl-3-methylimidazolium cation EMIM+ is studied in the presence of three different anions, namely, chloride Cl-, tetrafluoroborate B F4 , and bis((trifluoromethyl)sulfonyly)imide TF2N . In the focus of the present study are the static distributions of anions and cations around a cation as a function of anion size. It is found that the preferred positions of the anions change from being close to the imidazolium hydrogens to being above and below the imidazolium rings. Lifetimes of hydrogen bonds are calculated and found to be of the same order of magnitude as those of pure liquid water and of some small primary alcohols. Three kinds of short-range cation-cation orderings are studied, among which the offset stacking dominates in all of the investigated ionic liquids. The offset stacking becomes weaker from [EMIM][Cl] to [EMIM][BF4] to [EMIM][TF2N]. Further investigation of the dynamical behavior reveals that cations in [EMIM][TF2N] have a slower tumbling motion compared with those in [EMIM][Cl] and [EMIM][BF4] and that pure diffusive behavior can be observed after 1.5 ns for all three systems at temperatures 90 K above the corresponding melting temperatures.

I. Introduction Room-temperature ionic liquids (ILs) refer to the class of salts whose melting temperatures are below 100 °C.1-5 In recent decades, ILs have received increasing interest because of their attractive properties: they are liquid in a wide temperature range, are nonvolatile, are nonflammable, have high thermal stability and ionic conductivity, have relatively low viscosity, and so forth. ILs have been advertised as “green solvents” in chemical, catalytical, and electrochemical reaction media. The most common ILs are composed of nitrogen-based heterocyclic cations, for example, imidazole and pyridine, in combination with a variety of anions. Changing the cations or the anions allows a practical tuning of the physicochemical properties to meet specific requirements; therefore, ILs are also referred to as “designer solvents”. Moreover, a wide range of solutes display good solubility in ILs because of their polarity. ILs have already been used as alternative solvents in organic chemical synthesis, electrodepositions, battery electrolytes, and so forth. Many experiments have been performed to study the structure,6 the transport property,7-9 the conductivity,8-11 the viscosity,8,9 the dielectric response,12-15 the nanoscale segregation,16 and some other properties. Relationships between several of these properties have been proposed on the basis of computer simulation17 and theory.18 Recently, theories have been introduced to explain the low melting temperatures of ILs,19 to define ILs by their transport properties,20 and to study the interface of ILs and nonILs.21 In order to be able to simulate the family of ILs, several force fields were developed,22-28 among which that of Lopes et al. has already been used to study several different ILs.14,29-31 * Corresponding author. E-mail: [email protected]. † Johann Wolfgang Goethe-Universita ¨ t. ‡ Max Planck Institute for Polymer Research.

Even though many properties, for example, density, melting temperature, radial distribution function (RDF), mean-square displacement (MSD), and electrical conductivity, have been studied frequently, some other properties have been rarely investigated. Among these are, for example, orientation correlations of the cations, strength of the hydrogen bonds (HB), the dielectric response, and the red-edge effect.32 In the present work, we will focus on the orientational correlation of the cations and anions, the lifetimes of hydrogen bonds and rotational and translational dynamics. HB interactions are believed to be important to affect the microscopic structure of ILs. The crystal structure analysis and IR study supported HBs in [EMIM[Cl][AlCl3].33 Liu et al. reported the C-H‚‚‚F HBs in several systems by ab initio calculations at the Hartree-Fock (HF) level using a small polarized split-valence double-ζ basis set augmented by diffuse basis functions (6-31+G(d)).27 Dong et al. reported that, in single-atomic anion ILs, only single-hydrogen-bonded ion-pairs exist, while there are single- and double-hydrogen-bonded ionpairs in multi-atomic anion ILs at the hybrid density functional theory level (B3LYP) with the 6-31+G(d) basis set.34 HB networks have been reported in [EMIM][Cl]34 and [EMIM][TF2N]35 from density functional theory (DFT) calculations (B3LYP/6-31+G(d)). From ab initio molecular dynamics simulations with the gradient corrected PBE density functional, Del Po´polo et al. proposed that HB interaction in [MMIM][Cl] can be classified as moderate HB according to the geometric criteria of HB,36 while the lifetimes of HBs in ILs have to the best of our knowledge never been reported. Besides HB, the aromatic stacking interaction is also believed to play an important role in nature. Three different structures are proposed to result from the aromatic interaction, that is, edge-face, offset stacking, and face-face stacking (see Figure

10.1021/jp0759067 CCC: $40.75 © 2008 American Chemical Society Published on Web 01/19/2008

1744 J. Phys. Chem. B, Vol. 112, No. 6, 2008

Figure 1. Geometries of aromatic interactions: (a) edge-face orientation, (b) offset stacked orientation, (c) face-face stacked orientation.

Figure 2. Schematic representation of (a) 1-ethyl-3-methylimidazolium (EMIM+), (b) chloride (Cl-), (c) tetrafluorobarate (BF4 ), and (d) bis((trifluoromethyl)sulfonyl)imide (TF2N-).

1).37,38 The edge-face geometry (T-shaped, C-H‚‚‚π interaction) is commonly found between aromatic residues in proteins. The offset stacking (parallel displaced) is commonly found in proteins, and this is the geometry in the base stacking of DNA. The offset and the edge-face orientation are the most common arrangements and have similar energies. For ILs, however, orientations of the cations are rarely reported.6,26,39-43 X-ray diffraction studies give evidence for stacking in [C18MIM][Cl‚ H2O], [C18MIM][Cl‚H2O‚CHCl3], and [C14MIM][Cl‚H2O],39 while not in [BMIM][Cl].40 The structures determined by X-ray diffraction of [MMIM][TF2N] and 1,2,3-triethylimidazolium bis((trifluoromethyl)sulfonyl)imide showed that the cations are aligned in layers.41 On the other hand, a neutron diffraction study on the liquid [MMIM][TF2N] showed no evidence of stacking.6 By developing a united-atom model, Micaelo et al. reported the existence of offset stacking of the cations in both [BMIM][PF6] and [BMIM][NO3] and found that the offset stacking in [BMIM][PF6] is less pronounced than that in [BMIM][NO3].26 Liquid clathrates are reported in ionic liquid-aromatic mixtures.42,43 We would like to remark here that the stacked orientations discussed here are different from the nanoscale segregation of long alkyl chain ILs.16,44 In the present work, we apply the force field of Lopes et al., which has already been used by several groups to study different ILs. The focus of the present work is thus not the development of a new force field but instead to study systematically the influence of anions on the strength of hydrogen bonding and the stacked orientation of the cations within the framework of the force field of Lopes et al. For our systematic study, the anions, chloride (Cl-), tetrafluoroborate (BF4 ), and bis((trifluoromethyl)sulfonyly)imide ((CF3SO2)2N- or TF2N-), are used in combination with 1-ethyl-3-methylimidazolium (EMIM+) (see Figure 2). II. Simulation Details Classical molecular dynamics (MD) simulations were performed, using an all-atom model. The force field parameters

Qiao et al. were taken from a recent publication of Lopes and Pa´dua,28 which are based on the OPLS-AA/AMBER with several additional parameters obtained from ab initio calculations. The simulations were performed with the software package GROMACS 3.3.1.45,46 The leapfrog algorithm was used, with a time step of integration of 1 fs. The isobaric-isothermal ensemble (NPT) was utilized for all simulations. Temperature was scaled with the Berendsen thermostat (characteristic time τ ) 0.1 ps) in combination with the Berendsen barostat algorithm (characteristic time τ ) 1 ps, reference pressure 1 bar, compressibility 4.5 × 10-5 bar-1).47 Neighbor searching was done up to 1.2 nm and updated every five time steps. Lennard-Jones (LJ) interactions and forces were smoothly switched off between 0.9 and 1 nm.48 Electrostatic interactions were treated with a particle mesh Ewald (PME) summation49,50 with Fourier grid spacing of 0.12 nm, fourth-order interpolation, and a cutoff of 1.2 nm. The LJ 1-4 and Coulomb 1-4 interactions were scaled by 0.5. The long-range dispersion correction for the energy and pressure were applied. All covalent bonds involving hydrogens were constrained using the LINCS algorithm.51 Given that the melting temperature of [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N] is 360 K,52 288 K,8 and 257 K,8 respectively, the corresponding working temperature was chosen to be 450, 380, and 350 K, that is, about 90 K above the melting temperature. We found this necessary to get well-converged dynamical properties from MD runs53 that can be performed for several nanoseconds. Lower temperatures would have resulted in prohibitively long simulation times, as was checked in one example. All of the systems had 288 ion pairs. The equilibration was performed with a method similar to that of ref 53; that is, the initial configurations were equilibrated at higher temperatures for 2 ns and then cooled down stepwise and equilibrated again for 2 ns. For [EMIM][Cl], we use the temperatures, 600, 550, and 500 K; for [EMIM][BF4], we use 600, 550, 500, 450, and 400 K, and for [EMIM][TF2N], we use 600, 550, 500, 450, and 400 K. Afterward, the three systems were equilibrated for 2.5 ns at the corresponding working temperature, and then data were collected every 10 time steps for 5 ns. In the production stage, the potentials and the densities were monitored. The equilibrium density of the three systems was 1.07, 1.19, 1.53 g cm3, respectively. For the former two, the corresponding density at the working temperature, which is linearly extrapolated from the experimental data, is 1.04 g cm-3 for [EMIM][Cl]11 and 1.21 g cm-3 for [EMIM][BF4].8 For [EMIM][TF2N], the experimental value at 348 K is 1.47 g cm-3,54 and the simulation result at 343 K is 1.54 g cm-3.28 All of the deviations from experimental data are less than 4% and therefore within the expected accuracy. III. Results and Discussions A. Radial Distribution Function. The RDFs of the anions around the cations of the three systems are shown in Figure 3a. In the following calculations, the geometric centers of the imidazolium ring and the anion are used unless otherwise stated. The term imidazolium ring refers to the ring formed by the five atoms, that is, N1, C2, N3, C4, and C5 (see Figure 2a). The first peak of the three curves is located at (0.46 nm, 3.82), (0.50 nm, 2.99), and (0.54 nm, 2.03), respectively. The peak position is shifted to larger distance; the peak height decreases, and the curve becomes broader from [EMIM][Cl] to [EMIM][BF4] to [EMIM][TF2N]. On the RDF of [EMIM][Cl], there appears a shoulder at a distance around 0.7 nm. This shoulder becomes less pronounced in [EMIM][BF4] and disappears in [EMIM][TF2N]. The behavior of the RDFs, especially the

Effects of Anions

Figure 3. Radial distribution functions of (a) anion-cation, (b) cationcation and (c) anion-anion of [EMIM][Cl], [EMIMBF4], and [EMIM][TF2N]. The geometric centers of the anions and the imidazolium rings are used.

behavior of the first peaks and the shoulders, will be further analyzed in terms of spatial distribution function in the following section. The RDFs of the cations around the cations are shown in Figure 3b. The cation-cation RDF of [EMIM][Cl] exhibits a shoulder at the shorter distance side of the first peak (about 0.5 nm), while this shoulder becomes less apparent in [EMIM][BF4] and disappears in [EMIM][TF2N]. Andrade et al. showed that this shoulder exists in [EMIM][BF4] and [BMIM][BF4] but not in [EMIM][AlCl4] and [BMIM][AlCl4] by using the allatom force field developed by them.23 From the coarse-grained model, polarizable atomistic MD model and nonpolarizable atomistic MD model, Wang et al. showed that this feature has a temperature dependence: it is observed at 400 K but not at 700 K.25 The appearance of this shoulder can also be affected by the choice of force field and the way how the centers of the ions are defined. The shoulder on the RDF for [EMIM][NO3] was found by Wang et al.25 using an AMBER force field. On the contrary, in ref 53 using the AMBER force field and partial charges for EMIM+ from ref 55 and parameters for NO3- from ref 56, this shoulder was not observed. Similar to the above findings, we also found a large structural difference in the RDF that depends on the force field selection in [EMIM][BF4].57 From our point of view, this arises from the fact that when using different values for the partial charges, the interactions between anions and imidazolium hydrogens are different, which can therefore produce a noticeable effect on the distribution of anions and cation around cations. The imidazolium hydrogens refer to the three hydrogen atoms directly bonded to the imidazolium ring, that is, H2, H4 and H5 (see Figure 2a). The choice of the center, the geometric center of the imidazolium ring or the center of mass (COM) of the cation, affects whether this shoulder appears or not, for example, [BMIM][PF6].27 This arises from the fact that, because of the long side chain, butyl in [BMIM][PF6], the COM position of the cation can be substantially displaced from the location of the geometric center of the imidazolium ring. The further analysis in the following sections shows that this shoulder has specific physical meaning: it arises mainly from the offset stacking of the imidazolium rings. We will visualize this in terms of the spatial distribution functions in the following section. For consistency sake, the anion-anion RDFs are also presented in Figure 3c. From the anion-cation, cation-cation,

J. Phys. Chem. B, Vol. 112, No. 6, 2008 1745 and anion-anion RDFs, one can find that even at very large distances, all three systems show oscillations in the RDFs, which agrees with the idea that ILs are strongly ordered systems.53 B. Spatial Distribution Function. The liquid structure can be further investigated through the spatial distribution function (SDF), which describes the density distribution of a particle in three-dimensional space around a reference molecule. In the present work, the software package gOpenMol58 is used to visualize the SDFs. Figure 4a shows the SDF results of the three different anions around the cations (a.1 for Cl-, a.2 for BF4 , and a.3 for TF2N-). To study the first peaks on the anion-cation RDFs, the distance range 0-0.54 nm, 0-0.59 nm, and 0-0.54 nm is used for [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N], respectively. In Figure 4a, the SDFs are shown by the blue and red isosurfaces, which represent 10 and 5 times the corresponding average density. In [EMIM][Cl], the anions prefer to be located close to the imidazolium hydrogens, that is, H2, H4, and H5. In [EMIM][BF4] and [EMIM][TF2N], this trend becomes weaker and weaker, and the anions prefer positions above and below the imidazole ring. Computer simulations on [MMIM][F], [MMIM][Cl], and [MMIM][Br]52 and neutron diffraction experiments on [MMIM][Cl] and [MMIM][PF6]6 have also indicated that the regions above and below the imidazolium rings are preferred in case of large anions. The small anions feel the attraction of the local bond dipoles at the imidazolium ring. As the ion size increases, the anion is attracted by the delocalized total positive charge, and due to steric repulsion, the preferred position is above and below the imidazolium ring. The SDFs in the shoulder ranges, 0.54-0.66 nm and 0.590.74 nm, on the anion-cation RDF curves in [EMIM][Cl] and [EMIM][BF4] are also plotted in Figure 4a, shown by the green and yellow isosurfaces. The green and yellow isosurfaces represent 8 and 4 times the corresponding average density. Because we found no shoulder in the anion-cation RDF of [EMIM][TF2N], the corresponding SDF is not shown. It can be seen that the contributions for the shoulders come mainly from the Cl and F atoms that are close to the alkyl-, especially the methyl-hydrogen atoms. These distributions agree with the findings for [MMIM][Cl].29 To further study the cation-cation orientation in the short distance range (from 0 up to the shoulder left to the first peak on the cation-cation RDF curves; see caption of Figure 4), the SDF of the cations around the cations are calculated for [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N], and plotted in Figure 4b, shown by the blue and red isosurfaces. The blue and red isosurfaces represent 5 and 2.5 times the corresponding average density. The anion-cation SDFs are also plotted, denoted by black isosurfaces (6 times the corresponding average density), to display how the anions and cations coexist around the cations. The cation-cation SDFs of [EMIM][Cl] and [EMIM][BF4] display similar features: (a) In the higher density distribution region, the imidazolium rings are found in an offset stacked orientation. From [EMIM][Cl] to [EMIM][BF4], the strength of the offset stacking decreases. (b) In the lower one, the imidazolium rings show another type of orientation besides the offset stacking. In this orientation, the nearest neighbor cations prefer to fill the region bounded by the probable bonds of the central C4-H4 and C5-H5 atoms with the anions. However, in the cation-cation SDF of [EMIM][TF2N], no offset stacking of the imidazolium rings is found, which is consistent with the fact that there is no shoulder left to the first peak on the cationcation RDF curve of [EMIM][TF2N].

1746 J. Phys. Chem. B, Vol. 112, No. 6, 2008

Qiao et al.

Figure 4. Three-dimensional density distribution function (SDF) of (a) anions and (b) cations around the cations for [EMIM][Cl](a.1, b.1), [EMIM][BF4] (a.2, b.2), and [EMIM][TF2N] (a.3, b.3). The geometric centers of the anions and the imidazolium rings are used, and binwidth of the center-center distance r is 0.01 nm. The blue, red, green, yellow, and black isosurface is plotted at 10, 5, 8, 4, and 6 times the corresponding average density, while in (b.1, b.2, and b.3) the blue and red isosurface is plotted at 5 and 2.5 times the corresponding average density. (a.1) r ) 0-0.54 nm for blue and red isosurfaces, and 0.54-0.66 nm for green and yellow ones. (a.2) r ) 0-0.59 nm for blue and red ones, and 0.59-0.74 nm for green and yellow ones. (a.3) r ) 0-0.54 nm. (b.1) r ) 0-0.5 nm. The black isosurface is the anion-cation SDF in r ) 0-0.54 nm. (b.2) r ) 0-0.54 nm. The black isosurface is the anion-cation SDF in r ) 0-0.59 nm. (b.3) r ) 0-0.77 nm. The black isosurface is the anion-cation SDF in r ) 0-0.54 nm.

The anion-cation SDFs of smaller anions (Cl-, BF4 ) in the present work are consistent with the computer simulation results of some other ILs,27,29,36,52 that is, the anions prefer the positions close to imidazolium hydrogens. The authors of ref 53 calculated the cation-cation SDF in [EMIM][NO3] and found that the nearest neighbor cations are well-localized at about 0.6 nm above and below the imidazolium ring. However, part of the cation-cation orientations in our lower density distribution region was not displayed in ref 53. The SDFs from neutron diffraction6 are frequently mentioned. However, these experimental results, both anion-cation SDFs and cation-cation SDFs, are different from computer simulation results to some extent. To explain this discrepancy, Del Po´polo et al.36 have argued that this difference might be due to the limitation of the empirical potential structure refinement (EPSR) model in ref 6 when studying complex molecular liquids: the RDF contains not enough information about intermolecular distances, while these distances are used as precise input quantities by the EPSR model. The above-mentioned ion specific behavior, concerning the cation-cation and anion-cation orientation, should be discussed in terms of hydrogen bonding interactions. The structure of cation-cation is mainly determined by the anion-cation interactions which is to a large extent determined by the hydrogen bonding, differing from anion to anion. The following section will substantiate this connection in more details. C. Hydrogen Bonding. Figure 5 displays all of the sitesite RDFs between the hydrogen atoms and the negatively charged atoms in the anions. In [EMIM][TF2N], the oxygen atoms approach the imidazolium ring closest, which is consistent with the results in 1-N-octyl-3-methylpyridinium bis((trifluoromethyl)sulfonyl)imide.22 In the H-F and H-N RDFs of [EMIM][TF2N], even though there are shoulders at about 0.3

nm, the first apparent peaks exist at about 0.5 nm, which is much larger than the corresponding van der Waals distances.34,59,60 These H‚‚‚F and H‚‚‚N interactions are relatively weak and thus neglected in the following calculation. Here, we employ the geometric definition of HB that involves the donor-acceptor distance criterion and the hydrogen-donoracceptor angle criterion.61 For the donor-acceptor distance criterion, the donor-acceptor RDFs are also calculated, and the first minimum of the donor-acceptor RDFs are taken as the upper limit of the donor-acceptor distance criterion, as listed in Table 1. Zero is the lower limit of the donor-acceptor distance criterion. The hydrogen-donor-acceptor angle criterion is to select all distances that occur within an angle ranging from 0 to 30°. To our knowledge, lifetimes of HBs in ILs have never been simulated or experimentally measured. There are two methods available to calculate it: the continuous HB method and the interrupted HB method (first-order dynamics method).62,64 In the first method, a HB must exist continuously, while in the second method, a HB is allowed to break and re-form. The continuous HB method depends on the frequency of saving the coordinates.62 Therefore, we apply the latter method, which has already been implemented in GROMACS.62 In this method, the kinetics of HB breakage and re-formation can be obtained from a chemical dynamics analysis.63 k

Ah B k′

(1)

with k and k′ as the forward and backward rate constant for HB breakage and formation. A binary function h(t) is defined as 1 when a HB exists, and is 0 otherwise. c(t) is the correlation function of h(t), that is, the conditional probability that a HB

Effects of Anions

J. Phys. Chem. B, Vol. 112, No. 6, 2008 1747

Figure 5. Site-site radial distribution functions of (a) [EMIM][Cl], (b) [EMIM][BF4], and (c) [EMIM][TF2N].

TABLE 1: Lifetimes (Picoseconds) of Hydrogen Bondingsa C2-H2‚‚‚X C4-H4‚‚‚X C5-H5‚‚‚X C8-H8‚‚‚X C6-H6‚‚‚X C7-H7‚‚‚X

[EMIM][Cl]

[EMIM][BF4]

[EMIM][TF2N]

1.389 (0.047) [0.43] 1.710 (0.059) [0.42] 2.476 (0.106) [0.42] 0.466 (0.014) [0.46] 0.349 (0.006) [0.43] 0.304 (0.007) [0.45]

0.310 (0.006) [0.42] 0.350 (0.001) [0.40] 0.401 (0.004) [0.40] 0.161 (0.002) [0.45] 0.145 (0.001) [0.41] 0.124 (0.001) [0.42]

1.037 (0.024) [0.43] 1.000 (0.031) [0.43] 1.182 (0.059) [0.42] 0.406 (0.004) [0.46] 0.309 (0.007) [0.42] 0.286 (0.003) [0.43]

a To determine the hydrogen bondings, the acceptor-donor distance criterion is taken from 0 to the first minimum of the acceptor-donor RDF (the values listed in square brackets, in nm), the hydrogen-donor-acceptor angle criterion is taken from 0 to 30°. For [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N], X represents Cl, F, and O, respectively. Standard deviations are given in parentheses.

exists, given the HB exits at t ) 0. The rate of the relaxation to equilibrium is characterized by the reactive flux correlation function K(t)

K(t) ) -

dc(t) dt

(2)

and

K(t) ) kc(t) - k′n(t)

(3)

Figure 6. Trans conformation of the bis((trifluoromethyl)sulfonyl)imide.

where n(t) is the probability that the distance criterion of a HB holds and the angle criterion is broken, given the HB exits at t ) 0. The lifetime in this method is obtained from the forward rate constant

in BF4 is -0.49 e, and for oxygen in TF2N -0.53 e. (2) The ion pairs and HB networks. In ref 34, it is reported that, in [EMIM][Cl], there is only single-hydrogen-bonded ion pairs, while single- and double-hydrogen-bonded ion pairs are in [EMIM][BF4]. References 6 and 65 have shown that TF2Nadopts a trans conformation (see Figure 6) predominantly in liquid phase, in comparison with a cis conformation in solid phase.41 The oxygens in the trans conformation of TF2N- can form maximally double-hydrogen-bonded ion pairs with the hydrogens in one cation because of steric repulsion. While considering the HB networks,34,35 in multi-atomic anion ILs, it is easier to form a HB network, which can promote the stability of the HBs. Therefore, it can be concluded that the strengths of HBs mainly depend on the partial charges and the formation of the HB networks. From the values presented in Table 1, we

τ)

1 k

(4)

The calculated lifetimes are listed in Table 1. From the dielectric relaxation, Daguenet et al. estimated that the lifetime of the ion pairs in ILs should be shorter than 100 ps.12 The results of our work actually agree with their estimate. Table 1 shows that the strengths of the HBs follow the order: [EMIM][Cl] > [EMIM][TF2N] > [EMIM][BF4]. There are two reasons related to this result. (1) The value of the partial charge. The charge at chloride is -1 e, while the partial charge of fluorine

1748 J. Phys. Chem. B, Vol. 112, No. 6, 2008

Figure 7. Geometric representation of a stacked orientation between the imidazolium rings. The x axis is obtained from the geometric center of the imidazolium ring to the H2 atom. The y axis is obtained from the geometric center of the imidazolium ring to the ethyl (not shown here). The z axis is the normal vector of the imidazolium ring, obtained from the right-hand rule of x axis and y axis. r represents the distance between the geometric centers of the imidazolium rings. d represents the offset distance from the geometric center of the 2nd imidazolium ring to the normal vector of the 1st imidazolium ring. R is the angle between the two normal vectors. β is the angle between the normal vector of the reference imidazolium ring and the vector of the geometric centers of the two imidazolium rings.

can also infer that the HBs involving the three imidazolium hydrogens have longer lifetimes because of the higher positive charges of these atoms. Let us also remark here that the lifetimes of the HBs in ILs are of the same order of magnitude as those of some other pure liquids, like H2O (1.20 ps), MeOH (2.05 ps), EtOH (2.71 ps), and PrOH (2.56 ps) at room temperature.62 From the viewpoint of the strength of HBs, this comparison supports the idea that HBs in ILs are moderate HBs.36 D. Stacked Orientation of the Imidazolium Rings. In aromatic systems, the packing structure is determined by C-H‚‚‚π or π‚‚‚π interactions, while the orientation of the imidazolium rings in IL systems to each other is mainly based on HB interactions. Since the imidazolium rings carry positive charges, stacking structures can only be explained by taking the hydrogen bonding into account, which we have discussed before. In the following, we will focus on the structures in our ILs that occur because of hydrogen bonding. In order to investigate the stacked orientation of the imidazolium rings, we first define some parameters which will facilitate the analysis. As shown in Figure 7, r is the distance between the geometric centers of the two imidazolium rings, d is the offset distance from the geometric center of the second imidazolium ring to the normal vector of the first imidazolium ring, and R is the angle formed by the two normal vectors. Table 2 shows the measured stacking parameters in the short distance ranges (from 0 up to the shoulders of the cation-cation RDFs, see caption of Table 2), and the observed broad distributions indicate the weak nature of the orientation. Figure 8 displays the isolines of the probability of finding an orientation with certain r and cos(R), where we consider a distance range r ) 0-1 nm. In the short distance range, we observe that the parallel and antiparallel orientations are dominant in smaller anion ILs, which confirms the offset stacking of the imidazolium rings shown by the blue isosurface in Figure 4b. This is a consequence of purely goemetric arguments, since in the short distance regime steric interactions prevent other than stacked arrangements. In the middle distance range, the vertical orientation dominates to some extent in [EMIM][Cl] and [EMIM][BF4]. Micaelo et al. reported that the offset stacking in [BMIM][PF6] is less pronounced than that in [BMIM][NO3].26 From our results and those of Micaelo et al., it can thus be concluded

Qiao et al. that, with the increase of the size of the anions, the positions above and below the imidazolium rings become preferentially occupied by anions, as we have explained before, and this in turn decreases the possibility for the offset stacking of the imidazolium rings. For [MMIM][TF2N], single-crystal X-ray diffraction studies support the stacking of the cations in the crystal phase,41 while neutron diffraction studies find no evidence of the stacking in the liquid phase.6 These experimental results and the temperature dependence of the shoulder on cation-cation RDF (it exists at 400 K, but not at 700 K25) support that the short-range stacking of the cation is reminiscent of the stacking observed in the solid phase. We have calculated the isolines of the probability distribution to determine which orientation is dominant in our three systems. In order to substantiate our previous conclusion, we use the following orientational distribution function to study the average orientation of the imidazolium rings

gR(r) )

1 N(r)

∑ i 90°, 180 - R is calculated. c See ref 26. d See ref 38.

Figure 10. Rotational dynamics (eq 6) of [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N].

strengths, and the anion-cation and cation-cation orientations. In this and the following sections, we will focus on the dynamics. The rotational dynamics of the imidazolium ring can be described by the following function22 Figure 8. Isolines of the distribution probabilities of imidazolium rings in respect to r (nm) and cos(R) of (a) [EMIM][Cl], (b) [EMIM][BF4], and (c) [EMIM][TF2N] in the distance range 0-1 nm. r and R have been shown in Figure 7.

C(t) )

1

N

∑

N i)1

3cos2Ri(t) - 1 2

(6)

where Ri(t) is the angle between the normal vectors of the ith imidazolium ring at time t and time t ) 0; N is the total number of imidazolium rings. When the imidazolium rings loose their orientations, the function C(t) will decay to zero. Figure 10 shows the rotational dynamics of the three systems investigated. The rotational time constant is calculated by22

τ)

Figure 9. Orientational distribution functions (eq 5) of (a) [EMIM][Cl], (b) [EMIM][BF4], and (c) [EMIM][TF2N]. The curves of [EMIM][Cl] are shifted by 2. The curves of [EMIM][BF4] are shifted by 1. The cone region is defined as β e 60°, as shown in Figure 7.

with increasing size of the anions, the offset stacking becomes weaker and weaker. E. Rotational Dynamics of the Imidazolium Rings. In the previous sections, we studied static properties of ILs, the HB

∫0∞ C(t)dt

(7)

The result obtained is 70, 54, and 150 ps for [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N], respectively. The fact that [EMIM][BF4] rotates fastest is consistent with our previous observation that its HB strength is the weakest. The slow decay in [EMIM][TF2N] can be reasoned to occur because in this IL the cations are surrounded by the bulky anions, TF2N-, as shown in Figure 4a.3. For [EMIM][TF2N], Kelkar and Maginn reported a rotational time constant of the longest end-to-end vector of the cation as 293 ps at 343 K.22 By comparing their result to ours, we can conclude that the normal vector of the imidazolium ring relaxes faster than the longest end-to-end vector of the cation. This follows the fact that the cation can easily rotate around the axis formed by the methyl and the ethyl groups, which are relatively immobile because of their big size compared to the sizes of the imidazolium hydrogens.

1750 J. Phys. Chem. B, Vol. 112, No. 6, 2008

Qiao et al.

Figure 11. (a) The type of the motion (eq 8) and (b) non-Gaussian parameter (eq 9) of [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N]. The geometric centers of the imidazolium rings are used.

F. Diffusive Behavior of the Imidazolium Rings. In the final section, we investigate the diffusive behavior of the imidazolium rings. To this end, we calculate their mean-square displacement ∆r2. The mean-square displacement can be divided into three typical regions, described as22,32,53 2 〈|b(t) r - b(0)| r 〉 ) ∆r2 ∝ tβ

(8)

where the value of β indicates the type of the diffusion. In the initial subpicosecond ballistic region, β is equal to 2, while in the long-time tail when the system shows true diffusive behavior, we expect β ) 1. The true diffusive behavior is very important to accurately determine the transport properties. In the intermediate region, the ions can be trapped in local cages, and therefore β < 1. Another variant to determine deviations from Gaussian behavior is to look at the so-called non-Gaussian parameter. It is defined as22,32,53

γ(t) )

3〈|∆r(t)|4〉 5〈|∆r(t)|2〉2

-1

(9)

with the meaning that, when the system exhibits true diffusive motion, the random walk dynamics is recovered, and γ approaches zero. Figure 11 displays our calculated β and γ of the three systems investigated. It is found that the ballistic motions last for about 0.1 ps for all three investigated systems. Figure 11a shows that all three systems go to a plateau of β ) 1 after 1.5 ns. However, from Figure 11b, we infer that neither of the three systems reaches the limiting value γ ) 0. This slight discrepancy arises from the higher-order term in calculating the non-Gaussian parameter, which is more sensitive to the mean-square displacement. In ref 22, Kelkar and Maginn have argued that the [EMIM][TF2N] system will not show Gaussian motion in 1 ns at 343 K, which is consistent with our results. Figure 11b shows that the peak position of [EMIM][Cl], [EMIM][BF4], and [EMIM][TF2N] is located at 39, 58, and 100 ps, respectively, where the diffusive behavior of the cations has the maximum deviation from Gaussian behavior.32,53 This shows that the translational behavior can be directly linked to the size of the anions. IV. Summaries and Conclusions In the framework of the all-atom force field of Lopes et al., we have calculated the anion-cation and cation-cation orientations in three pure ILs, where we have always taken the same

EMIM+ cation but used three different counterions, Cl-, BF4, and TF2N-, respectively, in order to study the size effects of the anions occurring in ILs. With the increase of the size of the anion, from Cl- to BF4 to TF2N , the preferred positions of the anions change from close to the imidazolium hydrogens to above and below the imidazole rings, while the cations move from above and below the imidazolium rings to larger distances. In a second step, we determined the lifetimes of hydrogen bondings. We find that the calculated lifetimes are consistent with the estimated value from the dielectric relaxation method. Moreover, the strengths of the hydrogen bondings follow the order: [EMIM][Cl] > [EMIM][TF2N] > [EMIM][BF4], which can be explained by the distribution of partial charges and the formation of HB networks in multi-atomic anion ILs. As expected, the interactions of the anions and the imidazolium hydrogens are stronger than those with the side chain hydrogens. The lifetimes of the hydrogen bondings of the three ILs are in the same order of magnitude as those of pure water and of some small primary alcohols, which demonstrates that the HBs of these IL are moderate HBs. We find evidence that there exists offset stacking of the imidazolium rings in all the ILs. And the strengths of the offset stackings follow the order: [EMIM][Cl] > [EMIM][BF4] > [EMIM][TF2N]. Besides the offset stacked orientation, two other orientations, the T-shaped edge-face and the vertical orientation in the region of C4-H‚‚‚anion and C5-H‚‚‚anion, coexist in our three investigated systems. However, our calculations show that the latter two orientations are only very weakly populated. Our dynamic studies of the rotational and translational motions of the imidazolium rings show that the imidazolium rings in [EMIM][BF4] display the fastest rotational diffusion. This is in agreement with the calculated strength of the hydrogen bondings. In contrast, the imidazolium rings in [EMIM][TF2N] relax slowest. This, we interpret, is due to the fact that they are surrounded by bulky anions. Second, all our three investigated systems reach Gaussian diffusive behavior after 1.5 ns at the temperatures 90 K above the corresponding melting temperatures. In ref 8, the authors calculated the ratio of the molar conductivity from complex impedance measurements to that obtained from NMR diffusion coefficient measurements and found that the ionic association of [EMIM][TF2N] is stronger than that of [EMIM][BF4]. The diffusion coefficient of the cations and the viscosity of [EMIM][TF2N] and [EMIM][BF4] were also reported in ref 8: for the former, 2.4 × 10-10 m2 s-1 and 8.2 mPa s at 350 K; for the latter, 3.3 × 10-10 m2 s-1 and 5.2 mPa s at 380 K. In our study, we found that the strength of the hydrogen bondings follow the order: [EMIM][TF2N] > [EMIM][BF4] at the corresponding temperature. Therefore, the above experimental data can be rationalized and explained through the strengths of the hydrogen bondings. We emphasize again that the conclusions drawn within this work are based on computed results obtained within a given force field and are thus only valid to the extend that the employed force field is capable to describe the relevant underlying physics. This implicit assumption is currently being investigated by us in more detail. Acknowledgment. We gratefully acknowledge financial support through a stipend to QB by the MPG and through the DFG Grants Ho 1108/15-1, Be 3791/1-1, DE 1140/2-1. Computer time was obtained through the Frankfurt Center for Scientific Computing (CSC).

Effects of Anions References and Notes (1) Huddleston, J. G.; Visser, A. E.; Reichert, W. M.; Willauer, H. D.; Broker, G. A.; Rogers, R. D. Green Chem. 2001, 3, 156. (2) Wasserscheid, P.; Keim, W. Angew. Chem. Int. Ed. 2000, 39, 3772. (3) Hajiwara, R.; Ito, Y. J. Flu. Chem. 2000, 105, 221. (4) Oliver-Bourbigou, H.; Magna, L. J. Mol. Catal. A Chem. 2002, 182-183, 419. (5) Wikes, J. S. J. Mol. Catal. A Chem. 2004, 214, 11. (6) (a) Hardacre, C.; Holbrey, J. D.; MeMath, S. E. J.; Bowron, D. T.; Soper, A. K. J. Chem. Phys. 2003, 118, 273. (b) Hardacre, C.; MeMath, S. E. J.; Nieuwenhuyzen, M.; Bowron, D. T.; Soper, A. K. J. Phys. Condens. Matter 2003, 15, S159. (c) Deetlefs, M.; Hardacre, C.; Nieuwenhuyzen, M.; Padua, A. A. H.; Sheppard, O.; Soper, A. K. J. Phys. Chem. B 2006, 110, 12055. (7) Every, H. A.; Bishop, A. G.; MacFarlane, D. R.; Ora¨dd, G.; Forsyth, M. Phys. Chem. Chem. Phys. 2004, 6, 1758. (8) Noda, A.; Hayamizu, K.; Watanabe, M. J. Phys. Chem. B 2001, 105, 4603. (9) (a) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. J. Phys. Chem. B 2004, 108, 16593. (b) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. J. Phys. Chem. B 2005, 109, 6103. (c) Tokuda, H.; Ishii, K.; Susan, M. A. B. H.; Tsuzuki, S.; Hayamizu, K.; Watanabe, M. J. Phys. Chem. B 2006, 110, 2833. (d) Tokuda, H.; Tsuzuki, S.; Susan, M. A. B. H.; Hayamizu, K.; Watanabe, M. J. Phys. Chem. B 2006, 110, 19593. (10) Bonhoˆte, P.; Dias, A.; Papageorigiou, Kalyanasundaram, K.; Gra¨tzel, M. Inorg. Chem. 1996, 35, 1168. (11) Fannin, A. A., Jr.; Floreani, J. D. A.; King, L. A.; Landers, J. S.; Piersma, B. J.; Stech, D. J.; Vaughn, R. L.; Wilkes, J. S.; Williams, J. L. J. Phys. Chem. 1984, 88, 2614. (12) (a) Wakai, C.; Oleinkova, A.; Ott, M.; Weinga¨rtner, H. J. Phys. Chem. B 2005, 109, 17028. (b) Bright, F. V.; Baker, G. J. Phys. Chem. B 2005, 110, 5822. (c) Wakai, C.; Oleinikova, O.; Weinga¨rtner, H. J. Phys. Chem. B 2005, 110, 5824. (d) Daguenet, C.; Dyson, P. J.; Krossing, I.; Oleinikova, A.; Slattery, J.; Wakai, C. Weinga¨rtner, H. J. Phys. Chem. B 2006, 110, 12682. (13) Schro¨dle, S.; Annat, G.; MacFarlane, D. R.; Forsyth, M.; Buchner, R.; Hefter, G. Chem. Commun. 2006, 1748. (14) Schro¨der, C.; Wakai, C.; Weinga¨rtner, H.; Steinhauser, O. J. Chem. Phys. 2007, 126, 084511. (15) Weinga¨rtner, H.; Sasisanker, P.; Daguenet, C.; Dyson, P. J.; Krossing, I.; Slattery, J. M.; Schubert, T. J. Phys. Chem. B 2007, 111, 4775. (16) Triolo, A.; Russina, O.; Bleif, H.; Cola, E. D. J. Phys. Chem. B 2007, 111, 4641. (17) Rey-Castro, C.; Vega, L. F. J. Phys. Chem. B 2006, 110, 14426. (18) Salanne, M.; Simon, C.; Turq, P.; Madden, P. A. J. Phys. Chem. B 2007, 111, 4678. (19) Krossing, I.; Slattery, J. M.; Daguenet, C.; Dyson, P. J.; Oleinikova, A.; Weinga¨rtner, H. J. Am. Chem. Soc. 2006, 128, 13427. (20) Abbott, A. P.; Harris, R. C.; Ryder, K. S. J. Phys Chem. B 2007, 4910. (21) Aerov, A. A.; Khokhlov, A. R.; Potemkin, I. I. J. Phys Chem. B 2007, 111, 3462. (22) (a) Shah, J. K.; Frennecke, J. F.; Maginn, E. J. Green. Chem. 2002, 4, 112. (b) Morrow, T. I.; Maginn, E. J. J. Phys. Chem. B 2002, 106, 12807. (c) Cadena, C.; Zhao, Q.; Snull, R. Q.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 2821. (d) Cadena, C.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 18026. (e) Kelkar, M. S.; Maginn, E. J. J. Phys. Chem. B 2007, 111, 4867. (23) (a) de Andrade, J.; Boes, E. V.; Stassen, H. J. Phys. Chem. B 2002, 106, 3546. (b) de Andrade, J.; Boes, E. V.; Stassen, H. J. Phys. Chem. B 2002, 106, 13344. (24) Tsuzuki, S.; Tokuda, H.; Hayamizu, K.; Watanabe, M. J. Phys. Chem. B 2005, 109, 16474. (25) Wang, Y.; Izvekov, S.; Yan, T.; Voth, G. A. J. Phys. Chem. B 2006, 110, 3564. (26) Micaelo, N. M.; Baptista, A. M.; Soares, C. M. J. Phys. Chem. B 2006, 110, 14444. (27) (a) Liu, Z.; Huang, S.; Wang, W. J. Phys. Chem. B 2004, 108, 12978. (b) Liu, Z.; Huang, S.; Wang, W. Chem. Phys. Chem. Phys. 2006, 8, 1096. (28) (a) Lopes, J. N. C.; Deschamaps, J.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 2038. (b) Lopes, J. N. C.; Deschamaps, J.; Pa´dua, A. A. H.

J. Phys. Chem. B, Vol. 112, No. 6, 2008 1751 J. Phys. Chem. B 2004, 108, 11250. (c) Lopes, J. N. C.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 16893. (d) Lopes, J. N. C.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 7485. (e) Lopes, J. N. C.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 19586. (29) (a) Bhargava, B. L.; Balasubramanian, S. J. Chem. Phys. 2005, 123, 144505. (b) Bhargava, B. L.; Balasubramanian, S. J. Chem. Phys. 2006, 125, 219901. (c) Bhargava, B. L.; Balasubramanian, S. Chem. Phys. Lett. 2006, 417, 486. (d) Bhargava, B. L.; Balasubramanian, S. J. Phys. Chem. B 2007, 111, 4477. (30) Alavi, S.; Thompson, D. L. J. Chem. Phys. 2005, 122, 154704. (31) Schro¨der, C.; Rudas, T.; Steinhauser, O. J. Chem. Phys. 2006, 125, 244506. (32) Hu, Z.; Margulis, C. J. Proc. Natl. Acad. Sci. 2006, 103, 831. (33) Dymek, J. R.; C. J.; Grossie, D. A.; Fratini, A. V.; Adams, W. W. J. Mol. Struct. 1989, 213, 25. (34) Dong, K.; Zhang, S.; Yao, X. J. Phys. Chem. A 2006, 110, 9775. (35) Ko¨ddermann, T.; Wertz, C.; Heintz, A.; Ludwig, R. Chem. Phys. Chem. 2006, 7, 1944. (36) Del Po´polo, M. G.; Lynden-Bell, R. M.; Kohanoff, J. J. Phys. Chem. B 2005, 109, 5895. (37) Waters, M. L. Curr. Opin. Chem. Biol. 2002, 6, 736. (38) Sony, S. M. M.; Ponnuswamy, M. N. Cryst. Growth Des. 2006, 6, 736. (39) Downard, A.; Earle, M. J.; Hardacre, C.; Mcmath, S. E. J.; Nieuwenhuyzen, M.; Teat, S. J. Chem. Mater. 2004, 16, 43. (40) Saha, S.; Hayashi, S.; Kobayashi, A.; Hamaguchi, H. Chem. Lett. 2003, 32, 740. (41) Holbrey, J. D.; Reichert, W. M.; Rogers, R. D. Dalton Trans. 2004, 2267. (42) Holbrey, J. D.; Reichert, W. M.; Nieuweihuyzen, M.; Sheppard, O.; Hardcare, C.; Rogers, R. Chem. Commun. 2003, 476. (43) Lachwa, J.; Bento, I; Duarte, M. T.; Lopes, J. N. C.; Rebelo, L. P. N. Chem. Commun. 2006, 2445. (44) Lopes, J. N. C.; Pa´dua, A. A. H. J. Phys. Chem. 2006, 110, 3330. (45) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Comm. 1995, 91, 43. (46) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306. (47) Berendsen, H.; Postma, J.; DiNola, A.; Haak, J. J. Chem. Phys. 1984, 81, 3684. (48) Levitt, M.; Hirshberg, M.; Sharon, R.; Daggett, V. Comput. Phys. Commun. 1995, 91, 215. (49) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (50) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (51) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (52) (a) Urahata, S. M.; Ribeiro, M. C. C. J. Chem. Phys. 2004, 120, 1855. (b) Urahata, S. M.; Ribeiro, M. C. C. J. Chem. Phys. 2005, 122, 024511. (53) Del Po´polo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744. (54) Krummen, M.; Wasserscheid, P.; Gmehling, J. J. Chem. Eng. Data 2002, 47, 1411. (55) Hanke, C. G.; Prince, S. L.; Lynden-Bell, R. M. Mol. Phys. 2001, 99, 801. (56) Baaden, M.; Burgard, M.; Wipff, G. J. Phys. Chem. B 2001, 105, 11131. (57) Qiao, B.; Krekeler, C.; Berger, R.; Site L. D.; Holm, C., in preparation. (58) Laaksonen, L. J. Mol. Graphics 1992, 10, 33. (59) Bondi, A. J. Phys. Chem. 1964, 68, 441. (60) Hunt, P. A. J. Phys. Chem. B 2007, 111, 4844. (61) Luzar, A.; Chandler, D. J. Chem. Phys. 1993, 98 (10), 8160. (62) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Timneanu, N. J. Phys. Chem. B 2006, 110, 4393. (63) (a) Luzar, A.; Chandler, D. Nature 1996, 379, 55-57. (b) Luzar, A. J. Chem. Phys. 2000, 113, 10663. (64) Starr, F. W.; Nielsen, J. N.; Stanley, H. E. Phys. ReV. E 2000, 62, 579. (65) Fujii, K.; Fujimori, T.; Takamuku, T.; Kanzaki, R.; Umebayashi, Y.; Ishiguro, S. J. Phys. Chem. B 2006, 110, 8179.