Ab Initio Force Fields for Imidazolium-Based Ionic Liquids - The

Subscriber access provided by UNIV OF NEWCASTLE

Article

Ab Initio Force Fields for Imidazolium-Based Ionic Liquids Jesse Gatten McDaniel, Eunsong Choi, Chang Yun Son, Jordan R. Schmidt, and Arun Yethiraj J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b05328 • Publication Date (Web): 28 Jun 2016 Downloaded from http://pubs.acs.org on July 3, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Ab Initio Force Fields for Imidazolium-Based Ionic Liquids. Jesse G. McDaniel,† Eunsong Choi,‡ Chang Yun Son,† J. R. Schmidt,† and Arun Yethiraj∗,† Department of Chemistry, University of Wisconsin, Madison, Wisconsin, 53706, and Department of Physics,University of Wisconsin, Madison, Wisconsin, 53706 E-mail: [email protected]

∗ To

whom correspondence should be addressed of Chemistry, University of Wisconsin, Madison, Wisconsin, 53706 ‡ Department of Physics,University of Wisconsin, Madison, Wisconsin, 53706 † Department

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We develop ab initio force fields for alkylimidazolium based ionic liquids (ILs) that predict the density, heats of vaporization, diffusion, and conductivity that are in semi-quantitative agreement with experimental data. These predictions are useful in light of the scarcity of, and sometimes inconsistency in, experimental heats of vaporization and diffusion coefficients. We illuminate physical trends in the liquid cohesive energy with cation chain length and anion. These trends are different than those based on the experimental heats of vaporization. Molecular dynamics prediction of the room-temperature dynamics of such ILs is more difficult than is generally realized in the literature, due to large statistical uncertainties and sensitivity to subtle force field details. We believe that our developed force fields will be useful for correctly determining the physics responsible for the structure/property relationships in neat ILs.

1

Introduction

Room-temperature ionic liquids (ILs), normally composed of organic cations and anions, are potentially important solvents and electrolytes for use in a variety of different applications. 1 Due to the numerous combinations of different organic ions and functional variants thereof, ILs are potential “designer” solvents, as their properties can be adjusted by employing different constituent molecules. 2 It is therefore of practical importance to develop a molecular-level understanding of the properties of ILs in order to more effectively design and utilize them for present and future applications. 3 This has motivated a large number of theoretical investigations of these systems, usually employing either quantum-chemical calculations to analyze the electronic structure of ion pairs or clusters, or molecular dynamics (MD) simulations to study the properties of the bulk liquids. 4 A major purpose of these theoretical studies is to develop a thorough understanding of how the individual ion properties such as their size, shape, flexibility, charge delocalization, polarizability, and potential for charge transfer, effect the bulk liquid properties, namely the density, structure, cohesive energy, vapor pressure, viscosity, ion diffusion and conductivity; such insight would en2 ACS Paragon Plus Environment

Page 2 of 41

Page 3 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


able molecular-level design and tuning of new ILs. With the exception of some ab initio molecular dynamics (AIMD) studies, 5–8 almost all MD simulations rely on force fields to describe the intraand inter-molecular interactions among the ions. In line with the above comments, it is thus highly important for these force fields to not only be accurate, but also physically-meaningful, so that the influence of different physical interactions on the bulk properties can be correctly assessed from the simulation. There have been many previous force fields for ionic liquids but this remains an area of active research, i.e., there is no accepted standard force field. Non-polarizable force fields include those based on the OPLS-AA 9–13 and AMBER 14–20 framework. Other force fields include those of Chaban and co-workers 21–23 and Maginn and coworkers. 4,24–29 Borodin and coworkers have developed both explicitly polarizable atomistic force fields 30,31 as well as non-polarizable united-atom (UA) force fields, 32,33 for a variety of ILs, notably yielding good agreement with experimental enthalpies of vaporization and dynamic properties. 34 Wang and coworkers 35–38 have developed IL force fields using their multiscale coarse-graining approach. Yan and coworkers developed polarizable models for EMIM+ /NO− 3 , and present a thorough discussion of the influence of polarization on the structure and dynamics of this IL. 39–42 A number of other force fields have been empirically refined based on experimental data. 43–51 Notable exceptions to empirical parametrization include force matching or parameter refinement based on AIMD simulations. 52–54 For a further discussion of IL force field development, we refer the interested reader to recent reviews. 55,56 Most IL force fields have thus relied, to varying extent, on empirical parametrization. However, the development of accurate ab initio force fields for ILs is desirable for a number of reasons. First, experimental determination of certain properties (e.g. diffusion coefficients, heats of vaporization) is difficult, and this data is not available for all ILs; additionally, for the ILs that have been characterized, there is sometimes significant discrepancy among different experimental reports of these properties. Second, the magnitude of the cohesive energy, which is a physically intuitive quantity that has important implications for the IL’s properties, cannot be directly determined from experimentally observed heats of vaporization due to gas phase ion pair binding. We show that trends



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in the cohesive energy do not always follow trends in the heat of vaporization, meaning that accurate predictions of the cohesive energy from simulations are essential for correctly interpreting the physical properties of ILs.. Finally, predicting the properties of complex systems such as ILs using entirely first-principles force fields is, in its own right, an important challenge for the field of theoretical chemistry. In this work, we extend our previous efforts 57,58 to develop an entirely ab initio force field for ILs composed of imidazolium-based cations, EMIM+ (1-Ethyl-3-methylimidazolium), BMIM+ (1-Butyl-3-methylimidazolium), and C6 MIM+ (1-Hexyl-3-methylimidazolium), and anions BF− 4 and PF− 6 . Our resulting force field predicts the densities, heats of vaporization, diffusion coefficients, and conductivities of the ILs mostly within semi-quantitative agreement with experiment, at multiple temperatures, with no adjustable parameters. New physical insight is presented for the trends in the cohesive energy of the ILs, which cannot be correctly inferred based on previous characterization of the heats of vaporization. Additionally we carefully analyze the dynamics of these ILs at room temperature, characterizing the large statistical uncertainties and sensitivity to subtle force field details (e.g. alkyl chain flexibility of EMIM+ ).

2

Force Field Development

2.1

Evaluating Polarization and Charge Transfer

Before discussing the details of our force field development, we examine the presence (or absence) of charge transfer in these ILs, as previous research 53,54,59 has suggested that this effect may be important in certain ILs, which would have important implications for force field development. To examine this, we calculate electron density difference isosurfaces for interacting gas phase ion pairs at several local minima. − − − + + + For each type of ion pair (EMIM+ /BF− 4 , EMIM /PF6 , BMIM /BF4 , BMIM /PF6 , − + C6 MIM+ /BF− 4 , C6 MIM /PF6 ) ∼ 100 gas phase configurations are generated by randomly

sampling an even distribution of solid angles describing the cation/anion center-of-mass displace4 ACS Paragon Plus Environment

Page 4 of 41

Page 5 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ment (~rCOM ) vector, while keeping the internal molecular axes fixed. rCOM is then chosen such that the closest contact between ions occurs at ∼ 0.8 van der Waals contact surface. All configurations are then relaxed to the nearest local minima, at the AVDZ/PBE0 level of theory. There is much redundancy in the optimized geometries, and we define unique local minima by those separated by > 2.5 kJ/mol in energy. For each ion pair type we choose the 5-10 most energetically favorable unique local minima, and further optimize their geometry at the AVTZ/PBE0 level of theory. Electron density difference isosurfaces are then computed by subtracting the electron density of isolated ions from the density of the ion pair, employing a dimer-centered basis set for all calculations. Note that any charge transfer from the anion to cation would result in significant differences between these computed densities. The electron density difference isosurfaces (plotted for electron density values ±0.003 a.u.) for four representative ion pair configurations are shown in Figure 1. We find no significant net depletion (enhancement) of the total electron density of the anion (cation), indicating little to no charge transfer for these ion pairs. These observations are consistent in all other (not shown) ion pair configurations. We additionally compute Bader charges for all ion pairs, and find net charges of ± 1.00 ( ± 0.01) electrons for all cations and anions. The electron density difference isosurfaces do, however, show important polarization effects: The cation exhibits significant polarization on both the imidazolium ring as well as the alkyl groups, notably enhancing the dipole moment of the C-H groups closest to the anion; the polarization of the anion is in general less significant. These findings guide our subsequent physically-motivated force field development: Our force field explicitly incorporates polarization, and utilizes full (± 1) ion charges determined to reproduce the electrostatic fields of each isolated ion.

2.2

Force Field Functional Form

We employ the following functional form for our force field: improper proper Etot = Ebond + Eangle + Edihedral + Edihedral + Enonbond


(1)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Electron density difference maps of select ion pairs. Isosurface values of 0.003 a.u. (yellow) and -0.003 a.u. (red) correspond to enhancement and depletion respectively of the ion pair electron density relative to the isolated ions.


Page 6 of 41

Page 7 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The bond, angle, and improper dihedral terms are modeled using standard harmonic functional forms. The proper dihedral energy is modeled using the Ryckaert-Belleman functional form: proper = Edihedral

5

∑ ∑ Cni jkl (cos(φi jkl − 180◦))n

(2)

{φi jkl } n=0

where φi jkl is the angle between the planes defined by atoms i jk and jkl (φi jkl = 0 corresponds to cis configuration). For the non-bonded energy we employ the functional form, 60

Enonbond

ij qi q j Cn tot =∑ + ∑ Ai j exp(−Bi j ri j ) − ∑ fn(Bi j , ri j ) rn +Ushell ij i, j ri j i, j n=6,8,10,12

(3)

Where the summation ∑i, j runs over all atom pairs on different molecules, as well as all atom pairs on the same molecule separated by four or more bonds. 1-4 intra-molecular interactions employ the above interactions scaled by 0.5. The employed Tang-Toennies damping function is defined as, fn (λ , r) = 1 − e−λ r

n

(λ r)m m=0 m!

∑

(4)

The Drude oscillator energy, Ushell , which explicitly accounts for polarization, is calculated using standard procedures. 61 All intra-molecular shell-shell interactions are explicitly considered; those at 1-4 or closer are screened using a Thole-screening function, 62 giving an interaction, Ui j = T (ri j )

qi q j , ri j

T (ri j ) = 1 − 1 +

pri j −pri j /(αi α j )1/6 e 2(αi α j )1/6

(5)

q2shell where αi = is the polarizability of atom i. We set the spring constant k = 0.1 a.u. and k Thole parameter p = 2 as done previously. 60 At each time step, the Drude oscillator positions are self-consistently optimized to minimize the system energy, until the maximum force on each shell is less than 0.1 kJ/mol/nm. All force field parameters, with the exception of bond, angle, and improper dihedral potentials,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

are determined entirely based on ab initio calculations, as described in Sections 2.3 and 2.4

2.3

Intermolecular Terms

Intermolecular interaction parameters are fit based on intermolecular symmetry-adapted perturbation theory (SAPT) calculations, using the methodology of Schmidt and coworkers. 60,63 Briefly, atomic point charges are fit to a distributed-multipole analysis of the monomer electron densities. 64,65 Atomic polarizabilites and dispersion coefficients are fit to molecular linear response calculations based on an iterative extension of the Williams-Stone methodology. 60,66 These linear response calculations are computed with the Camcasp software, 67,68 at the coupled-perturbed Kohn-Sham level of theory. The short-range exponential terms accounting for exchange-repulsion and charge penetration effects are fit to homo-molecular DFT-SAPT 69–71 calculations, conducted using the Molpro software. 72 Charges are explicitly fit for every ion, while the remaining non− bonded parameters are explicitly fit only for the imidazolium ring, BF− 4 , and PF6 , and the

additional alkyl group parameters are taken from previous work. 60 We note the BMIM+ /BF− 4 non-bonded parameters are identical to those employed in our previous work. 57 The fidelity of the force field’s description of intermolecular interactions is examined through explicit comparison − 57 57 + to SAPT calculations of ion pairs: BF− 4 /BF4 (ref ), imidazolium/imidazolium (ref ), BMIM − − − − 57 + + /BF− 4 (ref ), PF6 /PF6 (Figure S5), EMIM /PF6 (Figure S6), and EMIM /BF4 (Figure 2).

We note that cation/anion interactions are described with similarly high fidelity for BF− 4 (Figure 2) and PF− 6 (Figure S6) anions. All intermolecular force field parameters are given in the supporting information.


Page 8 of 41

Page 9 of 41

-150

140

a)

b)

120

-200

EFF (kJ/mol)

EFF (kJ/mol)

100 80 60

-250

-300

40 -350

20 0

0

20

40

80 60 100 ESAPT (kJ/mol)

120

-400 -400

140

0

c)

-300 -250 ESAPT (kJ/mol)

-200

-150

-40


-10

0

-200

-150

d)

EFF (kJ/mol)

-10

-20

-30

-40

-50 -50

-350

0

-10

EFF (kJ/mol)

-20

-30

-40

-40

-30 -20 ESAPT (kJ/mol) -150

-10

0

-50 -50

e) -200

EFF (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-250

-300

-350

-400 -400

-350


Figure 2: Force field predicted energies (y-axis) compared to SAPT energies (x-axis) for EMIM+ /BF− 4 dimers pulled from MD simulation: a) exchange, b) electrostatic, c) induction+dhf, d) dispersion, and e) total interaction energy.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.4

Page 10 of 41

Intramolecular Terms

There are many intra-molecular interaction terms in the IL force field, and full ab initio parametrization of all of these terms is a tedious and laborious task. Therefore, we decided to focus our parametrization efforts on the proper dihedral potentials of the cation alkyl groups, as these are the major source of conformational flexibility. Intramolecular bond, angle, and improper dihedral parameters for the cations are therefore taken from Sambasivarao and Acevedo. 13 For the anions, 14 bond and angle force field parameters are taken from de Andrade et al. 18 for BF− 4 and Liu et al.

for PF− 6 . Preliminary results suggest that the simulation results are largely insensitive to the exact values of bond, angle, and improper force constants, and employing literature values for these terms is sufficient. Specifically, we test the sensitivity of the simulation results to the magnitude of the cation improper dihedral force constant(s) enforcing planar ring configurations; this is motivated by preliminary simulations, in which we noticed that the imidazolium rings exhibited significant deviations from planarity. We perform test calculations, artificially increasing these force constants by up to two orders of magnitude; while this changes the improper dihedral distributions of the rings, no significant changes to predicted thermodynamic or kinetic properties are observed. We fit dihedral potentials for the alkyl chains of EMIM+ , BMIM+ , and C6 MIM+ using the following procedure. For each cation, all possible dihedral angles of the alkyl chain are manually rotated at 18◦ increments, generating 20n configurations, where ”n” is the number of independent alkyl dihedral angles, including the terminal methyl hydrogen dihedral (n=2, 4, 6 for EMIM+ , BMIM+ , C6 MIM+ respectively). For BMIM+ and C6 MIM+ , consideration of all these configurations is intractable and thus we employ Monte Carlo sampling (based on non-bonded force field component energies) to omit high-energy, sterically repulsive configurations. For C6 MIM+ , the configurational space is still too large, and we further randomly parse the remaining configurations down to 30,000. Ab initio calculations at the PBE0/AVTZ level are conducted to compute the relative configurational energies; this enables a direct comparison of the relative dihedral energetics, as the configurations employ identical bonds and angles. The dihedral potentials are then fit to the 10 ACS Paragon Plus Environment

Page 11 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


residual energy difference between the relative ab initio energies and the force field non-bonded interactions. For every -CH2- group in the alkyl chain, dihedrals are only fit for heavy atoms, and R-R-C-H dihedrals are set to zero, to ensure that only one independent dihedral type is fit to each degree of freedom for the rigid configuration (an exception is for the terminal R-R-CH3 group, in which case the hydrogen atoms exhibit non-zero dihedral terms). The parametrized force field representation, ab initio, and absolute difference for the (relative) intra-molecular potential energy surface of EMIM+ as a function of the two alkyl dihedral angles are shown in Figure 3 a), b), and c) respectively (additionally, a different representation of this fit is given in Figure S2). We do not explicitly show the dihedral fits for BMIM+ and C6 MIM+ , due to the difficulty of clearly presenting this data in the enlarged phase space; these other fits are of similar quality to that shown for EMIM+ . We note that it is necessary to employ different dihedral potentials for the three cations because the charges on the alkyl chains are different for each cation, giving different non-bonded contributions to the configurational energies. All ab initio force field parameters are given in the supporting information.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3: Relative potential energy surface of EMIM+ as a function of the two alkyl dihedral angles computed with a) force field, b) ab initio, and c) difference between the two. We note that exact 3-fold symmetry along the N-C-C-H dihedral angle is broken due to asymmetry in the methyl group of the rigid configuration employed for fitting.


Page 12 of 41

Page 13 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3

Results & Discussion

− − + + Condensed phase properties of the six ILs ([EMIM+ ][BF− 4 ] , [EMIM ][PF6 ] , [BMIM ][BF4 ] − − + + , [BMIM+ ][PF− 6 ] , [C6 MIM ][BF4 ] , [C6 MIM ][PF6 ] ) are computed at 298 K and 353 K

from molecular dynamics (MD) simulations using the GROMACS 4.6 software 73 (an exception 74 and properties are only computed for is [EMIM+ ][PF− 6 ] , which is a solid at room temperature,

353 K). Equilibration is carried out in the NPT ensemble for 2 ns, followed by 50 ns production runs in the NVT ensemble, using a 2 fs time step; all simulations employ 200 ion pairs (unless otherwise noted). The Berendsen barostat 75 and Nose-Hoover thermostat 76,77 are used for pressure and temperature coupling, respectively. Particle mesh Ewald 78 is used for electrostatics, and van der Waals (VDW) interactions are shifted to zero at 1.4 nm (this introduces negligible error, vide infra). All GROMACs input files implementing our force field are included as supplementary material. As the dynamics of ILs at room temperature are extremely slow, prediction of these properties is very difficult, and subject to significant uncertainty. We therefore evaluate the statistical uncertainty of the computed dynamical properties (at 298 K) by conducting an additional eight independent 50 ns simulations for [BMIM+ ][BF− 4 ] at 298 K, starting from different uncorrelated snapshots of the initial NPT equilibration simulation. Additionally, we conduct a 50 ns simulation for a much larger [BMIM+ ][BF− 4 ] system composed of 1600 ion pairs; this allows for enhanced averaging and observation of any finite-size effects. We use the statistical uncertainty obtained from these additional simulations as an estimate of the general uncertainty in the 298 K dynamical properties for all IL systems studied. To the best of our knowledge, this is currently the longest total simulation time (∼ 500 ns) reported in the literature for estimating the uncertainty in IL dynamic properties, using an all-atom, polarizable force field.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.1

Page 14 of 41

Density

The computed IL densities are compared with experimental values in Table 1. For the BF− 4 based ILs, the predicted densities are generally ∼ 2 % lower than experimental values, while agreement − + with experiment for PF− 6 ILs is within 1 % (agreement is slightly worse for [EMIM ][BF4 ] and

[BMIM+ ][BF− 4 ] at 353 K). Based on the ab initio calculations, it is difficult to explain the slightly − poorer predictions for the density of ILs composed of BF− 4 (compared to PF6 ) anions. Despite

this, the overall agreement with experiment for the IL densities should be considered quite good, as there are no adjustable parameters. It is interesting to compare these ab initio predictions to similar studies of small-molecule organic liquids. Such organic liquids have comparatively low-cohesive energy density, and as such, long-range corrections to the VDW interactions as well as three-body dispersion significantly effect the predicted density. 79 However for ILs, the cohesive energy is much higher due to the electrostatics, and as such said physics may be neglected without incurring significant error. This is verified by additional simulations employing 1.0 and 1.8 nm cutoffs for the VDWs interactions, in which we observe respective changes in the density of ∼ 1 % and 0 % for [BMIM+ ][BF− 4] compared to the simulations employing a 1.4 nm cutoff. Table 1: Predicted and experimental densities for ILs at 298 K and 353 K. The experimental values are given as averages of available reference values, after omitting outliers. The range of experimental values and the corresponding references are given in brackets. Uncertainties in the simulation results are on the order of the last reported digit. a [EMIM+ ][PF− 6 ] is a solid at room temperature. IL MD − + [EMIM ][BF4 ] 1255 [EMIM+ ][PF− N/Aa 6] − + [BMIM ][BF4 ] 1180 [BMIM+ ][PF− 1378 6] − + [C6 MIM ][BF4 ] 1121 [C6 MIM+ ][PF− 6 ] 1299

ρ (kg/m3 ) 298 K Exp. 1280 [1240-1337 74,80–84 ] N/Aa 1203 [1202-1210 83,85–87 ] 1370 [1368-1371 85,86 ] 1146 [1145-1148 74,89 ] 1293 [1292-1295 74,87,89 ]

MD 1204 1432 1135 1327 1086 1259


ρ (kg/m3 ) 353 K Exp. 1246 [1239-1253 74,84 ] 1422 74 1162 [1161-1164 85,87,88 ] 1323 85 1108 74 1249 [1248-1250 74,87 ]

Page 15 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2

Heat of Vaporization and Cohesive Energy

For ILs, the cohesive energy of the liquid cannot be directly inferred from experimental heats of vaporization, due to gas phase ion pair binding, and trends in the former do not necessarily correspond to trends in the latter. The cohesive energy is a measure of how strongly the distinct components of a liquid are bound together, which has direct implication for all other liquid properties. The cohesive energy is most naturally defined as (normalized per ion pair):

Ecohesive (T ) =

1 Nion−pair

Eliquid (T ) − Ecation (T ) − Eanion (T )

(6)

where Ecation and Eanion are the energies of the isolated ions. This cohesive energy is distinct from the heat of vaporization, which is given by:

∆Hvap (T ) =

1 Nion−pair

Evapor (T ) − Eliquid (T ) + RT

(7)

Unlike common organic liquids, Evapor is an important, non-trivial contribution for ILs, which significantly reduces ∆Hvap relative to the liquid cohesive energy. This is because the ionic liquid vapor predominantly consists of neutral ion pairs. 90,91 Therefore, besides the “RT” term, Ecohesive and ∆Hvap differ by the gas phase ion pair binding energy, given by:

Eion−pair (T ) =

1 Nion−pair

Evapor (T ) − Ecation (T ) − Eanion (T )

(8)

The approximation that the vapor phase consists entirely of single, neutral ion pairs will be discussed later. While the above equations could be simplified by defining an energy reference Ecation = Eanion = 0, this is not a convenient choice, as these quantities are temperature dependent, due to intra-molecular interactions. The computed values of Eliquid , Evapor , Ecation , and Eanion for all ILs are given in the Supporting Information. Figure 4, 5, and 6 depict the cohesive energy, enthalpy of vaporization, and gas phase ion pair binding energy for the ILs (note we plot the magnitude of the cohesive and ion pair energies, the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

actual values are negative). There are pronounced trends in both the cohesive energy and gas phase ion pair binding energy that are rationalized based on chemical intuition. For a given anion, the magnitude of the cohesive energy decreases with increasing alkyl chain length of the cation; this is due to the corresponding decrease in total charge density of the system with the addition of “neutral” CH2 or CH3 groups. Also, the magnitude of the cohesive energy is lower for PF− 6 compared to BF− 4 ILs, again due to the change in total system charge density, this time because of the anion size. For the gas phase ion pair binding energies, the same trends hold; smaller ions, with more local charge density, have stronger binding energies, and thus the magnitude of the ion pair binding energies decreases with increasing cation alkyl chain length, and decreases with increasing − anion size, from BF− 4 to PF6 .

There is little (or even incorrect) physical interpretation from the heat of vaporization alone. For increasing cation alkyl chain length, the heat of vaporization also increases. This could be incorrectly interpreted to mean that interactions are stronger with increasing alkyl chain length. Rather the trend is an artifact of the subtle cancellation between liquid and gas phase ion energetics, and merely reflects that the magnitude of the gas phase ion binding energies decrease comparatively more than the magnitude of the liquid cohesive energy with increasing alkyl chain length. Due to this subtle cancellation, the trend in the heat of vaporization is not robust, and is − only seen for BF− 4 and not PF6 ILs (Figure 5). Similarly, the enthalpy of vaporization is higher − for PF− 6 compared to BF4 ILs. This observation could once again be misinterpreted to mean that

the PF− 6 ILs have stronger ionic interactions. Instead, this is due to the fact that there are larger changes in the gas phase ion pair energies than in the liquid cohesive energy with the two different anions. The fact that trends in the heat of vaporization may be easily misinterpreted is an important point, as considerable effort has been directed at rationalizing such trends. 92–94 Physical trends are much more apparent and interpretable when the cohesive energy of the liquid, rather than heat of vaporization, is considered, and this quantity is only directly obtained from computer simulations. Indeed, Borodin 95 has previously shown that the cohesive energy exhibited better correlation than the heat of vaporization with IL transport properties.


Page 16 of 41

Page 17 of 41

We finally note the different temperature dependence of the IL cohesive energy as a function of cation size (Figure 4). There is increasing temperature dependence of the cohesive energy with increasing alkyl chain length of the cation, e.g. the cohesive energy of [C6 MIM+ ][BF− 4 ] changes more with temperature than [EMIM+ ][BF− 4 ] ; this does not directly correlate with the temperature dependence of the density (Table 1). This may potentially be explained by spatial heterogeneity due to hydrophobic aggregation, which is more pronounced with increasing alkyl chain length; spatial heterogeneity has been generally reported for ILs, 35,96–98 and deserves a more extensive analysis that is beyond the scope of the present paper. 500

|Ecohesive| (kJ/mol)

490

298 K 353 K

BF4

PF6

480

470

460

450

[EMIM] [BMIM] [C6MIM] [EMIM] [BMIM] [C6MIM]

Figure 4: Magnitude of the cohesive energy (per ion pair) of ILs at 298 K and 353 K. The actual cohesive energy is negative.

140 135

∆Hvap (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


298 K 353 K

BF4

PF6

130 125 120 115 110


Figure 5: Enthalpy of vaporization of ILs at 298 K and 353 K. Data is the same as that given in Table 2.



380 370

|Eion-pair| (kJ/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 41

298 K BF4 353 K

PF6

360 350 340 330 320


Figure 6: Magnitude of gas phase ion pair binding energy of ILs at 298 K and 353 K. The actual ion pair binding energy is negative. Table 2: Predicted and experimental heats of vaporization for ILs at 298 K and 353 K. Uncertainties in the simulation results are ∼ 2-3 kJ/mol, which mainly come from prediction of the gas phase ion pair energies (see Supporting Information). We only list experimental data for one temperature, to avoid redundancy, as the data is extrapolated from higher temperature measurements. a [EMIM+ ][PF− 6 ] is a solid at room temperature. ∆Hvap (kJ/mol) 298 K IL MD Exp. − + [EMIM ][BF4 ] 121 139, 99 136 92 [EMIM+ ][PF− N/Aa N/Aa 6] − + 100 [BMIM ][BF4 ] 126 128, 141, 99 152 93 [BMIM+ ][PF− 134 155 100 6] [C6 MIM+ ][BF− 132 145 99 4] [C6 MIM+ ][PF− 137 140 100 6] − [C8 MIM+ ][BF4 ] N/A 122 100 , 162 101 − + [C8 MIM ][PF6 ] N/A 144 100 , 169 101

∆Hvap (kJ/mol) 353 K MD Exp. 117 — 131 138 92 120 — 129 — 126 — 131 — N/A — N/A —

In Table 2 we compare the predicted enthalpy of vaporization of the ILs to “experimental” results, the latter of which warrant further comment. Experimental ∆Hvap data are most widely available for imidazolium based ILs employing Bis(trifluoromethylsulfonyl)imide anions, [Cn MIM+ ][NTf− 2] − ; corresponding ∆Hvap data for BF− 4 and PF6 based ILs are relatively scarce. Verevkin and cowork-

ers 94 have thoroughly compiled and analyzed the numerous available experimentally reported 94 it is clear that there are significant dif∆Hvap for [Cn MIM+ ][NTf− 2 ] . From this compilation,

ferences among reported experimental values for identical ILs. Room temperature ∆Hvap values are determined by extrapolation from higher temperature experimental measurements, requiring 18 ACS Paragon Plus Environment

Page 19 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


knowledge of ∆glC p , the difference between the liquid and vapor phase heat capacity. ∆glC p must either be indirectly determined through experimental measurements of ∆Hvap at different temperatures, or determined through theoretical methods. It is suggested that the use of inaccurate values (e.g. a constant value for a wide range of ILs) of ∆glC p partially accounts for discrepancies among reported experimental values. 94 − + Experimental ∆Hvap data for [EMIM+ ][BF− 4 ] and [EMIM ][PF6 ] is given in Table 2 from

quartz crystal microbalance experiments of Zaitsau and coworkers. 92 We also report temperature programmed desorption measurement data from Deyko et al. 93 for [BMIM+ ][BF− 4 ] and for longer − + 101 all of these works use chain, [C8 MIM+ ][BF− 4 ] and [C8 MIM ][PF6 ] from Armstrong et al.;

an estimate for ∆glC p to extrapolate the experimental measurements to lower temperature. It is important to note that the rest of the reported experimental data 99,100 are not direct measurements; rather, an empirical correlation between ∆Hvap and IL surface tension was determined 100 based on a different IL series, namely [Cn MIM+ ][NTf− 2 ] , and this same correlation was employed to − estimate ∆Hvap values for the BF− 4 and PF6 ILs based on their surface tension values. Such values

cannot therefore be considered “exact”, as the transferability of this empirical correlation is an approximate assumption. Due to the inconsistency of the experimental data as well as the approximate extrapolation/correlation schemes, it is difficult to conclusively comment on the exact quantitative accuracy of our predictions for ∆Hvap ; any deviation of our predictions is mostly within the experimental scatter. However, based on overall trends of the entire ∆Hvap data, it is possible that our predictions for ∆Hvap are slightly too low, and this may reflect our approximation that the vaporized species are entirely composed of single, neutral ion pairs. While it is largely agreed that single ion pairs constitute the majority of vaporized species, the presence of other species would likely lead to higher heats of vaporization, 102 and some research indeed suggests a minority concentration of other species in the vapor. 90,91 One would additionally expect that any fraction of minority species in the vapor is temperature dependent, and thus dependent on the particular experimental technique used to measure ∆Hvap . Due to the entirety of the mentioned complications, we believe that our ab



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

initio predictions for the trends in both the enthalpy of vaporization and liquid cohesive energy are especially valuable for deriving correct insight about the structure/property relations of these quantities.

3.3

Diffusion Coefficients

The uncertainties for predicted dynamical properties of ILs at 298 K are significantly larger than is generally recognized in the literature. This is shown by Figure 7, where we plot both the cation and anion mean squared displacement (MSD) as computed from each of the 10 independent [BMIM+ ][BF− 4 ] trajectories. From the two trajectories bounding the data, one predicts diffusion coefficients (from the slopes) of D+ ∼ 0.008-0.015 (10−5 cm2 /s) and D− ∼ 0.006-0.01 (10−5 cm2 /s) for the cation and anion respectively; we note that differences in trajectories may partially be attributed to slightly different (< 1 %) system densities (Figure S10). The average values (mean and median are nearly identical) from the 10 simulations are D+ =0.011 (10−5 cm2 /s) and D− =0.008 (10−5 cm2 /s). We note that the predicted diffusion constants from the simulation of 1600 ion pairs, D+ =0.011 (10−5 cm2 /s) and D− =0.009 (10−5 cm2 /s), are very close to the average values, indicating minimal finite-size effects, and also that deviations of the independent simulations may be partially statistical in nature. It is important to note that 50 ns (for each trajectory) is a long simulation time relative to previously reported simulations in the IL literature, especially for polarizable force fields, and it is clear that such time scales are not long enough to ensure converged statistics at 298 K. Based on this analysis, we obtain a general estimate of the statistical uncertainty for the 298 K dynamics of ∼ ±30%. It is expected that uncertainties (at 298 K) are even higher for the ILs with slower dynamics, such as [C6 MIM+ ][PF− 6].


Page 20 of 41

Page 21 of 41

5

3

BMIM 2.5

4

BF4

2

3

2

MSD (nm )

2

MSD (nm )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

1.5 1

1

0

0

0.5

10

20

30

40

50

0

0

10

Time (ns)

20

30

40

50

Time (ns)

Figure 7: Mean squared displacement (MSD) for BMIM+ and BF− 4 from each independent 50 ns simulation at 298 K. The MSD from the large, 1600 ion pair system simulation is the bold black line. Table 3: Cation (D+ ) and anion (D− ) diffusion coefficients for ILs at 298 K. Uncertainties in the simulation results are at least ±30%, except for [BMIM+ ][BF− 4 ] where the uncertainties are much lower due to the multiple trajectories. a [EMIM+ ][PF− 6 ] is a solid at room temperature. IL D+ − + [EMIM ][BF4 ] [EMIM+ ][PF− 6] + [BMIM ][BF− 4] − + [BMIM ][PF6 ] [C6 MIM+ ][BF− 4] [C6 MIM+ ][PF− 6]

MD 298 K Exp. 298 K −5 2 −5 2 D− (10 cm /s) D+ (10 cm /s) D− (10−5 cm2 /s) 0.006 0.045, 103 0.05 104 0.035, 103 0.041 104 a a N/A N/A N/Aa 0.008 0.014, 85 0.011, 105 0.014 103 0.013, 85 0.013 103 0.003 0.006, 105 0.005, 105 0.007 85 0.004, 105 0.005 85 0.004 — — 0.001 — —

(10−5 cm2 /s) 0.015 N/Aa 0.011 0.005 0.005 0.003

The predicted diffusion coefficients for the six ILs are compared to experimental values in Table 3 for 298 K and Table 4 for 353 K. We note that many of the experimental values are extrapolated from the parametrized VFT equations, 85,104 and seemingly small uncertainties in the VFT parameters lead to large percentage differences in diffusion coefficients at 298 K. Within the significant uncertainties of both simulation and experiment, there is generally semi-quantitative agreement. The notable exception is for both the cation and anion diffusion coefficients of [EMIM+ ][BF− 4] at 298 K; our simulations significantly underestimate these values compared to experiment (deviations exist also at 353 K, but are not nearly as significant). This inconsistency unfortunately leads to difficulty in interpreting the effect of increasing alkyl chain length on diffusion; based on 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 4: Cation (D+ ) and anion (D− ) diffusion coefficients for ILs at 353 K. The uncertainties are significantly lower than those reported at 298 K (Figure 3). IL D+ [EMIM+ ][BF− ] 4 [EMIM+ ][PF− 6] + [BMIM ][BF− 4] − + [BMIM ][PF6 ] [C6 MIM+ ][BF− 4] − + [C6 MIM ][PF6 ]

MD 353 K Exp. 353 K −5 2 −5 D− (10 cm /s) D+ (10 cm2 /s) D− (10−5 cm2 /s) 0.102 0.2 104 0.18 104 0.035 — — 85 0.093 0.105 0.11 85 0.041 0.07 85 0.056 85 0.054 — — 0.024 — —

(10−5 cm2 /s) 0.155 0.060 0.113 0.059 0.051 0.029

experiment, increasing alkyl chain length dramatically slows down diffusion for both the cation and anion by a factor of ∼ 3 at 298 K and a factor of ∼ 2 at 353 K for [BMIM+ ][BF− 4 ] compared to [EMIM+ ][BF− 4 ] . Our simulations do predict that diffusion of both the cation and anion slows down with increasing alkyl chain length, EMIM+ > BMIM+ > C6 MIM+ , but to a much smaller extent. We therefore attempted to resolve this discrepancy, to better illuminate the physical mechanism(s) dictating diffusion in these ILs. An intuitive explanation for the enhanced dynamics of EMIM+ compared to BMIM+ is that due to the shorter alkyl chain of the former, there is less steric hindrance and greater mobility of the cation. The extent of the mobility enhancement will depend on the rotation barriers of the alkyl chain. Figure S8 shows that the ethyl chain dihedral of EMIM+ , θC−N−C−C , is largely confined to cis configurations, regardless of temperature. The PES in Figure 3, however, suggests that such a dihedral distribution may be a poor representation of the correct ensemble average. This distribution depends on the energetics of in-plane, steric interactions between the terminal methyl group and the imidazolium ring, which in turn is extremely sensitive to the methyl configuration (this interaction is henceforth termed “C6/ring” interaction, see atom labeling in supporting information). While the fit to the ab initio PES is semi-quantitatively accurate (Figure 3, Figure S2), the force field does predict too high barriers for methyl rotation of such eclipsed configurations (e.g. θC−N−C−C = 0). Additionally, the energy of the C6/ring interaction is very sensitive to the exact bond lengths and angles (as evidenced by the asymmetry in Figure 3), and we have not thoroughly investigated the coupling of these degrees of freedom. 22 ACS Paragon Plus Environment

Page 22 of 41

Page 23 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We therefore performed additional simulations to investigate the sensitivity of the [EMIM+ ][BF− 4] diffusion on the dihedral PES. For our first test, we uniformly scaled the (relative) ab initio PES of EMIM+ by a factor of 0.5, and subsequently refit the force field dihedral potentials to the scaled PES (Figure S3). As shown in Figure 8, the new dihedrals for EMIM+ lead to a smoother PES, however the large barriers for methyl rotation at eclipsed configurations (θC−N−C−C = 0◦ , ±180◦ ) still exist, due to the non-bonded interactions which remain the same. [EMIM+ ][BF− 4 ] simulations employing this new dihedral PES give negligible changes in the θC−N−C−C distribution, and correspondingly the predicted dynamics are unchanged. This result verifies that the C6/ring non-bonded interactions (which are unchanged) are a major determiner of the θC−N−C−C distribution. To examine the sensitivity of the dynamics on the energetics of the C6/ring interaction, we tested the limit in which the PES is independent of the C6 methyl group hydrogen configuration. To do this, we projected the 2-D PES of EMIM+ (Figure 3) onto a 1-D PES dependent only on the θC−N−C−C degree of freedom. The projection we employed is

E(θC−N−C−C ) = minθN−C−C−H E(θC−N−C−C , θN−C−C−H )

(9)

(note this PES is similar to that calculated and discussed by Tsuzuki 106 ). This projection, shown in Figure S4, gives a much smoother PES; to fit this new PES we exclude all intramolecular non-bonded interactions. Employing the new E(θC−N−C−C ) potential gives dramatically different θC−N−C−C distributions, with much greater configurational freedom (Figure S9). Because of this, the resulting diffusion coefficients of [EMIM+ ][BF− 4 ] are a factor of ∼ 2-3 times larger (D+ = 0.035 ∗ 10−5 cm2 /s, D− = 0.017 ∗ 10−5 cm2 /s) at 298 K; interestingly, the dynamics is much less effected at 353 K (D+ = 0.19 ∗ 10−5 cm2 /s, D− = 0.15 ∗ 10−5 cm2 /s). The two important conclusions are: 1) The dynamics of [EMIM+ ][BF− 4 ] at 298 K is strongly dependent on the θC−N−C−C distribution, which is significantly effected by the C6/ring interaction; possible errors in the predicted dependence of the C6/ring interaction on the methyl group configuration may explain some of the discrepancy between the predicted and experimental diffusion coefficients. 2)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 41

The change in diffusion at 298 K but not 353 K, suggests that the particular diffusion mechanisms are temperature-dependent. This latter conclusion is consistent with previous findings, 58 and will be analyzed in detail in future work

Figure 8: Potential energy surface (PES) of EMIM+ for force field fitted to full and scaled ab initio PES.

3.4

Conductivity

The conductivity of the ILs is calculated using the following relation 1 t→∞ 6tV kB T

σ = lim

n

∑

(qi [Ri (t) − Ri (0)]) · (q j [R j (t) − R j (0)])

(10)

i, j

Here, the sum runs over all ions i, j, with charges qi and positions Ri . To evaluate this limit, we calculate the slope from 0 < t < 6 ns; longer times exhibit poor statistical averaging. As for the diffusion coefficients, the uncertainty is estimated by calculating conductivities for the 10 independent [BMIM+ ][BF− 4 ] trajectories at 298 K. These calculated conductivities are bounded by 24 ACS Paragon Plus Environment

Page 25 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.20 < σ < 0.37 S/m , with an average value of 0.29 S/m (the 1200 ion pair simulation gives σ =0.26 S/m). Therefore we estimate a lower bound uncertainty of ∼ 30-40 % for the conductivity values at 298 K, with greater uncertainty expected for the ILs with slower dynamics (e.g. [C6 MIM+ ][PF− 6] ). Table 5: Predicted and experimental conductivities for ILs at 298 K and 353 K. Uncertainties in the simulation results are at least 30-40%. a [EMIM+ ][PF− 6 ] is a solid at room temperature. IL MD − + [EMIM ][BF4 ] 0.50 [EMIM+ ][PF− N/Aa 6] [BMIM+ ][BF− 0.29 4] − + [BMIM ][PF6 ] 0.12 [C6 MIM+ ][BF− 0.10 4] − + [C6 MIM ][PF6 ] 0.04

σ (S/m) 298 K Exp. 107 1.31, 1.38, 107 1.53 108 N/Aa 0.35, 85 0.44 109 0.15, 85 0.17 109 0.15, 109 0.12 108 0.06 109,110

MD 5.65 1.40 2.30 0.89 0.81 0.49

σ (S/m) 353 K Exp. 107 3.65, 4.35, 107 5.7 108 — 2.20, 85 2.17 108 1.24 85 1.00 108 0.64 110

The conductivity data at both 298 K and 353 K are compared to experimental values in Table 5. The trends in conductivity are similar to the corresponding trends in diffusion coefficients (Table 3 and Table 4); conductivity decreases with increasing cation alkyl chain length, and for a given − cation, the conductivity decreases going from BF− 4 to PF6 anions. Additionally, the agreement

with experiment is similar; there is semi-quantitative agreement with experiment for all ILs at both temperatures, except for the anomalous case of [EMIM+ ][BF− 4 ] at 298 K (employing the dihedral surface given by Equation 9 has a similar qualitative effect on the conductivity as it did on the diffusion of [EMIM+ ][BF− 4 ] ). We finally note that our current prediction of both the diffusion coefficients and conductivity for 57 using a very similar [BMIM+ ][BF− 4 ] differ somewhat from the values computed by Choi et al.,

force field. This is due to the fact that the earlier work utilized different proper dihedral potentials for the cation alkyl chain, and omitted 1-4 intramolecular non-bonded interactions (instead of scaling by 0.5), resulting in slightly faster dynamics.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4

Conclusion

Due to their complexity, it is very difficult to achieve a complete understanding of IL structure/property relationships from experiments alone. Therefore, the demonstrated accuracy of our entirely first principles force field is significant, as corresponding simulations enable direct analysis of the essential physics determining the IL properties. As we have shown, this is important for both correctly interpreting existing experimental data, as well as predicting properties that are not directly measurable by experiment (liquid cohesive energy), or investigating uncharacterized ILs. We have illustrated the importance of utilizing a physically-motivated force field for correctly determining structure/property relationships of ILs. We have shown that trends in experimentally measured heats of vaporization have little physical interpretation as they result from subtle cancellations in liquid and vapor phase energetics. Rather, physical trends are directly observed in the liquid phase cohesive energy, which may be best predicted from simulation, utilizing accurate force fields. Additionally, we suggest that conformational flexibility of the ethyl chain of EMIM+ is mechanistically important for the dynamics at 298 K, but is much less important for the dynamics at higher temperatures (353 K). Quantitative prediction of the dynamics of [EMIM+ ][BF− 4] at 298 K seemingly requires a very accurate description of in-plane, C6/ring steric interactions, and the sensitive dependence on the methyl group configuration. This analysis may have been obscured using empirically parametrized force fields, as force fields can predict similar condensed phase properties, while employing very different physics. 58 Lastly, we have highlighted several important issues that are general to the theoretical study of ILs. First, while our force fields incorporate many-body polarization, they explicitly neglect all other non-additive three-body interactions such as three-body exchange and dispersion. In spite of this, they are able to predict IL densities with quantitative accuracy; this would not be the case for simple organic liquids, where such interactions make significant contributions to the internal pressure. 79 For ILs, such interactions are obscured by the very large, electrostatically dominated cohesive energy. Additionally, we have extensively characterized the uncertainty in predicting dynamic properties of ILs at 298K, and we show that simulations of hundreds of nanoseconds are 26 ACS Paragon Plus Environment

Page 26 of 41

Page 27 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


necessary to achieve truly converged quantities; this is notably longer than typical simulation times reported in the literature. Finally, our characterization of the electron density of interacting gas phase, ion pairs reveals no noticeable charge transfer between the anion and cation. We therefore suggest caution in employing reduced/scaled charge models for ILs in general, as such physics may only be physically-motivated for specific ions (e.g. chloride 59 ).

5

Supporting information

All force field parameters are given in the supporting information. Additionally, all relevant GROMACs input files implementing our force fields are included as supplementary material. Additional information given in the supporting information includes dihedral force field fits for EMIM+ , PF− 6 − + / PF− 6 and EMIM / PF6 SAPT fits, EMIM dihedral distributions, equilibration verification, liquid

phase energies, ion pair energies, and single ion energies, and diffusion coefficients as a function of density for the 10 different [BMIM+ ][BF− 4 ] simulations.

6

Acknowledgement

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1111835. This work was partially supported by Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy, under Award DE-FG02-09ER16059. Computational resources were provided by the Center for High Throughput Computing (CHTC) at the University of Wisconsin. The CHTC is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the National Science Foundation, and is an active member of the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy’s Office of Science.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Galinski, M.; Lewandowski, A.; Stepniak, I. Ionic Liquids as Electrolytes. Electrochim. Acta 2006, 51, 5567 – 5580. (2) Kobrak, M. N. Adv. Chem. Phys.; John Wiley & Sons, Inc., 2008; pp 85–138. (3) Castner, E. W.; Wishart, J. F. Spotlight on Ionic Liquids. J. Chem. Phys. 2010, 132, 120901. (4) Maginn, E. J. Molecular Simulation of Ionic Liquids: Current Status and Future Opportunities. J. Phys. Condens. Matter 2009, 21, 373101. (5) Del Popolo, M. G.; Lynden-Bell, R. M.; Kohanoff, J. Ab Initio Molecular Dynamics Simulation of a Room Temperature Ionic Liquid. J. Phys. Chem. B 2005, 109, 5895–5902. (6) Buhl, M.; Chaumont, A.; Schurhammer, R.; Wipff, G. Ab Initio Molecular Dynamics of Liquid 1,3-Dimethylimidazolium Chloride. J. Phys. Chem. B 2005, 109, 18591–18599. (7) Bhargava, B. L.; Balasubramanian, S. Intermolecular Structure and Dynamics in an Ionic Liquid:

A Car-Parrinello Molecular Dynamics Simulation Study of 1,3-

Dimethylimidazolium Chloride. Chem. Phys. Lett. 2006, 417, 486–491. (8) Bhargava, B. L.; Balasubramanian, S.; Klein, M. L. Modelling Room Temperature Ionic Liquids. Chem. Commun. 2008, 3339–3351. (9) Canongia Lopes, J. N.; Deschamps, J.; Padua, A. A. H. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. J. Phys. Chem. B 2004, 108, 2038–2047. (10) Canongia Lopes, J. N.; Padua, A. A. H. Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. J. Phys. Chem. B 2004, 108, 16893–16898. (11) Canongia Lopes, J. N.; 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. 28 ACS Paragon Plus Environment

Page 28 of 41

Page 29 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(12) Canongia Lopes, J. N.; Padua, A. A. H.; Shimizu, K. Molecular Force Field for Ionic Liquids IV: Trialkylimidazolium and Alkoxycarbonyl-Imidazolium Cations; Alkylsulfonate and Alkylsulfate Anions. J. Phys. Chem. B 2008, 112, 5039–5046. (13) Sambasivarao, S. V.; Acevedo, O. Development of OPLS-AA Force Field Parameters for 68 Unique Ionic Liquids. J. Chem. Theory Comput. 2009, 5, 1038–1050. (14) Liu, Z.; Huang, S.; Wang, W. A Refined Force Field for Molecular Simulation of Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2004, 108, 12978–12989. (15) Liu, Z.; Wu, X.; Wang, W. A Novel United-Atom Force Field for Imidazolium-Based Ionic Liquids. Phys. Chem. Chem. Phys. 2006, 8, 1096–1104. (16) Liu, Z.; Chen, T.; Bell, A.; Smit, B. Improved United-Atom Force Field for 1-Alkyl-3methylimidazolium Chloride. J. Phys. Chem. B 2010, 114, 4572–4582. (17) Zhong, X.; Liu, Z.; Cao, D. Improved Classical United-Atom Force Field for ImidazoliumBased Ionic Liquids: Tetrafluoroborate, Hexafluorophosphate, Methylsulfate, Trifluoromethylsulfonate, Acetate, Trifluoroacetate, and Bis(trifluoromethylsulfonyl)amide. J. Phys. Chem. B 2011, 115, 10027–10040. (18) de Andrade, J.; Boes, E. S.; Stassen, H. Computational Study of Room Temperature Molten Salts Composed by 1-Alkyl-3-methylimidazolium Cations–Force-Field Proposal and Validation. J. Phys. Chem. B 2002, 106, 13344–13351. (19) Liu, X.; Zhang, S.; Zhou, G.; Wu, G.; Yuan, X.; Yao, X. New Force Field for Molecular Simulation of Guanidinium-Based Ionic Liquids. J. Phys. Chem. B 2006, 110, 12062–12071. (20) Zhou, G.; Liu, X.; Zhang, S.; Yu, G.; He, H. A Force Field for Molecular Simulation of Tetrabutylphosphonium Amino Acid Ionic Liquids. J. Phys. Chem. B 2007, 111, 7078– 7084.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(21) Chaban, V. Polarizability Versus Mobility: Atomistic Force Field for Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055–16062. (22) Chaban, V. V.; Prezhdo, O. V. A New Force Field Model of 1-Butyl-3-methylimidazolium Tetrafluoroborate Ionic Liquid and Acetonitrile Mixtures. Phys. Chem. Chem. Phys. 2011, 13, 19345–19354. (23) 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. (24) Shah, J. K.; Brennecke, J. F.; Maginn, E. J. Thermodynamic Properties of the Ionic Liquid 1n-Butyl-3-methylimidazolium Hexafluorophosphate from Monte Carlo Simulations. Green Chem. 2002, 4, 112–118. (25) Morrow, T. I.; Maginn, E. J. Molecular Dynamics Study of the Ionic Liquid 1-n-Butyl-3methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2002, 106, 12807–12813. (26) Cadena, C.; Maginn, E. J. Molecular Simulation Study of Some Thermophysical and Transport Properties of Triazolium-Based Ionic Liquids. J. Phys. Chem. B 2006, 110, 18026– 18039. (27) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. Molecular Modeling and Experimental Studies of the Thermodynamic and Transport Properties of Pyridinium-Based Ionic Liquids. J. Phys. Chem. B 2006, 110, 2821–2832. (28) Maginn, E. J. Atomistic Simulation of the Thermodynamic and Transport Properties of Ionic Liquids. Acc. Chem. Res. 2007, 40, 1200–1207. (29) Zhang, X.; Huo, F.; Liu, Z.; Wang, W.; Shi, W.; Maginn, E. J. Absorption of CO2 in the Ionic Liquid 1-n-Hexyl-3-methylimidazolium Tris(pentafluoroethyl)trifluorophosphate


Page 30 of 41

Page 31 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


([hmim][FEP]): A Molecular View by Computer Simulations. J. Phys. Chem. B 2009, 113, 7591–7598. (30) Borodin, O.; Smith, G. D. Structure and Dynamics of N-Methyl-N-propylpyrrolidinium Bis(trifluoromethanesulfonyl)imide Ionic Liquid from Molecular Dynamics Simulations. J. Phys. Chem. B 2006, 110, 11481–11490. (31) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113, 11463–11478. (32) Vatamanu, J.; Borodin, O.; Bedrov, D.; Smith, G. D. Molecular Dynamics Simulation Study of the Interfacial Structure and Differential Capacitance of Alkylimidazolium Bis(trifluoromethanesulfonyl)imide [Cnmim][TFSI] Ionic Liquids at Graphite Electrodes. J. Phys. Chem. C 2012, 116, 7940–7951. (33) Hu, Z.; Vatamanu, J.; Borodin, O.; Bedrov, D. A Molecular Dynamics Simulation Study of the Electric Double Layer and Capacitance of [BMIM][PF6] and [BMIM][BF4] Room Temperature Ionic Liquids Near Charged Surfaces. Phys. Chem. Chem. Phys. 2013, 15, 14234–14247. (34) Bedrov, D.; Borodin, O.; Li, Z.; Smith, G. D. Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 4984–4997. (35) Wang, Y.; Voth, G. A. Unique Spatial Heterogeneity in Ionic Liquids. J. Am. Chem. Soc. 2005, 127, 12192–12193. (36) Wang, Y.; Izvekov, S.; Yan, T.; Voth, G. A. Multiscale Coarse-Graining of Ionic Liquids. J. Phys. Chem. B 2006, 110, 3564–3575. (37) Wang, Y.; Jiang, W.; Yan, T.; Voth, G. A. Understanding Ionic Liquids through Atomistic



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and Coarse-Grained Molecular Dynamics Simulations. Acc. Chem. Res. 2007, 40, 1193– 1199. (38) Wang, Y.; Feng, S.; Voth, G. A. Transferable Coarse-Grained Models for Ionic Liquids. J. Chem. Theory Comput. 2009, 5, 1091–1098. (39) Yan, T.; Burnham, C. J.; Del Popolo, M. G.; Voth, G. A. Molecular Dynamics Simulation of Ionic Liquids: The Effect of Electronic Polarizability. J. Phys. Chem. B 2004, 108, 11877– 11881. (40) Yan, T.; Li, S.; Jiang, W.; Gao, X.; Xiang, B.; Voth, G. A. Structure of the Liquid-Vacuum Interface of Room-Temperature Ionic Liquids: A Molecular Dynamics Study. J. Phys. Chem. B 2006, 110, 1800–1806. (41) Yan, T.; Wang, Y.; Knox, C. On the Structure of Ionic Liquids: Comparisons between Electronically Polarizable and Nonpolarizable Models I. J. Phys. Chem. B 2010, 114, 6905– 6921. (42) Yan, T.; Wang, Y.; Knox, C. On the Dynamics of Ionic Liquids: Comparisons between Electronically Polarizable and Nonpolarizable Models II. J. Phys. Chem. B 2010, 114, 6886– 6904. (43) Urahata, S. M.;

Ribeiro, M. C. C. Structure of Ionic Liquids of 1-Alkyl-3-

methylimidazolium Cations: A Systematic Computer Simulation Study. J. Chem. Phys. 2004, 120, 1855–1863. (44) Micaelo, N. M.; Baptista, A. M.; Soares, C. M. Parametrization of 1-Butyl-3methylimidazolium Hexafluorophosphate/Nitrate Ionic Liquid for the GROMOS Force Field. J. Phys. Chem. B 2006, 110, 14444–14451. (45) Koddermann, T.; Paschek, D.; Ludwig, R. Molecular Dynamic Simulations of Ionic Liquids:


Page 32 of 41

Page 33 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A Reliable Description of Structure, Thermodynamics and Dynamics. ChemPhysChem 2007, 8, 2464–2470. (46) Koddermann, T.; Fumino, K.; Ludwig, R.; Canongia Lopes, J. N.; Padua, A. A. H. What Far-Infrared Spectra Can Contribute to the Development of Force Fields for Ionic Liquids Used in Molecular Dynamics Simulations. ChemPhysChem 2009, 10, 1181–1186. (47) Zhao, W.; Eslami, H.; Cavalcanti Welchy, L.; Muller-Plathe, F. A Refined All-Atom Model for the Ionic Liquid 1-n-Butyl 3-Methylimidazolium bis(Trifluoromethylsulfonyl)imide [bmim][Tf2N]. Z. Phys. Chem. 2007, 221, 1647–1662. (48) Wang, Y.; Pan, H.; Li, H.; Wang, C. Force Field of the TMGL Ionic Liquid and the Solubility of SO2 and CO2 in the TMGL from Molecular Dynamics Simulation. J. Phys. Chem. B 2007, 111, 10461–10467. (49) Tsuzuki, S.; Shinoda, W.; Saito, H.; Mikami, M.; Tokuda, H.; Watanabe, M. Molecular Dynamics Simulations of Ionic Liquids: Cation and Anion Dependence of Self-Diffusion Coefficients of Ions. J. Phys. Chem. B 2009, 113, 10641–10649. (50) Roy, D.; Maroncelli, M. An Improved Four-Site Ionic Liquid Model. J. Phys. Chem. B 2010, 114, 12629–12631. (51) Mondal, A.; Balasubramanian, S. Quantitative Prediction of Physical Properties of Imidazolium Based Room Temperature Ionic Liquids through Determination of Condensed Phase Site Charges: A Refined Force Field. J. Phys. Chem. B 2014, 118, 3409–3422. (52) Youngs, T. G. A.; Del Popolo, M. G.; Kohanoff, J. Development of Complex Classical Force Fields through Force Matching to ab Initio Data: Application to a Room-Temperature Ionic Liquid. J. Phys. Chem. B 2006, 110, 5697–5707. (53) Youngs, T. G. A.; Hardacre, C. Application of Static Charge Transfer within an Ionic-Liquid Force Field and Its Effect on Structure and Dynamics. ChemPhysChem 2008, 9, 1548–1558. 33 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(54) Bhargava, B. L.; Balasubramanian, S. Refined Potential Model for Atomistic Simulations of Ionic Liquid Bmim PF6. J. Chem. Phys. 2007, 127. (55) Dommert, F.; Wendler, K.; Berger, R.; Delle Site, L.; Holm, C. Force Fields for Studying the Structure and Dynamics of Ionic Liquids: A Critical Review of Recent Developments. ChemPhysChem 2012, 13, 1625–1637. (56) Salanne, M. Simulations of Room Temperature Ionic Liquids: From Polarizable to CoarseGrained Force Fields. Phys. Chem. Chem. Phys. 2015, 17, 14270–14279. (57) Choi, E.; McDaniel, J. G.; Schmidt, J. R.; Yethiraj, A. First-Principles, Physically Motivated Force Field for the Ionic Liquid [BMIM][BF4]. J. Phys. Chem. Lett. 2014, 5, 2670–2674. (58) Son, C. Y.; McDaniel, J. G.; Schmidt, J. R.; Cui, Q.; Yethiraj, A. First Principles United Atom Force Field for the Ionic Liquid [BMIM][BF4] : An Alternative to Charge Scaling. J. Phys. Chem. B 2016, 120, 3560–3568. (59) Hunt, P. A.; Kirchner, B.; Welton, T. Characterising the Electronic Structure of Ionic Liquids: An Examination of the 1-Butyl-3-methylimidazolium Chloride Ion Pair. Chem. Eur. J. 2006, 12, 6762–6775. (60) McDaniel, J. G.; Schmidt, J. R. Physically-Motivated Force Fields from Symmetry-Adapted Perturbation Theory. J. Phys. Chem. A 2013, 117, 2053–2066. (61) Rick, S. W.; Stuart, S. J. Potentials and Algorithms for Incorporating Polarizability in Computer Simulations. Rev. Comp. Chem. 2002, 18, 89–146. (62) Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6, 774–786. (63) Schmidt, J. R.; Yu, K.; McDaniel, J. G. Transferable Next-Generation Force Fields from Simple Liquids to Complex Materials. Acc. Chem. Res. 2015, 48, 548–556. 34 ACS Paragon Plus Environment

Page 34 of 41

Page 35 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(64) Stone, A. J. Distributed Multipole Analysis: Stability for Large Basis Sets. J. Chem. Theory Comput. 2005, 1, 1128–1132. (65) Ferenczy, G. G. Charges Derived from Distributed Multipole Series. J. Comput. Chem. 1991, 12, 913–917. (66) Williams, G. J.; Stone, A. J. Distributed Dispersion: A New Approach. J. Chem. Phys. 2003, 119, 4620–4628. (67) Misquitta, A. J.; Stone, A. J. Dispersion Energies for Small Organic Molecules: First Row Atoms. Mol. Phys. 2008, 106, 1631–1643. (68) Misquitta, A. J.; Stone, A. J.; Price, S. L. Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models Against SAPT(DFT) Energies. J. Chem. Theory Comput. 2008, 4, 19–32. (69) Hesselmann, A.; Jansen, G. Intermolecular Dispersion Energies from Time-Dependent Density Functional Theory. Chem. Phys. Lett. 2003, 367, 778–784. (70) Hesselmann, A.; Jansen, G.; Schutz, M. Density-Functional Theory-Symmetry-Adapted Intermolecular Perturbation Theory with Density Fitting: A New Efficient Method to Study Intermolecular Interaction Energies. J. Chem. Phys. 2005, 122, 14103–14119. (71) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Dispersion Energy from Density-Functional Theory Description of Monomers. Phys. Rev. Lett. 2003, 91, 33201–33204. (72) Werner, H.-J.; Knowles, P. J.; Lindh, R.; Manby, F. R.; Schutz, M.; Celani, P.; Korona, T.; Mitrushenkov, A.; Rauhut, G.; Adler, T. B. et al. MOLPRO, a Package of Ab Initio Programs, Version 2009.1; see http://www.molpro.net. (73) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D. et al. GROMACS 4.5: A High-Throughput



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845–854. (74) Taguchi, R.; Machida, H.; Sato, Y.; Smith, R. L. High-Pressure Densities of 1-Alkyl-3methylimidazolium Hexafluorophosphates and 1-Alkyl-3-methylimidazolium Tetrafluoroborates at Temperatures from (313 to 473) K and at Pressures up to 200 MPa. J. Chem. Eng. Data 2009, 54, 22–27. (75) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (76) Nose, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (77) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697, PRA. (78) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. (79) McDaniel, J. G.; Schmidt, J. R. First-Principles Many-Body Force Fields from the Gas Phase to Liquid: A “Universal” Approach. J. Phys. Chem. B 2014, 118, 8042–8053. (80) Shiflett, M. B.; Yokozeki, A. Liquid-Liquid Equilibria in Binary Mixtures of 1,3Propanediol + Ionic Liquids [bmim][PF6], [bmim][BF4], and [emim][BF4]. J. Chem. Eng. Data 2007, 52, 1302–1306. (81) Zhang, S.; Li, X.; Chen, H.; Wang, J.; Zhang, J.; Zhang, M. Determination of Physical Properties for the Binary System of 1-Ethyl-3-methylimidazolium Tetrafluoroborate + H2O. J. Chem. Eng. Data 2004, 49, 760–764.


Page 36 of 41

Page 37 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(82) Fuller, J.; Carlin, R. T.; Osteryoung, R. A. The Room Temperature Ionic Liquid 1-Ethyl-3methylimidazolium Tetrafluoroborate: Electrochemical Couples and Physical Properties. J. Electrochem. Soc. 1997, 144, 3881–3886. (83) Nishida, T.; Tashiro, Y.; Yamamoto, M. Physical and Electrochemical Properties of 1-Alkyl3-methylimidazolium Tetrafluoroborate for Electrolyte. J. Fluor. Chem. 2003, 120, 135– 141. (84) Shamsipur, M.; Beigi, A. A. M.; Teymouri, M.; Pourmortazavi, S. M.; Irandoust, M. Physical and Electrochemical Properties of Ionic Liquids 1-Ethyl-3-methylimidazolium Tetrafluoroborate, 1-Butyl-3-methylimidazolium Trifluoromethanesulfonate and 1-Butyl-1methylpyrrolidinium Bis(trifluoromethylsulfonyl)imide. J. Mol. Liq. 2010, 157, 43–50. (85) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 1. Variation of Anionic Species. J. Phys. Chem. B 2004, 108, 16593–16600. (86) Shiflett, M. B.; Yokozeki, A. Solubility and Diffusivity of Hydrofluorocarbons in RoomTemperature Ionic Liquids. AIChE 2006, 52, 1205–1219. (87) Gardas, R. L.; Freire, M. G.; Carvalho, P. J.; Marrucho, I. M.; Fonseca, I. M. A.; Ferreira, A. G. M.; Coutinho, J. A. P. High-Pressure Densities and Derived Thermodynamic Properties of Imidazolium-Based Ionic Liquids. J. Chem. Eng. Data 2007, 52, 80–88. (88) Klomfar, J.;

Souckova, M.;

Patek, J. Experimental P-ρ-T Data for 1-Butyl-3-

methylimidazolium Tetrafluoroborate at Temperatures from (240 to 353) K and at Pressures up to 60 MPa. J. Chem. Eng. Data 2011, 56, 426–436. (89) Letcher, T. M.; Reddy, P. Ternary (Liquid + Liquid) Equilibria for Mixtures of 1-Hexyl-3methylimidazolium (Tetrafluoroborate or Hexafluorophosphate) + Benzene + an Alkane at T=298.2 K and P=0.1 MPa. J. Chem. Thermodyn. 2005, 37, 415–421.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(90) Chaban, V. V.; Prezhdo, O. V. Ionic Vapor: What Does It Consist Of? J. Phys. Chem. Lett. 2012, 3, 1657–1662. (91) Neto, B. A. D.; Meurer, E. C.; Galaverna, R.; Bythell, B. J.; Dupont, J.; Cooks, R. G.; Eberlin, M. N. Vapors from Ionic Liquids: Reconciling Simulations with Mass Spectrometric Data. J. Phys. Chem. Lett. 2012, 3, 3435–3441. (92) Zaitsau, D. H.; Fumino, K.; Emel’yanenko, V. N.; Yermalayeu, A. V.; Ludwig, R.; Verevkin, S. P. Structure-Property Relationships in Ionic Liquids: A Study of the Anion Dependence in Vaporization Enthalpies of Imidazolium-Based Ionic Liquids. ChemPhysChem 2012, 13, 1868–1876. (93) Deyko, A.; Hessey, S. G.; Licence, P.; Chernikova, E. A.; Krasovskiy, V. G.; Kustov, L. M.; Jones, R. G. The Enthalpies of Vaporisation of Ionic Liquids: New Measurements and Predictions. Phys. Chem. Chem. Phys. 2012, 14, 3181–3193. (94) 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. (95) 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. (96) Li, S.; Banuelos, J. L.; Guo, J.; Anovitz, L.; Rother, G.; Shaw, R. W.; Hillesheim, P. C.; Dai, S.; Baker, G. A.; Cummings, P. T. Alkyl Chain Length and Temperature Effects on Structural Properties of Pyrrolidinium-Based Ionic Liquids: A Combined Atomistic Simulation and Small-Angle X-ray Scattering Study. J. Phys. Chem. Lett. 2011, 3, 125–130. (97) Russina, O.; Triolo, A.; Gontrani, L.; Caminiti, R. Mesoscopic Structural Heterogeneities in Room-Temperature Ionic Liquids. J. Phys. Chem. Lett. 2011, 3, 27–33.


Page 38 of 41

Page 39 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(98) Hettige, J. J.; Kashyap, H. K.; Margulis, C. J. Communication: Anomalous Temperature Dependence of the Intermediate Range Order in Phosphonium Ionic Liquids. J. Chem. Phys. 2014, 140, 111102. (99) Xu, W.-G.; Li, L.; Ma, X.-X.; Wei, J.; Duan, W.-B.; Guan, W.; Yang, J.-Z. Density, Surface Tension, and Refractive Index of Ionic Liquids Homologue of 1-Alkyl-3methylimidazolium Tetrafluoroborate [Cnmim][BF4] (n = 2,3,4,5,6). J. Chem. Eng. Data 2012, 57, 2177–2184. (100) 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. (101) Armstrong, J. P.; Hurst, C.; Jones, R. G.; Licence, P.; Lovelock, K. R. J.; Satterley, C. J.; Villar-Garcia, I. J. Vapourisation of Ionic Liquids. Phys. Chem. Chem. Phys. 2007, 9, 982– 990. (102) Kelkar, M. S.; Maginn, E. J. Calculating the Enthalpy of Vaporization for Ionic Liquid Clusters. J. Phys. Chem. B 2007, 111, 9424–9427. (103) Hayamizu, K.; Tsuzuki, S.; Seki, S.; Umebayashi, Y. Multinuclear NMR Studies on Translational and Rotational Motion for Two Ionic Liquids Composed of BF4 Anion. J. Phys. Chem. B 2012, 116, 11284–11291. (104) Noda, A.; Hayamizu, K.; Watanabe, M. Pulsed-Gradient Spin-Echo 1H and 19F NMR Ionic Diffusion Coefficient, Viscosity, and Ionic Conductivity of Non-Chloroaluminate RoomTemperature Ionic Liquids. J. Phys. Chem. B 2001, 105, 4603–4610. (105) Kruk, D.; Meier, R.; Rachocki, A.; Korpala, A.; Singh, R. K.; Rossler, E. A. Determining Diffusion Coefficients of Ionic Liquids by Means of Field Cycling Nuclear Magnetic Resonance Relaxometry. J. Chem. Phys. 2014, 140, 244509. 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(106) Tsuzuki, S. Factors Controlling the Diffusion of Ions in Ionic Liquids. ChemPhysChem 2012, 13, 1664–1670. (107) Zhang, S.; Sun, N.; He, X.; Lu, X.; Zhang, X. Physical Properties of Ionic Liquids: Database and Evaluation. J. Phys. Chem. Ref. Data 2006, 35, 1475–1517. (108) Stoppa, A.; Zech, O.; Kunz, W.; Buchner, R. The Conductivity of Imidazolium-Based Ionic Liquids from (-35 to 195) ◦ C. A. Variation of Cation’s Alkyl Chain. J. Chem. Eng. Data 2010, 55, 1768–1773. (109) Nakamura, K.; Shikata, T. Systematic Dielectric and NMR Study of the Ionic Liquid 1Alkyl-3-Methyl Imidazolium. ChemPhysChem 2010, 11, 285–294. (110) Ning, H.; Hou, M.; Mei, Q.; Liu, Y.; Yang, D.; Han, B. X. The Physicochemical Properties of Some Imidazolium-Based Ionic Liquids and Their Binary Mixtures. Sci. China Chem. 2012, 55, 1509–1518.


Page 40 of 41

Page 41 of 41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry


Ab Initio Force Fields for Imidazolium-Based Ionic Liquids - The

Recommend Documents