Mobility and association of ions in aqueous solutions: the case of imidazolium based ionic liquids †

The mobility and the mechanism of ion pairing of 1,1 electrolytes in aqueous solutions were investigated systematically on nine imidazolium based ionic liquids (ILs) from 1-methylimidazolium chloride, [MIM][Cl], to 1-dodecyl-3-methylimidazolium chloride, [1,3-DoMIM][Cl], with two isomers 1,2-dimethylimidazolium chloride, [1,2-MMIM][Cl], and 1,3-dimethylimidazolium chloride, [1,3-MMIM][Cl]. Molecular dynamics (MD) simulations, statistical mechanics calculations in the framework of the integral equation theory using one-dimensional (1D-) and three-dimensional (3D) reference interaction site model (RISM) approaches as well as conductivity measurements were applied. From experiment and MD simulations it was found that the mobility/diﬀusion coeﬃcients of cations in the limit of infinite dilution decrease with an increasing length of the cation alkyl chain, but not linearly. The aggregation tendency of cations with long alkyl chains at higher IL concentrations impedes their diﬀusivity. Binding free energies of imidazolium cations with the chloride anion estimated by RISM calculations, MD simulations and experiments reveal that the association of investigated ILs as model 1,1 electrolytes in water solutions is weak but evidently dependent on the molecular structure (alkyl chain length), which also strongly aﬀects the mobility of cations.


Introduction
Salt solutions are ubiquitous, and solvated ions greatly influence many processes, such as protein folding and conformational changes of nucleic acids, permeability, conductance, and electrostatic potential of cell membranes, micellization of surfactants and (in water) hydrophobic effects (also Hofmeister effects), and behaviour of enzymes, or have some other properties (e.g.acid-base properties) that influence the activity of biomolecules.Charged species also affect the mechanism and kinetics of chemical reactions consequently the mechanism of ion-exchange, the sol-gel transition, and many other phenomena.Thus, solvated ions were, not surprisingly, the subject of many experimental and theoretical studies.Recently, a systematic investigation of the interactions between oppositely charged ions in aqueous solutions of simple alkali halide salts has been carried out, based on conductance experiments and Monte Carlo simulations. 1 It has been found that ion pairs with a large difference in hydration free energy of cations and anions (large difference in size) have a relatively low association constant compared to those which exhibit similar hydration free energies and tend to associate in water.The microscopic structure was obtained by Monte Carlo simulations, using a simple two-dimensional water model.It has been shown that in all cases, the relative water density between the two ions increases, compared to the relative water density around two separate ions.In the case of two small ions, the two are associated due to stronger electrostatic interactions (in agreement with Collins' ''law of matching water affinities''), while in the case of two large ions, we attribute ion pairing to entropic reasons.On the other hand, in the case of two ions with very different water affinities, the small ion would be strongly hydrated and not release its water upon pairing, thereby leading to a small association constant.
In order to obtain more information about the mechanism of ion pairing in solutions, ionic liquids (ILs) can serve as an excellent model system, because they exist in diverse structures.For a long time, water was considered as an impurity in ILs.However, recent studies treat the mixtures of ILs and water as a novel class of solvents.These solvents are very attractive systems for different applications, in particular, they are used as the solvents for biopolymers [2][3][4] and therefore the body of literature examining IL-water mixtures has been growing steadily.
Molecular dynamic (MD) simulations have been used in combination with empirical force fields to study various ILs either in the bulk phase or in mixtures with water with varying IL concentrations. 5The association of ions in water has also been investigated using MD simulations directly.In ab initio MD simulations it was observed that small 1-ethyl-3-methylimidazolium ([1,3-EMIM] + ) cations and chloride anions remained associated in water throughout simulations, which was explained using hydrophobic effects. 6Ion association has also been observed in classical MD simulations, especially when hydrophobic ions are involved, as in the case of imidazolium-based cations with long alkyl chains paired with hydrophobic NTf 2 . 7When small hydrophilic iodide anions were used, however, ion association was observed with small dimethyl imidazolium cations but not with imidazolium cations that contained long octyl chains. 8Also the influence of the alkyl chain length of imidazolium-based cations on cation diffusion in water has been studied in a few cases.A linear decrease of diffusivity with increasing alkyl length was found. 9,10These values were obtained for large IL concentrations around 1.1-1.3M and 0.3-0.4M, respectively, where strong interactions among the ions determined selfdiffusion.However, the experimental values of diffusion coefficients for comparison were missing.
Among various experimental techniques for investigation of the mobility and the association properties of ionic systems in solutions, the electric conductivity measurement of dilute solutions is probably one of the most accurate techniques. 11,12ome attempts have been made to apply this method also on ILs diluted in aqueous solutions. 13,14Whereas for 1-butyl-3methylimidazolium chloride, [1,3-BMIM][Cl], at 298.15 K the value of K A is reported, 13 it is claimed that 1-ethyl-3-methylimidazolium tetrafluoroborate and 1-butyl-3-methylimidazolium tetrafluoroborate exist as free ions in aqueous solutions. 14fortunately there are no systematic investigations on diluted ILs in aqueous solutions which can provide more insights into their transport properties.
So we decided to carry out a systematic computational and experimental investigation of the influence of the ion size and its molecular structure on ion association and on ion mobility in water.Here we present the results for a series of imidazolium based chlorides with different alkyl chain lengths, from 1-methylimidazolium chloride, In this work, the experimental conductivity data will be analyzed by Barthel's low-concentration chemical model (lcCM) to obtain the association constants, K A (T), and limiting molar conductivities at infinite dilution, L N (T). 15L N (T) will be split into the ionic contributions with the help of known values of limiting conductivity of chloride ions, l N (T, Cl À ), enabling us to determine the diffusion coefficients D N (T, Cl À ).These diffusion coefficients will then be compared with values derived from MD simulations.In these simulations, the diffusivity of imidazolium cations will be derived in aqueous solution and close to the limit of infinite ion dilution to enable a straightforward comparison with D N .
To explain the obtained results of K A on a microscopic scale, we will apply statistical mechanics calculations in the framework of the integral equation theory using one-dimensional (1D-) and three-dimensional (3D-) reference interaction site model (RISM) approaches.The RISM integral equation theory is actively and successfully employed for the study of different liquids [16][17][18][19][20][21] (and references therein) including ILs. [22][23][24] Using the aforementioned approaches, one can obtain detailed information This journal is © the Owner Societies 2016 about atom-atom and molecular-atom interactions, including interactions among ions.Due to the latter it is thus possible to investigate ion association in the systems of interest.6][27][28][29][30][31][32][33][34][35] The PMF is a suitable property to study possible ion-binding, because the probability and the stability of ion pairs can be derived from it by estimating the depth of the corresponding PMF minima.

Molecular dynamics simulations
Imidazolium cations and chloride anions were described with the force field developed by Lopez et al. [36][37][38] However, partial charge on the atoms was adjusted to an aqueous environment.These charges were derived as described in Section 2.2 and are identical with those shown in Fig. B1-B9 in part B of the ESI.† Water was represented using the TIP5P model. 39D simulations were performed in the isothermal-isobaric ensemble (NPT) using GROMACS 4.5.5. 40Hydrogen was treated explicitly, while the lengths of covalent bonds that involved hydrogen were constrained using the LINC algorithm. 41Full period boundary conditions were used.The integration step size was 2 fs.The solution was simulated at 298 K using a velocity rescaling thermostat 42 and a coupling constant of 0.1 ps.The target pressure of 1 atm was regulated using a Berendsen barostat 43 with a coupling constant of 1 ps and a compressibility of 4.5 Â 10 À5 bar À1 .For long range electrostatic interactions beyond 1.5 nm the fast particle-mesh Ewald method with dispersion corrections was used. 44,45Lennard-Jones interactions were taken into account up to atom pair distances of 1.2 nm.The translation of the center-of-mass of the system was removed.
Six different systems were simulated that contained one cation respectively.These cations were solvated with 2050 water molecules, which correspond to an ion concentration of 0.0266 M. One chloride anion was added to ensure charge neutrality.After a minimization of the potential energy to remove initial close contacts between atoms, the systems were equilibrated for 1 ns.Subsequently, the systems were simulated for 50 ns for data analysis.The diffusion coefficients were derived from the resulting trajectories using the Einstein relation that connects the mean square displacement of a cation center-of-mass position, r com , with its diffusion coefficient, D: Mean square ion displacements were determined as a function of time intervals t, from which the diffusion coefficients were determined as the slope of the linear parts of these curves.The statistical errors were estimated as the difference between the diffusion coefficients derived from the first and second half of the curves.
In Fig. A1 in part A of the ESI † the mean square displacements of the center-of-mass of cations as a function of time interval, derived from MD simulations, are presented.Diffusion coefficients for the different cations were derived as slopes of these curves.
To examine the effect of cation aggregation on the diffusion coefficients, we also simulated ILs in aqueous solutions at high concentrations.In two MD simulations of 20 ns length, we solvated 90 ion pairs of [MIM][Cl] and [1,3-DoMIM][Cl] with 3500 water molecules, which resulted in IL concentrations of 1.0 M and 1.2 M, respectively.These very high concentrations are comparable to those that have been used before in ref. 9. Diffusion coefficients for the cations were derived as slopes of the curves presented in Fig. A2 in part A of the ESI.† Binding free energies, DG bind , of ion pairs were calculated using the thermodynamic integration method as implemented in GROMACS.In this process, DG bind is calculated by deriving the solvation free energy of Cl À at the cation binding site, while the ion pair is solvated in water, and the solvation free energy of Cl À in pure water without cations.The difference of these two solvation free energies is DG bind of the ion pair.This energy is also equivalent to the free energy change that occurred when Cl À is gradually displaced from its cation binding site into bulk water.The solvation free energy of Cl À is calculated in MD simulations by a step-wise switching off of all Coulomb interactions of the anion with its environment in the first phase.In the second phase, all Lennard-Jones interactions are gradually switched off, thereby fully decoupling the anion from the rest of the system.The coupling strength of Cl À to its environment is described by the coupling parameter l as shown in eqn (2): The parameter l is increased in steps from 0 to 1, so that the total energy, E(l), gradually changes from E 1 , the energy that contains all interactions of Cl À with its environment, to E 2 , where Cl À is fully decoupled.In our case l was increased with an increment of Dl = 0.1 from 0 to 1 for Coulomb and Lennard-Jones interactions, respectively.For each l-value, a MD simulation of 2 ns length was carried out (resulting in a total simulation length of 44 ns per system), in which after each 100 fs the derivative of the system enthalpy with respect to l, qH/ql, was written out.The results from the first 0.5 ns were discarded to account for equilibration.In the next step, values of qH/ql were averaged over each respective trajectory.These average hqH/qli values are plotted as a function of l in Fig. A3 and A4 in part A in the ESI.† A numerical integration of these curves over l yielded DG bind .In Lennard-Jones decoupling simulations where l is close to zero, i.e.where Cl À is almost fully decoupled from the system, it is necessary to use a soft-core potential.This additional potential avoids close contacts between atoms that would lead to singularities when the energy potentials are evaluated.For the soft-core potential we used an a-parameter of 0.5, a distance parameter of 0.3 nm and a l exponential of 1.The functional form of the soft-core potential and further details of the simulation procedure are given in the literature. 46e chose for Cl À binding sites at the cation that were found to be energetically most favorable according to results from 3D-RISM calculations.These 3D-RISM results were also in line with previously published work as discussed in the results part.To avoid diffusion of Cl À away from the cation binding site during free energy simulations, it is necessary to restrain the distance between Cl À and the cation hydrogen that according to 3D-RISM coordinates directly with Cl À .This distance restraint is particularly important taking into account that ion pair association was found to be very weak, thereby facilitating dissociation.The H-Cl distance was restrained with a harmonic potential with its minimum at the optimal distance and a force constant of 5000 kJ nm À2 mol À1 .The optimal H-Cl distances were found by minimizing the potential energy of the system.Otherwise, the same simulation parameters were used as described in the previous paragraph.Statistical errors were calculated using the block averaging method for each value of hqH/qli and error propagation was used to derive the error of DG bind , subsequently.

RISM calculations
The foundations of the 1D-and 3D-RISM theories were described in detail elsewhere 16,20,47,48 including our recent publications, [29][30][31][49][50][51] therefore it may suffice here to give only a brief outline of the aspects relevant to our study.
1D-RISM.The 1D-RISM approach is based on calculations of statistically averaged site-site (atom-atom) radial distribution functions (RDFs), g ab (r), via the numerical solution of the sitesite (atom-atom) Ornstein-Zernike integral equation 47 coupled with a corresponding 1D-closure relation.The RDFs describe the probability density of the distribution of sites (atoms), a and b, belonging to different particles (molecules, ions) of the system.The PMF between two ions in solution, W ab (r), is readily obtained as a logarithm of their RDF from the above mentioned coupled RISM equations for a solute-solvent mixture 25,52 as with k B being Boltzmann's constant and T the Kelvin temperature.In our case the ion-ion g ab (r) is the cation atom (a)-anion (b) RDF.The positions and values of the PMF at its first and second minima are corresponding, respectively, to a contact ion pair (CIP) and a solvent-separated ion pair (SSIP) (see Fig. 5 as an example).3D-RISM.This approach yields the three-dimensional (3D-) distribution of atoms (sites) of the solvent particles around a solute of arbitrary shape using molecular-atom (site) spatial distribution functions (SDFs).The ion binding of our interest can be described using the PMFs between two solutes (u) calculated from the corresponding SDFs g uu (r) as 28 W uu (r) = À k B T ln g uu (r). (4) The calculation of the SDFs is based on the solution of the 3D-RISM Ornstein-Zernike integral equation 28 coupled with a 3D-type closure relation.An important point is that the accuracy of both the 1Dand 3D-RISM PMFs depends significantly on the closure approximation used.It was established 53 that the most reasonable results are obtained using the Kovalenko-Hirata closure. 54orrespondingly, we used 1D- 55 and 3D- 28 versions of this closure in our calculations.As in this study, first of all, we were interested in the possibility of contact ion pairing (see below), and therefore we limited our discussion to the behavior of the PMF at its first minimum.
In the framework of the RISM calculations, the systems under study were considered as a mixture of water molecules and ILs with the common form of the interaction potential represented by the long-range electrostatic (Coulomb) and short-range Lennard-Jones (LJ) terms.The corresponding LJ parameters were taken from the Optimized Potentials for Liquid Simulations (OPLS) force field. 56For water the modified version of the SPC/E model (MSPC/E) was used. 57The parameters of the optimized (bond lengths and angles) structures of the IL cations under study were adopted by the DFT (Density Functional Theory) method at the B3LYP/6-311++G** level in combination with a polarizable continuum model (PCM) for the water solvent [58][59][60][61] as implemented in the ab initio code GAUSSIAN 09. 62The optimized geometries are presented in the ESI † in the form of atom coordinates (Tables in part B).The partial charges of IL atoms in water were derived from the electrostatic potential induced by the structurally optimized cations (ESP charges) using the CHELPG scheme. 63Charges on chemically equivalent atoms were averaged.The charges are shown in Fig. B1-B9 in part B of the ESI.† The RISM calculations were carried out using the rism1d and rism3d.snglpntcodes from the AmberTools package (version 14). 64The numerical solution of the 1D-and 3D-RISM integral equations was obtained using the MDIIS (Modified Direct Inversion in the Iterative Subspace) iterative scheme. 65The 1D-RISM equations were solved on a 1D grid of 16  .The residual tolerance of 10 À6 was also chosen for both 1D-and 3D-RISM calculations.All above mentioned parameters are large enough to accommodate the complex together with sufficient solvation space around it so that the obtained results are without significant numerical errors.All calculations were performed for ambient conditions and the same concentrations of ILs(aq) as in the conductivity experiments.

Conductivity measurements and data analysis
Materials.All compounds were purchased from Ionic Liquids Technologies (io-li-tec, Germany, the details are given in ESI † Part C) and were used without further purification.They were dried for 24 h at T E 313 K with a vacuum line ( p o5 Pa) and stored in a desiccator over P 2 O 5 before use.
Stock solutions were prepared by mass from the pure compounds and demineralized distilled water.Demineralized water was distilled two times in a quartz bidestillation apparatus (Destamat Bi 18E, Heraeus).The final product with the specific conductivity o6 Â 10 À7 S cm À1 was distilled into a flask permitting storage under an atmosphere of nitrogen permitting storage and transfer of water into the measuring cell under an atmosphere of nitrogen.The concentration of the stock solution was checked by titration with a standard solution of AgNO 3 potentiometric indication.
Conductivity measurements.The conductivities of the solutions were determined with the help of a three-electrode measuring cell, described elsewhere. 66The cell was calibrated with dilute potassium chloride solutions 67 and immersed in the high precision thermostat described previously. 68The temperature dependence of the cell constant was taken into account. 67The oil bath can be set to each temperature using a temperature program with a reproducibility of 0.005 K.The temperature in the precision thermostat bath was additionally checked using a calibrated Pt100 resistance thermometer (MPMI 1004/300 Merz) in connection with a Multimeter HP 3458A.The resistance measurements of the solutions in the cell were performed using a precision LCR Meter Agilent 4284A.The measuring procedure, including corrections and extrapolation of the sample conductivity to infinite frequency, has been described previously. 68he molar concentrations c were determined from the masses and the corresponding solution densities, d, of stock solutions and the last solutions in the cell.They were determined at 298.15 K using the method of Kratky et al. 69 C2 in part C of the ESI † together with molar conductivities, L(c), of all investigated systems.Molar conductivities are given as a function of electrolyte molality, m, which is related to the corresponding (temperaturedependent) molar concentration, c, via c = mÁd/(1 + M 2 Ám), where M 2 is the molar mass of the solute (given also in Table C2 in part C of the ESI †) and d is the density of solution.The estimated uncertainty of d is within 0.005 kg m À3 .The results for L(c) are shown in Fig. 2a for [1,3-EMIM][Cl] in water, whereas Fig. 2b compares the data for all investigated solutes at 298.15 K. Considering the sources of error (calibration, measurements, impurities), the values for L(c) are believed to be certain within 0.3%.
Data analysis.The obtained molar conductivities, L(c) (Table C2 in part C of the ESI † and Fig. 1 and 2), of dilute (c E r0.005 M) IL solutions were analyzed in the framework of Barthel's low-concentration chemical model (lcCM) 15 which describes successfully the thermodynamic and transport properties of diluted solutions.For the evaluation of molar conductivities, L, the lcCM uses the equations (5) where L N is the molar conductivity of the solute at infinite dilution, (1 À a) is the fraction of oppositely charged ions acting as ion pairs, and K A is the standard-state (infinite dilution) ion association constant.y AE 0 is the corresponding mean activity coefficient of the free ions, (y AE 0 ) 2 = y 0 + y À 0 , k D is the Debye parameter, e o is the proton charge, e is the relative permittivity of the solvent, and e o is the permittivity of a vacuum.Other symbols have their usual meaning.The coefficients of eqn (5) are explained in ref. 15.The parameters S and E are estimated with the help of the solvent parameters.The coefficients J 1 and J 2 additionally depend on the distance parameter R, which represents the distance up to which oppositely charged ions are counted as non-conducting ion pairs, and hence is the upper limit of ion-pair association represented by the upper limit of the association constant K A and the distance parameter of the mean activity coefficient y AE 0 (eqn ( 6) and ( 7)).
Analysis of the conductivity data was carried out by setting the coefficients S, E, J 1 and J 2 of eqn (5) to their calculated values 15 using the properties of water, d(T), 70 e(T), 71 Z(T) 72 (Table C1 in part C of the ESI †) and the distance parameter R. The lower limit a of the association integral is the distance of closest approach of the cation and the anion (contact distance), a = a + + a.It was calculated from the ionic radii of the chloride anion (a À = 0.181 nm) 15 and the same effective cation radius, a + = 0.325 nm for all investigated ILs.For the latter it was assumed that the anion mainly interacts with the acidic hydrogen at the C2 position of the imidazolium ring (for 1,3-ILs).Thus, the chosen a + corresponds to the molecular dimension of the imidazolium ring along the C2-H axis obtained from semiempirical MOPAC calculations and used also in our previous work on [1,3-BMIM][Cl] in water. 13rom extensive investigations on electrolyte solutions in different solvents it has been found 15 that the upper limit of association is most realistically defined as R = a + + a _ + nÁs, where s is the length of an orientated solvent molecule (located between the cation and the anion) and n is an integer.A value of s = d OH = 0.28 nm for water was taken into account.Assuming that only contact (CIP) and solvent-shared (SIP) ion pairs were likely to form, a value of n = 1 was used and R(J 2 ) was fixed to 0.786 nm.
Then two-parameter fits were used to obtain the limiting values of molar conductivity, L N , and the association constant, K A , by non-linear least squares iterations.All results are summarized in Table C3 in part C of the ESI.† Here it has to be stressed that the experimental quantification of weak association is not a trivial task.According to Marcus and Hefter 11 even the conductivity measurements, as one of the most established techniques, yield only a rough estimation for the association constant, K A r 10, what is the case in our work.However, in our previous work we have shown that this method is applicable at least for an excellent comparison of the associating behavior in diverse weak associating systems. 1,73,74rom Table C3 in the ESI † (part C) it can be seen that the same trend in K A was found at all temperatures.

Ion mobility
Using a literature value for l N À (Cl À ) 75 (Table C1 in part C of the ESI, † l N (Cl À , 298.15 K) = 76.35 S cm 2 mol À1 ), the limiting molar conductivity can be split into separate ionic contributions, l N i , where limiting ionic conductivities were estimated for all cations.Consequently, the diffusion coefficients D N i were calculated by using the relation where symbols have the same meaning as explained in eqn ( 6) and (7).Diffusion coefficients at infinite dilution are the characteristic properties of ionic transport, which are not affected by ionic interactions.Thus, for Cl À at 298.15 K, a value of 20.35 Â 10 À10 m 2 s À1 can be estimated.For all investigated cations, values at 298.15 K are gathered in Table 1 (together with values of L N and l N + at 298.15 K for all investigated solutes) and presented in Fig. 3, as a function of the number of C atoms in the alkyl side chain, together with values that were obtained from MD simulations that were carried out at the same temperature of 298 K.
As expected, the mobility of cations decreases with the increasing length of the alkyl chain (Table 1 and Fig. 3), but not linearly with its number of C atoms as it has been obtained in previous MD simulations of imidazolium-based cations in aqueous solution as reported by Bhargava and Klein. 9In these previous simulations high ion concentrations above 1 M were used in contrast to the very low ion concentration (c = 0.0266 M) and the lack of inter-ionic interactions in our case.For the medium sized cations [1,3-BMIM] + and [1,3-HMIM] + we observe a good agreement between calculated values in the present work and values reported in ref. 9 (about 9 Â 10 À10 m 2 s À1 and 7 Â 10 À10 m 2 s À1 , respectively).In the case of small [1,3-EMIM] + cations the previously reported value was substantially larger (15 Â 10 À10 m 2 s À1 ) and for larger [1,3-OMIM] + a slower diffusion (2 Â 10 À10 m 2 s À1 ) was predicted.These deviations are expected due to the presence of strong inter-ionic interactions that are found at high ion concentrations.For instance, it is well known that cations with long alkyl chains tend to aggregate to form structures similar to micelles. 76Even more, Goodchild and coworkers have found out that in water solutions of 1-alkyl-3methylimidazolium bromides there are no immediately apparent aggregates in the systems with n c = 2 and n c = 4, 77 but already at [1,3-HMIM][Br] small oblate aggregate forms were detected by SANS, increasing in size with increasing concentration.These bulky aggregates effectively lower the diffusivity of the involved cations.
This assumption was confirmed by the results of MD simulations at high concentrations.We observed that within less than 1 ns [1,3-DoMIM] + cations formed micelles that remained stable throughout the rest of the simulations as expected.In contrast, [MIM] + cations remained dispersed in water and did not show any signs of aggregation.The cation distribution in water is visualized for these two cases in Fig. 4. Diffusion coefficients of 10.4 and 1.1 Â 10 À10 m 2 s À1 were derived for [MIM] + and [1,3-DoMIM] + , which can be compared to 12.8 and 6.3 Â 10 À10 m 2 s À1 , respectively, close to the limit of infinite IL dilution.The drastic reduction of [1,3-DoMIM] + diffusivity is a result of the cation aggregation.[1,3-DoMIM] + cations diffuse as members of these micelle aggregates, which move much slower through the solvent than single cations because of the large aggregate sizes.

Ion association
As can be seen in Tables 2a and C3 (part C of the ESI †), K A values are small (B2.5 r K A r B6) but distinctly higher than those obtained using the same model recently used for alkali metal halides in water. 1 Before considering the RISM results to explain the experimental results for K A , we should give some information on features of ion-pairing in imidazolium-based ILs.The positive charge of their cations is located at the imidazolium ring 78 (see also Fig. B1-B9 in part B of the ESI †).Therefore, it is not a surprise that the interactions of anions with the non-polar cation alkyl chains are very weak. 95][86] Obviously, this location of the anion relative to the cation is a result of the strong electrostatic interaction mediated by the ring.5][86] 1).Thus, ion pairing is expected for species that can bind through specific intermolecular interactions such as hydrogen   bonding.Based on this information, we give the results of 1D-RISM PMF calculations for the above-named atoms (Table 2b and c).Fig. 5 presents the 1D-RISM PMF for H i -Cl as an example.
As it can be seen from the 1D-RISM PMF values (Tables 2b  and c), the hydrogen-bonded interactions between the anion and cations mediated by ring hydrogen atoms are weak.This result is in agreement with NMR/IR/Raman spectroscopy data 80,87 and MD simulations. 9,82,85It is also noteworthy that the obtained H i -Cl PMF values are rather similar.This fact suggests that the probabilities to form cation-anion pairs with the different available ring hydrogens are similar.However, earlier studies that used different methods 85,86,88 have shown that the H-bonding affinity with the anion, including Cl À , is highest for the acidic hydrogen (the H atom bonded to the C atom located between the two N atoms of the ring), compared to the other two ring hydrogens, i.e. it could be that the C1-H1Á Á ÁCl bond for our systems, with the exception of [1,2-MMIM] + having a structure different from other ILs (Fig. 1).
The analysis of the results of 3D-RISM calculations (Table 2d) has indicated that chloride ions prefer to interact with cations via hydrogen H7 (N2-H7Á Á ÁCl) in [MIM][Cl] aqueous solutions, via hydrogen H9 (N2-H9Á Á ÁCl) in [1,2-MMIM][Cl] aqueous solutions, and via atom H1 (C1-H1Á Á ÁCl) in all other IL aqueous solutions.The energetically most favourable position of the anion relative to the cations is shown in Fig. 6a-i.These results are also consistent with the charge distributions that were used for our IL models (see the Figures in part B of the ESI †).Thus, the valid values of 1D-RISM PMFs are for the above-named atoms of cations (Table 2b and c in bold) and should be chosen for the discussion.
As to the deeper minima of 1D-RISM PMFs (H2/H3-Cl À ) in comparison with (H1-Cl À ) for ILs ranging from ethyl to dodecyl (Table 2c), that, probably, they are overvalued due to the influence of additional charges of the alkyl chain atoms at the transition from [1,3-MMIM][Cl] to [1,3-EMIM][Cl] and to the next ILs in the set (Fig. B3-B9 in part B of the ESI †).
Thus, we have found the same type of cation-anion hydrogenbonded interaction (Cl-H1Á Á ÁAnion) for ILs ranging from dimethyl to dodecyl as it was shown in the literature previously. 85,86,88oreover, for the systems under study, the additional confirmation of our results follows from MD simulations of [1,3-EMIM][Cl](aq) 6,89 and [1,3-OMIM][Cl](aq). 89The data of NIR spectroscopy demonstrate a similar pairing of chloride anions with 1-alkyl-3-methylimidazolium cations ([1,3-AMIM] + ). 86ne can also see from Table 2 that the increase of alkyl chain length affects the interaction between hydrogen atoms and the chloride anion only slightly.
According to the results from the 3D-RISM, we found for [MIM] + and [1,2-MMIM] + a different energetically most favorable anion position relative to the cation compared to the other IL solutions.Chloride appears to interact most favorably in these systems with H7/H9 ([MIM] + /[1,2-MMIM] + ).This is not surprising, considering the charge distributions shown in Fig. B1 and B2 of the ESI.† The largest positive charges can be found on H7/H9, which resulted in stronger attractive electrostatic interactions with the anion.
Table 2b and c reveal that the cation-anion interaction strength, indeed, depends on the particular structure of the ions contained in the ILs.Whereas the order in could be explained by the change in 1D-RISM PMFs this is not the case for the other investigated systems.
For four different imidazolium cations, namely for [MIM] + , [1,3-BMIM] + , [1,3-HMIM] + and [1,3-DoMIM] + , the binding free energies with Cl À , DG MD bind , were also obtained from MD simulations, using thermodynamic integration as specified in Section 2.1.From measured values of K A binding free energies were calculated using the relation DG exp bind = ÀRTÁln K A .The values of DG MD bind and DG exp bind are presented in Table 2e.In the case of MIM + , DG MD bind seems to be overestimated, whereas for the other three cases with alkyl chains longer than n 4 0, the results from MD simulations agree with experimental values and deviations are within the statistical uncertainty.The somewhat overestimated binding in the case of [MIM][Cl] could be a result of a neglected transfer of charge from Cl À to [MIM] + in the MD model.Such a charge transfer is to be expected considering the tighter ion pairing of Cl À to the small cation, compared to other larger cations with more delocalized charge.Compared to the RISM, binding free energies derived using MD show a larger variation with the alkyl chain length, similar to what we observed in the experiment.The basically constant binding free energy (PMF) found in the RISM for chain lengths longer than or equal to n = 2 is likely a result of the fixed geometry of the cation that does not consider the flexibility of the cation alkyl chain.Overall, MD simulations confirm that the strength of ion association varies indeed with the cation alkyl chain length and that these variations are nonlinear and more generally that imidazolium-Cl association is very weak and that it depends on the specific imidazolium cation structure.
In Fig. 7 DG exp bind and DG MD bind (without the value for [MIM][Cl]) are compared with the minimal values of 1D-RISM-PMFs, calculated between chloride and those cation hydrogen atoms that were selected as described above, W RISM bind .The trend of W RISM bind values from the 1D-RISM suggests that the interaction between the anion and the cation is slightly decreasing with increasing alkyl chain length (n c 4 2).This can be rationalized with the fact that the cation charge becomes more delocalized with increasing cation size, thereby somewhat weakening the interaction strength with the anion.Our result is also confirmed by recent MD simulations, 89 where [1,3-MMIM] + was found to associate more strongly with Cl À than [1,3-OMIM] + .
Here it has to be stressed that the model of ion association applied in the data analysis is assuming only equilibrium between (neutral) ion pairs and free ions (Section 2.3) in diluted solutions (lcCM).For investigated ILs this assumption seems to be valid for systems with a shorter alkyl chain, where the association behavior can also be confirmed well by 1D-and 3D-RISM calculations.However, the results from MD simulation are in reasonable agreement with those from the experiment.
As it was already mentioned in Section 3.1, at longer chains the additional interactions, resulting in bulky aggregates, became inevitable. 77Despite the fact that experiments were carried out in diluted solutions (o0.005 mol dm À3 ) they may affect the possible ion equilibrium in the solutions and in this case the applied lcCM model would not be valid.Moreover, the investigated ILs should be regarded as typical 1,1-electrolytes, and with increasing length of the side alkyl chain hydrophobic interactions also become important and may lead to (additional) cation-cation interactions which are not taken into account in the here applied model.
From MD simulation, the diffusion coefficients close to infinite dilution, D N +,calc , were calculated and compared with the experimental values D N +,exp .The measured values were obtained from the cationic contributions, l N + , to the molar conductivity at infinite dilution, L N , and were estimated by applying the low concentration Chemical Model (lcCM) to experimental conductivity data.Experimental and calculated values from MD simulations were found to be in very good agreement, thereby confirming the conductivity measurements.
As expected, it was found that the mobility of cations is decreasing with increasing length of its alkyl chain, but not linearly as it has been reported previously. 9That seemingly linear decrease could be explained by the aggregation of cations with long alkyl chains and thus slower diffusivity, when larger cation concentrations, i.e. stronger inter-ionic interactions, are present.
The results of RISM calculations have shown a weak ion pairing between the anions and cations.It was found that their interactions are hydrogen-bonded and mediated by ring hydrogen atoms.From minimal values of 1D-RISM PMFs for atom pairs N2-Cl À ([MIM][Cl] and [1,2-MMIM][Cl]) or C1-Cl À (all other ILs) it was concluded that the ion association weakens first with the increasing length of the alkyl chain up to two C atoms but that beyond this length this association should not be affected any further.This assumption was partially confirmed by the values of association constants, K A , yielded also from experimental conductivity data by lcCM.Generally, estimated K A values are slightly higher than those obtained for 1,1-electrolytes in water, but the values for [1,3-BMIM][Cl] show distinctly a strong deviation with a value of K A = B5.2 at 298.15 K being still in reasonable agreement with a reported value of 6.4. 13This fact led to the assumption that there must be a difference in cation-anion interactions depending on the particular structure of the cations in the ILs, which may also affect the ion pairing process.
The most favorable anion position at cations was predicted by the 3D-RISM calculations.In particular, for the case of 1,3-ILs the chloride ions prefer to interact with cations via the acidic hydrogen of the IL ring.However, 3D-RISM data do not predict any strong influence of the alkyl chain length on the ion-pairing for n c Z 4.
Overall, we demonstrated that studied ILs serve as an excellent model system for investigating the influence of the ion structure on the mobility and ion pairing by combining MD simulations and RISM calculations with experiments, even in aqueous solutions, where the interionic interactions are weak.

Fig. 2
Fig. 2 (a) Molar conductivities of [1,3-EMIM][Cl] in water from 278.15 to 313.15 K in steps of 5 K in the concentration range of B0.0002 o c/M o B0.005;(b) molar conductivities of investigated IL chlorides at 25 1C in water; symbols denote the experiment and lines denote the lcCM calculations.

Fig. 3
Fig. 3 Diffusion coefficients of IL cations in water at 298.15 K as a function of the number of carbon atoms in the side chain; ( ) experiment; ( ) MD simulations.

Fig. 4
Fig. 4 Representative structures of [MIM][Cl] in solution (left hand side) and [1,3-DoMIM][Cl] (right hand side) taken from equilibrated MD simulations at very high IL concentrations (1.0 M and 1.2 M, respectively).It can be clearly seen that [1,3-DoMIM][Cl] forms micelles in contrast to [MIM][Cl], which remained dispersed.The display of water and hydrogens was omitted for clarity.Color code: C = cyan, N = blue and Cl = green.

Table 1
Limiting values of molar conductivity, L N , and limiting values for cations, l N + , together with a comparison of calculated D N +,calc , with experimental D N +,exp values of the cation diffusion coefficients at infinite dilution of all investigated ILs a This journal is © the Owner Societies 2016

Table 2
(a) Association constants, K A , at 298.15 K, the minimal values of 1D-RISM PMFs for atom pairs (b) X-Cl À and (c) H i -Cl À (as described in the text), (d) the minimal values of 3D-RISM PMFs and (e) binding free energies of imidazolium cations with Cl À as obtained from experiments and MD simulations.The valid 1D-RISM values and interacting cation atoms are given in bold.The atom numbering is given in Fig.1 a