Physical Nature of Intermolecular Interactions in [BMIM][PF6] Ionic

Feb 3, 2014 - Soraya Ebrahimi , Mohammad H. Kowsari ... Motahare Sadeghi Googheri , Mohammad Sadegh Sadeghi Googheri , Samira Hozhabr Araghi...
0 downloads 0 Views 344KB Size
Article pubs.acs.org/JPCB

Physical Nature of Intermolecular Interactions in [BMIM][PF6] Ionic Liquid Borys Szefczyk* and W. Andrzej Sokalski Institute of Physical and Theoretical Chemistry, Faculty of Chemistry, Wroclaw University of Technology, Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland ABSTRACT: The intermolecular interaction energy in a popular ionic liquid, [BMIM][PF6] is analyzed using the Hybrid Variation-Perturbation Theory approach. The analysis is performed on a sample of configurations from molecular dynamics simulation, instead of minimized structures. The interaction energy components are quantified, showing that the electrostatics is the dominating but not the only important term. It is found that two- and threebody electron delocalization components also contribute to the stabilization of the complexes; however, these interactions vanish beyond the first coordination sphere. The presented study shows a systematic way to obtain the amount of physically meaningful components of the interaction energy, which possibly could be related to macroscopic properties of ionic liquids (e.g., viscosity, melting point) or electron transfer in ionic liquids.



possible to find a relation between measured macroscopic properties and molecular properties calculated from first principles. In order to analyze intermolecular interactions in a systematic way and to assign them a physical meaning, several energy decomposition schemes have been devised. The prominent place among them is held by the SymmetryAdapted Perturbation Theory (SAPT).10,11 However, other schemes are used as well, for example, the Morokuma-Kitaura12 scheme or the Hybrid Variation-Perturbation Theory (HVPT) decomposition scheme.13,14 The latter one has been successfully applied to elucidate the intermolecular interactions (including many-body effects) in a series of large molecular complexes15−20 and will be used here. In the context of ionic liquids, several reports exist dealing with the interaction of an isolated pair at ab initio and density functional theory level7,21 also including the [BMIM][PF6] ionic liquid investigated here.1 Zahn et al. used the SAPT approach to evaluate the interactions in [MMIM][Cl].22 However, to the best of our knowledge, the many-body and cooperative effects are still not fully examined, except for a single work by Kossmann et al.23 They investigated cooperative effects in [MMIM][Cl], but without looking into the physical nature of the interactions (especially the nonadditive manybody effects). The present work aims at filling this gap, by calculating the interaction energy terms in a systematic way and by investigating the physical source of the cooperative (nonadditive) effects. The HVPT decomposition scheme is used to achieve this goal. Also, as opposed to the bottom-up

INTRODUCTION In the quest for sustainable and environmentally friendly technologies, one of the routes explored is replacement of “traditional” organic solvent with ionic liquids (IL). Formerly called molten salts, ILs are composed of ions, but contrary to inorganic salts, they are liquid at ambient temperature (below 100 °C). Due to the strong interaction between the ions (order of magnitude larger than in water),1,2 they are almost nonvolatile3 and may exert unexpected catalytic activity.4 Extremely low vapor pressure together with the broad range of anions that can be combined together to form an IL render them an attractive alternatives to organic solvents. This is despite the fact that ILs are more expensive; however, the idea is to recycle IL several times.5 An already existing library of possible cations and anions allows us to create a multitude of ILs. The choice of the ions determines the properties of the liquid, including the melting point (some are very viscous or even solid at room temperature).6 Consequently, there is great interest in the rational design and prediction of the physical properties of an IL, based on the choice of cation and anion. It has been shown that several properties of IL depend on the components of the interaction energy in cation−anion pairs. For instance, Bernard et al. have shown that the conductivity, viscosity, and melting point of ILs can be related to the dispersion interaction or the ratio of the total and dispersion interaction.7 Hayes et al. used neutron diffraction to show that the melting point correlates with the number, length, and direction of the hydrogen bonds formed in the cation−anion (CA) pair.8,9 It cannot be expected that the analysis of interactions between isolated ions will be sufficient to design IL with desired properties, due to the highly dynamic nature of the liquid and many effects (e.g., entropic) that cannot be captured that way. Nevertheless, for certain classes of IL it might be © 2014 American Chemical Society

Received: November 19, 2013 Revised: January 27, 2014 Published: February 3, 2014 2147

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

approach used by Kossmann et al.23 (they build up a complex by stacking monomers), here a top-down approach is applied: a classical MD is used to obtain snapshots and statistical analysis of the interactions is performed at high quantum-chemical level. In our study, we choose to investigate the IL composed of butylmethylimidazolium cation and PF6− anion, [BMIM][PF6]. This is a popular IL, often studied by means of experimental24,25 and theoretical methods.26,27 We analyze the components of interaction energy in this liquid, paying special attention to their amount as a function of the distance. This allows us to understand the difference in the interaction within and beyond the first coordination shell. This analysis is complemented by the study of nonadditive components of the interaction energy.

ΔEHL =

This interaction, thanks to the antisymmetrisation operator  , satisfies the Pauli exclusion principle. It can be further decomposed into the terms HL ΔEHL = ∈(10) +∈ex el (10) where ∈HL ex is the electron exchange component and ∈el is the electrostatic component calculated as the first-order perturbation in the polarization expansion. The latter term is identical in SAPT and HVPT approaches. Finally, one can use one of several atomic multipole expansions (DMA,30 CAMM,31 etc.) to obtain various components of the electrostatic interaction, like charge−charge, charge−dipole, dipole−dipole, and so on. The use of atomic instead of molecular multipole expansion allows us to obtain convergent results at intermolecular distances closer to equilibrium geometry. Distributed Multipole Analysis (DMA)30 is used in current HVPT implementation.16,28 The moment truncated sum of such component, called here the multipole interaction, ∈(1) el,MTP can be subtracted from the electrostatic term ∈(10) el , giving the electron penetration term



THEORETICAL BACKGROUND The interaction energy ΔE of an n-body complex is defined as n

ΔE(Ξ) = E123...n(Ξ) −

∑ Ei(Ξ) i=1

where E123...n(Ξ) is the total energy of the n-body complex and Ei(Ξ) is the energy of ith monomer in the unaltered geometry (Ξ) of the whole complex. The interaction energy is composed of additive (two-body) as well as the nonadditive (three-, fourbody, etc.) terms: ΔE =

(10) ∈(1) −∈(1) el,PEN = ∈el el,MTP

∑ ΔEi(12) + ∑ ΔEi(123) + ... + ΔE(123...n) i

̂ |Â ΨAΨB Â ΨAΨB|HAB Ψ |Ĥ |Ψ − A A A ΨA|ΨA Â ΨAΨB|Â ΨAΨB Ψ |Ĥ |Ψ − B B B ΨB|ΨB

i

If the atomic multipole series is truncated right after the monopole interaction, one obtains the Coulomb formula for the interaction of static point charges used by most force fields. This kind of interactions will be denoted as ΔEqq. Finally, at larger distances the Coulomb interaction can be “coarsegrained”  one can use simply single point molecular charge for each monomer (ΔEQQ). With the formal charges +Q and −Q for cation and anion, respectively, the interaction energy of a single pair is ΔEQQ = −Q2r−1. These interaction terms form the following hierarchy of approximations to the total energy of interaction:

Currently, the most rigorous and systematic scheme for describing the intermolecular interaction energies is probably the Symmetry-Adapted Perturbation Theory (SAPT).10,11 In this work, an alternative scheme is used, namely, the variationalperturbational energy decomposition scheme.13,14 The procedure has been implemented in the GAMESS-US package and was described in several papers.16,28 The main advantage of this scheme is the possibility to perform calculations of larger molecular systems; here, systems of up to 64 atoms are considered. According to this scheme, the interaction energy is decomposed into following terms:

HL HF (1) ΔE MP2 = ∈(1) el,MTP + ∈el,PEN + ∈ex +ΔEdel + ΔEcorr

HL HF (1) ΔEHF = ∈(1) el,MTP + ∈el,PEN + ∈ex +ΔEdel

ΔE ≈ ΔE MP2 HL HF (1) = ∈(1) el,MTP + ∈el, PEN + ∈ex +ΔEdel + ΔEcorr

HL (1) ΔEHL = ∈(1) el,MTP + ∈el,PEN +∈ex

The ΔEMP2 is the interaction energy calculated at the Møller− Plesset second level of theory. The electron correlation component, ΔEcorr, is nonadditive and is computed as a difference:

(1) ∈(10) = ∈(1) el el,MTP +∈el,PEN

ΔE qq =

∑ ∑ (qa·qb)rab−1 a

ΔEcorr = ΔE MP2 − ΔEHF

b

ΔEQQ = −Q 2r −1

where ΔEHF is the Hartree−Fock interaction energy. The electron delocalization term, ΔEHF del , can be interpreted as the interaction due to the charge redistribution, upon formation of the n-body complex from unperturbed charge distributions of the monomers. This term is nonadditive as well, and is calculated as the difference

Among them, only the multipole, penetration, and electrostatic interaction terms are additive, i.e., do not contribute to the many-body interactions. Since the ΔEcorr incorporates besides dispersion various physical phenomena, like the intramolecular correlation correction to the electrostatic interaction, further parts of electron exchange interaction and part of the dispersion interaction, it could be educative to estimate their amount. In this work, the second-order ∈(20) disp interaction is calculated in order to give a better estimate of the role of classical dispersion

HF ΔEdel = ΔEHF − ΔEHL

where ΔEHL is the Heitler-London therm, calculated according to the Löwdin formula:29 2148

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

counter-poise (CP) correction, in order to remove the basis set superposition error (BSSE).33 Two kinds of interaction energy calculations have been done. The first kind of calculations was aimed at computing the twobody interactions only, but for a large number of configurations. We have used 16 distinct snapshots from the MD simulations. Each snapshot contained a central ion pair (cation +anion) surrounded by 31 nearest cation+anion pairs. From each snapshot, the 32 configurations have been extracted, each one containing the central cation and one anion. In this way, we have obtained 16 × 32 configurations of a single cation+anion pair, differing in the distance between the ions. The distances were measured between the phosphorus atom in PF6− and geometrical center of imidazolium ring, spanning the range from 4.18 to 20.29 Å. The second kind of calculation was aimed at estimating the significance of many-body interactions. In this case, 4-mer complexes were analyzed, i.e., complexes of two cations and two anions. The same snapshots were used and the complexes were chosen in the following way: the central pair (cation +anion) was present in all complexes and the other cation +anion was selected from the first seven nearest pairs. By the “nearest pair” we mean the nearest cation and the nearest anion, which do not necessarily form a tightly bound pair (they may be located on the opposite sides of the central pair, for example). This choice is illustrated in Figure 1. In the case of

term in the investigated liquids. For that purpose, the following formula is used:32 ∈(20) disp (IJ ) = 4 ∑ ∑ ∑ ∑ i∈I j∈J a∈I b∈J

ij|ab ab|ij εi + εj − εa − εb

where i, j and a, b are occupied and virtual molecular orbitals, respectively, and ε refers to respective orbital energies. I, J indicate monomers. (20) Besides the additive ∈ disp component, the electron correlation interaction also contains second-order contribution to the exchange−delocalization interaction, ΔE(2) ex−del, which is nonadditive. It is calculated as a difference:16 MP2 (2) ΔEex − ΔE(12)(el) − ∈(20) disp − del = ΔE

where ΔE(12)(el) is the correlation correction to the electrostatic interaction. When considering a general system composed of n monomers (more than two) one has to be careful with the definition of the pair interactions. The total energy of interaction ΔE contains n(n − 1)/2 mutual (pairwise) interactions. The general formula for the number of k-body interactions of an n-mer is given by the binomial coefficient: ⎛n⎞ n! ⎜ ⎟ = ⎝k ⎠ k ! (n − k )!

The second issue that has to be always kept in mind is to compare like with like: for instance, in order to compare twobody interactions in a dimer and in a 4-mer, one should perform all calculations in the 4-mer basis set. Such an approach allows us to avoid the basis set superposition error (BSSE).33 In this respect, SAPT as well as the HVPT decomposition scheme provides elegant and systematic approach, because all the many-body components are automatically included and the BSSE corrections are properly handled.



METHODS The geometries used in calculations were obtained from classical molecular dynamics simulations (snapshots kindly provided by Dr. Gyorgy Hantal).34 Sixteen configurations were selected every 250 ps from a 5 ns NpT simulation at 298.15 K, using Parinello-Rahman barostat (pressure maintained at 1 bar) and velocity rescaling thermostat; the force field used for simulations was developed by Bhargava and Balasubramanian.35 All calculations of the interaction energy and its components have been done in a modified version of the GAMESS-US package.36 For most calculations a smaller basis set was used, namely, Lanl2DZ, with the core orbitals replaced by Effective Core Potentials (ECP).37 For a subset of 16 cation+anion configurations, calculations in larger basis sets have been performed as well, in order to check if the Lanl2DZ/ECP basis set offers sufficient reliability. These benchmark calculations have been done in cc-pVDZ and cc-pVTZ basis sets. Since the system contains negatively charged PF6− anions, tests were also performed with basis sets including diffuse functions; the Lanl2DZ/ECP basis set augmented with additional polarization and diffuse functions (hence, denoted Lanl2DZdp/ECP),38 and aug-cc-pVDZ, which is similar to cc-pVDZ, except that it also includes additional diffuse functions. All calculations of the interaction energy described in this paper are done with

Figure 1. Sample configuration used in the interaction energy calculations. The central pair C0:A0 is indicated in red. The choice for two-body calculations would be C0:A0, C0:A1, C0:A2, etc. The choice for four-body calculations would be C0:A0:C1:A1 (P1), C0:A0:C2:A2 (P2), etc.

the four-body interactions, only up to the seventh nearest pair was considered (designated P1, P2, ..., P7), i.e., up to the distance of ca. 8−9 Å, the first coordination shell. This is because the only components of the interaction energy that are nonadditive are also vanishing quickly, due to the dependence on the higher powers of 1/r (r being the distance). The nonadditive components have been indicated in the text, tables, and graphics with subscript, e.g., ∈HL ex,3 means three-body exchange component. The paper aims at investigating interactions at room conditions; therefore, the geometries extracted from the MD snapshots were not optimized at the quantum chemical level. Additionally, a single CA pair was optimized at the DFT level using B3LYP functional and aug-cc-pVDZ basis set and at the MP2 level using Lanl2DZdp basis with ECP. The analysis of 2149

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

Table 1. Energy of Interactions [kJ mol−1] of a Single CA Pair Calculated Using Different Basis Sets and Averaged over 16 Snapshotsa Component ∈(10) el ΔEHL ΔEHF ΔEMP2 ∈HL ex ΔEHF del ΔEcorr ∈(20) disp a

Lanl2DZ (ECP)

Lanl2DZdp (ECP)

cc-pVDZ

aug-cc-pVDZ

cc-pVTZ

−280. 14.4 −264. 5.9 −279. 8.1 −284. 7.1 16.65 11.94 −15.64 3.73 −4.97 2.49 −7.02 2.06

−280. 13.5 −262. 5.3 −282. 7.4 −291. 7.3 17.30 12.37 −19.23 4.25 −9.61 3.30 −13.93 3.99

−283. 14.7 −265. 6.0 −283. 8.7 −292. 8.6 17.82 12.46 −17.98 4.37 −8.87 3.00 −11.33 3.40

−280. 13.6 −262. 5.2 −283. 7.7 −298. 8.4 17.56 12.51 −21.03 4.61 −14.73 4.27 −19.87 5.43

−281. 13.7 −264. 5.4 −284. 7.7   17.52 12.49 −20.28 4.51    

Every second row contains the standard deviation (σ) of the above interaction term calculated among the snapshots.

The Lanl2DZdp/ECP basis set, which contains diffuse functions, is able to recover most of the delocalization interaction (91% compared against aug-cc-pVDZ), but the electron dispersion is still underestimated (70% of ∈(20) disp from aug-cc-pVDZ). These figures have to be taken into account in the following discussion that relies on data obtained using Lanl2DZ/ECP and Lanl2DZdp/ECP basis sets. Similar analysis has been extended for other possible CA pairs, extracted from the 16 snapshots. For each snapshot, the central cation was selected, and for each anion in this snapshot, the interaction energy components have been calculated. Subsequently, all results were plotted versus the distance between P in PF6− and the geometrical center of imidazolium ring. The resulting interaction energy at the MP2 level is presented as black dots in Figure 2. The interactions at other levels of theory have been calculated as well and, after smoothing with Bezier curve, were plotted as solid lines in Figure 2.

interaction energy components was subsequently performed for the optimized geometry, using two basis sets: Lanl2DZ/ECP and Lanl2DZdp/ECP.



RESULTS From 16 snapshots, a closely bound CA pair was selected and the energy decomposition was performed. Five basis sets were used: Lanl2DZ and Lanl2DZdp with ECP and three correlation-consistent basis sets, cc-pVDZ, aug-cc-pVDZ, and cc-pVTZ. These calculations were aimed at verifying the smaller basis set (Lanl2DZ) and the ECP’s, which were necessary due to the large number of complexes that were considered in this work. Table 1 shows the interaction energy calculated with these basis sets; except for cc-pVTZ, the calculations were carried out up to the MP2 level of theory. Not surprisingly, the major component of the interaction in these ionic pairs is the electrostatics. It amounts to −280 kJ· mol−1, which is 94% of the total energy, −298 kJ·mol−1 at the MP2/aug-cc-pVDZ level. However, due to the short intermolecular contacts, there is some contribution of the interaction energy components that depend on higher powers of 1/r: the electron correlation, about −15 kJ·mol−1 (5% of ΔE), delocalization, ca. −21 kJ·mol−1 (7% of ΔE), and exchange (about 6% of ΔE). The amount of dispersion interaction is ca. −20 kJ·mol−1 (7% of ΔE)  more than ΔEcorr, since it is partially canceled by the electron exchange contribution to the correlation energy. Despite the fact that the Heitler-London interaction is substantially lower than ΔEHF and ΔEMP2, the purely electrostatic component ∈(10) is very el close to the HF energy. This is due to the fact that electron delocalization ΔEHF del is always negative, while the exchange component ∈HL ex is always positive. Since the magnitude of these components is usually similar, they cancel each other out to a certain extent. The comparison of various basis sets shows that the Lanl2DZ/ECP combination provides a reasonable estimate of the total energy; however, the dispersion interactions (and therefore the electron correlation component) are underestimated. The Lanl2DZ basis set accounts for only 35% of the dispersion interaction, while the electron delocalization is reproduced in ca. 74% (i.e., ΔEHF del amounts to −16 and −21 kJ· mol−1, when basis sets Lanl2DZ and aug-cc-pVDZ are used).

Figure 2. Two-body energies of interactions at various levels of theory, calculated with Lanl2DZdp/ECP basis set. Black dots represent MP2 values; lines show interaction energies interpolated using Bezier curves. Dashed line shows the ΔEQQ interaction energy for |Q| = 1. 2150

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

mol−1 for the first neighbor, to −433 kJ mol−1 for the seventh one. The energies obtained for complexes with fourth to seventh neighbor are virtually the same (taking into account the error bars), because these neighbors are within the same coordination shell and the distance to the central pair is similar. However, one can notice that the sum of two-body interactions is not equal to the total energy of interaction; for instance, in the first complex, the total energy is −564 kJ mol−1, while the sum of two-body interactions is −585 kJ −1, i.e., slightly overestimated. The discrepancy is due to the nonadditive (cooperative) terms, electron exchange, delocalization, and correlation. Table 3 shows the amount of two-, three-, and four-body interactions. Although the four-body interactions are negligible, some of the three-body components have quite significant values. Namely, the delocalization component is ca. 21 kJ mol−1 and goes down to 6 kJ mol−1 for the seventh neighboring pair. The three-body electron correlation and exchange have rather insignificant magnitudes. These data are summarized in a concise way in Figure 4. It shows that the three-body delocalization is similar in magnitude to two-body electron correlation and about one-third of two-body exchange and delocalization. This is in line with the results obtained for ionic crystal, 3-pentenenitrile, 2-nitro-5-oxo, ion(−1), sodium.17

The data are noisy, especially until ca. 8−9 Å (first coordination shell), due to the fact that the cation is far from being spherical, and therefore the interaction depends strongly on the direction. Once again, in this graph one can observe that the electrostatic component agrees very well with the most accurate MP2 results. Interesting conclusions can be drawn from comparison with the r−1 curve (black dashed curve in Figure 2), which represents the interaction of two opposite unit charges in atomic units. From ca. 10 Å on, i.e., beyond the first coordination shell, the exact interaction can be described very well with this simple model. Naturally, this graph refers to two interacting monomers isolated from the rest of the system. The interactions with the environment are neglected here and the electric permittivity of vacuum is assumed. The empirical procedure of scaling down the charges to 0.8−0.9 would shift the dashed curve up, but it would not mean that the procedure is wrong. It would only be wrong for isolated monomers in vacuum, like in the case studied here. The deviation from the ideal r−1 interaction within the first coordination shell originates from the components of the interaction energy higher than electrostatic term. In Figure 3,



DISCUSSION One of the issues in the ongoing debate about the reason behind the special properties of ionic liquids is the nature of the interactions in the ionic pair. As pointed out by Zahn et al.,22 the interaction of an ionic pair in IL is weaker than in a typical salt and the nature of this interaction is that it is not completely electrostatic as one would expect. As suggested by Bernard et al.,7 the short-range interactions (e.g., dispersion) could be correlated with transport properties of ILs, such as the viscosity and diffusion coefficient. This is because the short-range interactions are the first that have to be overcome, when a molecule displaces from its original position. In context of this debate, a thorough investigation of the components of the energy of interaction in ILs emerges as a way to understand and possibly also predict the properties of IL. Quite apart from this, the study of the interaction energy may deliver answers to the technical problem of MD simulations of ILs with classical force fields: what should be the formal charge of the ions and is the scaling of the charges39 a correct approach? It is known that setting the charge of the ions to ±1 leads to underestimated dynamics of ions and various countermeasures have been undertaken: scaling down the charges,39,40 tuning the van der Waals parameters,41 or using polarizable force fields.42 When the charge transfer (or better, electron transfer, to distinguish from the charge transfer interaction) between ions forming IL is discussed, one has to rely on an arbitrary definition of the frontier between molecular fragments, partitioning the electron density between the cation and anion. Such partitioning can be done in a number of ways, for instance, based on atomic orbitals (population analyses, e.g., Mulliken population analysis43) or using topological analysis, like the one proposed by Bader (the AIM theory).44 Electron transfer values estimated from quantum chemical expectation values are closer to experimental values than values derived from any arbitrary population analysis.45 Electron delocalization, on the other hand, is much easier to define in terms of polarization interaction. The polarization is defined as charge redistribution due to the presence of the second, unperturbed monomer (first

Figure 3. Two-body components of the interaction energy  electron HF exchange ∈HL ex , delocalization ΔEdel and correlation ΔEcorr. The data were interpolated with Bezier curve.

the interpolated values of electron exchange, delocalization, and correlation components are shown. The largest and the only positive component is the electron exchange term. The largest negative component is the electron delocalization, which corresponds to the interaction due to the induced polarization (the “induced dipole” interaction). All these component vanish within the first coordination shell, because of their short-range character. The following analysis refers to energy of interactions in larger complexes of two cations and two anions (4-mers). This is in a certain way similar to the analysis proposed by Kossmann,23 in the sense that two CA pairs are considered. However, only the first (central) pair is closely bound  the other one is built from one cation and one anion of increasing distance from the central pair. Up to 7 nearest pairs are investigated in this way. Due to the size of these systems, the smaller Lanl2DZ/ECP basis set was used. Table 2 shows the total and two-body interaction in the 4mer complexes, averaged over 16 different complexes. All seven nearest neighbors are still in the first coordination shell; therefore, the interaction remains quite strong, from −564 kJ 2151

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

Table 2. Four-Body Complex of Two [BMIM][BF6] Units: Averaged Total Interaction Energy and Sum of Two-Body Interactionsa ΔEMP2

total two-body

ΔEHF

total two-body

ΔEHL

total two-body

ΔEcorr

total two-body

ΔEHF del

total two-body

P1

P2

P3

P4

P5

P6

P7

−564. 32.5 −585. 39.6 −540. 32.5 −560. 39.6 −500. 35.9 −500. 36.0 −24.1 2.63 −25. 2.63 −39.0 6.38 −59.3 5.55

−537. 38.9 −553. 44.5 −518. 37.2 −532. 42.4 −479. 37.8 −479. 37.8 −19.0 3.33 −20.0 3.52 −39.1 5.37 −53.7 6.54

−527. 31.2 −542. 34.5 −511. 31.3 −525. 34.5 −477. 32.3 −477. 32.3 −15.7 2.41 −16.6 2.48 −33.3 3.76 −47.3 3.90

−471. 31.0 −476. 33.7 −454. 30.0 −458. 32.6 −416. 30.4 −416. 30.5 −16.8 2.48 −17.1 2.57 −37.9 3.99 −42.7 4.73

−444. 47.9 −450. 53.2 −432. 46.5 −437. 51.7 −399. 47.5 −399. 47.5 −12.2 2.30 −12.5 2.40 −33.1 3.81 −38.2 4.98

−450. 52.9 −456. 57.2 −438. 52.0 −444. 56.1 −409. 52.6 −409. 52.6 −12.1 2.00 −12.4 2.15 −28.8 3.66 −35.4 4.67

−433. 38.6 −440. 42,7 −424. 38.2 −430. 42.1 −400. 37.7 −400. 37.8 −9.4 1.65 −9.6 1.74 −24.5 4.14 −30.6 5.29

a Energies in kJ·mol−1. The numbers in heading refer to the pair number (larger number means longer distance). Every second row indicates errors calculated at 95% confidence interval.

Table 3. Nonadditivity in Four-Body Complex of Two [BMIM][BF6] Units: Sum of Nonadditive Two-, Three-, and Four-Body Interaction Energy Termsa P1

P2

P3

P4

P5

P6

P7

42.00 9.79 −42.72 4.73 8.73 1.97

36.02 7.73 −38.17 4.98 7.74 1.69

33.66 7.32 −35.37 4.67 7.40 1.66

26.79 10.23 −30.65 5.29 5.92 2.13

−0.13 0.22 4.87 3.59 0.35 0.23

−0.05 0.10 5.26 5.85 0.33 0.24

−0.07 0.06 6.70 5.16 0.31 0.29

−0.02 0.04 6.15 4.40 0.28 0.23

0.00 0.00 −0.01 0.16 0.00 0.02

0.00 0.00 −0.16 0.15 −0.01 0.02

0.00 0.00 −0.10 0.16 −0.01 0.02

−0.00 0.00 −0.02 0.10 0.00 0.01

two-body ∈HL ex ΔEHF del ΔE(2) ex−del

60.52 11.92 −59.28 5.55 13.09 2.19

51.98 12.71 −53.70 6.54 11.08 2.95

44.11 9.73 −47.33 3.90 9.36 2.20

0.10 0.21 20.90 8.18 1.18 0.30

−0.22 0.27 14.79 5.98 1.07 0.39

0.10 0.12 14.13 3.95 0.91 0.18

0.00 0.00 −0.66 0.34 −0.06 0.05

0.00 0.00 −0.21 0.26 −0.03 0.04

0.00 0.00 −0.13 0.17 −0.01 0.02

three-body ∈HL ex ΔEHF del ΔE(2) ex−del

four-body ∈HL ex ΔEHF del ΔE(2) ex−del

The numbers in heading refer to the pair number (larger number means longer distance). Energies are given in kJ·mol−1. Every second row indicates errors calculated at 95% confidence interval. a

electron transfer than the electron transfer itself. Such an approach was applied based on other than HVPT decomposition schemes as well.46 As shown in Figure 3, the delocalization interaction vanishes beyond 10 Å radius, which is the size of the first coordination sphere. This does not answer the question of the amount of electron transfer, but it shows that these effects have local character  an observation which confirms findings by Schmidt et al.,47 based on comparison of

order) or as the charge redistribution due to the induced charge redistribution in the second monomer (second order). Taking the total energy at the ΔEHF level, for example, one can subtract the interaction of the unperturbed monomers (∈(10) el ) and the electron exchange (∈HL ), and the remaining part (ΔEHF ex del ) can be mostly attributed to the induction energy constituting a major part of the delocalization term. It is therefore easier to discuss the components of the interaction energy related to the 2152

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

Table 4. Energy of Interactions [kJ mol−1] of a Single, Optimized CA Complexa

Figure 4. Nonadditive components of the interaction energy in fourbody complexes of two [BMIM][BF6] units. Error bars were calculated at 95% confidence interval.

Geometry

charges calculated for isolated dimers and bulk liquid. The importance of electron delocalization was also confirmed by Prado et al.48 who have studied change in the dipole moment of cation and anion in [MMIM][Cl] due to the presence of the environment. Wendler et al. performed ab initio molecular dynamics simulations of a series of ILs and investigated the charge reduction on the ions.49 They have found that full reduction is achieved in clusters as small as eight ion pairs. Moreover, there is a correlation between dipole moment and center of charge vector in two ions that vanishes beyond 8 Å sphere. This is in line with the range of delocalization interaction presented here. The nature of interaction in IL was investigated by means of the SAPT method,1,22 however, only for single, isolated ionic pairs and rather small ones ([MMIM][Cl] and [EMIM][BF4]). Dong et al. report that only 71% of ΔE in [EMIM][BF4] is due to the electrostatics, the rest being attributed to electron exchange, induction, and dispersion terms. This rather low contribution of electrostatics is not supported by our results or by the data obtained by Zahn et al. from SAPT analysis.22 It has to be mentioned that electrostatic term definitions in SAPT and HVPT are identical. As shown in Table 1 and discussed in the Results section, the contribution of electron exchange, delocalization, and dispersion interactions is significantly lower than observed by Dong et al.1 The reason for this could be, for example, that they have used optimized structures in the minima of PES, while in our case snapshots from MD were used. The CA pair is closer bound in the minimum, hence shorter contact between the monomers and larger contribution from interaction components other than electrostatics. Therefore, one can argue that the real contribution of electrostatics might be underestimated when the global minimum is considered, while the contribution from delocalization and dispersion components could be overestimated. As a matter of fact, we have optimized a single CA complex at the DFT (B3LYP/aug-cc-pVDZ) and MP2 (MP2/Lanl2DZdp) levels and, subsequently, evaluated the interaction energy terms. The results are shown in Table 4. Both minima, DFT and MP2, are quite similar to each other, as can be seen in the graphics. When the Lanl2DZdp basis set is used to calculate the energy of the MP2 complex, the ∈(10) el term is 106% of the total (ΔEMP2) energy of interaction. The HF exchange, ∈HL ex constitutes 22% and ΔEdel constitutes 10% of

B3LYP/aug-cc-pVDZ

MP2/Lanl2DZdp

Energy

Lanl2DZ

Lanl2DZdp

Lanl2DZ

Lanl2DZdp

∈(10) el HL

−343. −281. −311. −313. 62.72 −30.60 −2.01 −13.41

−330. −278. −315. −325. 64.05 −37.20 −9.54 −26.25

−353. −282. −309. −316. 71.07 −26.86 −7.13 −18.12

−350. −277. −311. −329. 72.69 −33.86 −18.70 −36.11

ΔE ΔEHF ΔEMP2 ∈HL ex ΔEHF del ΔEcorr ∈(20) disp

“Geometry” indicates level used to optimize geometry, while “Energy” indicates basis set used to calculate energy terms. Lanl2DZ and Lanl2DZdp use ECP. The picture compares structure of the complex at the MP2 (green) and B3LYP (orange) level; hydrogen atoms are omitted in the graphics.

a

the total energy of interaction. On the other hand, when the snapshot geometries are considered (Table 4), the electrostatics HF is 96% of the total interaction, ∈HL ex is 6% and ΔEdel is 7%. Note that, in the equilibrium geometry, the electrostatic interaction is larger than the total interaction, while in the geometries from MD, it is smaller. This is in agreement with the results obtained by Zahn et al.22 for [MMIM][Cl] who also found that, at the minimum geometry, the electrostatic term is slightly larger than the total interaction energy and it drops below when the distance between monomers is increased. The amount of ∈HL ex and ΔEHF del are also similar. Kossmann et al.23 investigated up to nine-member linear complexes of [MMIM][Cl] pairs and found that the energy of interaction calculated per CA pair does not increase monotonically with the size of the complex. ΔE of a single pair is larger than half of ΔE in two-CA pair complex (four-body complex). It was also larger than 1/3 of ΔE in three-CA pair complex; however, for larger complexes the interaction increased above the value of a dimer. This could be compared with the data obtained here for single pair and two pairs in the 4-mer basis set (Table 5). In a four-body complex, there are six two-body interactions; since two monomers are cations and two are anions, one can easily expect to find four stabilizing and two destabilizing interactions. Table 5 follows the same approach taken in the whole paper  for a given central CA pair, a second cation and anion are added (not necessarily forming a second pair). For each snapshot, seven such complexes are selected, with varying distance of the second cation and anion from the central CA pair. It can be seen that the presence of the 2153

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

Table 5. Two- and Four-Body Total Interaction Energies [kJ mol−1] in the Four-Body Complexes, Calculated Using Lanl2DZ/ ECP Basis Seta Interaction/Complex C0:A0 A0:Cn C0:Cn Cn:An A0:An C0:An C0:A0:Cn:An a

P1

P2

P3

P4

P5

P6

P7

−284. 3.7 −221. 33.3 198. 8.0 −186. 30.2 174. 12.4 −266. 7.9 −564. 32.2

−284. 3.7 −203. 35.5 184. 8.6 −181. 32.8 190. 18.9 −259. 11.2 −537. 38.9

−284. 3.8 −209. 34.2 177. 7.2 −157. 26.8 171. 12.9 −240. 10.5 −527. 31.2

−284. 3.7 −145. 19.2 167. 10.7 −177. 30.7 187. 21.0 −223. 19.2 −471. 31.0

−284. 3.8 −161. 33.7 174. 8.9 −137. 31.8 162. 15.0 −203. 17.8 −444. 47.9

−284. 3.8 −171. 30.7 159. 12.7 −156. 35.0 157. 25.9 −161. 12.9 −450. 52.9

−284. 3.7 −176. 33.1 155. 7.6 −129. 21.9 144. 19.7 −150. 6.1 −433. 38.6

Every second row indicates error.

Therefore, in terms of the theory of intermolecular interactions it is the only “cooperative” effect that can be observed. The interaction energy components higher than electrostatics can be observed in the first coordination shell only. Beyond that shell, the interaction is perfectly described by the simple Coulomb formula applying formal charges of cation and anion (±1). The amount of the higer level interaction energy terms strongly depends on the distance. The presented analysis gives a slightly different picture compared to the results obtained by others.1,22 Contrary to those works, here, geometries from MD simulations are used instead of optimized complexes, resulting in a wider spectrum of distances between monomers and smaller contribution from electron exchange, delocalization, and dispersion interaction. One can postulate that this is closer to a “real” IL at room temperature. Based on the presented results and other works47−49 it should be possible to estimate the amount of charge reduction on the ion from quantum chemical calculations of smaller clusters, since the electron delocalization effects seem to be local. The components of the interaction energies could possibly correlated with properties like viscosity or melting point in a fashion proposed by Bernard.7

additional cation and anion has no influence whatsoever on the interaction in the central pair; it is constant and equals −284 kJ· mol−1. At the same time, the total interaction of the system drops down, because the distance increases. For the first 4-mer complex, the C0:A0 interaction is 50% of the total energy of interaction, while for the seventh complex, the C0:A0 interaction constitutes 66% of the total interaction energy. Kossmann et al.23 have found that the total interaction energy in their 4-mer complex was less than the 2-mer energy multiplied by two. This is most likely because, in their case, the monomers are more separated than in our four-body complex. However, when longer chains of dimers are considered, they find that the total interaction energy is larger than the 2-mer interaction multiplied by the number of monomers. In a linear structure, the attracting interactions are enhanced (each cation “sees” two anion next to it and vice versa), while the repulsive interactions are diminished. Dominant role of delocalization interaction in the three-body term is in agreement with pioneer study of Kolos et al. related to water−ion−water interactions.50 It may be tempting to extrapolate the calculated amounts of interaction energy components to other slightly different ILs, for example, from the [RMIM][PF6] family, where R stands for alkyl. However, caution must be exercised here. The properties of this class of IL have been found to change in a nonlinear fashion due to the nanostructuration, which becomes apparent for systems containing hexylmethylimidazolium and longer.51 At least two kinds of effects are expected to appear here: one is due to the heterogeneity of the system, the other is due to increasing dispersion interaction. Hoja et al. found that dispersion interaction due to the alkyl chain of varying length is present even in hydrogen bonded systems.52



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors would like to thank Dr. Gyorgy Hantal for providing the snapshots of the IL from MD simulations. The calculations have been carried out at Wroclaw Centre for Networking and Supercomputing (WCSS, Poland). This work has been financed by the National Science Centre (NCN, Poland) grant no. DEC-2011/01/B/ST5/06659.



CONCLUSIONS The analysis of interaction energy components in the series of snapshots from [BMIM][PF6] ionic liquid MD simulation confirms that the electrostatics is the dominating term; however, several other components should be considered as well. Calculations on a dimer system using basis set with polarized and diffuse function shows that electron exchange, delocalization, and dispersion each constitute a few percent of the total energy. In larger complexes, a contribution from threebody electron delocalization appears, being the only noticeable nonadditive component in agreement with earlier reports.



REFERENCES

(1) Dong, K.; Song, Y.; Liu, X.; Cheng, W.; Yao, X.; Zhang, S. Understanding Structures and Hydrogen Bonds of Ionic Liquids at the Electronic Level. J. Phys. Chem. B 2012, 116, 1007−1017. (2) Xantheas, S. S. On the Importance of the Fragment Relaxation Energy Terms in the Estimation of the Basis Set Superposition Error

2154

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

Correction to the Intermolecular Interaction Energy. J. Chem. Phys. 1996, 104, 8821−8824. (3) Earle, M. J.; Esperanca, J. M. S. S.; Gilea, M. A.; Lopes, J. N. C.; Rebelo, L. P. N.; Magee, J. W.; Seddon, K. R.; Widegren, J. A. The Distillation and Volatility of Ionic Liquids. Nature 2006, 439, 831− 824. (4) Chowdhury, S.; Mohan, R. S.; Scott, J. L. Reactivity of Ionic Liquids. Tetrahedron 2007, 63, 2363−2389. (5) Plechkova, N. V.; Seddon, K. R. Applications of Ionic Liquids in the Chemical Industry. Chem. Soc. Rev. 2008, 37, 123−150. (6) Rooney, D.; Jacquemin, J.; Gardas, R. L. In Ionic liquids; Kirchner, B., Ed.; Topics in Current Chemistry 290; Springer, 2009; pp 185− 212. (7) Bernard, U. L.; Izgorodina, E. I.; MacFarlane, D. R. New Insights into the Relationship between Ion-Pair Binding Energy and Thermodynamic and Transport Properties of Ionic Liquids. J. Phys. Chem. C 2010, 114, 20472−20478. (8) Fumino, K.; Peppel, T.; Geppert-Rybczyńska, M.; Zaitsau, D. H.; Lehmann, J. K.; Verevkin, S. P.; Köckerling, M.; Ludwig, R. The Influence of Hydrogen Bonding on the Physical Properties of Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 14064−14075. (9) Hayes, R.; Imberti, S.; Warr, G. G.; Atkin, R. The Nature of Hydrogen Bonding in Protic Ionic Liquids. Angew. Chem., Int. Ed. 2013, 52, 4623−4627. (10) Chałasiński, G.; Szcześniak, M. M. Origins of Structure and Energetics of van der Waals Clusters from ab initio Calculations. Chem. Rev. 1994, 94, 1723−1765. (11) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887−1930. (12) Kitaura, K.; Morokuma, K. New Energy Decomposition Scheme for Molecular-Interactions within Hartree-Fock Approximation. Int. J. Quantum Chem. 1976, 10, 325−340. (13) Sokalski, W. A. The Physical Nature of Catalytic Activity due to the Molecular Environment in Terms of Intermolecular Interaction Theory - Derivation of Simplified Models. J. Mol. Catal. 1985, 30, 395−410. (14) Sokalski, W. A.; Roszak, S.; Pecul, K. An Efficient Procedure for Decomposition of the SCF Interaction Energy into Components with Reduced basis set Dependence. Chem. Phys. Lett. 1988, 153, 153−159. (15) Dziekonski, P.; Sokalski, W. A.; Leszczynski, J. Physical Nature of Environmental Effects on Intermolecular Proton Transfer in (O2NOH ... NH3)(H2O)n and (ClH ... NH3)(H2O)n (n=1−3) Complexes. Chem. Phys. 2001, 272, 37−45. (16) Skwara, B.; Gora, R. W.; Bartkowiak, W. On the Influence of non-additive Interactions on the Optical Properties of the Selected Subsystems of Crystalline Urea. Chem. Phys. Lett. 2005, 406, 29−37. (17) Gora, R. W.; Sokalski, W. A.; Leszczynski, J.; Pett, V. B. The Nature of Interactions in the Ionic Crystal of 3-Pentenenitrile, 2-Nitro5-oxo, Ion(−1), Sodium. J. Phys. Chem. B 2005, 109, 2027−2033. (18) Dyguda, E.; Grembecka, J.; Sokalski, W. A.; Leszczynski, J. Origins of the Activity of PAL and LAP Enzyme Inhibitors: Toward Ab Initio Binding Affinity Prediction. J. Am. Chem. Soc. 2005, 127, 1658−1659. (19) Szefczyk, B.; Mulholland, A. J.; Ranaghan, K. E.; Sokalski, W. A. Differential Transition-State Stabilization in Enzyme Catalysis: Quantum Chemical Analysis of Interactions in the Chorismate Mutase Reaction and Prediction of the Optimal Catalytic Field. J. Am. Chem. Soc. 2004, 126, 16148−16159. (20) Szarek, P.; Dyguda-Kazimierowicz, E.; Tachibana, A.; Sokalski, W. A. Physical Nature of Intermolecular Interactions within cAMPDependent Protein Kinase Active Site: Differential Transition State Stabilization in Phosphoryl Transfer Reaction. J. Phys. Chem. B 2008, 112, 11819−11826. (21) Fernandes, A. M.; Rocha, M. A. A.; Freire, M. G.; Marrucho, I. M.; Coutinho, J. A. P.; Santos, L. M. N. B. F. Evaluation of Cation− Anion Interaction Strength in Ionic Liquids. J. Phys. Chem. B 2011, 115, 4033−4041.

(22) Zahn, S.; Uhlig, F.; Thar, J.; Spikermann, C.; Kirchner, B. Intermolecular Forces in an Ionic Liquid ([Mmim][Cl]) versus Those in a Typical Salt (NaCl). Angew. Chem., Int. Ed. 2008, 47, 3639−3641. (23) Kossmann, S.; Thar, J.; Kirchner, B.; Hunt, P. A.; Welton, T. Cooperativity in ionic liquids. J. Chem. Phys. 2006, 124, 174506. (24) Law, G.; Watson, P. R. Surface Tension Measurements of NAlkylimidazolium Ionic Liquids. Langmuir 2001, 17, 6138−6141. (25) Tokuda, H.; Hayamizu, K.; Ishii, K.; Abu Bin Hasan Susan, M.; Watanabe, M. Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 1. Variation of Anionic Species. J. Phys. Chem. B 2004, 108, 16593−16600. (26) Morrow, T. I.; Maginn, E. J. Molecular Dynamics Study of the Ionic Liquid 1-n-Butyl-3-methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2002, 106, 12807−12813. (27) Margulis, C. J.; Stern, H. A.; Berne, B. J. Computer Simulation of a “Green Chemistry” Room-Temperature Ionic Solvent. J. Phys. Chem. B 2002, 106, 12017−12021. (28) Gora, R. W.; Bartkowiak, W.; Roszak, S.; Leszczynski, J. A New Theoretical Insight into the Nature of Intermolecular Interactions in the Molecular Crystal of Urea. J. Chem. Phys. 2002, 117, 1031−1039. (29) Löwdin, P. O. Quantum Theory of Cohesive Properties of Solids. Adv. Phys. 1956, 5, 1−172. (30) Stone, A. J. Distributed Multipole Analysis, or How to Describe a Molecular Charge-distribution. Chem. Phys. Lett. 1981, 83, 233−239. (31) Sokalski, W. A.; Poirier, R. A. Cumulative Atomic Multipole Representation of the Molecular Charge Distribution and its Basis set Dependence. Chem. Phys. Lett. 1983, 98, 86−92. (32) Jeziorski, B.; van Hemert, M. C. Variation-Perturbation Treatment of Hydrogen-Bond between Water Molecules. Mol. Phys. 1976, 31, 713−729. (33) Boys, S. F.; Bernardi, F. Calculation of small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (34) Hantal, G.; Cordeiro, M. N. D. S.; Jorge, M. What Does an Ionic Liquid Surface Really Look Like? Unprecedented Details from Molecular Simulations. Phys. Chem. Chem. Phys. 2011, 13, 21230− 21232. (35) Bhargava, B. L.; Balasubramanian, S. Refined Potential Model for Atomistic Simulations of Ionic Liquid [bmim][PF6]. J. Chem. Phys. 2007, 127, 114510. (36) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; et al. General Atomic and Molecular Electronic-structure System. J. Comput. Chem. 1993, 14, 1347−1363. (37) Wadt, W. R.; Hay, P. J. Ab-Initio Effective Core Potentials for Molecular Calculations - Potentials for Main Group Elements Na to Bi. J. Chem. Phys. 1985, 82, 284−298. (38) Check, C. E.; Faust, T. O.; Bailey, J. M.; Wright, B. J.; Gilbert, T. M.; Sunderlin, L. S. Addition of Polarization and Diffuse Functions to the LANL2DZ Basis Set for P-Block Elements. J. Phys. Chem. A 2001, 105, 8111−8116. (39) Youngs, T. G. A.; Hardacre, C. Application of Static Charge Transfer within an Ionic-Liquid Force Field and Its Effect on Structure and Dynamics. Chem. Phys. Chem. 2008, 9, 1548−1558. (40) Zhong, X.; Liu, Z.; Cao, D. Improved Classical United-Atom Force Field for Imidazolium-Based Ionic Liquids: Tetrafluoroborate, Hexafluorophosphate, Methylsulfate, Trifluoromethylsulfonate, Acetate, Trifluoroacetate, and Bis(trifluoromethylsulfonyl)amide. J. Phys. Chem. B 2011, 115, 10027−10040. (41) Köddermann, T.; Paschek, D.; Ludwig, R. Molecular Dynamics Simulations of Ionic Liquids: A Reliable Description of Structure, Thermodynamics and Dynamics. Chem. Phys. Chem. 2007, 8, 2464− 2470. (42) 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. (43) Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave functions. J. Chem. Phys. 1955, 23, 1833−1840. 2155

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156

The Journal of Physical Chemistry B

Article

(44) Bader, R. F. W. A Quantum Theory of Molecular Structure and its Applications. Chem. Rev. 1991, 91, 893−928. (45) Szefczyk, B.; Sokalski, W. A.; Leszczynski, J. Optimal Methods for Calculation of the Amount of Intermolecular Electron Transfer. J. Chem. Phys. 2002, 117, 6952−6958. (46) Khaliullin, R. Z.; Bell, A. T.; Head-Gordon, M. Analysis of Charge Transfer Effects in Molecular Complexes Based on Absolutely Localized Molecular Orbitals. J. Chem. Phys. 2008, 128, 184112. (47) Schmidt, J.; Krekeler, C.; Dommert, F.; Zhao, Y.; Berger, R.; Delle Site, L.; Holm, C. Ionic Charge Reduction and Atomic Partial Charges from First-Principles Calculations of 1,3-Dimethylimidazolium Chloride. J. Phys. Chem. B 2010, 114, 6150−6155. (48) Prado, C. E. R.; Del Pópolo, M. G.; Youngs, T. G. A.; Kohanoff, J.; Lyndel-Bell, R. M. Molecular Electrostatic Properties of Ions in an Ionic Liquid. Mol. Phys. 2006, 104, 2477−2483. (49) Wendler, K.; Zahn, S.; Dommert, F.; Berger, R.; Holm, C.; Kirchner, B.; Delle Site, L. Locality and Fluctuations: Trends in Imidazolium-Based Ionic Liquids and Beyond. J. Chem. Theory Comput. 2011, 7, 3040−3044. (50) Clementi, E.; Kistenmacher, H.; Kolos, W.; Romano, S. Nonadditivity in Water-Ion-Water Interactions. Theor. Chim. Acta 1980, 55, 257−266. (51) Rocha, M. A. A.; Ribeiro, F. M. S.; Ferreira, A. I. M. C. L.; Coutinho, J. A. P.; Santos, L. M. N. B. F. Thermophysical Properties of [CN−1C1im][PF6] Ionic Liquids. J. Mol. Liq. 2013, 188, 196−202. (52) Hoja, J.; Sax, A. F.; Szalewicz, K. Chem.Eur. J., DOI: 10.1002/ chem.201303528.

2156

dx.doi.org/10.1021/jp411363d | J. Phys. Chem. B 2014, 118, 2147−2156