Polarizable molecular dynamics simulations of ionic liquids: Influence of temperature control

Ionic liquids are an interesting class of soft matter with viscosities of one or two orders of magnitude higher than water. Unfortunately, classical, non-polarizable MD simulations of ionic liquids result in too slow dynamics and demonstrate the need for explicit inclusion of polarizability. The inclusion of polarizability, here via the Drude oscillator model, requires amendments to the employed thermostat, where we consider a dual Nosé-Hoover thermostat, as well as a dual Langevin thermostat. We investigate the effects of the choice of a thermostat, as well as the underlying parameters such as the masses and force constants of the Drude particles on static and dynamic properties of ionic liquids. Here, we show that Langevin thermostats are not suitable for the investigation of the dynamics of ionic liquids. Since polarizable MD simulations are associated with high computational costs, we employed a self-developed GPU enhanced code within the MD program CHARMM to keep the overall computational effort reasonable.


I Introduction
Aqueous concentrated salt solutions and ionic liquids, i.e. salts which are liquid at ambient temperatures, have numerous applications in chemistry, physics and energy storage/ conversion. 1,2 Often molecular dynamics (MD) simulations are employed to gain fundamental insights into their structural and transport properties. The (molecular) ions generate strong electric fields polarizing solvent molecules as well as neighboring ions. Consequently, the prediction of physico-chemical properties of such systems employing MD simulations often requires explicit inclusion of polarizable forces to damp the Coulomb interactions between the ions. 3 Although non-polarizable MD simulations are able to reproduce structural experimental data, they usually result in dynamics which are almost one order of magnitude too slow compared to experiment. [4][5][6] In contrast, polarizable ionic liquid simulations were shown to produce dielectric spectra, 7 conductivities, 4,7 diffusion coefficients, 4 viscosities 4 and solvation dynamics 8,9 in close agreement to experiment. However, despite their high accuracy, the increased computational effort of polarizable simulations often limits their use. Due to the very small self-diffusion constants of ionic liquids, realistic diffusive behavior can only be observed by simulating several dozens to hundreds of nanoseconds. 10 The need for long simulation times, together with the slow implementation of polarizability in common MD packages, renders polarizable simulations of ionic liquids difficult and time-consuming.
In general, different approaches to account for explicit inclusion of polarization exist, with the most popular ones being the classical Drude oscillator model, 11,12 the induced dipole model 13,14 and the fluctuating charge model. 15 In this study, we focus on the Drude model, where a charged particle is attached via a harmonic spring to each nucleus, mimicking the electron cloud. As explained in detail in section 2, the determination of the positions of the Drude particles is timeconsuming, since the displacement vector of each Drude particle depends on all other Drude particles, and an exact solution can thus only be achieved in a self-consistent field (SCF) approach.
A convenient alternative to the SCF procedure is the so-called extended Lagrangian method, 12,16,17 where the induced electronic polarization at each timestep is obtained via Newton's equations of motion, i.e. by assigning the Drude particles a small mass and treating them as any other physical particle, instead of optimizing the positions of the Drude particles via an SCF procedure. However, the Drude particles need special treatment concerning the thermostat (and thus for NVT and NPT ensembles), since the temperature of the Drude degrees of freedom needs to be low (few degrees above absolute zero), whereas other atoms are kept at the desired temperature, e.g. room temperature. Thus, the extended Lagrangian formalism requires the use of a dual thermostat, which narrows down the possible choices of a thermostat considerably. To the best of our knowledge there are only two implementations of dual thermostats currently available, which are summarized in the following. A dual Nosé-Hoover thermostat 18,19 was implemented in the program package CHARMM, 12,20 to enable polarizable simulations based on the Drude oscillator model. The Nosé-Hoover thermostat is known to yield accurate structural and dynamic properties of liquids. 12 Furthermore, a dual-Langevin thermostat 21 was implemented in NAMD and is also available in CHARMM. A drawback of the Langevin thermostat is the dependency of the dynamic properties of a system on the chosen friction coefficients of the dual heat baths, 22,23 so that often the Langevin thermostat is only used to generate structural properties. Very low friction coefficients were used by Jiang et al. to calculate dynamic properties of a salt solution; 21 however, it is unknown whether this method translates well to ionic liquids.
Both implementations of dual thermostats suffer from high computational costs. Fast GPUbased variants such as DOMDEC CHARMM 24 are available for nonpolarizable simulations. Also, several other program codes 25,26 offer GPU enhancement for polarizable simulations with various thermostats and barostats. However, there are currently no options to speed up polarizable simulations in CHARMM. Therefore, we developed a GPU enhanced polarizable CHARMM code. All thermostats implemented in CHARMM are available in the newly developed code, and we test both a dual Nosé-Hoover, as well as a dual Langevin thermostat.
The accuracy of properties obtained from MD simulations can depend on the employed temperature control parameters, such as the magnitude of the mass assigned to the Drude particles, the force constants of the harmonic Drude oscillators, or the friction coefficients of the Drude degrees of freedom in the case of the Langevin thermostat. In this study, we examine the effects of those parameters on the structural and dynamic properties of two ionic liquid systems, as well as compare the performance and obtained properties between the Langevin and the Nosé-Hoover dual thermostats.
The remainder of this article is organized as follows. In section 2 details of the Drude oscillator model and its implementations via self-consistent field, dual Langevin and dual Nosé-Hoover thermostat are given. In section 3, details of the GPU enhanced CHARMM code are presented, and the computational setup of this study is described. A detailed analysis of the timescales and properties of various simulations of the ionic liquids 1ethyl-3-methylimidazolium dicyanamide and 1-ethyl-3-methylimidazolium trifluoromethanesulfonate with and without GPU enhancement, different thermostats, and different simulation parameters is given in section 3. Concluding remarks are presented in section 4.

II Theory
In the Drude oscillator model a pair of Drude particles representing the induced dipole is attached to a polarizable atom β of molecule i with partial charge q iβ and an atomic polarizability α iβ . The mobile Drude particle has a partial charge of q iβ D , which is usually negative, and moves around the polarizable atom. It is attached to the nucleus via a harmonic spring. The harmonic potential is centered around a dislocation of zero, where k iβ D is the force constant and d iβ (t) the vector from the nucleus to the mobile Drude particle. We note that within the program package CHARMM 20 the respective parameter for specifying the force constant is KDRUDE, where KDRUDE = k iβ D /2. The distance d iβ (t) to the central atom depends on the local electric field, i.e. a strong electric field leads to a large dislocation and thus a large induced dipole via μ iβ ind (t) = q iβ D ⋅ d iβ (t) . (2) The second Drude particle is always located at the position of the polarizable atom and has a partial charge of −q iβ D .
The harmonic force constant k iβ D and the Drude charge q iβ D are correlated via where ε 0 is the vacuum permittivity. As a consequence, usually only the force constant k D is given in polarizable force fields with a standard value of 1000 or 2000kcal/mol/Å 2 for all atoms iβ. 3,12,27 The Drude charges then arise from k D and the atomic polarizability α iβ , where a larger force constant leads to larger Drude charges at smaller displacements, and thus the same induced dipole moment. In other words, the induced dipoles should in principle not depend on the chosen q D or k D but only on α iβ . In a previous study 28  polarizable atoms. We noticed that the average strength of the induced dipoles is slightly higher for q D =−0.5 e than for the other Drude charges. 3,28 Also its distribution is significantly broader in this case. In addition, the static dielectric constant was higher compared to the polarizable simulations using value of |q D | > 1 e. 28 Consequently, the Drude force constant of 2000kcal/mol/Å 2 is superior to 1000kcal/mol/Å 2 as it results in higher Drude charges which are closer to the polarizable atoms.

A Self-consistent field minimization
The induced dipole μ iβ ind (t) depends on the local electric field E iβ (t) depends on the local electric field E iβ (t) and the atomic polarizability as where the local electric field is with the distance r iβjγ (t) between the atoms β of molecule i and atoms γ of molecule j at time t. Since all the surrounding induced dipoles μ jγ ind (t) of atoms γ of molecule j enter the local field E iβ acting on atom β, the induced dipoles have to be determined in a selfconsistent field manner 3,12 to find the optimal displacement d iβ which originates from the balance of the electrostatic forces of the local electric field E iβ and the forces of the harmonic spring constant D iβ D between the nucleus and the mobile Drude particle: The traditional self-consistent field procedure is limited and computationally very expensive, 12 since many iterations may be necessary depending on the convergence criterion δ ≃ 0, i.e. ∂U/∂d iβ < δ. If the energy minimization is not completely converged, it will introduce systematic drag forces affecting the stability of the simulations as well as the physicochemical properties of the investigated particles. 3

B Extended Lagrangian thermostat
An alternative to the self-consistent field procedure is extending the Lagrangian of the system to treat k iβ D as additional dynamical degrees of freedom. 12,16,17 The mobile Drude particles are attributed a small mass m D which is subtracted from the mass m iβ of the polarizable atom. Thus the Drude particles can be treated according to Newtonian mechanics similar to the nuclei. However, this treatment requires the use of a dual thermostat, i.e. different thermostatting of the non-Drude degrees of freedom and the Drude degrees of freedom at different temperatures. In the following, the dual Nosé-Hoover thermostat and the dual Langevin thermostat are explained in detail.
1 Nosé-Hoover thermostat-Lamoureux and Roux 12 implemented a dual Nosé-Hoover thermostat 18,19 method in the molecular dynamics package CHARMM. 20 The first thermostat acts on each center-of-mass R iβ of the Drude-atom systems (8) which is close to the positions of the respective nucleus since the Drude mass m D and the displacement d iβ (t) of the mobile Drude particle are small. Increasing the Drude mass m D corresponds to a shift of the center-of-mass R iβ (t) of the Drude-atom system in direction of the mobile Drude particle.
The thermostat keeps the non-Drude atoms at the system temperature, usually around T =300 K via The second thermostat is coupled to the thermal fluctuations of the mobile Drude particle around the nucleus, i.e. ∂d i β (t)/∂t.
with the reduced mass m iβ D = m iβ − m D ⋅ m D /m iβ ⋅ 12 The forces F iβ (t) = − ∂U(t)/ ∂R iβ (t) and F iβ D (t) = − ∂U(t)/ ∂d iβ (t) act on the center-of-mass and the reduced mass, respectively. The parameter η and η D are friction coefficients of the thermostats Heid et al. Page 5 and are correlated with the number of atoms n and number of Drudes n D , the relaxation time of the thermostats, τ and τ D , as well as the respective thermostat temperatures T and T D .
The integration of Eq. 9 and 10 is usually split into multiple steps, namely, into a one-step "Velocity-Verlet" propagation where the atomic forces are evaluated once and a multistep "Nosé-Hoover" propagation, as explained in Ref. 12. A high number of intermediate steps (in CHARMM termed 'NSTEP') reduces systematic integration errors on η and η D , and thus increases the accuracy of the algorithm. 12 In order to be close to the self-consistent field procedure, the second thermostat should act on the vibration of the mobile Drude particle with the nucleus rather than the motion of the Drude particle directly. 12,21 The temperature, usually around 1 K, should be small enough to leave almost no kinetic energy to these vibrations but large enough to respond to the roomtemperature motions of the nuclei. Consequently, the mass m D and the force constant k iβ D are important parameters for the thermostatting of the Drude pairs. Also, a low m D requires a small timestep, and might thus affect the stability of a simulation. 12,[29][30][31] 2 Langevin thermostat-In the MD package NAMD 32 Jiang et al. 21 implemented a dual stochastic Langevin thermostat to improve the performance of polarizable MD simulations on parallel infrastructures. The dual Langevin thermostat is also available in CHARMM. The equations of motions read 14) and are decoupled in contrast to eq. (11 and (12 of the Nosé-Hoover thermostat, where in the differential equation for each atom the masses and coordinates of all other atoms (summation in the evaluation of the friction coefficients) need to be known. Within the Langevin formalism each equation can be treated separately and independently (no summation over all masses or coordinates), so that the Langevin thermostat is particularly fit for highly-parallel computing platforms. 21 f iβ (t) and f iβ D (t) are fluctuating random forces obeying the fluctuation-dissipation theorem where s(t) and s D (t) are univariate Gaussian random processes. The parameters γ and γ D are the frictions strengths on the center of mass of a Drude pair and the relative motion of the Drude pair, respectively. 21 Similar to the dual Nosé-Hoover thermostat, the dual Langevin thermostat keeps the Drude degrees of freedom at a different temperature T D , usually around 1 K, than the non-Drude atoms, which are kept at T, usually room temperature.

A GPU enhanced CHARMM
The most time consuming step in MD simulations of liquids and biomolecules is the calculation of the non-bonded interactions. 33 To keep the programming effort low, as well as to keep changes to the existing CHARMM code minimal, we decided to accelerate only the (real-space) non-bonded interactions on the GPU following the approach taken by Elber and coworkers. 33 Using NVIDIA's CUDA platform, CPU and GPU operate independently. This leads to the following, trivial strategy: while the GPU computes nonbonded interactions, the CPU takes care of the bonded interactions and particle mesh Ewald summation. 34,35 Integration of the equation of motion, application of constraints (e.g., SHAKE/RATTLE 36 ) are kept on the CPU. The normal parallel code of CHARMM is used for the CPU part; only a single GPU is supported. As a result, given a processor with n CPU cores plus GPU, the code will exploit all n CPU cores, as well as the GPU.
Given this strategy, only the calculation of Lennard-Jones and real-space Coulomb interactions needs to be ported to the GPU. Also, since for explicit solvent simulations only atom-based cutoffs are used, the group based cutoff methods existing in regular CHARMM were not supported.
The implementation is based on the work described in "Particle Simulation using CUDA" 37 and the corresponding example code from the CUDA toolkit. The primary changes to this reference code were threefold: (1) the simple repulsive potential of the example was replaced by the functional forms for Lennard-Jones and electrostatic interactions used in CHARMM, and provisions for different particle types (Lennard-Jones parameters, partial charges) were added. (2) The example code uses a spatial decomposition approach, in which one considers only interactions between particles within the same cell or in neighboring cells. Our implementation classifies interactions in three categories, i.e., the interaction of a particle with particles within the same cell, particles in the nearest neighbor cells, and particles in next to nearest neighbor cells. The spatial decomposition is carried out for each energy / force calculation, the cell-width is half the cut-off distance. While this scheme involves a bit more book-keeping, one does not have to worry about nonbonded exclusions during the calculation of particles with partners in next to nearest neighbor cells, the largest number of interactions, simplifying the code and the performance of these kernels. (3) The code takes into account non-bonded exclusions (cf. above) and special treatment of 1-4 interactions.
The performance of double precision floating point arithmetics in consumer grade graphics cards from NVIDIA is crippled. While the regular CHARMM code is all double precision, for performance reasons most calculations on the GPU have to be carried out in single precision. To avoid round-off errors, resulting in failure to conserve energy, care is needed during the accumulation of forces and energies. We used ideas from a paper by Levesque and Verlet, 38 in particular scaling of coordinates and forces by the box size. The accumulation of energies and forces is done in 32-bit fixed-point arithmetic, 39 a technique which was also adopted by AMBER. 40 Potential energy contributions are accumulated by particle; hence, even in this case, 32bit integers are enough. The total Lennard-Jones and (real-space) Coulomb interactions are obtained by a reduce operation at the very end; this is done in double precision and the performance cost is negligible. A similar strategy is used to obtain the components of the virial tensor. The CPU/GPU mixed-precision code described here conserves energy under NVE conditions almost as good as the regular, all double precision CHARMM code (loss in total energy is less than 0.5 kcal/mol over 100 ns, no noticeable temperature drift). We furthermore tested the differences in energies and forces between the conventional CPU code and the CPU/GPU mixed-precision code for the largest system in this study (1000 ion pairs, 37000 atoms). The difference in energy amounted to 2.0 × 10 −4 kcal/mol in absolute terms, and -2.2 × 10 −9 relative to the CPU energy. The largest deviation in forces on the atoms was found to be 8.3 × 10 −3 kcal/mol/Å (which corresponds to 5.8 × 10 −13 N). The mean absolute error averaged over all atoms amounted to 1.3 × 10 3 kcal/mol/Å (which corresponds to 9.0 × 10 −14 N).

B Computational setup
Polarizable Bonds involving hydrogens were kept at fixed length using the SHAKE/RATTLE algorithm. 36 We employed a timestep of 0.5 fs with different methods of temperature control, as specified in the following. Furthermore, the use of the TP control utility necessary for a dual Langevin or dual Nosé-Hoover thermostat in CHARMM requires the use of a Velocity-Verlet integrator.
Benchmark simulations with and without GPU enhancement were conducted for 500, 750 or 1000 ion pairs of [C 2 mim]DCA, where the ions were packed into cubic boxes using the program PACKMOL. 43 A Drude mass of m D =0.2 amu and a force constant of k D =1000 kcal/mol/Å 2 (corresponding to KDRUDE=500 in CHARMM) was employed. After an NPT equilibration period of 0.25 ns and an NVT equilibration period of 0.1ns, a 1ns NVT production run was calculated on a node equipped with an 8-core Xeon processor (Xeon W-2145) and an NVIDIA GTX1080 graphics processing unit using the following setup. For the systems containing 500, 750 and 1000 ion pairs, a CPU and GPU production run was calculated with a dual Nosé-Hoover thermostat. 18,19 The non-Drude particles were kept at 300 K, and the Drude particles at 1 K with response times of 0.1 ps and 0.005 ps respectively. In the following, we denote simulations with the Nosé-Hoover thermostat as 'NH'. For 1000 ion pairs, a CPU and GPU production run was produced additionally without Drude particles/polarizability. Furthermore, for 1000 ion pairs a GPU production run was calculated with a dual Langevin thermostat 21 with different friction coefficients. First, the friction strength on the center of mass of a Drude pair was set to 5.0 ps −1 , and the friction strength on the relative motion of the Drude pair to 10.0 ps −1 (default parameters in CHARMM documentation, similar to Ref. [44][45][46][47]. In the following, this setup is referred to as 'L1'. Second, the friction strengths were set to 1.0 ps −1 and 2.0 ps −1 respectively, denoted as 'L2', to estimate the influence of the chosen parameters. In biomolecular simulations, friction coefficients of 1-10 ps −1 are frequently used. 23 Within both setups, the temperature was kept at 300 K for the non-Drude particles and at 1 K for the Drude degrees of freedom. During the production run, no coordinates were saved. We note that equilibration periods would need to be much longer if any properties were to be evaluated from the produced trajectory. Since we were only interested in computation times, the equilibration periods were kept as short as possible.  Table I lists the computation times for a 1 ns MD simulation employing a different number of ion pairs, as well as different thermostats. Fig. 1 shows a graphical representation of the computation times. Regardless of the system size, the GPU enhanced code speeds up polarizable simulations by a factor of nearly five, and nonpolarizable simulations by a factor of nearly four. On both CPU and GPU a polarizable simulation takes about double the computation time than an equivalent nonpolarizable simulation. This is equivalent to polarizable simulations in NAMD, where the addition of polarizability also approximately doubles the simulation cost. 21 We note, however, that the GPU based CHARMM DOMDEC 24 speeds up nonpolarizable simulations even more, so that our comparison is only valid for conventional CHARMM. Furthermore, on both CPU and GPU the computation time increases by a factor of two if the system size is doubled. The dual Langevin and dual Nosé-Hoover thermostats perform equally well.   (17) where the first term corresponds to the positively charged Drude particle at the nucleus, and the second term to the negatively charged Drude particle which can move around the nucleus.  Fig. 4), it is systematic, and occurs in both ionic liquids. We also calculated the distributions for each atomic induced dipole moment (data shown in the Supporting Information), where we find that for k D = 1000kcal/mol/Å 2 larger atomic induced dipoles are more abundant than for k D = 2000 kcal/mol/Å 2 . Constraining the maximum distance of the Drude particle respective to the atom it is attached to (via the 'hardwall' command, maximum distance 0.2 Å) did not affect our finding at all, in fact, during the whole simulation 'hardwall' was not invoked a single time. We thus find very small but systematic differences in the atomic dipole moment distribution depending on the choice of k D , but not on the overall molecular dipole moment.

C Diffusion coefficients and static conductivities
So far, all structural properties were largely independent of the choice of m D , k D and the thermostat. However, dynamic properties are more sensitive towards the simulation setup, and are analyzed in the following. First, the diffusion coefficients of the cations, D ⊕ , and the anions, D ⊖ , were calculated from the derivative of the mean squared displacement where ΔD is a finite size correction. Due to the limited size of the simulation box of length L, long-range hydrodynamic interactions restrict diffusion. 48  where T is the temperature of the system, k B the Boltzmann constant and η the experimental viscosity of 14.54 mPas. 50 Furthermore, the static conductivity σ(0) was calculated from the mean-squared displacement of the translational part M j (t) of the collective dipole moment (20) where the charge of the molecule q i is +1 for cations or -1 for anions, and r i is the center of mass of molecule i. The conductivity then arises from the slope of the mean-squared where V is the volume of the simulation box. Table II lists  would further increase the conductivity and thus decrease the quality of some of the results, especially for the system m D = 0.2 amu and k D =1000 kcal/mol/Å 2 . Also, decreasing the friction would further increase the diffusion coefficients, which are already too large compared to the NH thermostat in most systems. Also, too low friction coefficients may cause a poor temperature control. 22 The optimal friction coefficients thus differ for the analyzed diffusive properties. We furthermore note that the choice of m D and k D largely affects the observed results, especially in the case of the L2 thermostat. Whereas different parameters for the NH thermostat produce consistent diffusion coefficients and conductivities, which differ only due to statistical uncertainties, the Langevin thermostat yields largely different diffusion constants and conductivities dependent on m D and k D . While it might be possible to find a combination of m D , k D and friction coefficients which resembles diffusive properties within reasonably small deviations from the true values, this combination will likely not carry over to a different system. We, therefore, do not recommend a dual Langevin thermostat to determine diffusive properties of ionic liquids.

D Thermostatting of the Drude degrees of freedom
To test whether the choice of parameters affects the temperature of the Drude degrees of freedom, and thus the thermostatting of the Drude particles, we calculated the temperature of each Drude particle T iβ D as Heid et al. Page 12 where v iβ D is the relative motion of the Drude particle to the motion of the attached atom, i.e.
where the motion of the main atom, which corresponds to a temperature of 300 K, was subtracted from the motion of the Drude particle to yield the relative Drude motion which should correspond to a temperature of 1 K. The temperature was then averaged over all atoms β of all molecules i in the system. For all choices of thermostats (NH, L1 and L2), m D (0.2 amu, 0.4 amu) and k D (1000 kcal/mol/Å 2 , 2000 kcal/mol/Å 2 ) the resulting averaged temperature was below 1 K in both [C 2 mim]DCA and [C 2 mim]OTf, so that sufficient thermostatting can be assumed. However, inspecting the cations and anions separately, the combination 0.4 amu with 1000 kcal/mol/Å 2 leads to slightly too large cation Drude temperatures (data shown in the Supporting Information). We thus assume that the deviations in the diffusive properties caused by this combination of m D and k D are the result of poor thermostatting of the imidazolium cation, leading to slightly overestimated induced atomic dipoles.

E Structural and diffusive properties: CPU versus GPU
For the system [C 2 mim]DCA, m D =0.2amu with k D =1000kcal/mol/Å 2 and the dual Nosé-Hoover thermostat, a 15 ns trajectory was also calculated on CPU. The radial distribution functions, dipolar distributions, diffusion constants, and static conductivities were found to be, within statistical uncertainties, identical to the data obtained from simulations calculated on GPU. Diffusive properties are given in Table II, structural properties in the Supporting Information. We note, that due to the poor statistics of the static conductivity calculation, the deviation is quite large; however, the diffusion coefficients, which can be evaluated with much better statistics, agree well between the two different codes.

F Optimality of Drude positions
To check whether the choice of parameters and thermostats affects the ability to recover the true (self-consistent field) Drude positions, we randomly selected 20 configurations from each [C 2 mim]DCA trajectory, omitted movement of all non-Drude particles and minimized the energy by a steepest descent approach followed by a Newton-Raphson method. In a previous study, we found that changes in the total energy of up to 1% still produced trajectories in close agreement with the self-consistent field limit. 28 Table III shows the percentage of decrease in the Drude bond energy, which is largely counteracted by an increase in electrostatics, so that the change in total energy is very small. In general, the dual Langevin thermostat produces Drude positions slightly closer to the self-consistent field limit than the Nosé-Hoover thermostat. Also, a k D of 2000 kcal/mol/Å 2 performs slightly better than a k D of 1000 kcal/mol/Å 2 . However, all combinations of m D , k D and thermostats are sufficiently close to self-consistency.

V Conclusion
We have shown that structural and diffusive properties of ionic liquids can be simulated accurately with different combinations of the Drude mass m D and the Drude force constant k D when employing a dual Nosé-Hoover thermostat. However, we do not recommend the combination of high Drude masses with low force constants, as recommended in the CHARMM documentation, for the simulation of ionic liquids since it led to slight deviations in diffusive properties of [C 2 mim]DCA, and, too less extent, of [C 2 mim]OTf. Also, slightly too large induced atomic dipoles, and slightly too large temperatures of the cationic Drude particles were found, especially for the atoms in the aromatic ring of the [C 2 mim] cation. The dual Langevin thermostat correctly recovered structural properties for both ionic liquids regardless of the choice of m D , k D or the Langevin friction coefficients, but failed for diffusive properties such as ion self-diffusion coefficients or conductivities. The dual Nosé-Hoover thermostat should thus be preferred for the simulation of ionic liquids dynamics.
The high computational effort for polarizable MD simulations of ionic liquids was overcome by using a GPU enhanced CHARMM code, which accelerated the simulations by a factor of five regardless of system size. The CPU/GPU mixed-precision code produced the same results for static and dynamic observables within statistical uncertainties and thus provides a quick and convenient opportunity to use polarizable MD simulation also for large and viscous systems, for example, the simulation of ionic liquids.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. and different thermostats (NH, L1, L2). Inset: Enlarged view on the first peak of selected radial distribution functions.    a 95% confidence intervals were estimated as mean(x) ± stdev(x). t t/ n with t being the Student-t factor. The slope of (Δr 2 (t)) was fit at intervals