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

Abstract

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/diffusion coefficients 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 diffusivity. 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 affects the mobility of cations.

---

1. 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.¹

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–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.⁵ The 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.⁶ Ion 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₂.⁷ When 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.⁸ Also 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,10 These values were obtained for large IL concentrations around 1.1–1.3 M and 0.3–0.4 M, respectively, where strong interactions among the ions determined self-diffusion. 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,12 Some attempts have been made to apply this method also on ILs diluted in aqueous solutions.^13,14 Whereas for 1-butyl-3-methylimidazolium chloride, [1,3-BMIM][Cl], at 298.15 K the value of K_A is reported,¹³ it is claimed that 1-ethyl-3-methylimidazolium tetrafluoroborate and 1-butyl-3-methylimidazolium tetrafluoroborate exist as free ions in aqueous solutions.¹⁴ Unfortunately 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, [MIM][Cl] to 1-dodecyl-3-methylimidazolium chloride, [1,3-DoMIM][Cl]. Among them there are also two isomers 1,2-dimethylimidazolium chloride, [1,2-MMIM][Cl] and 1,3-dimethylimidazolium chloride [1,3-MMIM][Cl]. The structures of all investigated systems are presented in Fig. 1.

Fig. 1 Structures of investigated IL cations with corresponding atom numeration: (a) 1-methylimidazolium, [MIM]⁺; (b) 1,2-dimethylimidazolium, [1,2-MMIM]⁺; (c) 1,3-dimethylimidazolium, [1,3-MMIM]⁺; (d) 1-ethyl-3-methylimidazolium, [1,3-EMIM]⁺; (e) 1-butyl-3-methylimidazolium, [1,3-BMIM]⁺, x = 1, y = 8, and z = 13; (f) 1-hexyl-3-methylimidazolium, [1,3-HMIM]⁺, x = 3, y = 10, and z = 17; (g) 1-octyl-3-methylimidazolium, [1,3-OMIM]⁺, x = 5, y = 12, and z = 21; (h) 1-decyl-3-methylimidazolium, [1,3-DeMIM]⁺, x = 7, y = 14, and z = 25; (i) 1-dodecyl-3-methylimidazolium, [1,3-DoMIM]⁺, x = 9, y = 16, and z = 29.

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, Λ^∞(T).¹⁵Λ^∞(T) will be split into the ionic contributions with the help of known values of limiting conductivity of chloride ions, λ^∞(T, Cl⁻), enabling us to determine the diffusion coefficients D^∞(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^∞.

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–21 (and references therein) including ILs.^22–24 Using the aforementioned approaches, one can obtain detailed information 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. For our purpose the free energy of ion pairing will be estimated based on the RISM potential of mean force (PMF).^25–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.

2. Methods

2.1. Molecular dynamics simulations

Imidazolium cations and chloride anions were described with the force field developed by Lopez et al.^36–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.³⁹

MD simulations were performed in the isothermal–isobaric ensemble (NPT) using GROMACS 4.5.5.⁴⁰ Hydrogen was treated explicitly, while the lengths of covalent bonds that involved hydrogen were constrained using the LINC algorithm.⁴¹ Full period boundary conditions were used. The integration step size was 2 fs. The solution was simulated at 298 K using a velocity rescaling thermostat⁴² and a coupling constant of 0.1 ps. The target pressure of 1 atm was regulated using a Berendsen barostat⁴³ with a coupling constant of 1 ps and a compressibility of 4.5 × 10⁻⁵ bar⁻¹. For long range electrostatic interactions beyond 1.5 nm the fast particle-mesh Ewald method with dispersion corrections was used.^44,45 Lennard-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 [MIM]⁺, [1,3-EMIM]⁺, [1,3-BMIM]⁺, [1,3-HMIM]⁺, [1,3-OMIM]⁺, [1,3-DeMIM]⁺ or [1,3-DoMIM]⁺, 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, ΔG_bind, of ion pairs were calculated using the thermodynamic integration method as implemented in GROMACS. In this process, ΔG_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 ΔG_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 λ as shown in eqn (2):

E(λ) = (1 − λ)E₁ + λE₂ (2)

The parameter λ is increased in steps from 0 to 1, so that the total energy, E(λ), gradually changes from E₁, the energy that contains all interactions of Cl⁻ with its environment, to E₂, where Cl⁻ is fully decoupled. In our case λ was increased with an increment of Δλ = 0.1 from 0 to 1 for Coulomb and Lennard-Jones interactions, respectively. For each λ-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 λ, ∂H/∂λ, was written out. The results from the first 0.5 ns were discarded to account for equilibration. In the next step, values of ∂H/∂λ were averaged over each respective trajectory. These average 〈∂H/∂λ〉 values are plotted as a function of λ in Fig. A3 and A4 in part A in the ESI.† A numerical integration of these curves over λ yielded ΔG_bind. In Lennard-Jones decoupling simulations where λ 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 α-parameter of 0.5, a distance parameter of 0.3 nm and a λ exponential of 1. The functional form of the soft-core potential and further details of the simulation procedure are given in the literature.⁴⁶ We 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⁻² mol⁻¹. 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 〈∂H/∂λ〉 and error propagation was used to derive the error of ΔG_bind, subsequently.

2.2. 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–31,49–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_αβ(r), via the numerical solution of the site–site (atom–atom) Ornstein–Zernike integral equation⁴⁷ coupled with a corresponding 1D-closure relation. The RDFs describe the probability density of the distribution of sites (atoms), α and β, belonging to different particles (molecules, ions) of the system. The PMF between two ions in solution, W_αβ(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

W_αβ(r) = k_BT ln g_αβ(r) (3)

with k_B being Boltzmann's constant and T the Kelvin temperature. In our case the ion–ion g_αβ(r) is the cation atom (α)–anion (β) 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²⁸

W^uu(r) = − k_BT ln g^uu(r). (4)

The calculation of the SDFs is based on the solution of the 3D-RISM Ornstein–Zernike integral equation²⁸ coupled with a 3D-type closure relation.

An important point is that the accuracy of both the 1D- and 3D-RISM PMFs depends significantly on the closure approximation used. It was established⁵³ that the most reasonable results are obtained using the Kovalenko–Hirata closure.⁵⁴ Correspondingly, we used 1D-⁵⁵ and 3D-²⁸ 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.⁵⁶ For water the modified version of the SPC/E model (MSPC/E) was used.⁵⁷ The 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–61 as implemented in the ab initio code GAUSSIAN 09.⁶² The 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.⁶³ Charges 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.snglpnt codes from the AmberTools package (version 14).⁶⁴ The numerical solution of the 1D- and 3D-RISM integral equations was obtained using the MDIIS (Modified Direct Inversion in the Iterative Subspace) iterative scheme.⁶⁵ The 1D-RISM equations were solved on a 1D grid of 16 384 points with a spacing of 2.5 × 10⁻³ nm and 10 MDIIS vectors. A residual tolerance of 10⁻⁶ was chosen. The 3D-RISM equations were solved on a three-dimensional grid of 256 × 270 × 250 points for [MIM][Cl] and [1,3-MMIM][Cl], 270 × 270 × 250 points for [1,2-MMIM][Cl] and [1,3-EMIM][Cl], 270 × 280 × 250 points for [1,3-BMIM][Cl], 280 × 288 × 252 points for [1,3-HMIM][Cl], 288 × 288 × 250 points for [1,3-OMIM][Cl], 294 × 300 × 250 points for [1,3-DeMIM][Cl] and 300 × 320 × 250 points for [1,3-DoDMIM][Cl] with 4 MDIIS vectors and a spacing of 0.025 nm in a parallelepiped cell. These parameters correspond to a cell of size (6.38 × 6.75 × 6.25) nm³ for [MIM][Cl] and [1,3-MMIM][Cl], (6.75 × 6.75 × 6.25) nm³ for [1,2-MMIM][Cl] and [1,3-EMIM][Cl], (6.75 × 7.00 × 6.25) nm³ for [1,3-BMIM][Cl], (7.00 × 7.20 × 6.30) nm³ for [1,3-HMIM][Cl], (7.20 × 7.20 × 6.25) nm³ for [1,3-OMIM][Cl], (7.35 × 7.50 × 6.25) nm³ for [1,3-DeMIM][Cl], and (7.50 × 8.00 × 7.25) nm³ for [1,3-DoDMIM][Cl]. The residual tolerance of 10⁻⁶ 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.

2.3. 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 ≈ 313 K with a vacuum line (p <5 Pa) and stored in a desiccator over P₂O₅ 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 <6 × 10⁻⁷ S cm⁻¹ 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₃ potentiometric indication.

Conductivity measurements. The conductivities of the solutions were determined with the help of a three-electrode measuring cell, described elsewhere.⁶⁶ The cell was calibrated with dilute potassium chloride solutions⁶⁷ and immersed in the high precision thermostat described previously.⁶⁸ The temperature dependence of the cell constant was taken into account.⁶⁷ The 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.⁶⁸

The 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.⁶⁹ using a Paar densimeter (DMA 5000) combined with a precision thermostat. A linear change of d with increasing salt content for diluted solutions was assumed, i.e. d(T) = d₀(T) + bm, where d₀(T) is the density of water (Table C1 in part C of the ESI†), m is the molality, and the b-coefficients are assumed as independent of temperature. They are given in Table C2 in part C of the ESI† together with molar conductivities, Λ(c), of all investigated systems. Molar conductivities are given as a function of electrolyte molality, m, which is related to the corresponding (temperature-dependent) molar concentration, c, via c = m·d/(1 + M₂·m), where M₂ 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⁻³. The results for Λ(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 Λ(c) are believed to be certain within 0.3%.

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 ∼0.0002 < c/M < ∼0.005; (b) molar conductivities of investigated IL chlorides at 25 °C in water; symbols denote the experiment and lines denote the lcCM calculations.

Data analysis. The obtained molar conductivities, Λ(c) (Table C2 in part C of the ESI† and Fig. 1 and 2), of dilute (c ≈ ≤0.005 M) IL solutions were analyzed in the framework of Barthel's low-concentration chemical model (lcCM)¹⁵ which describes successfully the thermodynamic and transport properties of diluted solutions. For the evaluation of molar conductivities, Λ, the lcCM uses the equationswhere Λ^∞ is the molar conductivity of the solute at infinite dilution, (1 − α) is the fraction of oppositely charged ions acting as ion pairs, and K_A is the standard-state (infinite dilution) ion association constant. y_±′ is the corresponding mean activity coefficient of the free ions, (y_±′)² = y′ + y₋′, κ_D is the Debye parameter, e_o is the proton charge, ε is the relative permittivity of the solvent, and ε_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₁ and J₂ 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_±′ (eqn (6) and (7)).

Analysis of the conductivity data was carried out by setting the coefficients S, E, J₁ and J₂ of eqn (5) to their calculated values¹⁵ using the properties of water, d(T),⁷⁰ε(T),⁷¹η(T)⁷² (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)¹⁵ 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.¹³

From extensive investigations on electrolyte solutions in different solvents it has been found¹⁵ 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₂) was fixed to 0.786 nm.

Then two-parameter fits were used to obtain the limiting values of molar conductivity, Λ^∞, 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¹¹ even the conductivity measurements, as one of the most established techniques, yield only a rough estimation for the association constant, K_A ≤ 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,74 From Table C3 in the ESI† (part C) it can be seen that the same trend in K_A was found at all temperatures.

3. Results and discussion

3.1. Ion mobility

Using a literature value for λ^∞₋(Cl⁻)⁷⁵ (Table C1 in part C of the ESI,†λ^∞(Cl⁻, 298.15 K) = 76.35 S cm² mol⁻¹), the limiting molar conductivity can be split into separate ionic contributions, λ^∞_i, where limiting ionic conductivities were estimated for all cations. Consequently, the diffusion coefficients D^∞_i were calculated by using the relationwhere 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⁻¹⁰ m² s⁻¹ can be estimated. For all investigated cations, values at 298.15 K are gathered in Table 1 (together with values of Λ^∞ and λ^∞₊ 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.

Table 1 Limiting values of molar conductivity, Λ^∞, and limiting values for cations, λ^∞₊, together with a comparison of calculated D^∞_+,calc, with experimental D^∞_+,exp values of the cation diffusion coefficients at infinite dilution of all investigated ILs^a

Cation	n c	Λ ^∞	λ ^∞+	D ^∞+,calc	D ^∞+,exp
a Units: Λ^∞, λ^∞+, S cm^2 mol^−1; D^∞+,calc, D^∞+,exp, 10^−10 m^2 s^−1. b Ref. 13.
[MIM]^+	0	124.06 ± 0.10	47.71 ± 0.05	12.8 ± 0.6	12.72 ± 0.05
[1,2-MMIM]^+	1	119.94 ± 0.05	43.59 ± 0.03	—	11.62 ± 0.03
[1,3-MMIM]^+	1	119.13 ± 0.01	42.78 ± 0.01	—	11.40 ± 0.01
[1,3-EMIM]^+	2	115.05 ± 0.05	38.70 ± 0.03	11.5 ± 0.3	10.31 ± 0.03
[1,3-BMIM]^+	4	109.34 ± 0.09 108.64^b	32.99 ± 0.05	9.2 ± 0.1	8.79 ± 0.05
[1,3-HMIM]^+	6	105.07 ± 0.10	28.72 ± 0.05	8.3 ± 0.8	7.65 ± 0.05
[1,3-OMIM]^+	8	102.71 ± 0.09	26.36 ± 0.05	7.5 ± 0.2	7.03 ± 0.05
[1,3-DeMIM]^+	10	100.73 ± 0.09	24.38 ± 0.05	5.8 ± 0.1	6.50 ± 0.05
[1,3-DoMIM]^+	12	99.90 ± 0.14	23.55 ± 0.07	6.3 ± 0.4	6.28 ± 0.07
Cl^−	—	—	—	—	20.35

---

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.

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.⁹ In 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⁻¹⁰ m² s⁻¹ and 7 × 10⁻¹⁰ m² s⁻¹, respectively). In the case of small [1,3-EMIM]⁺ cations the previously reported value was substantially larger (15 × 10⁻¹⁰ m² s⁻¹) and for larger [1,3-OMIM]⁺ a slower diffusion (2 × 10⁻¹⁰ m² s⁻¹) 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.⁷⁶ Even more, Goodchild and coworkers have found out that in water solutions of 1-alkyl-3-methylimidazolium bromides there are no immediately apparent aggregates in the systems with n_c = 2 and n_c = 4,⁷⁷ 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⁻¹⁰ m² s⁻¹ were derived for [MIM]⁺ and [1,3-DoMIM]⁺, which can be compared to 12.8 and 6.3 × 10⁻¹⁰ m² s⁻¹, 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.

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.

Fig. 5 1D-RISM PMFs between Cl⁻ and H7 ([MIM][Cl]), H9 ([1,2-MMIM][Cl]) and H1 for all other ILs.

3.2. Ion association

As can be seen in Tables 2a and C3 (part C of the ESI†), K_A values are small (∼2.5 ≤ K_A ≤ ∼6) but distinctly higher than those obtained using the same model recently used for alkali metal halides in water.¹

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

Cation	n c	(a)	(b)		(c)						(d)		(e)
		Exp	1D-RISM PMF X–Cl^−		1D-RISM PMF Hi–Cl^−						3D-RISM PMF		ΔG^expbind^d	ΔG^MDbind
		K A	Interacting cation atom X		Interacting cation atom						Cation's atom closest to the anion
					H1	H2	H3	H7	H8	H9
a Units: PMF, ΔG^expbind, ΔG^MDbind; kcal mol^−1. b Estimated errors are given in Table C3 in part C in the ESI. c Ref. 13. d Errors are denoted with error bars in Fig. 7.
[MIM]^+	0	4.9	−0.84 −0.87	C1 N2	−0.82	−0.71	−0.81	−0.70			−1.44	H7	−0.94	−2.7 ± 0.2
[1,2-MMIM]^+	1	2.5	−0.78	N2				−0.68	−0.69	−0.62	−1.51	H9	−0.54
[1,3-MMIM]^+	1	1.6	−0.65	C1	−0.70	−0.69	−0.69				−1.46	H1	−0.28
[1,3-EMIM]^+	2	3.1	−0.52	C1	−0.55	−0.68	−0.69				−1.38	H1	−0.67
[1,3-BMIM]^+	4	5.2 6.4^c	−0.53	C1	−0.54	−0.70	−0.70				−1.39	H1	−0.98	−0.8 ± 0.2
[1,3-HMIM]^+	6	2.6	−0.52	C1	−0.54	−0.69	−0.69				−1.40	H1	−0.57	−0.3 ± 0.1
[1,3-OMIM]^+	8	3.8	−0.51	C1	−0.53	−0.69	−0.69				−1.39	H1	−0.79
[1,3-DeMIM]^+	10	3.5	−0.51	C1	−0.53	−0.69	−0.69				−1.39	H1	−0.74
[1,3-DoMIM]^+	12	4.9	−0.51	C1	−0.52	−0.69	−0.68				−1.39	H1	−0.94	−0.9 ± 0.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⁷⁸ (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.⁹ At the same time, however, the anions are found to be present in the vicinity of the five-membered ring, both in pure ILs (see, for instance, ref. 79–83) and in their mixtures with water.^80,84–86 Obviously, this location of the anion relative to the cation is a result of the strong electrostatic interaction mediated by the ring. Moreover, similar to the fact that the polar heads of pure imidazolium-based ILs are known to form H-bonded ion pairs (see, for instance, ref. 79–83), it was found for their aqueous mixtures that the anions are also interacting with cations via hydrogen bonds mediated by ring hydrogen atoms.^6,84–86 Thus, we encounter C–H⋯Cl interactions with imidazolium ring protons acting as hydrogen-bond donors to the anion. Note that for [MIM][Cl] and [1,2-MMIM][Cl] one may also consider N–H⋯Cl interactions. More precisely, because of steric hindrances, only N2–H7⋯Cl for [MIM][Cl] and N2–H9⋯Cl for [1,2-MMIM][Cl] could occur (Fig. 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,85 It 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.

Fig. 6 2D-maps of the PMF between the cation and the anion for [MIM][Cl] (a), [1,2-MMIM][Cl] (b), [1,3-MMIM][Cl] (c), [1,3-EMIM][Cl] (d), [1,3-BMIM][Cl] (e), [1,3-HMIM][Cl] (f), [1,3-OMIM][Cl] (g), [1,3-DeMIM][Cl] (h) and [1,3-DoDMIM][Cl] (i). The most probable locations of the anion are indicated by red circles.

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 hydrogen-bonded interaction (Cl–H1⋯Anion) for ILs ranging from dimethyl to dodecyl as it was shown in the literature previously.^85,86,88 Moreover, 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).⁸⁹ The data of NIR spectroscopy demonstrate a similar pairing of chloride anions with 1-alkyl-3-methylimidazolium cations ([1,3-AMIM]⁺).⁸⁶ One 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 K_A for [MIM][Cl] > [1,2-MMIM][Cl] > [1,3-MMIM][Cl] 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⁻, ΔG^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 ΔG^exp_bind = −RT·ln K_A. The values of ΔG^MD_bind and ΔG^exp_bind are presented in Table 2e.

In the case of MIM⁺, ΔG^MD_bind seems to be overestimated, whereas for the other three cases with alkyl chains longer than n > 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 ΔG^exp_bind and ΔG^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 > 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,⁸⁹ where [1,3-MMIM]⁺ was found to associate more strongly with Cl⁻ than [1,3-OMIM]⁺.

Fig. 7 Comparison of binding free energies as obtained from the conductivity experiment, ΔG^exp_bind, and from MD simulations, ΔG^MD_bind, with the minimal values of 1D-RISM PMFs (for atom pairs N2–Cl⁻ (for [MIM]⁺ and [1,2-MMIM]⁺) or C1–Cl⁻ for all other IL solutions), W^RISM_bind.

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.⁷⁷ Despite the fact that experiments were carried out in diluted solutions (<0.005 mol dm⁻³) 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.

4. Conclusion

The ion mobility and association of dilute solutions of nine imidazolium based 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,2-MMIM][Cl] in water, were systematically studied using MD simulations, statistical mechanics calculations in the framework of the integral equation theory with one-dimensional (1D-) and three-dimensional (3D-) reference interaction site model (RISM) approaches as well as conductivity measurements.

From MD simulation, the diffusion coefficients close to infinite dilution, D^∞_+,calc, were calculated and compared with the experimental values D^∞_+,exp. The measured values were obtained from the cationic contributions, λ^∞₊, to the molar conductivity at infinite dilution, Λ^∞, 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.⁹ That 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 = ∼5.2 at 298.15 K being still in reasonable agreement with a reported value of 6.4.¹³ This 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 ≥ 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.

Acknowledgements

The financial support by the Slovenian Research Agency through Grant No. P1-0201 and the Bilateral project SLO-Ru (2014–2015) is gratefully acknowledged. This work was partially supported by COST Action CM1206. Access to the A*STAR Computational Resource Centre in Singapore, where the MD simulations were performed, is greatly appreciated. M. B.-R. wishes to express her sincere thanks to Prof. Dr J. Mavri for his kind support of this work.

References

1. J. Gujt, M. Bešter-Rogač and B. Hribar-Lee, J. Mol. Liq., 2014, 190, 34–41.
2. Q. Yang, H. Xing, B. Su, K. Yu, Z. Bao, Y. Yang and Q. Ren, Chem. Eng. J., 2012, 181–182, 334–342.
3. N. D. Khupse and A. Kumar, J. Phys. Chem. A, 2011, 115, 10211–10217.
4. V. V. Chaban and V. O. Prezhdo, J. Phys. Chem. Lett., 2013, 4, 1423–1431.
5. S. Feng and G. A. Voth, Fluid Phase Equilib., 2010, 294, 148–156.
6. C. Spickermann, J. Thar, S. B. C. Lehmann, S. Zahn, J. Hunger, R. Buchner, P. A. Hunt, T. Welton and B. Kirchner B, J. Chem. Phys., 2008, 129, 104505.
7. P. Yee, J. K. Shah and E. J. Maginn, J. Phys. Chem. B, 2013, 117, 12556–12566.
8. H. V. R Annapureddy and L. X. Dang, J. Phys. Chem. B, 2013, 117, 8555–8560.
9. B. L. Bhargava and M. L. Klein, Soft Matter, 2009, 5, 3475–3480.
10. S. Palchowdhury and B. L. Bhargava, J. Phys. Chem. B, 2014, 118, 6241–6249.
11. Y. Marcus and G. Hefter, Chem. Rev., 2006, 106, 4585–4621.
12. G. Hefter, Pure Appl. Chem., 2006, 78, 1571–1586.
13. M. Bešter-Rogač, A. Stoppa, J. Hunger, G. Hefter and R. Buchner, Phys. Chem. Chem. Phys., 2011, 13, 17588–17598.
14. A. Boruń, C. Fernandez and A. Bald, Int. J. Electrochem. Sci., 2015, 10, 2120–2129.
15. J. M. G. Barthel, H. Krienke and W. Kunz, Physical Chemistry of Electrolyte Solutions: Modern Aspects, Springer, New York, 1998.
16. Molecular Theory of Solvation, ed. F. Hirata, Kluwer Academic Publishers, Dordrecht, 2003.
17. S. A. Shah, Y. L. Chen, S. Ramakrishnan, K. S. Schweizer and C. F. Zukoski, J. Phys.: Condens. Matter, 2003, 15, 4751–4778.
18. A. Ben-Naim, Molecular Theory of Solutions, University Press, Oxford, 2006.
19. L. Sarkisov and P. R. Van Tassel, J. Phys.: Condens. Matter, 2008, 20, 333101.
20. M. V. Fedotova and M. F. Holovko, in Theoretical and Experimental Methods of Solution Chemistry (Monograph in Russian), Prospect, Moscow, 2011, pp. 68–152.
21. T. Luchko, I. S. Joung and D. A. Case, in Innovations in Biomolecular Modeling and Simulations, ed. T. Schlick, Royal Society of Chemistry, London, 2012, pp. 51–86.
22. S. Bruzzone, M. Malvaldia and C. A. Chiappe, Phys. Chem. Chem. Phys., 2007, 9, 5576–5581.
23. M. Malvaldi, S. Bruzzone, C. Chiappe, S. Gusarov and A. Kovalenko, J. Phys. Chem. B, 2009, 113, 3536–3542.
24. H. Nakano, J. Noguchi, T. Mochida and H. Sato, J. Phys. Chem. A, 2015, 119, 5181–5188.
25. B. M. Pettitt and P. J. Rossky, J. Chem. Phys., 1986, 84, 5836–5844.
26. S. W. W. Chen and P. J. Rossky, J. Phys. Chem., 1993, 97, 6078–6082.
27. K. Koga, X. C. Zeng and H. Tanaka, Fluid Phase Equilib., 1998, 144, 377–385.
28. A. Kovalenko and F. Hirata, J. Phys. Chem. B, 1999, 103, 7942–7957.
29. M. V. Fedotova and S. E. Kruchinin, J. Mol. Liq., 2012, 169, 1–7.
30. M. V. Fedotova and S. E. Kruchinin, J. Mol. Liq., 2013, 179, 27–33.
31. M. V. Fedotova and S. E. Kruchinin, Biophys. Chem., 2014, 190–191, 25–31.
32. M. V. Fedotova and O. A. Dmitrieva, Amino Acids, 2015, 47, 1015–1023.
33. M. V. Fedotova and O. A. Dmitrieva, New J. Chem., 2015, 39, 8594–8601.
34. I. Terekhova, E. Chibunova, R. Kumeev, S. Kruchinin, M. Fedotova, M. Kozbiał, M. Wszelaka-Rylik and P. Gierycz, Chem. Eng. Sci., 2015, 122, 97–103.
35. A. Eiberweiser, A. Nazet, S. E. Kruchinin, M. V. Fedotova and R. Buchner, J. Phys. Chem. B, 2015, 119, 15203–15211.
36. J. N. C Lopes, A. A. H. Padua and K. Shimizu, J. Phys. Chem. B, 2008, 112, 5039–5046.
37. C. A. N de Castro, E. Langa, A. L. Morais, M. L. S. M. Lopes, M. J. V. Lourenco, F. J. V. Santos, M. S. C. S Santos, J. N. C. Lopes, H. I. M. Veiga and M. Macatrao, et al. , Fluid Phase Equilib., 2010, 294, 157–179.
38. K. Shimizu, D. Almantariotis, M. F. Costa Gomes, A. A. Padua and J. N. Canongia Lopes, J. Phys. Chem. B, 2010, 114, 3592–3600.
39. M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys., 2000, 112, 8910–8922.
40. D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
41. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
42. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
43. H. J. C Berendsen, J. P. M. Postma, W. F. v. Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
44. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
45. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
46. A. Seduraman, P. Wu and M. Klähn, J. Phys. Chem. B, 2011, 116, 296–304.
47. D. Chandler and H. C. Andersen, J. Chem. Phys., 1972, 57, 1930–1937.
48. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, New York, 2006.
49. M. V. Fedotova and S. E. Kruchinin, J. Mol. Liq., 2013, 186, 90–97.
50. M. V. Fedotova and S. E. Kruchinin, J. Mol. Liq., 2011, 164, 201–206.
51. G. N. Chuev, M. Valiev and M. V. Fedotova, J. Chem. Theory Comput., 2012, 8, 1246–1254.
52. H. Hirata, P. J. Rossky and B. M. Pettitt, J. Chem. Phys., 1983, 78, 4133–4144.
53. A. Kovalenko and F. Hirata, J. Chem. Phys., 2000, 112, 10403–10417.
54. A. Kovalenko and F. Hirata, J. Chem. Phys., 2000, 113, 2793–2805.
55. A. Kovalenko and F. Hirata, J. Chem. Phys., 1999, 110, 10095–10112.
56. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
57. L. Lue and D. Blankschtein, J. Phys. Chem., 1992, 92, 8582–8594.
58. S. Miertuš, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117–129.
59. R. Cammi and J. Tomasi, J. Comput. Chem., 1995, 16, 1449–1458.
60. M. Cossi, V. Barone, R. Cammi and J. Tomasi, Chem. Phys. Lett., 1996, 255, 327–335.
61. C. Amovilli, V. Barone, R. Cammi, E. Cancès, M. Cossi, B. Mennucci, C. S. Pomelli and J. Tomasi, in Advances in Quantum Chemistry, ed. L. Per-Olov, Academic Press, New York, 1998, pp. 227–261.
62. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalman, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota,R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz and J. Cioslowski and D. J. Fox, Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford CT, 2010.
63. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361–373.
64. D. A. Case, V. Babin, J. T. Berryman, R. M. Betz, Q. Cai, D. S. Cerutti, T. E. Cheatham, T. A. Darden, R. E. Duke, H. Gohlke, A. W. Goetz, S. Gusarov, N. Homeyer, P. Janowski, J. Kaus, I. Kolossváry, A. Kovalenko, T. S. Lee, S. LeGrand, T. Luchko, R. Luo, B. Madej, K. M. Merz, F. Paesani, D. R. Roe, A. Roitber, C. Sagui, R. Salomon-Ferrer, G. Seabra, C. L. Simmerling, W. Smith, J. Swails, R. C. Walker, J. Wang, R. M. Wolf, X. Wu and P. A. Kollman, Amber 14, University of California, San Francisco, 2014.
65. A. Kovalenko, S. Ten-no and F. Hirata, J. Comput. Chem., 1999, 20, 928–936.
66. J. Barthel, R. Wachter and H.-J. Gores, in Modern Aspects of Electrochemistry, ed. B. E. Conway and J. O. M. Bockris, Plenum Press, New York, 1979, pp. 1–79.
67. J. Barthel, F. Feuerlein, R. Neueder and R. Wachter, J. Solution Chem., 1980, 9, 209–219.
68. M. Bešter-Rogač and D. Habe, Acta Chim. Slov., 2006, 53, 391–395.
69. O. Kratky, H. Leopold and H. N. Stabinger, Z. Angew. Phys., 1969, 27, 273–277.
70. E. F. G. Herington, Pure Appl. Chem., 1976, 48, 1–9.
71. B. B. Owen, R. C. Miller, C. E. Milner and H. L. Cogan, J. Phys. Chem., 1961, 65, 2065–2070.
72. L. Korson, W. Drost-Hansen and F. J. Millero, J. Phys. Chem., 1969, 73, 34–39.
73. D. Rudan-Tasic, T. Župec, C. Klofutar and M. Bešter-Rogač, J. Solution Chem., 2005, 34, 631–644.
74. M. Bešter-Rogač, C. Klofutar and D. Rudan-Tasic, J. Mol. Liq., 2019, 156, 82–88.
75. H. S. Harned and B. B. Owen, The Physical Chemistry of Electrolytic Solution, Reinhold, New York, 1958.
76. K. R. Ramya, P. Kumar, A. Kumar and A. Venkatnathan, J. Phys. Chem. B, 2014, 118, 8839–8847.
77. I. Goodchild, L. Collier, S. L. Millar, I. Prokeš, J. C. D. Lord, C. P. Butts, J. Bowers, J. R. P. Webster and R. K. Heenan, J. Colloid Interface Sci., 2007, 307, 455–468.
78. J. N. Canongia, J. Deschamps and A. A. H Pádua, J. Phys. Chem. B, 2004, 108, 2038–2047.
79. J. Ross and J. Xiao, Chem. – Eur. J., 2003, 9, 4900–4906.
80. A. Mele, C. D. Tran and S. H. De Paoli Lacerda, Angew. Chem., 2003, 115, 4500–4502.
81. C. Hardacre, J. D. Holbrey, S. E. J. McMath, D. T. Bowron and A. K. Soper, J. Chem. Phys., 2003, 118, 273–278.
82. B. L. Bhargava and S. Balasubramanian, Chem. Phys. Lett., 2006, 417, 486–491.
83. T. Köddermann, C. Wertz, A. Heintz and R. Ludwig, ChemPhysChem, 2006, 7, 1944–1949.
84. B. L. Bhargava and M. L. Klein, J. Phys. Chem. A, 2009, 113, 1898–1904.
85. M. Moreno, F. Castiglione, A. Mele, C. Pasqui and G. Raos, J. Phys. Chem. B, 2008, 112, 7826–7836.
86. B. Wu B, Y. Liu, Y. Zhang and H. Wang, Chem. – Eur. J., 2009, 15, 6889–6893.
87. B. Fazio, A. Triolo and G. Di Marco, J. Raman Spectrosc., 2008, 39, 233–237.
88. T. Takamuku, Y. Kyoshoin, T. Shimomura, S. Kittaka and T. Yamaguchi, J. Phys. Chem. B, 2009, 113, 10817–10824.
89. H. V. R. Annapureddy and L. X. Dang, J. Phys. Chem. B, 2013, 117, 8555–8560.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6cp05010g

---