# Charge transfer and polarisability in ionic liquids: a case study

Frederik Philippi, Kateryna Goloviznina, Zheng Gong, Sascha Gehrke, Barbara Kirchner, Agílio Pádua, Patricia Hunt

► To cite this version:

Frederik Philippi, Kateryna Goloviznina, Zheng Gong, Sascha Gehrke, Barbara Kirchner, et al. Charge transfer and polarisability in ionic liquids: a case study. Physical Chemistry Chemical Physics, 2022, 24 (5), pp.3144-3162. 10.1039/d1cp04592j. hal-03644382

HAL Id: hal-03644382

https://hal.science/hal-03644382v1

Submitted on 10 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
ARTICLE

Charge Transfer and Polarisability in Ionic Liquids: a Case Study

Frederik Philippi,ª Kateryna Goloviznina,ᵇ Zheng Gong,ᵇ Sascha Gehrke,ᶜ,ᵈ Barbara Kirchner,ᶜ Agílio A. H. Pádua,ᵇ Patricia A. Hunt*ª,ᵉ

Received 00th January 20xx,
Accepted 00th January 20xx

DOI: 10.1039/x0xx00000x

The practical use of ionic liquids (ILs) is benefiting from a growing understanding of the underpinning structural and dynamic properties, facilitated through classical molecular dynamics (MD) simulations. The predictive and explanatory power of a classical MD simulation is inextricably linked to the underlying force field. A key aspect of the forcefield for ILs is the ability to recover charge based interactions. Our focus in this paper is on the description and recovery of charge transfer and polarisability effects, demonstrated through MD simulations of the widely used 1-butyl-3-methylimidazolium bis(trifluoromethanesulfonyl)imide [C₄C₂im][NTf₂] IL. We study the charge distributions generated by a range of ab-initio methods, and present an interpolation method for determining atom-wise scaled partial charges. Two novel methods for determining the mean field (total) charge transfer from anion to cation are presented. The impact of using different charge models and different partial charge scaling (unscaled, uniformly scaled, atom-wise scaled) are compared to fully polarisable simulations. We study a range of Drude particle explicitly polarisable potentials and shed light on the performance of current approaches to counter known problems. It is demonstrated that small changes in the charge description and MD methodology can have a significant impact; biasing some properties, while leaving others unaffected within the structural and dynamic domains.

Introduction

The molecular level interactions that lead to the macroscopic
properties of ionic liquids (ILs) span orders of magnitude in both
length and time scale.¹⁻⁴ For classical molecular dynamics (MD)
simulations the limiting factor is the quality of the force field
which depends strongly on approximations (at the atomic level)
in the underlying model. Classical MD simulations have
primarily involved fixed partial charges placed on individual (or
groups of) atoms and MD simulations of this type have provided
valuable insight into the structure and dynamics of ILs.⁵⁻⁷

In contrast to classical MD, *ab-initio* simulations avoid the use
of arbitrary fixed partial charges, however at the expense of

computational efficiency. Nevertheless, *ab-initio* simulations
are becoming more feasible and are employed within the ionic
liquids community.⁸⁻¹² *Ab-initio* level simulations can also be
used in a complementary manner as a benchmark for classical
MD simulations.⁸,¹³

A deep-seated problem in classical MD potential development
is identifying the best choice for atomic partial charges. ¹⁴
Moreover, the validity and efficacy of fixed partial charges,
particularly for the study of ILs is in question. ¹⁵⁻¹⁸ Local
interactions of the ions in an ionic liquid (IL) can lead to charges
smaller than those formally assigned (typically ±1). ¹,¹⁰,¹⁹⁻²²

Unravelling the underpinning contributions to charge
movement within and between ions in an IL is challenging.
Conceptually different types of charge transfer can be
considered. "Molecular charge transfer" occurs when electron
density is moved from one ion to another due to molecular
orbital overlap between adjacent ions. Molecular charge
transfer can fluctuate as a result of transient molecular orbital
overlap due to ion movement within the liquid. In addition, if a
polyatomic ion loses or receives charge, this is not distributed
evenly over the whole ion, some atoms will lose or receive more
charge than others.

"Environmental charge transfer" is an apparent change in
charge due to changes in the local electrostatic potential or
dielectric environment. For example, a test charge approaching
an ion of charge +Q in a low dielectric "sees" the full charge,
but in a high dielectric the ion will be shielded/screened and the
test charge "sees" a reduced charge. As ions move, the local
electrostatic screening potential will change and charge
fluctuation will occur.

Electronic Supplementary Information (ESI) available: LAMMPS data and input files, atom labelling, atom-wise charge scaling details, simulation densities and box sizes, thermostat details, energy drift, CL&Pol implementation details, coverage of conformational space, solvation model influence on charges, derivation of theoretical charge transfer, atomic charges and histograms, ESP reproducibility, interpolation of atomic charges, visualisation of charges, fukui functions and orbitals, Coulomb interaction energy (derivation, plot, fit parameters), critical points of radial distribution functions, Nernst-Einstein equation and charge transfer, Portuguese flat plot, dihedral angle occurrence and dihedral potential energy surfaces, charge arm distribution histograms, charge arm lengths as a function of dihedral angles, diffusion coefficients and mean squared displacements, performance tests. See DOI: 10.1039/x0xx00000x

<footnote>ª Department of Chemistry, Molecular Sciences Research Hub, Imperial College London, White City Campus, London W12 0BZ, UK.</footnote>

<footnote>ᵇLaboratoire de Chimie, École Normale Supérieure de Lyon & CNRS, 69364 Lyon,
France.</footnote>

<footnote>ᶜ Mulliken Center for Theoretical Chemistry, University of Bonn, Beringstr. 4+6, D-53115 Bonn, Germany.</footnote>

<footnote>d. Department of Physics and Astronomy, University College London, London, WC1E 6BT United Kingdom.</footnote>

<footnote>e. School of Chemical and Physical Sciences, Victoria University of Wellington, Wellington, New Zealand.</footnote>
Moreover, the total charge on the ion can remain unchanged while charge within the ion is redistributed due to "polarisation" induced by local electric fields originating from neighbouring ions. Induction or polarisation is the induced formation of a dipole (or higher multipole) on the (whole polyatomic) ion due to neighbouring ions. In polyatomic ions, polarisation can be responsible for atomic charge transfer, as the whole ion is polarised charge is increased on some atoms and reduced on others. In general, an asymmetric polyatomic ion will already have a dipole, thus the inductive effect will be to change the dipole relative to the dipole present in the absence of an electric field.

Where ions are connected through orbital overlap or H-bonding, a larger intermolecular interaction network is present and polarisation across the network can result in real charge transfer from one ion to another.²¹ Orbital overlap can also be facilitated through orbital (electron) polarisation, and thus charge transfer and polarisation are inter-twined.

There is an ongoing investigation into the importance and effects of charge transfer vs charge polarisation in ILs.¹⁹⁻²¹ The extent of charge transfer and/or polarisation is highly dependent on the specific cation and anion. A strong H-bonding anion will exhibit a larger extended network and increased polarisation induced charge transfer.²¹ A highly symmetrical and diffuse anion will exhibit stronger molecular polarisation but little charge transfer.

We assign as "static" charge effects those that are recovered by fixed partial charge potential models. In static partial charge potentials, atomic partial charges cannot be redistributed (increase or decrease). Moreover, fixed charges cannot move off the atomic centre, nor can they be augmented by local multipoles, for example, an atomic dipole cannot form. Important charge contributions are neglected in MD simulations based on static partial charge potentials, leading to an over representation of electrostatic interactions and systematically biased results.¹⁸,²³,²⁴ Consequently a wide range of mechanisms have been employed to address the effects of charge transfer/polarisation.²³

One straightforward way of recovering some of the charge transfer and polarisation contributions is to adapt the atomic charges to emulate the *average effects* of charge transfer and polarisation.²⁵⁻²⁹ The Coulomb potential $U_{Coulomb}$ is defined in Eqn (1) where $q^C$ is the charge on the cation, $q^A$ is the charge on the anion.¹⁸ Typically $q^C$ and $q^A$ are formally assigned based on the molecular species, ie as ±1 ions. Charges can then be scaled uniformly, for example with a factor of $1/\sqrt{\epsilon_r} \approx 1/n$. For [C₄C₁im][NTf₂] $n^{25°C} \approx 1.44$, which would correspond to a scaling factor 0.69. Scaling factors around 0.8 are commonly used in MD simulations of ionic liquids.¹⁵,¹⁹,²⁰,²³,²⁷,²⁸,³⁰⁻³³ The rationale behind the scaling is that the charge transfer and polarisation effects are encapsulated within the empirical dielectric constant $\epsilon_r$. $q^C$ and $q^A$ are now "effective" charges, taking into account a range of underlying processes. Charge scaling has been shown to achieve the desired acceleration of dynamics, and in particular an increase in diffusion. However, recent work has shown that scaled charge models fail on a local

level, where a model with explicit atomic level polarisability delivers more accurate results.¹⁸,²³,³⁴

$$U_{Coulomb} = \frac{q^C q^A}{4\pi\epsilon_0\epsilon_r r} \qquad (1)$$

A range of methods exist to include explicit atomic level recovery of polarisation including; fluctuating point charges, induced point dipoles and Drude-dipoles.¹⁸,²⁰,²³,³⁵⁻³⁷ Fluctuating point charges are rarely implemented in molecular dynamics code, and the treatment of the interactions involving mathematical point dipoles under periodic conditions adds considerable complexity to a simulation. The interested reader is referred to a recent review on this topic.²³ We will focus on the Drude-dipole methodology, which can be added onto existing fixed charge non-polarisable models and which is already implemented in many efficient software packages.³⁸,³⁹ Heavy atoms (non-H) are split into two particles, a Drude core (DC) fixed on the atomic centre and a mobile Drude particle (DP), connected to the Drude core by a harmonic potential. The Drude core and particle together have the charge and mass of the original atom. The Drude particle is considerably lighter and oppositely charged with respect to its DC. The Drude particle can move around the atomic centre creating an atomic dipole, and hence atomic inductive polarisation is recovered. Polarisation components that require dynamic changes of the partial charge on specific atoms during a simulation are not recovered.

Average as well as explicit models for polarisability in molecular dynamics simulations are not only important for an accurate description of pure ionic liquids, but also related systems such as ionic solutions,⁴⁰,⁴¹ electrode surfaces in contact with ionic liquids,⁴² and miscibility studies of ionic liquids in electrolytes.⁴³ Combining potentials developed independently (i.e. for IL vs aqueous systems), needs to be evaluated carefully before use.⁴¹ It is already well established that polarizable salt potentials are strongly preferred when water is present,⁴⁰,⁴⁴,⁴⁵ a preference that transfers over to water/IL systems.⁴⁰,⁴⁶

Our aim in this work is to explore the impact of different charge descriptions, first we consider different static partial charge assignments and the impact of charge scaling, then we probe the recovery and impact of atomic charge transfer and polarisation effects. Classical MD simulations are carried out on a representative IL system: [C₄C₁im][NTf₂]. We use the standard fixed partial charge force field developed by Canongia Lopes and Pádua as a reference, denoted by CL&P.³⁹,⁴⁷ The CL&P force field has been highly successful in modelling a wide range of ILs.⁴⁷,⁴⁸ We systematically change the atomic charges to those obtained from the CHELPG,⁴⁹,⁵⁰ ADCH,⁵¹ MSK,⁵²,⁵³ and NPA⁵⁴ schemes. CHELPG and MSK are based on the electrostatic potential while ADCH and NPA are based on the electron density distribution. The impact of charge scaling and choice of scaling parameter is explored.

We also explore the effects of introducing polarisability via Drude induced dipoles, we employ a potential implemented by Clark *et al.* that includes Drude particles without further modification of the underlying CL&P force field, denoted
CL&P+Drude.⁵⁵,⁵⁶ However, for a consistent treatment of the interactions, the original force field Lennard-Jones (LJ) parameters need to be modified to avoid double counting of (the Drude particle recovered) polarisation effects. Thus, we employ the recently published CL&P force field which includes both Drude particles and LJ-parameter adjustments.⁵⁷⁻⁵⁹ A newly implemented temperature-grouped Nose-Hoover thermostat (TG-NH) method is also investigated.⁶⁰ Table 1 summarises the methods explored, the methodology for each method is explained further below. Each (scaled) force field with explicitly mentioned charge scheme is scaled atom-wise, i.e. scaling individual atoms, except for the original CL&P force field which is scaled uniformly. The same molecular level scaling factor, 0.86, is used throughout, the method and choice of scaling factor are explained in detail in the paper.

<table><caption>Table 1: Summary of the methods explored within this work.</caption><tr><td>CL&P</td><td>original potential and charges (developed using CHELPG)</td></tr><tr><td>CL&P scaled</td><td>CL&P with uniformly scaled charges</td></tr><tr><td>CL&P-MSK</td><td>CL&P replacing charges with MSK charges (unscaled)</td></tr><tr><td>CL&P-MSK scaled</td><td>CL&P replacing charges with atom-wise scaled MSK charges</td></tr><tr><td>CL&P-CHELPG</td><td>CL&P replacing charges with CHELPG charges (unscaled)</td></tr><tr><td>CL&P-CHELPG scaled</td><td>CL&P replacing charges with atom-wise scaled CHELPG charges</td></tr><tr><td>CL&P-ADCH</td><td>CL&P replacing charges with ADCH charges (unscaled)</td></tr><tr><td>CL&P-ADCH scaled</td><td>CL&P replacing charges with atom-wise scaled ADCH charges</td></tr><tr><td>CL&P+Drude/NH</td><td>CL&P with Drude particles</td></tr><tr><td>CL&P+Drude/TG-NH</td><td>CL&P with Drude particles, TG-NH applied</td></tr><tr><td>CL&Pol/TG-NH</td><td>CL&Pol with Drude particles, scaled LJ and TG-NH applied</td></tr></table>

# Method

## Molecular Dynamics Simulation

A range of classical MD potential implementations based around the CL&P force field have been examined for [C₄C₁Im][NTf₂]. The standard CL&P force fields for [C₄C₁Im]⁺⁴⁸ and [NTf₂]⁻⁶¹,⁶² are employed, with the atomic nomenclature given in the ESI Section 1. A range of scaling scenarios have been explored which are discussed in the body of the paper, a summary of the equations for determining atom-wise scaled charges is provided in the ESI Section 2. Simulations with 512 cations and 512 anions (box sizes and densities are given in the ESI Section 3) were run at 333.15 K (60°C) using the LAMMPS software package with periodic boundary conditions.⁶³ Details of the equilibration procedure are outlined shortly.

The temperature of 60°C is well above the melting point of [C₄C₁Im][NTf₂] which is -4°C.⁶⁴ Bonds and angles are harmonic and an OPLS type dihedral is used. C-H bonds are constrained using the SHAKE⁶⁵ algorithm with a tolerance of 10⁻⁴ and a

maximum of 20 iterations. The impact of constraining the C-H bonds for the static partial charge CL&P potential has been evaluated, ESI Section 4 and ESI Section 5. A timestep of 1.0 fs was used unless otherwise stated.

The Martyna-Tuckerman-Klein correction was used for the NPT simulations.⁶⁶ Equations of motion in the microcanonical (NVE) ensemble were integrated using the velocity Verlet algorithm. Simulations in the canonical (NVT) and isobaric-isothermal (NPT) ensemble were performed using a Nosé-Hoover (NH) chain thermostat/barostat with a chain length of three for both thermostat and barostat at a pressure of 1 atmosphere. The temperature damping parameter of the NH thermostat was 100 fs for real particles and 25 fs for the internal motions of the Drude core-particle pairs. The impact of the thermostat parameters has been tested for the non-polarized CL&P potential, and we have shown that the NH thermostat reproduces the dynamics of an NVE ensemble, ESI Section 4. The damping parameter for the barostat was 1 ps. The centre of mass velocity of the whole box is subtracted every 1000 timesteps, in the microcanonical trajectories this energy loss was accounted for by rescaling.

Short-range Coulombic pairwise interactions and 12/6 LJ interactions were calculated with a cut-off of 12 Å and a tail correction for energy and pressure. A particle-particle particle-mesh solver was used to evaluate long-range coulombic interactions with a relative RMS error in the per-atom forces of 10⁻⁶ for simulations with Drude particles, and 10⁻⁵ for all other simulations. Pairwise Coulombic and LJ interactions of atoms that are separated by one or two bonds are set to zero, and weighed by 0.5 for those separated by three bonds. For simulations without Drude particles, LJ parameters for the interaction of different atom types were obtained by geometric mixing.

The polarisable potentials employed the same parameters given above, with some modifications as outlined in the following. The USER-DRUDE package was used to perform the polarisable simulations in LAMMPS.³⁹ Simulations with Drude particles employed a timestep of 0.5 fs. All atoms except hydrogen atoms were considered as polarisable, distinct H-atoms are retained, but the polarisability of the hydrogen atoms is summed into the heavy atoms they are bounded to. The potential implemented by Clark *et al.* also includes no Drude particles on the F-atoms, the polarisability is summed into the linked C-atoms.

The mass of the Drude particle (DP) is set to 0.4 a.u. (subtracted from the parent atom) and the force constant $k_{DP}$ between the Drude core and Drude particle is set to 4184 kJ/mol.⁶⁷ The partial charges $q_{DP}$ of the Drude particles were subsequently evaluated using Eqn (2) with atomic polarisabilities $\alpha$ determined from the literature.⁶⁸ The total charge of the Drude core and particle is equal to the initial charge of the nonpolarisable atom. To avoid over polarisation of polarisable atoms in close contact a Thole screening with parameter $\alpha$ = 2.6 is employed to damp short range dipole-dipole electrostatic interactions.³⁹,⁶⁹
$$q_{DP} = \sqrt{\alpha \cdot k_{DP}} \qquad (2)$$

In the CL&Pol potential the LJ potential is modified to avoid double counting of the polarisation effects recovered by the Drude-particle model. LJ parameter interaction well depths inherited from the CL&P database are scaled down as described in ref 57,59, further details can be found in the ESI Section 6.

The motion of the atoms (A), which consist of the Drude core-particle pairs, is thermostatted at T_A=333.15 K. As the electrons are able to instantaneously adapt to nuclear motions, motion of Drude particles with respect to their cores has a thermostat set to 1 K.⁶⁷ The Drude polarisable simulations can suffer from heat flow into the Drude core-particle pair motion, which leads to a violation of the equipartition theorem. The newly implemented temperature-grouped Nose-Hoover thermostat (TG-NH) method solves this problem by using a separate thermostat for the molecular ion (M) CoM motions, this thermostat is also set to T_M = 333.15 K,⁶⁰ further details can be found in the ESI Section 7.

The typical procedure for preparing the simulations was as follows. The initial LAMMPS data files were created using fftool (https://github.com/paduagroup/fftool) and packmol.⁷⁰ The energy was minimised using a conjugate gradient algorithm for a maximum of 100 iterations with a stopping tolerance of 10⁻⁴ in the energy. Atoms were then assigned random velocities consistent with a temperature of 333 K. The system was kept at 333 K for 50 ps, heated (linear ramp in temperature) to 600 K over 50 ps, kept at 600 K for 100 ps, and cooled to 333 K over 50 ps. These four steps were performed in NPT. Subsequently, the system was allowed to relax for 25 ps in NVT, followed by an NPT run for 150 ps, during which the box dimensions were stored every 5 fs. Both parts (NVT followed by NPT) were repeated at least 10 times, and the box size was checked for the absence of drift. The average box size was then taken as the target volume and the system compressed into this average box size over 25 ps, and allowed to further equilibrate in NVT for at least 4 ns. Individual simulations were initiated using this procedure, with the following two exceptions. First, the simulations with different sets of (scaled and unscaled) charges were all started using the positions and velocities from the equilibrated simulation with the original CL&P force field. From these starting points, after assigning the new atomic charges, the simulation was run for 250 ps in NPT, followed by 250 ps in NVE for initial relaxation. Then, the system was allowed to relax for 300 ps in NVT, followed by an NPT run for 600 ps, during which the box dimensions were stored every 10 fs. The box was then compressed to the obtained (local) average box size in NVT over 100 ps. This sequence (NVT, NPT, NVT) was repeated at least 10 times, and the local averages of the box size were checked for the absence of drift. The box was then compressed to the overall (global) average box size over 100 ps, followed by equilibration for 10 ns and a production run of 10 ns total length. Second, the 100 simulations employing CL&P for the determination of the cross- and self-diffusion coefficients were started from the unscaled CL&P simulation after the production run, first running 10m steps in the microcanonical ensemble,

where m is the index of the simulation from 1 to 100. This was followed by a short run (200 ps) in the NVT ensemble, allowing the simulations to diverge (Lyapunov instability), and finally a 10 ns production run was carried out. Data files and example input files for the LAMMPS simulations are provided as supporting information.

Power spectra, SAXS patterns, diffusion coefficients and radial distribution functions (RDFs) have been calculated using the TRAVIS software package, code version Jan 01 2019.⁷¹⁻⁷³ The prealpha software package (source code, executable, input files and manual available on github, https://github.com/FPhilippi/prealpha) was used to calculate velocity correlation functions, mean squared displacements (Einstein Relation), temperatures, total Coulomb interaction energy, charge arm distributions, conformer distributions and to separate out the motion of the Drude particles. Velocity correlation functions and mean squared displacements were calculated using the centre of mass.

## Quantum Chemical Computation

Ab initio calculations have been performed using the Gaussian 09 software package (Revision D.01).⁷⁴ Calculations were carried out on isolated ions (gas-phase) employing a frozen-core UMP2/cc-pVTZ or (RO)B3LYP/6-311+G(2d,p) method with a tightened convergence criteria of 10⁻¹⁰ on the density matrix and 10⁻⁸ on the energy.

512 individual anion and 512 individual cation structures were obtained from a single timestep and "box" of the equilibrated CL&P force field MD. These structures provide an ensemble of conformers mimicking the behaviour of the bulk liquid. The original CL&P force field was used because these structures are required to generate input data (such as the atomic charges) for the other MD simulations (with altered charges). We demonstrate that the 512 structures per ion type reasonably cover the full conformational space of both the anions and cations in ESI Section 8.

CHELPG charges were obtained within G09 from fitting electrostatic potential (ESP) and constrained to reproduce the molecular dipole moment. The CHELPG charges were also evaluated in the presence of a polarisable continuum solvent environment, and did not show significant differences with respect to those evaluated in the gas-phase, ESI Section 9. The Natural Population Analysis (NPA, part of the NBO scheme) and the Atomic Dipole Corrected Hirshfeld (ADCH) charges were obtained from the UMP2 density (dipole moment referenced to the center of mass of the molecule). NPA charges were obtained from NBO 3.1 internal to G09, the ADCH charges were determined using the Multiwfn software package.⁵¹,⁷⁵ For each of the cation and anion structures, the chemical potential of the electrons µ and the hardness η have been calculated at the B3LYP/6-311+G(2d,p) level and evaluated as described in the ESI Section 10.⁷⁶
## Results and Discussion

This paper has two key sections, the first section investigates and discusses different partial charge models, different implementations of charge scaling and different ways of evaluating the charge scaling parameter. MD simulations employing the various charge scaled models are compared with methods that incorporate explicit polarisation and areas of concern are highlighted. In the second section, polarisable models are the focus; key implementation parameters are investigated. All the methods (charge scaled and polarisable) explored are compared for the ability to robustly evaluate diffusion coefficients.

**Choice of partial charge method.** A key concern in the development of MD potentials is the best choice for atomic partial charges. In the following we compare four charge schemes CHELPG, MSK, ADCH and NBO. 512 individual cation and anions were extracted from the CL&P MD simulation and the electronic structure evaluated at the UMP2/cc-pVTZ level of theory. The 512 structures (both cation and anion) represent the accessible conformational space, i.e. the Boltzmann weighed MD PES which is sampled during the simulation. Partial charges were determined using the UMP2 ESP for the CHELPG and MSK schemes, and the UMP2 density for the ADCH and NBO schemes, ESI section 11. Histograms of the atomic partial charges for [NTf₂]⁻ are shown in Figure 1 (MSK is not shown because the charges are very similar to the CHELPG). Charge distribution histograms for CHELPG, ADCH and NBO for [NTf₂]⁻, the corresponding neutral radicals and for a selection of atoms in the [C₄C₁Im]⁺ cation can be found in the ESI section 12.

The CHELPG scheme is widely used; it is the basis for the charges in the CL&P force field. The CHELPG scheme fits charges to the ESP outside of the Van der Waals surface of the molecule, a key problem is the determination of inner atom partial charges which can vary widely without impacting the fit.⁷⁷ Another

disadvantage of the CHELPG scheme is the wide distribution of charges for each atom with a strong dependence on ion conformation, as evidenced in Figure 1a. The least polarisable atoms (fluorine, oxygen) lead to rather sharp distributions, the most polarisable atoms (sulphur, carbon) lead to broader distributions.⁷⁸ In contrast, the distributions of the ADCH scheme and especially the NPA scheme are significantly narrower, Figure 1b,c.

The NBO and ADCH schemes are based on the electronic density and thus evaluate a quantity that is inside the van der Waals radius of the atoms. In the NBO scheme natural atomic orbitals are determined, these diagonalise the 1-electron density matrix, are orthogonal and are formed with a maximum occupancy requirement. Atomic charges are found from summing the populations of the natural atomic orbitals for each atom. No dipole moment constraints are applied in the NBO scheme, typically leading to significantly larger charges, Figure 1c.

The ADCH scheme is based on the Hirshfeld method, which treats the total density as a pro-molecule density composed of the spherically averaged ground state densities of neutral atoms and a deformation density. Atomic charges are based on the weighted deformation density, Figure 1b. Variation is introduced in the spherical averaging of the reference atomic density and by the choice of van der Waals radii used in the evaluation. The ADCH method includes an additional constraint that the molecular dipole is reproduced, leading to charges that are similar to the CHELPG method. There are small differences in the ADCH compared to the CHELPG charges; charges on the nitrogen and oxygen atoms differ more in the CHELPG scheme, and charges on the sulphur atom are considerably lower. In general, electron density based charges do not well reproduce the ESP, and the ESP charges do not well reproduce the electronic density. However, the density based, dipole fitted

<figure><img src="image_2.png" /><figcaption>Figure 1: Histograms of atomic partial charges evaluated at the (RO)B3LYP/6-311+G(2d,p) level for the [NTf₂]⁻ anion for the (a) CHELPG, (b) ADCH and (c) NBO charge fitting schemes, normalised by total number of atoms used for creating the histogram.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Interpolation of ADCH atomic partial charges. a) Between the negatively charged anion ($\Delta q_{total} = 0$) and the neutral radical ($\Delta q_{total} = 1$), b) between the positively charged cation and the neutral radical. The shaded blue region corresponds to $\Delta q_{total} \approx 0.2$.</figcaption></figure>

ADCH charges are an exception, producing a reasonable fit to the ESP (ESI Section 13).

Overall, the ADCH charges provide a good compromise; they are density based but not as polarised as the NBO charges, the ADCH charges also reproduce the ESP well, but are not as sensitive to conformer changes as the CHELPG charges. Thus, the ADCH charges look to be a reasonable choice for use in MD potentials.

**Atom-wise charge scaling.** The charge of an ion can change due to charge transfer from one ion to the other, in a liquid the amount of charge transfer will change as ions associate and dissociate over time. Changes in the local screening environment can alter the "observed" charge and the movement of local ions can generate changing polarisation effects. Thus the "average" charge (over different ion conformers and over time) is relevant in the context of MD simulations. Some of these effects can be recovered in an approximate way by scaling the partial charges obtained for isolated ions.

Typically charge scaling is achieved by uniformly scaling all the charges ($q_i$) on the individual atoms (i) within a molecule. The fixed scaling factor is $(1 - \Delta q_{total})$ where $\Delta q_{total}$ is determined from the *total molecular charge transfer*. However, a better model would allow charge transfer to vary for different atoms ($\Delta q_i$) within a molecular ion. In the following the variable amount of charge transfer exhibited by $[C_4C_1im]^+$ and $[NTf_2]^-$ is demonstrated, and a procedure for obtaining atom-wise scaled charges is outlined.

The neutral analogues (fixed at the ion geometry) of the 512 structures taken directly from the MD simulation, i.e. with one electron added to the cation, or one electron removed from the anion, have been computed. The neutral species are radicals with an un-paired electron, requiring the use of unrestricted quantum chemical methods. To retain consistency, the spin

paired systems have also been evaluated using unrestricted methods. Atomic partial charges have been obtained for the ions and neutral analogues using the CHELPG, NBO, ADCH and MSK methods. The charges on the neutral species represent a vertical ionisation or vertical electron affinity for the removal or addition of one electron. Symmetry equivalent atoms from all 512 structures have been averaged, and the variation given as a standard deviation. For example, for the two sulphur atoms in $[NTf_2]^-$, 1024 atomic charges were averaged. The averages and standard deviations for the two ions and two radicals are given in the ESI Section 11.

The charge on each atom can be determined for both the molecular ion ($q_i^{charged}$, $\Delta q_{total} = 0$) and neutral analogue ($q_i^{neutral}$, $\Delta q_{total} = 1$). Thus, two end-points for zero and one electron of charge transfer in the molecular ion have been obtained, which can be connected by a straight line as shown for the density based ADCH charges of $[C_4C_1im]^+$ and $[NTf_2]^-$ in Figure 2. A similar graph for the ESP based CHELPG charges is in the ESI Section 14. From Figure 2, it is clear that some atoms exhibit larger charge transfer than others; for example the N atom in $[NTf_2]^-$ exhibits strong charge transfer while the S atom does not. The change in charge is counter-intuitive based on relative electronegativity of the atoms (S is a soft atom and N is a hard atom). However, the change in atomic charge is reflecting the impact on the electronic density of removing of an electron from the highest occupied molecular orbital (HOMO) of $[NTf_2]^-$, or correspondingly adding an electron to the lowest occupied molecular orbital of $[C_4C_1im]^+$ (ESI Section 15). The value of the scaled charge $q_i$ on each atom for an intermediate total molecular charge transfer $\Delta q_{total}$ can now be evaluated via interpolation of the straight line. For example, the shaded region blue in Figure 2 identifies $\Delta q_{total} \approx 0.2e$, related to a scaling factor of $\approx 0.8$, which is commonly used for ILs.<sup>15,19,20,27,28,30-33</sup> Using a simple linear relationship $y = mx +$
<figure><img src="image_4.png" /><figcaption>Figure 3: CHELPG atomic charges for a) positively charged [C₄C₁im]⁺ and b) the corresponding neutral radical are combined to yield the condensed finite difference Fukui function c). d) Magnitude of the three-dimensional Fukui function mapped on the 0.02 isosurface of the cation electron density. All calculations at the UMP2/cc-pVTZ level of theory. The geometry shown here is representative and was arbitrarily selected for illustrative purposes.</figcaption></figure>

c where $x = \Delta q_{total}$, the slope $(q_i^{neutral} - q_i^{charged})/\Delta q_{total}$ (for $\Delta q_{total} = 1$) and intercept $q_i^{charged}$ (for $\Delta q_{total} = 0$) lead to Eqn (3).

$$q_i = \Delta q_{total} (q_i^{neutral} - q_i^{charged}) + q_i^{charged} \quad (3)$$

Another way of accessing information on the impact of charge transfer is to consider Fukui functions. The Fukui function $f^+(r)$ is the density difference between the cation and its neutral analogue (with one electron added), Eqn (4).‡ For the anion the analogous Fukui function $f^-(r)$ is the density difference between the anion and its neutral analogue (with one electron removed). The condensed Fukui function resolves the 3-dimensions density difference into single point charges on each atom Eqn (5). However, while the electron density used in Eqn (4) is an observable, and thus $f^+(r)$ is well defined, this is not the case for $f_i^+$, which depends on the choice of the atomic partial charge scheme.

$$f^+(r) = \rho^{neutral}(r) - \rho^{charged}(r) \quad (4)$$

$$f_i^+ = q_i^{neutral} - q_i^{charged} \quad (5)$$

MD potentials are (nearly always) parametrised using ESP based charges and to be consistent with existing parametrisations the same charge schemes should be employed. The condensed Fukui functions can be approximated using ESP derived charges, taking the difference between the atomic partial charges of the charged and neutral ions.⁷⁹⁻⁸¹ The vertical electron density distribution and ESP based charges have been evaluated at the

UMP2/cc-pVTZ level of theory, illustrated in Figure 3 [C₄C₁im]⁺, for a representative example of one of the 512 structures. The full sets of averaged atomic charges for the two ions and their neutral counterparts as well as the condensed Fukui function values and the interpolated atom-wise scaled charges are given in ESI Section 11. Figure 3a and 3b illustrate the colour mapped CHELPG charges for the [C₄C₁im]⁺ cation and its neutral radical analogue respectively. Figure 3c is the CHELPG based condensed Fukui function. Figure 3d shows the ESP based Fukui function (Magnitude of the three-dimensional Fukui function mapped on the 0.02 isosurface of the cation electron density). After the addition of an electron the cation has become less positive overall, however some atoms have become more positive (i.e. lost negative charge), in particular the N atoms are more positive (green in Figure 3c). The imidazolium ring atoms are impacted more than the alkyl side chain atoms, with the C2 carbon atom becoming significantly more negative (red in Figure 3c), see also Figure 2b. Again counter-intuitive results (N atoms are more positive) strongly reflect changes in density, in this case due to occupation of the LUMO orbital, (ESI Section 15). An analogous analysis has been carried out for [NTf₂]⁻ and is described in ESI Section 15.

Thus, it can be seen that some parts of the ion are more significantly affected by charge transfer than others, and that changes in atomic partial charge can be the reverse of those that would be applied via a uniform scaling parameter. However, there is also an underlying issue that the Fukui function is density based and the charge analysis is ESP based, and that ESP based charge schemes do not well reflect density
<figure><img src="image_5.png" /><figcaption>Figure 4: Distribution of Δq<sub>ion</sub> determined from 512 cations and 512 anions computed at the R(O)B3LYP/6-311+G(2d,p) level of theory. The blue arrow indicates the charge transfer estimate obtained from the two (individually) optimised structures (shown in the top left).</figcaption></figure>

changes.<sup>82–84</sup> In this context the ADCH charges being density based, and producing a good ESP fit, are a more optimal choice in terms of potential development.

**Determining the magnitude of charge scaling.** A value for the molecular level charge transfer ($\Delta q_{total}$) is required to determine the individual atomic charge transfers ($\Delta q_i$). $\Delta q_{total} = \Delta q_{ion}$ can be estimated using total energy derivatives via Eqn (6).<sup>85,86</sup> Here, $\Delta q_{total}$ is the charge transfer experienced by the whole molecular ion, i.e. not atom-wise. $C$ denotes cations, $A$ denotes anions, $\mu$ is the chemical potential of the electrons within the molecular ion, and $\eta$ is the hardness of the molecular ion, further details are in ESI Section 10.

$$ \Delta q_{total} = \Delta q_{ion} = \frac{\mu^{C} - \mu^{A} - \frac{\alpha e^{2}}{2\pi\varepsilon_{0}d}}{2\eta^{C} + 2\eta^{A} - \frac{\alpha e^{2}}{2\pi\varepsilon_{0}d}} \quad (6) $$

$$ \mu = \frac{E_{N_0-1} - E_{N_0+1}}{2} \quad (7) $$

$$ \eta = \frac{E_{N_0-1} + E_{N_0+1} - 2E_{N_0}}{2} \quad (8) $$

The amount of charge transfer depends both on the ability (cost in energy) of the anion to release charge and the cation to receive charge. Chemical potential is the gradient of energy $E$ with respect to charge, a low gradient will indicate there is little change in energy for a change in charge and thus a large magnitude of charge transfer will be easy. Hardness is the curvature or second derivative of $E$ with respect to charge, low curvature is indicative of an extended low gradient and "soft" ions. The chemical potential and hardness were computed using Eqn (7) and Eqn (8) (derivations in the ESI, section 10): $E_{N_0}$ is the energy of the singly charged ion, $E_{N_0-1}$ is the energy of the radical with one electron removed, and $E_{N_0+1}$ is the energy of the radical with one electron added.

A Coulomb correction can be included in Eqn (6) to account for losses in electrostatic attraction due to charge transfer, ESI, section 10 and 16. $\alpha$ is a Madelung Constant computed from the Coulomb interaction energy. For $\alpha = 0$ no correction is made, for $\alpha = 1$, a simple unshielded Coulomb correction at distance $d$ is included. To the best of our knowledge, Eqn (6) has only been used for ILs with $\alpha = 0$ and $\alpha = 1$.<sup>76,86</sup> $\alpha$ is therefore a critical parameter which we will show has a pronounced impact on the predicted charge transfer. As $\Delta q_{ion}$ is based on isolated ions, it does not fit the definition of "environmental" or "molecular" charge transfer. Instead, it is purely based on the energetics of hypothetical non-integer electron transfer between isolated ions.

The 512 individual cation and anions previously extracted from the CL&P MD simulation were evaluated at the (gas phase) restricted open-shell (RO)B3LYP/6-311+G(2d,p) level. The chemical potential and hardness were computed using Eqn (7) and Eqn (8) and the total charge transfer estimated according to Eqn (6) with $\alpha = 0$, the distribution of $\Delta q_{ion}$ is given in Figure 4. An average $\Delta q_{ion} = 0.36 \pm 0.01$ e is determined. The charge transfer distribution has a breadth of $\approx 0.025$ e which is relatively narrow indicating that $\Delta q_{ion}$ varies minimally with respect to changes in ion conformation. For comparison, the lowest energy (optimized) B3LYP/6-311+G(2d,p) geometry of each ion gives a $\Delta q_{ion} = 0.35$ e.

A computed $\Delta q_{ion} = 0.36$ e corresponds to a scaling factor of $(1 - \Delta q) = 0.64$ e. So far $\Delta q_{ion}$ has been computed using information from "isolated" ions ($\alpha = 0$), the Coulomb correction depends on $\alpha$ and the distance $d$ between the ions. As a distance estimate, we use 5.5 Å, which is the position of the first peak in the cation-anion radial distribution function (RDF) for the MD simulation of $[C_4C_1im][NTf_2]$ using the original CL&P potential (RDF critical points reported in the ESI Section 17). For the point charge ion pair Coulomb correction with $\alpha = 1$ and $d = 5.5$ Å, $\Delta q_{ion} = 0.14 \pm 0.02$ e.

Eqn (6) with $\alpha = 0$, represents the electron population transfer for the hypothetical partial charge transfer between isolated ions, without Coulombic interactions or even orbital overlap. $\alpha = 1$ includes a correction for Coulombic interactions in a cation-anion pair, calculated as point charges. Beyond the ion pair approximation, the Coulomb interaction potential $U = U^C + U^A$ Eqn (9) can be evaluated from MD simulations and used to determine $\alpha$. $U^A$ is the Coulomb interaction energy of all the other ions with a central reference anion (A) and $U^C$ is the corresponding energy for the cation (C).

$$ U = U^C + U^A = \frac{-\alpha e^2}{4\pi\varepsilon_0 d} \quad (9) $$

$U^C$ and $U^A$ are based on the particle (charge) density and pairwise RDFs, a detailed derivation can be found in the ESI Section 16. The screened Coulomb potential has been computed using $\pm 1$ point charges at the centre of charge (for each cation and anion) for MD trajectories employing variations of the CL&P force field; unscaled, uniformly scaled and polarisable, Figure 5. A full graph including ADCH, CHELPG and
<figure><img src="image_6.png" /><figcaption>Figure 5: Distribution of total ion charge transfer $\Delta q_{ion}$ in [C$_4$C$_1$im][NTf$_2$] predicted by Eqn 3 using different values for the coulomb correction term. The histogram for $\alpha = 0$ is also shown in Figure 4.</figcaption></figure>

MK unscaled and atom-wise scaled can be found in the ESI Section 18.

Pronounced oscillations are evident in the Coulomb interaction energy $U$, i.e. the charge of the reference ion is over screened.$^{87-89}$ Moreover, $U$ converges relatively slowly and shows considerable oscillation even 30 Å from the central ion. The Coulomb interaction energy has been fit to an exponentially damped sine function for each force field, ESI Section 19. The decay constant of the exponential envelope function can be interpreted as the Debye length (how far the electrostatic effect of an ion persists in solution).

The asymptotic long-distance limit value of $U$ ($E_0$) has been used in Eqn (9) with $d$ determined from the first maxima in the cation-anion RDF for each force field (ESI Section 17), to determine a stable value of $\alpha$, details can be found in the ESI Section 19. The average $\alpha$ is 1.45(1) across the simulations. Inorganic salts typically exhibit a Madelung constant greater than 1, for example that of NaCl is 1.7476, indicating enhancement of the Coulomb interactions (beyond a simple ion-pair) due to the surrounding ions. This value for the Madelung constant of [C$_4$C$_1$im][NTf$_2$] is comparable to the known Madelung constants of other "ionic liquid" like singly charged organic salts, for example [C$_4$mpyr][NTf$_2$] $\alpha =1.2854$.$^{90,91}$ The value of $\alpha =1.45$ results in a negative $\Delta q_{ion}$ molecular charge transfer. For the full range of MD potentials investigated $\Delta q_{ion}$ varies from -0.01 to -0.09, ESI Section 19. A negative charge transfer implies transfer of charge from the cation to the anion, which is unphysical.

Figure 5 demonstrates the sensitivity of $\Delta q_{ion}$ to the choice of $\alpha$. Without Coulomb correction ($\alpha = 0$) $\Delta q_{ion}$ is at a maximum, with ion pair based Coulomb correction ($\alpha = 1$) $\Delta q_{ion}$ has been significantly reduced. Using $\alpha = 1.45$ results in an unphysical negative $\Delta q_{ion}$. Moreover, using the experimental refractive index of [C$_4$C$_1$im][NTf$_2$] 1.4271,$^{64}$ gives rise to a scaling factor of $(1 - \Delta q) = 1/n = 0.70$ or $\Delta q_{total} = 0.30$ e. These results lead to the conclusion that Eqn 3 is not fully recovering key physical / chemical quantities. The $\mu$ and $\eta$ employed in Eqn 3 are based

<figure><img src="image_7.png" /><figcaption>Figure 6: Coulomb interaction energy for cation+anion in [C$_4$C$_1$im][NTf$_2$] as a function of distance, calculated using ±1 point charges at the centre of charge. The green dashed line shows an exponential envelope. The blue shaded areas correspond to the positions of the first maxima in the cation-anion RDFs (around 5.5 Å) and the cation-cation (C-C) and anion-anion (A-A) RDFs (around 8.5 Å)</figcaption></figure>

on isolated ions, suggesting that intermolecular interactions are a crucial prerequisite for charge transfer, as opposed to the approach purely based on energetics. The Madelung constant term may need to be included in a more subtle manner, incorporated through moderation of the chemical potential and hardness due to the surrounding environment.

**Nernst-Einstein based determination of charge transfer.** The Nernst-Einstein equation is used to describe ionic conduction within low concentration salts, and it is well recognised that this equation does not recover the full physics of concentrated electrolytes. Explaining deviations from the Nernst-Einstein Equation is the subject of an ongoing debate where three major interpretations co-exist.$^{92,93}$ Deviations can be explained by transient or static ion pairing or ion aggregation, reducing the number of charge carriers and/or introducing larger clusters with a reduced charge.$^{94-96}$ Deviations can also be related to a breakdown in the assumption that diffusion is a sum of self-(individual ion) diffusion coefficients. Cross contributions related to correlations between ions (cation-cation, anion-anion, cation-anion) can be introduced.$^{97,98}$ Both the aggregation and the cross contribution models assume integer charges on the ions. A third approach rationalises deviations from the Nernst-Einstein equation as due to charge transfer, and a reduction in the charge carried by each molecular ion.$^{21,76}$ Experimental and MD data can be used as an alternative route to obtain an estimate $\Delta q_{Bulk}$ for the total charge transfer $\Delta q_{total}$. Eqn 10 (derivation details are in the ESI Section 20). Eqn 10 encompasses both the cross contribution and charge transfer approaches. A major advantage of Eqn 10 is that only relative quantities are used, which reduces the impact of systematic errors. $\Delta q_{Bulk}$ is the average charge transfer, and $\pm (1-\Delta q_{Bulk})$ is the average charge transported by an ion over time. Thus,
$\Delta q_{\text{Bulk}}$ is an effective quantity, at any point in time a real ion might receive or donate more or less charge.

$$ \Delta q_{\text{total}} = \Delta q_{\text{Bulk}} = 1 - \sqrt{\frac{1 - \Delta_{\text{NE}}}{1 + \frac{\Lambda_{\text{cross}}}{\Lambda_{\text{self}}}}} \quad (10) $$

Here, $\Delta_{NE}$ is the deviation from the Nernst-Einstein ionic conductivity, which is accessible experimentally (via pulsed-field gradient NMR and impedance spectroscopy). A $\Delta_{NE} \approx 0.366$ has been experimentally evaluated (at the simulation temperature of 333 K).⁹⁹ Experimental measurements will include the effects of molecular and environmental charge transfer as well as those originating in the polarisability of the ions. The ratio of cross and self contributions to the molar electrical conductivity ($\Lambda_{cross}$ and $\Lambda_{self}$) can be obtained from MD simulations.

Diffusion coefficients (Table 2) have been determined using the Green-Kubo expression as well as the Einstein expression. 100 separate simulations using the CL&P potential, of 10 ns each (thus a total of 1000 ns of simulation time) have been evaluated, full fitting equations can be found in the ESI section 20. Employing Eqn 10, an average total charge transfer $\Delta q_{total} = 0.12 \pm 0.04e$ is determined, Table 2. The few experimental estimates of charge transfer within ILs accordingly support a small, but non-zero charge transfer.¹⁰⁰,¹⁰¹ These small values conceptually fall into the category of “molecular” charge transfer. In contrast, the larger experimental charge transfer values up to $\Delta q_{total} \approx 0.3e$ obtained from the refractive index include “environmental” charge transfer. Future work could explore the (box) size dependence of the simulated $\Lambda_{cross}/\Lambda_{self}$,¹⁰² with scaled or polarisable force fields (which would provide better absolute dynamics).

<table><caption>Table 2: Self- and cross-diffusion coefficients and $\Lambda_{cross}/\Lambda_{self}$ from MD simulations, with a total of 1 µs simulation time in 100 separate trajectories. Simulations were run in NVE, with the exception of the results in the last column, which were obtained from an NVT ensemble. GK=Green-Kubo, ES=Einstein evaluation. Estimated uncertainty in diffusion is 5%. Estimated uncertainty in D is the 95% confidence interval evaluated for the 100 separate simulations. Fit equations: ESI Section 20.</caption><thead><tr><th>Fit type</th><th>NVE, GK<br>Biexp</th><th>NVE, GK<br>Monoexp</th><th>NVE, GK<br>Stretche<br>d exp</th><th>NVE, ES<br>Linear</th><th>NVT, ES<br>Linear</th></tr></thead><tbody><tr><td>$D_S^C$ / $10^{-8}$<br>cm<sup>2</sup>/s</td><td>10.0(5)</td><td>10.2(5)</td><td>9.7(5)</td><td>10.1</td><td>9.6</td></tr><tr><td>$D_S^A$ / $10^{-8}$<br>cm<sup>2</sup>/s</td><td>7.2(4)</td><td>7.3(4)</td><td>6.9(3)</td><td>7.2</td><td>6.9</td></tr><tr><td>$D_{RMMVCF}^A$ /<br>$10^{-8}$ cm<sup>2</sup>/s</td><td>7.0(3)</td><td>7.0(4)</td><td>7.0(4)</td><td>7.3</td><td>6.2</td></tr><tr><td>$\frac{\Lambda_{cross}}{\Lambda_{self}}$</td><td>-0.19(5)</td><td>-0.20(5)</td><td>-0.16(5)</td><td>-0.16</td><td>-0.25</td></tr><tr><td>$\Delta q$</td><td>0.12(4)</td><td>0.11(4)</td><td>0.13(4)</td><td>0.13</td><td>0.08</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: a) Cation-Cation and b) Anion-Anion RDFs obtained from different MD simulations with and without charge scaling compared to MD simulations employing polarisable potentials.</figcaption></figure>

The impact of altering partial charges on the liquid structure. One of the shortcomings of using fixed charges without explicit treatment of polarisability is the overemphasis of long-range electrostatic interactions. In order to evaluate the extent of this effect, MD simulations were carried out using the CL&P potential as a base with the charges altered in the following ways: CL&P partial charges uniformly scaled (CL&P-scaled), ADCH charges unscaled (CI&P-ADCH) and ADCH charges atom-wise scaled (CI&P-ADCH scaled). In particular a uniform scaling has not been applied to the ADCH charges, rather the charges for each atom have been individually scaled through the interpolation process identified in Figure 2, using a charge transfer value of $\Delta q_{total} = 0.14e$. These fixed charge models have been compared with two CL&P based polarisable potentials, one employing Drude particles CL&P+Drude/TG-NH and the other using Drude particles with LJ scaling (CL&Pol/TG-NH).

Ion-ion (cation-cation or anion-anion) RDFs are given in Figure 7, the position of RDF critical points and number integrals can be found in the ESI, Section 17. The first peak in the like-charged ion-ion RDFs characterises the second solvation shell of an ion,
<figure><img src="image_9.png" /><figcaption>Figure 8: Simulated SAXS patterns generated for different MD simulations with and without charge scaling and using polarisable potentials.</figcaption></figure>

this peak is too pronounced for the non-scaled charges (dashed lines). Moreover, although the changes in the second peak of the ion-ion RDF appear small, a much larger number of ions are affected. The first peak in the ion-ion RDFs corresponds to approximately 17-22 neighbours, whereas the second peak corresponds to approximately 45-60 neighbours. The ADCH charges in combination with the atom-wise scaling yields ion-ion RDFs which are very close to those obtained using the polarisable force fields.

Simulated small-angle X-ray scattering (SAXS) patterns have been generated, Figure 8. SAXS patterns of ILs typically have two or three peaks at wave vectors smaller than 20 nm⁻¹ (or conversely distances greater than 3.1 Å). The peak at 14 nm⁻¹ (4.5 Å) arises from correlations of atoms close to each other (inter- and intra-molecular adjacency correlations)⁶ and shows little variation between the different simulation types. In contrast, the "charge" peak at 8-9 nm⁻¹ (7.9 Å) can be attributed to (same sign) ion-ion correlations and is representative of charge ordering in ILs.⁶ In Figure 8 the charge peak intensity decreases significantly with unscaled < scaled < polarisable simulations. The smaller charge peak obtained for the polarisable simulations is consistent with the less pronounced structuring indicated by the lower intensity of the first maxima and disordering of the second maxima in the ion-ion RDFs of Figure 7.⁵⁷ The atom-wise scaled ADCH approach follows most closely the polarisable potential simulated SAXS patterns.

Surprisingly, compared to the other potentials, the CL&Pol/TG-NH simulation shows a qualitatively distinct behaviour at low q (below 6 nm⁻¹, 15-20 Å in real space). Peaks in this region arise from structural heterogeneity in ILs, i.e. the formation of polar and non-polar domains through the aggregation of alkyl chains.⁶,¹⁰³,¹⁰⁴ This peak (or shoulder) is not observed in the experimental SAXS pattern for [C₄C₁im][NTf₂].¹⁰⁵⁻¹⁰⁷ Nevertheless, it has been suggested that this peak/shoulder could be washed out due to the tail of the charge peak dominating in this region.¹⁰³ Deconvolution of the scattering

pattern for [C₄C₁im][NTf₂] did reveal a peak at d ≈ 12.3 Å,¹⁰⁶ close to the value of 11 Å predicted from the extrapolation of [CₓC₁im][NTf₂] for x = 6-10.¹⁰⁷ A prediction was made for the onset of nano-aggregation at x=3.¹⁰⁷ Thus the shoulder visible in Figure 8 could indicate the onset of domain formation. However, visualisation of the molecules in the simulation does not show the expected aggregation of non-polar groups, ESI Section 21. Thus, the existence of nano-scale segregation for x=4 is not resolved, and if a structural feature does exist within the IL it may have unique or different characteristics compared to the nano-aggregation of alkyl domains observed for x≥6. This aspect could be further explored in future work through more detailed analysis of MD simulations with the help of partial SAXS patterns.

The CL&Pol/TG-NH and CL&L+Drude/TG-NH potentials differ in the adjustment of the LJ parameters and in the treatment of the F-atoms. The CL&Pol force field uses Drude particles on fluorine atoms, whereas Clark et al. account for the polarisability of fluorine *via* the neighbouring carbon atom.⁵⁵,⁵⁶ The shoulder might be an artefact of the different fluorine atom treatment, together with the dominance of the anion contribution to the SAXS patterns.¹⁰⁴,¹⁰⁸ It is important to exercise caution when interpreting SAXS patterns since the peaks often are the result of many different contributions and cancellations.⁶ Further investigation would be needed to clarify this issue, which is beyond the scope of this article. Nevertheless, this result highlights the importance of the fluorine atom treatment in obtaining an accurate representation of nano-scale aggregation.

**The impact of altering charge on the conformational space spanned by the MD potential.** We have shown that the choice of atomic charges and (atom-wise) charge scaling can improve the long-range structure within a simulation, mimicking the effects of polarisability. In this work the partial atomic charges of the original CL&P force field have been altered. This approach is common for simulations where charges are uniformly scaled, but also carries the danger of introducing bias into a simulation. Biases arise because charge based non-bonded interactions are altered without concomitant adjustment of relevant dihedral (torsion) angle parameters, which can reasonably be expected to vary with the magnitude of atomic partial charges.²⁵,⁴⁷ As a result, the conformational space sampled by a simulation can change.

In the gas and liquid phase, [NTf₂]⁻ adopts two different geometries, cis and trans.¹⁰⁹ The conformer distribution is best visualised by plotting the probability of finding a certain geometry as a contour plot in the space spanned by the C-S-N-S and S-N-S-C dihedral angles, as shown in Figure 9. Figure 9a is the Boltzmann weighted potential energy surface evaluated at the MP2/cc-pVTZ//RB3LYP-GD3BJ/6-311+G(d,p) level of theory.¹¹⁰ Figure 9b and 9c are MD conformer distributions showing the number of occurrences of a geometry in the whole trajectory. Distributions for MD simulations using the other potentials are shown in the ESI Section 22.
<figure><img src="image_10.png" /><figcaption>Figure 9: Occurrence of the two backbone dihedrals in the [NTf₂]⁻ anion. a) Boltzmann-weighed ab initio surface¹¹⁰, b) classical MD simulation with original charges of the CL&P force field and c) with ADCH charges.</figcaption></figure>

The charge arm distribution for the cation, Figure 10a, shows normalised frequency

the original CL&P force field. However, in contrast, the
conformer distribution obtained with the CL&P force field
employing unscaled ADCH charges is significantly distorted,
Figure 9c. There is minimal difference between the relaxed
potential energy surfaces (ab-initio, CL&P and CL&Pol) for the
torsional angles of [NTf₂]⁻, ESI Section 23, demonstrating that
the Drude particles do not bias the dihedral parameterisation.
Coulombic forces act on the centre of charge of a molecular ion.
The ion will experience torque when the centre of mass and
centre of charge don't coincide. Thus, the distance between
centre of charge and centre of mass, i.e. the charge arm, can be
used to probe the extent to which an ion will interact with
electric fields.⁵⁵,¹¹¹,¹¹² The charge arm can be interpreted as a
dipole moment. We will use the term 'charge arm' to distinguish
this quantity from the dipole moment of neutral molecules.¹¹³

only minor differences between scaled and unscaled charges.
The charge distribution used in the CL&P force field differs
slightly from the optimal CHELPG fitted charges for [C₄C₁im] and
[NTf₂], due to constraints required to ensure the transferability
of the force field.⁴⁸ This leads to a minimal shift by 0.2 Å towards
larger charge arms for the CL&P force field (compare the dashed
red and solid black lines in Figure 10). A fuller Figure including
results from all the charge model MD simulations can be found
in ESI, section 24.

The chloride anion with a greater concentration of negative
charge can H-bond with the alkyl groups and may influence
more significantly alkyl chain rotations and the charge arm. The
[NTf₂]⁻ anion is generally considered to only have weak H-bond
interactions with the cation and stronger steric constraints,
consistent with the shifts in charge arm length determined here

<figure><img src="image_11.png" /><figcaption>Figure 10: Charge arm distribution for a) cations and b) anions for different MD simulations.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: a) Power spectral density, offset vertically for clarity. The red curve only shows internal motions of the Drude particles where the intensity is multiplied by a factor of 500 compared to the other curves. The black curve, CL&P, is the only one where C-H bonds were not kept fixed, showing the high-frequency modes at 3000 cm⁻¹. b) Temperature of the real particles T<sub>R</sub> and the centre of mass motions T<sub>CM</sub>, CL&P+Drude/NH. The black dashed line is the specified thermostat temperature.</figcaption></figure>

as the cis/trans balance of the [NTf₂]⁻ conformers changes. Ab-initio calculations have shown a non-systematic impact on the rotational profile of the three alkyl chain torsion angles for an ion-pair; up to ≈5 kJ/mol and shifting minima by up to 40° (on a shallow PES).¹⁴ Thus, we could expect charge scaling to have an impact on the electrostatic H-bond and hence chain orientation in MD simulations, however we rationalise that the impact at room temperature or above is likely be subtle.

The cis and trans conformers of [NTf₂]⁻ exhibit significantly different charge arm lengths, Figure 10b and ESI section 24. Despite a relatively small scaling factor, the scaled CL&P force field exhibits a significant distortion of the distribution, favouring the cis conformer. Uniform scaling of charges does not change the charge arm length as the centre of charge remains the same, and thus the shift towards a larger charge arm is a clear indication of a biased conformer distribution, cf Figure 9 and ESI section 22 and 24. A biased conformer distribution (shifted to the cis conformer) is also directly evident for the ADCH atom-wise scaled partial charge simulations. Thus, the ADCH atom-wise scaled partial charges compared to the uniformly scaled CL&P force field deliver a better liquid structure interpreted through RDFs (Figure 7) but a worse anion conformer distribution (Figure 10).

The charge arm distribution broadens and shifts to slightly larger values for both cation and anion when polarisability is introduced, (black in Figure 10a and 10b). An advantage of Drude based methods over those employing scaled charges is that the dihedral components of the force field do not need to be adjusted, ESI section 23. The charge on each atomic site (sum of Drude-core and Drude-particle charge) is the same as that used for the parameterisation of the dihedral components.

## Polarisable Models and Diffusion Coefficients

For Drude-particle simulations the use of dual thermostats is normal, one for the motion of the atoms (Drude core-particle pair) T<sub>R</sub> = 333 K and one for the Drude particle-core motions T<sub>D</sub> = 1 K. However, heat flow can occur from the atomic motions into the Drude core-particle motions potentially leading to a violation of the equipartition theorem. Heat flow is most significant from the fast molecular-ion C-H vibrations which effectively couple with the fast Drude vibrations, Figure 11a shows the strong spectral overlap around 3000 cm⁻¹. Consequently, constraining the C-H bond vibrations is an effective measure in reducing energy flow to the Drude particles.

Heat flow from atomic motions into the Drude core-particle creates a distorted kinetic energy distribution. As the average atomic temperature is held constant via the NH thermostat, to compensate for the low temperature in the high-frequency modes there is an increase in the temperature of the low-frequency motions, resulting in a molecular CoM motion that is too high (= too hot).⁶⁰ The effective temperature of the molecular CoM (T<sub>CM</sub>) for a trajectory is shown in Figure 11b (blue), T<sub>CM</sub> is too high by ≈70 K. After C-H vibrations are constrained (at 10ps) T<sub>CM</sub> drops to an average of 360 K. The effective (average) temperature over all atoms T<sub>R</sub> stays at the target temperature with only minor fluctuations in Figure 11b (red).

Despite constrained C-H vibrations, the translational temperature of the ions in Figure 11b is still ≈30 K above the target temperature. The temperature grouped Nosé-Hoover thermostat (TG-NH) scheme solves this problem by using separate thermostats, for the molecular CoM motions and atomic intramolecular motions, each independently thermalised at 333K, see also ESI
<figure><img src="image_13.png" /><figcaption>Figure 12: Log-log plot of the mean squared displacements of the Anions for different simulations. The bold green line shows the expected behaviour assuming normal diffusion using experimental data for the diffusion coefficients at the simulation temperature.⁹⁹</figcaption></figure>

Section 4.⁶⁰,¹¹⁴ The relative motion of Drude particles with respect to their DCs is thermalised at 1K.

Figure 12 shows the mean squared displacement (MSD) $\langle r^2(t) \rangle$ of the anion for a range of simulations employing non-polarised and polarised potentials under varying conditions (C-H constrained, TG-NH approximation). Full graphs (including results from all the charge model MD simulations) for the cation and anion can be found in the ESI Section 25. ILs are usually viscous and heterogeneous, which means that long simulation times are required to reach the diffusive regime. For the simulations of $[C_4C_1im][NTf_2]$ in this work, the (normal) diffusive regime is reached after approximately 1 ns.

Diffusion coefficients have been determined from the slope of a linear fit of $\langle r^2(t) \rangle$ against t from 3 to 7 ns. This fitted region corresponds to 40% of the total trajectory length. We use the Einstein relation, Eqn 12, to calculate the diffusion coefficients from the mean squared displacements. D is the diffusion coefficient and r the distance that a tagged particle travels in a time t. The cornered brackets denote an ensemble average.

Deviations of the parameter a from 1 indicate anomalous diffusion.¹¹⁵

$$ \langle r^2(t) \rangle = 6Dt^a \qquad (11) $$

$$ \log\langle r^2(t) \rangle = \log 6D + a \log t \qquad (12) $$

The computed diffusion coefficients span an order of magnitude, Figure 13, see the ESI Section 25 for numeric data. The simulations with unscaled charges show the slowest diffusion. Charge scaling significantly accelerates the dynamics but only reaches about 50% of the experimental diffusion coefficient with the scaling factors employed in this work. The choice of the charge fitting scheme does not have a pronounced impact on the diffusion (despite the impact on the conformational space sampled by the ions).

The polarisable simulations are sensitive to the specifics for each model. The addition of Drude particles without fixed C-H bonds (CL&P+Drude/NH) leads to a more than 15-fold increase in diffusion, giving higher diffusion coefficients than found experimentally. The reason for this increase is the too high molecular CoM motion due to energy flow into the Drude particles. Consequently, constraining C-H bonds (CL&P+Drude/NH (C-H fix)) reduces heat flow, slowing the molecule CoM motion and reducing the ion diffusion coefficients. Using a temperature grouped thermostat while keeping the C-H bonds constrained (CL&P+Drude/TG-NH) further lowers the diffusion coefficients. The drawback of the CL&P+Drude/TG-NH potential is double counting of the polarisation effects implicitly through the non-bonded attractive empirical LJ potential and explicitly through induced dipoles. The "Pol" force field rescales the LJ parameters to avoid this double counting of inductive polarisation affects.⁵⁷,⁵⁸ The CL&Pol force field with the temperature grouped thermostat and constrained C-H bonds has the best agreement with experimental diffusion data. Corrections such as the Yeh-Hummer finite size correction are often used to further improve predictions.⁶⁰,¹¹⁶,¹¹⁷ Yeh-Hummer corrections are indicated (for the CL&Pol force field) on Figure 13 by the black arrows, the

<figure><img src="image_14.png" /><figcaption>Figure 11: Comparison of the self-diffusion coefficients for a) cation and b) anion as well as c) their ratio. For comparison, the experimental value is shown as dashed grey line.⁹⁹ The arrows indicate the change using the Yeh-Hummer¹¹⁶ finite size correction.</figcaption></figure>
corrections are small compared to the variation with the
simulation setup.

The significant changes in dynamics observed for subtle changes
in the force field or technical details of the simulation are a
current problem for the predictive power of classical MD
simulations. However, the ratio of the cation and anion
diffusion coefficients has been evaluated and found to be
similar for all simulations, Figure 13c. This is a remarkable result,
since some of the simulations in this work were (by choice)
heavily biased, especially the scaled simulations without
reparameterisation of dihedral angles and the CL&P+Drude
simulation without constrained C-H bonds. This observation has
motivated the use of $\Lambda_{cross}/\Lambda_{self}$ in Eqn (10) instead of an
absolute measure such as $\Lambda_{cross}$.

Conclusions

In this work, we have shed light on the interplay between
charge transfer and polarisability in ILs. The importance of these
two phenomena is clear, but the relation between them, and
their realisation in classical force fields, remains the subject of
ongoing inquiry. Charge transfer can be interpreted in a number
of ways, and there are areas where the concepts of charge
transfer and polarisation overlap. Charge transfer can occur due
to molecular orbital overlap and/or become apparent due to
electrostatic screening. Polarisation can also occur; the electric
field of local ions leading to atomic partial charge transfer within
a (fixed charge) molecular ion and molecular charge transfer
over a larger collection of interacting ions.

The two common approaches to improve dynamics in ILs,
charge scaling and explicit polarisability, are closely linked to the
interpretations of charge transfer/polarisability above. In
general, the choice of fixed partial atomic charges gives the user
freedom to include mean field effects of phenomena such as
polarisation. This approach is often chosen over models that
explicitly include polarisability for its simplicity and low
computational cost.

However, we found that scaling of atomic partial charges
neglects the heterogeneity of the organic/inorganic ions used in
ILs, i.e. the fact that some parts of a molecule are more
susceptible to charge transfer than others. Atomic partial
charge transfer can be in the opposite direction to the total
molecular charge transfer, this is due to the dependence of
density changes on the character of the HOMO and LUMO
orbitals (Fukui functions). We thus introduced atom-wise scaled
charges and evaluated their effectiveness in $[C_4C_1im][NTf_2]$ MD
simulations.

We evaluated the magnitude of the total molecular ion charge
transfer ($\Delta q_{total} = \Delta q_{ion}$) predicted by DFT descriptor theory
and explored the variation that can be introduced by including
a correction for local environment Coulomb effects (the
response of local ions providing charge screening) through the
Madelung constant $a$. In the absence of charge screening
$\Delta q_{ion} = 0.36$ e, which decreases to $\Delta q_{ion} = 0.14$ e for an ion
pair correction. $\Delta q_{ion}$ vanishes when charge screening is fully
switched on. This shows that charge transfer does not take

place purely for energetic reasons, and that other mechanisms
such as molecular orbital overlap are necessary.

The presence of charge transfer and/or polarisability will leave
a signature in the macroscopic electrical conductivity.
Deviations from the electrical conductivity predicted by the
Nernst-Einstein equation can be directly interpreted as charge
transfer. However, this approach neglects the effects of cross-
correlation, which have also been shown to be relevant for ILs.
We thus developed a new measure for charge transfer which
accounts for these effects and also remains valid for the limiting
cases of no charge transfer and no cross correlation. Using this
new approach, we determined $\Delta q_{total} = \Delta q_{bulk}$ based on
empirical (NMR) diffusion measurements and individual and
cross-correlation ion diffusion coefficients obtained from MD
simulations (CL&P potential). This led to a $\Delta q_{bulk} = 0.12$ e.

We compared structural and dynamic properties obtained from
MD simulations using different charge schemes (unscaled,
uniform scaling, atom-wise scaling and polarisable). In general,
scaled charges produced structural (Coulomb interactions,
RDFs, SAXS patterns) and dynamic (diffusion coefficients)
properties closer to those of the models with explicit
polarisation. Notably, the ADCH scheme showed the best
correspondence with the polarisable simulations and produced
narrower distributions of the predicted atomic charges
compared to CHELPG. However, the potential energy surface of
the $[NTf_2]_-$ anion showed a pronounced dependence on the
choice of atomic charges. It is thus recommended to re-
parameterise the backbone dihedrals when the atomic partial
charges are changed (for those of another method) or scaled.
Particularly for local properties (such as conformations), it may
be better to stay with poor charges rather than to choose new
charges without re-parameterising the dihedral space.

We believe the best strategy is to re-parametrise the dihedral
angles of the $[NTf_2]_-$ anion for the ADCH atom-wise scaled
charges. Such a potential is expected to deliver; a better
cis/trans distribution, a good density/ESP model, a narrow
charge distribution and a better approximation to the fully
polarizable potential. Variation in the diffusion coefficients or
density between the different charge scaling methods is
negligible. Moreover, cation-anion relative diffusion
coefficients should be analysed where possible due to the
systematic cancellation of errors. For future work, it would be
interesting to explicitly compare atom-wise scaled ADCH and
uniformly scaled CL&P charges with different scaling factors,
but with both having (independently) properly re-
parameterised dihedral spaces. The ADCH charges appear to be
a competitive alternative to CHELPG charges.

We discussed in detail several potential pitfalls of using Drude
induced dipoles for explicit polarisability. As opposed to charge
scaling schemes, a reparameterisation of backbone dihedral
angles is not necessary. However, care must be taken to avoid
(and identify) biases introduced by the thermostat, the violation
of the equipartition theorem, and double counting of induction
effects. Our results furthermore highlighted the importance of
the fluorine atom treatment in obtaining an accurate
representation of nano-scale aggregation.
The diffusion coefficients, which represent the bulk dynamics, showed a pronounced method dependence. Charge scaling and inclusion of explicit polarisability significantly increased diffusion, whereas the difference between different sets of charges with the same scaling factor was negligible. By far the best correspondence between simulation and experiment was found for the CL&Pol force field with a temperature-grouped thermostat, which includes all recent improvements to the well-known CL&P force field. The temperature-grouped NH thermostat ensures equipartition and scaling of the Lennard-Jones ε parameters is used to remove double counting of induction effects.

Although the absolute diffusion coefficients were found to vary significantly across all the simulations, the ratio of the cation and anion diffusion coefficients is more constant. This is a remarkable result, since some of the simulations in this work were (by choice) heavily biased. Thus, there is an advantage to using relative quantities or internal references whenever possible.

Polarisability is a phenomenon of central importance to ILs which needs to be considered to yield accurate classical force fields. Careful choice of atomic charges – including charge scaling – can provide a good approximation to the mean field effects of polarisability. Models with explicit consideration of polarisability can yield accurate dynamic properties, but at the expense of considerable added complexity. Based on the results of this paper, where it is not feasible to employ explicitly polarisable models, it appears reasonable to choose larger scaling factors to emulate environmental charge transfer (including polarizability and screening effects), even though some measures of charge transfer might yield smaller or even negligible values.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

FP thanks Oldamur Hollóczki (Mulliken Center for Theoretical Chemistry) for insightful discussions of charge transfer, Fernando Bresme for advice regarding velocity autocorrelation functions, and Tom Welton for excellent supervision and general advice.

KG thanks IDEX Lyon Fellowship (ANR-16-IDEX-005) and ZG thanks ANR project LIQUI2D (ANR-18-CE09-0018) for the financial support.

This work is funded by the Imperial President's PhD Scholarship. Computational resources provided by the Imperial College Research Computing Service are gratefully acknowledged. DOI: 10.14469/hpc/2232

## Notes and references

‡ The Fukui functions $f^+(r)$ and $f_i^+$ identify regions which are susceptible to nucleophilic attack.<sup>81,118</sup> For the anion, the relevant Fukui function is that for electrophilic attack, defined as

$f^-(r) = \rho^{\text{charged}}(r) - \rho^{\text{neutral}}(r)$. However, Eqn (3) remains the same for the anion, thus the sign of the Fukui function needs to be reverted for the anion if the original definition is followed.

1. B. Kirchner, O. Hollóczki, J. N. Canongia Lopes and A. A. H. Pádua, *Wiley Interdiscip. Rev. Comput. Mol. Sci.*, 2015, **5**, 202–214.

2. K. Dong, X. Liu, H. Dong, X. Zhang and S. Zhang, *Chem. Rev.*, 2017, **117**, 6636–6695.

3. R. Hayes, G. G. Warr and R. Atkin, *Chem. Rev.*, 2015, **115**, 6357–6426.

4. Q. Berrod, F. Ferdeghini, J.-M. Zanotti, P. Judeinstein, D. Lairez, V. García Sakai, O. Czakkel, P. Fouquet and D. Constantin, *Sci. Rep.*, 2017, **7**, 2241.

5. E. J. Maginn, *J. Phys. Condens. Matter*, 2009, **21**, 373101.

6. J. C. Araque, J. J. Hettige and C. J. Margulis, *J. Phys. Chem. B*, 2015, **119**, 12727–12740.

7. Y. Wang, W. Jiang, T. Yan and G. A. Voth, *Acc. Chem. Res.*, 2007, **40**, 1193–1199.

8. M. G. Del Pópolo, R. M. Lynden-Bell and J. Kohanoff, *J. Phys. Chem. B*, 2005, **109**, 5895–5902.

9. E. I. Izgorodina, *Phys. Chem. Chem. Phys.*, 2011, **13**, 4189.

10. M. Bühl, A. Chaumont, R. Schurhammer and G. Wipff, *J. Phys. Chem. B*, 2005, **109**, 18591–18599.

11. J. K. Shah, *Ab Initio Molecular Dynamics Simulations of Ionic Liquids*, Elsevier B.V., 1st edn., 2018, vol. 14.

12. M. Salanne, *Phys. Chem. Chem. Phys.*, 2015, **17**, 14270–14279.

13. A. Szabadi, R. Elfgen, R. Macchieraldo, F. L. Kearns, H. Lee Woodcock, B. Kirchner and C. Schröder, *J. Mol. Liq.*, 2021, 116521.

14. P. A. Hunt and I. R. Gould, *J. Phys. Chem. A*, 2006, **110**, 2269–2282.

15. M. Kohagen, M. Brehm, J. Thar, W. Zhao, F. Müller-Plathe, B. Kirchner, F. Müller-Plathe and B. Kirchner, *J. Phys. Chem. B*, 2011, **115**, 693–702.

16. M. M. Francl and L. E. Chirlian, in *Reviews in Computational Chemistry*, 2007, vol. 14, pp. 1–31.

17. S. Riniker, *J. Chem. Inf. Model.*, 2018, **58**, 565–578.

18. C. Schröder, *Phys. Chem. Chem. Phys.*, 2012, **14**, 3089.

19. T. G. A. Youngs and C. Hardacre, *ChemPhysChem*, 2008, **9**, 1548–1558.

20. C. Schröder, A. Lyons and S. W. Rick, *Phys. Chem. Chem. Phys.*, 2020, 49–52.

21. O. Hollóczki, F. Malberg, T. Welton and B. Kirchner, *Phys. Chem. Chem. Phys.*, 2014, **16**, 16880–16890.

22. K. Wendler, F. Dommert, Y. Y. Zhao, R. Berger, C. Holm and L. Delle Site, *Faraday Discuss.*, 2012, **154**, 111–132.

23. D. Bedrov, J. P. Piquemal, O. Borodin, A. D. MacKerell, B. Roux and C. Schröder, *Chem. Rev.*, 2019, **119**, 7940–7995.

24. T. Yan, C. J. Burnham, M. G. Del Pópolo and G. A. Voth, *J. Phys. Chem. B*, 2004, **108**, 11877–11881.

25. F. Dommert and C. Holm, *Phys. Chem. Chem. Phys.*, 2013, **15**, 2037–2049.

26. V. V. Chaban and I. V. Voroshylova, *J. Phys. Chem. B*, 2015, **119**, 6242–6249.

27. B. Doherty, X. Zhong, S. Gathiaka, B. Li and O. Acevedo, *J. Chem. Theory Comput.*, 2017, **13**, 6131–6145.
28 K. Cui, A. Yethiraj and J. R. Schmidt, *J. Phys. Chem. B*, 2019, **123**, 9222–9229.

29 F. Dommert, K. Wendler, B. Qiao, L. Delle Site and C. Holm, *J. Mol. Liq.*, 2014, **192**, 32–37.

30 M. H. Kowsari and S. Ebrahimi, *Phys. Chem. Chem. Phys.*, 2018, **20**, 13379–13393.

31 X. Zhong, Z. Liu and D. Cao, *J. Phys. Chem. B*, 2011, **115**, 10027–10040.

32 B. L. Bhargava and S. Balasubramanian, *J. Chem. Phys.*, 2007, **127**, 114510.

33 V. V. Chaban, I. V. Voroshylova and O. N. Kalugin, *Phys. Chem. Chem. Phys.*, 2011, **13**, 7910.

34 J. G. McDaniel and A. Yethiraj, *J. Phys. Chem. Lett.*, 2018, **9**, 4765–4770.

35 O. Borodin, *J. Phys. Chem. B*, 2009, **113**, 11463–11478.

36 M. Schmollngruber, V. Lesch, C. Schröder, A. Heuer and O. Steinhauser, *Phys. Chem. Chem. Phys.*, 2015, **17**, 14297–14306.

37 S. Patel and C. L. Brooks, *Mol. Simul.*, 2006, **32**, 231–249.

38 C. Schröder and O. Steinhauser, *J. Chem. Phys.*, 2010, **133**, 154511.

39 A. Dequidt, J. Devémy and A. A. H. Pádua, *J. Chem. Inf. Model.*, 2016, **56**, 260–268.

40 E. Marin-Rimoldi, J. K. Shah and E. J. Maginn, *Fluid Phase Equilib.*, 2016, **407**, 117–125.

41 D. Dhabal and T. Patra, *Phys. Chem. Chem. Phys.*, 2020, **22**, 28325–28338.

42 M. Buraschi, S. Sansotta and D. Zahn, *J. Phys. Chem. C*, 2020, **124**, 2002–2007.

43 J. G. McDaniel, *J. Phys. Chem. B*, 2018, **122**, 4345–4355.

44 G. L. Warren and S. Patel, *J. Phys. Chem. B*, 2008, **112**, 11679–11693.

45 P. B. Petersen and R. J. Saykally, *Annu. Rev. Phys. Chem.*, 2006, **57**, 333–364.

46 G. Chevrot, R. Schurhammer and G. Wipff, *Phys. Chem. Chem. Phys.*, 2006, **8**, 4166.

47 J. N. Canongia Lopes and A. A. H. Pádua, *Theor. Chem. Acc.*, 2012, **131**, 1129.

48 J. N. Canongia Lopes, J. Deschamps and A. A. H. Pádua, *J. Phys. Chem. B*, 2004, **108**, 2038–2047.

49 L. E. Chirlian and M. M. Francl, *J. Comput. Chem.*, 1987, **8**, 894–905.

50 C. M. Breneman and K. B. Wiberg, *J. Comput. Chem.*, 1990, **11**, 361–373.

51 T. LU and F. CHEN, *J. Theor. Comput. Chem.*, 2012, **11**, 163–183.

52 U. C. Singh and P. A. Kollman, *J. Comput. Chem.*, 1984, **5**, 129–145.

53 B. H. Besler, K. M. Merz and P. A. Kollman, *J. Comput. Chem.*, 1990, **11**, 431–439.

54 A. E. Reed, R. B. Weinstock and F. Weinhold, *J. Chem. Phys.*, 1985, **83**, 735–746.

55 R. Clark, M. von Domaros, A. J. S. McIntosh, A. Luzar, B. Kirchner and T. Welton, *J. Chem. Phys.*, 2019, **151**, 164503.

56 C. E. S. Bernardes, K. Shimizu, J. N. C. Lopes, P. Marquetand, E. Heid, O. Steinhauser and C. Schröder, *Phys. Chem. Chem. Phys.*, 2016, **18**, 1665–1670.

57 K. Goloviznina, J. N. Canongia Lopes, M. Costa Gomes and A. A. H. Pádua, *J. Chem. Theory Comput.*, 2019, **15**, 5858–5871.

58 A. A. H. Pádua, *J. Chem. Phys.*, 2017, **146**, 204501.

59 K. Goloviznina, Z. Gong and A. A. H. Pádua, *WIREs Comput. Mol. Sci.*, 2021, 1–16.

60 C. Y. Son, J. G. McDaniel, Q. Cui and A. Yethiraj, *J. Phys. Chem. Lett.*, 2019, **10**, 7523–7530.

61 J. N. Canongia Lopes and A. A. H. Pádua, *J. Phys. Chem. B*, 2004, **108**, 16893–16898.

62 A. S. L. Gouveia, C. E. S. Bernardes, L. C. Tomé, E. I. Lozinskaya, Y. S. Vygodskii, A. S. Shaplov, J. N. C. Lopes and I. M. Marrucho, *Phys. Chem. Chem. Phys.*, 2017, **19**, 29617–29624.

63 S. Plimpton, *J. Comput. Phys.*, 1995, **117**, 1–19.

64 P. Bonhôte, A.-P. Dias, M. Armand, N. Papageorgiou, K. Kalyanasundaram and M. Grätzel, *Inorg. Chem.*, 1996, **35**, 1168–1178.

65 J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, *J. Comput. Phys.*, 1977, **23**, 327–341.

66 G. J. Martyna, D. J. Tobias and M. L. Klein, *J. Chem. Phys.*, 1994, **101**, 4177–4189.

67 G. Lamoureux and B. Roux, *J. Chem. Phys.*, 2003, **119**, 3025–3039.

68 E. Heid, A. Szabadi and C. Schröder, *Phys. Chem. Chem. Phys.*, 2018, **20**, 10992–10996.

69 B. T. Thole, *Chem. Phys.*, 1981, **59**, 341–350.

70 L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, *J. Comput. Chem.*, 2009, **30**, 2157–2164.

71 M. Brehm and B. Kirchner, *J. Chem. Inf. Model.*, 2011, **51**, 2007–2023.

72 M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, *J. Chem. Phys.*, 2020, **152**, 164105.

73 O. Hollóczki, M. Macchiagodena, H. Weber, M. Thomas, M. Brehm, A. Stark, O. Russina, A. Triolo and B. Kirchner, *ChemPhysChem*, 2015, **16**, 3325–3333.

74 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone,

75 B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J.

76 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. J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V.

77 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. Kleene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.

78 E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G.

79 A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A.

80 D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2009.

81 T. Lu and F. Chen, *J. Comput. Chem.*, 2012, **33**, 580–592.

82 F. Philippi, D. Rauber, M. Springborg and R. Hempelmann, *J. Phys. Chem. A*, 2019, **123**, 851–861.
77 C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, *J. Phys. Chem.*, 1993, **97**, 10269–10280.

78 P. Schwerdtfeger and J. K. Nagle, *Mol. Phys.*, 2019, **117**, 1200–1225.

79 R. Balawender and L. Komorowski, *J. Chem. Phys.*, 1998, **109**, 5203–5211.

80 P. Fuentealba, P. Pérez and R. Contreras, *J. Chem. Phys.*, 2000, **113**, 2544–2551.

81 Y. Li and J. N. S. Evans, *J. Am. Chem. Soc.*, 1995, **117**, 7756–7759.

82 P. A. Hunt, *Mol. Simul.*, 2006, **32**, 1–10.

83 P. A. Hunt, B. Kirchner and T. Welton, *Chem. - A Eur. J.*, 2006, **12**, 6762–6775.

84 R. M. Fogarty, R. Rowe, R. P. Matthews, M. T. Clough, C. R. Ashworth, A. Brandt, P. J. Corbett, R. G. Palgrave, E. F. Smith, R. A. Bourne, T. W. Chamberlain, P. B. J. Thompson, P. A. Hunt and K. R. J. Lovelock, *Faraday Discuss.*, 2018, **206**, 183–201.

85 P. Geerlings, F. De Proft and W. Langenaeker, *Chem. Rev.*, 2003, **103**, 1793–1874.

86 F. Mississippi, A. Quinten, D. Rauber, M. Springborg and R. Hempelmann, *J. Phys. Chem. A*, 2019, **123**, 4188–4200.

87 J. G. McDaniel and A. Yethiraj, *J. Phys. Chem. B*, 2019, **123**, 3499–3512.

88 P. Keblinski, J. Eggebrecht, D. Wolf and S. R. Phillpot, *J. Chem. Phys.*, 2000, **113**, 282–291.

89 R. M. Lynden-Bell, *Phys. Chem. Chem. Phys.*, 2010, **12**, 1733–1740.

90 E. I. Izgorodina, U. L. Bernard, P. M. Dean, J. M. Pringle and D. R. MacFarlane, *Cryst. Growth Des.*, 2009, **9**, 4834–4839.

91 D. R. MacFarlane, M. Kar and J. M. Pringle, *Fundamentals of Ionic Liquids*, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2017.

92 M. Kar, N. V. Plechkova, K. R. Seddon, J. M. Pringle and D. R. MacFarlane, *Aust. J. Chem.*, 2019, **72**, 3.

93 D. R. MacFarlane, M. Forsyth, E. I. Izgorodina, A. P. Abbott, G. Annat and K. Fraser, *Phys. Chem. Chem. Phys.*, 2009, **11**, 4962–4967.

94 G. W. Driver, Y. Huang, A. Laaksonen, T. Sparrman, Y.-L. Wang and P.-O. Westlund, *Phys. Chem. Chem. Phys.*, 2017, **19**, 4975–4988.

95 A. Rupp, N. Roznyatovskaya, H. Scherer, W. Beichel, P. Klose, C. Sturm, A. Hoffmann, J. Tübke, T. Koslowski and I. Krossing, *Chem. - A Eur. J.*, 2014, **20**, 9794–9804.

96 K. Ueno, H. Tokuda and M. Watanabe, *Phys. Chem. Chem. Phys.*, 2010, **12**, 1649–1658.

97 H. K. Kashyap, H. V. R. Annapureddy, F. O. Raineri and C. J. Margulis, *J. Phys. Chem. B*, 2011, **115**, 13212–13221.

98 K. R. Harris, *J. Phys. Chem. B*, 2016, **120**, 12135–12147.

99 K. R. Harris and M. Kanakubo, *Phys. Chem. Chem. Phys.*, 2015, **17**, 23977–23993.

100 T. Cremer, C. Kolbeck, K. R. J. Lovelock, N. Paape, R. Wölfel, P. S. Schulz, P. Wasserscheid, H. Weber, J. Thar, B. Kirchner, F. Maier and H.-P. Steinrück, *Chem. - A Eur. J.*, 2010, **16**, 9018–9033.

101 W. Beichel, N. Trapp, C. Hauf, O. Kohler, G. Eickerling, W. Scherer and I. Krossing, *Angew. Chemie Int. Ed.*, 2014, **53**,

3143–3146.

102 Y. Shao, K. Shigenobu, M. Watanabe and C. Zhang, *J. Phys. Chem. B*, 2020, **124**, 4774–4780.

103 H. K. Kashyap, C. S. Santos, R. P. Daly, J. J. Hettige, N. S. Murthy, H. Shirota, E. W. Castner and C. J. Margulis, *J. Phys. Chem. B*, 2013, **117**, 1130–1135.

104 H. V. R. Annapureddy, H. K. Kashyap, P. M. De Biase and C. J. Margulis, *J. Phys. Chem. B*, 2010, **114**, 16838–16846.

105 S. T. Keaveney, T. L. Greaves, D. F. Kennedy and J. B. Harper, *J. Phys. Chem. B*, 2016, **120**, 12687–12699.

106 N. J. Brooks, F. Castiglione, C. M. Doherty, A. Dolan, A. J. Hill, P. A. Hunt, R. P. Matthews, M. Mauri, A. Mele, R. Simonutti, I. J. Villar-Garcia, C. C. Weber and T. Welton, *Chem. Sci.*, 2017, **8**, 6359–6374.

107 O. Russina, A. Triolo, L. Gontrani, R. Caminiti, D. Xiao, L. G. Hines Jr, R. A. Bartsch, E. L. Quitevis, N. Pleckhova and K. R. Seddon, *J. Phys. Condens. Matter*, 2009, **21**, 424121.

108 J. J. Hettige, H. K. Kashyap, H. V. R. Annapureddy and C. J. Margulis, *J. Phys. Chem. Lett.*, 2013, **4**, 105–110.

109 J. N. Canongia Lopes, K. Shimizu, A. A. H. Pádua, Y. Umebayashi, S. Fukuda, K. Fujii and S. Ishiguro, *J. Phys. Chem. B*, 2008, **112**, 1465–1472.

110 F. Mississippi, D. Pugh, D. Rauber, T. Welton and P. A. Hunt, *Chem. Sci.*, 2020, **11**, 6405–6422.

111 M. N. Kobrak and N. Sandalow, *ECS Proc. Vol.*, 2004, **2004-24**, 417–425.

112 H. Li, M. Ibrahim, I. Agberemi and M. N. Kobrak, *J. Chem. Phys.*, 2008, **129**, 124507.

113 M. N. Kobrak and H. Li, *Phys. Chem. Chem. Phys.*, 2010, **12**, 1922.

114 K. Goloviznina, Z. Gong, M. F. Costa Gomes and A. A. H. Pádua, *J. Chem. Theory Comput.*, 2021, **17**, 1606–1617.

115 G. Gradenigo, A. Sarracino, D. Villamaina and A. Vulpiani, *Acta Phys. Pol. B*, 2013, **44**, 899.

116 I.-C. Yeh and G. Hummer, *J. Phys. Chem. B*, 2004, **108**, 15873–15879.

117 E. J. Maginn, R. A. Messerly, D. J. Carlson, D. R. Roe and J. R. Elliott, *Living J. Comput. Mol. Sci.*, 2019, **1**, 1–20.

118 F. Jensen, *Introduction to Computational Chemistry*, John Wiley & Sons, Inc., West Sussex, 2nd edn., 2007.