Computationally Efficient Prediction of Ionic Liquid Properties

Letter pubs.acs.org/JPCL

Computationally Efficient Prediction of Ionic Liquid Properties Vitaly V. Chaban*,†,‡ and Oleg V. Prezhdo‡ †

MEMPHYS − Center for Biomembrane Physics, Syddansk Universitet, Odense M. 5230, Kingdom of Denmark Department of Chemistry, University of Rochester, Rochester, New York 14627, United States

‡

ABSTRACT: Due to fundamental differences, room-temperature ionic liquids (RTIL) are significantly more viscous than conventional molecular liquids and require long simulation times. At the same time, RTILs remain in the liquid state over a much broader temperature range than the ordinary liquids. We exploit the ability of RTILs to stay liquid at several hundred degrees Celsius and introduce a straightforward and computationally efficient method for predicting RTIL properties at ambient temperature. RTILs do not alter phase behavior at 600−800 K. Therefore, their properties can be smoothly extrapolated down to ambient temperatures. We numerically prove the validity of the proposed concept for density and ionic diffusion of four different RTILs. This simple method enhances the computational efficiency of the existing simulation approaches as applied to RTILs by more than an order of magnitude. SECTION: Liquids; Chemical and Dynamical Processes in Solution

A

is a serious challenge for laboratory work with ionic liquids. It also presents problems for computational studies7,8 as viscous liquids normally exhibit a very complicated potential energy surface. Extensive sampling efforts are necessary to derive reliable structure and dynamical properties and, especially, to conduct free Gibbs energy simulations. While 1 ns long sampling imposed on a system near the free-energy minimum provides accurate results for the majority of ordinary organic solvents, RTILs require tens and often hundreds of nanoseconds to get analogous physical insights. Thus, investigation of RTILs is more expensive in both laboratory and mathematical modeling.11,12 In this Letter, we demonstrate a simple and computationally efficient methodology for predicting properties of RTILs. The method reduces the simulation effort by over an order of magnitude. We make use of the fact that RTILs maintain the same liquid phase over a very large temperature interval. The normal evaporation temperatures of RTILs are above 800 K and can be much higher depending on the ion sizes. Further, many RTILs decompose below the normal evaporation point and far below the critical point. Aqueous vapor is believed to play an important role in thermally driven decomposition of RTILs, while ionic liquids are notoriously stable in vacuum. Problematic in experimental studies, decomposition can be forbidden in atomistic simulations. Therefore, a complete phase diagram of a RTIL can be recorded at a low computational expense using the proposed technique. For a representative set of RTILs, including 1-ethyl-3methylimidazolium tetrafluoroborate, [C2C1IM][BF4], 1-ethyl3-methylimidazolium bis(trifluoromethylsulfonyl)imide,

room-temperature ionic liquid (RTIL) is a specific kind of salt in which bulky asymmetric cations and primarily multiatom anions are poorly coordinated.1−6 Such a structure ensures that these compounds remain liquid below 100 °C or even at room temperature. The strong intermolecular interactions between the charged particles increase the viscosity and extend the temperature range over which RTILs exist in the liquid state. Delocalized electronic charge is another key feature of many RTILs.7 At least one ion is organic, that is, it does not contain metallic elements. This favors specific valence electronic structure and prevents RTILs from forming an energetically stable crystal lattice at lower temperatures.5,7,8 The purely organic methylimidazolium, [C1IM]+, and pyridinium, [PY]+, cations are excellent starting points for construction of a great variety of stable RTILs. Multiple organic and inorganic anions, such as tetrafluoroborate, [BF 4 ] − , and bis(trifluoromethylsulfonyl)imide, [TFSI]−, can diversify physical properties of these two classes of RTILs. A particularly important feature, many RTILs decompose into harmless products and are considered environmentally friendly. Note that recent studies cast some doubts regarding the extent of their green nature.9 RTIL properties, such as melting and evaporation points, electric conductivity, shear viscosity, and density, can be controlled by substituents in the organic component and, largely, by the counterion.5 Many ionic liquids have already been developed for specific synthetic problems and robust separation applications. RTILs are worthily termed “designer solvents”. Although RTILs at ambient conditions maintain arrangement characteristic to liquids, the viscosity of these substances is much higher than that of conventional liquids. The viscosity of certain RTILs reaches over 100 cP (100 times more than the viscosity of liquid water). The viscosity is dramatically increased by addition of lyophobic chains to the cation.8,10 Huge viscosity © XXXX American Chemical Society

Received: April 9, 2014 Accepted: May 20, 2014

1973

dx.doi.org/10.1021/jz5007127 | J. Phys. Chem. Lett. 2014, 5, 1973−1977

The Journal of Physical Chemistry Letters

Letter

[C 2 C 1 IM][TFSI], N-butylpyridinium tetrafluoroborate, [C4PY][BF4], and N-butylpyridinium bis(trifluoromethylsulfonyl)imide, [C4PY][TFSI], we show that self-diffusion at ambient conditions can be predicted based on the simulations at significantly higher temperatures. High-temperature simulations take exponentially less time because the required sampling is achieved fast. The viscosity is linked to the average self-diffusion constant of the system via the Stokes−Einstein equation, and viscosity rapidly (exponentially) decays with increasing temperature. To validate the force field models, we use density (Figure 1) as this property is easily measurable and

Figure 1. Mass density of [C2C1IM][BF4], [C2C1IM][TFSI], [C4PY][BF4], and [C4PY][TFSI] simulated at 600−800 K and extrapolated to ambient temperatures. All extrapolated values agree well with the experimental data. The standard deviations of the computed densities are 1−2% of their converged values.

Figure 2. (Top panel) Percent deviation between the simulated and experimental mass densities of [C2C1IM][BF4], [C2C1IM][TFSI], [C4PY][BF4], and [C4PY][TFSI]. (Bottom panel) Point-by-point comparison of the experimental (right column) and simulated (left column) mass densities. The maximum deviation of the extrapolated values from the experimental values amounts to 4%. Note that throughout the work, we compare the extrapolated values to the experimental ones and do not compare these values to the simulations at room temperature. One should distinguish between inaccuracies from the atomistic models and inaccuracies from the mathematically processed computer simulation data.

is not susceptible to large fluctuations in molecular dynamics (MD) simulations. It is important for the proposed approach that a RTIL remains in the liquid state at the selected simulation temperatures. We performed MD simulations of [C 2 C 1 IM][BF 4 ], [C2C1IM][TFSI], [C4PY][BF4], and [C4PY][TFSI] at 600, 650, 700, 750, and 800 K to obtain mass densities (Figure 1). Each density was obtained using an independent simulation, where constant temperature was controlled. The obtained temperature dependence of density is perfectly described by the first-order polynomial; the correlation coefficients exceed 0.99 for all four systems. Extrapolated to room temperature, the simulation provides 1236 kg m−3 for [C2C1IM][BF4], 1579 kg m−3 for [C2C1IM][TFSI], 1167 kg m−3 for [C4PY][TFSI], and 1494 kg m−3 for [C4PY][TFSI] (Figure 2). The corresponding experimental densities13 are 1275 kg m−3 for [C2C1IM][BF4], 1512 kg m−3 for [C2C1IM][TFSI], 1216 kg m−3 for [C4PY][TFSI], and 1444 kg m−3 for [C4PY][TFSI]. All simulated and subsequently extrapolated densities are in good agreement with the experimental data. Interestingly, the densities are slightly overestimated for the RTILs containing [TFSI]− anions and are slightly underestimated for the two RTILs containing tetrafluoroborate. The largest deviation is 4%, observed for [C2C1IM][TFSI]. The relative values of the densities of the studied systems are accurately reproduced. It is important to note that for the sake of generality, the applied force field has been derived exclusively on the basis of ab initio quantum chemical data for isolated cations and anions.14 Certain interest is credited to zero-temperature densities. Hypothetical values, these densities are determined under the assumption that a RTIL is brought down to 0 K without a liquid−solid phase transition. Experience shows that the difference between the solid and liquid densities of the same compound is smaller in many cases than the difference in the densities of different compounds. Therefore, this measure can provide important insights about both thermal compressibility

coefficients, and the most energetically favorable packing of various ions. The zero-temperature density of the investigated RTILs amounts to 1430 kg m−3 for [C 2C 1IM][BF4], 1844 kg m−3 for [C2C1IM][TFSI], 1353 kg m−3 for [C4PY][BF4], and 1747 kg m−3 for [C4PY][TFSI]. RTILs containing [TFSI] exhibit higher densities mostly because [TFSI]− is constructed of heavier elements than tetrafluoroborate. The same considerations apply for comparison between [C4PY]+ and [C2C1IM]+. Furthermore, the greater affinity of [BF4]− to the imidazole ring than to the pyridine ring is easy to interpret from the chemical perspective. The cationic and anionic self-diffusion constants are summarized in Figure 3. As all transport properties, diffusion accelerates exponentially with an increase in the kinetic energy of a system. The temperature dependence can be approximated using a power function, such as D = A × 10−(B/T), where D is the self-diffusion constant, defined separately for the cation and the anion, T is temperature in Kelvin, and A and B are constants that reflect the physical nature of each RTIL. The power function can be rewritten for simplicity as log D = log A − (B/T) (Figure 3). The two-parameter function is sufficient to describe the five simulated self-diffusion constant values for each RTIL. The correlation coefficients of the fits were 0.98 or greater, which can be regarded as an excellent approximation. Several arguments motivate us to use the simple D = A × 10−(B/T) expression instead of the frequently applied13,15−18 and more sophisticated Vogel−Tamman−Fulcher (VTF) equation. The VTF kinetics is followed around the glass transition temperature. For the RTILs investigated in this work, the glass 1974

dx.doi.org/10.1021/jz5007127 | J. Phys. Chem. Lett. 2014, 5, 1973−1977

The Journal of Physical Chemistry Letters

Letter

Figure 4. Point-to-point comparison of self-diffusion constants obtained from pulsed-gradient spin−echo 1H and 19F NMR (left columns)13 and through extrapolation of the simulated data at high temperatures (right columns).

Figure 3. Self-diffusion constants (D) simulated at 600−800 K. See the legends for description. The dotted lines were obtained by fitting of the simulated numbers to the proposed mathematical function. The standard deviations of simulated diffusion constants are 10−20% of their original values.

transition temperatures are 184 ([C 2C1IM][BF4]), 186 ([C2C1IM][TFSI]), and 202 K (C4PY][BF4]).13 At higher temperatures, the VTF regime reduces to Arrhenius kinetics. In this case, fitting to the two-parameter exponential expression is preferable to the use of the three-parameter VTF equation. Caution must be taken when applying this method to study low-temperature transport properties. Diffusion below the melting temperature is particularly slow, corresponding to shear viscosity above 1000 cP. Usually, it presents little practical interest. The transition between the VTF and Arrhenius behaviors deserves particular attention. Figure 4 provides a point-to-point comparison between the simulated data and the results of direct experimental measurement using nuclear magnetic resonance (NMR).13 The diffusion constants at room conditions are notoriously small, complicating high-precision investigation both in simulations and experimentally. The largest self-diffusion constant is observed for the [C2C1IM]+ cation in [C2C1IM][BF4]. Interestingly, this conclusion is only partially confirmed by the experiment, where [C2C1IM]+ is the fastest ion, but only in [C2C1IM][TFSI]. [C2C1IM][TFSI] exhibits the highest zerotemperature density, which rather implies reduced ionic mobility. According to NMR,13 TFSI is not only more mobile than [BF4]− but also makes the corresponding RTILs more mobile as a whole. This is in spite of its large size and bulky shape. [C4PY]+ is less mobile than [C2C1IM]+, regardless of the anion used. The discussed differences are, however, negligibly small in absolute units, keeping the trends exclusively of theoretical interest. The applied model allows us to separate the self-diffusion constant of RTILs into the temperature-dependent, B/T, and temperature-independent, log A, terms (see Figure 5). The temperature-independent term, A, is similar for all investigated RTILs, while the B values vary more significantly. Larger B

Figure 5. Temperature-independent (log A) and temperature-dependent (B) factors determining ionic self-diffusion constants in [C2C1IM][BF4], [C2C1IM][TFSI], [C4PY][TFSI], and [C4PY][TFSI] ionic liquids. The left columns correspond to cations, and the right columns correspond to anions of the marked RTILs.

means weaker dependence of the diffusion constant on temperature. Both anions possess larger B, regardless of the cation used. Cations appear more mobile in general, with the only exclusion of [C4PY][BF4] at some elevated temperatures (Figure 3). The computed values can be used to forecast the phase transition points of the RTILs. The remainder of the Letter is devoted to comparison between the sampling efforts necessary at room and elevated temperatures (Figure 6). The energy landscape of liquids near the glass transition point is enormously complicated, making it impossible to accumulate adequate statistics without techniques for enhanced sampling, such as replica exchange. The average self-diffusion constant of [C4PY][TFSI] at 800 K equals 3.3 × 10−9 m2 s−1. This value is somewhat larger than D in room-temperature water and 3 orders of magnitude 1975

dx.doi.org/10.1021/jz5007127 | J. Phys. Chem. Lett. 2014, 5, 1973−1977

The Journal of Physical Chemistry Letters

Letter

is required. Diffusion constants between room and vaporization/decomposition temperatures of each RTIL can be promptly derived. Electronic polarization effects7,8,19 between the cation and anion are another problem in most RTILs at ambient conditions. The problem is circumvented using the newly introduced method. At 600−800 K, polarization decreases dramatically, and nonpolarizable models can be used, providing the same level of accuracy as the polarizable force fields.20,21 Neglect of polarization at room temperature would have produced systematically underestimated values of diffusion constants. Performance of the Arrhenius-type fitting equations does not depend on the type of interaction potential. Different force field models may be needed for various temperature ranges. The fitting will provide correct predictions as long as the diffusion constant at each considered temperature is simulated with a reasonable accuracy. Using nonpolarizable models, we are able to evaluate accurately and efficiently the RTIL properties at the elevated temperatures and obtain by extrapolation the low-temperature data in good agreement with experiment.

■

METHODOLOGY The presented analysis is based on MD simulation of 21 systems (Table 1) involving four RTILs. These four RTILs

Figure 6. Deviation of the average diffusion constant, D, obtained using different sampling times, from the converged value. The converged D was obtained using 40 ns trajectories at 300 and 800 K. All ions were used for the calculation simultaneously. Therefore, the depicted D represents an average system D. The analysis is performed with the [C4PY][TFSI] system (150 ion pairs), as an example.

Table 1. List of the Systems Simulated in the Present Worka

larger than D of [C4PY][TFSI] at 300 K. On the basis of our convergence analysis (Figure 6), the diffusion constant at 800 K can be reliably estimated based on very short trajectories, well below 1 ns. This is in complete agreement with the sampling requirements for small solvent molecules, such as water, methanol, acetonitrile, and so forth, simulated at 300 K or above. However, simulations of [C4PY][TFSI] diffusion at 300 K require at least 10 ns of sampling. [C4PY][TFSI] follows a subdiffusion regime during a few nanoseconds. Thus, the resulting D is overestimated if proper sampling of mean-square displacements is not achieved. Furthermore, equilibration at 300 K requires additional 10 ns of MD simulation. To recapitulate, we proposed and numerically validated a very straightforward method to predict RTIL transport properties, such as ionic self-diffusion constants, at ambient conditions using MD simulations at substantially elevated temperatures. In comparison, our method provides the diffusion constant for [C4PY][TFSI] using only 0.5−1 ns (five simulations at 600, 650, 700, 750, and 800 K using 100− 200 ps per MD run), whereas the direct approach requires at least 10 ns; plus, similar time is needed to achieve the freeenergy minimum state. The computational cost is decreased by over 10 times. Even if certain RTILs decompose below 600− 800 K in experiments, for example, due to the presence of water, we can still perform the simulations using the same protocol and achieve the order of magnitude speedup. With high reliability, the obtained properties are extrapolated to the room-temperature range. The proposed method can be used to predict RTIL properties over the entire temperature range of the RTIL liquid phase, provided that the temperature dependence of the given property is well represented by the chosen analytic expression, such as the Arrhenius and VTF equations. The introduced method can be particularly useful in cases where screening of a large set of RTIL-containing systems

#

RTIL

ion pairs

interaction centers

sampling time, ns

1−5 6−10 11−15 16−21

[C2C1IM][BF4] [C2C1IM][TFSI] [C4PY][BF4] [C4PY][TFSI]

250 200 250 150

6000 6800 7250 5850

50 40 75 40

a

Each RTIL was simulated at 600, 650, 700, 750, and 800 K. Additionally, the MD system containing [C4PY][TFSI] was simulated at 300 K in order to record convergence of ionic diffusion constants with simulation time.

represent two independent families of compounds. A 0.001 ps MD integration time step was used. The simulation time for each system is indicated in Table 1. Each ionic composition was simulated at 600, 650, 700, 750, and 800 K. Additionally, the [C4PY][TFSI] system was simulated at 300 K for 40 ns in order to obtain the data for Figure 6. The number of ions in each system was determined by the cation and anion sizes and the total number of atoms. All systems contain approximately the same number of interaction centers to ensure similar computational efforts. The first 1 ns (at 600−800 K) was disregarded as equilibration, and the equilibrium properties were obtained based on the remaining trajectory part. The coordinates were saved every 0.01 ps, that is, every 10 time steps. The MD trajectories were propagated using the GROMACS simulation suite.22−24 The computation of density and diffusion was carried out using the GROMACS standard analysis tools.22−24 The mathematical processing of data was performed using in-home scripts. The force field models for the [C2C1IM][BF4], [C2C1IM][TFSI], [C4PY][TFSI], and [C4PY][TFSI] ionic liquids were developed by Lopes and Padua recently.14 We did not modify these models. The comparison versus the experimental data (Figures 2 and 4) indicates a good to excellent performance of the force fields. 1976

dx.doi.org/10.1021/jz5007127 | J. Phys. Chem. Lett. 2014, 5, 1973−1977

The Journal of Physical Chemistry Letters

Letter

(10) Li, H. L.; Ibrahim, M.; Agberemi, I.; Kobrak, M. N. The Relationship between Ionic Structure and Viscosity in RoomTemperature Ionic Liquids. J. Chem. Phys. 2008, 129, 124507. (11) Spohr, H. V.; Patey, G. N. Structural and Dynamical Properties of Ionic Liquids: The Influence of Charge Location. J. Chem. Phys. 2009, 130, 104506. (12) Spohr, H. V.; Patey, G. N. Structural and Dynamical Properties of Ionic Liquids: The Influence of Ion Size Disparity. J. Chem. Phys. 2008, 129, 064517. (13) Noda, A.; Hayamizu, K.; Watanabe, M. Pulsed-Gradient SpinEcho 1H and 19F NMR Ionic Diffusion Coefficient, Viscosity, and Ionic Conductivity of Non-Chloroaluminate Room-Temperature Ionic Liquids. J. Phys. Chem. B 2001, 105, 4603−4610. (14) Lopes, J. N. C.; Padua, A. A. H. Molecular Force Field for Ionic Liquids III: Imidazolium, Pyridinium, and Phosphonium Cations; Chloride, Bromide, and Dicyanamide Anions. J. Phys. Chem. B 2006, 110, 19586−19592. (15) Tokuda, H.; Hayamizu, K.; Ishii, K.; Abu Bin Hasan Susan, M.; Watanabe, M. Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 1. Variation of Anionic Species. J. Phys. Chem. B 2004, 108, 16593−16600. (16) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 2. Variation of Alkyl Chain Length in Imidazolium Cation. J. Phys. Chem. B 2005, 109, 6103−6110. (17) Tokuda, H.; Ishii, K.; Susan, M. A. B. H.; Tsuzuki, S.; Hayamizu, K.; Watanabe, M. Physicochemical Properties and Structures of RoomTemperature Ionic Liquids. 3. Variation of Cationic Structures. J. Phys. Chem. B 2006, 110, 2833−2839. (18) Borodin, O.; Gorecki, W.; Smith, G. D.; Armand, M. Molecular Dynamics Simulation and Pulsed-Field Gradient NMR Studies of Bis(fluorosulfonyl)imide (FSI) and Bis[(trifluoromethyl)sulfonyl]imide (TFSI)-Based Ionic Liquids. J. Phys. Chem. B 2010, 114, 6786−6798. (19) Znamenskiy, V.; Kobrak, M. N. Molecular Dynamics Study of Polarity in Room-Temperature Ionic Liquids. J. Phys. Chem. B 2004, 108, 1072−1079. (20) Borodin, O. Relation between Heat of Vaporization, Ion Transport, Molar Volume, and Cation−Anion Binding Energy for Ionic Liquids. J. Phys. Chem. B 2009, 113, 12353−12357. (21) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113, 11463−11478. (22) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Gromacs  A Message-Passing Parallel Molecular-Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43. (23) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (24) Lindahl, E.; Hess, B.; van der Spoel, D. Gromacs 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306. (25) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (26) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182. (27) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089.

All simulations were conducted in the constant temperature, constant pressure ensemble (N, P, T). The systems were maintained at the requested temperature using the velocity rescaling Bussi−Parrinello thermostat.25 The time constant of 1.0 ps was applied to ensure weak temperature coupling. The Parrinello−Rahman barostat26 was responsible for constant atmospheric pressure with the time constant of 2.0 ps and the compressibility constant of 4.5 × 10−5 bar−1. The cutoff distance of 1.2 nm for the Lennard-Jones potential was employed in conjunction with shifted force modification between 1.1 and 1.2 nm. The electrostatic interactions were computed using the direct pairwise Coulomb potential at separations smaller than 1.5 nm and using the particle mesh Ewald scheme27 for all distances beyond the cutoff. The neighbor list was updated every 0.01 ps within 1.5 nm. Periodic boundary conditions in the three Cartesian directions were simulated. All computations were conducted using 12-core nodes of the Exciton cluster (University of Rochester). Domain decomposition23 was used to distribute computational load over the nodes.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] or [email protected]. Tel.: +45 6550 9198. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS V.V.C. thanks M. Watanabe (Yokohama National University) for insightful discussion regarding NMR measurements of diffusion constants and their thermal dependences. The research was supported by Grant CHE-1300118 from the U.S. National Science Foundation. MEMPHYS is the Danish National Center of Excellence for Biomembrane Physics. The Center is supported by the Danish National Research Foundation.

■

REFERENCES

(1) Lin, R. Y.; Taberna, P. L.; Fantini, S.; Presser, V.; Perez, C. R.; Malbosc, F.; Rupesinghe, N. L.; Teo, K. B. K.; Gogotsi, Y.; Simon, P. Capacitive Energy Storage from −50 to 100 °C Using an Ionic Liquid Electrolyte. J. Phys. Chem. Lett. 2011, 2, 2396−2401. (2) Maginn, E. J. What to Do with CO2. J. Phys. Chem. Lett. 2010, 1, 3478−3479. (3) Urahata, S. M.; Ribeiro, M. C. C. Unraveling Dynamical Heterogeneity in the Ionic Liquid 1-Butyl-3-methylimidazolium Chloride. J. Phys. Chem. Lett. 2010, 1, 1738−1742. (4) Wishart, J. F. Ionic Liquids and Ionizing Radiation: Reactivity of Highly Energetic Species. J. Phys. Chem. Lett. 2010, 1, 3225−3231. (5) Chaban, V. V.; Prezhdo, O. V. Ionic and Molecular Liquids: Working Together for Robust Engineering. J. Phys. Chem. Lett. 2013, 4, 1423−1431. (6) Maciel, C.; Fileti, E. E. Molecular Interactions between Fullerene C-60 and Ionic Liquids. Chem. Phys. Lett. 2013, 568, 75−79. (7) Chaban, V. Polarizability versus Mobility: Atomistic Force Field for Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055−16062. (8) Chaban, V. V.; Voroshylova, I. V.; Kalugin, O. N. A New Force Field Model for the Simulation of Transport Properties of Imidazolium-Based Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 7910−7920. (9) Davey, S. ‘Ecosystem in a Box’ Strategy to Determine the Effects of Structure on Toxicity  Testing the Toxicity of Ionic Liquids. Phys. Chem. Chem. Phys. 2007, 9, C75−C75. 1977

dx.doi.org/10.1021/jz5007127 | J. Phys. Chem. Lett. 2014, 5, 1973−1977