Article pubs.acs.org/JPCB

Saturation Properties of 1‑Alkyl-3-methylimidazolium Based Ionic Liquids Kaustubh S. Rane and Jeffrey R. Errington* Department of Chemical and Biological Engineering, University at Buffalo, The State University of New York, Buffalo, New York 14260-4200, United States S Supporting Information *

ABSTRACT: We study the liquid−vapor saturation properties of room temperature ionic liquids (RTILs) belonging to the homologous series 1alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([Cnmim][NTf2]) using Monte Carlo simulation. We examine the effect of temperature and cation alkyl chain length n on the saturated densities, vapor pressures, and enthalpies of vaporization. These properties are explicitly calculated for temperatures spanning from 280 to 1000 K for RTILs with n = 2, 4, 6, 8, 10, and 12. We also explore how the identity of the anion influences saturation properties. Specifically, we compare results for [C4mim][NTf2] with those for 1-butyl-3-methylimidazolium tetrafluoroborate ([C4mim][BF4]) and 1-butyl-3-methylimidazolium hexafluorophosphate ([C4mim][PF6]). Simulations are completed with a recently developed realistic united-atom force field. A combination of direct grand canonical and isothermal−isobaric temperature expanded ensemble simulations are used to construct phase diagrams. Our results are compared with experimental data and Gibbs ensemble simulation data. Overall, we find good agreement between our results and those measured experimentally. We find that the vapor pressures and enthalpies of vaporization show a strong dependence on the size of the alkyl chain at low temperatures, whereas no particular trend is observed at high temperatures. Finally, we also discuss the effect of temperature on liquid phase nanodomains observed in RTILs with large hydrophobic groups. We do not observe a drastic change in liquid phase structure upon variation of the temperature, which suggests there is not a sharp phase transition between a nanostructured and homogeneous liquid, as has been suggested in earlier studies.

1. INTRODUCTION Room temperature ionic liquids (RTILs) are being pursued as potential solvents in numerous chemical processes. A good understanding of liquid−vapor equilibria (VLE) is necessary to transfer these technologies from the laboratory scale to the production scale. RTILs are considered “green solvents” because of their low volatility, which limits their spread in the atmosphere. However, recent studies have shown that they have appreciable toxicity.1 Therefore, knowledge of saturation properties can help in quantifying the volatility of RTILs and understanding their environmental impact. Moreover, VLE studies can play an important role in developing processes to recover these fluids from waste streams. Finally, studying the bulk liquid−vapor phase behavior of RTILs represents an important first step toward understanding the interfacial behavior of these compounds, which is important to many technologies. Molecular simulations have played an important role in studying the thermophysical properties of RTILs.2 Given the highly tunable nature of these fluids, simulations can aid in their “rational design” for specific applications. The liquid−vapor saturation properties of RTILs have recently received considerable attention. Rai and Maginn used Gibbs ensemble © XXXX American Chemical Society

Monte Carlo (GEMC) simulation to compute saturated densities and vapor pressures of RTILs at temperatures spanning from 850 K to the critical region.3,4 The lower temperature limit is related to sampling difficulties generally encountered in GEMC simulations. Saturation pressures have also been measured experimentally.5,6 However, the temperature range considered within these studies spans from approximately 400 to 600 K, with the upper limit related to the thermal instability of RTILs at high temperatures. Therefore, in order to compare the vapor pressures obtained from molecular simulations with those from experiments, one requires a simulation strategy to study VLE at low temperatures. Although RTILs break down at temperatures well below the critical point, there is significant interest in understanding how critical properties evolve with chemical identity. Knowledge of these parameters is important for various corresponding state theories and equations of state that are employed to predict the bulk and interfacial properties of RTILs.7−9 Rebelo et al. used Received: April 26, 2014 Revised: June 24, 2014

A

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Several studies have reported the formation of polar and nonpolar “nano-domains” within the liquid phase of RTILs with sufficiently long alkyl side chains.21−23 At sufficiently low temperature, the alkyl chains and ionic moieties organize into distinct nanometer-scale domains, with the regions forming a semiregular three-dimensional network. At relatively high temperature, this local ordering breaks down. Wang and Voth identified these domains within configurations generated from molecular dynamics simulations involving a coarse-grained model of an imidazolium based ionic liquid containing relatively long side chains.21 These authors later used similar models to study the temperature dependence of the inhomogeneous liquid phase.24 They noted that the temperature evolution of the specific heat and a structural order metric that captures the extent of nanodomain formation exhibited signatures consistent with a structural phase transition between a nanodomain-dominated structure at low temperature and a uniform structure at high temperature. Moreover, room temperature molecular dynamics simulations performed by Lopes and Pádua with all-atom models of imidazolium based ionic liquids containing long alkyl chains also supported the presence of these nanometer-scale domains.22 Triolo et al. provided the first experimental evidence for these structural heterogeneities in RTILs on the basis of X-ray scattering studies.23 The presence of these nanodomains can significantly influence the potential applications of RTILs as solvents and as a media for self-assembly and photochemical reactions.25,26 However, except for the earliest studies completed using coarse-grained models, very little information is available on the temperature dependence of the liquid structure. In this work, we use isothermal−isobaric ensemble (NPT) simulations to probe the temperature dependence of the liquid structure of [C8mim][NTf2]. The paper is organized as follows. Within section 2, we describe the force field employed within our calculations and specify details related to the simulations. Our results are presented in section 3. Here, we compare simulation data for saturation properties with those from experiment and discuss the structure of the liquid phase. Within section 4, we provide a few concluding remarks and comment on future directions for this work.

the Guggenheim and Eötvos equations to estimate the critical properties of different RTILs from experimentally determined surface tensions and liquid densities.10 For the [Cnmim][NTf2] family, their results indicate that the critical temperature decreases with an increase in the length of the alkyl chain attached to the cation. Similar trends were observed by Rai and Maginn, who obtained critical properties from Gibbs ensemble simulations.3,4 Additionally, they observed that the magnitude of change in critical temperature per methylene group decreases with an increase in the length of the alkyl chain. They reasoned that this trend may be the result of a saturation of electrostatic screening associated with the alkyl chains. In this work, we use the Monte Carlo simulation based approach introduced by Rane et al.11,12 to compute saturation properties at temperatures ranging from room temperature to those above the hypothetical boiling point of a given ionic liquid. Within this approach, one first completes a series of direct grand canonical simulations at relatively high temperatures and subsequently employs a temperature expanded ensemble scheme to trace saturation curves to relatively low temperature. This strategy helps us in overcoming difficulties often encountered at low temperatures in sampling inhomogeneous states characterized by densities intermediate between those of the vapor and liquid phases. The method has been successfully applied to study liquid−vapor phase equilibria of various complex molecules, including nonionic molecules12 and the ionic liquid 1,3-dimethylimidazolium tetrafluoroborate ([C1mim][BF4]).11 Here, we study six RTILs belonging to the homologous series 1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)amide ([Cnmim][NTf2]). We also provide results for two additional RTILs, namely, 1-butyl-3methylimidazolium tetrafluoroborate ([C4mim][BF4]) and 1butyl-3-methylimidazolium hexafluorophosphate ([C4mim][PF6]), in an effort to understand how anion identity influences the phase behavior. We compare our results with those from earlier simulation3 and experimental studies.5,6,13−17 The enthalpy of vaporization ΔHvap features prominently in many experimental and computational studies. Knowledge of ΔHvap is useful in understanding the influence of interactions on the vapor and liquid phase properties. For example, in 1alkyl-3-methylimidazolium based ionic liquids, the dependence of ΔHvap on the alkyl chain length shows how the competition between electrostatic and dispersion interactions influences coexistence behavior.3 Additionally, ΔHvap is also used to calculate the Hildebrand solubility parameter, which is useful for predicting the solubility of solutes in RTILs.18 Moreover, ΔHvap is often used to validate force fields employed within molecular simulations.19,20 Experimental data for ΔHvap from different sources show significant discrepancies, as was highlighted recently by Verevkin et al.13 They attribute most of these deviations to the values of specific heats used to calculate ΔHvap at standard conditions from measurements performed at higher temperatures. The alkyl chain length dependence of ΔHvap for the homologous series [Cnmim][NTf2] has also proven to be sensitive to the temperature range examined. The relatively high temperature Gibbs ensemble results from Rai and Maginn3 suggest that ΔHvap generally decreases with increasing n, whereas results obtained from experiments and molecular dynamics (MD) simulations at relatively low temperature indicate the opposite. In the current work, we compute ΔHvap over a wide range of temperatures for each of the RTILs examined here and further probe these issues.

2. FORCE FIELD AND SIMULATION DETAILS Molecular Models. In this work, we use a united-atom model developed by Zhong et al.19 for imidazolium based ionic liquids. In this model, the imidazolium ring is treated explicitly, while all other groups like CH2, CH3, CF3, etc., are treated as united atoms. The functional form of this force field is similar to that of AMBER/OPLS. Most of the van der Waals (vdW) parameters are taken from the general AMBER force field (GAFF).27 The partial charges on different groups stem from a fitting procedure involving the electrostatic potential generated from quantum calculations performed on ion pairs. The bond and angle parameters are fitted to reproduce the bond length and angle distributions obtained from all-atom simulations. Some of the dihedral angle parameters are fitted to reproduce the energy distributions obtained from quantum calculations, while others are taken from the TraPPE united-atom model with some refitting. This force field has been shown to satisfactorily reproduce liquid densities measured from experiments. Moreover, the enthalpies of vaporization computed for the homologous series [Cnmim][NTf2] using molecular dynamics (MD) simulations show good agreement with B

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

experimental results. Additionally, this force field was previously used by Rai and Maginn3,4 to complete Gibbs ensemble liquid− vapor equilibrium calculations for two homologous series of ionic liquids. We also used this force field to perform a series of benchmark calculations featuring [C1mim][BF4].11 Methods. We use a combination of direct grand canonical (GC) Monte Carlo simulations and isothermal−isobaric temperature expanded ensemble (NPT-TEE) simulations to trace liquid−vapor saturation curves over the desired range of temperatures. A detailed explanation of these methods is provided in our earlier work, where we examined the performance of these techniques in determining the liquid− vapor saturation properties of the restricted primitive model (RPM) and the RTIL [C1mim][BF4].11 Here, we provide a brief overview of our approach. We first perform GC simulations at two relatively high temperatures and compute saturation properties at these temperatures with the assistance of histogram reweighting. We use the vapor pressures at these two temperatures to generate a guess for the pressure− temperature relationship along the saturation line. This relationship is then used within independent NPT-TEE simulations focused on the liquid and vapor phases. These simulations allow us to track the Gibbs free energy of the two phases along the proposed (guess) saturation curve. We then use a histogram reweighting based scheme to develop a revised estimate for the pressure−temperature relationship along the saturation line, which serves as input for a subsequent set of NPT-TEE simulations. The above process is repeated until the vapor pressure curve satisfies a convergence criterion. Generally, three iterations are sufficient to achieve convergence for the RTILs studied here. We use the Binder cumulant approach28 to locate critical points. The technique involves analysis of the temperature and system size dependence of the fourth order cumulant UL = 1 − ⟨ρ4⟩/(3⟨ρ2⟩2), where ρ is the fluid density. This cumulant is calculated from analysis of a density probability distribution that spans the vapor and liquid domains. One evaluates UL at liquid−vapor saturation conditions over a range of near-critical temperatures for select system sizes. At the critical temperature, UL is nearly independent of system size. Therefore, the critical temperature corresponds to the point at which the Binder cumulant curves corresponding to two or more system sizes intersect. Examples of UL curves for [C2mim][NTf2] at three system sizes are provided in Figure 1. From a theoretical perspective, these system-size-specific curves are expected to intersect at a single temperature only if the corrections to finitesize scaling are negligible.28 However, as can be seen in Figure 1, these distinct intersection points are relatively close to each other and their average provides a satisfactory estimate of the critical temperature. The molecule number probability distribution required for this analysis is obtained by performing a direct grand canonical simulation12 at an approximate value of the critical temperature, which we obtain from a law of rectilinear diameters analysis.29 For each molecule number sampled, we also collect the configurational energy probability distribution. Histogram reweighting30,31 is then used to obtain probability distributions at nearby temperatures and to establish coexistence conditions. Here, we consider the maximum in compressibility as an indicator of liquid−vapor coexistence. The critical density is obtained from the average number of molecules present in the largest simulation box at the critical condition. Finally, we note that the calculation of UL requires a probability distribution that

Figure 1. Variation of UL with temperature for [C2mim][NTf2] for three system sizes. Solid black, red, and green lines represent data obtained using simulation box volumes of 1 × 105, 1.5 × 105, and 2.0 × 105 Å3, respectively. The dashed circle denotes the approximate location of the critical point.

is sufficiently smooth, such that the free energy barrier between liquid and vapor domains is clearly identified. This condition was not met for n = 12 because of sampling challenges associated with the long alkyl chain. Therefore, we do not present critical parameters for n = 12. Growth expanded ensemble and distance-bias insertion and deletion techniques are used to improve the efficiency of the molecule addition and removal process within GC simulations. The sampling of microstates in GC and NPT-TEE simulations is facilitated using reservoir based moves,32 hybrid moves,33 and aggregation-volume-biased Monte Carlo moves.34,35 The probability distributions (e.g., density or temperature probability distributions) are calculated using a transition matrix Monte Carlo (TMMC) based scheme.12 Simulation Details. Direct GC simulations are performed at temperatures of 900 and 1000 K for all RTILs within the [Cnmim][NTf2] series and for [C4mim][PF6]. For [C4mim][BF4], we perform GC simulations at 850 and 1000 K. GC simulations are performed using a cubic simulation cell of volume 1.0 × 105 Å3 for smaller ionic liquids (n < 6) and of volume 2.0 × 105 Å3 for larger ones. The MC move mix includes ion displacements and rotations, ion-pair displacements, HMC, HMC-One, and molecular insertions (or growth) and deletions (or reductions). The details of these moves are provided in our earlier work.11 Insertion/deletion moves are selected with 60% probability, and the displacement, rotation, HMC, and HMC-One moves are partitioned such that the computational time invested in each move is roughly equal. The Lennard-Jones interactions are truncated using a spherical cutoff distance of 14 Å, and long-range interactions are accounted for with a mean-field correction. The Ewald summation with tinfoil boundary conditions is used to calculate electrostatic interactions.36 NPT-TEE simulations span the temperature range defined by Tmin ≈ 280 K to Tmax ≈ 1000 K. Within the expanded ensemble scheme, subensembles are arranged with a constant shift in the inverse temperature β = 1/kBT of Δβ = 2.0 × 1017 J−1. NPT-TEE simulations are performed with 100 ion pairs. The MC move mix includes ion displacements and rotations, C

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

increasing n over the entire temperature range. Our results for n = 4 and 6 are in good agreement with the high-temperature Gibbs ensemble results of Rai and Maginn.3 However, for n = 2, our density values are low relative to their results. Further studies are necessary to understand the reason for this difference. We also provide liquid densities at low temperature from the molecular dynamics study of Zhong et al.19 These results were obtained from isothermal−isobaric simulations at a pressure of 1 bar. Since the vapor pressures at these temperatures are very low (refer to Figure 4), we expect that the liquid densities are negligibly affected by this difference in pressure, and therefore, we compare our results with them. In all cases, there is very good agreement between our results and those of Zhong et al. Figure 2 also shows experimental density values.14−17 These measurements were performed near room temperature at a pressure of 1 bar. For n = 2 and 8, we show data from Gardas et al.,15 whereas results from Fredlake et al.,16 Azevedo et al.,17 and Tomé et al.14 are provided for n = 4, 6, and 10, respectively. Our simulation results are in good agreement with those from experimental studies. Figure 3 provides saturated vapor densities as a function of temperature. Unlike the liquid densities, the trend with respect

ion-pair displacements, aggregation-volume-bias ion and ionpair displacements, HMC and HMC-One moves, volume change attempts, and subensemble change moves. Volume changes are selected with a probability of approximately 1/N, subensemble change moves are attempted with a probability of 0.01, and the remaining moves are distributed using the “equal computational cost” approach noted above. The spherical cutoff distance for the Lennard-Jones and realspace Ewald contribution is taken to be a constant fraction of the simulation cell length, usually 50%. Critical properties are obtained from GC simulations performed at 1200, 1210, 1190, 1180, and 1170 K for n = 2, 4, 6, 8, and 10, respectively. For each component, we perform simulations with three system sizes. The volumes used are 1.0 × 105, 1.5 × 105, and 2.0 × 105 Å3 for n = 2, 4, 6, and 8. For n = 10, we use simulation cells of volume 1.5 × 105, 2.0 × 105, and 2.5 × 105 Å3. Isothermal−isobaric ensemble simulations are used to probe the structure of liquid phase nanodomains within [C8mim][NTf2] at 36 temperatures. After the system is equilibrated, configuration snapshots are collected every 100 MC cycles, and these configurations are used to construct radial distribution functions of interest. Statistical uncertainties for all properties are taken as the standard deviation of results from four independent simulation runs.

3. RESULTS AND DISCUSSION Saturated Densities. Figure 2 contains saturated liquid densities for RTILs belonging to the homologous series [Cnmim][NTf2]. We observe that densities decrease with

Figure 3. Saturated vapor densities of RTILs belonging to the homologous series [Cnmim][NTf2] as a function of temperature. Symbol and color schemes are the same as those used in Figure 2. Note that the densities are plotted using a logarithmic scale. Insets a and b show magnified profiles at relatively high and low temperatures, respectively. Figure 2. Saturated liquid densities of RTILs belonging to the homologous series [Cnmim][NTf2] as a function of temperature. Black, red, green, blue, violet, and cyan solid lines represent NPT-TEE simulation results for n = 2, 4, 6, 8, 10, and 12, respectively. The NPTTEE results consist of approximately 890 data points. Uncertainties are provided at select temperatures. The black, red, and green squares denote GEMC simulation results from Rai and Maginn.3 Black, red, green, blue, and violet pluses denote results from molecular dynamics simulations of Zhong et al.19 Black and blue triangles represent the experimental data of Gardas et al.15 Red circles, green down triangles, and violet squares correspond to experimental data from Fredlake et al.,16 Azevedo et al.,17 and Tomé et al.,14 respectively. The color scheme used for symbols matches that for solid lines.

to n is not consistent over the entire temperature range. Within inset a of Figure 3, we have magnified the plot to show data at the highest temperatures studied. We note that RTILs generally degrade at temperatures above 500−600 K.37 Therefore, the properties calculated at these temperatures are hypothetical in nature. However, such high-temperature saturation properties are useful in estimating critical properties, which help in developing equations of state and corresponding state theories.7 Here, we do not observe a discernible trend in saturated vapor densities. Our values are systematically lower than those reported by Rai and Maginn.3 Similar differences between our D

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

noted that, for RTILs with n < 4, the force field employed here shows an n dependence for enthalpies of vaporization that is stronger than that observed in experiments.19 For RTILs with n > 1, the force field uses the same set of partial charges for the constituents of the alkyl chain. Perhaps this charge distribution requires modification for RTILs with low n. The vapor pressures decrease sharply with n at low temperatures. This n dependence becomes increasingly weaker with increasing temperature, and no particular trend is observed at relatively high temperatures. In other words, increasing n makes it more difficult to transfer an ion pair from the liquid to vapor phase at low temperatures but not at high temperatures. This trend is also reflected in the normal boiling points shown in Figure 5. Here, we compare the n dependence of Tb obtained

vapor phase results and those of Rai and Maginn for 1,3dimethylimidazolium tetrafluoroborate ([C1mim][BF4]) were previously noted.11 However, at this point, we are unable to provide an explanation for this deviation. Inset b of Figure 3 shows magnified profiles at the lowest temperatures investigated in this work. At these temperatures, we observe the vapor density decreases with an increase in n. Vapor Pressures. Figure 4 provides vapor pressures psat for [Cnmim][NTf2] ionic liquids over a wide range of temperature.

Figure 4. Vapor pressures of RTILs belonging to the homologous series [Cnmim][NTf2] as a function of reciprocal temperature. Symbols are defined in the same manner as for Figure 2. Note that the pressures are plotted using a logarithmic scale. Inset a provides data at the highest temperatures studied. Here, squares have the same meaning as in Figure 2. Inset b shows vapor pressures at temperatures ranging from 450 to 500 K. Here, the black, red, green, blue, violet, and cyan diamonds represent experimental data from Rocha et al.5 for n = 2, 4, 6, 8, 10, and 12, respectively.

Figure 5. Boiling temperatures Tb as a function of n for RTILs belonging to the homologous series [Cnmim][NTf2]. Blue circles represent results from NPT-TEE simulations. Red squares correspond to results from the GEMC simulations of Rai and Maginn.3 Green diamonds represent values calculated using vapor pressures extrapolated from the experimental results of Zaitsau et al.6 The dashed lines are provided as a guide to the eye.

Insets a and b of Figure 4 show magnified profiles at high and moderate temperatures, respectively. We have also plotted the Gibbs ensemble results of Rai and Maginn.3 Their saturation pressures are systematically higher than our values (see inset a of Figure 4). As noted above, similar disagreement was observed in a previous study focused on [C1mim][BF4].11 We note that the approaches used to calculate the vapor pressure are rather different. Rai and Maginn employ the Gibbs ensemble and compute the pressure via the molecular virial. In our work, we use the NPT-TEE approach wherein the vapor pressure is set directly within the simulation and saturation conditions are established by equating Gibbs free energies. Hummer et al. have previously noted differences between the virial and thermodynamic pressures for water systems,38 and it appears that this issue should be investigated with ionic liquids in future work. Inset b of Figure 4 shows vapor pressures within the moderate temperature range 450−500 K. Here, we compare our values with those determined experimentally by Rocha et al.5 We find very good agreement between simulation and experimental vapor pressure data for n ≥ 4. In contrast, simulation data for the n = 2 case exceed experimental values by approximately a factor of 2. In a related study, Zhong et al.

from the current study with the simulation results of Rai and Maginn3 and values extrapolated from experimental measurements completed at temperatures ranging from 450 to 500 K.6 Our results indicate that the boiling point does not vary significantly with n. The difference between our results and those of Rai and Maginn is due to the difference in vapor pressures, as was discussed earlier. The boiling points based upon experimental data are obtained by assuming a linear relationship between ln psat and 1/T, as is discussed by Zaitsau et al.6 These extrapolated values are appreciably lower than our results and show a strong n dependence. This result is expected, as our simulation results suggest that ln psat varies nonlinearly with 1/T at relatively high temperatures. It is also clear from the data provided in Figure 4 that linear extrapolation of ln psat from the experimentally accessible temperature range will result in underestimation of the boiling temperature. Critical Properties. Figure 6 shows the dependence of the critical temperature and density on the length of the alkyl chain attached to the cation. Interestingly, we observe a maximum in the critical temperature at n = 4. In contrast, Rai and Maginn report a monotonic decrease in Tc with increasing n. We compute a lower value of Tc than Rai and Maginn for n = 2, E

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 6. Cation alkyl chain length dependence of the critical temperature and density for the [Cnmim][NTf2] homologous series. The upper and lower panels provide critical temperatures and densities, respectively. Blue symbols denote results obtained from this work using the Binder cumulant method. Violet symbols correspond to results obtained by Rai and Maginn3 using the law of rectilinear diameters approach.

Figure 7. Enthalpies of vaporization of RTILs belonging to the homologous series [Cnmim][NTf2] as a function of temperature. Symbol and color schemes are the same as those used in Figure 2. Insets a and b show the n dependence of ΔHvap at three different temperatures. Inset a: Blue circles and blue squares represent data from NPT-TEE simulations at 298 and 500 K, respectively. Red circles and squares represent data from molecular dynamics simulations of Zhong et al.19 at 298 and 500. K, respectively. Inset b: Blue circles denote values obtained at 900 K from NPT-TEE simulations, and red circles denote values obtained at 900 K from the Gibbs ensemble simulations of Rai and Maginn.3

whereas Tc values for n = 4 and 6 are relatively close. The discrepancy in Tc values for n = 2 is likely related to the differences in saturated liquid densities discussed earlier. The critical density monotonically decreases with increasing n. Our values for ρc are low relative to those reported by Rai and Maginn, with the difference between data sets reducing with increasing n. These differences are likely related to the different techniques used to determine the critical properties. Enthalpies of Vaporization. Figure 7 shows the temperature dependence of the enthalpy of vaporization ΔHvap for the homologous series [Cnmim][NTf2]. We calculate ΔHvap as follows ΔH vap = (U v − U l) + psat (V v − V l) i

We first compare our ΔHvap values to those from other simulation studies wherein the same force field was employed. Inset a of Figure 7 shows the n dependence of ΔHvap at 298 and 500 K. Our data are generally consistent with the results of Zhong et al.19 at 298 K, while our results at 500 K are systematically around 6 kJ/mol lower than the molecular dynamics results. Zhong et al. calculated ΔHvap via eq 1 and approximated the vapor phase as an ideal gas consisting of ion pairs. Our simulation results (e.g., compressibility factors, residual energies shown in Figure 10, cluster size distributions11) suggest that this is a reasonable approximation at these temperatures. Inset b of Figure 7 shows the n dependence of ΔHvap at 900 K. Here, our results differ from those of Rai and Maginn, with the two data sets suggesting qualitatively different trends regarding the variation in ΔHvap with n. We now compare our ΔHvap results with those from experimental studies. The procedures used to obtain ΔHvap from experimental data have drawn considerable attention within the recent literature.5,13 Typically, one collects some measure of volatility over a relatively narrow temperature range and uses these data to infer ΔHvap at the average collection temperature. Estimates for ΔHvap at a reference temperature of 298 K are then deduced from the relation

(1)

i

where U and V correspond to the average saturated molar configurational energy and volume of phase i, respectively. The superscripts l and v denote the liquid and vapor phase, respectively. All relevant quantities are obtained from NPTTEE simulations at liquid−vapor saturation conditions. We observe that ΔHvap varies strongly with n at low temperatures, whereas no particular trend is observed at high temperatures. This outcome is expected given the vapor pressures reported above, as ΔHvap and psat are related via the well-known Clausius−Clapeyron relation d ln psat d(1/T )

=−

ΔH vap psat ΔVvap/T

≅−

ΔH vap R

(2)

where ΔVvap = V − V represents the difference between the saturated molar volumes of the vapor and liquid phases. The second equality in eq 2 is applicable at temperatures sufficiently below the critical point, where one can approximate ΔVvap ≈ Vv and take the vapor to be an ideal gas, for which psatVv = RT, where R is the ideal gas constant. We note that estimates for ΔHvap generated via numerical implementation of eq 2 are consistent with those presented in Figure 7. We use eq 1 as our working equation simply because we find that it provides ΔHvap estimates that are more precise than those obtained from eq 2. v

l

ΔH vap(T ) = ΔH vap(T0) + ΔCP ,vap(T − T0)

(3)

where ΔCP,vap = − is the difference in the constant pressure molar heat capacity of the saturated vapor and liquid phases. In several earlier studies ΔCP,vap was simply taken to be a constant value of −100 J/mol-K.13,39 This approach resulted in considerable scatter in the estimates for ΔHvap at 298 K. CvP

F

C1P

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

and Rocha et al. data sets give an increase of 2.7, 4.2, and 3.9 kJ/mol per carbon atom at 298 K, respectively. Figure 9 shows the n dependence of ΔCP,vap. We obtain estimates for this quantity from a fit to ΔHvap data that span the

More recently, attention has been placed on determining accurate values for ΔCP,vap, such that reliable extrapolations to 298 K can be completed. We begin by focusing on ΔHvap comparisons at three temperatures that match the average temperature explored within the aforementioned volatility measurements, thereby minimizing the impact of the uncertainty in ΔCP,vap. Figure 8

Figure 9. Difference in saturated liquid and vapor constant pressure heat capacities as a function of the length of the alkyl chain attached to the cation. Blue circles represent values calculated from NPT-TEE simulations. Red squares denote experimental values reported by Verevkin et al.13 Green diamonds correspond to values calculated by Rocha at al.5 using an empirical relationship. Figure 8. Enthalpies of vaporization plotted as a function of n for the RTILs belonging to the homologous series [Cnmim][NTf2]. Insets a, b, c, and d provide values at 380, 475, 572, and 298 K, respectively. Blue circles represent results from NPT-TEE simulations. Green crosses and diamonds correspond to results from Rocha et al.5 at 475 and 298 K, respectively. Red squares, up triangles, and down triangles denote the results from Verevkin et al. at 298, 380, and 572 K, respectively.13

temperature range 380−572 K. We selected this temperature range to be consistent with the approach used by Verevkin et al.13 The three data sets are qualitatively similar, indicating that ΔCP,vap increases in magnitude with increasing n. Our data suggest that ΔCP,vap is relatively linear in n, with a decrement in ΔCP,vap of approximately 4.2 J/mol-K per carbon atom. We also note that our values are systematically lower than the frequently used value of −100 J/mol-K. Our values are relatively consistent with the ΔCP,vap(n) function deduced by Rocha et al. In contrast, the Verevkin et al. results suggest a stronger variation of ΔCP,vap with n relative to simulation data. We now further explore the origin of the chain length dependence of ΔHvap. Consideration of the magnitudes of the energy and pressure−volume terms in eq 1 reveals that the energy term dominates. For example, at 400 K, the energy term accounts for more than 97% of ΔHvap. Therefore, we focus on this term. Figure 10 provides the temperature dependence of the residual vapor and liquid energies of the [Cnmim][NTf2] series. Here, the residual energy is defined as the difference between the total and ideal gas energies, where we take the ideal gas to consist of ion pairs. We compute this ideal gas energy by determining the ensemble average energy of an isolated ion pair. The residual energies of the vapor phase collapse to zero at low to moderate temperature, indicating that this quantity does not vary with n. In contrast, the liquid phase residual energies clearly show an n dependence similar to that for ΔHvap, with the variation in the residual energy increasing with decreasing temperature. Effect of Anions. We now briefly examine the impact the anion has on liquid−vapor saturation properties. Saturated densities and vapor pressures for [C4mim][NTf2], [C4mim][BF4], and [C4mim][PF6] are provided in Figures 11 and 12, respectively. The liquid density decreases in the order [NTf2] >

provides data from quartz crystal microbalance measurements of Verevkin et al. at 380 K (inset a),13 Knudsen effusion measurements of Rocha et al. at 475 K (inset b),5 and thermogravimetry measurements of Verevkin et al. at 572 K (inset c).13 First, we note that our data vary linearly with the chain length at all temperatures. At 380 K, experimental values are systematically 5−10 kJ/mol lower than ours. At 475 K, experimental and simulation results are in reasonably good agreement. At 572 K, experimental data show a nonlinear dependence with respect to n, with good agreement between experimental and simulation data at moderate n and relatively large deviations at higher n. Figure 8d provides our ΔHvap values at 298 K as well as those extrapolated from experimental data. The Verevkin et al. results are generated using eq 3, with ΔCP,vap deduced from the ΔHvap data collected at 380 and 572 K. We find that the level of agreement between our data and the Verevkin et al. data increases with increasing n. Rocha et al. employed eq 3 and an n-dependent ΔCP,vap relation deduced from a fit of available experimental data (see ref 5 for details) to generate estimates for ΔHvap at 298 K. Here, we find that the level of agreement between our data and the Rocha et al. data decreases with increasing n. Both experimental data sets suggest a stronger variation in ΔHvap with n than the results obtained from simulation. More specifically, the simulation, Verevkin et al., G

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 10. Liquid and vapor phase residual energies of RTILs belonging to the homologous series [Cnmim][NTf2] as a function of temperature. Details regarding the method used to calculate residual energies are provided within the text. Symbol and color schemes are the same as those used in Figure 2.

Figure 12. Influence of the anion identity on the vapor pressure. Note that the vapor pressures are plotted using a logarithmic scale. Inset a shows vapor pressures at relatively high temperatures. Gibbs ensemble simulation results of Rai and Maginn3,4 are denoted by squares. Inset b shows enthalpies of vaporization as a function of temperature. Symbol and color schemes are the same as those used in Figure 11.

while for high temperatures they are similar. The results for vapor pressures show trends similar to those for vapor densities. Inset a of Figure 12 provides a comparison between our vapor pressure results and those from the Gibbs ensemble study of Rai and Maginn. Inset b of Figure 12 provides the enthalpies of vaporization for [C 4 mim][NTf 2 ], [C 4 mim][BF 4 ], and [C4mim][PF6]. We observe that ΔHvap is significantly influenced by the nature of the anion. The profile for [C4mim][NTf2] exhibits considerable curvature at low temperature. Using the approach outlined above, we arrive at estimates of ΔCP,vap = −111, −87, and −94 J/mol-K for [C4mim][NTf2], [C4mim][BF4], and [C4mim][PF6], respectively. Overall, these results show that the nature of anion has a significant impact on liquid−vapor saturation properties. Moreover, the influence of the anion identity does not vary significantly with temperature. Liquid Phase Structure. Previous studies indicate that the liquid phase of RTILs with relatively long alkyl chains organizes into polar and nonpolar “nano-domains”.21−23 For the homologous series [Cnmim][NTf2], such structures are generally observed for n > 4.22 In this work, we focus on the temperature dependence of these domains. It is generally expected that nanodomains will feature prominently at low temperature and that entropic driving forces will frustrate such organization at high temperature. Their presence is usually characterized by a pronounced first peak in the radial distribution function (rdf) associated with the terminal methyl groups of the alkyl chain attached to the imidazolium cation. In the top panel of Figure 13, we plot rdf’s for the ionic liquid [C8mim][NTf2] at four points along the liquid−vapor saturation line. One finds that the height of the first peak decreases with temperature, with the peak vanishing at T = 950 K, suggesting the liquid structure is relatively homogeneous at this temperature. In the bottom panel of Figure 13, we plot the height of the first peak of the rdf as a function of temperature.

Figure 11. Influence of anion identity on the saturated liquid (top panel) and vapor (lower panel) densities. Solid black, red, and green lines represent NPT-TEE simulation data for [C4mim][NTf2], [C4mim][BF4], and [C4mim][PF6], respectively. The inset in the lower panel shows the magnified vapor density profiles at relatively low temperatures. Note that the vapor densities are plotted using a logarithmic scale. Squares denote GEMC simulation results of Rai and Maginn,3,4 and triangles denote experimental results from Jacquemin et al.40 The color scheme used for symbols is the same as that for solid lines.

[PF6] > [BF4] over the entire temperature range. We also observe a large difference between the liquid densities of [C4mim][BF4] and those of the other two RTILs at lower temperatures. The saturated vapor densities of [C4mim][NTf2] are much higher than those of the other two RTILs, especially at high temperatures. The vapor densities of the other two RTILs are in the order [PF6] < [BF4] at low temperatures, H

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

temperatures. We observed that this property increases with the cation alkyl chain length at fixed temperature, which is qualitatively consistent with trends reported within experimental studies. Quantitative agreement between simulation data and experimental results was also reasonable. Our results indicate that the alkyl chain length dependence of vapor pressures and enthalpies of vaporization weakens with increasing temperature. The ln psat versus 1/T vapor pressure relation exhibits considerable curvature, which limits the extent to which one can use low-temperature data to estimate boiling temperatures. Finally, we briefly looked at the influence of temperature on the nanoscale structure of the liquid phase of RTILs with relatively long alkyl chains. We did not find evidence for a sharp structural phase transition while traversing from a nanostructured liquid at low temperatures to a more homogeneous liquid phase at high temperatures. The current work demonstrates the utility of Monte Carlo simulation and histogram based strategies for computing the liquid−vapor saturation properties of realistic room temperature ionic liquids at experimentally relevant conditions. The approach pursued here provides pseudo-continuous saturation curves over a broad range of temperature with a reasonably high level of precision. These results enable one to identify subtle trends in thermophysical properties of interest, and to obtain differential properties such as the specific heat. We are currently extending these techniques to develop strategies41 to study the interfacial properties of ionic liquids. These methods require as input bulk saturation properties.

Figure 13. Influence of temperature on the liquid phase structure of [C8mim][NTf2]. The top panel shows the intermolecular radial distribution functions (rdf’s) associated with the terminal methyl groups of the alkyl chain attached to the cation. Solid black, red, green, and blue lines represent the results at 298, 450, 700, and 950 K, respectively. The lower panel provides the height of the first peak of the above-mentioned rdf as a function of temperature.

■

We observe a gradual decrease in peak height, which indicates that the distribution of alkyl chains evolves gradually with temperature. In a previous study, Wang and Voth used a coarse-grained model of an imidazolium based ionic liquid to study the temperature dependence of liquid nanodomains.24 They observed the characteristics of a phase transition while going from low to high temperature. Specifically, they identified a discontinuity in the liquid phase energies and a pronounced peak in the constant volume specific heat. They proposed that this transition is related to the formation of nanodomains at low temperature. In contrast, we do not observe a drastic change in the liquid phase energy while moving from low to high temperature (see Figure 10). Collectively, the configurational energy and structural metrics examined here suggest that [C8mim][NTf2] does not exhibit a nanodomain related phase transition upon variation of the temperature. Further studies are needed to understand the generality of this result.

ASSOCIATED CONTENT

S Supporting Information *

The Xmgrace source files for Figures 1−13 are provided in a zip archive. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: jerring@buffalo.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge the financial support of the National Science Foundation (Grant No. CHE-1012356). Computational resources were provided in part by the University at Buffalo Center for Computational Research and the Rensselaer Polytechnic Institute Computational Center for Nanotechnology Innovations.

4. CONCLUSIONS We employed Monte Carlo simulation to calculate the liquid− vapor saturation properties of imidazolium based RTILs over a range of temperatures spanning from 280 to 1000 K. We focused on understanding the influence of the cation alkyl chain length on saturation properties of RTILs belonging to the [Cnmim][NTf2] homologous series. We also provided results for [C4mim][BF4] and [C4mim][PF6], which helped us to understand the effect of anion identity. The saturation properties were obtained via a histogram based scheme that was introduced in a previous study. Our results for saturated densities and vapor pressures were compared with those from experimental and simulation studies. Overall, the computed properties show good agreement with those obtained from experiment, suggesting that the Zhong et al. model provides a reasonable description of saturation properties. We also calculated enthalpies of vaporization over a wide range of

■

REFERENCES

(1) Cvjetko Bubalo, M.; Radošević, K.; Radojčić Redovniković, I.; Halambek, J.; Gaurina Srček, V. A Brief Overview of the Potential Environmental Hazards of Ionic Liquids. Ecotoxicol. Environ. Saf. 2014, 99, 1−12. (2) Maginn, E. J. Molecular Simulation of Ionic Liquids: Current Status and Future Opportunities. J. Phys.: Condens. Matter 2009, 21, 373101−373117. (3) Rai, N.; Maginn, E. J. Critical Behaviour and Vapour-Liquid Coexistence of 1-Alkyl-3-Methylimidazolium Bis(Trifluoromethylsulfonyl)Amide Ionic Liquids Via Monte Carlo Simulations. Faraday Discuss. 2012, 154, 53−69.

I

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

(22) Canongia Lopes, J. N. A.; Pádua, A. A. H. Nanostructural Organization in Ionic Liquids. J. Phys. Chem. B 2006, 110, 3330−3335. (23) Triolo, A.; Russina, O.; Bleif, H.-J.; Di, C. E. Nanoscale Segregation in Room Temperature Ionic Liquids. J. Phys. Chem. B 2007, 111, 4641−4644. (24) Wang, Y.; Voth, G. A. Tail Aggregation and Domain Diffusion in Ionic Liquids. J. Phys. Chem. B 2006, 110, 18601−18608. (25) Greaves, T. L.; Drummond, C. J. Solvent Nanostructure, the Solvophobic Effect and Amphiphile Self-Assembly in Ionic Liquids. Chem. Soc. Rev. 2013, 42, 1096−1120. (26) Castner, E. W., Jr.; Margulis, C. J.; Maroncelli, M.; Wishart, J. F. Ionic Liquids: Structure and Photochemical Reactions. Annu. Rev. Phys. Chem. 2011, 62, 85−105. (27) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (28) Binder, K. Computer Simulations of Critical Phenomena and Phase Behaviour of Fluids. Mol. Phys. 2010, 108, 1797−1815. (29) Zollweg, J. A.; Mulholland, G. W. On the Law of the Rectilinear Diameter. J. Chem. Phys. 1972, 57, 1021−1025. (30) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo Technique for Studying Phase Transitions. Phys. Rev. Lett. 1988, 61, 2635−2638. (31) Errington, J. R. Direct Calculation of Liquid−Vapor Phase Equilibria from Transition Matrix Monte Carlo Simulation. J. Chem. Phys. 2003, 118, 9915−9925. (32) Errington, J. R.; Panagiotopoulos, A. Z. New Intermolecular Potential Models for Benzene and Cyclohexane. J. Chem. Phys. 1999, 111, 9731−9738. (33) Mehlig, B.; Heermann, D. W.; Forrest, B. M. Hybrid Monte Carlo Method for Condensed-Matter Systems. Phys. Rev. B 1992, 45, 679. (34) Chen, B.; Siepmann, J. I. A Novel Monte Carlo Algorithm for Simulating Strongly Associating Fluids: Applications to Water, Hydrogen Fluoride, and Acetic Acid. J. Phys. Chem. B 2000, 104, 8725−8734. (35) Chen, B.; Siepmann, J. I. Improving the Efficiency of the Aggregation−Volume−Bias Monte Carlo Algorithm. J. Phys. Chem. B 2001, 105, 11275−11282. (36) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic Press: London, 2001; Vol. 1. (37) Huddleston, J. G.; Visser, A. E.; Reichert, W. M.; Willauer, H. D.; Broker, G. A.; Rogers, R. D. Characterization and Comparison of Hydrophilic and Hydrophobic Room Temperature Ionic Liquids Incorporating the Imidazolium Cation. Green Chem. 2001, 3, 156− 164. (38) Hummer, G.; Gronbech-Jensen, N. Pressure Calculation in Polar and Charged Systems Using Ewald Summation: Results for the Extended Simple Point Charge Model of Water. J. Chem. Phys. 1998, 109, 2791. (39) M. S. S. Esperança, J.; Canongia Lopes, J. N.; Tariq, M.; Santos, L. s. M. N. B. F.; Magee, J. W.; Rebelo, L. s. P. N. Volatility of Aprotic Ionic Liquids  a Review. J. Chem. Eng. Data 2009, 55, 3−12. (40) Jacquemin, J.; Husson, P.; Mayer, V.; Cibulka, I. High-Pressure Volumetric Properties of Imidazolium-Based Ionic Liquids: Effect of the Anion. J. Chem. Eng. Data 2007, 52, 2204−2211. (41) Rane, K. S.; Kumar, V.; Errington, J. R. Monte Carlo Simulation Methods for Computing the Wetting and Drying Properties of Model Systems. J. Chem. Phys. 2011, 135, 234102−234113.

(4) Rai, N.; Maginn, E. J. Vapor−Liquid Coexistence and Critical Behavior of Ionic Liquids Via Molecular Simulations. J. Phys. Chem. Lett. 2011, 2, 1439−1443. (5) Rocha, M. A. A.; Lima, C. F. R. A. C.; Gomes, L. R.; Schroder, B.; Coutinho, J. A. P.; Marrucho, I. M.; Esperanca, J. M. S. S.; Rebelo, L. P. N.; Shimizu, K.; Lopes, J. N. C.; Santos, L. M. N. B. F. High-Accuracy Vapor Pressure Data of the Extended [Cnc1im][Ntf2] Ionic Liquid Series: Trend Changes and Structural Shifts. J. Phys. Chem. B 2011, 115, 10919−10926. (6) Zaitsau, D. H.; Kabo, G. J.; Strechan, A. A.; Paulechka, Y. U.; Tschersich, A.; Verevkin, S. P.; Heintz, A. Experimental Vapor Pressures of 1-Alkyl-3-Methylimidazolium Bis(Trifluoromethylsulfonyl)Imides and a Correlation Scheme for Estimation of Vaporization Enthalpies of Ionic Liquids. J. Phys. Chem. A 2006, 110, 7303−7306. (7) Valderrama, J. O.; Zarricueta, K. A Simple and Generalized Model for Predicting the Density of Ionic Liquids. Fluid Phase Equilib. 2009, 275, 145−151. (8) Shiflett, M. B.; Yokozeki, A. Solubilities and Diffusivities of Carbon Dioxide in Ionic Liquids: [Bmim][Pf6] and [Bmim][Bf4]. Ind. Eng. Chem. Res. 2005, 44, 4453−4464. (9) Weiss, V. C.; Heggen, B.; Muller-Plathe, F. Critical Parameters and Surface Tension of the Room Temperature Ionic Liquid [Bmim][Pf6]: A Corresponding-States Analysis of Experimental and New Simulation Data. J. Phys. Chem. C 2010, 114, 3599−3608. (10) Rebelo, L. P. N.; Lopes, J. N. C.; Esperanca, J. M. S. S.; Filipe, E. On the Critical Temperature, Normal Boiling Point, and Vapor Pressure of Ionic Liquids. J. Phys. Chem. B 2005, 109, 6040−6043. (11) Rane, K. S.; Errington, J. R. Using Monte Carlo Simulation to Compute Liquid−Vapor Saturation Properties of Ionic Liquids. J. Phys. Chem. B 2013, 117, 8018−8030. (12) Rane, K. S.; Murali, S.; Errington, J. R. Monte Carlo Simulation Methods for Computing Liquid−Vapor Saturation Properties of Model Systems. J. Chem. Theory Comput. 2013, 9, 2552−2566. (13) Verevkin, S. P.; Zaitsau, D. H.; Emel’yanenko, V. N.; Yermalayeu, A. V.; Schick, C.; Liu, H.; Maginn, E. J.; Bulut, S.; Krossing, I.; Kalb, R. Making Sense of Enthalpy of Vaporization Trends for Ionic Liquids: New Experimental and Simulation Data Show a Simple Linear Relationship and Help Reconcile Previous Data. J. Phys. Chem. B 2013, 117, 6473−6486. (14) Tomé, L. I. N.; Carvalho, P. J.; Freire, M. G.; Marrucho, I. M.; Fonseca, I. M. A.; Ferreira, A. G. M.; Coutinho, J. o. A. P.; Gardas, R. L. Measurements and Correlation of High-Pressure Densities of Imidazolium-Based Ionic Liquids. J. Chem. Eng. Data 2008, 53, 1914− 1921. (15) Gardas, R. L.; Freire, M. G.; Carvalho, P. J.; Marrucho, I. M.; Fonseca, I. M. A.; Ferreira, A. G. M.; Coutinho, J. A. P. Pρt Measurements of Imidazolium-Based Ionic Liquids. J. Chem. Eng. Data 2007, 52, 1881−1888. (16) Fredlake, C. P.; Crosthwaite, J. M.; Hert, D. G.; Aki, S. N. V. K.; Brennecke, J. F. Thermophysical Properties of Imidazolium-Based Ionic Liquids. J. Chem. Eng. Data 2004, 49, 954−964. (17) Gomes de Azevedo, R.; Esperança, J. M. S. S.; Szydlowski, J.; Visak, Z. P.; Pires, P. F.; Guedes, H. J. R.; Rebelo, L. P. N. Thermophysical and Thermodynamic Properties of Ionic Liquids over an Extended Pressure Range: [Bmim][Ntf2] and [Hmim][Ntf2]. J. Chem. Thermodyn. 2005, 37, 888−899. (18) Yoo, B.; Afzal, W.; Prausnitz, J. M. Solubility Parameters for Nine Ionic Liquids. Ind. Eng. Chem. Res. 2012, 51, 9913−9917. (19) Zhong, X.; Liu, Z.; Cao, D. Improved Classical United-Atom Force Field for Imidazolium-Based Ionic Liquids: Tetrafluoroborate, Hexafluorophosphate, Methylsulfate, Trifluoromethylsulfonate, Acetate, Trifluoroacetate, and Bis(Trifluoromethylsulfonyl)Amide. J. Phys. Chem. B 2011, 115, 10027−10040. (20) Ponder, J. W.; Case, D. A. Force Fields for Protein Simulations. In Advances in Protein Chemistry, 1st ed.; Valerie, D., Ed.; Academic Press: Waltham, MA, 2003; Vol. 66, pp 27−85. (21) Wang, Y.; Voth, G. A. Unique Spatial Heterogeneity in Ionic Liquids. J. Am. Chem. Soc. 2005, 127, 12192−12193. J

dx.doi.org/10.1021/jp504085t | J. Phys. Chem. B XXXX, XXX, XXX−XXX