Interaction of Water with the Model Ionic Liquid [bmim][BF4]: Molecular

Jun 7, 2008 - At high water content, the ionic network is somewhat disrupted or ... Ionic Liquids As Tunable Toxicity Storage Media for Sustainable Ch...
0 downloads 5 Views 925KB Size
Download PDF
7826

J. Phys. Chem. B 2008, 112, 7826–7836

Interaction of Water with the Model Ionic Liquid [bmim][BF4]: Molecular Dynamics Simulations and Comparison with NMR Data Margherita Moreno, Franca Castiglione, Andrea Mele,* Carlo Pasqui, and Guido Raos* Dipartimento di Chimica, Materiali e Ingegneria Chimica “G. Natta”, Politecnico di Milano, 20131 Milano, Italy ReceiVed: January 15, 2008; ReVised Manuscript ReceiVed: April 3, 2008

We report on molecular dynamics simulations of the ionic liquid [bmim][BF4] and its mixtures with water, from zero up to 0.5 mol fraction of water. All of the simulations are carried out with two published force fields. The results are compared with each other and with published as well as new NMR data on the same mixtures, whenever possible. We perform extensive analyses of structural quantities, such as pair correlation functions, nearest-neighbor analysis and size distribution of the water clusters formed at higher concentrations. We show that the water clusters are formed almost exclusively by linear chains of hydrogen-bonded molecules. There is a nanoscale structuring of the mixtures but no macroscopic phase separation among the components, in agreement with experiment. Roughly, we identify two solvation regimes. At low water content, the ions are selectively coordinated by individual water molecules, but their ionic network is largely unperturbed. At high water content, the ionic network is somewhat disrupted or swollen in a nonspecific way by the water clusters. 1. Introduction A widely used and studied class of room-temperature ionic liquids (RTILs) is based on 1-alkyl-3-methylimidazolium (Cnmim) salts with a large repertoire of anions, PF6-, (CF3SO2)2N-, and BF4- being the most popular. The need for a better understanding of physicochemical processes in RTILs has stimulated many studies on the structure and the solvation mechanisms of molecular solutes in these systems.1–7 Some kind of structural organization in these liquids is expected on the basis of simple but perhaps also simplistic considerations. The imidazolium cation can be sketched as a polar head (the ring) with an apolar tail (the alkyl chain, Cn). The alkyl tails will tend to aggregate when they exceed some critical length, in analogy with other amphiphilic systems. This aggregation is not only evident in aqueous solutions but also evident in the neat ionic liquids, where a smectic liquid-crystalline phase has been reported for Cnmim-based RTILs with n g 10.8 This tendency to nanoscale segregation must be balanced against the requirement of even interdispersion of positive and negative charges within the liquid. Finally, these broad features may also be modulated by other specific interactions, such as hydrogen bonding between suitable sites of the cations, the anions, and possible solutes. Understanding the complex interplay among these factors represents an intriguing structural problem. A number of investigations by neutron and X-ray scattering,9–12 IR-Raman spectroscopy,13–16 NMR17,18 and molecular dynamics (MD) simulation19–22 have confirmed the presence of some mesoscale order or nanostructural organization in these systems. Information on the local organization of bulk ionic liquids was recently obtained by NMR methods based on the detection of intermolecular nuclear Overhauser enhancement in both rotating frame (ROE) and laboratory frame (NOE).23–25 NMR studies were carried out on the model system 1-butyl-3methylimidazolium tetrafluoroborate [bmim][BF4] (see Scheme * Corresponding author. E-mail: [email protected] (A.M.) and [email protected] (G.R.).

SCHEME 1: Chemical Structure and Atom Numbering of the Ions of [bmim][BF4]a

a Within the text, label Hn indicates the H atom(s) covalently bound to carbon atom n (e.g., H10 for three hydrogens bonded to C10). The oxygen and hydrogen atoms of water are denoted by OW and HW, respectively. The greek letters indicate the n-butyl torsion angles.

1 for molecula formula and atom numbering), either neat or with variable amounts of added water. Water was studied as an important representative model of a polar solute dissolved in RTILs. The collected data concerned cation-cation interactions (obtained by the method of Osteryoung and co-workers),26 water-cation, water-anion and cation-anion interactions. The first attempt of semiquantitative interpretation of intermolecular NOEs24 was carried out under the reasonable assumption of long-lived aggregates and thus neglecting, at this stage of approximation, long-range and nonspecific contributions to relaxation arising from translational diffusion.27–29 NOE data provided (a) distance restraints for cation-cation organization, (b) details of the head-to-tail interactions, and (c) information on the mode of interaction of the cations with the water molecules. The results of our NMR studies may be summarized as follows: (i) There is significant cation-cation aggregation (headto-head) in the pure liquid. Some head-to-head intermolecular distances are comparable to those found in the solid state, at least in those cases where X-ray crystal structures are available. (ii) Ion pairing is mediated by C-H···F hydrogen bonds. (iii) Up to 1:1 molar fraction, added water binds both to the anion and to the cation. The imidazolium-imidazolium interactions are loosened by the presence of water, which acts as hydrogenbond competitor: C-H···F interactions are progressively replaced

10.1021/jp800383g CCC: $40.75  2008 American Chemical Society Published on Web 06/07/2008

Simulation of Water in [bmim][BF4] by C-H···O hydrogen bonds. At low water content, such interaction is selectively localized around the H2-H4-H5 atoms of the cation, thus confirming the role played by C-H···O type hydrogen bonds. (iv) Apolar groups (the n-butyl chains) interact with the polar imidazolium cores, so that they cannot be fully segregated in a micelle-type domain. The NMR measurement did not allow us to investigate the organization (if any) of the apolar domains and the state of water dissolved in IL. Furthermore, the local structural motifs derived from the NOE distance restraints were only approximate. The present contribution has three main objectives. First, our MD simulations on [bmim][BF4] and its mixtures with water are aimed at providing a picture of the short-to-medium range organization of the IL consistent with the experimental data, investigating the structure of the apolar domains and role of hydrogen bonding within the bulk liquid, using water as probe for this type of interaction. Second, we provide a detailed comparison of simulation results from two distinct force fields (FF) for the representation of intra- and intermolecular forces, developed independently by different groups.30–33 We are not aware of similar systematic comparisons in this context, despite several valuable MD studies on these and related systems.19,20,34–37 The value of such systematic comparisons has become appreciated in the biomolecular context.38 By highlighting the relative strengths and weaknesses of different parametrizations, they have shown the way for their improvement (e.g., subsequent generations of the CHARMM, AMBER, and OPLS FFs for proteins) and at times demonstrated the need to overcome the FF paradigm by the introduction of polarization terms or by a suitable electronic structure method (semiempirical or ab initio Hartree-Fock and density functional theories). Finally, we have also reanalyzed some of our original NMR data in order to sharpen our comparison with the MD simulation results. 2. Experimental Methods and Computational Approach Materials and NMR Measurements. Samples of [bmim][BF4] were prepared and purified according to standard literature methods.39 The samples of [bmim][BF4] with water were prepared by mixing weighted amounts of [bmim][BF4] with HPLC grade water. The water quantity was measured by microsyringe, assuming d(H2O) ) 1.00 g mL-1. The mixtures were tranferred to 5 mm NMR tubes equipped with a coaxial capillary containing d6-DMSO with tetramethylsilane (TMS) as internal chemical shift reference. 1H NMR spectra were collected on a Bruker Avance 500 spectrometer operating at 500 MHz proton frequency. All of the spectra were carried out at 305.0 ( 0.1 K. The concentration range explored was water molar fractions χw ) 0; 0.09; 0.17; 0.27; 0.36; 0.45; 0.52, corresponding to the following concentrations of compound [bmim][BF4] (mM): 5354, 5302, 5251, 5171, 5077, 4968, and 4844, respectively. The conformational analysis of the butyl side chain was carried out by considering the torsion angles β and γ, defined as the torsions across the C(6)-C(7) and C(7)-C(8) chemical bonds, respectively (see Scheme 1 for atom numbering). These angles, in turn, are related to the values of NMR coupling constants 3JH6-H7 and 3JH7-H8. The populations of the trans (t) and gauche (g) rotamers for both torsion angles were obtained according to the classic approach proposed by Abraham et al.40 Simulations. The interaction of water with the ionic liquid was investigated by the molecular dynamics (MD) method. All simulations were performed with the TINKER 4.2 molecular modeling package.41 We carried out two independent sets of simulations, using two different force fields. The first force field,

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7827 TABLE 1: Composition of the Simulated Systems system pure liquid 1:8 mixture 1:4 mixture 1:2 mixture 1:1 mixture

no. of no. of water water mole fraction (χW) ion pairs molecules 0.000 0.111 0.200 0.333 0.500

125 125 125 125 125

0 15 31 62 125

total no. of atoms in the simulation box 3750 3795 3843 3936 4125

named FF1 in this paper, is based on the AMBER-type parametrization by Wang and co-workers30 for both the [bmim]+ cation and the BF4- anion. The second force field, FF2, is based on the mixed OPLS/AMBER-style parametrization for the cation by Canongia Lopes and co-workers,31,32 and on AMBER-style anion parameters from Stassen and co-workers.33 In both cases, water molecules were simulated using the TIP3P model.42b Simulations were carried out under periodic boundary conditions, computing the long-range electrostatic interactions with the particle mesh Ewald (PME) method.43 Both the cutoff for van der Waals interactions and the switching distance for electrostatic interactions were fixed at the value of 12.0 Å. We prepared five model systems, corresponding to five different compositions. These are given in Table 1, both as water mole fraction and as water/IL ratio. Our simulation strategy consisted of two steps. First, we equilibrated the systems by constant pressure and temperature MD runs at T ) 310 K and P ) 1.0 atm, which were controlled by the Berendsen method44 (the coupling times are 0.2 ps for the pressure and 2.0 ps for the temperature). The MD equations were integrated with the velocity-Verlet scheme and a 2.0 fs time step, with fixed C-H bond lengths. Each equilibration run started from a low-density condition and typically lasted for 0.7-0.9 ns, until attainment of a roughly constant density. The final densities are compared with the available experimental values45,46 in Table 2 (see also ref 47 for further data). The second step consisted of MD runs at constant volume and temperature (310 K), starting from the equilibrated systems. The length of each of these production runs was between 8 and 10 ns. Atomic coordinates were saved to disk every 10.0 ps for later analysis. 3. Results and Discussion In our discussion of the simulation results and comparison with experiment, we shall start from basic physicochemical properties (e.g., densities and heats of mixing). We shall then consider increasingly complex and detailed features, such as the analysis of molecular conformation and specific intermolecular interactions. Most of these properties provide a statistical characterization of the equilibrium system configuration. However, we shall also discuss a series of simple but important dynamical quantities, namely, the diffusion coefficients of the ions and water and the time correlation functions defining the average lifetime of the water clusters. For the sake of conciseness, several plots and tables contain the results from a single FF. We report the data from both FFs only in the presence of significant differences or for the sake of comparison. 3.1. Thermodynamic Properties. The densities of the pure ionic liquids and their water mixtures are plotted in Figure 1 (see also Table 2). FF1 systematically predicts lower densities than FF2. These differences are not large, being smaller than 1%. In fact, both force fields predict somewhat higher densities than the experimental ones at 308 K45 (see again Table 2). The overall trends for the mixtures are identical for both force fields:

7828 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Moreno et al.

TABLE 2: Densities of the MD Simulations (in g/cm3) and Comparison with Experimental Data simulationsa

experimentb

water mole fraction (χW)

FF1

FF2

data from ref 45, except where indicated

0.000

1.1886 ( 0.0003

1.1974 ( 0.0008

0.111 0.200 0.333 0.500

1.1901 ( 0.0002 1.1852 ( 0.0004 1.1846 ( 0.0003 1.1795 ( 0.0004

1.1939 ( 0.0003 1.1955 ( 0.0004 1.1924 ( 0.0005 1.1871 ( 0.0003

1.1748 (308), 1.1940 ( 0.0014 (313)c, 1.19423 (313)d, 1.1822 (298) 1.19115 (313)d, 1.1792 (298) 1.18782 (313)d, 1.1761 (298) 1.18124 (313)d, 1.1701 (298) 1.17025 (313)d, 1.1601 (298)

a Calculated values are averages over the last part of the equilibration runs. The uncertainties correspond to twice the standard error. Experimental temperatures are reported in brackets. The data at χW > 0 have been obtained by a parabolic fit to the data tabulated by Li et al.45 and Rebelo et al.48 c Brennecke and co-workers.46 d Rebelo et al.48 b

Figure 1. Computed densities of the pure liquids and mixtures, as a function of water mole fraction (black and red symbols). The continuous lines are parabolic fits to the numerical data. The blue diamonds and line represent experimental data at 313 K from ref 48. The dashed lines represent the “ideal” mixture results, obtained under the assumption of volume additivity. The bulk density of the TIP3P water model has been obtained by an independent MD simulation, and is 0.9886 g/cm3 at 310 K and 1 atm (expt value 0.9933 g/cm3; see ref 49).

Figure 3. Representative snapshots from the MD simulations of [bmim][BF4]/water systems using FF1: (a) 1:8 mixture, (b) 1:4 mixture, (c) 1:2 mixture, (d) 1:1 mixture. Color code: water in green, cations in red, anions in blue. Hydrogens are not shown, for clarity. The figures have been generated with VMD.52

Figure 2. Computed energies of mixing for the mixtures [bmim][BF4]/ water, as a function of the water mole fraction.

the densities decrease in an approximately parabolic fashion with χW, but they are larger than the values calculated assuming volume additivity (dashed lines in Figure 1). Thus, the simulations predict negative excess molar volumes for the mixtures. This disagrees with the experimental measurements by both Li et al.45 and Rebelo et al.,48 who reported a positive volume change upon mixing. The excess internal energies of mixing obtained from the MD simulations of the pure liquids and mixtures are plotted in Figure 2. The results from FF1 and FF2 are almost identical. In this case, there is good agreement with the experimental measurements by Rebelo et al.48 Their mixing enthalpies are also positive, and they exhibit a maximum of about 2.2 kJ/mol at χW = 0.6. Our maximum is slightly displaced, being located somewhere around χW ) 0.4, but the value of the maximum mixing energy is almost quantitatively correct (for condensed phases, the PV term is negligible compared with the internal energy U, so that the latter may be safely identified with the enthalpy H).

In summary, the volumes and internal energies of mixing provide simple but very valuable tests of the quality of the force fields and of the overall simulation strategy. Our results are encouraging for the energies but disappointing for the volumes. It is unlikely that the results would be entirely reversed by carrying out longer equilibration runs. Rather, we believe that volume change upon mixing is a sensitive probe of the force field quality, and the present results point to the possibility (and the need) of further fine-tuning of their parameters. In fact, both the RTIL and the water force fields were obtained and tested by simulations of the pure components. There is probably room for improvement in the treatment of water-IL cross interactions. In particular, polarization effects22 may be important, and the Lennard-Jones parameters do not necessarily follow any simple set of “mixing rules”, neither AMBER-style nor OPLS-style.38,50 Before embarking on such a task, however, it would be wise to consider also other water models. Examples are TIP4P and TIP5P,42,51 which provide improved estimates of both static and dynamic properties of bulk water (including the temperatureand pressure-dependence of its density) at a relatively modest computational cost. We did not do this in this work, since our main focus was the comparison of the RTIL force fields. 3.2. Pair Distribution Functions and Coordination Numbers. Figure 3 contains snapshots of the final system configurations, obtained from the MD production runs of the mixtures. They show a clear tendency of the water molecules to aggregate

Simulation of Water in [bmim][BF4]

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7829

Figure 4. Radial distribution functions between the anion (B) and the ring hydrogens (H2, H4, H5) of [bmim][BF4]. Data are from the simulations with FF1 of the neat RTIL (continuous black lines) and of the 1:1 mixture (red dashed lines).

at higher concentrations. Still, the water clusters remain limited in size, and they do not demix to form a two-phase system, in agreement with the hydrophilic nature of [bmim][BF4]. The radial (pair) distribution functions gAB(r) are a classic tool to characterize the structure of the liquids and their mixture, through the effective pairwise interaction between atoms A and B. Because of the large number of possible atom pairs, one is almost forced to select a relatively small subset among them. Figure 4 shows that there is a certain weakening of the interaction between the anion and the cation (represented by H2, H4, H5, the three hydrogen bond donors on the fivemembered ring) on going from the neat RTIL to the water mixtures. The first maximum of these pair distribution functions is shifted to slightly smaller values at slightly larger distances. The second maximum is almost unaffected, meaning that water perturbs the RTIL network only at short range. Hydrogen H2 is the one which appears to interact more closely with the anion. This is consistent with their charge distributions within the FF. Also the ab initio calculations (MP2/6-311G**) on selected gas-phase ion pairs by Tsuzuki et al.53 confirm that H2 interacts more strongly than H4 or H5. Conversely, the B...H4 and B...H5 interactions are weaker and also more liable to be affected by water. Figure 5 contains further plots of some gAB(r) from the FF1 simulation of the 1:1 sample. To a smaller degree, the interaction between the rings hydrogens and the oxygen atom of water (OW) follows the previous trends for cation-anion interactions (H2 > H4 ≈ H5). Thus, these groups have indeed a propensity to become engaged as hydrogen bond donors with the anion and water. However, a water molecule is less selective than the anion in choosing its interaction partner, probably also because of its smaller size. Water can also act as a hydrogen bond donor toward the anion, as shown by the OW...B pair distribution function (we have plotted this in place of HW...F to avoid multiple complicated oscillations in g(r), arising from the fourfold and twofold “degeneracies” of F and HW). Finally, the large and well-defined short-range peak of the OW...OW pair distribution function provides a first quantitative confirmation of the existence of the water clusters seen in Figure 3. These clusters will be further characterized below. To quantify better and gain an overall view of the structural trends produced by water in all the simulated systems, we have computed a series of coordination numbers. These represent the average number of atoms within a prescribed cutoff distance from the given reference atom. They are obtained by numerical integration of the pair distribution functions. After proper normalization, gXY(r) yields both NX(Y) (number of X-type atoms surrounding Y) and NY(X) (number of Y-type atoms surrounding X). These coincide only when the total number of Xs and Ys are identical. For simplicity and ease of comparison,

Figure 5. Selected radial distribution functions of the 1:1 mixture, obtained from the MD simulations with force field FF1. (a) OW...H2 continuous black line, OW...H4 dashed red line, OW...H5 dash-dot blue line. (b) OW...OW continuous black line, OW...B dashed red line.

we have chosen the cutoff of 4.5 Å for all atom pairs, but in some cases, we have also performed additional analyses with different cutoffs. Figure 6 shows the systematic variations of NH2(B)[)NB(H2)] as the most important cation-anion interaction, of NB(OW) and NH2(OW) as representatives of water-ion interactions, and NOW(OW) for water-water interactions. The graphs in part (a) of the figure clearly show that, for both force fields, there is a loosening of the cation-anion interactions on increasing the water content. The plots obtained with the 5.0 Å cutoff, which corresponds to the first minimum of the gBH2(r) distribution function in Figure 4, shows that this trend is robust. NH4(B) and NH5(B) (not shown) are always lower than NH2(B), but they have the same behavior. The value of about 2 for this coordination number in neat [bmim][BF4] is consistent with other simulations and NMR experiments.36 Naively, one might expect an increase of water-cation and water-anion coordination numbers, to compensate for the decrease in the anion-cation coordination numbers. Figure 6b shows that this is not the case: both NB(OW) and NH2(OW), representing the average number of ions surrounding a water molecule, decreases upon going from the 1:8 to the 1:1 mixture (we have not plotted NOW(B) nor NOW(H2) because these increase almost trivially with χW). The steep increase in NOW(OW) demonstrates once more that the water molecules tend to aggregate with increasing volume fraction. 3.3. Clustering of Water Molecules. We have pointed out that visual inspection of the simulation snapshots (Figure 3) and analysis of the pair correlations hints at the formation of water clusters at increasing χW. These clusters were identified and analyzed quantitatively by suitable modifications of the TINKER routines. Two water molecules were considered to be connected whenever the O...O distance fell below a cutoff radius of 3.5 Å (fully beyond the sharp peak in their pair distribution function; see Figure 5). The list of intermolecular bonds was then converted into a list of connected clusters. Let p(n) be the probability of existence of a cluster of n water molecules (ncluster for short). This was calculated by counting and then

7830 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Moreno et al.

Figure 8. Equilibrium constants for the formation of water clusters of size n, obtained from the analysis of the MD simulations of the 1:1 mixtures.

1:1 mixture, there is a small but appreciable probability of finding sizable water clusters (even larger than the 10-molecule cutoff used in the figure). Somewhere around the 1:4 composition, there is a transition, from a regime with water monomers and dimers to a regime where larger aggregates become increasingly likely. It is also of interest to investigate the structure of these clusters, that is, whether they resemble roughly spherical and compact droplets, linearly connected chains, or some intermediate situation between these extremes. This has been done through the following connectivity index: ∞

C)

〈N 〉

∑ p(n) n -b 1n

(2)

n)3

Figure 6. Selected coordination numbers (see the main text for the definition), obtained from the pair distribution functions of [bmim][BF4] and its mixtures with both FF1 and FF2. The cutoffs are 4.5 Å, except where indicated.

Figure 7. Normalized histograms of P(n), obtained from the MD simulations with FF1.

normalizing the number of their occurrences over a whole MD run. A related quantity is the probability P(n) that a given water molecule belongs to an n-cluster:

P(n) )

∑

np(n) ∞ mp(m) m)1

(1)

Compared with p(n), P(n) has the advantage of being more sensitive to the formation of large clusters. The P(n) probabilities are plotted in Figure 7, for each water mole fraction. The simulations with FF1 (shown in the Figure) and with FF2 yield similar results. The probability of isolated water molecules (1-clusters) decreases with increasing water content, while the probabilities of larger clusters grow. In the

In eq 2, 〈Nb〉n is the average number of intermolecular “bonds” among water molecules within an n-cluster. The C index is designed to yield different values for different cluster connectivities. C ) 1 corresponds to a linear association of water molecules (two bonds in a trimer, three in a tetramer, etc.). A circularly connected trimer would have C ) 3/2; a circular tetramer would have C ) 4/3, etc. The upper limit for C is obtained for an infinite cluster of tetrahedrally coordinated water molecules (as in bulk water), where C ) 2. Note that clusters with n ) 2 have been excluded from the overall average, since they are trivially linear. The calculated C index was always extremely close to unity with both FFs (1:1 mixture FF1, C ) 1.00000; 1:1 mixture FF2, C ) 1.00008), meaning that the formation of linear clusters is highly favored compared with other types of aggregation. This conclusion is consistent with the water coordination numbers shown in Figure 6, which are of the order of unity even in the 1:1 mixture. The cluster probabilities p(n) are directly proportional to the clusters’ molar concentrations [(H2O)n], which enter the phenomenological equilibrium constants Kn (n ) 1, 2, 3, ...) defined by eq 3a,b:

(H2O)n + H2O h (H2O)n+1 Kn )

[(H2O)n+1] [(H2O)n][H2O]

(3.a) (3.b)

Thus, K1 is the constant for formation of a dimer from two isolated molecules; K2 is the constant for formation of the trimers from a dimer and a molecule, and so forth. The results of the Kn calculations from both FFs are summarized in Figure 8. They are obtained from the simulations with χW ) 0.5, which provide the most reliable results: the large number of n-clusters with n e 5 allows the collection of good statistics in this case. Within the error bars, the equilibrium constants for formation of the trimers, tetramers, pentamers, and so forth are all identical. Their n independence signifies that the free energy changes on going

Simulation of Water in [bmim][BF4]

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7831

Figure 9. Normalized histogram of nearest-neighbor interactions experienced by the anion (B atom). Left-hand side: force field FF1. Right-hand side: force field FF2.

Figure 10. Normalized histogram of nearest-neighbor interactions experienced by water (OW atom). Left-hand side: force field FF1. Right-hand side: force field FF2.

TABLE 3: Results of Linear Regression (y ) mx + q) for Selected Pair of Protons of [bmim][BF4] at Increasing Water Molar Fractionsa entry

y

x

m

q

R2

1 2 3 4 5 6 7

δH2 δH2 δH4 δH2 δH2 δH9 δH10

δH4 δH5 δH5 δH6 δH10 δH8 δH9

2.38 1.98 0.83 7.73 11.24 0.99 -1.33

-8.30 -5.59 1.13 -19.49 -28.92 -0.43 3.53

0.999 0.999 >0.999 0.998 0.996 0.992 0.916

a

See section 2 for concentration of the seven samples.

from dimer to trimer to tetramer and so forth are roughly identical. This is consistent with the linear connectivity of the clusters: the “incoming” molecule experiences a single interaction with the “growing end” of the cluster. The only exception to this scenario is represented by dimer. The lower value of K1 is most likely due to a larger entropy loss in this case. The dynamics of cluster formation and breakup will be discussed in section 3.6. The results guarantee that the equilibria (eq 3.a) are effectively achieved over the time scale of our simulations. 3.4. Statistics of Nearest-Neighbor Contacts and Comparison with NMR Data. The statistics of nearest-neighbor contacts are a sensitive structural indicator related to NMR NOE intensity, as the set of nearest-neighbors of a given nucleus includes the subset of nuclei whose distances are below the threshold for observable NOE, in the approximation of intermolecular cross-relaxation in long-lived associated systems (see above). The nearest-neighbor analysis consists of determining how many times a certain atom is the nearest-neighbor of a given reference one. To achieve that we developed a special purpose routine within the TINKER package. The frequency of atom X as nearest-neighbor of a reference atom Y was

normalized with respect to the total number of Y atoms within a single frame and the total number of frames. For simplicity, the reference atom Y was always chosen among “backbone” atoms (B, O, C, N). Figure 9 displays the normalized histograms of the nearestneighbors of the B atom with the H atoms of water and [bmim]+ at increasing water molar fraction, computed with both FF1 and FF2. Thus, they represent the anion-water and anion-cation interactions, and they can be compared with the 1H-1H and 1H-19F NOE data. The first set of histogram bars (labeled as b-hw) clearly indicates again water-anion interactions mediated by O-H...F-B hydrogen bonds. This is in good agreement with the observed heteronuclear NOE between the water hydrogens and the F nuclei of BF4-, as reported for the sample at χW ) 0.36 in ref 23. The percentage of HW atoms which are first neighbors of B rapidly increases with increasing χW, as expected on the basis of the strong affinity of water for the anion.14 The contact statistics of the B atom with the [bmim]+ protons (categories b-h2 to b-h10) accounts for the formation of ion pairs and shows their stability in the presence of water. Noticeably, the data clearly indicate that the anion contacts selectively the H2-H4-H5-H10 subset of [bmim]+, which belong to the polar end of the cation. Both FF1 and FF2 lead to the conclusion that the protons of n-butyl chain are quite insensitive to the presence of BF4-, while H2, H4, and H5 are the major sites of interaction with the anion, in good agreement with the H-bond donor properties of C(sp2)-Hs of the imidazolium ring and the stabilization of ion pairs through C(sp2)-H...F-B hydrogen bonds. Even in the presence of equimolar amounts of water, the histograms show significant aggregation of anion and cation, thus confirming the attitude of the ion pairs to survive also in the presence of a strong H-bond donor/acceptor species such as water. These conclusions

7832 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Moreno et al.

TABLE 4: Conformational Preferences of n-Butyl Side Chain in 1-Butyl-3-methylimidazolium and 1-Butyl-2,3-dimethylimidazolium Ions in the Solid State system R(°) β(°) γ(°) Conf. typee ref.

bmimCla bmimClb bmimBr bmimBrc bmimPF6 79 68 180 (u, g+, t) 59

-83 -174 -175 (d, t, t) 60

78 66 180 (u, g+, t) 61

-124 171 194 (d, t, t) 61

104 -61 178 (u, g-, t) 62

bmimOTfd 102 (-106) -68 (66) 178 (180) (u, g-, t), (d, g+, t) 62

bdmimBF4 bdmimCl bdmimPF6 bdmimSbF6 bdmimHSO4 -97 -175 -180 (d, t, t) 63

75 66 180 (u, g+, t) 64

80 62 180 (u, g+, t) 63

-81 -63 176 (d, g-, t) 63

80 58 180 (u, g+, t) 64

Orthorhombic. b Monoclinic. c Cocrystallized with 1,4-phenylendiamine. d OTf ) Trifluoromethanesulphonate. Two independent [bmim]+ cations in the crystal cell. e The conformational states are labeled as follows: for torsion angle R, the arbitrary labels “up” and “down” (u and d respectively) have been used for +90° and -90° angles. For β and γ torsion, the conformational states trans and gauche have been identified as t, g+, and g-. a

Figure 11. Dihedral angle distributions for the alkyl chains, calculated from MD of the pure ionic liquid (continuous black line) and 1:1 water mixtures (dashed red line) with both FF1 (top) and FF2 (bottom).

provide a detailed atomistic model fully consistent with that deduced on the basis of observed 1H-19F NOE.23 The histograms of Figure 10 collect the nearest neighbor statistics for the water molecules (OW atoms). The first set of bars (category ow-hw) shows that water clustering increases rapidly with χW. While the number of water molecules surrounding the anion increases with χW (Figure 9, category b-hw), the average number of anions surrounding a given water molecule follows the opposite trend (Figure 10, category owf). Thus, the MD simulations successfully reproduce competitive processes such as formation of water aggregates and water-anion adducts. The ow-h categories in Figure 10 are related to water-cation interactions. The intensities for water-cation contacts (from ow-h2 to ow-h10) clearly show that H2 and, to a lesser extent, H4 and H5 are the preferential sites for water binding. This confirms our observation of specific water-cation interactions mediated by C(sp2)-H...O-H hydrogen bonds, obtained by following the NOE intensities as a function of water content. Also the selectivity of the interaction parallels the selectivity of the previously reported intermolecular NOEs, at low water molar fractions (up to χW ) 0.2). While the two sets of MD simulations reproduce the repertoire of hydrogen bond interactions in [bmim][BF4] and their modulation by water, significant differences between FF1 and FF2 are clearly evidenced by comparing the nearest neighbor statistics of H2 with those of H4 and H5. Figure 9 shows that the B...H2 contacts predicted by FF1 are nearly 4 times more frequent than those of B...H4 and B...H5 pairs, thus suggesting

a stronger capability of H2 in coordinating the anion. On the other hand, FF2 largely smooths out the difference between the H2 and the other aromatic hydrogens. Similarly, Figure 10 shows subtle differences among the force fields in handling the interaction of water with the aromatic rings. The FF1 simulations apparently show that O...H2 contacts are roughly twice as frequent as O...H4 and O...H5, while these differences are not evident from the FF2 simulations. Conversely, FF2 seems to emphasize the role of water-anion contacts compared with FF1 (category ow-f in Figure 10). In view of these differences among the two force fields, we decided to compare their predictions with the 1H chemical shift of [bmim][BF4] protons at different χW (unpublished experimental values, listed in Table S1 of Supporting Information) that may give insights into the different capability of the imidazolium ring H atoms to interact with water molecules. The data are analyzed following the approach by Headley et al.,54,55 in order to detect the sensitivity of each proton pair to changes in the solvation medium produced by water. The results are summarized in Table 3. Entries 1 and 2 indicate that the shielding effect on H2 due to the molecular environment is roughly twice that experienced by H4 and H5, as indicated by the slope of the regression lines. These data are consistent with water molecules preferentially interacting with H2, thus indicating that the C(sp2)-H...OH interaction involving H2 is stronger than those involving the other aromatic protons, H4 and H5. The latter, in turn, shows similar sensitivity to solvation effects, as stated by the slope of δH4 versus δH5 regression line close

Simulation of Water in [bmim][BF4] to unit (entry 3). The different response to water content of H2 and of H4 and H5 agrees nicely with the simulations, expecially those using FF1. Perusal of Table 3 reveals that δH2 can also be linearly correlated to δH6 and δH10 (entries 4 and 5). The slope of the regression lines are positive, confirming that the nature of the interaction of H6 and H10 with water is of the same type of that observed for H2, H4, and H5. This must be due not only to their spatial proximity but also to the fact that C(sp3)-H bonds can act as hydrogen bond donors when the carbon is bonded to a positively charged nitrogen.56 Note also the linear correlation of δH9 with δH8 (entry 6) with slope close to unity, meaning that protons on the hydrophobic tail of [bmim]+ undergo the same effect with increasing water content. Instead, entry 7 demonstrates that the two methyl groups of [bmim]+ (H10 and H9, respectively belonging to its polar head and apolar tail) undergo opposite chemical shift variations with increasing χW. This is consistent with incipient aggregation phenomena, with possible formation of apolar domains within the RTIL-water mixtures, also discussed by Singh and Kumar in the case of a more dilute aqueous solution of [bmim][BF4].57 3.5. Conformation of the Alkyl Chains. Another point of interest for the family of 1-alkyl-3-methylimidazolium based ionic liquids is the conformation of the alkyl chain, because of the possible influence of conformational transitions on the physicochemical properties of the liquid.58 The conformation of the n-butyl chain can be described in terms of three dihedral angles: R (defined as C2-N1-C6-C7, see Scheme 1), β (N1-C6-C7-C8) and γ (C6-C7-C8-C9). The analysis of available literature data of solid state structures containing [bmim]+ and [bdmim]+ (the latter being the 1-butyl-2,3dimethylimidazolium cation) shows that only a limited number of conformations are populated in the solid state. Table 4 summarizes the experimental angles measured in the solid state via single crystal X-ray diffraction. The following general trends emerge: (i) the R angle takes values close to (90°, with maximum variation (15%; (ii) preferred values for β are close to 180° (t conformation) and (60° (g+ and g- conformation); (iii) preferred values for γ are close to 180° (t conformation) only. Accordingly, the populated conformational states can be labeled with the triplet of symbols related to dihedral angles R, β, and γ, as described in Table 4. The distribution of the R angle obtained from MD trajectories is displayed in Figure 11. In both cases, the frequency distribution of R torsion angles displays two maxima, centered on (90° for FF1 and (120° for FF2. The calculated distributions agree with the observed values in the solid state, as reported in Table 4. Unfortunately, the R angle cannot be easily obtained in the liquid via NMR methods. On the contrary, the conformational preferences corresponding to the β and γ torsion angles can be obtained at different χW by analysis of the vicinal coupling constants 3J(H7-H8) and 3J(H8-H9). The percentage of t conformation calculated from the experimental coupling constants are reported in Table 5, together with the MD values obtained by integration of the central peaks in Figure 11. Both the experimental and the simulated populations are almost insensitive to the presence of water, up to χW ) 0.5. According to the NMR measurements, the t populations are about 50% for β and 60% for γ. The t population of γ calculated with FF1 is in excellent agreement with the experimental value, whereas the t conformation of β is slightly overpopulated. The simulations based on FF2 underestimate the t population of β and overestimate that of γ.

J. Phys. Chem. B, Vol. 112, No. 26, 2008 7833 TABLE 5: Comparison of the Experimental and Calculated Population (%) of t Conformation for Torsion Angles β and γ for Both Force Fields % t rotamer β γ β γ β γ β γ β γ

experimental

calcd FF1

χW) 0 49 59 χW) 0.09 50 60 χW) 0.17 51 59 χW ) 0.36 51 59 χW ) 0.52 52 60

χW) 0 62 64 χW) 0.11 65 63 χW ) 0.20 67 62 χW ) 0.33 65 62 χW ) 0.50 64 62

calcd FF2 34 74 33 74 31 74 34 74 34 73

The differences in the side-chain conformations predicted by the two FFs are further illustrated by Figure 12, reporting the distributions of their end-to-end distances (from N1 to C9) for the pure liquid and the 1:1 water mixture. There are at least three sets of conformations, corresponding to a broad peak around 3.8 Å (folded conformation with two consecutive g states), a second peak at 4.5 Å (one g torsion), and a final one around 5.0 Å (extended all-t conformation). By integration of the area of the outermost peak in Figure 12, we find a 43% and a 32% probability for the all-t conformation respectively with FF1 and FF2. These may be compared with a recent estimate of 28%, based on Raman data.58e In summary, our NMR data in Table 5 are consistent with the Raman spectra of Berg et al.,58b and they tend to support the FF1 parametrization. However, FF258d compares better with more recent experiments.58e Further studies are necessary to achieve a quantitative understanding of the conformational equilibria in this system. 3.6. Dynamical Properties. The diffusion coefficients of the ions and water were calculated from the MD trajectories by application of the Einstein formula

Di ) lim tf∞

〈∆ri2(t)〉 6t

(4)

where 〈∆r2i (t)〉 is the mean-square displacement of species i after a time t. Because of the high viscosity and overall “sluggishness” of [bmim][BF4] and other RTILs, it is important36 to carry out relatively long MD runs (8-10 ns or more) in order to reach the true diffusive regime (〈∆r2i 〉 ∝ t) and thus estimate these quantities in a reliable way.65 Indeed, after a whole simulation run, we typically find rms values (〈∆r2i 〉1/2) for both ions of about 3 Å in the pure ionic liquid and 9 Å in the 1:1 mixture. Thus, even in the most unfavorable case, the average ion displacements are at least comparable with their size. Plots illustrating the composition dependence of our diffusion coefficients are given in Figure 13. The ion diffusion coefficients from the MD simulations of the pure liquid can be compared with published NMR data.36,47 At 310 K, Tokuda et al.47 reported that Dbmim ) 0.25 × 10-10 m2 s-1 for the cation and DBF4 ) 0.24 × 10-10 m2 s-1 for the anion. Our numbers from MD are lower by 1 order of magnitude. These results agree with those from Saielli and co-workers,36 who performed a 5 ns MD run on neat [bmim][BF4] at different temperatures. Wang and coworkers30 reported values apparently in better agreement with experiment, but these were computed from a 100 ps MD run,

7834 J. Phys. Chem. B, Vol. 112, No. 26, 2008

Moreno et al.

Figure 12. Distribution of the end-to-end distances of the alkyl chains (N1 to C9) from the simulations of the pure RTIL (continuous black line) and of the 1:1 mixtures (dashed red line) with both FF1 (left-hand side) and FF2 (right-hand side).

Figure 13. Diffusion coefficients of the ions (left) and of the water molecules (right), from the MD simulations with both FF1 and FF2.

definitely too short a time to reach the diffusive regime. Although this might be disappointing at first sight, it is also important and encouraging that the most important trend is reproduced by the simulations: in all systems, including the mixtures with water, the cation and anion diffusion coefficients are virtually identical. As also pointed out by others,66,67 this observation points to highly coordinated motion of oppositely changed ions, which may possibly form long-lived ion pairs. This holds even in the 1:1 mixture, but eventually it will break at higher water contents (remember that [bmim][BF4] and water form homogeneous mixtures over the whole concentration range). Note that the differences between the predictions of two force fields are relatively minor in this case. Water diffusion coefficients are higher by 1 order of magnitude, compared with the ions. Our diffusion coefficient for the pure TIP3P water model (DH2O ) 6.7 × 10-9 m2 s-1 at 310 K) is in line with other published simulations.42,51 This value is known to be higher than experiment by a factor of 2. More elaborate water models (TIP4P and more recently TIP5P51) have corrected this deficiency, to the point of becoming almost quantitatively correct. We believe that, also in the case of the ionic liquid, improving the comparison of computational and experimental data will require both longer simulation runs and adjustments in the force field parameters. All of the diffusion coefficients of the cation, the anion, and water are roughly constant up to χW = 0.2, after which they increase steeply. It is tempting to interpret this increase by a decrease in the overall solution viscosity, according to the Stokes-Einstein relationship Di ∝ T/(ηRi) (where T is the absolute temperature, η is the macroscopic shear viscosity of the medium, and Ri is an effective radius of species i). However, the macroscopic Stokes-Einstein picture cannot be generally valid, for at least two reasons. First of all, we have already pointed out that bmim+ and BF4- have identical diffusion properties, despite their different sizes. Second, experimental viscosity data have been collected at 298 K by Li et al.,45 who reported an exponential decrease of the viscosity from 116 mPa·s

for neat [bmim][BF4] to 60 mPa·s for the 1:4 mixture, to 21 mPa·s for the 1:1 mixture. Thus, the decrease in macroscopic viscosity does not correspond to an enhanced diffusion of the molecular components, for low water contents. Instead, we believe that the nanoscale organization of the RTIL-water mixtures may contribute to produce an anomalous diffusion behavior. We recall that the previous analysis of the structural data of the mixtures has allowed us to identify two distinct solvation regimes: a dilute regime for χW < 0.2 (1:4 mixture) with individual water molecules and a few dimers, and another regime at higher concentrations with formation of larger water clusters (see for example the discussion of Figure 7). The RTIL network is perturbed relatively little by water in the first case, while it is significantly disrupted in the second one. This, in our opinion, is also the reason for the two concentration regimes seen in Figure 13. The dynamics of formation and breakup of the water clusters can be characterized by computing the following time correlation function:68

C(t - t0) )

∑ i