Polarizable molecular dynamics simulations on the conductivity of pure 1-methylimidazolium acetate systems

The protic ionic liquid 1-methylimidazolium acetate is in equilibrium with its neutral species 1-methylimidazole and acetic acid. Although several experimental data indicate that the equilibrium favors the neutral species, the system exhibits a significant conductivity. We developed a polarizable force field to describe the ionic liquid accurately and applied it to several mixtures of the neutral and charged species. In addition to comparing single values, such as density, diﬀusion coeﬃcients, and conductivity, with experimental data, the complete frequency-dependent dielectric spectrum ranging from several MHz to THz can be used to determine the equilibrium composition of the reaction mentioned above. In the this to provide further insight and a better understanding of the interplay


Introduction
High demands are imposed on current battery and fuel cell generations, thus implying continuous research of more efficient and improved components. A promising approach is the application of ionic liquids (ILs), which are basically molten salts. [1][2][3] They have gained broad scientific interest during the last decades, not only because of countless ways of designing different candidates but also due to their specific and adjustable properties such as non-flammability, high thermal and electrochemical windows, as well as low volatility. 4 Among the multitude of different ILs, a subclass called protic ionic liquids (PILs) is especially interesting in terms of electrochemistry applications. 2 This family of substances is formed by a proton transfer from a Brønsted acid (AH) to a Brønsted base (B).
HA + B " A À + BH + (1) The available proton and ability to form hydrogen-bonded networks, 5,6 distinguishes PILs from aprotic ILs. Additionally, they usually have higher conductivities and lower viscosities.
Strictly speaking, only if the equilibrium in eqn (1) is shifted far to the right, the substance can be properly named an IL. The more the equilibrium shifts to the left, the more neutral species are in the mixture. Still, a remarkable difference to common ionic solutions remains, as the neutral species are in direct chemical relation to the ions. 7,8 Since it is not very likely that the proton transfer is fully complete, all four different species of eqn (1) may be present, resulting in the possible formation of (non-)neutral aggregates. 6 The Walden rule is often employed to access the ionicity of ILs, which states that the product of the molar conductivity and viscosity is constant. 9 On a log-log plot of these two properties, also known as the Walden plot, 1 M KCl solution is used as an ideal reference. 10 Many ILs lie below the reference line, indicating incomplete proton transfer and poor conductivity. Still, some ILs are remarkably close to the ideal line, although they are quite different from a dilute electrolyte solution in terms of ion correlation. 11 Angell and co-workers concluded that ILs with a DpK a value above 10 are good ionic liquids, where DpK a = pK a (BH + ) À pK a (AH). Regarding the main transport mechanisms, vehicle and Grotthus transport, the KCl-line in the Walden plot resembles vehicle transport by the fully dissociated species. On the other hand, Grotthus transport could exceed vehicle transport, as it is less hindered by diffusion, thus playing an important role, especially for ILs. 12,13 The DpK a for the IL 1-methylimidazolium (Im 1 H + ) acetate (OAc À ) of 2.2 indicates that the equilibrium favors the neutral species acetic acid (HOAc) and 1-methylimidazole (Im 1 ). 14 This finding is supported by several spectroscopic data also reporting on the dominance of the neutral species. 11,15 Nevertheless, a significant conductivity of about 3.3 mS cm À1 to 4.4 mS cm À1 was reported by several groups which argues for a significant amount of charged species. [15][16][17][18] In the light of these remaining ambiguities, this study aims to provide further insight and a better understanding of the interplay between high conductivity in PILs on the one hand and the dominance of neutral species on the other hand, by investigating the PIL [Im 1 H]OAc using polarizable molecular dynamics (PMD) simulations. While some methods, such as constant pH simulations, [19][20][21] exist to emulate proton transfer in MD simulations, these are not suitable for the studied PILs where several  hundreds of cations and anions are subject to (de-)protonation. Therefore, we model the equilibrium of HOAc + Im 1 " OAc À + Im 1 H + by sampling various ratios of charged and neutral molecules and comparing them to experimental data. In addition to comparing single values like density and conductivity, a more profound test is the juxtaposition of the frequencydependent dielectric spectrum ranging from several MHz to THz and covering all local and collective motions of the liquid.

Methods
Classical force fields have fixed charges located at the nuclei positions. These force fields are additive since the total electrostatic energy emerges from the sum of the pairwise interactions. While this description is often sufficient for reproducing structural properties, they fail at describing dynamical properties of ILs, 22,23 as flexible electronic charge distributions are essential to properly capture the response of the charged system to fluctuations in the local environment. [23][24][25] The induced dipoles implemented in polarizable molecular dynamics (MD) simulation, however, may model this essential electronic flexibility. 23,26

Parametrization process
The starting point for the parametrization was the CHARMM General Force Field (CGenFF). 27 First, the initial parameters of the molecules were obtained by using the PARAMCHEM website. 28 Molecules, which are not already fully featured in the CGenFF, are built based on analogy. Subsequently, the recently published FFParam package 29 was used to convert the atom types of the molecules to Drude force field compatible ones since polarizable MD was used throughout the whole study. The optimization procedure of FFParam covers the inspection of structural similarities and fitting of electrostatic as well as bonded parameters in the given order. After each step, checks were made to ensure the quality of the parameters by comparison to the reference data. The whole procedure was done in a self-consistent-field-like approach.
The quantum mechanical (QM) reference data was generated using Gaussian09 30 as well as Psi4 31 for dipoles and polarizabilities since the more efficient RI-MP2 algorithm is implemented. According to the Drude Generalized Force Field (DGenFF) standard CHARMM procedure, all intermolecular potentials were parametrized using water. These water interactions were calculated at a MP2/aug-cc-pVDZ, dipoles and polarizabilities at MP2/cc-pVQZ, and PES scans of important bonded parameters at a MP2/aug-cc-pVDZ level of theory.
On top of that, force matching was applied to optimize force constants and equilibrium distances further. The total force -F MM of a force field is the sum of intra-and intermolecular interactions. The superscript MM denotes that the forces are from a classical molecular mechanics force field, and Y represents the physical parameters.
The goal of force matching, as portrayed in ref. 32, is to fit intramolecular parameters to high level forces -F QM , e.g. derived from QM calculations. This is done via a Bayesian formalism which minimizes the negative log-likelihood and thus the force residuals: The last term subtracts the intermolecular interactions from the QM-derived forces before fitting, as the difference between intramolecular forces of the empirical force field and the QMderived forces should be minimized. z t( j) is an estimate of the mean square error of the force for atom type t( j).
The generation of force data was done by extracting conformations from a CHARMM 33 trajectory at different time steps and calculating the forces using GAUSSIAN09. The calculations were done at PBE1PBE/cc-pVDZ level of theory and using the FORCE keyword to request the calculation of forces. Since the results were still not satisfying, additional charge fitting was done. Instead of the standard Mulliken method for calculating partial atomic charges, calculation of the charges from electrostatic potentials using a grid-based method (CHELPG) was utilized. They depend less on the underlying level of theory. 34,35 Due to the introduction of Drude particles to a MD simulation, van der Waals interactions between induced dipoles are modeled explicitly. To account for this fact in the simulation, the depth of the Lennard-Jones potential, e LJ b is scaled. Otherwise, London forces would be counted twice, as they are also covered in the general Lennard-Jones term of the force field. This adjustment was made using the formula of Vlugt and coworkers, 36 which scales the well depth for each atom b in dependence of the polarizability a b . Atoms with the highest polarizability, max (a b ), are scaled by the value of l and all other atoms between l and 1. Using this formula, it is also ensured that nonpolarizable atoms are unaffected by the adjustment. The final bonded and non-bonded parameters of the different molecules can be found in the ESI † (Tables S1-S5).

Computational setup
As our goal is to determine the composition of the equilibrium depicted in Fig. 1, the nine systems in Table 1 were studied.
All systems studied contained a total of 1000 molecules and fulfilled charge neutrality, as always the same amount of cations and anions was used. The relative content f IL of ionic/neutral pairs varied between the pure systems containing 500 pairs of only neutral molecules Im 1 and HOAc to 500 ion pairs of only Im 1 H + and OAc À . Three replicas of each system were generated to increase the statistical accuracy of our results.
As Drude oscillators emulate the induced forces in our polarizable simulations, only non-hydrogens can be made polarizable. Nevertheless, the polarizabilities of the hydrogens are added to the next atom to which they are bonded. Drude particles were assigned a mass of 0.4 u and a force constant k = 1000 kcal mol À1 Å À2 . To obtain stable simulation, the maximum distance for the mobile Drudes was set to 0.2 Å. Lennard-Jones interactions were reduced according to eqn (4) using l values between 0.2 and 0.5. A previous study 37 showed that a l value of 0.4 yielded best results. However, in this work, best results are recorded for a l of 0.25.
Trajectories were generated by first getting the initial configuration of the system using Packmol. 38 After an initial minimization using CHARMM 33 a timestep of 0.1 fs for 20 ps at 300 K was used to obtain stable simulations. Afterwards the timestep was changed to 0.5 fs. The temperature was increased to 500 K to ensure proper mixing. After 2 ns the temperature was decreased to 300 K again for another 3 ns. A Monte-Carlo barostat at 1.0 atm was used for these NPT runs. The average boxlength of the last nanosecond was then used for the subsequent 50 ns NVT production run. All simulations were done using periodic boundary conditions with an initial boxlength of 48 Å. The final boxlengths are given in Table S7 (ESI †).
The OpenMM package 39 was used for generating trajectories with a time step of 0.5 fs. Integration was done using the velocity Verlet algorithm with the temperature-grouped Dual-Nosé-Hoover thermostat as described in ref. 40 and 41, which implements an additional temperature group for the center-ofmass translations. The non-Drude particles were kept at 300 K and the Drudes at 1 K. The collision frequency was set to 10 ps À1 for the atoms and to 200 ps À1 for the Drude particles. Keeping the vibrations of the Drude particle at a extremely low temperature results in an effective energy minimization on the fly, i.e. finding the optimal positions of the Drude particles.
Electrostatic interactions were treated using the Particle Mesh Ewald algorithm with a cut-off distance of 11 Å and a switch distance of 10 Å to truncate the interactions smoothly. Bonds containing hydrogens were fixed using the constraints = HBonds option. The error tolerance d was set to 0.0001. From the given error tolerance and cut-off, the number of mesh nodes is specified internally, via #mesh = 2/3add À1/5 with d being the box length in each dimension, and a ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi À logð2dÞ p r cutoff . Simulations were run on the CUDA platform using the default setting of mixed precision, which means that forces are calculated in single precision, and integration is done in double precision.
The input files for OpenMM were provided as CHARMM psf and crd files, as well as the Drude General Force Field (DGenFF) 42

Nernst-Einstein conductivity r NE
Charge transport depends on the number of charge carriers and their mobility. Quite generally, the mobility of a molecular species k is characterized by its diffusion coefficient i ðtÞ Àr i ð0Þ j jusing an index i for each molecule belonging to species k. The number of charge carriers correlates with the number density r k = N k /V. As a result the conductivity of a system can be estimated via the Nernst-Einstein equation Here, q k is the molecular charge of the species k. 45 However, this evaluation neglects all correlations between the motions of the charged species and consequently overestimates the static conductivity s(0). This is usually corrected by a factor (1 À D), called Haven ratio. However, our results show that this fudge Fig. 1 Equilibrium between the IL 1-methylimidazolium acetate [Im 1 H]OAc and its neutral analogues 1-methylimidazole Im 1 and acetic acid HOAc. Table 1 Systems under investigation. All systems contain a total of 1000 molecules distributed between the species Im 1 H + , OAc À , Im 1 , HOAc. For each system three replica were simulated. Please note that these tabulated systems belong to one specific l-value and have to be repeated for each additional l-value. The ratio f IL is the number of the charged molecules divided by the total number of molecules factor is a function of the composition f IL rendering this correction even more questionable.

Frequency-dependent conductivity r(x)
Consequently, the static conductivity is usually computed via the collective translational dipole momentM J ðtÞ ¼ P i q i Ár i ðtÞ which describes an intermolecular dipole moment using the molecular charges q i and the center-of-mass positions r i (t) of the molecules i. Applying an Einstein-Helfand approach, 46 the static conductivity yields with Actually, we showed 47 that this correlation function correlates with the mean-squared displacement of the collective translational dipole moment: since the electric current -J(t) = d -M J (t)/dt is the time derivative of the collective translational dipole moment. In principle, the static conductivity can also be obtained from eqn (9) but with less statistical accuracy.

Spectrum of the generalized dielectric constant R 0 (x)
The charged species also contribute to the frequencydependent spectrum of the generalized dielectric constant S 0 (o), in particular to the frequency-dependent dielectric conductivity W 0 (o) which is defined as However, this property is experimentally not directly accessible and should not be confused with the conductivity s(o) in eqn (9). The dielectric conductivity W 0 (o) usually characterizes cage librations or any intermolecular vibrational motion between charged species and not the directed motion of the ions, which is characterized by the static conductivity s(0). 47 The second significant contribution to the generalized dielectric constant stems from the permittivity e(o) and describes the collective rotation of the molecular dipole moments. The collective rotational dipole moment -M D (t) is the sum of molecular dipole momentsm i (t). As -M D (t) can be computed for each species k separately, a decomposition of the permittivity into molecular contribution is feasible, i.e. eðoÞ ¼ P k e k ðoÞ.
The comparison between the computational and experimental spectra are done on the level of the generalized dielectric constant S 0 (o), 47,48 which consists of the overlapping contributions from the permittivity e(o) and dielectric conductivity W 0 (o) and additionally the high-frequency limit of the dielectric constant e N .
Please keep in mind that experimenter uses the frequency n = o/(2p) instead of the pulsatance o.

Voronoi-shell resolved potential of mean force
The potential of mean force PMF kl is a measure for the aggregation of a species l in the proximity of a species k. Commonly, it is defined as a function of the distance from the central molecule. However, as our molecules are not spherical but have significant anisotropic shapes, we apply here a shell-resolved PMF kl (s) which can be defined for each non-spherical solvation shell s: 49 The ratio of the concentration c l (s) of species l in solvation shell s and the bulk concentration c l resembles the distance resolved radial distribution function g(r) which is also a ratio of local and global mass density.
Here, we restrict ourselves to the first solvation shell s = 1, which can be unambiguously determined by Voronoi tesselation. 50,51 The coordination number N kl = N lk counts the number of molecules of species l in the first solvation shell of molecules of species k. The volume V k (s = 1) is the averaged sum of molecular volumes of all molecules in the first solvation shell of a molecule of species k. N l is the total number of molecules of species l in the simulation box, and V is the box's volume.

Results and discussion
The top panel of Fig. 2 depicts the overall mass density r as a function of f IL . The mass density r(f IL ) increases linearly with increasing content of charged molecules. The Coulomb interaction between the oppositely charged ions is stronger, leading to closer distances between the ions. As a result, the overall density of molecules is increased. The density line r(f IL ) crosses the experimental density 52 of r = 1.07 g cm À3 at f IL = 45%. Lowering the Lennard-Jones scaling factor l (see eqn (4)) has only marginal consequences for the density.

Diffusion
The diffusion coefficients of all four species are calculated using the Einstein relation in eqn (5)  The simulation systems become more viscous with an increasing number of charged molecules because of the overall increased Coulombic interactions. Therefore, the diffusion coefficients of the charged species Im 1 H + (filled pentagons) and OAc À (filled stars) as well as the neutral species Im 1 (empty pentagons) and HOAc (empty stars) decrease exponentially with f IL . Im 1 and HOAc have higher diffusion coefficients than their charged analogs which can also be explained by the additional Coulomb interactions between the molecular charge and the surrounding dipoles of the ions.
Since the proton transfer HOAc + Im 1 " OAc À + Im 1 H + is very fast compared to the NMR measurements, experimental diffusion coefficients always contain contributions from the charged and neutral species. 53 However, simulation mixtures containing f IL = 10% to 30% charged molecules are close to the experimental data of ref. 18. Reducing the l-value from 0.40 (green curves) to l = 0.25 (orange curves) improves the agreement with the experiment. Applying smaller l-values than 0.25 does not improve the results but may increase stability issues.
As already mentioned in the Theory section, the conductivity is determined by the number of charge carriers and their mobility. In Fig. 3, the Nernst-Einstein conductivity s NE based on eqn (7) is displayed as a function of f IL . The bell-shaped behavior can be described by an empirical fit: The corresponding fit parameters for l = 0.25 and l = 0.40 are given in Table 2. This fit can also be used to determine the maximum conductivity max(s NE ) and the respective Although the maximum conductivity for l = 0.25 is higher than that for l = 0.40, the respective f IL coincides.
Up to a fraction f IL of roughly 0.2, the increased number of charge carriers promotes the conductivity of the system. Beyond f IL = 0.2, the Nernst-Einstein conductivity s NE drops down significantly as the diffusion coefficients of species decrease exponentially with f IL . This is the first indication that only a small number of charged molecules (E20%) are sufficient to realize a significant conductivity in the system. However, the Nernst-Einstein approach neglects all correlations between the motion of the molecular species. Consequently, we continue our analysis with the collective static conductivity s(0).

Conductivity
The Einstein-Helfand approach of eqn (8) can be applied to the unfolded trajectories of MD simulations straightforwardly. As the mean-squared displacement of the collective translational dipole moment is usually noisy, we averaged the corresponding hDM J 2 (t)i of the three replicas before determining the slope between 500 ps and 2 ns. The obtained s(0) as a function of f IL   is shown in Fig. 3 in comparison to the s NE (f IL ). The bellshaped curves are again fitted by eqn (16). As visible from the fit parameter in Table 2 the maximum values of s(0) are roughly a quarter from the s NE -values arguing for a significant D-value in eqn (7). This indicates a strong correlation between the (charged) species. Additionally, this overestimation on the basis of s NE prohibits the discrimination which l-value suits better the experimental data. Moreover, the position of the maximum has shifted to f IL (max(s(0))) = 0.25-0.28. This shift casts additional doubts on the usefulness of s NE for the interpretation of conductivity phenomena in ionic liquids as not only the absolute values are overestimated by a factor (1 À D) but also this factor D(f IL ) changes drastically with the composition f IL (see ESI † for further information).
The maximum conductivity s(0) = 3.1 mS cm À1 can be obtained for the trajectories with l = 0.25 at f IL = 0.28. This value is still lower than 3.3 mS cm À1 to 4.4 mS cm À1 (light green area in Fig. 3) found in experiment. [15][16][17][18] On the one hand, small amounts of water increase mobility and, hence, conductivity in experiments. Various water and impurity contents may also explain the variance of the experimental s(0)-values indicated by the green area in Fig. 3. These impurities are not present in the simulations. On the other hand, our simulations do not allow proton transfer. As a result, our polarizable MD simulations mainly monitor vehicle transport. The simulation's lack of explicit Grotthus transport may also explain the slightly lower conductivity than experimental data.
The conductivity near the maximum can also be described by the Bahe-Varela theory 54-56 with x IL = f IL /f IL (max (s(0))). Using the values of Table 2 results in the black dash-dotted parabola in Fig. 3. This finding indicates that the equilibrium between the charged and neutral species behaves like a solution of the charged species dissolved in the neutral species. However, the inclusion of quadratic and cubic terms of (1 À x IL ) as suggested by ref. 56 leads to a prediction of a minimum conductivity at f IL close to unity (see ESI † for further information). Consequently, we stick to our empirical fit function in eqn (16).

Dielectric spectroscopy
Based on our findings for the conductivity and density, a mixture of roughly 30% charged Im 1 H + + OAc À and 70% neutral Im 1 + HOAc seems most probable. This finding is also in agreement with AIMD simulations, which found concentrations of 70% to 80% neutral species in a [IM 1 H]OAc system. 57 However, this optimal ratio is found by comparing single values of the diffusion coefficients, the conductivity, and the density in our study. A more profound check can be made on the basis of the dielectric spectrum S 0 00 (o) as it covers the rotational and translational motion of the system over several orders of magnitude in frequency o.   Table 1. A continuous shift to higher frequencies for increasing content of neutral species is observed. The amplitude increases as well, except for the neutral system. The mixtures with f IL = 10% to 30% match the experimental spectrum best. In particular, the f IL = 0.30 mixture agrees very well in terms of peak height as well as the position of the maximum. It can be clearly seen that neither the pure IL nor the neutral species reproduce the experimental spectrum satisfactorily. It is also impossible to reproduce the experimental spectrum with a linear combination of the pure neutral and the pure ionic system. This finding clearly demonstrates that neutral and charged species interact to a vast extent, changing the individual components' properties in the mixture. Fig. 4(b) demonstrates the exhaustive overlap between the permittivity e(o) and dielectric conductivity W 0 (o) for f IL = 0.30. In addition, the conductivity parabola s(0)/o for the experimental data and the mixture at f IL is shown in Fig. 4(b) as dashdotted line. As already deduced from Fig. 3, the computational static conductivity s(0) at f IL = 0.30 matches the experimental one best. Shifting to lower f IL would increase the gap between the orange dash-dotted line (computational) and the gray dashdotted line (experimental conductivity). The static dielectric constant S 0 (0) of 31.3 (see Table 3) in our simulations at f IL = 0.30 agrees well to the experimentally found value of 33.3. At f IL = 0.20, the static dielectric constant is slightly overestimated with 36.7.
The permittivity e(o) in Fig. 4(c) emerges from the collective rotation of the molecular ions. The maximum amplitude does not change very much with varying mole fractions of the ions. Only the shift to lower frequencies with increasing ion content (as already noticed for S 0 (o)) is found again. Taking a closer look at e(o) for f IL = 0.30 in Fig. 4(d) reveals that all molecular species contribute to the permittivity. The largest share comes from the neutral 1-methylimidazole (red area, e Im1 (0) = 8.2, Table 3) which has also the broadest peak. Interestingly, the second most important contribution stems from the acetic acid molecules (gray area, e HOAc (0) = 4.4). This is counter-intuitive as acetic acid has a low molecular dipole moment m HOAc of 1.7 D (see ESI †). Furthermore, pure acetic acid is known to form dimers with hydrogen bonds between the hydroxy-hydrogen of one acetic acid molecule with the carbonyl oxygen of the other acetic acid molecule (see Fig. 5).
Thus, the dipoles of the two hydrogen-bonded acetic acid molecules should cancel out, which results in a dielectric constant of 6.2 in pure acetic acid. However, in our mixtures, the acetic acid prefers to interact with the 1-methylimidazole as shown by the first shell potential of mean force PMF(s=1) = À0.36 kJ mol À1 . In contrast, the first shell potential of mean force PMF(s=1) between two acetic acid molecules is +0.43 kJ mol À1 at f IL = 0.30, meaning an overall depletion of acetic acid molecules in the first solvation shell around acetic acid. Additionally, the formation of hydrogen bonds and especially dimers of HOAc molecules is reduced. About 5.6% of the acetic acid molecules have at least one hydrogen bond to another acetic acid molecule, of which only 2.2% form HOAc dimers at f IL = 0.30. The maximum distance between the hydrogen donor and oxygen acceptor atom was set to 3 Å to define a hydrogen bond in our analysis.
Obviously, as shown by the overlap of e(o) and W 0 (o) in Fig. 4(b), the amplitude S k of the peak cannot be used to evaluate an effective dipole moment m eff,k based on the Cavell equation 58 S k p c k Á m eff,k 2 (18) The maximum amplitude at o = 2.5 GHz for the f IL = 0.30 mixture comes from equal contributions of e(o) and W 0 (o). The latter only depends on the translation of the ions but not their dipoles. However, a decomposition into permittivity contributions e k (0) of the molecular species k is possible for the MD trajectory. The respective contributions are listed in Table 3. Thus, we are not only able to remove the translational contribution but analyze each species separately. Nevertheless, there is no correlation between e k (0) and S k computed from the Cavell equation (see ESI †). The dielectric conductivity W 0 (o) in Fig. 4(e) shows an opposite trend to higher frequencies with increasing f IL but this trend is less pronounced compared to that of the permittivities. Interestingly, the amplitude decreases with increasing ion content f IL which is surprising at first sight as only the ions Im 1 H + and OAc À can contribute. However, the interactions between a central ion and its neighboring ions are stronger at lower f IL as characterized by the first shell potential of mean force PMF(s=1). This PMF(s=1) of a central 1-methylimidazolium (blue filled pentagons) and a central acetate (orange filled circles) in Fig. 6 is most negative at low f IL showing an aggregation of oppositely charged ions around a central ion. These strong interactions go along with librational motion between the ions leading to significant contributions to W 0 (o). With increasing ionic concentrations, these ion clusters transform into an ionic network in which the neutral molecules penetrate. As a result, the local concentration c l (s = 1) of the oppositely charged ions converges to the bulk density c l of these ions, and fewer particular ion clusters are found. The dashed lines in Fig. 6 correspond to a Gaussian fit f (f IL ) = aÁ exp(Àbf IL 2 ) + c of the PMF. The corresponding fit parameters are also given in Fig. 6. Interestingly, the exponential b-factor of 10.7 and 11.4 is very close to 10.9 found for the Nernst-Einstein conductivity as a function of f IL in Table 2.

Conclusions
The position of the equilibrium HOAc + Im 1 " OAc À + Im 1 H + is not so easy to determine as various experiments came to different conclusions. By comparing density, diffusion coefficients and conductivity of various mixtures modelled by polarizable molecular dynamics with the experimental data, a ratio of 20% to 30% ionic and 80% to 70% neutral molecules seems most probable. In addition to comparing these single values, we also computed the frequency-dependent dielectric spectrum of all these mixtures and compared the computational results with the experimental spectrum over several orders of magnitude in frequency. It turned out that the mixture with 30% ions fits the experimental data best. In addition, the dielectric spectrum cannot be reproduced by a linear combination of systems containing only ionic or only neutral molecules. Consequently, there have to be strong correlations between all species altering their individual dynamic properties. For example, the shell-resolved potential of mean force PMF(s=1) decays Gaussian-like with increasing ion content f IL for the interaction of a central ion and its oppositely charged neighboring ions.

Conflicts of interest
There are no conflicts to declare.