Ind. Eng. Chem. Res. 2008, 47, 9115–9126
9115
Determining the Accuracy of Classical Force Fields for Ionic Liquids: Atomistic Simulation of the Thermodynamic and Transport Properties of 1-Ethyl-3-methylimidazolium Ethylsulfate ([emim][EtSO4]) and Its Mixtures with Water Manish S. Kelkar,† Wei Shi, and Edward J. Maginn* Department of Chemical and Biomolecular Engineering, UniVersity of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556-5637
The ability of simple classical force fields to predict the structure and density of ionic liquids is now wellestablished. However, it is less clear how accurate such force fields are for a range of other pure and mixture properties of ionic liquids. In this work, a single classical force field is used to compute a wide range of thermodynamic and transport properties for the ionic liquid 1-ethyl-3-methylimidazolium ethylsulfate ([emim][EtSO4]). In addition to liquid densities, the volumetric expansivity, heat capacity, enthalpy of vaporization, rotational relaxation time, self-diffusivity, shear viscosity, and thermal conductivity are computed at various temperatures for the pure ionic liquid. The density, excess molar volume, enthalpy of mixing, partial molar enthalpy, water solubility as a function of partial pressure, rotational relaxation time, self-diffusivity, shear viscosity, and thermal conductivity are also computed for mixtures that contain different concentrations of water at various temperatures. The agreement between simulations and experiment is fair for most properties, although deviations in enthalpy of mixing, viscosity, and self-diffusivity are often large. It is shown that much of the error for mixtures with water likely is due to neglect of the water polarizability, which results in too strong of an attraction between water and the [EtSO4] anion. 1. Introduction Ionic liquidsssalts that are molten close to ambient conditionsshave come under intense study recently, because of their unique properties. Many ionic liquids have negligible vapor pressure, high thermal stability, excellent solvation power, and a broad liquidus range. They have found use in many industrial applications, such as gas storage, acid scavenging, hydrosilyation, and electroplating.1 One particularly promising ionic liquid is 1-ethyl-3-methylimidazolium ethylsulfate ([emim][EtSO4]). It is easy to prepare in a halide-free manner from diethyl sulfate2 and, therefore, is comparatively inexpensive. It is water-miscible, air-stable, and has a relatively low viscosity. Preliminary studies suggest that it also has low toxicity.3 To use this or any other ionic liquid in an industrial application, a wide range of thermodynamic and transport properties must be known. A reasonable amount of experimental physical property data exists for [emim][EtSO4],4,5 although only under a limited range of conditions. Much effort is being devoted to the measurement and tabulation of the properties of ionic liquids,6 but the sheer number of possible compounds that can be made into ionic liquids means that not all candidates that could be effective in a given application will be examined. Moreover, it is often difficult to relate how properties are dependent on the underlying structural and chemical features of an ionic liquid from experiments alone, thus frustrating the search for new ionic liquids. Classical atomistic-level simulations have been shown to be an effective means for estimating the physical properties of ionic liquids.7 Simulations have been used to compute pure liquid properties such as liquid and crystalline density,8,9 liquid struc* To whom correspondence should be addressed. Tel.: (574) 6315687. Fax: (574) 631-8366. E-mail address:
[email protected]. † Present address: DuPont Central Research and Development, Experimental Station, Wilmington, DE 19880.
ture,10-13 heat capacity,14 enthalpy of vaporization,15-17 compressibility and volumetric expansivity,18,19 self-diffusivity 20-23 and viscosity.23-26 Several mixture properties have also been computed, although to a lesser degree.27-33 In most studies, analytic potential models (“force fields”) are proposed for a particular ionic liquid and then a handful of properties are computed. Validation of the force field is done by comparing the computed properties against available experimental data, most often liquid densities. Unfortunately, it has been shown that several different parameterizations of a force field can yield essentially the same liquid density, making this property not a particularly good choice for validation studies. The goal of the present study is to focus on one particular ionic liquid, [emim][EtSO4], and compute a wide range of thermodynamic and transport properties with a single force field. Because this ionic liquid has had many of its properties measured experimentally, an overall assessment of the accuracy of force field-based simulations for the prediction of several different ionic liquid properties can thus be obtained. To conduct this study, we use three different simulation techniques: equilibrium molecular dynamics (EMD), reverse nonequilibrium molecular dynamics (RNEMD), and continuous fractional component Monte Carlo (CFC MC). A new force field is also proposed for [emim][EtSO4]. For the neat ionic liquid, the following properties are computed: density, volumetric expansivity, heat capacity, enthalpy of vaporization, rotational relaxation time, self-diffusivity, shear viscosity, and thermal conductivity. In addition, the following properties are computed for water/[emim][EtSO4] mixtures: density, excess molar volume, enthalpy of mixing, partial molar enthalpy, water solubility as a function of partial pressure, rotational relaxation time, selfdiffusivity, shear viscosity, and thermal conductivity. To our knowledge, no other simulation study has been conducted in which so many different properties have been computed for a
10.1021/ie800843u CCC: $40.75 2008 American Chemical Society Published on Web 09/09/2008
9116 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008
Figure 1. Schematic representation of the [EtSO4] anion.
single ionic liquid with a single force field. Whenever possible, direct comparison with experimental data is made to assess the predictive power of the force field and methods. 2. Simulation Details 2.1. Force Field. A standard classical force field is used34 that accounts for bond stretching, angle bending, dihedral or torsion angle rotation, and improper angle bending. In addition, intramolecular and intermolecular van der Waals and Coulombic interactions are modeled. The functional form is
Figure 2. Comparison of the two torsional potentials, in terms of total energy of the configuration: (- - -) eq 3 (dashed curve) and (s) eq 4 (solid curve); (b) energies computed using MP2 and (O) energies computed using DFT.
Utotal ) Ubond + Uangle + Utorsion + Uimproper + Uvdw + Uqq (1) Utotal )
∑ k (r - r )
2
b
0
bonds
+
∑ k (θ - θ )
2
θ
0
+
angles
∑ ∑ k [1 + cos(nφ - δ)] + ∑ φ
dihedrals n
kψ(ψ - ψ0)2 +
impropers
N-1 N
∑∑ i)1 j>1
{ [( ) ( ) ] 4ij
σij rij
12
-
σij rij
6
+
}
qiqj (2) 40rij
The force field parameters for the [emim] cation were developed starting from the force field originally proposed by Cadena and Maginn.35 New nominal bond angles (θ0) and bond lengths (r0) were determined by performing quantum mechanical geometry optimization on the isolated cation at the B3LYP/ 6-311++G** level of theory using Gaussian 03.36 New partial atomic charges were derived using the CHELPG method.37 All other parameters were the same as those used by Cadena and Maginn.35 For the [EtSO4] anion, shown in Figure 1, nominal bond lengths, angles, and partial charges were derived by performing geometry optimization and CHELPG calculations on an isolated anion at the B3LYP/6-311++G** level of theory. Other parameters were taken from the work of Huige and Altona,38 who extended the CHARMM force field39 for sulfates and sulfamates. Note that, after the present work was completed, a force field for alkylsulfonate and alkylsulfate anions appeared in the literature.40 The [EtSO4] force field was tested by performing a short geometry optimization on an isolated anion and comparing the resulting conformation with that obtained from the quantum calculations. It was found that the force field and quantum calculations did not yield the same minimized structure. Specifically, the parametrization of the C1-C2-O3-S4 dihedral reported by Huige and Altona38 showed subtantial deviation from the quantum calculations. To rectify this, new parameters for this dihedral angle were determined in the following manner. The C1-C2-O3-S4 angle φ was scanned through 360°, in steps of 10°, and the total energy of each relaxed structure was computed at the B3LYP/6-311++G** level. As a check, the energies of these configurations were also computed using second-order Møller-Plesset (MP2) calculations. As shown below, there was very little difference between the B3LYP and MP2 calculation results. The resulting energy profile, EQM(φ), was compared against the energies obtained for the same configurations from the force field with all terms included,
except the dihedral energy for C1-C2-O3-S4. The resulting molecular mechanics energies, EMM(φ), were then compared against EQM(φ), using a common reference at φ ) 0. The difference between the two energy profiles was attributed to the dihedral energy (i.e., Utorsion(φC1-C2-O3-S4) ) EQM(φ) EMM(φ)) and was fit to two different functional forms, viz the standard OPLS potential with a zero term Utorsion )
∑
c0 + c1[1 + cos(φ)] + c2[1 - cos(2φ)]
dihedrals
+ c3[1 + cos(3φ)](3) and a three-term CHARMM dihedral potential Utorsion )
∑
k1,φ[1 + cos(φ - δ1)] +
dihedrals
k2,φ[1 + cos(2φ - δ2)] + k3,φ[1 + cos(3φ - δ3)](4) Figure 2 shows a comparison of the total molecular energy as a function of dihedral angle for the two force fields and the quantum calculations. The CHARMM (eq 4) and OPLS (eq 3) forms, fit to the MP2 calculations, yield almost-identical results. Therefore, either of the two models are recommended for the [EtSO4] anion; the parameters for these models are listed in the Supporting Information. Because of technical issues that are related to the different software used for the calculations presented in this work, the OPLS form of the dihedral potential was used for viscosity and thermal conductivity calculations, while the CHARMM form was used for all other calculations. The force-field parameters for water correspond to the widely used simple point charge (SPC) model,41 with flexibility added to the bonds and angles. The SPC water model has been shown to give reliable vapor-liquid equilibria results,42 and therefore was chosen for this work. Intramolecular nonbonded energies were neglected for atoms separated by fewer than three bonds and were scaled by 50% for atoms separated by exactly three bonds. For all other cases, full interactions were used. Lorentz-Berthelot combining rules were used to compute the cross-interactions. A complete list of force field parameters can be found in the Supporting Information. 2.2. Simulation Methods. Several different simulation methods, algorithms, and programs were utilized to enable the wide range of properties to be computed in this work. Specifically, some properties were computed with EMD simulations using
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9117 43
the code NAMD, whereas two different locally developed codes were used to calculate other properties with RNEMD and Monte Carlo. As such, there are some variations in the simulation details used for computing properties with the different methods. Whenever possible, however, validation tests were conducted to ensure that these variations in simulation details did not result in inconsistencies among the various computed properties. In other words, despite the differences in some of the details used in the various methods, the computed properties are expected to be internally consistent. For example, densities and volumetric expansivities were computed from isothermal-isobaric EMD simulations run at a pressure of 1 bar and various temperatures. Other properties were computed using simulations in the canonical (NVT) ensemble at the density determined from the isothermal-isobaric simulations. Rotational relaxation times were computed from molecular end-to-end vector autocorrelation functions during NVT EMD simulations. Self-diffusivity (Ds) values were also obtained from NVT EMD via application of the well-known Einstein relation to the mean square displacement: Ds ) lim tf∞
1〈 |r(t) - r(0)|2 〉 6t
(5)
Heat capacities were calculated using a combination of NVT EMD and ab initio calculations, as described in previous work.22 Viscosities and thermal conductivities were computed using the RNEMD method44-46 in the NVT ensemble. We have used this procedure previously to compute the viscosity of a series of polyhydric alcohols47 and an ionic liquid.25 Here, we utilize the RNEMD method for the calculation of ionic liquid thermal conductivities. The method works as follows. The simulation box is divided into N ) 14 small bins in one particular direction (e.g., the “y” direction). The method is most effective for elongated boxes, so, for this reason, the simulation boxes for all molecular dynamic runs were orthogonal, with an aspect ratio of 2.5. To compute viscosity, a stress or momentum flux in the “y” direction is imposed by exchanging the x-component of momentum between an atom having the most negative “x” momentum in bin n ) 1 with another atom in the central bin nc ) (N/2) + 1 that has the largest momentum in the positive “x”direction. The momenta swaps are conducted periodically, so that the total exchanged momentum is ptotal ) ∑(px,nc - px,n1)
(6)
Note that the frequency of the swaps determines the total imposed flux on the system. After any time t, the imposed momentum flux is given by jy(px) )
ptotal 2tLxLz
(7)
where Lx and Lz are the lengths of the simulation box in the “x”- and “z”-directions, respectively. The system responds to the nonphysical imposed momentum flux by establishing a real momentum flow in the opposite direction. This sets up a linear velocity profile in the upper and lower halves of the simulation box. From the slopes of the velocity profiles (which are easily determined) and the imposed momentum flux, the shear viscosity may be determined. To compute the thermal conductivity, the same basic algorithm is used, with the only difference being that, instead of atoms with maximum positive or negative “x”-directional momenta, atoms are selected with maximum and minimum kinetic energies. All three components of the velocity are swapped between molecules, resulting in the establishment of
a heat flux Jheat. This nonphysical heat flux generates a temperature gradient in the system,∇T, which is easily measured. The thermal conductivity λT is then determined via the application of Fourier’s law: Jheat ) -λT ∇ T
(8)
Water solubilities were computed using the continuous fractional component Monte Carlo (CFC MC) procedure in the osmotic ensemble. Details of this method have been given in recent publications.32,48 In this method, a series of Monte Carlo moves are constructed to satisfy the thermodynamic constraints of constant temperature, pressure, and solute fugacity that are imposed on the system. The resulting average liquid density and solute composition are computed. Configurations are changed using hybrid Monte Carlo to obtain thermal equilibrium. The volume of the box is changed to ensure that the imposed pressure is maintained. To satisfy the imposed water fugacity (partial pressure), water molecules are added to and removed from the cubic simulation box. Because the probability of inserting and deleting entire water molecules into this dense, highly charged liquid is extremely small, the CFC MC algorithm is used to perform the insertions and deletions gradually with an optimized bias potential. Approximately 60% of the moves were thermal equilibration moves, ∼38% were CFC insertion/ deletion moves, and ∼2% were volume change moves. The insertions were performed using a thermally equilibrated reservoir of 500 water molecules. The Lennard-Jones potential was multiplied by a switching function of the form
{
S(|b r ij|) ) (if |b r ij| e rs)
1 r ij| ) (rc + 2|b r ij| - 3rs ) (rc - |b 2
2 2
2
2
(rc2 - rs2)3 0
2
(if rs < |b r ij| < rc) (9) (if rc e |b r ij|)
during EMD and CFC MC simulations. The cutoff distance was rc ) 12 Å in all cases. The onset distance was rs ) 10 Å for the EMD simulations and rs ) 10.5 Å for the Monte Carlo calculations. In the RNEMD simulations, a cutoff of 12 Å was used without a switching function but long-range corrections were applied. Electrostatic interactions were computed in the MD simulations with the particle-mesh Ewald algorithm49 with real space cutoff distance of 12 Å. A previously validated damped shifted force procedure32,50 was used in the Monte Carlo simulations. A reversible multiple time-stepping algorithm (r-RESPA)51 was used in the MD simulations with inner time step (δ) of 0.5 fs and an outer time step (τ) of 2.0 fs. A four-chain Nose´-Hoover thermostat52 was used to control the temperature, whereas a Langevin damping method was used to control pressure in the isothermal-isobaric simulations. 2.3. State Points and System Size. For the pure ionic liquid simulations, 200 ion pairs were studied at T ) 308, 328, 348, and 368 K. To examine the effect of water content on transport properties, three different water mole fractions (xwater ) 0.26, 0.5, and 0.75) were studied at 348 K. To understand the effect of water content on the equilibrium properties such as density and heat of mixing, nine additional compositions (xwater ) 0.8-0.96) were studied at 348 K. A pure water simulation was also conducted with 2000 SPC water molecules at 348 K. Water solubility calculations were conducted using 160 ion pairs at temperatures of 313, 348, and 473 K and pressures in the range of 0.0001-0.5 bar.
9118 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 Table 2. Results of Solubility Simulations for Water in [emim][EtSO4]a
Table 1. Computed Liquid Densities at Different Temperatures, Compared with Experimental Dataa Density (g/cm3) temperature, T (K)
Fexp
Fsim
308 328 348 368
1.2306 1.2172 1.2040 N/Ac
1.2175 1.2036 1.1896 1.1746
b
deviation (%) -1.1 -1.2 -1.2 N/Ac
a Data taken from ref 4. b Subscripts indicate the statistical uncertainties in the final digit. c The experimental density was not available at 368 K.
3. Results 3.1. Thermodynamic Properties. 3.1.1. Pure Liquid Density. EMD simulations were conducted for 2 ns in the isothermal-isobaric ensemble to compute liquid densities. The results are compared against experimental values in Table 1. Computed densities are systematically lower than the experimental values, but the agreement is within 1.2% at all temperatures. It has been noted in previous simulation studies22 that ionic liquids have a tendency to exhibit sluggish dynamics, which can pose a problem for EMD simulations. To test for nonergodic behavior, up to three independent simulations were performed to test whether the computed density was repeatable. In all cases, each replicate was consistent with the others, giving an indication that the simulations were performed long enough to obtain a reliable density. 3.1.2. Volumetric Expansivity. The volumetric expansivity was computed from the temperature dependence of the molar volumes via the following relation: RP )
( )
1 ∂〈V〉 〈V〉 ∂T
P
(10)
The experimental value4 for this quantity is 5.37 × 104 K-1, whereas the simulations yield a value of 6.26 × 104 K-1, which is 16.5% higher than the experimental value. This level of agreement is typical for ionic liquid simulations; derivative properties such as volumetric expansivity are typically less accurate than single-point properties. 3.1.3. Heat Capacity. The constant pressure heat capacity (CP) was computed for [emim][EtSO4] at 348 K by following the procedure that was described by Cadena et al.22 This procedure involves determining the residual contribution to the heat capacity from the temperature dependence of the intermolecular portion of the enthalpy in classical simulations. This residual contribution is added to the ideal-gas portion of the heat capacity obtained from a frequency analysis of the quantum mechanically optimized cation and anion conformations. The ideal-gas contribution was computed at 298 K and 1 bar and scaled to the proper temperature. An additional scale factor of 0.97 was used to account for the basis set.53 The estimated heat capacity at 348 K was 454 J/mol K, which is ∼12% higher than the experimental quantity reported by Ficke et al.5 3.1.4. Enthalpy of Vaporization. The enthalpy of vaporization (∆Hvap) is given as ∆Hvap(T, P) ) 〈Ugas(T) 〉 - 〈Uliq(T, P) 〉 + RT - PVliq
(11)
A procedure for how to compute the enthalpy of vaporization of ion pairs has been presented elsewhere.16 The computed enthalpy of vaporization for [emim][EtSO4] is 178 kJ/mol. The experimental value for this property is not yet available, but the number compares well with the enthalpy of vaporization of other ionic liquids determined computationally15,16 and experi-
temperature, T (K)
pressure, P (bar)
xH2Oa
γ
313 313 313 313 348 348 348 348 348 473 473 473
0.0001 0.0005 0.001 0.002 0.001 0.0025 0.005 0.02 0.04 0.074 0.2 0.5
0.133 (0.001) 0.204 (0.003) 0.227 (0.001) 0.259 (0.002) 0.077 (0.002) 0.148 (0.004) 0.170 (0.002) 0.322 (0.003) 0.539 (0.002) 0.028 (0.003) 0.086 (0.002) 0.209 (0.002)
0.019 0.063 0.113 0.198 0.071 0.092 0.161 0.339 0.406 0.234 0.206 0.212
a Uncertainties are computed as standard deviations of block averages and are shown in parentheses.
mentally.15,54 As the temperature T increases, the enthalpy of vaporization ∆Hvap decreases at a rate of ∼0.1 kJ/(mol K). 3.1.5. Water Isotherms. The compound [emim][EtSO4] is known to be completely miscible with water at room temperature and, like most ionic liquids, has a tendency to absorb large amounts of water, even under low humidity conditions. Therefore, it is important to understand how the properties of the ionic liquid change as a function of water content, because it is unlikely to ever be used in a pure state. CFC MC simulations were run to compute water isotherms at 313, 348, and 473 K at various fugacities, with pressures computed from the Peng-Robinson equation of state. Given the low pressures, the water vapor phase was essentially ideal in all cases; therefore, the results are not sensitive to the choice of equation of state. Results are given in Table 2. The activity coefficient (γ) was computed as γ ) P/(xPsat), where x is the computed mole fraction of water in the ionic liquid and Psat is the saturation pressure for the SPC water model. The activity coefficients are all less than unity, indicating negative deviations from Raoult’s law, which is consistent with favorable interactions between water and the ionic liquid. Note that the saturation pressures were computed from separate Gibbs ensemble simulations,55 with values of 0.039 ( 0.002 bar at 313 K, 0.183 ( 0.006 bar at 348 K, and 11.3 ( 0.5 bar at 473 K. These saturation pressures differ from the experimental values by 50%-70%, which is a well-known artifact of the SPC model, and a manifestation of the difficulty in modeling the phase behavior of water with simple fixed-charge models. Although there are more than 46 distinct water models from which to chose,56 the SPC model was used in the present work, because it does a reasonable job matching many properties of water,is simple and used widely. Other water models could be used and comparison be made to the SPC model, but such a study was beyond the scope of the present work. To compare the computed isotherms with the experimental results of Heintz and co-workers,57 the solubility was plotted as a function of reduced pressure, where the computed (experimental) vapor pressure was used as the reducing parameter for the simulations (experiments). Figure 3 shows the results. The agreement between the simulations and experiments is quite good, although there are no experimental data in the low concentration (x < 0.2) regime. (This part of the computed isotherm has been expanded in the inset of Figure 3.) Note that, for both the simulations and experiments, the solubility values almost fall on a single curve. This is similar to the behavior observed for SO2 in an imidazolium- and pyridinium-based ionic liquid58 and for water in three different “hydrophobic” ionic liquids59 and indicates that, for these systems, differences in vapor pressure are the key determination of solubility.
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9119
Figure 3. Computed and experimental mole fraction versus reduced pressure P/Psat for H2O in [emim][EtSO4]. Computed values are obtained at (9) 313 K, (b) 348 K, and ([) 473 K. (Error bars are smaller than the symbol size.) Experimental values are shown at (4) 303 K, (0) 313 K, and (O) 323 K. Inset shows the computed isotherm at low pressure.
Henry’s law constants were estimated from the slope of a straight-line fit to the low-pressure simulation results. Values obtained were 0.002 ( 0.001 bar at 313 K, 0.017 ( 0.002 bar at 348 K, and 2.37 ( 0.04 bar at 473 K, respectively. From the van’t Hoff equation, the partial molar enthalpy at infinite dilution was estimated to be -9 ( 5 kJ/mol (using liquid water as a reference). The relatively large error bar is indicative of the fact that only three temperatures were examined. However, the result is consistent with the EMD mixture calculations described below. The high-temperature simulation results were extrapolated to 298 K and the Henry’s law constant at this temperature was estimated to be 0.0007 bar. We have already performed simulations for water absorption in 1-n-hexyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([hmim][Tf2N]) at 298 K.32 The computed Henry’s law constant for water in this ionic liquid at 298 K was 0.04 ( 0.01 bar, which demonstrates how much more soluble water is in [emim][EtSO4] than in [hmim][Tf2N] and that the simulations correctly capture this known trend. The difference in water solubility is mainly due to difference between the “hydrophillic” [EtSO4] anion and the “hydrophobic” [Tf2N] anion, but also because the longer alkyl chain on the [hmim] cation decreases water solubility. Recall that the SPC water model gives a saturation pressure that is less than the experimental value. Therefore, although the simulations and experiments yield the same solubility at the same reduced pressure, the model predicts higher water solubility than the experiment at the same total pressure. This suggests that the interaction between water and the ionic liquid is too strong. As shown below, this may be due to the fact that the SPC water model has an exaggerated dipole moment to match the properties of pure water. In the presence of a lowdielectric material such as an ionic liquid, the actual dipole moment of water should be lower, yielding weaker electrostatic interactions between water and the ionic liquid. Despite these problems with the SPC model, it is widely used by many groups and therefore we decided to test its predictive ability in modeling additional mixture thermodynamic and transport properties. 3.2. Densities as a Function of Water Content. The effect of water content on liquid density was examined across the entire composition range by performing isothermal-isobaric EMD simulations at a temperature of 348 K and pressure of 1 bar. At a common water loading, the densities obtained from the EMD simulations are consistent with the CFC MC calculations, but it is easier to span a given composition range with EMD. A comparison of the simulation results and experimental data is
Figure 4. Comparison of the experimental4 (open circles) and simulated (MD (filled symbols) and MC (crosses)) densities, as a function of water content (xwater) at 348 K; q1, q2, q3, and qSPC correspond to different sets of partial charges (see text). The dashed and dotted lines are guides to the eyes for the simulation and experimental results, respectively.
given in Figure 4. The agreement between simulation and experiment is within 1.5% over the entire composition range, which is remarkable, given that no adjustments were made to the force field parameters to match any experimental data (except the SPC model has been fit to pure water density). The error bars for the MD simulation results were estimated from the standard deviation of the density during the 2-ns production run. For the CFC MC calculations, the error bars were estimated from a standard error propagation, using fluctuations in the mixture molar volume and mole fraction during the production run. Similar to many other fixed-charge models for water, the SPC model uses partial charges that result in a dipole moment greater than what occurs in the gas phase. For example, the partial charge on oxygen is -0.82 (with the charge on each hydrogen being +0.41), whereas ab initio calculations on a single water molecule under vacuum result in charges of -0.72 and +0.36. These larger partial charges account for the fact that, in the high dielectric medium of liquid water, significant polarization occurs that would not normally be captured with a fixed-charge model. In fact, the computed density of liquid water is far too low if vacuum partial charges are used with the SPC water model. However, it is known that the dielectric constant of ionic liquids is small (typically ∼10).60,61 Based on this, it is possible that water will not be nearly as polarized when dissolved in an ionic liquid as it is when surrounded by other water molecules. To test whether polarizability plays a large role in mixture properties, simulations were performed at xwater ) 0.26, using three additional sets of partial charges:q1represents a “vacuum” condition, with an oxygen charge of -0.72; q2 represents charges obtained from ab initio calculations of water in a dielectric continuum similar to liquid CCl4, resulting in an oxygen charge of -0.77498; and q3 represents charges typical of water in a more polar solvent, with an oxygen charge of -0.80. For pure liquid water simulations, only the standard SPC charges were used, because these are optimized for the pure water environment. As shown in Figure 4, the various charges have little effect on the density of the mixture at xwater ) 0.26. This suggests that polarizability should have little influence on density, which is consistent with previous findings.23 E 3.2.1. Excess Molar Volume (Vm ). The excess molar E volume, Vm, is an extremely difficult property to obtain accurately from a simulation, because very small differences
9120 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008
Figure 5. Comparison of the experimental4 (open circles) and simulation (filled symbols) excess molar volumes, as a function of water content (xwater) at 348 K. The dashed and dotted lines are guides to the eyes for the SPC and experimental results, respectively. (See text for an explanation of the different charge sets q.)
in density can result in large differences in excess molar volume. As such, it is a stringent test of a force field. The excess molar volume was computed as VEm ) Vm-[xILVIL - xwaterVwater]
(12)
where xi is the mole fraction of component i in the mixture and VIL and Vwater are the pure component molar volumes of the ionic liquid (IL) and water, respectively. Vm is the molar volume of the mixture. The quantity given in the square brackets is the E molar volume of an ideal mixture. One also can compute Vm from the density as VEm )
∑ x M [ F1 i
i
i
m
-
1 Fi
]
(13)
where Fi is the density of the pure component, Fm the density of the mixture, and Mi the molecular weight of the pure component. Both of these formulas gave the same results. The computed excess molar volumes are compared with the experimental results of Rodriguez and Brennecke4 in Figure 5. While both the simulations and experiments indicate that the excess molar volume has a tendency to be small and slightly negative, the simulations using the SPC model overestimate the magnitude of excess molar volumes by ∼1 cm3/mol. This is a relatively small absolute difference, however, which is indicative of the sensitivity of the excess molar volume to the shape of the density/composition curve. The simulations using the “least polarized” model for water, charge set q1, agree best with the experimental data. As the magnitude of the partial charges decreases, the excess molar volume becomes less negative. Given the magnitude of the error bars, however, there is essentially no difference between the different charge sets. Note that the experimental data of Rodriguez and Brennecke4 are in good agreement with a second set of experimental data from Lu et al.,62 and that a slightly positive excess molar volume is observed at very high water concentrations. There was no E indication of the reversal in sign of Vm in the simulations with the SPC water model. 3.2.2. Enthalpy of Mixing. The enthalpy of mixing (∆hmix) was computed as ∆hmix ) hm - [xILhIL + xwaterhwater]
(14)
where hm is the molar enthalpy of the system, based on the total number of moles in the mixture, whereas hIL and hwater are
Figure 6. Comparison of the experimental5 (open circles) and simulated (filled symbols) enthalpies of mixing, as a function of water content (xwater) at 348 K. The dashed lines are polynomial fits to the results, and the dotted lines are van Laar model fits.
the molar enthalpies of the pure ionic liquid and water, respectively. The quantity in the square brackets denotes the ideal enthalpy of mixing. The enthalpies of mixing computed from the simulations are compared with the experimental data of Ficke et al.5 in Figure 6. The experimental and simulation results were fit to a van Laar model (denoted by the dotted lines in Figure 6) with the functional form ∆hmix ) x(1 - x)[A1x + A2(1 - x)]
(15)
as well as to a fourth-order polynomial (denoted by the dashed line in Figure 6). While both the experiments and simulations indicate a slightly negative enthalpy of mixing, the simulations using the SPC model consistently yield an enthalpy of mixing that is 1-7 kJ/mol greater in magnitude than the experimental data. In absolute terms, this is a relatively small energy difference (certainly within the range of what can reasonably be expected of a simple classical force field), but it is curious that the computed energies were systematically lower than the experimental values. The simulations that were performed at xwater ) 0.26 and using different charge sets indicate that much of the discrepancy between experiments and simulations is due to a neglect of polarizability. The enthalpy of mixing becomes less negative as the magnitude of the charges on water decreases from that used with the SPC model to that of a vacuum (q1). This is consistent with the fact that the dielectric constant of water is much larger than that of the ionic liquid, and it suggests that the inclusion of polarizability would reduce the enthalpy of mixing. Note that the modified charges were only used in computing hm; for hwater, the normal SPC water charges were used, because this parameter set is most appropriate for pure water simulations. These results indicate that, although densities are insensitive to the inclusion of polarizability, mixture enthalpies for fluids that have very different dielectric constants (such as water and ionic liquids) are sensitive to polarization, and this should be included in models if a high degree of accuracy is required. We would anticipate that the error associated with neglecting polarizability will be smaller for solutes whose dielectric constants are smaller than water, but this remains to be verified. 3.2.3. Partial Molar Enthalpy. The partial molar enthalpy of water (h˜water) is defined as ˜ hwater ) ∆hmix - xIL
∂∆hmix + hwater ∂xIL
(16)
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9121 Table 3. Temperature Dependence of the Rotational Relaxation Times for the [emim] Cation and [EtSO4] Anion in the Pure Liquid Rotational Relaxation Time (ps) τcation
temperature, T (K)
τanion
a
308 328 348 368
2321a 454 263 149
1777 623 361 216
a Value is approximate, because of the fact that these times are on the same order as the simulation time.
Table 4. Rotational Time Constants for the [emim] Cation, the [EtSO4] Anion, and Water, as a Function of Water Content at 348 K Rotational Relaxation Time (ps) water content, xwater
τcation
τanion
0.00 0.26 0.50 0.75 1.00
361 352 252 152
263 300 272a 140
Figure 7. Comprison of experimental and calculated partial molar enthalpy of water in [emim][EtSO4] at 348 K.
where hwater is the pure component molar enthalpy of liquid water. All the terms in the aforementioned equation are known from the simulations except ∂∆hmix/∂xIL. The partial derivatives were obtained analytically from the van Laar and polynomial fits of the data. The infinite dilution partial molar enthalpy of water can be compared with that obtained from a van’t Hoff plot of the Henry’s law constants. As can be seen from Figure 7, however, the results are quite sensitive to the fitting procedure at low water concentrations, with the simulations showing the largest deviation between the two fitting procedures. The fact that there are so few simulation points in the low water concentration regime makes estimates of the infinite dilution partial molar enthalpy extremely imprecise. Additional simulations could be run at lower water compositions, but with a fixed number of IL ion pairs (set at 200), the accuracy of the calculations becomes questionable when only a few water molecules are modeled. Alternatively, a much larger system could be modeled such that a representative sampling of water molecules is contained in the system. However, as the system size becomes large, the computational effort required to perform the simulation becomes unacceptably large. Consistent with the enthalpy of mixing results, the SPC water model yields partial molar enthalpies that are larger in magnitude (more negative) than the experiments. The simulations predict that the infinte dilution partial molar enthalpy ranges from -8 kJ/mol to -17 kJ/mol, which agrees with the result obtained from the Monte Carlo solubility simulations of -9 ( 5 kJ/mol. On the other hand, the experimental results vary from -1.8 to -3 kJ/mol, depending on the fitting procedure. Again, we anticipate that the use of a polarizable water model would result in much better agreement with the experiment. 3.3. Dynamic Properties. In addition to thermodynamic properties of the pure ionic liquid and its mixtures with water, dynamical properties may also be computed from force-fieldbasedsimulations.Startingfromtheequilibratedisothermal-isobaric simulation results, additional equilibration runs of 5 ns were conducted in the canonical ensemble at 308, 328, 348, and 368 K. These equilibrated conformations were then used as initial configurations for the RNEMD simulations, from which viscosities and thermal conductivities were computed. Other dynamic properties such as rotational relaxation times and self-diffusivities were computed from the 5-ns canonical ensemble simulations. 3.3.1. Rotational Relaxation Times. The rotational relaxation times of the cation and anion were computed as a function
a
τwater 95 81 37 1
This estimate of rotational relaxation is inaccurate.
of temperature for the pure ionic liquid. These time constants were determined by first calculating the correlation function C(t): C(t) ) (e1(0) · e1(t))
(17)
where e1 is the unit vector of the longest axis of the species in question. At the lowest temperatures, the correlation function does not decay to zero over the length of a single simulation, so estimates of the time constants were obtained by first fitting the correlation function to two stretched exponentials of the following form:
[ ( )]
C(t) ) a1 exp -
t τ0
β1
[ ( )]
+ a2 exp -
t τ0
β2
(18)
The first term was fit to the data from 0-0.5 ns, whereas the second was fit to the remainder of the data. The rotational time constant was then found from a numerical integration of C(t): τ)
∫
0.5
0
[ ( )]
a1 exp -
t τ0
β1
dt +
∫
∞
0.5
[( )]
a2 exp -
t τ0
β2
dt (19)
Table 3 summarizes the estimated rotational relaxation times computed for the pure ionic liquid ions at different temperatures. The relaxation times increase dramatically as the temperature decreases, which is why very long simulations are required to properly compute properties at low temperature. Note that, at the lowest temperature, the time constants are of the same order of magnitude as the actual MD simulations; these time constants should be taken as only estimates. Generally, the anion rotates slightly faster than the cation, but the differences are not as large as what is observed for systems with small, spherical anions.8 Rotational relaxation times for the ions and water were calculated at 348 K and water mole fractions of xwater ) 0, 0.26, 0.50, 0.75, and 1.0. The mobility of the ions and water was great enough that direct numerical integration of C(t) could be used to determine τ. Table 4 summarizes the findings. The rotational time constant for the cation decreases as water content increases, and the time constant for water increases as the ionic liquid content increases. This is expected behavior, because highly mobile water has a tendency to accelerate the dynamics of the cation and the ionic liquid has a tendency to inhibit the motion of the water. The results shown here are consistent with the recent computational findings of Annapureddy and co-
9122 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 Table 5. Simulation Results for Self-Diffusivities of the [emim] Cation and the [EtSO4] Anion at Different Temperatures Self-diffusivity (× 1012 m2/s) temperature, T (K) 308 328 348 368 a
Figure 8. Snapshot from a bulk simulation showing water-forming hydrogen-bonding “bridges” between two [EtSO4] anions. This behavior may be responsible for the slowing of anion rotational motion at low water concentrations.
workers,33 who observed increased rotational dynamics for the 1-hexyl-3-methylimidazolium cation as water was added to the system. Interestingly, however, the trend is different for the [EtSO4] anion. At low water concentrations, the rotational motion of the anion actually slows somewhat, from 263 ps for the neat liquid up to 300 ps at xwater ) 0.26. At xwater ) 0.5, the anion correlation function showed unexplained fluctuations for up to 5 ns, making the integration subject to large error. The actual time constant likely is larger than the value reported here. At xwater ) 0.75, the rotational time constant of the anion decreases again, to 140 ps. One possible reason for this behavior is that water preferentially associates with the [EtSO4] anions in the simulations, forming hydrogen-bonded “bridges” between two [EtSO4] anions, as shown in Figure 8. These hydrogenbond bridges have a tendency to strengthen the interaction between two anions, effectively linking the anions together and hindering their rotational motion. Above xwater ) 0.5, there is more than one water molecule per anion, so the ions begin to be fully solvated and the rotational dynamics of the anion increases again. 3.3.2. Self-Diffusivity. Because of the sluggish dynamics of ionic liquids, it is important to ensure that the simulations are long enough to probe diffusive behavior. For the present simulations, 5-ns production runs were used. These simulations are longer than any ion rotational relaxation time, although at 308 K, the rotational and simulation times approach this time. One way to test if the simulations are long enough to probe diffusive behavior is to compute a linearity parameter, β(t): d log(∆r2(t)) (20) d log(t) where ∆r2(t) is the mean square displacement at time t. By definition, β(t) must approach unity in the diffusive regime. At values below β(t) ) 1, the system is in the subdiffusive regime and any apparent self-diffusivity will underestimate the true selfdiffusivity. Table 5 lists the computed self-diffusivities for the pure [emim][EtSO4], whereas Figure 9 compares these results with experimental pulsed field gradient nuclear magnetic resonance (PFG NMR) values obtained by Vasenkov and co-workers.63 Even after 5 ns, the simulations at 308 K had not reached the
D-
D+ a
4.9 15.4 21.4 39.3
1.9a 7.2 9.5a 20.0
This value has a β(t) value of
