Letter pubs.acs.org/JPCL
First-Principles, Physically Motivated Force Field for the Ionic Liquid [BMIM][BF4] Eunsong Choi,†,§ Jesse G. McDaniel,‡,§ J. R. Schmidt,‡ and Arun Yethiraj*,‡ †
Department of Physics, University of Wisconsin, 1150 University Avenue, Madison, Wisconsin 53706, United States Department of Chemistry, University of Wisconsin, 1101 University Avenue, Madison, Wisconsin 53706, United States
‡
S Supporting Information *
ABSTRACT: Molecular simulations play an important role in establishing structure−property relations in complex fluids such as room-temperature ionic liquids. Classical force fields are the starting point when large systems or long times are of interest. These force fields must be not only accurate but also transferable. In this work, we report a physically motivated force field for the ionic liquid 1-butyl-3-methylimidazolium tetrafluoroborate ([BMIM][BF4]) based on symmetryadapted perturbation theory. The predictions (from molecular dynamics simulations) of the liquid density, enthalpy of vaporization, diffusion coefficients, viscosity, and conductivity are in excellent agreement with experiment, with no adjustable parameters. The explicit energy decomposition inherent in the force field enables a quantitative analysis of the important physical interactions in these systems. We find that polarization is crucial and there is little evidence of charge transfer. We also argue that the often used procedure of scaling down charges in molecular simulations of ionic liquids is unphysical for [BMIM][BF4]. Because all intermolecular interactions in the force field are parametrized from first-principles, we anticipate good transferability to other ionic liquid systems and physical conditions. SECTION: Liquids; Chemical and Dynamical Processes in Solution
I
theory (SAPT).11 All inter-molecular interactions are parametrized entirely from ab initio calculations, with no adjustable parameters. Intramolecular bond, angle, and dihedral potentials are taken from the literature12,13 for convenience; in principle, these terms could also be parametrized from ab initio calculations. We find that our new force field is highly accurate, with the predicted liquid properties (liquid density, enthalpy of vaporization, diffusion coefficients, viscosity, and conductivity) exhibiting excellent agreement with experiment. The computational expense of our force field is comparable to other polarizable models (∼10 times slower than nonpolarizable force field as implemented in GROMACS;14 other implementations, i.e., extended Lagrangian15 may be faster). The force field contains a one-to-one correspondence with the inherent energy decomposition of SAPT, allowing us to obtain new insight into the physical interactions important in ionic liquids and to shed light on existing empirical molecular dynamics force fields for these systems.9 The predictions of the force field for various physical properties are in good agreement with experiment. Table 1 compares simulation results for the density, enthalpy of vaporization, self-diffusion coefficient of anions and cations, viscosity, and ionic conductivity to experimental data at 300 K. The enthalpy of vaporization and dynamic properties of ionic liquids are often very difficult to reproduce using conventional
onic liquids are an important class of compounds that exhibit extremely low volatility, nonflammability, high thermal and chemical stability, high conductivity, and a wide liquid temperature range, resulting in extensive potential applications, including: solvents for synthesis and catalysis,1,2 media for gasseparation processes,3−5 electrolytes in fuel cells and batteries,6 and storage media for solar thermal electric power systems.7 The large number of possible cation−anion combinations yields the potential for tuning the properties for specific applications.8 To fully exploit this tunability, we must develop rational design principles that relate the properties of the individual ions to the macroscopic static and dynamic properties of the ionic liquid. Computer simulations provide a means of obtaining molecular-level insight that could be used to develop such design principles. There have been a large number of computational studies that have explored the microscopic mechanisms dictating the structure and dynamics of ionic liquids.9 The majority of these studies rely on force fields to describe the intra- and intermolecular interactions of the ions. There have been a few ab initio molecular dynamics (AIMD) simulations,10 but these are typically limited to tens of picoseconds due to their computational expense. Because the relaxation times in common ionic liquids can extend well into the nanosecond regime (or longer), obtaining accurate and transferable force fields is therefore of considerable importance. In this work, we obtain a first-principles-based force field for the ionic liquid 1-butyl-3-methylimidazolium tetrafluoroborate ([BMIM][BF4]) based on symmetry-adapted perturbation © 2014 American Chemical Society
Received: May 30, 2014 Accepted: July 22, 2014 Published: July 22, 2014 2670
dx.doi.org/10.1021/jz5010945 | J. Phys. Chem. Lett. 2014, 5, 2670−2674
The Journal of Physical Chemistry Letters
Letter
Table 1. Comparison of Predictions for Physical Properties of [BMIM][BF4] to Experiment at 300 K ρ (kg/m3) MD exp.
1181 120216
ΔH (kJ mol−1) 128 12817c
D+ (10−11 m2 s−1) b
a
1.8 (2.2 ) 1.6016
D− (10−11 m2 s−1)
η (cP)
σ (S/m)
1.4b (1.5a) 1.4916
93b (98a) 90.416
0.4b (0.5a) 0.3916
a
Simulated data reported in parentheses are from explicit MD simulations at 300 K and are subject to significant statistical uncertainty as described in the paper. bValues extrapolated from the VFT fit line. (See Figures 1 and 2.) cEstimated using the empirical correlation with the surface tension at T = 298 K.
Figure 1. Computed self-diffusion coefficients as a function of temperature. The VFT line is fit to MD results computed at the temperatures T = 350 K and higher. The inset figures show the same plots on a logarithmic scale.
Figure 2. Computed (a) shear viscosity and (b) ionic conductivity as a function of temperature. ηMD/σMD indicates viscosity/conductivity directly computed from the simulations using the Einstein relation, ηEXP/σEXP, which indicates experimental viscosity/conductivity.16 The VFT line is fit to MD results computed at the temperatures T = 350 K and higher in both cases. The inset figures show the same plots on a logarithmic scale.
nonpolarizable force fields.18−21 Borodin19 has developed a polarizable model for various ionic liquids including [BMIM][BF4], which yields accurate values for the enthalpy of vaporization and self-diffusion coefficients. However, that model was empirically parametrized to fit experimental density and dynamic properties. Chaban and coworkers18 introduced a reduced charge (nonpolarizable) force field for [EMIM][BF4] and [BMIM][BF4] developed based on Liu et al.’s model22 and reported good agreement with experiments. The diffusion coefficients, viscosity, and ionic conductivity are in good agreement with experiment over a range of temperatures. Because the dynamics of the ionic liquid is extremely sluggish, computing accurate transport properties at room temperature is a challenging task.21 We therefore fit the Vogel−Fulcher−Tamman (VFT)23−25 equation, f(T) = A exp((B)/(T − T0)), where A, B, and T0 are fit parameters, to the dynamic properties computed at temperature greater than or equal to 350 K, to estimate the room temperature result from the fit line. In Figures 1 and 2, extrapolated room-
temperature results are provided along with the directly computed results. The energy decomposition inherent to SAPT allows us to obtain physical insight into the important interactions in these systems. Such an analysis has been carried out for gas-phase dimers26−28 and clusters,29 and we investigate the decomposition for the liquid, in terms of energy contribution to the entire liquid and contributions to the energy of ion pairs in the liquid. From a 10 ns portion of the (equilibrated) MD simulation of [BMIM][BF4] at constant temperature of T = 300 K, we compute the average intermolecular exchange, electrostatic, induction, and dispersion interactions, employing the one-to-one correspondence between distinct force-field terms and SAPT interaction energies. (See the Supporting Information.) In addition to the liquid, we also calculate the energy decomposition of cation/anion pairs (termed ion complexes), which we define as those species bound by at least −300 kJ/mol, which is 85% of the strongest pairwise binding. Figure 3 depicts the different contributions to the 2671
dx.doi.org/10.1021/jz5010945 | J. Phys. Chem. Lett. 2014, 5, 2670−2674
The Journal of Physical Chemistry Letters
Letter
results31,33 and one should be skeptical toward conclusions based on a single methodology. Finally, we comment briefly on the commonly employed reduced/scaled charge models of [BMIM][BF4]. Such models are important because they can exhibit very good accuracy for the properties of neat ionic liquids at comparatively low computational cost.18,34,35 However, we suggest that explicitly polarizable models are essential to capture the correct physics in [BMIM][BF4], using the calculation of enthalpy of vaporization as an example. Because vaporizing an ionic liquid results in a gas phase ion pair, the enthalpy of vaporization is affected by differences in the polarization energy of the gas and liquid phase. The polarization energy is strictly attractive in both phases but is “optimized” for the close contact ion pair in the gas phase (i.e., the many-body contribution in the liquid is repulsive29), leading to a lower enthalpy of vaporization by ∼30 kJ/mol. (See Supporting Information.) In contrast, reduced/ scaled charge models may also predict lower (and thus more accurate18) enthalpies of vaporization, but in this case unphysically, by artificially reducing the overall cohesive energy of the liquid. Explicit polarization and scaled charges exhibit opposite effects on the cohesive energy of both the liquid and ion complexes; explicit polarization is strictly attractive, leading to lower (more negative) ion complex and liquid energies, while reduced/scaled charge models dramatically underestimate the electrostatic interactions. While it is common to employ scalings of 0.6 to 0.9 of the unit ion charge,9 our energy decomposition analysis shows that such scaling underestimates the cohesive energy of both the ion complex and the liquid by over 100 kJ/mol per ion pair! Scaled charge force fields may therefore be problematic when computing solute−solvent interaction energies and related properties. In summary, we have developed a first-principles-based, polarizable force field for [BMIM][BF4] from SAPT calculations. The MD simulations using our force-field show excellent agreement with available experimental data including density, enthalpy of vaporization, self-diffusion coefficients, shear viscosity, and conductivity over a wide range of temperatures. While popular scaled/reduced charge models have been successful at accurately predicting properties of neat ionic liquids, we propose that they cannot entirely capture the physical effects of explicit polarization and cannot be motivated by charge transfer for the case of [BMIM][BF4]. We expect that the advantages of utilizing a force field that captures all of the correct physical interactions may be more apparent when studying ionic liquid-based mixtures and interfaces or ionic liquid systems under abnormal conditions.
energy in the liquid and ion complexes, plotted as a fraction of the total.
Figure 3. Comparison of energy decomposition for entire liquid to strongly bound cation/anion pairs from a 10 ns trajectory of [BMIM][BF4] at T = 300 K and fixed volume.
There are several interesting features of the decomposition. (1) The relative contribution of electrostatics is approximately the same for the ion complexes and the liquid. This is not a priori obvious because the total electrostatic energy of the liquid is a complex sum of short-range and long-range pairwise ion−ion interactions. (2) The induction terms are a small but significant fraction for both the liquid (5%) and the ion complexes (9%). (3) The dispersion and exchange terms are small for the ion complexes (7−10%) but quite significant for the liquid (28−33%). However, these two terms are roughly equal in magnitude and of opposite sign and cancel for both the ion complexes and the liquid. To a good approximation the total energy mirrors the electrostatic energy. An examination of the induction energy sheds light on the importance of polarization and charge transfer in [BMIM][BF4]. Within the context of SAPT, the induction energy includes both polarization and charge-transfer contributions. However, we propose that polarization is the dominant component of the induction energy in [BMIM][BF4] and that charge transfer is negligible. We support this claim by calculating the ab initio (SAPT) induction energy for 1300 [BMIM][BF4] ion complexes from the liquid. This ab initio induction energy is extremely well-reproduced by the Drude oscillator polarization energy (itself parametrized to molecular response calculations), suggesting no substantial contribution from charge transfer. In addition, molecular charges for these ion complexes using three different charge distribution schemes (DMA, NAO, Bader) yield cation and anion charges very close to +1.0 and −1.0. (See the Supporting Information.) While we find negligible charge transfer for [BMIM][BF4], the existence of charge transfer in ionic liquids may depend on the specific ions. Ionic liquids with small anions such as Cl− seem to exhibit some charge transfer, similar to ionic melts.30 For example, the electronic structure of [BMIM][Cl] suggests that there is modest charge transfer (0.15e) from the anion to the cation.31 Using a very different charge partitioning scheme based on projections of plane-wave electron densities, Wendler et al. find significantly greater charge transfer of 0.4e for ion pairs 1,3-dimethylimidazolium chloride ([MMIM][Cl]), 1ethyl-3-methylimidazolium thiocyanate ([EMIM][SCN]), and 1-ethyl-3-methylimidazolium dicyanamide ([EMIM][DCA]).32 Because there is no unique way to partition molecule charges, different partitioning schemes may lead to very different
■
METHODS The method for constructing a force field based on SAPT calculations is identical to previous work, and we refer the reader to that paper for parametrization details.36 In this work (with the exception of charges, which are fit explicitly for BMIM and BF4), force-field parameters are derived for the imidazolium ring and the BF4 anion. Alkyl group parameters are taken from our previous force field,36 with the parameters of the α carbon atoms (adjacent to imidazolium ring) slightly modified to reflect the charge transfer to the ring; see the Supporting Information. In brief, atomic static and frequencydependent polarizabilites, used to obtain Drude oscillator charges and dispersion coefficients, are obtained using linear response calculations at the coupled-perturbed Kohn−Sham (CPKS) level of theory,37 as implemented in the CAMCASP 2672
dx.doi.org/10.1021/jz5010945 | J. Phys. Chem. Lett. 2014, 5, 2670−2674
The Journal of Physical Chemistry Letters
Letter
software package.38−42 Point charges are derived by fitting43 atomic charges to distributed multipole analysis (DMA)44,45 of the BMIM cation and BF4 anion separately. (The cation and anion exhibit full charges of +1.0 and −1.0, respectively.) Shortrange, charge penetration terms are fit to DFT-SAPT46 calculations of homomolecular dimers, namely, imidazolium dimers and BF4 dimers. All DFT-SAPT calculations are conducted using the MolPro 2009 software package.47 Deriving our force field requires several hundred DFT-SAPT calculations (each
