WWW.C-CHEM.ORG

FULL PAPER

Stability of Ion Triplets in Ionic Liquid/Lithium Salt Solutions: Insights from Implicit and Explicit Solvent Models and Molecular Dynamics Simulations Andrzej Eilmes* and Piotr Kubisiak Binding energies of ion triplets formed in ionic liquids by Li1 with two anions have been studied using quantum-chemical calculations with implicit and explicit solvent supplemented by molecular dynamics (MD) simulations. Explicit solvent approach confirms variation of solute-ionic liquid interactions at distances up to 2 nm, resulting from structure of solvation shells induced by electric field of the solute. Binding energies computed in explicit solvent and from the polarizable continuum model approach differ largely, even in sign, but relative

values generally agree between these two models. Stabilities of ion triplets obtained in quantum-chemical calculations for some systems disagree with MD results; the discrepancy is attributed to the difference between static optimized geometries used in quantum chemical modeling and dynamic strucC 2015 Wiley Periodicals, tures of triplets in MD simulations. V Inc.

Introduction

effect of the medium. Common methods of computationally inexpensive treatment of solvent are implicit solvent models (COSMO, conductor-like solvent model[12] or PCM, polarizable continuum model[13]) yielding corrections to energies due to dielectric continuum. Description of the ionic liquid as a dielectric medium in a model like PCM may be, however, questionable. One may argue that the static dielectric constant is not a well suited descriptor for a liquid composed entirely of ions.[14] Moreover, experimental data on dyes in ILs suggest that solvatochromic shifts of probe molecules induced by ionic liquids are larger than expected from the values of their static dielectric constant.[14,15] Though there are continuum solvent models designed specially for ionic liquids,[16] still the most widely available are the models rooted from the description of typical molecular solvents, for example, PCM. It is unclear whether PCM-like models are adequate for calculations of interaction energies in ionic liquids and, if the answer is positive, which value of dielectric permittivity should be used in such cases. In this article, we would like to investigate the interactions between ions in some typical ionic liquids and to compare predictions of continuous solvent approach with the results of explicit solvent modeling to obtain some clues on the applicability of both methods. A study of this kind was performed for ion pairs of EMIM-based ionic liquid using clusters up to

Ionic liquids (ILs) are organic salts with melting points below 100 C. Widespread interest in these systems is motivated by their properties, including high conductivity, high thermal stability, nonflammability and low vapor pressure at room temperatures. ILs are also good solvents for a wide range of chemicals. Owing to these properties ionic liquids found a wide range of applications not only as electrolytes in batteries, fuel cells and electrochemical devices but also as alternative solvents in organic synthesis and separation technology.[1–4] An important advantage of ILs is the possibility of tailoring liquid properties by changing anion/cation combination leading to a vast variety of solvents.[5] At microscopic level ILs are composed of cations and anions only, what makes them different from molecular liquids. Strong electrostatic interactions between ions result in unique properties of these systems.[6] Appropriate modeling of ion–ion interactions in IL is therefore a key point for understanding of their properties and rational design of new systems.[7–11] A large amount of relevant data can be obtained from molecular dynamics (MD) simulations; however, such calculations are time consuming and require proper setup, quite often development of force-field parameterization, and complicated calculation schemes to obtain for example, Gibbs free energies of solvation. Conversely, quantum-chemical (QC) calculations can easily and relatively fast provide structures of ions and their complexes and the interaction energies between ions or molecules in an ionic liquid and these data are often of great help in the understanding of underlying physics. Various methods of quantum-chemistry were applied to such computations. Although valuable information can be obtained from QC calculations for systems in vacuum, advanced modeling should take into account solvent effects changing the interaction energies between dissolved species because of the screening

DOI: 10.1002/jcc.23853

A. Eilmes, P. Kubisiak Jagiellonian University, Faculty of Chemistry, Department of Computational Methods in Chemistry, Ingardena 3, 30-060, Krak ow, Poland E-mail: [email protected] Contract grant sponsor: National Science Centre (Poland); Contract grant number: 2012/07/B/ST4/00573; Contract grant sponsor: PL-Grid computational infrastructure; Contract grant sponsor: European Regional Development Fund in the framework of the Polish Innovation Economy Operational Program (Part of the equipment used in calculations); Contract grant number: POIG.02.01.00-12-023/08 C 2015 Wiley Periodicals, Inc. V

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

1

FULL PAPER

WWW.C-CHEM.ORG

130 atoms.[17] For our work we chose some systems which have been already studied through quantum-chemical calculations. Both experimental and theoretical results suggest that in solutions of lithium salts in ionic liquids charge-carrying species are the ion triplets built of Li1 and two anions,[18–20] as a consequence of preferred fourfold coordination of lithium cation. When the anion of the IL (An1) differs from the anion of lithium salt (An2) homogenous [Li(An1)2]2 or mixed [Li(An1)(An2)]2 triplets may be formed depending on their relative stability. In a recent work[21] binding energies of such complexes in five ionic liquids were studied using quantumchemical calculations in vacuum or with one pair of ionic liquid ions. Here we will reexamine these systems using different solvent models and molecular dynamics simulations and we will discuss the correspondence between different approaches to see whether explicit and implicit solvent can yield similar predictions and whether the results of energy calculations agree with stabilities observed in MD.

Methodology In this work we studied five ionic liquids based on 1-ethyl-3methylimidazolium (EMIM) cation and five anions: BðCNÞ2 4 2 2 (bison), BF2 4 , PF6 , bis(fluorosulfonyl)imide (FSI ) and bis(trifluoromethane)sulfonimide (TFSI2). This is the same set of liquids and anions as used in a computational study on stability of ion triplets.[21] Gaussian 09 rev. A.02[22] was used to find optimal geometries of ion pairs and triplets and their energies. Geometry optimizations were performed in vacuum for LiAn pairs and [Li(An)2]2 triplets (where An2 is one of five anions considered) using MP2 method and aug-cc-pVDZ basis set. Additional sets of optimized structures of all 15 [LiAn1An2]2 triplets with all possible An1 and An2 combinations were obtained at the B3LYP/6-311G* and B3LYP/aug-cc-pVDZ levels. We intended to use quantum-chemical calculations for ion triplets solvated in a number of IL ion pairs, therefore we conducted some initial tests on performance of different methods. Although there were no problems for ion pairs or triplets isolated in vacuum, we noticed that for clusters containing an ion aggregate solvated in 3 to 15 EMIM1 – An2 pairs convergence of the SCF procedure depended on the density functional used. There were no convergence problems in Hartree–Fock calculations or in density functional theory calculations with local functionals. Conversely, when hybrid functionals were used, convergence was poor and for some systems calculations did not converge at all. The problem is apparently related to ionic nature of the system. This suggests a possible way to alleviate it by use of range-corrected fuctionals, and indeed we did not encounter problems in CAM-B3LYP or wB97X calculations. Therefore we chose CAM-B3LYP (with 6-311G* basis set) for quantumchemical studies of larger systems. In such cases TeraChem v. 1.5[23] running on machines equipped with Nvidia K20 or Nvidia M2090 cards was used to apply hardware acceleration on Graphics Processing Units. Additional Gaussian 09 energy calculations employing Polarizable Solvent Model were performed for ion pairs and triplets at vacuum geometries using 2

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

MP2 or CAM-B3LYP method and 6-311G* or aug-cc-pVDZ basis sets. Classical polarizable force field was used in molecular dynamics simulations. The majority of parameters (bonds, angles, torsions, charges, polarizabilities, and initial values of van der Waals parameters) for EMIM cation and all five anions were taken from the APPLE&P (Atomistic Polarizable Potential for Liquids, Electrolytes & Polymers) force field.[24] Few missing bonded parameters were adapted from Refs. [25] and [26]. Original APPLE&P force field uses Buckingham potential to describe van der Waals interactions. For our purposes we had to convert these parameters to Lennard-Jones-type 6–12 potential. This might be done by inspection of the Buckingham potential curve to find corresponding values of atomic radius and depth of the minimum. However, for given values of these parameters, the 6–12 potential is more repulsive at short distances than Buckingham potential and simple reuse of the Buckingham APPLE&P parameters in Lennard-Jones potential leads to an overestimate of repulsive interactions and too low density of the liquid. Therefore original APPLE&P atomic radii were reduced to reproduce in NPT simulations experimental densities of neat liquids with accuracy of about 0.01 g/cm3 (experimental and calculated densities are given in Table S1 in Supporting Information). We determined scaling factors by trial-and-error methods, aiming to introduce as small corrections as possible. In the case of EMIM cation atomic vdW radii were reduced by 15%; for anions the reduction factors depended on the anion and were between 7.5, and 15%; individual values are listed in Table S2 in Supporting Information. A separate fitting of Lennard-Jones potential was performed for pairwise interactions of Li1 with atoms of IL anions. Van der Waals interaction parameters were fitted for best reproduction of both the geometries and interaction energies for Li1-An2 pairs and [Li(An)2]2 triplets obtained from the quantum-chemical calculations at the MP2/aug-ccpVDZ level. Initial geometries for MD simulations were prepared with Packmol package.[27] Ion pairs or triplets (in B3LYP/aug-ccpVDZ geometry) were solvated in 100 EMIM-anion pairs of the IL in the cubic simulation box. In simulations of ion triplets or individual anions one extra EMIM cation was added to ensure the neutrality of the system. Size of the box was set according to the density of given liquid. MD calculations were conducted using Tinker v. 5.1 molecular modeling software.[28] Most of the simulations were performed in NVT ensemble at T 5 298 K using the Berendsen thermostat. The exception were the systems in EMIM-PF6 in which cases temperature was raised to 353 K because of higher melting point of this IL (335 K). Default Tinker value of 0.1 ps was set as the coupling time for temperature bath. Beeman algorithm with time step of 1 fs was used to integrate the equations of motion. Cutoff of 9 A˚ was applied to van der Waals interactions; electrostatic interactions were treated by Ewald summation. Induced dipoles were calculated iteratively assuming mutual interactions between polarizable units to achieve self-convergence; the Thole scheme of short-range polarization damping[29] was applied. In a set of initial simulations, the geometry of the LiAn ion pair WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

or the [LiAn1An2]2 triplet was kept frozen at the geometry optimized in quantum-chemical calculations. We also performed additional simulations for solvated individual ions (Li1 and 5 anions) in fixed geometry. After about 4–6 ns of relaxation, approximately 400 frames from the subsequent part of the trajectory (about 2 ns long) were recorded in 5 ps intervals and used in further analysis. Simultaneously, for simulations of ion pairs and triplets, five individual frames of each trajectory were selected as starting points for new simulations in which geometry of the central ion pair or triplet was allowed to change. Trajectories from these MD runs were recorded up to 10 ns.

Quantum-Chemical Calculations Relative stabilities of ion pairs and triplets investigated in this paper were already computed and analyzed in[21] using DFT methodology and B3LYP functional. Vacuum values are the starting point for further modeling; therefore, we needed to repeat such calculations to obtain reference data for our forcefield calculations and for QC methods used in this work. Here we briefly summarize the results. Stabilization energies of ion pairs LiAn and ion triplets [LiAn1An2]2 in vacuum were calculated in a standard way as the difference between the energy of the complex and the energies of its constituents: 0 0 0 DELiAn 5ELiAn 2ðELi0 1EAn Þ

(1a)

0 0 0 0 DELiAn1An2 5ELiAn1An2 2ðELi0 1EAn1 1EAn2 Þ

(1b)

where EX0 stands for the energy of species X calculated in vacuum. In Table S3 (Supporting Information), we compare binding energies of LiAn ion pairs obtained from MP2 calculations (in the procedure of FF fitting) with CAM-B3LYP calculations and with the values calculated using our force field. Data obtained by different methods agree within 5 kcal/mol with force fieldcalculated values usually failing between MP2 and DFT results. Binding energies decrease in the order BðCNÞ2 >> 4 2 2 2 PF2 > FSI > TFSI > BF (note that these values are negative, 6 4 therefore, strength of the interaction increases in the above sequence). This is the same sequence as reported in the literature for B3LYP calculations.[21] In all cases Li1 is complexed in a bidentate manner, Table S4 in Supporting Information collects distances to coordinating atoms obtained in DFT calculations and in the structures optimized with force field. Binding energies for [LiAn1An2]2 ion triplets calculated at force-field-optimized geometries in vacuum using force field are collected in Table 1. To facilitate comparison with quantum-chemical results, values obtained from CAM-B3LYP/6311G* and MP2/6-311G* calculations (at B3LYP geometries) are also presented. We included also B3LYP/6-3111G* results of former theoretical work.[21] Depending on the method, absolute energies may differ up to 10 kcal/mol, but general pattern of stabilization energies for triplets follows the binding energies of LiAn pairs, that is, binding energies in Table 1 decrease within rows or within columns in the sequence

Table 1. Binding energies (in kcal/mol) of [LiAn1An2]2 ion triplets calculated in vacuum using force field (FF) and quantum chemical methods. IL anion Salt anion

BðCNÞ2 4

FF BðCNÞ2 2150.9 4 2166.0 PF2 6 FSI2 2166.6 2169.4 TFSI2 BF2 2175.1 4 MP2/6-311G* BðCNÞ2 2149.2 4 PF2 2164.6 6 2166.0 FSI2 TFSI2 2168.8 2171.3 BF2 4 CAM-B3LYP/6-311G* BðCNÞ2 2151.3 4 PF2 2167.8 6 FSI2 2170.5 2173.4 TFSI2 BF2 2174.5 4 B3LYP/6-3111G*, Ref. [21] 2 BðCNÞ4 2145.8 2162.8 PF2 6 FSI2 2165.2 TFSI2 2167.5 2170.7 BF2 4

PF2 6

FSI2

TFSI2

BF2 4

2166.0 2181.9 2180.9 2183.4 2190.4

2166.6 2180.9 2180.8 2183.3 2189.3

2169.4 2183.4 2183.3 2185.5 2191.8

2175.1 2190.4 2189.3 2191.8 2198.7

2164.6 2182.5 2182.3 2184.5 2188.4

2166.0 2182.3 2183.0 2185.7 2188.1

2168.8 2184.5 2185.7 2188.8 2190.1

2171.3 2188.4 2188.1 2190.1 2194.2

2167.8 2186.6 2187.7 2190.1 2192.4

2170.5 2187.7 2189.4 2191.9 2193.4

2173.4 2190.1 2191.9 2194.5 2195.7

2174.5 2192.4 2193.4 2195.7 2198.1

2162.8 2180.0 2182.4 2184.5 2186.9

2165.2 2182.4 2184.0 2186.2 2189.3

2167.5 2184.5 2186.2 2188.3 2191.2

2170.7 2186.9 2189.3 2191.2 2194.8

Literature data from Ref. [21] are included for easy comparison. Colors denote stability of mixed triplets with respect to homogenous triplets; see text for details.

shown above (in accord with Fig. 2 from Ref. [21]). Small devia2 tion is observed in FF results for PF2 anions—in 6 and FSI 2 some cases triplets with PF6 are slightly more stable than corresponding systems with FSI2. Listed values correspond to vacuum results, therefore offdiagonal entries in Table 1 (mixed anion triplets) are symmetric with respect to the diagonal. However, when analyzing stabilities of ion triplets one should note that in real salt/IL systems concentrations of salt anion and the anion of the ionic liquid are different thus the lower and the upper triangle of the table are not equivalent. In particular, homogenous triplets can be formed only from the IL anion. Relative stability of ion triplets should be therefore analyzed within columns, that is, so that different salts in given IL are compared. Mixed triplets can be stable only when their binding energy is lower than that for homogenous triplets, otherwise they are supposed to break up, for example, [Li(BF4)(B(CN)4)]2 triplet is stable in EMIM-B(CN)4 but not in EMIM-BF4 in which case formation of [Li(BF4)2]2 is expected. We labeled mixed triplets in Table 1 according to the stability predicted from vacuum calculations. Different QC methods (also those not shown in Table 1) give slightly different binding energies but their relative values usually do not differ between methods more than 4 kcal/mol, therefore we arbitrarily used this value as the uncertainty threshold. All mixed triplets with binding energy at least 4 kcal/mol weaker than homogenous triplet are considered as unstable (marked red), conversely, mixed triplet is stable (green entries in Table 1) Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

3

FULL PAPER

WWW.C-CHEM.ORG

EMIM-BF4 should therefore form [Li(BF4)2]2 triplets; the sole 2 uncertain case with BF2 4 anion is the [Li(BF4)(TFSI)] triplet. Having obtained some clues about stability of LiAn pairs and [LiAn1An2]2 ion triplets from quantum-chemical calculations in vacuum we turned our attention to the main goal of our work – analysis how the solvent effect of an ionic liquid influences the stabilization energies of ion triplets and whether the predictions of explicit and implicit solvent models will be comparable. In a liquid, binding energy of an ion aggregate is measured with respect to the energies of individual ions solvated in bulk solvent. To obtain necessary quantities we used the trajectories recorded in MD simulations with frozen geometry of ions. Geometries of individual ions, ion pairs or triplets surrounded by increasing number of ion pairs of the IL were extracted from the frames of the trajectory using the Trajectory Sculptor tool.[30] Energies of these systems (solvated ions 1 IL pairs) were calculated using our force field and then averaged over whole trajectory. Owing to the speed of FF energy calculations we were able to easily use up to 60 solvating IL pairs and about 300 geometries for each system type and size. One possibility of estimating binding energy of an ion aggregate solvated in N ion pairs of the IL would be a simple modification of formulas (1a,b) used in vacuum: n2 N N DELiAn 5ELiAn 2ðELin1 1EAn Þ

(2a)

n2 n3 N N 5ELiAn1An2 2ðELin1 1EAn1 1EAn2 Þ DELiAn1An2

(2b)

or

where EXn is the energy of species X in n pairs of IL ions and X

ni 5N

i

Figure 1. Solute–ionic liquid EMIM-TFSI interaction energies for ions a), ion pairs b), and ion triplets c) calculated from the force-field for increasing N according to eq. (4).

when it is bound at least 4 kcal/mol stronger than the homogenous one. Triplets with binding energy falling within 64 kcal/ mol interval around the stabilization energy of homogenous triplet are classified as unclear cases (yellow entries). Color patterns arising for different methods of energy calculations are similar with uncertain classification obtained mainly for triplets 2 2 consisting of PF2 6 , FSI and TFSI anions. Because of relatively 1 weak interaction between Li and BðCNÞ2 4 , mixed triplets containing bison anion are always unstable except in EMIM1 B(CN)4. Conversely, BF2 4 binds strongly to Li and triplets with tetrafluoroborate are almost always stable. Most salts in 4

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

so that the number of IL pairs on both sides of the formula is the same. Binding energy in bulk liquid could be estimated from an extrapolation to an infinite N. Sample plot of binding energy values versus N calculated for ion triplets in EMIM-TFSI for ni 5 N/3 is shown in Figure S1 in Supporting Information. Inapplicability of such approach becomes obvious: because the number of ionic liquid ions solvating the aggregate is three times larger than that for individual ions, a dominant contribution to the calculated value is not the binding energy of the triplet but the cohesion energy of the liquid and therefore DE N constantly decreases with increasing N. We used therefore an improved scheme (based on a thermodynamic cycle) int ðnÞ

int ðnÞ

N 0 DELiAn 5DELiAn 1DELiAn 2ðDELi int ðnÞ

int ðnÞ

N 0 DELiAn1An2 5DELiAn1An2 1DELiAn1An2 2ðDELi

int ðnÞ

1DEAn

Þ

int ðnÞ

(3a) int ðnÞ

1DEAn1 1DEAn2 Þ (3b)

int ðnÞ

with DEX being the interaction energy of species X with n IL solvent pairs calculated as WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

Figure 2. Binding energies of ion pairs a) and triplets b) in EMIM-TFSI calculated from the force-field for increasing N according to eqs. (3a,b). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] int ðnÞ

DEX

5EXn 2ðEX0 1Esolv2X Þ

(4)

n or 0 specify the number of IL pairs solvating X (thus EX0 is the energy calculated in vacuum) and Esolv2X is the energy of n solvent pairs in the geometry of the system [X 1 solvent] (therefore for each system [X 1 solvent] an additional energy calculation is required after X is removed from the “droplet”). With this definition scatter of calculated binding energies is greatly reduced because energy fluctuations inevitably arising from different arrangements of solvent ions partially cancel in (4). In Figure 1, we displayed an example how the interaction energies with the liquid obtained for different ions solvated in EMIM-TFSI change with increasing amount of the solvent. When there are only a few solvating ion pairs of the IL (small N values) each next ion of the liquid is being added to the cluster close to the solvated ion and the interaction energy decreases rapidly with increasing N. When the first solvation layer of the ion, pair or triplet is completed, DE int starts to oscillate and its changes reflect the structure of the solvation shells differing in orientation of effective dipole moments of solvating IL ions with respect to the central charge (or dipole in case of solvated ion pair). Electric field of smaller ions is stronger; inducing more ordered structure of the solvent and

FULL PAPER

therefore more pronounced oscillations of interaction energy. Accordingly, the effect is the largest for Li1 cation, in which case also the interaction energy with the solvent has the most negative value. Long-range nature of the solvent-structuring effect is readily noticeable; oscillations of the interaction energy are still observable for 60 EMIM-TFSI ion pairs, corre˚ distance from the central ion(s) to the sponding to about 18 A outermost ion pair of the droplet. Positions of the maxima and minima of interaction energy depend on the size of the central solute as seen in Figure 1 but also on the size of the anion An1 of the EMIM-An1 IL. Such oscillatory behavior of electrostatic interactions have been shown in studies of screening effects in ionic liquids [31–33] and the variations of solvent potential have been found to persist at distances up to 2 nm.[34] Plots of interaction energy vs. N obtained for single ions and ion triplets are similar, with the major difference that the minima shift to larger N for triplets because of their larger size. Oscillations of DE int are also visible for solvated ion pairs, although trends exhibited for large N for different ions are less regular than in the case or ions or triplets. One may note that DE int calculated for ion pairs decreases very rapidly for small N values, large interaction energy is reached even for 2–3 IL pairs and, unlike two other types of system, the first maximum is deeper than the second. Conversely, for ions and triplets there is a decrease of interaction energy between the first two minima. This behavior may be rationalized by the charge of the central solute. In our calculations size of the [solute – IL] aggregate was increased in steps of one IL pair. For solvated Li – anion pair (electrically neutral) and small N values, anions and cations of the IL can be arranged in a way that optimizes their interactions with the dipole moment of the solute. For isolated ions or triplets, the solute is charged and for small number of IL pairs at least one of IL ions will be always in unfavorable position, therefore, the strength of stabilizing interaction with the solvent is reduced. This also explains why the trends for different ions and triplets are similar – oscillations of DE int reflect the ordering of solvation shells resulting from the charge of the solute, therefore they do not depend much of the IL and solute (apart from the effect of ion size). Conversely, such structuring of the liquid caused by a dipole of Li-An2 ion pair is weaker than for a central charge and other interactions may become comparable with electrostatics causing the trends for different ion pairs to look quite differently. As readily seen in Figure 1 even for 60 pairs of IL ions we did not reach the convergence of DE int . As the differences of these interaction energies are used to calculate the binding energy in the solvent according to the formula (3a,b), also the latter quantity did not saturate for largest N used in our calculations. Indeed, oscillatory character of binding energy is visible in Figure 2 which shows the DE N calculated for different ion pairs and triplets in EMIM-TFSI. Changes of binding energy DE N reflect the DE int dependence on N. For small numbers of IL ion pairs DE N increases rapidly from the value of binding energy calculated in vacuum (lower than 2150 kcal/mol) to positive values. After oscillation around the first maximum Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

5

FULL PAPER

WWW.C-CHEM.ORG

Figure 3. Comparison of solute–liquid interaction energies in EMIM-TFSI a) and binding energies of ion triplets b) calculated using force field (continuous lines) and CAM-B3LYP/6-311G* calculations (points and broken lines). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

(appearing at about 10 pairs of IL), changes of DE N become much smaller. Lines obtained for different triplets in EMIM-TFSI above N 5 30 are more-less parallel, therefore although absolute binding energies of triplets change with N, their relative energies remain approximately constant. This behavior is observed for most IL investigated here (Cf. Supporting Information). This is usually not the case for DE N of ion pairs because of rather irregular changes of their interaction energies as shown in Figure 1b. Relative energies of ion pairs in the solvent still depend much on N even for N 5 50 – 60 (Cf. also plots of DE N for pairs in EMIM-PF6 or EMIM-BF4 presented in Supporting Information). In force field based calculations we did not reach convergence of binding energies even for largest droplets of explicit solvent. System sizes attainable in quantum-chemical calculations are much smaller (even with hardware acceleration). Nevertheless, we performed some checks of CAM-B3LYP/6311G* energy calculations for selected ions and ionic liquids up to N 5 15. An example is presented in Figure 3. There is a difference between QC and FF-based DE int values – for example, DFT calculations yield stronger interaction of Li1 with EMIM-TSFI liquid than molecular mechanics, conversely QCcalculated interactions between IL and anions are weaker. Likewise, binding energies of triplets obtained from the two approaches differ. However, the dependences of DE int versus N are essentially the same apart from a rather constant shift. Clearly, QC values for ions follow the FF-based curve and the picture for binding energies DE N is similar, though it is more difficult to be noticed, because in Figure 3 we are in the range of most rapid changes of DE N, and also the averaging of QC data is significantly worse (20 values per point instead of 300 used in FF calculations). The conclusion from Figure 3 is that the binding energies computed from the force field reproduce the trends observed in quantum-chemical calculations. Before discussing binding energies in more details, we present the results obtained within continuous solvent model. We 6

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

applied the Polarizable Continuum Model either with default Gaussian 09 parameters used to construct the molecular cavity or with the atomic radii derived from the SMD[35] model. MP2 or CAM-B3LYP methodology and 6-311G* basis set were used. We tested a range of values of static dielectric constant, up to e 5 80 to investigate the general trend although the values of dielectric permittivity of our ILs are much smaller: 12.3 for EMIM-TFSI or 13.6 for EMIM-BF4[36] and for most methylimidazolium-based ILs e is between 12 and 15.[37] Binding energies in the PCM model are obtained in analogy to eqs. (1a,b) with the sole modification that the energies of individual species are calculated in the solvent rather than in vacuum. To check the effect of the basis set on stabilization energies we repeated energy calculations in vacuum and for two values of static dielectric constant (e 5 15 and 80) in augcc-pVDZ basis. Results are compared in Figure S2 in Supporting Information. Typically, increase of the basis set decreases stabilities of ion aggregates in MP2 calculations and increases them in CAM-B3LYP computations. Differences between energies obtained in these two basis sets do not exceed 3 kcal/ mol. Accordingly, relative stabilization energies for different anions are affected by less than 3 kcal/mol; in most cases the difference is between 1 and 2 kcal/mol. For some ionic liquids dependencies of binding energy of ion pairs on the dielectric permittivity of the solvent calculated in implicit solvent have been already presented in recent papers,[38,39] therefore, in Figure 4, we show the sample MP2 results obtained for ion triplets containing TSFI anion (plots for ion pairs and other triplets are available in the Supporting Information). The picture is similar to that reported in the case of ion pairs: when e increases from the vacuum value, initially DE changes very fast, then the slope of the dependence decreases and above e 5 50 further increase of dielectric constant has little effect on binding energy. The effect of different atomic radii (default vs. SMD) is mainly a vertical shift of all curves: SMD-based energies are more negative than the values WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 2. Solute–ionic liquid interaction energies (in kcal/mol) obtained from the force field calculations for ions, ion pairs and triplets in explicit EMIM-An ionic liquid and from the PCM quantum-chemical approach (MP2/6-311G*) with default or SMD-derived cavity. Solvent

BðCNÞ2 4

Ions 273.3 EMIM-B(CN)4 EMIM-FSI 268.5 EMIM-PF6 262.8 EMIM-TFSI 268.8 EMIM-BF4 261.7 MP2/PCM-def, e 5 15 240.8 MP2/PCM-SMD, e 5 15 229.9 MP2/PCM-def, e 5 80 243.2 MP2/PCM-SMD, e 5 80 232.5 ion pairs Li1 – An2 EMIM-B(CN)4 253.3 EMIM-FSI 260.3 EMIM-PF6 257.7 EMIM-TFSI 264.4 EMIM-BF4 263.2 MP2/PCM-def, e 5 15 233.6 MP2/PCM-SMD, e 5 15 211.5 MP2/PCM-def, e 5 80 236.8 MP2/PCM-SMD, e 5 80 213.6 ion triplets [Li(TFSI)An]2 EMIM-TFSI 270.3 MP2/PCM-def, e 5 15 237.4 MP2/PCM-SMD, e 5 15 218.2 MP2/PCM-def, e 5 80 240.4 MP2/PCM-SMD, e 5 80 221.1

Figure 4. Binding energies of [Li(TFSI)An]2 ion triplets in solvents with increasing dielectric constant calculated at the MP2/6-311G* level using PCM approach with default a) and SMD b) atomic radii. Shaded area marks the range of e typical for EMIM ionic liquids. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

obtained with standard molecular cavity. Binding energies of all ion triplets calculated for SMD cavity are negative even for highest values of the dielectric constant. For default cavity binding energy of some rather weakly bound triplets (these with BðCNÞ2 4 ion) becomes slightly positive in strongly polar solvents. Similar results were obtained in CAM-B3LYP calculations. At this point we can compare the ion–solvent interaction energies calculated within explicit and implicit solvent models. Table 2 collects the interaction energies obtained from force field computations with N 5 60 IL ion pairs for individual ions and ion pairs in all five IL liquids and for TFSI-containing triplets in EMIM-TFSI. The PCM-based MP2 interaction energies for two values of dielectric constant (e 5 15—a typical value for ILs and e 5 80 as a case of highly polar solvent) are given for default and SMD atomic radii. Interaction energies calculated from the force field for the same system in different IL liquids typically differ by 10–15 kcal/mol–these values may reflect a real change in DE int , however, a large part of the difference comes from the oscillatory behavior of the interaction energy still pronounced at N 5 60. Though the FF-based DE int for different anions vary with the IL, some general trends may be observed. For the five anions

FSI2

PF2 6

TFSI2

BF2 4

Li1

279.0 273.2 267.2 274.7 265.8 243.7 233.7 246.3 236.5

278.1 213.7 268.1 275.8 265.7 250.0 251.5 253.0 254.5

280.0 216.8 267.4 279.7 269.1 242.4 234.1 245.1 237.1

282.9 278.9 269.4 283.8 271.7 254.9 256.8 258.1 260.1

2128.2 2136.9 2155.8 2144.7 2173.0 2115.0 285.1 2121.6 290.1

248.6 251.3 252.5 254.0 261.5 227.3 26.3 230.0 27.7

247.3 249.5 251.6 253.3 262.1 232.9 219.7 235.8 221.3

252.7 254.0 252.9 252.5 268.5 225.7 27.1 228.4 28.7

244.4 244.7 252.3 248.2 256.8 232.1 219.8 234.8 221.4

274.0 235.6 217.8 238.1 220.4

271.4 238.6 230.0 241.4 232.7

276.5 235.3 219.2 237.9 222.0

271.7 239.8 231.8 242.7 234.6

considered here interaction energies change in the order 2 2 2 BF2  PF2 4 < TFSI < FSI 6 < BðCNÞ4 . Much stronger interaction with the solvent is observed for the Li1 ion with the DE int ranging from 2128 kcal/mol in EMIM-B(CN)4 to as much as 2173 kcal/mol in EMIM-BF4. In this case large differences between ionic liquids are not surprising as the electric field of small lithium cation leads to strong reorganization of the solvent in the solvation shell. Apparently, strength of the Li1 interaction with ionic liquid increases (i.e., DE int becomes more negative) with decreasing size of the IL anion. Values of DE int obtained from the PCM model decrease as expected with increasing solvent polarity. In general, strength of the interaction with implicit solvent is smaller (less negative values) than predicted in explicit solvent calculations, even for largest values of dielectric permittivity. The difference increases even more for SMD cavity – results calculated for SMD radii are less negative than obtained for default cavity, the difference is the largest for Li1 in which case it amounts to about 30 kcal/mol. In the PCM model with SMD radii also the differences between different solutes in the same solvent are larger. PCM solvent apparently yields different sequence of interaction energies, for example, for the anions DE int increases from 2 tetrafluoroborate to bison anion in the order BF2 4 < PF6 < 2 2 2 2 TFSI  FSI < BðCNÞ4 (PF6 anion changed its position in the sequence compared to the FF results). Differences between explicit and PCM results explain the difference in the sign of binding energies. The effect is mainly due to Li1 interaction energy, which is subtracted in (3) from the value of binding energy in vacuum to get the binding Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

7

FULL PAPER

WWW.C-CHEM.ORG

energy in the solvent. Strong stabilization of lithium cation in explicit solvent results in positive binding energies. In PCM results with SMD cavity this stabilization is much smaller, therefore even in polar liquids binding energies are negative except for systems with relatively weakly bound bison anion. In the middle between these two cases is DE int calculated with default radii, accordingly the default PCM binding energies may have positive or negative sign depending on the anion and polarity of the liquid. As we have shown above, values of binding energy obtained from the implicit solvent model and calculated using explicit solvent and force field differ even in sign. Part of this difference is related to the solvent–ion interaction energies giving a uniform shift of DE N. It may be therefore reasonable to compare not the absolute values but relative binding energies of ion triplets or pairs. From the data presented in Figure 2a it is difficult to decide which size N may be used to estimate binding energies for ion pairs because of different behavior of curves for different ions. Conversely, DE N versus N dependences for triplets are similar and relative energies do not change much above certain size of the droplet (as indicated by approximately parallel lines in Fig. 2b), therefore we decided to present relative binding energies only for triplets. In Figure 5, we collected relative values of DE N for [LiAn1An2]2 triplets in all five EMIM-An1 ionic liquids. In most cases the triplet [LiAn1(TFSI)]2 has the most negative binding energy, therefore it was set as zero of relative energy scale. In the plot we included vacuum data, values obtained from the PCM model with default and SMD cavity calculated at e 5 15 or e 5 80, and the results of force field energy calculations in 60 ion pairs of explicit solvent. In vacuum the triplet containing BF2 4 anion is always the most stable (1–5 kcal/mol lower than the triplet with TFSI2) and the systems with BðCNÞ2 4 are remarkably higher in energy (about 10 kcal/mol) than the other four triplets. These features look different in the case of results obtained from the solvent model (implicit or explicit). The BF2 4 triplet is never the most stable one and the gap between BðCNÞ2 4 triplets and other systems is reduced. The energy of the triplet with PF2 6 in the solvent, unlike vacuum calculation, is always noticeably higher than the energies of FSI2 and TFSI2 triplets. The main difference between relative binding energies calculated in the PCM model and in the explicit solvent is that in two ILs the latter predicts that the FSI2 triplet is more stable than the one with TFSI2. Otherwise the agreement between relative energies in continuous and explicit solvent is surprisingly good (in view of the fact that the absolute values may even differ in sign between the two models) and major differences with respect to vacuum data are well captured. Let us discuss the problem of the sign of calculated binding energies. All values calculated for N 5 60 in explicit ionic liquid are positive for all ion pairs and triplets. From energetic point of view this means that the LiAn2 salt dissociates in EMIM-An1 liquid, therefore binding energies just agree with known solubility of lithium salts in ionic liquids. For the salt sharing common anion with IL (i.e., LiAn in EMIM-An) positive binding energy means that the Li1-An2 pair will be present in solution 8

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

Figure 5. Relative binding energies for [Li(An1)(An2)]2 ion triplets calculated at the MP2/6-311G* level in vacuum, within PCM approach with default or SMD radii and two values of the dielectric constant and in explicit EMIM-An1 ionic liquid (data for N 5 60). Color codes for the An2 anion as in Figure 4. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

(because there is no other anion to coordinate the lithium cation) but the cation will dynamically exchange anions in its solvation shell. Using such reasoning for ion triplets we may expect that all triplets will exchange the An1 anion of the ionic liquid (because all are unstable judging from positive DE N ) and may also swap the salt An22 anion for the An12 anion of the IL in cases where it lowers the energy (i.e., leads to formation of a triplet with less positive binding energy). Implicit solvent models with dielectric constant of the solvent set in the range corresponding to liquids investigated here (e 5 12–15) give negative binding energies of ion pairs and triplets. SMD model yields negative values even for much more polar solvents. This is the common pitfall of using continuous solvent models to calculate interactions between WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

FULL PAPER

ions – in most cases they predict stability of ion aggregates suggesting that the salt does not dissociate in solution in disagreement with experiment.[40] Inclusion of few explicit solvent molecules usually changes the sign of binding energy.[21,40,41] Comparing the relative stabilities of triplets presented in Figure 5 with the last column of Table 2 in Ref. [21] one can note that the predictions of our explicit solvent model are generally the same as obtained in[21] from a simple model with one EMIM – anion pair. Our explicit solvent may be seen as an extreme case of such approach, and, as shown above, suggests dissociation of ion pairs or triplets. Relative binding energies obtained from the PCM model agree with explicit solvent data despite of the sign problem, therefore the PCM calculations seem to estimate correctly relative strength of interactions of different ion aggregates in given ionic liquid. Our results confirm that either applying of the PCM model or inclusion of one ion pair of the IL can correct the vacuum values of the binding energy and the relative stabilities of triplets are similar to that obtained from the analysis of large clusters.

Molecular Dynamics Simulations At the end we tried to confront the predictions of energy calculations in the solvent with molecular dynamics. To check the stability of ion triplets in MD we performed NVT simulations (without any restrictions on the geometry of initial triplet) of [LiAn1An2]2 triplets in EMIM-An1 liquids for all possible combinations of An1 and An2 anions. The simulation box contained only one ion triplet to avoid interactions between solvation shells of different Li1 ions and to make analysis of the trajectory easier. To collect some statistics five independent MD runs were performed for each system. First we discuss structures of selected Li1 – An2 complexes observed in MD simulations. Figure 6 presents radial distribution functions (RDF) for Li-X pairs (where X is an atom belonging to coordinating anion) and distance-dependent coordination numbers for Li-FSI in EMIM-FSI and Li-TFSI interactions in EMIM-TFSI. Data for mixed structure LiTFSI in EMIMFSI and LiBF4, LiB(CN)4 and LiPF6 in corresponding liquids are available as Figure S7 in Supporting Information. The maximum in RDF for Li-O interaction in LiTFSI and LiFSI ˚ . Two maxima appear in systems is located slightly below 2 A ˚ Li-N RDF: at about 3.5 and 4.5 A for LiTFSI or at 3.4 and 4.4 A˚ for LiFSI. The first of these peaks corresponds to a bidentate coordination of Li1 by the (T)FSI2 anion, the other to the monodentate orientation. Distance-dependent coordination numbers show that the cation is coordinated in average by 4.3 and 4.5 oxygen atoms in systems with FSI2 and TFSI2, respec˚ is tively. The average number of nitrogen atoms up to 5 A about 3.4 for LiFSI and 3.0 for LiTFSI. This suggests that instead of double-bidentate coordination to two anions, lithium cation quite often interacts with three anions, usually coordinating two oxygen atoms from one anion and one atom from each of two other anions. Such behavior was already reported for experimental data and MD simulations[42–45] on Li1 coordination in ionic liquids with FSI2 or TFSI2 anions. In some cases two anions are coordinated as

Figure 6. Li-X radial distribution functions and distance-dependent coordination numbers n(r) for Li1 interactions with TFSI2 and FSI2 anions obtained from the MD simulations. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

bidentate and the third as monodentate increasing the average number of O atoms in the first coordination shell above 4. Similar behavior is observed for mixed Li1-FSI2-TFSI2 interactions (Fig. S7). TFSI2 anions coordinating Li1 ion in our MD simulations prefer cis conformation; cis and trans conformers of anions are more equally populated in the case of LiFSI. Preferential cis conformation of TFSI2 anion has been reported in Raman spectroscopy and DFT studies of Li1 solvation in TFSIcontaining ILs.[46,47] Conversely, analysis in Ref. [18] suggests mixed cis/trans conformation of TFSI2 in the vicinity of lithium cation. However, the energy differences between the two conformations of LiTFSI are about 1 kcal/mol[46,47] and therefore our MD simulations may not capture well the ratio of the two conformers, because force field was not parametrized to reproduce such a small effect. The first maximum in Li-F RDFs for Li-BF4 and Li-PF6 interactions is located about 1.9 A˚ (Supporting Information Fig. S7). ˚ is about 4.9 and 4.0 Integrated number of F atoms up to 2.5 A 2 2 for BF4 and PF6 , respectively. The average number of B or P ˚ from the Li1 cation is 3.2 and 3.5, suggesting atoms up to 4 A 2 that most PF6 anions are coordinated as monodentate while the probability of bidentate coordination is significantly larger for BF2 4 (though monodentate orientation is populated as Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

9

FULL PAPER

WWW.C-CHEM.ORG

Figure 7. Selected snapshots from MD trajectories for [Li(TFSI)(B(CN)4)]– in EMIM-TFSI a–d), [Li(TFSI)(BF4)]2 in EMIM-TFSI e–f ) and [Li(FSI)(TFSI)]2 in EMIM-FSI g–h).

well). Accordingly, in Li-B RDF for Li-BF4 the peak at about 2.5 ˚ corresponding to bidentate coordination is higher than the A ˚ , resulting from monodentate configuration. other above 3 A Stability of ion triplets in MD simulations was checked by visual inspection of the first solvation shell of the Li1 cation. Ion triplet was considered unstable if we observed separation of the An22 (salt) anion from the lithium ion lasting more than 0.1 ns; in such a case the simulation was terminated. If no long-enough separation event occurred during 10 ns, [LiAn1An2]2 triplet was counted as stable. In Figure 7, three typical examples are presented. For some salt/IL combinations we observed that initial mixed triplets disintegrated after relatively short time: salt anion (An2) coordinating the Li1 ion was exchanged for the anion of the 10

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

liquid (An1) and the homogenous [Li(An1)2]2 was formed. An example is [Li(TFSI)(B(CN)4)]2 triplet in EMIM-TFSI liquid shown in Figure 7a-d. After about 0.5 ns the BðCNÞ2 4 ion separated from Li1 and the lithium cation is coordinated to 4 oxygen atoms from TFSI2 anions. Figures 7ef and 7gh present two cases of stable mixed triplets: [Li(TFSI)(BF4)]2 in EMIM-TFSI and [Li(FSI)(TFSI)]2 in EMIM-FSI. In both cases after 10 ns of MD simulation the Li1 ion is coordinated to the An22 anion from the lithium salt indicating stability of the [LiAn1An2]2 triplet. We may note, however, that it does not necessary mean that there is no exchange of ions – in the second example (Fig. 7h) at the end of the trajectory the triplet is still composed of FSI2 and TFSI2 anions but the FSI2 ion originally bound to Li1 has been replaced by another FSI2 ion from the liquid. Therefore, even “stable” triplets, that is, preserving original composition, may dynamically exchange ions. This behavior is observed for almost all of stable triplets, exceptions are mainly systems in EMIM-TFSI, for example, for [Li(TFSI)(BF4)]2 in EMIM-TFSI only in one of five trajectories an exchange of the TFSI2 ion was observed during 10 ns of simulation. Dynamical changes in the triplet structure agree with positive binding energies of triplets from explicit solvent calculations indicating that these aggregates tend to break. With this respect explicit solvent model performs better than PCM model yielding negative binding energies for most systems over almost whole range of values of static dielectric constant. Summary of the triplet stability observed in MD simulations is presented in Table 3. Average life-time is given for unstable triplets. A quick comparison with the predictions of stability obtained from FF calculations in vacuum (Table 1) reveals a general (and rather surprising) agreement between both sets of results. All mixed triplets were stable in EMIM-B(CN)4 during 10 ns MD run because exchange of the salt anion for weakly binding BðCNÞ2 4 is energetically unfavorable. Conversely, in EMIM-BF4 homogenous [Li(BF4)2]2 triplet is the most stable species and mixed triplets decayed in average during first 2 ns of simulation. Turning on to the triplets classified in Table 1 as unclear we may see that there is noticeable preference for 2 binding the PF2 and TFSI2. This 6 anion in systems with FSI may be partially related to the force field used in simulations – relative stability of triplets with PF2 6 obtained in vacuum from FF is slightly larger than resulting from QC calculations. For all systems in which mixed triplets changed into homogenous triplets in all 5 independent realizations of MD simulation (red entries in Table 3), average decay time of the triplet was relatively short—no longer than 2 ns. This indicates that for low concentrations of lithium salt in these systems the Li1 is most likely to form triplets with two identical anions provided by the ionic liquid. Conversely, in salt/IL combinations marked deep green in Table 3 mixed triplet was stable. Although the length of simulation was limited to 10 ns, thus one can not exclude possibility that changes may occur at longer times, repeatability of the result together with vacuum stabilization energies suggest the mixed triplet as the probable charge carrier in these systems. In three remaining cases with PF2 6 , where only in a fraction of simulated systems WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

Table 3. Summary of stability of [LiAn1An2]2 ion triplets in EMIM-An1 ionic liquid obtained from five independent MD simulations for each system. IL anion Salt anion

BðCNÞ2 4

BðCNÞ2 4 PF2 6 FSI TFSI BF2 4

PF2 6

FSI2

TFSI2

1.0 (0.5)

1.4 (0.9) 6.3

1.0 (0.2) 3.4 (3.5) 2.1 (1.4)

1.0 (0.3) 4.6 (2.9)

BF2 4 0.8 1.6 1.5 1.7

(0.3) (0.7) (1.1) (1.0)

0

1

2

3

not happen for the energies of triplets with TFSI2 and FSI2 which remain too low because of the idealized geometry. We did not try to calculate binding energies corresponding to structures seen in MD simulations, because the fourfold coordination of Li1 is reached by interactions either with two or three anions thus the energies are not directly comparable and the calculation scheme in explicit solvent would be more complicated. We conclude that the differences between MD simulations and predictions of explicit solvent model can be attributed to the structure of studied systems.

Conclusions

Legend No. of separation events

FULL PAPER

4

5

Colors decode the number of observed separations of the An2 anion from the triplet. Numbers are the average lifetimes of unstable triplets and their standard deviations (in parentheses) in nanoseconds.

coordination of Li1 has changed, one may expect that there will be a population of both kinds of triplets (mixed and homogenous) containing Li1 cation. It is however obvious that the results of MD simulations do not agree with relative binding energies calculated in the solvent (either implicit or explicit). According to calculated binding energies, triplets with FSI2 and TFSI2 anions should be 2 the most stable, binding with BF2 4 is less preferable and PF6 is the second last. This is not observed in simulations. While one may argue that the stability of [Li(FSI)(BF4)]2 in EMIM-FSI and [Li(TFSI)(BF4)]2 in EMIM-TFSI is an artifact caused by too short time of simulations, fast conversion of [Li(BF4)(FSI)]2 and [Li(BF)4(TFSI)]2 in EMIM-BF4 into [Li(BF4)2]2 suggests that bind2 ing to BF2 or TFSI2. 4 is preferred over interaction with FSI 2 Likewise, [Li(PF6)(FSI)] in EMIM-PF6 decays quite fast and the [Li(PF6)2]2 is formed. As mentioned earlier, MD simulation results are closer to relative binding energies in vacuum rather than to the values obtained in the solvent. Nevertheless we think that the explicit solvent model gives in principle better description of stability of ion aggregates, but it may fail when the geometry of the system used in energy calculations differs too much from the structure in “real” solvent. In our calculations in explicit solvent we used ion triplets fixed in the geometry of the optimum obtained from quantum-chemical calculations. Li1 was coordinated to FSI2 or TFSI2 in bidentate geometry (and double-bidentate when it was complexed with two anions of such type). As explained at the beginning of this section in “unfrozen” MD simulations we observed that Li1 cation usually interacts with three FSI2 or TFSI2 anions with mixed bi- and monodentate orientation. Such coordination leads to smaller stabilization per Li1-O interaction, therefore average strength of cation–anion interaction for real geometry is smaller than obtained for our ideal, frozen structure. The probable origin of the discrepancy is that in explicit solvent binding energy of the triplet with BF2 4 gets closer to binding energy of the triplet containing BðCNÞ2 4 but this does

We have analyzed relative stabilities of ion aggregates in five ionic liquids using quantum-chemical calculations supported by Molecular Dynamics simulations. Stability pattern of ion triplets in vacuum results mainly from the stabilities of LiAn pairs, with strongest binding of BF2 4 anion and weakest interaction with BðCNÞ2 . Explicit solvent force field calculations per4 formed for large solvent clusters yield different relative stabilities: triplets containing TFSI2 or FSI2 anions are more stable than triplets with tetrafluoroborate. Explicit solvent model proved also that reorganization of the liquid in the field of the solute leads to variations of interaction energies observable at distances up to 2 nm. Binding energies obtained from the Polarizable Continuum Model are in most cases negative (as opposed to the values in explicit liquid), thus overestimating stability of ion pairs and triplets, but the relative stabilities of ion triplets in ionic liquids predicted by both approaches are in general agreement and capture most of the difference with respect to vacuum calculations. Solute–solvent interactions calculated in continuous model are significantly weaker than in explicit solvent and the negative binding energies in PCM can be largely attributed to smaller stabilization of Li1 in the liquid. Stabilities of ion triplets observed in Molecular Dynamics for triplets with TFSI2 or FSI2 do not agree with the predictions of solvent models. In MD simulations apparently most stable are the systems with BF2 4 anion. We suggest that this discrepancy results from the structure of the triplets used in explicit solvent calculations frozen at the optimal geometry, while the real coordination of Li1 ion is different. We conclude therefore that calculations of binding energies of ions in ionic liquids within explicit model require large system sizes because of the nature of long-range correlations in the structure of the liquid. Conversely, estimates of relative values obtained in this work agree well with predictions of PCM approach or results obtained in[21] for one explicit pair of IL ions. This suggests that computationally cheap quantumchemical calculations with implicit solvent or QC computations with explicit solvent reduced to one pair of ionic liquid ions can provide quite reliable estimates of relative binding energies. However, absolute values obtained from the two approaches differ greatly and may even have opposite sign. Moreover, results of quantum-chemical computations for optimal geometries of isolated ion aggregates should be always taken with care—for some systems simulations of dynamics Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

11

FULL PAPER

WWW.C-CHEM.ORG

do not agree with static calculations. We hope that these results and conclusions will provide some more insight into ion–ion interactions in ionic liquids and will help in proper choice of modeling approach in future research. Keywords: ionic liquids  ion aggregates  explicit solvent model  continuous solvent

How to cite this article: A. Eilmes, P. Kubisiak J. Comput. Chem. 2015, DOI: 10.1002/jcc.23853

]

Additional Supporting Information may be found in the online version of this article.

[1] J. P. Hallett, T. Welton, Chem. Rev. 2011, 111, 3508. [2] A. Rehman, X. Zeng, Acc. Chem. Res. 2012, 45, 1667. [3] S. J. Zhang, J. Sun, X.C. Zhang, J.Y. Xin, Q.Q. Miao, J.J. Wang, Chem. Soc. Rev. 2014, 43, 7838. [4] M. V. Fedorov, A.A .Kornyshev, Chem. Rev. 2014, 114, 2978. [5] C. Chiappe, C.S. Pomelli, Eur. J. Org. Chem. 2014, 28, 6120. [6] H. Weing€artner, Angew. Chem. Int. Ed. 2008, 47, 654. [7] E. J. Maginn, Acc. Chem. Res. 2007, 40, 1200. [8] Y. Wang, W. Jiang, T. Yan, G.A. Voth, Acc. Chem. Res. 2007, 40, 1193. [9] B. Kirchner, Top. Curr. Chem. 2009, 290, 213. [10] E. I. Izgorodina, Phys. Chem. Chem. Phys. 2011, 13, 4189. [11] S. Zahn, D. R. MacFarlane, E. I. Izgorodina, Phys. Chem. Chem. Phys. 2013, 15, 13664. €rmann, J. Chem. Soc. Perkin. Trans. 2, 1993, 5, 799. [12] A. Klamt, G. Sch€ uu [13] S. Miertus, E. Scrocco, J. Tomasi, Chem. Phys. 1981, 55, 117. [14] M. N. Kobrak, H. Li, Phys. Chem. Chem. Phys. 2010, 12, 1922. [15] C. Wakai, A. Oleinikova, M. Ott, H. Weing€artner, J. Phys. Chem. B 2005, 109, 17028. [16] V. S. Bernales, A. V. Marenich, R. Contreras, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. B. 2012, 116, 9122. [17] K. Angenendt, P. Johansson, J. Phys. Chem. C 2010, 114, 20577. [18] J.-C. Lasse`gues, J. Grondin, C. Aupetit, P. Johansson, J. Phys. Chem. A 2009, 113, 305. [19] A. M. Fernandes, M.A.A. Rocha, M. G. Freire, I. M. Marrucho, J. A. P. Coutinho, L. M. N. B. F. Santos, J. Phys. Chem. B 2011, 115, 4033. [20] Y. Umebayashi, H. Hamano, S. Seki, B. Minofar, K. Fujii, K. Hayamizu, S. Tsuzuki, Y. Kameda, S. Kohara, M. Watanabe, J. Phys. Chem. B 2011, 115, 12179. [21] K. Angenendt, P. Johansson, J. Phys. Chem. B 2011, 115, 7808. [22] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J.

12

Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23853

[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]

Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford CT, 2009. TeraChem v. 1.5, PetaChem, LLC, Available at: http://www.petachem. com/. Accessed February 9, 2015. O. Borodin, J. Phys. Chem. B 2009, 113, 11463. J. N. C. Lopes, J. Deschampes, A.A.H Padua, J Phys. Chem. B 2004, 108, 2038. J. N. C. Lopes, A.A.H Padua, J. Phys. Chem. B 2004, 108, 16893. L. Martınez, R. Andrade, E.G. Birgin, J.M. Martınez, J. Comput. Chem. 2009, 30, 2157. Tinker Molecular Modeling Package v. 5.1, Available at: http://dasher. wustl.edu/tinker/. Accessed February 9, 2015. B. T. Thole, Chem. Phys. 1981, 59, 341. A. Eilmes, M. Sterzel, T. Szepieniec, J. Kocot, K. Noga, M. Golik, In eScience on Distributed Computing Infrastructure; M. Bubak, J. Kitowski, K. Wiatr, Eds.; Lecture Notes in Computer Science, v. 8500, Springer International Publishing, Switzerland, 2014; pp. 250–263. R. M. Lynden-Bell, J. Chem. Phys. 2008, 129, 204503. B. L. Bhargava, M. L. Klein, S. Balasubramanian, ChemPhysChem 2008, 9, 67. R. M. Lynden-Bell, Phys. Chem. Chem. Phys. 2010, 12, 1733. R. M. Lynden-Bell, J. Phys. Chem. B 2007, 111, 10800. A. V. Marenich, C.J. Cramer, D. G. Truhlar, J. Phys. Chem. B 2009, 113, 6378. K. Nakamura, T. Shikata, ChemPhysChem 2010, 11, 285. M.-M. Huang, Y. Jiang, P. Sasisanker, G. W. Driver, H. Weing€artner, J. Chem. Eng. Data 2011, 56, 1494. X. Hu, Q. Lin, J. Gao, Y. Wu, Z. Zhang, Chem. Phys. Lett. 2011, 516, 35. H. Roohi, S. Khyrkhah, Comp. Theor. Chem. 2014, 1037, 70. P. Johansson, Phys. Chem. Chem. Phys. 2007, 9, 1493. A. Eilmes, P. Kubisiak, Electrochim. Acta 2011, 56, 3219. Z. Li, G. D. Smith, D. Bedrov, J. Phys. Chem. B 2012, 116, 12801. H. Liu, E. Maginn, J. Chem. Phys. 2013, 139, 114508. K. Fujii, H. Hamano, H. Doi, X. Song, S. Tsuzuki, K. Hayamizu, S. Seki, Y. Kameda, K. Dokko, M. Watanabe, Y. Umebayashi, J. Phys. Chem. C 2013, 117, 19314. C.J.F. Solano, S. Jeremias, E. Paillard, D. Beljonne, R. Lazzaroni, J. Chem. Phys. 2013, 139, 034502. Y. Umebayashi, T. Mitsugi, S. Fukuda, T. Fujimori, K. Fujii, R. Kanazaki, M. Takeuchi, S. Ishiguro, J. Phys. Chem. B 2007, 111, 13028. Y. Umebayashi, S. Mori, K. Fujii, S. Tsuzuki, S. Seki, K. Hayamizu, S. Ishiguro, J. Phys. Chem. B 2010, 114, 6513.

Received: 9 December 2014 Revised: 14 January 2015 Accepted: 17 January 2015 Published online on 00 Month 2015

WWW.CHEMISTRYVIEWS.COM