Polarizable MD simulations of ionic liquids: How does additional charge transfer change the dynamics?

Abstract

Both experimental and computational evidence exist that Coulomb interactions between the molecular ions in ionic liquids are significantly damped by almost a factor of two. This circumstance is often used to justify charge scaling. However, as polarizable MD simulations are also capable of explaining the reduced Coulomb interaction between the ionic liquid ions [C. Schröder, Phys. Chem. Chem. Phys., 2012, 14, 3089], the question arises, if the reduced Coulomb interactions are due to a charge transfer between the molecules or due to an overall effect of induced dipolar interactions. We aim to contribute to this discussion using polarizable MD simulations of 1-ethyl-3-methylimidazolium tetrafluoroborate including a new model for treating charge transfer between the cations and anions. The diffusion time scales are not changed significantly with the inclusion of charge transfer, but individual ions show a strong dependence on charge transfer amounts. Ions which have transferred more charge, and have a charge with a smaller magnitude, diffuse slower. The charge transfer model shows a slightly larger conductivity, despite having smaller charges, and shows a much stronger contribution of the anions to the conductivity. With charge transfer, the anions become the dominant species for charge transport, while the polarizable models show a roughly equal contribution from the anions and the cations.

---

Room temperature ionic liquids (ILs), composed of an organic cation and an organic or inorganic anion, are an interesting class of solvents, with applications including electrolyte solutions in batteries and capacitors, lubricants, and solvents for synthesis and catalysis.^1,2 Significant effort has gone into the simulation of these liquids and a variety of potentials have been developed for ILs.^3–16 Non-polarizable models, similar to what is commonly used in the simulation of aqueous solutions and proteins, can accurately reproduce a number of properties,^4–6 but tend to underestimate diffusion rates.^7,10,16–18 Polarizable models for ILs using point inducible dipoles,^7,9,17 fluctuating charges,^13,14 and Drude oscillators¹¹ give better dynamical properties and it is not surprising given the large electric fields created by the ions that polarizability is important.¹⁹ A group of non-polarizable potentials uses reduced charges on the ions, so that each ion has a charge of ±0.7 e to ±0.9 e.^{4,8,12,18,20,21} These models, which will have weaker Coulombic interactions, have faster dynamics, closer to the polarizable models and experimental results. This charge reduction is taken as an adjustable parameter^8,18,20 or as an easy method to account for polarizability.²² A comparison of charge scaled to polarizable models has shown that the charge scaled approach does not completely capture the dynamics of ILs, especially at the local level.^18,19

Scaled charges have been used in models for aqueous ions.^22–25 For both, ILs^{4,12,21,26–28} and aqueous ions,^29–33 the reduction in charge in the scaled charge models is consistent with the amount of charge transfer (CT) between particles determined from ab initio calculations. Computational studies have shown that the charge transfer between particles influences the dynamics of water in salt solutions^34–36 and ILs.³⁷ Experimental studies with X-ray photoelectron spectroscopy (XPS)²⁸ and two-dimensional infrared (2D-IR) spectroscopy³⁸ indicate the interactions between the ions in ionic liquids are influenced by CT. The importance of CT for ion-water and water–water interactions are seen in IR,^39–41 Raman,^42–44 and X-ray absorption⁴⁵ spectroscopy.

The experimental and ab initio studies both indicate the significance of CT on the properties on ILs. In this paper, we describe a new model for treating CT effects that builds on earlier work developing CT models for water^46–48 and ions.^32,49–52 What makes this CT approach different from the charge scaling methods for ionic liquids is that our method transfers charge between particle pairs based on their separation. The charge of a particle will depend on its local environment and will fluctuate. The amount of charge transfer is determined from ab initio calculations. The model is also polarizable, using the Drude oscillator approach. Using this model, we will examine how charge transfer influences the structure and dynamics of the IL with the 1-ethyl-3-methylimidazolium cation, C₂mim⁺, and the tetrafluoroborate anion, BF₄⁻.

1 Methods

1.1 The charge transfer model

Charge transfer in the model is treated using a method adapted from the Discrete Charge Transfer (DCT) method developed previously.⁴⁶ Previous DCT applications have been applied to molecules with interacting through hydrogen bonds, e.g. water,^46–48 and the hydroxide or hydronium ions,⁵¹ or monatomic ions.^32,50 Here, we extend the approach to polyatomic ions. In this approach, charge is transferred between two ions, i and j, based on the distance between atoms on the two ions,where the first sum is over all atoms β on the molecular ion i, the second sum is over all atoms γ on the molecular ion j, and r_βγ is the distance between atoms β and γ. Total charge is conserved, so Δq_ij equals −Δq_ji.

Following the other DCT models, the charge transferred, δq_βγ, is taken to be a constant, δq₀, at small separations, and to smoothly go to zero at large separations, as given by

and depicted in Fig. 1.

Fig. 1 Charge transfer scheme between anions and cations (including definition of the atom types).

The total charge q_i(t) on a molecular ion i is then the formal charge of that ion, q⁰_i, plus all the charge that is transferred to other ions,

There are a number of ways the charge could be distributed among the atoms on the ion. We choose to distribute the charge equally among all the atoms, so that the charge of atom β on ion i, iswhere q⁰_iβ is the charge of the atom before charge transfer and N_i is the number of atoms on ion β. By distributing the charge evenly, an ion's dipole moment is not changed by charge transfer and the polarizability and charge transfer effects are distinct. An ion's dipole will depend on the origin of the coordinates used but with the ion's center as the origin, the dipole moment before and after charge transfer will be the same, fromwhere δq_i is the total charge transferred to each atom of ion i and μ⁰_i is the dipole moment of the molecule before charge transfer. Using the ion's center as the origin means equals zero, so μ_i will equal μ⁰_i.

Charge is only transferred between ions of opposite charge and it transferred from the anion to the cation. To keep things simple, only heavy atoms on the C₂mim⁺ cation and fluorine atoms on BF₄⁻ anion are used in eqn (1) and the parameters (δq⁰, r₁, and r₂) are taken to be the same for all atom pairs. (The charges of all atoms are modified through eqn (4).)

In addition to the change in the electrostatic interactions, there is a charge transfer energy E_CT, for each amount of charge transferred, given by

E_CT(t) = −μ_CT·|δq_βγ(t)| (6)

This term is energetically favorable, corresponding to moving electron density to lower energy unoccupied orbitals,⁵³ and compensates for the increased energy due the reduction in the Coulombic energy from charge transfer.

This CT method introduces four parameters (δq⁰, r₁, r₂, and μ_CT). Density functional theory (DFT) calculations with the M06-2X functional⁵⁴ and aug-cc-pVTZ basis sets, using the program Gaussian,⁵⁵ determined the distance dependence of the amount of charge transfer between pairs.

Pair geometries were generated by varying the distance along a coordinate in the plane of the imidazolium ring in the direction outward from the CR atom and in a direction perpendicular to the ring, moving outward from the ring center (Fig. 1). The acidic CH points directly towards the B atom along the bisector of the F–B–F angle. Other pair geometries were taken from liquid simulations. Charge was determined both with the Bader atoms in molecules (AIM) method⁵⁶ using the AIMALL program,⁵⁷ and the CHELPG electrostatic potential fitting method,⁵⁸ as implemented by Gaussian. Comparing to these results it was found that an upper limit of charge transfer, r₂, of 4.8 Å is consistent with the distance at which the DFT calculations find that the charge transfer between pairs goes to zero. The charge transfer parameters δq⁰ and r₁ were chosen to reproduce the liquid states charges as calculated by Mondal and Balasubramanian.²¹ Using configurations of the liquid generated from molecular dynamics simulations, charges were determined, using the density-derived electrostatic and chemical (DDEC)⁵⁹ charge partitioning method, which is based on AIM. This method gave ions with average charges equal to ±0.79 e.²¹ Values of δq⁰ = 0.0050 e and r₁ equal to 4.0 Å were found to reproduce the Mondal and Balasubramanian results, giving an average charge of ±0.77 e. The final parameter, μ_CT, was adjusted to give the correct liquid density at 298 K. This gives a value of μ_CT equal to 162.0 kcal mol⁻¹e⁻¹.

1.2 Polarizable model using Drude oscillators

The Drude model treats polarizability using the Drude oscillator method, in which charge are placed on harmonic springs attached to heavy atoms.¹⁹ The Lennard-Jones parameters, partial charges q^FF_iβ, and bonded parameters for bond stretches, bends, and torsions, are all taken from the polarizable force field of ref. 60 which is based on the non-polarizable force field of Canongia Lopes and Pádua.⁵ In the polarizable model, the dipoles are induced by the charges, as well as the other dipoles. Since our CT model has smaller charges, the induced dipoles will also be smaller. To compensate for this, the polarizabilities α_iβ were increased by a factor of s = 1.04 by increasing the Drude charge. With this change, the dipole moment distributions of the models with and without charge transfer almost coincide (see Fig. 2), making comparisons between the two models easier. For additional comparison, the dipole distributions of a scaled charge model are also depicted. For the dipole calculation, the charge scaling was removed as this procedure can be considered as a poor man's way to take an average polarizability into account.^18,22 The corresponding shift of the maxima of the dipole distributions between using the scaled/full charges for analysis are shown as horizontal black arrows in Fig. 2.

Fig. 2 The distribution of dipole moments for the pure polarizable model (ref. 60, solid line), the polarizable CT model (dashed line) and the non-polarizable scaled charge model (dotted line). In case of the non-polarizable scaled charge model, the black arrows depict the shift due to using full charges for analysis resulting in the dotted blue and orange distributions for the cations and anions, respectively.

All mobile Drude particles have mass of 0.2 amu, which was subtracted from the mass of the heavy atom associated with that particle. The partial charge q^D_iβ of the mobile Drude particle is defined by:

and listed for the various atom types in Table 1 for both the pure polarizable model and the polarizable CT model. The CT model has polarizabilities a factor of 1.04 larger than the polarizable model (vide supra). The Drude charge q^D_iβ for the polarizable atom iβ is subtracted from the partial charge q^FF_iβ given in the force field, i.e. q⁰_iβ = q^FF_iβ − q^D_iβ. Thus, the Drude pair of the mobile particle and the one located at the position of the polarizable atom generates an induced dipole for that atom and the total charge is conserved. The force constant k^D for all mobile Drude particles is 500 kcal mol⁻¹ Å⁻² and interactions between induced dipoles of atoms connected by bonds or angles are neglected in analogy to the treatment of Coulomb interactions.

Table 1 Partial charges q⁰_iβ and Drude charge q^D_iβ for the polarizable and CT models. The polarizability of the hydrogens is merged to that of the heavy atom to which they are bonded. The atom types are displayed in Fig. 1

	Atom type	q ^FF iβ [e]	α iβ [Å^3]	q ^D iβ [e]	α iβ [Å^3]	q ^D iβ [e]
			Pure pol.		Pol. CT
BF4^−	Boron	0.8276	1.186	−1.336	1.233	−1.363
	Fluorine	−0.4569	0.400	−0.776	0.416	−0.792
C2mim^+	NA	0.150	1.123	−1.300	1.168	−1.326
	CR	−0.110	1.578	−1.542	1.641	−1.572
	CW	−0.130	1.578	−1.542	1.641	−1.572
	HA	0.210
	C1 (methylene)	−0.170	1.788	−1.641	1.860	−1.674
	C1 (methyl)	−0.170	1.998	−1.735	2.078	−1.769
	H1	0.130
	CE	−0.050	1.998	−1.735	2.078	−1.769
	HC	0.060

---

As the hydrogen cannot be made polarizable, their polarizability is added to the heavy atom to which they are attached to. Consequently, the polarizability for the atom with the atom type C1 depends on the number of bonded hydrogen. The terminal C1 has three hydrogen and consequently a slightly higher polarizability than the C1 carbon of the CH₂-group.

The non-polarizable scaled charge simulations were performed using a constant charge scaling factor of 0.77. As the command “SCALAR CHARGE MULT 0.77 SELE ALL END” causes dubious simulation artifacts in CHARMM, the partial charges were already scaled in the topology file.

1.3 Computational setup

All molecular dynamics simulations were performed using the program CHARMM,⁶¹ version 38b1, which we modified to include the charge transfer algorithm. The modified routines will be shared on request with other CHARMM users. The Lennard-Jones interactions smoothly were truncated between a distance of 10 to 11 Å and the particle mesh Ewald (PME) method was used for long-ranged electrostatic interactions with a real-space cut-off of 14 Å and a k-space damping constant of 0.41 Å. All bonds were constrained with SHAKE.

All simulations were run with a time step of 0.5 fs at a temperature of 300 K, using a Nosé–Hoover thermostat, with a relaxation time equal to 0.1 ps. In case of the polarizable force fields, the Drude positions were thermostated at a temperature of 1 K with a relaxation time equal to 5 fs. Simulations for parameter optimization were run for at least 1 ns.

The initial configuration of 500 ion pairs in a periodic box of an initial length of 52 Å was generated by PACKMOL⁶² and equilibrated after a short steepest descent energy minimization for 500 ps at constant pressure of 1 atm using an Andersen barostat with a relaxation time equal to 0.2 ps yielding the densities in Table 2. Dynamical properties were generated in the canonical, NVT ensemble during production runs of 4 ns and 8 ns for polarizable CT and the pure polarizable as well as the non-polarizable scaled charge simulations, respectively. Structures for ensemble averages and time correlation functions were saved every 200 steps.

Table 2 Liquid density, ρ, heat of vaporization, ΔH_vap, average dipole moments 〈μ^⊕〉 and 〈μ^⊖〉, as well as diffusion constants, D^⊕ and D^⊖, at a temperature of 300 K

	Pure pol.	Pol. CT	Scaled charges	Experiment
a These simulations were performed at constant volume with the density of the pure polarizable simulations to avoid additional density effects. b The dipole distributions in Fig. 2 were fitted by a Gaussian and its standard deviation. Selected properties of all models are given on Table 2.
ρ [g cm^−3]	1.232	1.279	1.232^a	1.280^63
ΔHvap [kcal mol^−1]	44.7	43.8	32.2	35.6^64
〈μ^⊕〉^b [Debye]	2.4 ± 0.7	2.3 ± 0.7	2.1 ± 0.2
〈μ^⊖〉^b [Debye]	0.4 ± 0.2	0.4 ± 0.2	0.3 ± 0.1
D ^⊕ [10^−8 cm^2 s^−1]	24	28	26	53^65
D ^⊖ [10^−8 cm^2 s^−1]	19	22	18	44^65
σ(0) [S m^−1]	0.43	0.55	1.1	1.55^66

---

2 Results

Selected properties of all models are given on Table 2. By parameterization, the density of the charge transfer model is in good agreement with experiment.⁶³ As the pure polarizable model⁶⁰ was not optimized in this respect, a small discrepancy of 3% was detected. For the heat of vaporization, the agreement is not as good but reasonable. Please note that the heat of vaporization is difficult to measure experimentally, since the vapor pressure is so low, and different methods, including surface tension, Knudsen, and temperature programmed desorption, can lead to values that differ by 10 kcal mol⁻¹ or more.^67,68 Simulations find a range of values, from 27.8 kcal mol⁻¹,⁶⁹ 29.4 kcal mol⁻¹,¹² 32.3 kcal mol⁻¹,⁹ to 43.5 kcal mol⁻¹.⁷⁰

2.1 Electrostatic properties

The model was parameterized to have similar dipoles to the pure polarizable model. The distributions of dipole moments for the ionic liquid as given by the charge transfer and the pure polarizable models are shown in Fig. 2 and look very similar. The distributions of molecular charges |q_i| is shown in Fig. 3. All cations and anions have an absolute charge value less than one with a standard deviation of about 0.08 e. Interestingly, the cationic distribution looks like a Gaussian. In case of the anions, the distribution is skewed allowing for very low molecular anionic charges of q_i < −0.7 e. The average charge was parameterized to be ±0.77 e, in agreement with the DFT-derived charges of Mondal and Balasubramanian.²¹ The distribution of charges q_i(t), resulting from differing local geometries, is one obvious difference from the non-polarizable, scaled charge models, which would have the same charge q_i at all times.

Fig. 3 The distribution of ion charges q_i for the charge transfer model. The vertical dashed line corresponds to the average charge 〈q_C2mim〉 = − 〈q_BF4〉 = 0.77 e.

The charges can deviate from their average values for relatively long periods of time, reflecting the long time scales to change an ion's local environment (see Fig. 4a and b). There are short time-scale fluctuations, on the order of 0.25 ps, and a longer time relaxation of the charges. The time scale for this relaxation can be determined by calculating the charge–charge autocorrelation function,

where q_i(t) is the charge of the molecular ion i at time t. From Fig. 4c, the fast and the longer exponential relaxation times can be seen with exponential decay constants of 50 ps for C₂mim⁺ and of 110 ps for BF₄⁻.

Fig. 4 Typical time series for the molecular charge q_i(t) of the cations (a) and the anions (b). The dashed line corresponds to the average of Fig. 3. The respective charge-charge autocorrelation function for the cations (blue) and the anions (orange) is plotted in (c).

2.2 Structure of the ionic liquid

The radial distribution functions g(r) between the centers-of-mass of the ions in Fig. 5 using the polarizable charge transfer model (orange dashed line), the non-polarizable scaled charge model (black dotted line) and the pure polarizable model (gray solid line, ref. 60) are fairly similar. In case of the radial distribution functions of ions of like charge, i.e. g₊₊(r) and g₋₋(r), the pure polarizable model shows minor differences with the CT model in the first solvation shell probably due to π–π stacking of the imidazolium rings. For a more detailed analysis, we considered imidazoliums with a center-of-mass/center-of-mass distance of less than 5 Å and to have significant π–π-interactions. The angle θ(t) was computed via vectors which are perpendicular to the imidazolium rings respectively. A second Legendre polynomial of cos θ(t) larger than 0.90 corresponds to angles 0° ≤ θ(t) ≤ 15° (above the reference imidazolium ring) or 165° ≤ θ(t) ≤ 180° (below the reference imidazolium ring). The numbers in Fig. 5b represent the average number of parallel imidazolium pairs per time frame and the total number of imidazolium pairs per time frame which are 5 Å apart. As already visible in g₊₊(r) the total number in case of the pure polarizable force field is higher. Overall roughly 20% of the imidazolium pairs fulfill our criterion on parallel rings in all force field models.

Fig. 5 Radial distribution functions g(r) between the ion center-of-masses comparing pure polarizable model (solid line), the polarizable CT (dashed line) and the scaled charge model (dotted line).

The increased number of pairs in the pure polarizable model may be due to the fact that the repulsion between ions of the same species is realized by an inner solvent effect.^18,19,22 Here, the induced dipoles act as an additional solvent to screen the partial charges and thereby reducing the repulsion. However, this dielectric continuum theory works for larger distances which may be violated by ions which are direct neighbors. In case of the polarizable CT and the non-polarizable scaled charge model the reduced repulsion is simply governed by the sub-integer charges. In principle, the π–π stacking is not parametrized in the force field. As the partial charges of the reference imidazolium ring atoms are higher on average in the pure polarizable system than for the other two force fields, the induced dipoles of nearby imidazolium ring atoms are more attracted in this case which may enhance the overall number of pairs and hence the number of π–π stacking pairs.

2.3 Transport properties of the ionic liquid

The diffusion constants in Table 2 are calculated from the slope of the mean-square displacement of the center-of-mass of each ion 〈Δr_i(t)²〉 after 900 ps which should correspond to the diffusive regime. A system size correction^71,72 is made, which adds a term to the diffusion constant equal to 2.837 k_BT/(6πη L), where k_B is Boltzmann's constant, L is the simulation box length, and η is the viscosity. We use the experimental value of η equal to 29.7 mPa s.⁶⁵ With an average box length equal to 50.5 Å, this gives a correction equal to 4.2 × 10⁻⁸ cm² s⁻¹. The computational diffusion coefficients are about a factor of two too small compared to the experimental values.⁶⁵ This may be due to the fact that the Lennard-Jones parameter of the polarizable force field have not been changed when adding polarizability. This results in a double counting of dispersion¹⁹ as the dispersion was already parametrized for the non-polarizable force field and the interaction of the induced dipoles further increases the dispersion interactions. A scaling scheme for the Lennard-Jones well depth has been used for polarizable simulations of ionic liquids, which lead to increased the conductivity (and hence the overall dynamics) by a factor of almost four in case of C₂mim triflate and C₂mim dicyanamide.⁷³ However, in order to compare the non-polarizable scaled charge model, the polarizable CT and the pure polarizable simulation in this work, we use the same Lennard-Jones parameters for all simulations so that all differences are due to electrostatic interactions. Non-polarizable full-charge models are about an order of magnitude too small in dynamics.^{7,10,16,18,19,74}

As visible in Table 2 we find for all models in agreement with experiment that the C₂mim⁺ ions diffuse faster, despite being heavier. Reliable estimates of the error bars for the diffusion constants would require longer simulations than carried out here.⁷⁵ However, the pure polarizable model gives diffusion constants that are slightly smaller than the polarizable CT and non-polarizable scaled charge model but within the range of the statistical uncertainty. Additional information about translational motion is given by the velocity–velocity autocorrelation functions, as defined by

where v_i(t) is the center-of-mass velocity of ion i and t is time. Both ions show a mean collision time near 0.2 ps. Since the charges have a much slower relaxation time than the mean collision time, the velocity autocorrelation functions can be calculated separately for ions with bigger charges from those with smaller charges. Ions are considered to have a bigger charge if the magnitude of the charge is greater than 0.81 e at time zero and to have a smaller charge if the charge magnitude is less than 0.74 e at time zero. These charge values are about one standard deviation away from the mean (Fig. 3).

There is not a big difference between various models as visible in Fig. 6a and b. From the short-time dependence of C_vv(t), the mean-squared-forces 〈F²〉 can be found according to

where m_i is the mass of the ion.^76,77 Fitting eqn (10) to the first 0.05 ps of the velocity autocorrelation functions gives the results shown in Table 3.

Fig. 6 Velocity–velocity autocorrelation functions for: (a) C₂mim⁺ (b) BF₄⁻. (c) high charged (dark orange line, q_i ≥ 0.81 e) and low (light orange line, q_i ≤ 0.74 e) charged polarizable CT C₂mim⁺. (d) High charged (dark orange line, q_i ≤ −0.81 e) and low (light orange line, q_i ≥ −0.74 e) charged polarizable CT BF₄⁻.

Table 3 The mean-squared forces on the ions and diffusion coefficients of the ions in the different models

	|qi| [e]	〈F^2〉 [(10^−6 N)^2]		D [10^−8 cm^2 s^−1]
		C2mim^+	BF4^−	C2mim^+	BF4^−
Pure polarizable	= 1.00	0.50	0.36	24	19
Polarizable CT
All ions	<1.00	0.53	0.39	28	22
High charge	>0.81	0.49	0.32
Low charge	<0.74	0.62	0.44
Scaled charges	= 0.77	0.47	0.33	26	18

---

As expected, the 〈F²〉 for the pure polarizable model and the polarizable CT model with absolute molecular charges close to 1 e reasonably agree with each other. Also the forces acting on the scaled charge ions are similar to the pure polarizable model which means that the inner solvent screening of the induced dipoles can be reproduced by scaled charges. However, 〈F²〉 acting on polarizable CT ions with low absolute molecular charge differ significantly from those acting on the scaled charge ions. This underlines the fact that non-polarizable scaled charge models are not appropriate to model charge transfer.

For the CT model, the ions with larger charges diffuse faster, with less backscattering, and experience smaller mean forces (Table 3) than those with smaller charges. This is counter-intuitive to what might be expected. Ions with decreased charges should have smaller forces and diffuse faster, as is found in the scaled charge models. The high charge ions diffuse faster due to the mechanism of charge transfer. These ions have larger charges (so closer to ±1 e) because they are transferring less charge to neighboring, oppositely-charged ions. Therefore, the ions with larger charges are in a less confined environment and those with smaller charges (in magnitude) are in a more cage-like geometry. This is consistent with a strong correlation between local environments with high number density and slow diffusion, as described recently.⁷⁸

According to eqn (10) stronger mean-square forces on the molecules lead to faster decay of the corresponding velocity autocorrelation functions. One interpretation involves bounce effects of a strong surrounding ion cage.⁷⁶ Increased interactions with the ion cage – a larger 〈F²〉 – would lead to increased librational motion and hindered translations of the ion in its cage. Consequently, the velocity autocorrelation function relaxes very fast. However, this tight ion cage will hamper the diffusion of the central ion and hence a decrease in diffusion coefficients should be observed. This is the case for polarizable simulations of Lithium doped 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide.⁷⁶ However, all our simulations find that the C₂mim⁺ ions feel stronger forces than BF₄⁻ despite the fact that the cations diffuse faster. Interestingly, the ratio of the mean-square-forces of cations and anions is similar to the ratio of diffusion coefficients.

Another way to interpret the decay of 〈v_i(0)·v_i(t)〉 is a loss of memory effect. Memory can be of a collective nature as ionic liquids are known to be organized in spatial and dynamically heterogeneous domains. Applying hole theory to describe the conductivity in ionic liquids,^79,80 the emergence of a suitable hole for the jumping ion is not only a function of the ion but also of the collective motion of the surrounding ions. From the point of view of the jumping ion, the drag forces towards the hole are random and consequently the memory of its former velocity decays rapidly. In this new local environment, new partners for charge transfer exist driving the ion in a new direction.

Another mechanism may also explain the coexistence of high mean-square forces and high mobility of the ions.^75,81 In a coarse-grained simulation of the mixture of [C₂mim][BF₄] and [C₂mim] triflate, the mobile ions were not randomly distributed, but associated into clusters of mobile particles, resulting in much higher radial distribution functions for like-charged ions for example. Here, the electrostatic repulsion of the low charge ions in Table 3 is weaker resulting in a stronger force keeping this clusters together.

2.4 Conductivity

Information about collective motion of the ions can be given by the charge current and total charge dynamics.^82,83 The electric current J(t) is given bywith the sum over all ions of the system. The molecular charges are only functions of time for the polarizable CT simulations. In case of the pure polarizable as well as the scaled charge simulations, the molecular charge is constant, i.e. q_i(t) = q_i = ±1 e or ±0.77 e. In principle, the electric current autocorrelation functioncan be used to compute the electrical conductivity σ(0).^83,84 At first sight, C_JJ(t) relaxes to zero within a few hundred femtoseconds (see Fig. 7a). However, the integration of C_JJ(t) has to be performed for several dozens of picoseconds to yield reliable conductivities as long-time processes keep C_JJ(t) slightly below zero.⁸⁵

Fig. 7 (a) Electric current time autocorrelation function. (b) Mean-squared displacement of the collective translational dipole moment.

Interestingly, the negative “rebound” regime in Fig. 7a of the collective C_JJ(t) is much more pronounced than for the molecular C_vv(t) with a deeper first minima, around 0.16 ps.

The static conductivity σ(0) can be obtained from the mean-square displacement of the collective translational dipole,⁸³ which is statistically more reliable compared to the integration of the current-current correlation function. The collective translational dipole of a particular species k is defined by

where r_i(t) is the time-dependent position of the center-of-mass of ion i of species k. The overall collective translational dipole moment M_J(t) is the sum of all components belonging to a particular species:⁷⁵

At times beyond the correlation time of the C_JJ(t), approximately 100 ps in this work, the mean-square displacement of M_J(t) becomes linear⁸³

〈ΔM_J²(t)〉 = 〈(M_J(0) − M_J(t))²〉 (15)

∝ 6 V k_BTσ(0)·t. (16)

and its slope correlates with the static conductivity σ(0). If one restricts the collective translational dipole moment M^k_J(t) to the cations (k = ⊕) or anions (k = ⊖), their contribution to the overall conductivity can be computed:⁷⁵

〈ΔM_J²(t)〉 = 〈ΔM^⊕_J(t)·ΔM_J(t)〉 + 〈ΔM^⊖_J(t)·ΔM_J(t)〉 (17)

〈ΔM^⊕_J(t)·ΔM_J(t)〉 ∝ 6 V k_BTσ^⊕(0)·t (18)

〈ΔM^⊖_J(t)·ΔM_J(t)〉 ∝ 6 V k_BTσ^⊖(0)·t (19)

Both summands in eqn (17) show the typical behavior of a mean-square displacement. Consequently, their slopes can be attributed to a conductivity contribution of that ion type. These values are also given is Table 4. Within the statistical uncertainty σ(0)-values of all our computational models are around 0.5 S m⁻¹ which is significantly lower than the experimental value of 1.55 S m⁻¹.⁶⁶ The discrepancy is probably due to the double counting of dispersion effects in our MD simulations, as discussed above in connection to the diffusion constants results.

Table 4 Contributions to the static conductivity σ(0). Please note that σ(0) = σ^⊕(0) + σ^⊖(0) (see eqn (17))

	σ(0) [S m^−1]	σ ^⊕(0) [S m^−1]	σ ^⊖(0) [S m^−1]	σ NE [S m^−1]
Pure pol.	0.43	0.19	0.24	1.00
Pol. CT	0.55	0.06	0.49	1.20
Scaled charges	0.36	0.16	0.20	1.06

---

So far, we have seen an amazing coincidence of the structural and dynamical results of the pure polarizable, the polarizable CT and the scaled charge model. The even better agreement between the pure polarizable and the scaled charge model points out that: First, one should use full charges for analysis as the scaling only mimics an averaged polarizability effect and not charge transfer. Second, as long as one is interested in mesoscopic properties like the static conductivity, charge scaled models are capable of reproducing results from polarizable simulations as noted previously.¹⁸ However, we expected local electrostatic interactions to be different in polarizable and charge scaled simulations since the dielectric continuum theory^18,22 is not valid anymore at short distances between the interacting charges.

Does charge transfer change the dynamics? At first glance, the answer is no. Pure polarizable simulations seem to be capable of reproducing the structure and dynamics of the CT simulations, although the molecular charges differ significantly, i.e. q^pol_i = ±1 e and 〈q^CT_i(t)〉 = ±0.77 e. The different response of the induced dipoles in both cases leads in the end to similar results.

However, at second glance, there is a difference between polarizable charge transfer and pure polarizable simulations as visible in Table 4. The partial conductivities σ^⊕(0) and σ^⊖(0) are noticeably different among the models. The difference in the cation and anion conductivities are much more pronounced for the CT model, for which σ^⊕(0) is much smaller than σ^⊖(0).

Comparing the individual contributions to the overall conductivity σ(0), the pure polarizable and the scaled charge simulations predict only a minor enhancement for the anions. Similar results are found for tetrafluoroborate and dicyanamide.⁷⁵ However, in the polarizable CT simulations the anions become the dominant species for charge transport. Here, the cations play only a minor role. Since the Nernst–Einstein conductivity σ_NE

is almost twice the value of the collective static conductivity σ(0), long-living (non-)neutral ion aggregates are expected. These aggregates are a prerequisite for charge transfer. Even more, reducing the cationic molecular charge shifts the complex balance between repulsive Coulomb and attractive Lennard-Jones interactions in favor of the latter. As a result, aggregates consisting of two or more cations may stay longer together. Compared to the larger cations with more atoms, the Lennard-Jones interactions of the anions are expected to be weaker. Hence, aggregates containing several anions may not be as long-living as the aggregates dominated by the cations. Consequently, the charge transport of the anions is more efficient.

3 Conclusions

Our pure polarizable, charge scaled and polarizable charge transfer MD simulations resulted in very similar structure and dynamics of the ionic liquid 1-ethyl-3-methylimidazolium tetrafluoroborate. The agreement between the pure polarizable and the charge scaled simulations using full charges for analysis demonstrates that the charge scaling mimics an average polarizability rather than charge transfer in ionic liquids. Charge transfer effects in ionic liquids seem to be very subtle. Given that charge transfer amounts depend on the ion pairs,^21,28,37 these effects may be different for other ILs. Charge transfer seems to depend more on the identity of the anion than the cation,^21,28 and is larger for monatomic ions like chloride than for polyatomic ions,^21,37 so these effects could be greater for other ILs. One major discrepancy between the pure polarizable and the polarizable charge transfer simulations was detected for the partial conductivities. Ions can have charges that deviate from average values for about 100 ps (Fig. 4). On this time scale, ions which have transferred more charge, and have an overall charge with a smaller magnitude, diffuse slower than those with a larger charge.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Science Foundation under contract number CHE-0611679.

Notes and references

1. P. Wasserscheid, Nature, 2006, 439, 797.
2. R. D. Rogers, Nature, 2007, 447, 917–918.
3. C. J. Margulis, H. A. Stern and B. J. Berne, J. Phys. Chem. B, 2002, 106, 12017–12021.
4. T. I. Morrow and E. J. Maginn, J. Phys. Chem. B, 2002, 106, 12807–12813.
5. J. N. Canongia Lopes, J. Deschamps and A. A. H. Pádua, J. Phys. Chem. B, 2004, 108, 2038–2047.
6. M. G. Del Pópolo and G. A. Voth, J. Phys. Chem. B, 2004, 108, 1744–1752.
7. T. Yan, C. J. Burnham, M. G. Del Pópolo and G. A. Voth, J. Phys. Chem. B, 2004, 108, 11877–11881.
8. T. G. A. Youngs and C. Hardacre, ChemPhysChem, 2008, 9, 1548–1558.
9. O. Borodin, J. Phys. Chem. B, 2009, 113, 11463–11478.
10. E. J. Maginn, J. Phys.: Condens. Matter, 2009, 21, 373101.
11. C. Schröder and O. Steinhauser, J. Chem. Phys., 2010, 133, 154511.
12. V. V. Chaban, I. V. Voroshylova and O. N. Kalugin, Phys. Chem. Chem. Phys., 2011, 13, 7910.
13. A. O. Cavalcante, M. C. C. Ribeiro and M. S. Skaf, J. Chem. Phys., 2014, 140, 144108.
14. Y. Wu, Y. Li, N. Hu and M. Hong, Phys. Chem. Chem. Phys., 2014, 16, 2674.
15. M. Campetella, L. Gontrani, F. Leonelli, L. Bencivenni and R. Caminiti, Chem. Phys. Chem., 2014, 16, 197–203.
16. M. Salanne, Phys. Chem. Chem. Phys., 2015, 17, 14270–14279.
17. O. Borodin and G. D. Smith, J. Phys. Chem., 2006, 110, 11481–11490.
18. C. Schröder, Phys. Chem. Chem. Phys., 2012, 14, 3089.
19. D. Bedrov, J.-P. Piquemal, O. Borodin, A. D. MacKerell Jr., B. Roux and C. Schröder, Chem. Rev., 2019, 119, 7940–7995.
20. M. Kohagen, M. Brehm, Y. Lingscheid, R. Giernoth, J. Sangoro, F. Kremer, S. Naumov, C. Iacob, J. Kärger, R. Valiullin and B. Kirchner, J. Phys. Chem. B, 2011, 115, 15280–15288.
21. A. Mondal and S. Balasubramanian, J. Phys. Chem. B, 2014, 118, 3409–3422.
22. I. Leontyev and A. Stuchebrukhov, Phys. Chem. Chem. Phys., 2011, 13, 2613–2626.
23. L. Pegado, O. Marsalek, P. Jungwirth and E. Wernersson, Phys. Chem. Chem. Phys., 2012, 14, 10248–10257.
24. Z. R. Kann and J. L. Skinner, J. Chem. Phys., 2014, 141, 104507.
25. J. Li and F. Wang, J. Chem. Phys., 2015, 143, 194505.
26. M. Bühl, A. Chaumont, R. Schurhammer and G. Wipff, J. Phys. Chem. B, 2005, 109, 18591–18599.
27. J. Schmidt, C. Krekeler, F. Dommert, Y. Zhao, R. Berger, L. Delle Site and C. Holm, J. Phys. Chem. B, 2010, 114, 6150–6155.
28. T. Cremer, C. Kolbeck, K. R. J. Lovelock, N. Paape, R. Wölfel, P. S. Schulz, P. S. Wasserscheid, H. Weber, J. Thar, B. Kirchner, F. Maier and H. P. Steinrück, Chem. – Eur. J., 2010, 16, 9018–9033.
29. M. Dal Peraro, S. Raugei, P. Carloni and M. L. Klein, Chem. Phys. Chem., 2005, 6, 1715–1718.
30. Z. Zhao, D. M. Rogers and T. L. Beck, J. Chem. Phys., 2010, 132, 014502.
31. S. Varma and S. B. Rempe, Biophys. J., 2010, 99, 3394–3401.
32. M. Soniat and S. W. Rick, J. Chem. Phys., 2012, 137, 044511.
33. B. Sellner, M. Valiev and S. M. Kathmann, J. Phys. Chem. B, 2013, 117, 10869–10882.
34. Y. Yao, Y. Kanai and M. L. Berkowitz, J. Phys. Chem. Lett., 2014, 5, 2711–2716.
35. Y. Yao, M. L. Berkowitz and Y. Kanai, J. Chem. Phys., 2015, 143, 241101.
36. M. Nguyen and S. W. Rick, J. Chem. Phys., 2018, 148, 222803.
37. O. Hollóczki, F. Malberg, T. Welton and B. Kirchner, Phys. Chem. Chem. Phys., 2014, 16, 16880.
38. T. Brinzer, E. J. Berquist, Z. Ren, S. Dutta, C. A. Johnson, C. S. Krishner, D. S. Lambrecht and S. Garrett-Roe, J. Chem. Phys., 2015, 142, 212425.
39. P. Ayotte, C. G. Bailey, G. H. Weddle and M. A. Johnson, J. Phys. Chem. A, 1998, 102, 3067–3071.
40. W. H. Thompson and J. T. Hynes, J. Am. Chem. Soc., 2000, 122, 6278–6286.
41. S. G. Ramesh, S. Re and J. T. Hynes, J. Phys. Chem. A, 2008, 112, 3391–3398.
42. D. E. Wu, S. Duan, X. M. Liu, Y. C. Xu, Y. X. Jiang, B. Ren, X. Xu, S. H. Lin and Z. Q. Tian, J. Phys. Chem. A, 2008, 112, 1313–1321.
43. Q. Wan, L. Spanu, G. A. Galli and F. Gygi, J. Chem. Theory Comput., 2013, 9, 4124–4130.
44. D. Sidler, M. Meuwly and P. Hamm, J. Chem. Phys., 2018, 148, 244504.
45. C. D. Cappa, J. D. Smith, B. M. Messer, R. C. Cohen and R. J. Saykally, J. Phys. Chem. B, 2006, 110, 5301–5309.
46. A. J. Lee and S. W. Rick, J. Chem. Phys., 2011, 134, 184507.
47. C. D. Wick, A. J. Lee and S. W. Rick, J. Chem. Phys., 2012, 137, 154701.
48. S. W. Rick, J. Comput. Chem., 2016, 37, 2060–2066.
49. M. Soniat and S. W. Rick, J. Chem. Phys., 2014, 140, 184703.
50. M. Soniat, L. Hartman and S. W. Rick, J. Chem. Theory Comput., 2015, 11, 1658–1667.
51. M. Soniat, R. Kumar and S. W. Rick, J. Chem. Phys., 2015, 143, 044702.
52. M. Soniat, G. Pool, L. Franklin and S. W. Rick, Fluid Phase Equilib., 2015, 407, 31–38.
53. J. P. Piquemal, A. Marquez, O. Parisel and C. Giessner-Prettre, J. Comput. Chem., 2005, 26, 1052–1062.
54. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
55. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. 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, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09 Revision A.1, Gaussian Inc., Wallingford CT, 2009.
56. R. F. W. Bader, Acc. Chem. Res., 1985, 18, 9–15.
57. T. G. S. Todd and A. Keith, AIMAll (Version 14.06.21), Overland Park KS, USA, (aim.tkgristmill.com), 2014.
58. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361.
59. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2012, 8, 2844–2867.
60. M. Schmollngruber, C. Schröder and O. Steinhauser, J. Chem. Phys., 2013, 138, 204504.
61. B. R. Brooks, C. L. Brooks III, A. D. MacKerell Jr., L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York and M. Karplus, J. Comput. Chem., 2009, 30, 1545.
62. L. Martínez, R. Andrade, G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157.
63. M. B. Shiflett and A. Yokozeki, J. Chem. Eng. Data, 2007, 52, 1302–1306.
64. S. P. Verevkin, Angew. Chem., Int. Ed., 2008, 47, 5071–5074.
65. A. Noda, K. Hayamizu and M. Watanabe, J. Phys. Chem. B, 2001, 105, 4603–4610.
66. A. Stoppa, O. Zech, W. Kunz and R. Buchner, J. Chem. Eng. Data, 2010, 55, 1768–1773.
67. J. P. Armstrong, C. Hurst, R. G. Jones, P. Licence, K. R. J. Lovelock, C. J. Satterley and I. J. Villar-Garcia, Phys. Chem. Chem. Phys., 2007, 9, 982–990.
68. S. P. Verevkin, D. H. Zaitsau, V. N. Emelyanenko, A. V. Yermalayeu, C. Schick, H. Liu, E. J. Maginn, S. Bulut, I. Krossing and R. Kalb, J. Phys. Chem. B, 2013, 117, 6473–6486.
69. S. V. Sambasivarao and O. Acevedo, J. Chem. Theory Comput., 2009, 5, 1038–1050.
70. G. Raabe and J. Köhler, J. Chem. Phys., 2008, 128, 154509.
71. B. Dünweg and K. Kremer, J. Chem. Phys., 1993, 99, 6983–6997.
72. I.-C. Yeh and G. Hummer, J. Phys. Chem. B, 2004, 108, 15873–15879.
73. E. Heid, B. Docampo-Àlvarez, L. M. Varela, K. Prosenz, O. Steinhauser and C. Schröder, Phys. Chem. Chem. Phys., 2018, 20, 15106–15117.
74. S. Tsuzuki, W. Shinoda, H. Saito, M. Mikami, H. Tokuda and M. Watanabe, J. Phys. Chem. B, 2009, 113, 10641–10649.
75. V. Zeindlhofer, L. Zehetner, W. Paschinger, A. Bismarck and C. Schröder, J. Mol. Liq., 2019, 288, 110993.
76. V. Lesch, H. Montes-Campos, T. Méndez-Molales, L. J. Gallego, A. Heuer, C. Schröder and L. M. Varela, J. Chem. Phys., 2016, 145, 204507.
77. T. Katö, K. Machida, M. Oobatake and S. Hayashi, J. Chem. Phys., 1988, 89, 3211–3221.
78. R. P. Daly, J. C. Araque and C. J. Margulis, J. Chem. Phys., 2017, 147, 061102.
79. A. P. Abbott, ChemPhysChem, 2005, 6, 2502–2505.
80. H. Zhao, Z.-C. Liang and F. Li, J. Mol. Liq., 2009, 149, 55–59.
81. W. Kob, C. Donati, S. J. Plimpton, P. H. Poole and S. C. Glotzer, Phys. Rev. Lett., 1997, 79, 2827–2830.
82. M. Salanne, C. Simon, P. Turq and P. A. Madden, J. Phys. Chem. B, 2007, 111, 4678–4684.
83. C. Schröder, M. Haberler and O. Steinhauser, J. Chem. Phys., 2008, 128, 134501.
84. B. L. Bhargava and S. Balasubramanian, J. Chem. Phys., 2005, 123, 144505.
85. C. Schröder, J. Chem. Phys., 2011, 135, 024502.

---