Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields

Many applications in chemistry, biology, and energy storage/conversion research rely on molecular simulations to provide fundamental insight into structural and transport properties of materials with high ionic concentrations. Whether the system is comprised entirely of ions, like ionic liquids, or is a mixture of a polar solvent with a salt, e.g., liquid electrolytes for battery applications, the presence of ions in these materials results in strong local electric fields polarizing solvent molecules and large ions. To predict properties of such systems from molecular simulations often requires either explicit or mean-field inclusion of the influence of polarization on electrostatic interactions. In this manuscript, we review the pros and cons of different treatments of polarization ranging from the mean-field approaches to the most popular explicit polarization models in molecular dynamics simulations of ionic materials. For each method, we discuss their advantages and disadvantages and emphasize key assumptions as well as their adjustable parameters. Strategies for the development of polarizable models are presented with a specific focus on extracting atomic polarizabilities. Finally, we compare simulations using polarizable and nonpolarizable models for several classes of ionic systems, discussing the underlying physics that each approach includes or ignores, implications for implementation and computational efficiency, and the accuracy of properties predicted by these methods compared to experiments.


INTRODUCTION
The increase in computational power and accessibility of massively parallel architectures combined with the maturation of advanced modeling techniques and force fields allowed atomistic molecular dynamics (MD) simulations to transform into an essential tool for providing molecular scale insight into the structure−property relationships and virtual design of novel materials. Condensed phase ionic systems have attracted attention of the modeling and simulation community due to their applications in chemistry, biology, and energy storage research. 1,2 Room temperature ionic liquids (ILs) and solventin-salt concentrated electrolytes are particularly interesting due to their potential applications in batteries, nanoreactors, and separation mediums. 1,2 Examples of typical ionic liquid cations and anions are depicted in Figure 1, including their commonly accepted abbreviations. High ionic concentration in these materials results in local electric fields polarizing solvent and large ions that often requires either explicit or mean-field inclusion of polarization effects in order for molecular models to be accurate and predictive, especially for sampling far-from-equilibrium structures such as electric double layers or nanoconfinement.

Nonpolarizable Simulations of Ionic Systems
Initial attempts to model these ionic systems using traditional nonpolarizable potential energy functions or force fields, which successfully worked for a variety of nonpolar and polar systems with low salt concentrations, uncovered challenges for the accurate prediction of transport and thermodynamic properties. Numerous studies indicated that additional attention should be paid to polarization treatment in systems with high ionic concentrations. 3−41 While simulations using nonpolarizable models provided important insight into molecular level correlations and structure in ionic systems, 42,43 they also demonstrated that thermodynamic and transport properties predicted from MD simulations are often inconsistent with experiments showing much slower dynamics, with the agreement becoming worse at higher salt concentrations. 43−46 For example, recent simulations by Rajupt et al. of 3 M lithium bis(trifluoro-methanesulfonyl)imide (referred as NTf 2 or TFSI or TFSA) solution in 1,3-dioxolane and 1,2-dimethoxyethane (DOL:DME) mixture predicted Li + and NTf 2 ion diffusion coefficients more than 100 times slower than experiments. 43 The self-diffusion coefficients obtained from simulations using the Canongia-Lopes et al. force field 16,18−21 were an order of magnitude smaller than experimental values 29 for 1-ethyl-3methylimidazolium bis(trifluoro-methanesulfonyl)imide [C 2 mim][NTf 2 ] and significantly slower than in experiment for 1,3-dimethylimidazolium chloride [C 1 mim][Cl] and 1butyl-3-methyl-imidazolium hexafluorophosphate [C 4 mim]- [PF 6 ]. 27 Similarly, a sluggish ion transport has been found in simulations of [ 30 MD simulations of alkylpyridinium-based ILs predicted apparent self-diffusivities that are roughly 10 times lower than experimental values. 33 In simulations of 13 different ionic liquids by Tsuzuki et al. 47 using a modified OPLS force field, selfdiffusion coefficients with deviations from experiments ranging from a factor of almost 10 up to as much as a factor of 40 were obtained. This trend is illustrated in Figure 2, where a correlation between ion self-diffusion coefficients (average for cation and In the IL research community, NTf 2 is the common reference for this anion, while in the electrochemical community working with battery electrolytes, TFSI is the more common notation. Taking this into account, throughout the paper we will use both notations to be consistent with the most common usage in the discussed corresponding application. anion) gained from MD simulations using nonpolarizable force fields and experiments is summarized for various ILs. Similarly, large deviations were observed for ionic conductivity and viscosity indicating sluggish dynamics from simulations using the nonpolarizable force fields without charge scaling.
The systematic trend observed for a variety of ILs, as well as concentrated organic solvent electrolytes clearly indicates that these discrepancies are not a problem of a particular force field or a specific system, but rather that the important physics or interactions are missing in the nonpolarizable simulations. Not surprisingly, these interactions are related to the induced polarization which plays a crucial role in molecular ionic systems as expected. The inclusion of induced atomic polarization, where in addition to fixed partial atomic charges each atom/ molecule has a fluctuating induced dipole responding to the local environment, noticeably improved the description of transport and thermodynamic properties of ionic systems. For example, Figure 2 shows a significantly enhanced correlation between simulation vs. experiment data for the self-diffusion coefficients obtained from simulations using a transferable many-body polarizable force field (APPLE&P). 44

Mean-Field Treatment of Polarization Effects
While the concept of polarizable molecular models existed since the late 1970s, the practical application of this approach was not straightforward due to the lack of generic polarizable force fields, limited availability in popular simulation packages, and a substantial (factor of 3−10) increase in computational costs. Therefore, while simulations using polarizable models were regarded as more accurate compared to those with nonpolarizable models, alternative approaches were investigated to bypass the complexity of explicit induced polarization in MD simulations, aiming at a reduction in the computational cost and ensuring the stability of numerical integration. 48 One of the approaches is to account for polarization implicitly by modifying the parameters of van der Waals dispersion interactions. For example, Koddermann et al. 29 modified the Lennard-Jones parameters from the original Canongia-Lopes et al. force field 16,18−21 to match the description of dynamical and thermodynamic properties of the [C n mim][NTf 2 ] series with n = 1, 2, 4, 6, and 8. Taking into account that the interaction between two induced dipoles separated by distance r scales as r −6 , i.e., the same scaling as the dispersion term in the Lennard-Jones potential, such an approach represents an effective meanfield approximation of induced polarization effects. However, it does not take into account the directionality of interactions with and between induced dipoles and requires a significant modification of the Lennard-Jones parameters, which makes such a force field system specific and less transferable. In another approach, the effective repulsion−dispersion parameters of nonpolarizable models with united atom representation of sp 3 carbons were adjusted to reproduce the data (density, heat of vaporization, and self-diffusion coefficients) predicted from fully atomistic polarizable MD simulations of ILs or electrolytes. 49−51 Unlike the original atomistic polarizable force field that utilized essentially a universal set of repulsion−dispersion parameters for most of the atom types (i.e., only one oxygen type of interaction for all ether compounds, carbonate solvent molecules, and anions), therefore ensuring transferability and predictive capabilities for novel materials, the effective two-body nonpolarizable force field had to introduce numerous chemistryspecific additional repulsion−dispersion parameters leading to limited transferability. The approach of adjusting the van der Waals interactions in nonpolarizable models was taken even further by using the force matching approach to fit numerical, system-specific two-body potential functions based on instantaneous atomic forces predicted from polarizable simulations in the condensed phase. 52 While such approximations of polarization by the two-body terms noticeably improved the description of structural properties (compared to just omitting polarization), it was not possible to completely reproduce the ion transport, highlighting the challenges with incorporating the many-body terms into effective two-body interactions.
Another mean-field approach to effectively take into account induced polarization is to scale ionic charges. This approach is one of the most popular in application to ILs and electrolytes. It was motivated by several ab initio calculations on ion pairs and clusters, which suggested that the net charges on ions have to be reduced due to charge transfer or/and polarization effects. 6,53−58 To define an atomic or molecular charge, a specific approximation for correlating the electron density or wave function distributions obtained from DFT or ab initio calculations with the point charge assignments must be employed (e.g., restrained electrostatic potential (RESP), charge from electrostatic potential (CHELPG), natural population analysis (NPA), etc.). Taking into account that strong polarization effects can lead to significant spatial overlap of electron orbitals, it is expected that without considering polarization explicitly (i.e., imposing the constraint of nonpolarizable models) the partial atomic charge distribution and molecular charges extracted from such calculations will be effectively reduced. However, if the same data from DFT or ab initio calculations are approximated assuming a polarizable model, the effective ionic charges can be kept close to unity while the weakening of interactions between ionic species is captured by explicit induced polarization interactions. 45 Chaban and Voroshylova reported that the cluster size of 1-ethyl-1methylpyrrolidinium chloride and dicyanamide had no significant impact on the charge scaling factor. 59 In contrast, Dommert et al. showed that the charge distributions obtained from fitting gas-phase first-principles data are significantly (qualitatively) different from those obtained from fitting the data from condensed phase DFT calculations. 48 In the absence of explicit polarization, the scaling of ionic charges is the most effective option which has been employed for the adjustment of various force fields for a wide range of ionic systems. Below we list several popular strategies for including polarization in a mean-field sense by modifying Coulomb interactions. The exact strategies for how to implement the scaling also vary: (1) Chaban has uniformly scaled the atomic charges by different factors to match various experimental thermodynamic properties 60 such as density, viscosity, conductivity, and heat of vaporization. Although these properties have no obvious correlation, a uniform charge scaling factor of 0.7−0.8 was sufficient to reproduce these properties in nonpolarizable simulations. Interestingly, this scaling factor is quite close to the inverse of the refractive index.
(2) Schroder reported on a linear correlation between the charge scaling factor of 1-ethyl-3-methylimidazolium trifluoromethanesulfonate [C 2 mim][OTf] and the polarizability on the basis of the effective Coulomb energy. 45 (3) Muller-Plathe and co-workers used an effective dielectric constant ϵ ∞ of 1.8 to scale down the Coulomb interaction. 61 This ansatz corresponds to an electronic continuum correction (ECC) or dielectric continuum model to reproduce the effect of the induced dipoles. 45,62,63 Following this model, one may use the inverse of the refractive index n 1 1 D = ϵ ∞ as a charge scaling factor (i.e., yielding a scaling factor of 0.75). 45 In addition to the ECC correction that weakens the ion− solvent interaction, one needs to decrease the ion− solvent repulsion parameters to avoid unrealistically high solvent−cation coordination numbers and overestimation of the ion−solvent packing. 64−66 Examination of 29 combinations of nonpolarizable ion−water force fields based on SPC/E, TIP4P, and TIP4Pew water models for the LiCl/water system at four different salt concentrations showed that none of the investigated models yield satisfactory results for all tested properties such as density, static dielectric constants, self-diffusion coefficients, and structure factor, but the force field using the ECC correction resulted in by far better prediction of transport properties. 46 (4) Unlike the ECC correction that often results in a weaker (compared to the full charge model) ion−solvent binding for solvent molecules in the first coordination shell of the ion, a mean-field polarizable model was proposed. This model added an effective polarizable term to enhance the ion−solvent attraction at short distances for solvent molecules within the first ion solvation shell. The effective polarization term was scaled to zero beyond the first coordination shell using a distance dependent dielectric constant approach. Thus, it effectively accounted for the increased dipole for the solvents directly coordinating small ions (e.g., Li + ) as has been observed in DFT-based MD simulations. 67 It roughly took into account the ion− solvent polarization screening due to electronic polarizability within the first solvation shell and due to dipole (and multipole) contributions beyond the first coordination shell. 68,69 Such an approach, however, overestimated the size of the first ion solvation shell and resulted in the slower ion dynamics compared to MD simulations with an explicit inclusion of polarization. 70 (5) Fileti and Chaban recommended that prior to the determination of partial charges the structure should be derived from DFT functionals including dispersion, e.g., ωB97XD, 71,72 followed by partial charge assignment via the electrostatic potential and Møller−Plesset secondorder perturbation theory to be more compatible with CHARMM force fields. (6) Schmidt et al. used the Blochl method to fit ionic charge distributions based on DFT data for the condensed phase (bulk) conditions. 58 Similarly, Mondal and Balasubramanian used crystalline and liquid phase DFT data to extract reduced charges for imidazolium-based ILs, thus implicitly accounting for polarization. 73,74 They also pointed out that utilization of different charge fitting schemes does not result in a consistent set of charges. A more extended discussion of earlier charge scaling strategies and approaches can be found in ref 48. However, such an approach might be to some extent counterintuitive for organic solvent electrolytes because DFT calculations revealed that the solvent dipole moment is higher for the molecules coordinating small cations than in a bulk solvent. Thus, the solvent−cation interactions in the first solvation shell should be stronger and not weaker due to polarization effects. 67 Chen et al. 75 compared three different nonpolarizable force fields commonly used for ILs for simulations of [C 4 mim] [Cl] and its mixtures with ethanol. For neat IL with nonscaled charges, i.e., ion q = ±1.0e, the self-diffusion coefficients, conductivity, and viscosity were about 2 orders of magnitude off from experimental values. Simulations of [C 4 mim][Cl] with scaled charges (ion q = ±0.85e) showed noticeable improvement of transport property predictions; however, they were still a factor of 5−20 off from experiments. For the [C 4 mim] [Cl] mixtures with ethanol the accuracy of predicted thermodynamic properties, such as the ethanol activity coefficient was found to be even more complicated. At dilute ethanol concentrations, the force fields with full ionic charges yielded values in reasonable agreement with experimental correlations, while simulations with scaled charges increased the activity coefficient by about a factor of 2. On the other hand, at high ethanol concentrations the trend reversed, showing a better description from simulations using the scaled charges while simulations with full charges overestimated the solubility of ethanol compared to experiments. Therefore, if someone is interested in investigating Chemical Reviews this IL/solvent mixture in the entire composition range, the selection of an appropriate force field/model becomes a taunting task.
The lack of transferability of nonpolarizable force fields with scaled charges specifically tuned to describe one system (e.g., pure IL) to mixtures of several ILs, IL with solvents or polymers, or organic solvent electrolytes, is one of the major disadvantages and limitations of the scaled charge approach. For example, if one mixes two ionic liquids sharing the same cationic species but with different charge scaling factors for the pure systems, the question arises which charge scaling factor should be applied to the cations? If one sticks to the factors for the pure systems, cations of the very same type may have different total charges which may lead to spurious artifacts in the simulations. Similarly, the alternate possibility to determine the uniform cationic scaling factor as a function of the mole fraction is not satisfactory as cations and the two anionic species can have different scaling factors. Furthermore, although the anionic scaling factors can be kept constant when changing the mixture composition, the cationic scaling factor has to be re-evaluated.
Choi and Yethiraj 76 reported that the scaled charge model for [C 4 mim][BF 4 ] failed to predict the phase separation of this IL with poly(ethylene oxide) chains at any temperatures, which is in disagreement with experiments. It was argued that this failure is a generic phenomenon, primarily due to significant and artificial underestimation of IL cohesive energy in simulations with scaled charges, therefore leading to poor predictions of phase behavior in IL mixtures. 77 This conclusion is further supported by recent work of McDaniel that states "Due to the important contribution of polarization, we f ind that non-polarizable force f ields qualitatively fail to predict mixing of ionic liquids with low dielectric solvents, predicting phase-separation instead!" after investigating the phase behavior of [C 4 mim][BF 4 ] mixtures with 1,2-dichloroethane, acetone, acetonitrile, and water using polarizable and nonpolarizable force fields. 78 McDaniel and Yethiraj also found that scaling of charges cannot remedy the artificially enhanced long-range ion−ion correlations in ionic liquids. 79 They demonstrated that this artifact is apparent to simulations with any nonpolarizable force field due to inability to capture the infinite frequency dielectric response ε ∞ (and hence electrostatic screening) by nonpolarizable models. The authors showed that the long-range electrostatic interactions, usually handled by one of the Ewald-type summation techniques, are fundamentally altered by electronic polarization. Polarizable simulations predict ε ∞ ≈ 2, while all nonpolarizable models (whether with scaled ionic charges or not) give ε ∞ = 1, therefore leading to overestimation of electrostatic interactions in the small wavevector limit. This work clearly illustrates the fundamental difference between polarizable and nonpolarizable models that attempt to approximate the influence of polarization.
In this review, we discuss the pros and cons of different treatments of polarization ranging from the mean-field approaches to the most popular explicit polarization models in MD simulations of highly ionic materials. Most of molecular simulation softwares now have an option to include induced polarization interactions. Therefore, each user (whether a modeling expert or a novice user) has to make choices deciding which model to select for his/her particular system of interest. We believe this review will serve as a practical guide both to expert and nonexpert users of molecular simulations.
The goal of this review is, first, to briefly introduce the basic concepts of induced polarization interactions in molecular systems and to describe several strategies for implementing these effects into molecular simulations such as fluctuating charge models, classical Drude oscillator, and point dipoles. For each method, we discuss the advantages and disadvantages as well as key assumptions and adjustable parameters (the latter are often taken as default universal parameters, which is not always the case). Furthermore, we compare the outcomes using one or the other polarizable model and discuss the origin of the observed discrepancies.
Then, we review the methods/strategies for the parametrization of polarizable models, specifically focusing on extracting atomic polarizabilities (the key adjustable parameter in polarizable force fields) and universal scaling behavior that has been discovered recently. Next, we compare simulation predictions obtained using polarizable and nonpolarizable models. We discuss the underlying physics that each approach includes or ignores, implications for implementation and computational efficiency, accuracy, and transferability. These comparisons are made for several classes of bulk systems that include: dilute and concentrated electrolytes with the focus on emerging battery applications, bulk ILs, and electrolytes at charged surfaces, where the isotropic approximations are no longer valid and explicit inclusion of polarization of electrolyte and the electrode can be important.

Fluctuating Charge Model
Polarizability refers to the response of the electron density due to a local electric field. In classical MD simulations, the electron density is accounted by the partial charges q iβ of atom β in molecule i of the system, allowing them to change during the simulation, i.e., q iβ (t), mimicking polarizability effects. 80−82 This approach has a big advantage that no new interaction types have to be implemented as polarizable and permanent electrostatic interactions are handled using Coulomb's law with the same monopoles. This method is referred to as the fluctuating charge, charge equilibration, or chemical potential equilibration model. While it has been successfully implemented in a number of programs, one has to be careful because a possible overestimation of the polarizability as significant charge flows within the molecule (or even between them) may happen at small energetic costs. 81 The energy U(q iβ (t)) required to create a charge q iβ (t) on an isolated atom β of molecule i can be described by a Taylor series with the amount of charge flow based on the Mulliken electronegativity χ iβ 0 and the hardness J iβ 0 . Going beyond the quadratic term possibly increases the accuracy of the model but causes problems for treating the charge dynamics as discussed below. However, in molecules, this electrostatic energy is augmented by the interaction with the partial charges of the other atoms γ of the same molecule i or another molecule j: Here, {q(t)} and {r⃗ (t)} represent the sets of partial charges and coordinates of all atoms at time t. The second-order coefficient J iβjγ (r⃗ iβjγ (t)) depends on the distance between the atoms iβ and jγ and should become 1 / 4 πϵ 0 r iβjγ at very long distances. At shorter distances, the screening of the electrostatic interaction between the atoms due to delocalized charge distributions ρ(r⃗ (t)) can be computed by the integral Using these mixed hardness parameters ensures the correct limiting behavior as 1/r iβjγ (t) in eq 2.1.2 for separations greater than 2.5 Å. Consequently, eq 2.1.4 applies to induced charge interaction of atoms connected by bonds, angles, and dihedrals.
The molecular polarizability 3 × 3 tensor α⃡ (t) can be evaluated from the inverse of the n × n hardness matrix containing the elements J iβjγ (t): where Δr⃗ (t) are the atomic coordinates relative to the center of geometry of the n atoms under investigation. 2.1.1. Charge Flux and Electronegativity Equalization. The current partial charges q iβ (t) on each atom β of molecule i are obtained by minimizing the electrostatic energy using the Lagrange multiplier method with one of the following charge restraining conditions f(q): (1) Bond charges b βγ between atoms β and γ connected by a covalent bond cancel each other: 84−87 The unambiguous assignment of bond charges b βγ is always possible for neutral molecules but needs additional rules for distributing the net charge in the case of ions. This has the considerable advantage that partial charges cannot be accumulated at a particular site of the molecule. This might happen in ionic systems as the redistribution stops once the charge has arrived at the site with the closest distance to the counterion, as the electrostatic potential between unlike charges is very steep at short distances.
(2) The molecular charge q i is constant: The molecular charge can be set to zero for neutral species, +1e and −1e for the monovalent ionic liquid cations and anions, respectively, or to subinteger values as inspired by the work of Morrow and Maginn. 6 This way, intermolecular charge transfer processes are prevented. This is the most common option, in particular if only one ionic species is made polarizable. 88,89 (3) Charge transfer is allowed between an ion pair: 90 However, the complete assignment which cation i belongs to which anion j is ambiguous as a particular anion may have the shortest distance from several cations. Nevertheless, temporary charge transfer between cations and anions can be realized this way.
(4) The total charge of the system is zero which should always be the case, even if the charge flow is restricted otherwise. This option guarantees free charge flow in the system. In the case of pure ionic liquids, this may jeopardize the stability of the polarizable MD simulation. The energy gradients are electronegativities χ iβ (t) of the corresponding atoms by Mulliken's definition 91 12) and result in Because eq 2.1.13 equals a constant λ for all atoms sharing the respective charge condition described above, their electronegativities are equal as well. In other words, the electronegativity χ iβ (t) of each atom β in molecule i equals the average value χ i (t) of the molecule i if the molecular charge is constant (see second charge condition eq 2.1.9). If the charge flow is allowed between all atoms in a system, i.e., only the total charge of the system is zero, all atomic electronegativities χ iβ (t) at a particular time have the very same value. The fluctuating charge model is often also called electronegativity equalization method or chemical potential equalization method because the chemical potential μ iβ (t) equals the negative electronegativity χ iβ (t): Using the Taylor expansion of the electrostatic energy up to the second term (see eq 2.1.1) and the minimization condition in eq 2.1.13 results in a coupled set of (n −1) linear independent equations of charges and the charge condition f(q). Solving these coupled equations by Cramer's rule yields all partial charges q iβ (t).
The partial charge solution can also be obtained by means of an extended Lagrangian dynamics method. Here, a fictitious uniform mass, m q , is assigned to each charge. Time evolution of the charges are computed by Newton's equation of motion with the average electronegativity χ i (t) of the molecule i in case of second charge condition eq 2.1.9 As a consequence, eq 2.1.15 describes forces experienced by each atom β driving the atomic electronegativity χ iβ (t) toward the molecular average χ i (t). In this view, this approach mimics the process of moving electrons from atoms with low electronegativity (= high chemical potential) to highly electronegative atoms having a lower chemical potential.  88,90 However, most of these simulations are performed with self-written MD codes as only LAMMPS, earlier versions of CHARMM, 83 and nonofficial versions of GROMACS have implemented a fluctuating charge algorithm. This explains the far less common use of these polarizable method compared to Drude oscillators or induced point dipoles.
Of course, the atomic electronegativity χ iβ (t) of atom β in a molecule i differs from the value of the isolated atom χ iβ 0 as the chemical environment and, in particular, neighboring atoms directly bonded to atom β influence the behavior of the charge flow. Hence, the parameter values χ iβ 0 and J iβ 0 for an atom type in eq 2.1.2 are usually not computed from the ionization potential and electron affinity as proposed by Mulliken 91 but are subject to optimization routines. 82,100 However, they may reveal some insights in the electronic structure of that atom, e.g., the oxidation state.
In Figure 3, the fluctuating charge parameters for methylimidazolium C n mim from ref 88 are displayed. The anions were made nonpolarizable in the reported systems and the charge flow was restricted to each cation corresponding to the molecular charge restraining condition (eq 2.1.9). Using these parameters, the partial charges of the cations resemble the Mulliken charges of these molecules. 88 However, the parameters do not seem to be transferable as prolonging the alkyl chain has a significant impact on χ 0 and J 0 . In particular, the electronegativity of the imidazolium nitrogens jumps to higher values from C 2 mim to C 4 mim. The values for the acidic ring carbon C2 increases with increasing alkyl chain length. This is also true for the carbon adjacent to the nitrogens. Interestingly, the electronegativity of the C5 carbon for C 2 mim is much higher than for C 1 mim and C 4 mim. The charge flow within the molecules seems to be strong as the electronegativity of the terminal methyl carbon C6 changes by a factor of 25 between C 1 mim and C 4 mim. The electronegativity of the alkyl carbons increases with the distance from the imidazolium ring. Quite counterintuitively, the fluctuating charge electronegativities do not correlate with the acidity of the atoms because the alkyl carbons possess higher values than the ring carbons. The acidity of the C2 carbon should be highest, but the corresponding electronegativity has the lowest value. The reported hardness values J 0 show less variation among the atom types.
Standard charge equalization schemes work reasonably well for structures of neutral molecules close to equilibrium but may have problems for nonequilibrium situations or charged systems, 81,82,99 resulting in nonphysical charge flows. For example, increasing the bond length of a sodium chloride ion pair should result in a vanishing partial charge of the participating atoms at infinite distance. However, as visible in Figure 4, the partial charge of the sodium as a function of the distance to the chlorine does not approach zero (red line) for the standard equalization scheme although Mulliken charges of CAS(8,5)/3-21G calculations (blue dashed line) went to zero  The fluctuating charge formalism could also be used in conjunction with the Drude or atom polarizable models to describe the intermolecular charge transfer (CT). 102 It is a complementary approach to account for CT via the explicit introduction of additional terms in the potential energy function, 103 like it is done in the sum of intermolecular fragment ab initio (SIBFA) model. 104,105 In the fluctuating charge model of CT, parameters in eq 2.1.12 model only the CT part. The negative sign of the Mulliken electronegativity, χ iβ 0 , represents the tendency of an atom to attract electrons in intermolecular CT. The hardness represents the atoms' resistance to losing electrons in intermolecular CT. There is also a maximum amount of CT for each pair that is dependent on the distance and decays to zero at large separations in order to limit the unphysical low-energy long-range transfer. 102 Here, the first, second, and third exponential term determines the single, double, and triple bond order, respectively. The parameters c 1···6 and the equilibrium distance r σ , r π , and r ππ have been modeled to agree with corresponding quantum-chemical results at a distance of r iβiγ (t) between the atoms β and γ. Angle and torsional potentials depend on the respective bond orders.
The Coulomb interaction between two atoms is shielded at shorter distances with the shielding parameter s iβjγ , which is the arithmetic 110 or geometric mean 111 of the respective atomic contributions. The time-dependent partial charges q iβ (t) and q jγ (t) are obtained by an electronegativity equilibration scheme as described in section 2.1.1. When the charge flow between molecules is not restricted, the spurious long-range charge transfer occurs in the standard ReaxFF. To minimize this artifact, the atom-condensed Kohn− Sham density functional (ACKS2) charge calculation scheme was developed by adding quadratic energy terms to eq 2.1.2 in order to control the range over which the charge is allowed to delocalize. 112 The reference charges were introduced in ACKS2 in order to distinguish between ions, neutral molecules, and to correctly describe limiting charge transfer. In the ACKS2, every atom must have a reference charge and the total charge of an isolated molecule is always equal to the sum of these reference charges. These rules are essential for the correct dissociation limits and linear response properties of ACKS2. In addition, the hardness parameters are dependent on the interatomic distance controlling the change of charge transfer, similar in spirit to the approach by Martinez and co-workers discussed above and shown in Figure 4.

Classical Drude Oscillator Model
Another technique to allow for a polarizable response of molecules to an electric field is the classical Drude oscillator model, which is also known as the "charge-on-a-spring" or "shell" model. 113−115 In contrast to the fluctuating charge model, 82,115 the Drude oscillator model does not modify partial charges due to changes in the local electric field during the simulation, but adds additional particles (the oscillators) that mimic physical dipoles on each polarizable atom to model the corresponding distortion of the electron density. The first particle comprising the dipole is located at the position of the nuclei of the polarizable atom β of molecule i and usually its partial charge, −q D that contributes to the dipole of the atom is merged with the corresponding partial atomic charge q iβ of the polarizable atom. The second, mobile Drude particle or oscillator carries a partial charge of q D and is tethered by a harmonic spring to the atomic nuclei and moves around the polarizable atom as depicted in Figure 5. Because the Drude charges q D are set to negative values, the mobile Drude particle is interpreted as a representation of the electron cloud of the polarizable atom. This picture should not be taken literally because the Drude pair is a simple method to introduce a dipole on each atom that can be handled with the same or similar algorithms for atoms in molecular dynamics simulations, e.g., reaction fields, Ewald summations, and particle mesh Ewald techniques 114−117 as only additional charged particles are introduced. Furthermore, QM/MM with polarizable forces are easier to implement with Drude oscillators 118 than induced point dipoles. As the Drude particles may point in any direction in three-dimensional space, they are not restricted to the dimensionality of the underlying molecule. For example, the polarizability of a fluctuating charge model of imidazole is more or less restricted to the two-dimensional plane in which the atoms lie. However, as the Drude oscillators may point below and above the imidazole plane, the polarizability is truly threedimensional. In addition, while recent studies have shown that the Drude oscillator model may be considered equivalent to the induced dipole model, 119 it has the advantage that a van der Waals term may be included on the oscillator, thereby offering steric effects associated with distortion of the electron cloud. 82,120,121 Each Drude pair results in a physical atomic dipole μ⃗ iβ ind shown as transparent arrows in Figure 5: The displacement d ⃗ iβ (t) of the mobile Drude particle from the position of the polarizable nucleus r⃗ iβ (t) originates from the balance of electrostatic forces characterized by the local electric field E ⃗ iβ (t) and forces of the harmonic spring between the nucleus and Drude particle with the force constant k iβ D . Usually, these displacements are quite short; in fact, they are largely exaggerated in Figure 5, with typical lengths |d ⃗ iβ (t)| of less than 0.1 Å. 122 A hardwall restraint has been introduced to enforce an upper limit for this distance, e.g., |d ⃗ iβ (t)| < 0.2 Å, 123,124 to prevent instabilities in the simulation from overpolarization of the Drude particles.
In practice, increasing the Drude charge q iβ D while greatly increasing the force constant lowers the average displacement of the mobile Drude particle to keep the induced dipole μ⃗ iβ ind (t) almost constant. However, this holds only true for Drude charges |q iβ D | > 1.0e as visible in Figure 6 for 1-ethyl-3methylimidazolium C 2 mim + . 122 Van Gunsteren and co-workers used q D = −8.0e to yield a high force constant on the nucleus− Drude pair, thereby keeping the displacement at a minimum. 117 This approach, versus smaller Drude charges and associated force constants, yields physical dipoles that more closely mimic a pure induced dipole model. However, very high values for the Drude charge may cause problems for the dynamics of the mobile Drude particle as discussed below.
The energy U pol (t) due to polarizable forces consists of three contributions: 45,122 The harmonic potentials of the Drude springs the Coulombic interaction between the nucleus and Drude particle charges associated with the dipoles on the polarizable atoms (see Figure 5)    The prefactor 4πϵ 0 includes the vacuum permittivity of 8.85 × 10 −12 As/Vm. The dipole−dipole tensor T r t r r r t ( ( )) 1 4 depends on distance vector r⃗ iβjγ (t) = r⃗ jγ (t) − r⃗ iβ (t) from the atom β of molecule i to atom γ of molecule j. 125 ... − + + at small Drude distances d ⃗ iβ (t) and d ⃗ jγ (t) and show that the electrostatic interactions of the Drude pair model resembles those of mathematical induced dipoles described in the next chapter concerning induced point dipoles. Although the last equations suggest a more or less complete analogy between Drude oscillators and induced-point dipoles differing only in technical details, the polarizable model of Drude particles is based on charge−charge interactions and offers additional features, e.g., modeling steric effects of the distortion of the electron cloud by introducing van der Waals interactions with the mobile Drude particle 120 or QM/MM mixed approaches. 126 130 DLPOLY, 131 ESPResSo, 132 GROMACS, 123 LAMMPS, 133 OPenMM, 134 and NAMD. 116 These implementations may be performed in various fashions, as described below.
(1) All Drude charges are set to a uniform value, i.e., q iβ D = q D . 45,122,135−139 Consequently, the harmonic force constant k iβ D depends on the isotropic polarizability α iβ of the respective atom β of molecule i: High polarizabilities weaken the spring of the respective mobile Drude particle allowing for larger displacements and consequently higher induced dipoles. Figure 6 shows that uniform polarizabilities α iβ and hence uniform k iβ D for carbons do not necessarily result in comparable induced dipole moment |μ⃗ ind | as these dipole depends on the local field E ⃗ iβ (t) which is stronger for C2 and C8 because the anions approach the imidazolium cations from this direction.
(2) All force constants for the harmonic Drude spring are set to a uniform value, i.e., k iβ D = k D . 140 and are depicted in Figure 7 for common atoms in ionic liquids. A harmonic force constant k D of 1000 kcal mol −1 Å −2 ensures that the Drude charges are above the limit to yield induced dipoles independent of the actual value of q iβ D . For example, |q iβ D | for nitrogens and carbons is approximately 2.0e and hence the induced dipole of those atoms should not vary in the simulations if one increases k D .
(3) In the polarizable, coarse-grained MARTINI force field 143,144 two mobile Drude particles carrying opposite charges (q ≅ ±0.46e) are attached to the polarizable atom. The bond length of each mobile Drude particle to the atom is fixed to l = 1.4 Å and the angle between the two bonds is modeled by an harmonic potential (k θ = 4.2 kJ mol −1 rad −2 ,θ eq = 0°). The Coulomb and Lennard-Jones interaction between the two Drude particles are disregarded. Consequently, a zero angle represents the nonpolarized situation. The maximum induced dipole is μ ind = 2 l q D . 2.2.2. Dynamics of Drude Particles. After the movement of all atoms in the simulation at each MD step according to the potentials given in the force field, the positions of the mobile Drude particles are determined by the electric field generated by the remaining atomic dipoles and the permanent charges of the atoms. This corresponds to a relaxation of the electronic degrees of freedom immediately upon any change in the nuclear configuration in analogy to the Born−Oppenheimer approximation. A self-consistent field approach may be applied to find the minimum total electrostatic energy for the Drude particles according to for all mobile Drude particles of polarizable atoms β of molecules i. 146,147 However, this procedure is computationally demanding due to multiple force evaluations at each minimization step. As a Car−Parrinello alternative to the SCF minimization, an extended Lagrangian approach was developed by Lamoureux and Roux based on a dual-thermostat approach. The method has been implemented in CHARMM, 130 GROMACS, 123 LAMMPS, 133 OpenMM, 134 and NAMD. 116 Because of the speedup of computational demands over the SCF methods, the

Chemical Reviews
Review dual-thermostat extended Lagrangian is the standard for MD simulations using the Drude force field. 114,123 In the extended Lagrangian approach, the first thermostat keeps the atoms at the desired simulation temperature T. The thermostated forces 114,123,133 act on the center-of-mass of the atomic nucleus iβ and its corresponding mobile Drude particle. Here, r⃗ iβ (t) and r⃗ iβ D (t) = r⃗ iβ (t) + d ⃗ iβ (t) are the positions of the atomic nucleus and the mobile Drude particle, respectively.
The second thermostat serves to dynamically maintain the nucleus−Drude oscillator at a very low temperature to a polarization energy U pol (t) close to the SCF value "on the fly". 130 In practice, each mobile Drude particle is assigned a small mass m D that is subtracted from the mass of the polarizable atom to keep the total atomic mass of the polarizable atom constant.
The implementation of the model using a mass m D of 0.4 amu and a force constant on the nucleus−Drude pairs of 500 kcal mol −1 Å −2 allows for an integration time step of 1.0 fs in MD simulations. 114,124,130,141,148−152 Polarizable simulations of ionic liquids have been carried out with a Drude mass m D of 0.1 amu, 45,122,138,139 although using too small values for Drude mass might also affect the integration time step that can be used. Because of the finite mass ascribed to the mobile Drude particles, treating the polarizability of light atoms such as hydrogens using an extended Lagrangian may also encounter some challenges, although this is not an issue if the model is propagated via a SCF approach.

Induced Point Dipoles
The polarization response of a molecule to an electric field can also be modeled via the induction of atomic point dipoles. In contrast to the Drude oscillators, these induced point dipoles are mathematical dipoles without additional particles having a mass. Hence, in principle, they can be placed anywhere in the molecule. Typically, induced point dipoles are placed at the center of atoms but can also be added to the massless force centers situated off the atomic sites in order to better represent anisotropy and spatial distribution of the polarization response.
In the latter case, the forces then are transferred to the basis of three atoms defining the position of the massless center using chain rule differentiation. Induced point dipoles are implemented in AMBER, 153  ind at the force center i is proportional to the total electric field E ⃗ iβ at this point that is composed of the electric field due to permanent charges and multipoles E ⃗ iβ q and a field due to induced dipoles E ⃗ iβ The prefactor αî β is a 3 × 3 atomic polarizability tensor. If the polarizability tensor is isotropic (i.e., diagonal and α iβ xx = α iβ yy = α iβ zz ), it can be replaced by a scalar value α iβ which is one-third of the trace of αî β . Atomic polarizability terms can be straightforwardly included in polarization at both atom and off-atom massless force centers. Anisotropic polarizabilities lead to a torque on the force center that needs to be distributed among the bonded atoms. In practice, models mostly use isotropic polarizabilities but some force fields that actually embody anisotropic dipole polarizabilities enable a more detailed representation of the polarization response. Here we only give an illustration of the polarization equations based on models limited to point charge electrostatics, but it can be extended to the more general case of multipole interactions, leading to more complex equations that have been shown to be fully extended up to periodic boundary conditions in Ewald summations. 160,161 The polarization energy U pol = U μq + U μμ + U self can be decomposed into contributions from the interaction of the induced dipoles with permanent charges (dipole−field interaction, cf U qD (t) in eq 2.2.7) and the self-polarization energy (see eq 2.2.3) The interaction energies U μq (t) and U μμ (t) stem from the interaction with the external electric field, whereas the last contribution represents the required work to create the induced dipoles. 82,147,162 Equation 2.3.5 defines a many-body problem for finding the induced dipoles i k j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z of the N atoms. Consequently, the polarization energy U pol is Using a variational method, e.g. the iterative atomic multipole optimized energetics for biomolecular simulation (AMOEBA) approach, 164 the induced dipoles are obtained by minimizing the residual The dominant computational cost is the repeated evaluation of Rμ⃗ ind during the self-consistent field computation. Of course, one way to reduce the computational cost is to use a direct polarization method (iAMOEBA), 165 where the coupling (= off diagonal) terms are neglected. This can be rationalized by the fact that the interaction between close induced dipoles will be damped (see section 2.4.3). If the interacting induced dipoles are more distant, the corresponding interaction approaches zero as a function of r −3 , making the respective off-diagonal elements in R very small. Hence, the computational bottleneck of the direct method is shifted to the evaluation of the electric fields and their derivatives with respect to the coordinates. 165 However, despite the speed gain, for most systems including ionic ones, inclusion of fully converged induced dipoles matter and one would like to stick to the full resolution of the polarization equations. Fortunately, such problems are mathematically well-defined and various highperformance strategies can be then employed to speed up the process. We will discuss these advanced techniques in section 2.4.
To conclude, it is also important to highlight that due to its self-consistent nature, the point dipole model is well-suited for hybrid QM/MM simulations 166 where the induced dipoles can be fully coupled to the electronic density, i.e., a full selfconsistent relaxation of both the AMOEBA induced dipoles and the DFT electron density at each MD step is then possible allowing for embedded DFT Born−Oppenheimer/AMOEBA simulations. 167 In the same spirit, force fields such as AMOEBA can be self-consistently coupled to polarizable solvation methods to perform MD simulations. 168 2.3.2. Comparison to Drude Oscillator Model. Both the induced point dipole and the Drude oscillator approaches try to model the electronic degrees of freedom by an induced dipole. In a recent publication, 135 Schroder and co-workers checked in the case of 1-ethyl-3-methylimidazolium trifluoromethylsulfonate [C 2 mim][OTf] if both approaches result in similar dynamic properties. This is particularly important as the Drude oscillator approach using the extended Lagrangian does not allow for polarizable hydrogens. Consequently, the effect of merging hydrogen polarizabilities with the polarizability of the atom to which they are attached (implicit H) versus the full atomistic polarizable representation (explicit H) was tested and compared to the induced point dipoles in AMBER. 153 As visible in Figure 8, the radial distribution functions g 000 (r) and g 110 (r) for cation−cation, cation−anion, and anion−anion of the implicit (violet) and explicit (orange) hydrogen polarization are in excellent agreement. The radial distribution function g 110 (r) weights the respective g 000 (r) with the cosine of the angle between the total molecular dipoles of the species. Negative values (e.g., for the cation−cation orientation) indicate that antiparallel alignment of the dipoles is preferred. However,  also for dynamic properties such as the rotational relaxation constant (via the relaxation of the dipole−dipole autocorrelation function) as well as the diffusion coefficients (via the meansquare displacement), the perfect agreement between implicit and explicit hydrogen polarizability using induced point dipoles (IPD) holds. Neglecting the hydrogen polarizability (IPD no H (dashed black line)) in Figure 8 results in lower molecular polarizabilities and hence lower translational and rotational mobility of the ions. Nonpolarizable simulations yield much lower diffusion coefficients and longer rotational relaxation times. Because the force field used in CHARMM (DRU) and in AMBER (IPD) is the same, results from the nonpolarizable simulations coincide (blue and green lines). However, small discrepancies are observed between the polarizable Drude oscillators in CHARMM (DRU no H (red lines)) and induced point dipoles in AMBER (IPD no H (black dashed line)) for the transport properties. This becomes more obvious for the implicit hydrogen polarization (gray and violet lines in Figure  8). This can be due to the different thermostats applied to the polarizable model. In CHARMM, the temperature of the mobile Drude particles was close to 1 K, whereas the induced point dipoles in AMBER showed higher temperatures up to 40 K. The discrepancies between DRU and IPD is more pronounced for the anions which possess atoms with higher polarizabilities, in particular, the sulfur atoms are highly polarizable. For these highly polarizable atoms, the representation of the induced dipole by a pair of charged Drude particles may not be accurate anymore. However, the biggest drawback of nonpolarizable hydrogens in polarizable simulations using Drude oscillators and Lagrangian thermostats can be solved by adding hydrogen polarizabilities to the atoms they are attached to.
To conclude on such model comparisons, it is worth noting that an automatic strategy allowing mapping of the Drude polarizable force field onto a multipole and induced dipoles model is currently developed to enable the direct use of Drude models into induced dipoles codes. 119

Beyond Induced Dipoles
At this point, the following question can then be raised: are induced dipole models enough to deal with many-body interactions in complex systems? Indeed, one can always try to compare the induction energy values extracted from modern ab initio energy decomposition analysis such as symmetry adapted perturbation theory (SAPT) 169 to force field estimates, and as a matter of fact, usually they do not match. This comes from various reasons, the first one being the definition itself of the polarization energy. In ab initio theory, the induction energy is more general than the classical polarization energy obtained from a point dipole approximation. Indeed, such a term appears at second order in a Rayleigh−Schrodinger perturbation expansion of the total intermolecular interactions. 125,169 Physically, the induction contribution is the energy of interaction of the permanent multipole moments of one molecule interacting with the induced multipole moments of another. Force field models do approximate the long-range behavior of the induction energy (i.e., the polarization) but usually fail to give a good approximation at the short-range which embodies both charge transfer and other nonclassical effects. Various solutions can be found to overcome such difficulties.
2.4.1. Higher-Order Induced Moments. As a matter of fact, the long-range, i.e., classical, part of the induction could be modeled by expanding electrostatics to high-rank multipoles to compute accurately permanent electric fields and by using higher-order polarizabilities. Higher-order terms for the polarization computation appear as the induced point dipole model is only the truncation of the total response. Indeed, in practice, if a dipole moment is by far the largest contribution to the response, it is induced along a long series of higher-order electric induced moments by introducing dependence not only to the electric field but also to the field gradients and so on. 125 Effects of the higher-order electric induced moments are discussed since the 1960s 170 and include dipole−quadrupole, quadrupole−quadrupole polarizabilities as well as the first hyperpolarizability. They have competing effects but lead to a modulation of the final induced dipoles values. 170 Over all of these higher-order induced moments, the induced quadrupoles were reported to have the most noticeable effects 171 and lead to a non-negligible polarization contribution in the case of metal and heavy metal cations. 172 Net benefits of a generalized inclusion of higherorder polarizabilities are yet to be demonstrated in real-life simulations, but accurate approaches offering direct evaluation of such quantities are now available. 170,173−175 However, the corresponding quantities are not available in common MD programs.
As a first, cheap remedy to implement anisotropic field gradients in the Drude model, anisotropic polarizabilities can be defined in CHARMM, 113,114,176,177 i.e., the induced dipole μ⃗ iβ ind (t) may point in a different direction than the local electric i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z The isotropic force constant k iβ D is expanded to a 3 × 3 tensor, with zero off-diagonal and the diagonal elements k xx D , k yy D , and k zz D allowing for anisotropic displacement with respect to an intramolecular reference frame located at that atom (see Figure  5). This feature was introduced in CHARMMs Drude 2013 force field to more accurately describe hydrogen-bond acceptors 177 and may be in particular valuable for protic ionic liquids. Furthermore, mobile Drude particles are allowed to have Lennard-Jones parameters in CHARMM 120 to model the anisotropic van der Waals surface of atoms. Anisotropic distributed polarizabilities are, by definition, also used in point induced polarizable force field such as the sum of interaction between fragment ab initio (SIBFA). 175 When distributed on nonatomic centers such as lone pairs, anisotropic polarizabilities were shown to provide a closer agreement with the ab initio reference data, both in terms of polarization energy and in terms of dipole moment. 178 2.4.2. Modeling Short-Range Induction. As we discussed, point dipole models usually encounter difficulties to model the short-range contribution of induction. Indeed, when Chemical Reviews Since the relative importance of short-range inductions differs between variational energy decomposition analysis 179,180 and perturbation approaches 181,182 such as SAPT, its exact magnitude is still a matter of debate. 125 It was furthermore shown that the short-range induction should not be included within classical polarization contribution and should preferably be either incorporated into the pairwise van der Waals contribution or treated explicitly. 183 If explicitly taken into account, charge transfer is usually modeled by a simple exponential function 184 according to its known exponential decay. 125 However, as part of the induction, charge transfer exhibits a many-body behavior 105,125 and many-body force field approximations exist in the framework of SIBFA 104,105,175 or effective fragment potential (EFP) 185 polarizable force fields and are particularly useful for the modeling of metal ions where donation and back-donation become non-negligible contributions. 175,185 The absence of the Pauli repulsion in the induced dipole model leads to potential polarization catastrophe at the short range. The earlier Drude oscillators and induced dipole models for water tended to overpolarize if the gas phase polarizability was used. This was due to lack of screening resulting in too high dielectric permittivity. A typical solution was to reduce the water model polarizability from 1.44 to 1.04 Å and scaling polarizabilities for other solvents as done in CHARMM. Such adjustments limit overpolarization and allow for the correct treatment of the dielectric constant in the liquid.

Thole Screening Functions.
When molecules closely approach each other, the molecular orbitals overlap and atom-based charge models no longer describe the electrostatic potential adequately. 186,187 Quite often, the pointcharge model predicts that the electrostatic potential continues to increase/decrease in the immediate vicinity of atoms, while QM calculations clearly indicate that the interaction should taper off or be damped. As a consequence, such inadequate electrostatic potentials often lead to overestimation of the many-body polarization energy at short distances (see eq 2.3.5). This shortcoming could be partially corrected by increasing repulsive interactions between force centers. However, this leads to shifting the energy from the many-body nonadditive polarization response to the two-body additive repulsion terms, thus degrading the description of the potential energy. Four approaches can be applied to improve the modeling of the electrostatic potential and field around a molecule or ion: (1) Simply add additional charges on off-atom force centers.
(2) Introduce a screening function s k (r) to correct for the penetration energy. 163 (3) Add atomic or bond centered dipoles, quadrupoles, etc. (4) Combine the first two approaches and use screened multipoles such as Gaussian multipoles. 186 Bearing analogy to the screening of partial charges, the screening functions between induced dipoles at short distances do not only improve the description of the interaction energies but are necessary to obtain convergence of the self-consistent interative process under certain conditions. The most popular approaches for damping the induced dipole−induced dipole interactions is (1) to either utilize the screened induced dipoles in conjunction to the distributed charges using Gaussian, Slater charge density distributions, 188 or point charges or (2) to combine Thole screening 163,189 with point charges and multipoles. While Thole originally examined numerous scaling functions, the exponential screening of the induced dipoles became the most widely adopted perhaps due to the similarity with Slater orbitals. In 1981, Thole 163 developed screening functions on the basis of smeared charges qĩ β of atom iβ instead of point charges q iβ using a dimensionless distance u r /( ) and the smearing function ρ k (u) (see Figure  9). The first derivative with respect to u of the screening function s k ′(u) is as depicted in Figure 9. In the following, the vector r⃗ is defined by r⃗ (t) = r⃗ iβ (t) − r⃗ jγ (t). This charge smearing has consequences for the local fields emerging from the charges (see eq 2.
and from other induced dipoles (see eq 2. as the corresponding electrostatic interaction tensors 115,162 T r t s r t r t r t ( ( )) ( ( )) 4 T r t s r t r t r t r t s r t are screened by the first and second derivative of the damping function s k (r(t)). Here, Îis a 3 × 3 unity matrix. The moderate effect of these functions on structure and dynamics of [C 2 mim][OTf] was reported in ref 190. The exponential Thole functions s 1 and s 2 seem to slightly reduce the molecular polarizability by shifting the distances between the ions but keeping mutual positions and orientations. This effect can be reduced for s 1 by using a higher default radius a. In fact, the default value is increased from 2.089 in AMBER to 2.600 in CHARMM. More recent versions of the CHARMM Drude model apply an atom-specific s 1 value making the term a parameter that my be optimized as part of the parametrization process. 191 As not all screening functions s k are available in standard MD programs like TINKER, CHARMM, DLPOLY, GROMACS, and LAMMPS, one may replace them by 1−2, 1−3 exclusions for the interaction of induced dipoles like already used for Coulomb and Lennard-Jones interactions. Here, the induced interaction between atoms connected by bonds and angles are not computed. This prevents an intramolecular "polarization catastrophe," and the effective molecular polarizability gets close to the sum of all atomic polarizabilities which might not be the case using the shape functions. However, intermolecular induced dipole interactions are not damped by the exclusions. Therefore, one has to be careful during the force field parametrization process as only large enough Lennard-Jones spheres prevent the intermolecular "polarization catastrophe". However, polarizing the Canongia-Lopes and Padua force field 192 as reported in refs 45, 122, 135−138, 190, and 193 caused no problem using the exclusions instead of the Thole screening functions.
2.4.4. Other Polarizability Damping Approaches. Nonlinear short-range induction effects are linked to orbital overlap, and can be modeled by damping models such as Tholetype models, or can be traced back to the value of the electric field as obtained from point charges or multipole classical approximations. Kunz and van Gunsteren 194 suggested an alternative approach by introducing nonlinear effects that would limit the maximum electric field in either Drude or induced dipole models as given by with the adjustable parameters p and E 0 . Here, α and E are the original polarizability and the electric field, respectively. The corresponding induced dipole of triflate [OTf] as a function of the electric field E is depicted in Figure 10 as dashed lines with p = 8 and E 0 = 2 V/Å. At this strong electric field, nonlinear effects are expected. 194 AIM (= atoms in molecules) calculations of Schroder and co-workers also show nonlinear behavior beyond 2 V/Å but instead of leveling off with increasing field strength, the  polarizabilities diverge. However, the underlying QM calculations start to have convergence problems and hence the actual values at these high field strengths should not be trusted too much. On the other hand, the evaluation of the polarizability α iβ is constant over a large electric field regime. In practice, they determine their polarizabilities at an electric field strength of 0.0008 au = 0.041 V/Å as indicated by the dashed line in Figure  10. 195−197 Discrepancies from the linear behavior can also be modeled by the hyperpolarizability βÊ In CHARMM, 176 additional potentials U hyper = k hyper (d − d 0 ) n may be applied to account for the nonlinear hyperpolarizability effects. Here, d is the distance of the Drude particle from the corresponding polarizable atom. The default value of the "equilibrium distance" d 0 is 0.2 Å. The potential is only calculated at distances larger than d 0 . Consequently, the force resulting from this potential is reducing the distance d and thereby weakening the induced dipole μ ind . The default value for the exponent n is four.
Models such as SIBFA 175 or EFP 198 also used damping functions of the electric field to maintain a linear regime in the polarization evaluation. Indeed, classical/multipolar electric fields clearly differ from their ab initio counterpart missing some screening effects. 199 Such behavior was noticed by Chelli et al. 200 They showed that nonlinear effects cannot be simply linked to the lack of hyperpolarization in the polarizable models but depend on the strength of the electric field. In contrast to the linear response at weak fields, strong fields enforce an intramolecular charge redistribution resulting in a nonlinear response. In other words, having some QM description of the fields is important and another strategy evolved starting with the replacement of point charges by Gaussian charges, 201 followed by an inclusion of polarizability via atomic dipoles or Drude particles and finally damping electrostatics via Gaussian distributed dipoles and polarizabilities. 202,203 Some remarkable manifestation of the importance of screening of the short-range electrostatic and polarization response was that the polarizable model with the Gaussian charges was able to accurately reproduce not only properties of water at ambient temperature but also at liquid vapor equilibrium, while previous attempt to achieve this proved unsuccessful. The developed water model satisfied the water monomer and dimer properties and simultaneously yielded very accurate predictions of dielectric, structural, vapor−liquid equilibria, and transport properties over the entire fluid range. 201 Beyond this first-generation Gaussian model, more evolved models fitting the electron density itself appeared. On the basis of Hermite Gaussians 204,205 (or Slatertype functions 206 ) that are closely related to distributed multipoles, they reproduce very accurately reference ab initio surfaces by enabling an extremely precise density-based evaluation of other contributions such as electrostatic (including penetration effects) exchange repulsion. The extension to particle Mesh Ewald of Hermite Gaussian treatment 205 make such models available in molecular dynamics. 207 Although still in development, 208 such methods are able to reproduce quasiexact quantum permanent electric fields in standard QM quality. 199 They recently allowed to understand a little bit further the dual level of nonlinear effects discussed by Chelli et al. 200 It was reported by Piquemal and co-workers, 199 that in the case of high electric fields generated by metal cations, high-level Gaussian models were not able to fully recover the ab initio polarization but still required some damping function. In other words, Gaussian electrostatics (exact electric fields) can deal with the linear regime discussed by Chelli but necessitates damping or Gaussian dipole screening to recover the second intramolecular nonlinear effect which can be traced back to exchange-polarization effects 199 due to the neglect of Pauli repulsion. Such observations led recently to simple modifications of the Thole model. 209 In the new framework, the Thole parameter for the direct (permanent) field was chosen to be different from the current Thole damping value used for the mutual induction, which leads to a significant improvement of results by separating the two nonlinear effects.

Parametrization of Polarizable Force Fields
A suitable force field model balances the needed accuracy, computational costs, transferability of parameters, and available software supporting particular potentials. Simple models that are readily available in simulation packages often dominate until the community learns about their drawbacks promoting more accurate but often computationally more expensive models to be implemented in common MD codes. For example, water models progressed from simple three site to 4, 5, and 6 sites and/or multipoles, 210 or Gaussian screened charges 211 in order to improve the description of the electrostatic potential followed by the inclusion of polarizability via induced point dipoles or Drude particles. Further improvements included damped electrostatics using Gaussian distributed dipoles and polarizabilities. 202,203 More sophisticated polarizable models such as SIBFA or EFP are now fast enough to be applied to ionic liquids. 212 Also, multipolar models originating from the biological applications, such as AMOEBA, propose new parametrizations. 213−215 An extended discussion on force field parametrization is out of the scope of this review. Hence, we will focus primarily on the parametrization of electrostatic interaction parameters, e.g., partial charges and polarizabilities, which should be most important for modeling ions, and their balance with the corresponding van der Waals parameters.
One approach to develop polarizable force fields (e.g., for ILs) is to use the existing nonpolarizable force fields as a starting point for the parametrization as depicted in Figure 11. Molecular polarizabilities α i can be deduced from experimental refractive indices n D or calculated quantum mechanically. However, the decomposition into atomic polarizabilities α iβ is not straightforward. Section 2.5.2 describes statistical approaches (shown in green in Figure 11) as well as QM procedures (shown in orange) to get atomic α iβ values. Some of these methods (relay matrix optimizations and the electrostatic grid-based approach 216 ) already include the Thole screening functions (see section 2.4.3). 163,189,190,217,218 After addition of the atomic polarizabilities, the intramolecular potentials have to be readjusted because the induced dipoles may affect internal torsions as well as the intramolecular geometries and general vibrational properties.
Importantly, induced polarization has a direct impact on intermolecular interactions. Because nonpolarizable force fields already include average dispersion between the molecules, existing Lennard-Jones parameters need to be reparametrized. The protocol for this optimization is described in section 2.5.3. Finally, condensed phase MD simulations are performed to finetune and validate the parameters of the force field assembled according to the above procedure. Depending on the agreement with experimental data and the computational setup, a number

Chemical Reviews
Review of parameters may be slightly adjusted to fine-tune the force field following the protocol summarized in section 2.5.4. Please note that the approach illustrated in Figure 11 practically relies on the assumption that an existing nonpolarizable force field is a good starting point, which is not necessarily the case, even if the simulations using a nonpolarizable force field were tuned to accurately describe some specific properties of the system of interest. There is no guarantee that in the original empirical adjustment of the nonpolarizable force field the approximation of average polarization effects has been uniformly distributed into van der Waals interactions and no artificial imbalance of repulsion−dispersion interactions has been introduced. Thus, it is generally better to start development of the repulsion− dispersion interactions for the polarizable force field from scratch by fitting quantum chemistry data and experimental measurement for density, heats of vaporization, and solvation energies to improve the force field transferability.
2.5.1. Intramolecular Potentials. In contrast to simple atomic ions like halides, ILs are charged molecules. Consequently, care must be taken to optimize the intramolecular parameters to accurately treat the geometries, vibrational properties, and relative conformational energies. The geometries are strongly influenced by the bond and angle equilibrium parameters and dihedral multiplicities and phase shifts, while the vibrations are governed by the force constants. The conformational energies are dominated by the dihedral parameters. Importantly, it must be remembered that the intramolecular parameters, especially the dihedral force constants, are coupled with the nonbonded parameters. Accordingly, when any parameter of a force field is changed, it is necessary to check all aspects of the model with respect to reproduction of the target data, with additional optimizations performed as required. Generally, such an approach requires only one to two iterations, especially when the initial guess parameters are of high quality such as those from the Canongia-Lopes et al. force field (from OPLS), 192 CGenFF, 219,220 or APPLE&P. 44 As intramolecular parameter optimization procedures have been described in detail elsewhere, 114,221 we only present them briefly here. Equilibrium, multiplicity, and phase term optimization typically targets QM geometries obtained at the MP2/6-31G(d) or higher level chemistry. When available, information from experimental data, such as microwave spectra and crystal structures, may be used as target data. When using crystal structures, ideally the geometries are obtained from MD simulations of the molecule in the crystal environment at the experimental temperature. A nice alternative is the use of crystal survey data that can be taken advantage of when large numbers of structures containing the molecular connectivity of interest are available. Force constants are optimized by targeting QM vibrational spectra, although experimental spectra may be used when available. When performing such optimizations, it is important to reproduce the contribution of the intramolecular degrees of freedom to the individual frequencies (i.e., the potential energy distribution) and apply the appropriate scale

Chemical Reviews
Review factor for the QM frequencies associated with the model chemistry used. 222,223 Targeting of the vibrations is used to optimize the bond and angle force constants as well as those of nonrotatable dihedrals and dihedrals terminated with hydrogens, excluding those terminated by hydroxyl or sulfhydryl groups. Optimization of dihedral force constants along with the determination of correct multiplicities are usually performed on the basis of QM potential energy scans. QM potential energy surfaces are typically calculated using electron correlation (e.g., MP2/6-31G(d) optimized geometries with single point energies at the MP2/aug-cc-pVTZ or MP2/cc-pVQZ model chemistry). Accurate treatment of the dihedral parameters of rotatable bonds is quite important for obtaining the correct conformational energies, although all the intramolecular terms contribute. Indeed, accurate optimization of the intramolecular parameters is essential to ensure proper treatment of the intramolecular distortions molecules undergo during MD simulations.
2.5.2. Electrostatic Interactions. In addition to the permanent charges q iβ already present in nonpolarizable MD force fields (an excellent review is given by Holm and coworkers 48 ), atomic polarizabilities α iβ are responsible for the strength of the induced dipoles emerging from the local electric fields at the position of the atom β of molecule i. In a first attempt to optimize these parameters, one may use the corresponding polarizabilities α iβ reported for neutral molecules. 189,224−227 However, it is typically preferable to derive new parameters for both neutral and ionic species, either from experimental data or QM calculations as depicted in Figure 12.
Further optimization and validation of electrostatic parameters can take advantage of a range of condensed phase experimental data. This general philosophy of exploiting experimental data to empirically improve the model has been applied successfully to optimize the van der Waals parameters in the case of additive nonpolarizable force fields. 228 On the basis of a large database of experimental densities ρ and refractive indices n D of a range of structurally diverse ILs averaged molecular volumes ⟨V i ⟩ and molecular polarizabilities ⟨α i ⟩ of the cations and anions can be determined assuming no particular correlation between the molecules (see Figure 12). 145,229 As the molecular composition of the ion pairs is known, a "designed regression"-analysis yields atomic volumes ⟨V iβ ⟩ and polarizabilities ⟨α iβ ⟩ of each chemical element involved, assuming again that there is no particular correlation between the atoms with respect to these properties. These atomic values can be used to predict molecular polarizabilities of ion pairs which were not part of the initial database. Notably, these predictions match QM calculations at the MP2/aug-cc-pVTZ level 229,230 despite the crude approximations made. Furthermore, the predicted density can be used as a starting point for the simulation box length in a constant pressure equilibration run. In addition, this approach is superior compared to classical quantitative structure−property relationship models 231 as physically mean-   145,195,196,218,a designed regression relay matrix AIM

Chemical Reviews
Review ingful descriptors are used in the designed regression approach without losing accuracy in the predicted refractive indices. As each chemical element was assigned an averaged atomic polarizability and volume, it does not matter if this atom is part of the anion or cation. Also, the position of methyl groups within the cations depicted in Figure 13 does not change the predicted results as the composition stays the same. However, the predicted refractive indices n D and mass densities ρ agree very well with all experimental values with a maximum deviation of 0.2% and 2%, respectively.
The biggest discrepancy to experimental refractive indices is found for dicyanamide based ILs. 229 This is due to the sp hybridization state of the carbon and nitrogen of the cyanogroups. During the designed regression analysis in ref 229, all carbons were treated the same way irrespective of the hybridization state. Because in ionic liquids most nitrogens are sp 2 hybridized in the imidazolium rings and carbons are sp 2 for the aromatic rings or sp 3 for the aliphatic chains, the different electronic environments for triple bonded carbons and nitrogens are not well represented. In a subsequent paper, 145 the hybridization state of the carbons was explicitly taken into account, resulting in much better predictions for dicyanamide based ILs. Originally, the hybridization state of nitrogens and oxygens was also incorporated. The results did not improve very much, and it was very difficult to assign the correct hybridization state based on the chemical formula for some compounds. For example, sp 2 and sp 3 hybridization state of nitrogen had to be detected sometimes by planar or tetrahedral configurations obtained from QM calculations. However, if one has to perform these calculations, the additional effort to determine the atomic polarizabilities, e.g., by means of AIM models (see below), is negligible.
The resulting atomic polarizabilities of ref 145 do not differ substantially from values derived from noncharged species 189,218,224,225 and are depicted in Table 1. The polarizability of bromine and iodine in Table 1 was obtained from the refractive index and density of imidazolium-based ionic liquids containing mixed polyhalides using the atomic polarizabilities of ref 145 for the nonhalide atoms.
On the basis of the ⟨α iβ ⟩ in Table 1 In a first attempt, this excess polarizability or its lack is distributed equally over the corresponding atoms of the molecule q Taking α e into account, the atomic polarizability of sp 3 carbons differs if they are part of an anion or cation. In the latter case, the polarizability is slightly less. The polarizability of sp 3 carbon also decreases comparing C 8 mim and C 2 mim due to the respective number of atoms.
As observed in Figure 14, the polarizability seems to be a linear function of the respective volume. Consequently, Uhlig et al. 232   analysis 234,235 of single ions or ion pairs. These AIM volumes differ from those V iβ computed by designed regression as only intramolecular space is attributed to the atoms. In contrast, the designed regression values contain also space between the ions in the liquid phase. As a consequence, the linear correlation between the atomic polarizability and atomic volume in eq 2.5.3 may be violated. 196 So far, no interaction between the induced dipoles is assumed. However, the local electric field E ⃗ iβ of atom iβ also contains a contribution from all other atoms jγ as shown in eq 2.3.5 connected via the distant-dependent dipole−dipole tensor T. The inverse of the relay matrix R in eq 2.3.6 consists of N 2 3 × 3 submatrices, which have to be summed up to get the molecular polarizability tensor α i  Table 1. Again, the net charge of the molecule has an impact on the atomic polarizabilities α iβ . Hydrogens, nitrogens, oxygens, chlorines, and bromines follow the expected trend and have higher average polarizabilities in the anions. However, boron and phosphorus show the opposite behavior.
The statistical approaches described so far rely on experimental data (see Figure 12). In principle, QM calculations of the molecular polarizability are readily performed 230,236 and use Stark's relation which links the energy U to an applied external electric field E ⃗ . The subscripts a and b denote the x-, y-, and z-directions. Either the second derivative of U or the first derivative of the dipole moment μ with respect to the electric field yields the components of the polarizability tensor α⃡ The average molecular polarizability α i is defined as the third of the trace of the polarizability tensor. The total molecular dipole moment μ can be evaluated relative to an arbitrary reference point. However, for charged molecules i, it turned out that the center of mass r⃗ i is an appropriate choice as a reference site. 237,238 Hence, the molecular dipole moment μ⃗ i reads where Z iβ is the nuclear charge of the atom β at the nuclear position r⃗ iβ and ρ(r) is the electron density. Following an AIM approach, 195,234,239 nonoverlapping atomic integration basins Ω iβ can be defined for each atom iβ as shown in Figure 15. The decomposition of the molecular dipole moment into atomic contributions i i μ μ ⃗ = ∑ ⃗ β β is then realized via atomic charges q iβ and atomic dipoles: The atomic charge q iβ includes the nuclear charge and some localized surrounding electron density. Consequently, μ⃗ iβ C describes the charge contribution to the dipole moment. Within the atomic basin Ω iβ , the electron density can polarize and contribute to an atomic polarizable dipole moment μ⃗ iβ P . Applying an external electric field E ⃗ perturbs the charge distribution of the molecule. The electron density can be transferred from one atom to another, giving rise to Δμ⃗ iβ C . The electron density around an atomic site also responds to E ⃗ changing the polarization within the atomic basins ΔΩ iβ and giving rise to Δμ⃗ iβ P . 195,234,239 To make the calculation of the polarizability (see eq 2.5.6) independent of the origin, the charge distribution contributing to the dipole moment μ⃗ iβ C can be converted to a sum of surrounding bond charges b βγ (2.5.9) For example, the partial charge of the C2 atom in Figure 15 equals the sum of the bond charges q C2 = b C2,H2 + b C2,N1 + b C2,N3 . Because these bond charges describe directed contributions, the  Figure 15. To determine the b βγ values, a nonunique set is chosen where γ > β fulfill eq 2.5.9 for each atom in the molecule. 195 As hydrogens are connected to other atoms via a single bond, the respective bond charges are set to the value of the partial charge of the hydrogen, b H,γ = q H . Finding the best solution for the bond charges yields atomic contributions to the charge term in eq 2.5. (2.5.11) with the position of the bond charge midway between the atoms,  Table 2 using M06-2X/Sadlej pVTZ for the QM calculation and an electric field strength of 0.0008 au = 0.041 V/Å. The agreement with the corresponding designed regression values is quite reasonable (taking into account that different functionals/basis sets alter the results slightly) for the sulfur and oxygen in the triflate and the carbons in the dicyanamide, with an exception being the acidic H2 of the imidazolium ring. This indicates that the averaged designed regression values do not completely resemble the current situation of the ionic liquid. In the case of imidazolium H2, the acidity and hydrogen bond capability are known to be enhanced compared to H4 or H5. The electron density of triflate in the bond charge model and consequently its ease to be distorted shifts from the sulfur to the oxygens. In the case of the dicyanamide, the shift is from the carbons to the nitrogens. Atomic polarizabilities α iβ (gained by the AIM approach described above) and partial charges q iβ derived from RESP calculations of several thousand molecules were also used to train a linear increment scheme as well as a neural net. 242 On the basis of these data, fast prediction of q iβ and α iβ with average errors of <0.02e/0.03 Å 3 and <0.07e/0.07 Å 3 are possible using the neural net and linear increment system, respectively. This method to determine the electrostatic parameters may be handy for large molecules where QM calculations might be too expensive computationally.
However, for small molecules, atomic polarizabilities α iβ , corresponding Thole screening parameters and atomic charges q iβ can be optimized in the standard CHARMM procedure on the basis of a QM calculation of the electrostatic potential (ESP). To determine α iβ , a test charge was placed on various positions of a grid and a series of perturbed ESPs is generated for each charge location around the molecule and the electrostatic parameters were optimized to reproduce the perturbed ESPs. 243 This procedure was implemented by Roux and co-worker 216,244 in a general force field generation tool general automated atomic model parametrization (GAAMP) and was used during APPLE&P force field parametrization for battery electrolytes. 44,245,246 Here, this procedure was followed to determine q iβ , polarizabilities α iβ , and the corresponding Thole screening parameters (s 2 ) to reproduce the perturbed ESPs, with the restraints included to prevent nonphysical values. The resulting atomic polarizabilities are depicted in Table 1 and differ slightly from those values obtained from designed regression 145 and the AIM 196 analysis in the case of the cations. However, terminal oxygens and fluorine in trifluoromethylsulfonate have significant lower polarizabilities as well as the terminal nitrogens of the dicyanamide anion.
2.5.3. van der Waals Interactions. Nonpolarizable force fields already contain Lennard-Jones parameters characterizing the exchange repulsion and dispersion of the molecules. Simply adding induced dipoles to these force fields results in an overestimation of the attraction between the molecules. 122,[136][137][138]190 This leads to a significant increase in density for polar neutral solvents, like methanol and N-methylacetamide (NMA), and a moderate increase for various ionic liquids as visible in Figure 16. This was also observed for other ionic liquids. 52 Apolar liquids, e.g., hydrocarbons, do not show an increase in density because the partial charges on the atoms of hydrocarbons are low and consequently the local electric field is negligible. Hence, the induced dipoles in these systems are weak (see eq 2.5.5). In contrast to polar, neutral liquids, the attraction in ionic liquid is dominated by Coulomb interactions between the cations and anions. Although the strong induced dipoles increase the density in these ionic liquid systems, the effect is counteracted by a reduction of the electrostatic interactions as the partial charges are now immersed in an "inner solvent" 45,61,62,190  This reduction in Coulomb interaction is stronger for the ionic liquids compared to the neutral polar solvents and consequently the increase in density is weaker.
For the design of a polarizable force field, several strategies have evolved: (1) All Lennard-Jones parameters are reparametrized after adding the polarizable forces to the simulation.
(2) The ratio between the interaction of the induced dipoles and dispersion can be determined by DFT calculations using symmetry-adapted perturbation theory. 247 As shown in Table 3, dispersion entirely dominates the interaction between hydrocarbons, 248 whereas the contribution of the induced dipoles become almost an equal partner in case of the interaction between the ions. The two values for the NTf 2 -anion discriminate between the oxygens or the fluorine of the NTf 2 being closest to the hydrocarbon.  Here, α max is the highest atomic polarizability in the system and Δα is the difference between α max and the polarizability of the current atom α iβ . The scaling parameter λ varies between zero (disregarding ϵ iβ for the atom with the highest polarizability) and one (no scaling at all). It can be determined by comparison of computational and experimental data of the mass density and conductivity. 140 (4) The molecular dispersion coefficients C 6 can also be determined via the isotropic dynamic polarizabilities: 232,250 16) and then scaled to the atomic property by the respective volumes: 232,251,252 Another possibility to obtain C 6 coefficients is from maximally localized Wannier functions. 253 A choice of combining rules for the repulsion−dispersion interactions is also important but is often hard to assess due to the presence of electrostatic and polarizable interactions in polar or charged molecules. Therefore, simulations of noble gases and mixtures of relatively uncharged hydrocarbon molecules can be good candidates for the evaluation of combining rules. In noble gases, standard arithmetic and geometric combining rules perform poorly. 254 Combining rules used in the OPLS-AA force field fail to predict the enthalpy and volume of mixing for nalkanes and fluoroalkanes while modified Waldman−Hagler combining rules for the repulsion−dispersion parameters used in APPLE&P accurately predict these mixing properties. Therefore, the application of the Waldman−Hagler combining rules in simulations of ionic system should provide an accurate description of the fluorinated and nonfluorinated parts of solvents and anions that are of high interest in battery electrolytes. 255 The repulsion−dispersion parameters in polarizable force fields are quite transferable. For example, in the APPLE&P, 47 only one set of oxygen, hydrogen, and carbon repulsion−dispersion parameters is used with a minor exception of slightly different parameters for the terminal methyl group. Yet the densities and heat of vaporization were accurately predicted for more than 30 common ILs, a wide range of solvents (alkanes, fluoroalkanes, ethers, carbonates, sulfones, phosphates), and a wide range of electrolytes from low to high salt concentrations indicating good transferability of a selected set of repulsion−dispersion parameters and combining rules. 44,256−261 2.5.4. Fine Tuning. To improve the computational efficiency, polarizable Drude MD simulations are propagated using a dual-thermostat extended Lagrangian approach to approximate the time-consuming self-consistent field calculations of the induced dipoles. 130 The approach involves assigning all Drude particles a mass of 0.4 amu which is subtracted from the mass of the respective polarizable atom. Hydrogens are not polarizable in these Drude-MD simulations to reduce the number of Drude particles. In principle, it is possible to polarize the hydrogens as well using Drudes, although this requires a recalibration of the Drude mass within the extended Lagrangian approach. Alternatively, it could be propagated via slow selfconsistent field optimization of the Drude particles. 146 As the dynamics of the ionic liquids depend on the molecular polarizability, 45,135 neglecting all hydrogen polarizability can be disadvantageous, in particular for cations with long alkyl chains. A possible byway is the merging of the hydrogen polarizability with the polarizability of the atom to which they are attached to. 45,122,135,137,138,140,190 Alternatively, the GAAMP procedure 216,243 is able to assign an α iβ to these non-hydrogen atoms to reproduce the ESPs without the need to polarize the hydrogens as well. As a result, the molecular polarizability is preserved. In contrast to the Drude model, induced point dipoles are mathematical dipoles having no additional particle with an artificial mass. As a result, polarizable hydrogens pose no problems using Lagrangian thermostats. In practice, induced point dipole simulations of [C 2 mim][OTf] showed very similar dynamics if all atoms were made polarizable or the hydrogen polarizabilities merged to the corresponding carbons. 135 In case of the polarizable, coarse-grained MARTINI force field for monovalent ions 262 in aqueous solution 143,144 the ion beads carry two Drude particles (see section 2.2.1). The partial charges and van der Waals parameters were determined by running a manifold of coarse-grained simulations with incremental changes in these parameters to reproduce the density and the dielectric constant as a function of the ion concentration as close as possible.

Computational Efficiency and Benchmarking
One key aspect of polarizable force fields in term of efficiency compared to classical simulations is the mandatory additional computational cost associated with the evaluation of polarization energy. Over the years, considerable work has been performed to overcome this computational bottleneck without compromising the accuracy. Solving the polarization equations using point dipoles usually costs more than half the total cost of an MD step (depending on the force field and simulation setup). To reduce that cost and to keep accuracy, various strategies are possible. Nevertheless, some limitations exist due to the imperfect time reversibility and volume preservation that they may imply. Furthermore, the ability to parallelize the method efficiently also influences the choice of the optimal method and therefore the final efficiency.
2.6.1. The Speed versus Accuracy Dilemma. In practice, the evaluation of the polarization equations for point dipole models can be seen as the resolution of a large set of linear equations therefore requiring a matrix inversion. 263,264 Usually polarizable simulations deal with thousands to hundred of thousands of atoms, consequently, as the polarization matrix size depends on such a large number of polarizable sites, "exact" direct matrix inversion approaches such as LU or Cholesky factorization are unfeasible. Therefore, one has to resort to stateof-the-art mathematical iterative methods. 265 Iterative methods used in molecular dynamics 264,265 can be grouped in two main families: stationary methods, like Jacobi Iterations, or the Jacobi over-relaxation method, and Krylov subspace methods, such as the conjugate gradient (CG) or the Jacobi/direct inversion of the iterative subspace (JI/DIIS). Historically the first set of methods was used in the community. For example, the Jacobi over-relaxation approach was used in the context of the AMOEBA force field. 164 Gradually, as convergence issues were known in the mathematics community and because they cannot be recovered by adding more iterations, 147 they were abandoned and replaced by Krylov approaches. 147,263,264 Efficiency of CG and JI/DIIS are similar and only potentially differ when a large number of cores is used within massively parallel implementations. 147 By definition, iterative techniques are nonexact inversion approaches, so they have to embody two qualities: a low computational cost and a high accuracy on both energy and forces. Of course, the devil being hidden in the details, the standard way of computing the forces assumes that the dipoles are fully converged. Unfortunately, to enforce the quality of the nonanalytical forces, a very tight convergence threshold of 10 −5 −10 −8 D on the dipoles is mandatory but also associated with a lot more iterations leading to a slowdown of the simulation. 147,263,264 In practice, such setup being rarely chosen, the dipoles are not fully converged and thus the forces are not the exact opposite of the gradient of the polarization energy, generating errors that accumulate leading to energy drifts. This degrades the computational efficiency of the solvers and limits the use of molecular dynamics with point dipole polarizable force fields, however several strategies have recently been developed to solve this issue.
2.6.2. Fast and Accurate Algorithms for Point Dipole Models. As we just discussed, the more advanced Krylov solvers are the most suited to be employed to solve the polarization equations being ensured of a guaranteed mathematical convergence. 265 To reduce the computational cost associated with reached convergence, several techniques have been developed and aim at reducing the number of necessary iterations to do so. For example, in the context of the conjugate gradient approach, it is possible to use a preconditioner. 263,264 It consists in choosing a matrix P such that P −1 is close to the inverse relay matrix R −1 (see eq 2.3.6), and in applying the iterative method to the modified linear system where the matrix and the right-hand side are multiplied by P −1 . The convergence of the solver is then accelerated because of the clustering of the eigenvalues of the matrix P −1 ·R.
Efficient preconditioners exist with various associated complexities, 263,264 and some of them have been designed for the polarization problem, such as the ones proposed by Wang and Skeel. 263 All provide a reduction in the number of iterations required to reach convergence up to 10−20%, depending on the nature of the chemical system which, of course, impacts the condition number of the matrix that one needs to invert. If preconditioning is not suitable for Jacobi/Direct Inversion of the Iterative Subspace (JI/DIIS) solvers, solutions also exist for this approach. For example, Nocito and Beran recently introduced the faster Divide and Conquer Block-Jacobi/DIIS method 266,267 that solves the polarization equations by partitioning the molecular system into a set of smaller clusters treated at the JI/DIIS level, leading to a 10−20% speedup compared to the standard approach.
Another solution that can be added on top of the CG and JI/ DIIS strategies is to choose an initial "predictor" guess as close as possible to the actual solution of the linear equations. This guess can be constructed using information from one or a few of the past values of the dipoles. The simplest choice is therefore to take the "previous guess", i.e., the value of the dipoles at the previous time step, but one can go for more elaborate and efficient strategies such as Kolafa's Always Stable Predictor Corrector (ASPC) 268 or Skeel's Least Square Predictor Corrector (LSPC). 263 These two advanced predictor/corrector approaches reduce the number of iterations by a factor two in a standard production simulation context. One should note that these techniques lose their efficiency when one uses larger time steps. In the case of the Reversible Reference System Propagator Algorithm (RESPA) multiple time step 269 integrator instabilities occur when such predictors are used with time steps larger than 2 fs.
Close in spirit, the last strategy, derived from ab initio MD, has been shown to be successful, leading to convergence acceleration through the addition of an extended Lagrangian scheme to propagate a set of dipoles that are used as an initial guess to standard iterative solvers (iEL/SCF or Extended Lagrangian Self Consistent Field). 270 This approach offers comparable performances and time step capability as the ASPC predictor but requires an additional thermostat in order to prevent energy flows between the degrees of freedom.
At this point of the discussion, the presented acceleration strategies do speed up the computations, but all ultimately suffer from the drawback of the presence of nonanalytical forces which lead to accumulating errors. Therefore, to ensure stability of very long time scale simulations toward microseconds, they should all employ a tighter dipole convergence criterion leading to a higher number of iterations than usually discussed in benchmarks for short "10 −5 D-like" simulations. This uncomfortable diagnostic led the community to consider approaches offering analytical formulas for the polarization energy.
Analogous to the iEl-SCF approach, 270 one can consider the actual induced dipoles as new degrees of freedom and build an extended Lagrangian. It has been shown that such an approach could be extended by defining the way to propagate the dipoles during the dynamics without any SCF cycles. The first results using this nonempirical strategy are promising, and the resulting method named iEl/0-SCF 271 does not require any iteration, therefore offering higher energy conservation compared to standard iterative approach. On the performance side, the method does presently offer performance in line with iterative techniques but appears limited to standard time steps due to its use of extended Lagrangian.
Another family of methods, closer to the mathematical ideas governing the matrix inversion techniques, were introduced by Wang and Skeel. 263,272 Indeed, they introduced 13 years ago the initial idea of "analytical forces" for polarizable forces and proposed a method relying on Chebyshev polynomials enabling the simultaneous analytical expression of both the energy and of its derivatives. Despite this conceptual advance removing any source of energy drifts, the approach was not applicable to production simulation as the resulting final energies were too degraded compared to fully converged reference ones. However, a few years ago, Simmonett et al. revisited the concept and improved it by proposing the extrapolated perturbation theory (ExPT). 273 ExPT can be seen as a truncation of the Jacobi iterative method at a predetermined order combined with the use of a few parameters. Initially, as the empirical and difficult choice of parameters limited the full applicability of the approach to all kinds of systems, the authors proposed an evolution of the method. Now denoted as OPT3 (OPT = orders of perturbation theory), 274 the strategy is pushed to a higher order of perturbation, and although it still involves parameters, OPT3 provides a systematic way for the parametrization, extending the applicability of the method. On the computational point of view, the analytical aspect of the evaluation of derivatives also reduces the cost of the method compared to the best iterative approaches by roughly a factor 2, making it attractive.
To overcome the EXPT/OPT3 limitations, a nonempirical and noniterative strategy denoted the truncated conjugate gradient (TCG) has been proposed. 275,276 It is derived by explicitly writing down all numerical operations of a finite number of CG cycles of iteration. The level of TCG can be userchosen: TCG-n, n = 1,3. For a chosen TCG level, the number of operations is fixed once and for all, and it is then possible to derive an exact analytical expression of the gradient of the energy exactly in the same spirit as in the Skeel or ExPT/OPT3 approaches. By construction, it avoids any energy drift in microcanonical simulations. TCG remains a Krylov CG approach, therefore its error is monotonically reduced at each cycle: the higher the TCG level is, the higher its accuracy is. Other advantages exist as the CG-method being mathematically optimal at each iteration provides "on the fly" optimal coefficients that do not need to be parametrized as in ExpT/ OPT and therefore guarantees that the number of the required matrix-vector products (1 per iteration in any iterative approach) is reduced to a minimum compared to any other methods. In practice, the TCG accuracy can be improved at negligible costs: (1) by using preconditioners as previously introduced leading to the truncated preconditioned conjugate gradient (TPCG); (2) by using the residue of the final CG step, available without any additional cost, to perform an additional "peek" pseudoiteration, equivalent to one step of Jacobi Over Relaxation with a relaxation parameter which can be found adaptively. As TPGC3 is virtually exact, TPCG2 can be used coupled to a peek step as a production method for any type of systems ranging from biological systems to ionic liquids at a cost comparable to OPT3. Finally, being analytical, TCG does not rely on history (no predictor-corrector) and can be applied to larger timesteps for the same fixed computational cost. It can be coupled to efficient multi-timestep integrators to provide a very fast evaluation of short-range polarization. In practice, an integrator such as BAOAB-RESPA1 435 uses TCG-1 at short-range and provides up to a 7-fold acceleration of polarizable point dipole simulations compared to standard 1fs approaches without loss of accuracy of the dynamics.
2.6.3. Evaluation of Polarization: Other Sources of Acceleration. All discussed polarization solvers have the advantage to ensure convergence and to be compatible with massively parallel implementations (see for example ref 154), but in practice, to be able to achieve scalable performance on modern supercomputers, they need to be coupled to two categories of algorithms. The first category is a linear scaling algorithm that will enable the fast evaluation of electric fields that are required to evaluate the polarizable energy and forces. The most common choice is the Smooth Particle Mesh Ewald (SPME) that is a well suited (n log(n)) fast periodic boundary conditions approach that was extended to the use of point dipole models as well as the alternative P 3 M (particle−particle particle−mesh). If one wants to use distributed multipoles, especially at high angular momenta, specific recursions or tree code techniques should be used to ensure a fast evaluation of electric fields. 160,277,278 Finally, to tackle very large systems, one should rely on techniques such as fast multipoles 279,280 or limit communications using 3D decomposition techniques such as the midpoint approach that is at the heart of the Tinker-HP code. 154,281

COMPARATIVE CASE STUDIES
In this section, we focus on the discussion of influence of polarization (or the lack of thereof) on the prediction of properties for several prominent classes of ionic systems. Taking into account that each application emphasizes the importance of a specific subset of properties, this discussion is organized based on different applications and the corresponding key physical phenomena. Specifically, first we discuss aqueous solutions focusing on the properties related to biological and biomedical systems. Then, we review MD simulations of electrolytes for rechargeable battery applications, which primarily focus on the ability of polarizable and nonpolarizable models to adequately predict small metal cation solvation and transport. Bulk ILs are another class of ionic systems where interactions due to induced polarization impact important structural and dynamic characteristics. Finally, electrolytes near charged surfaces represent another challenging case, where polarization of both electrolyte and electrodes can influence the mechanisms of charge separation and storage.

Pure Solvents and Dilute Aqueous Solutions
An important quality of a force field is its ability to accurately model pure solvents and dilute aqueous solutions as these are representative of the conditions typically present in biological systems. For example, the interior of lipid bilayers or certain domains of proteins have characteristics of pure solvents, while the environments around biological molecules are dilute aqueous solutions. In addition, such systems have often been subjected to extensive experimental analysis from which thermodynamics parameters are available that may be used to optimize or validate force fields. 282−284 In this section, we will briefly summarize the various types of pure solvents and dilute aqueous solutions used in the polarizable force field development with the emphasis on Drude oscillators and the AMOEBA induced point dipole force field.
3.1.1. Polarizable Water Models. Pure solvents played a central role in the development of force fields designed for condensed phase simulations. The most widely studied pure solvent is water, and there have been numerous comprehensive reviews of this topic. 285−287 Drude Model. When the development of the Drude force field was initiated, the initial focus was on a water model. Given the success of the TIP4P 288 and related additive four-point water models, 289 an analogous four-point model was developed in which the charge sites were the "MW" (see Figure 17) or fourth off-center site and the hydrogens while the Lennard-Jones parameters and polarizability were placed on the oxygen. 148 Target data included high level QM data of the water dimer and a range of experimental pure solvent water properties. A systematic search of parameter space in which the polarizability was fixed to the gas experiment value, yielded a number of models that produced good agreement with a range of target data at room temperature, such as density and heat of vaporization, but yielded a dielectric constant systematically larger than the experimental value. This led to reconsideration of the approximation of the fixed polarizability, with subsequent optimization allowing the polarizability to vary as part of the fitting process, yielding a model with the experimental gas phase polarizability scaled by 0.72 that was in good agreement with the dielectric constant as well as a range of other properties. The model was named SWM4-DP, standing for "simple water model with four sites and Drude polarizability". The physical justification for the need for scaling the polarizability is still a point of debate, although a decreased ability of the electron cloud to distort in the condensed phase and the approximation that the polarization of the water model (or of a polarizable atom) is based on the electric field at the nucleus rather than being integrated over the entire volume 291 appear to be an important factor. 114 A modified version of SWM4-AD, 292 where AD stands for isotropic atomic dipole polarizability was also developed and showed an equivalent performance to SWM4-DP. However, while the SWM4-DP model did reproduce a range of experimental data, the model was developed with the Drude oscillator carrying a positive charge. As the Drude particle is meant to represent the electronic degrees of freedom, a second water model was developed, SWM4-NDP, where NDP stands for negative Drude polarizability. 149 This model was optimized in a manner similar to that used for the SWM4-DP model, yielding good agreement with a range of experimental data. However, because the impact of system size on the calculated diffusion coefficient, 293 (an issue that was not identified until the period during which SWM4 models were being developed) was not considered during the optimization process, both models actually overestimate the self-diffusion constant. This is due to the fact that the apparent self-diffusion constant was in good agreement with experiment while the size correction leads to a value that is too high. More recently, the SWM6 water model included two lone pairs along with the atomic and MW sites as depicted in Figure 17. During the development of this model, it was determined that six rather than just five sites (MW-site omitted) were required to model the full range of gas phase (e.g., quadrupole moment) as well as condensed phase properties. 294 The SWM6 model yielded improved agreement with respect to water clusters and for the diffusion constant relative to the SWM4-NDP model, although at an additional computational cost. Finally, it should be noted that the temperature-density profile of the SWM4 and SWM6 water models are not in good agreement with experiment, a problem that ongoing optimization efforts are addressing.

Chemical Reviews
As an alternative to the polarizable water models in CHARMM (SWM4-NDP, SWM6), several other Drude-like models exist in literature, for example, the COS/G3 model 295 developed by van Gunsteren and co-workers or the polarizable model of Dang and Chang 297 used by Salanne and co-workers for aqueous electrolytes. 298−300 The COS/G3 water model is the oldest and uses a tetrahedral angle of 109.47°as shown in Table 4, leading to a distorted water structure. POL3 is the standard polarizable water model for AMBER force fields. 301 The polarizable water model of Dang and Chang 297 is parametrized to reproduce structure and thermodynamic data of aqueous clusters, bulk phase, and liquid/vapor interface of water. However, a frequency-dependent dielectric analysis of several nonpolarizable and polarizable water models 290 showed that SWM4-DP does not only reproduce the static permittivity ε but also comes closest to the experimental relaxation time. The various nonpolarizable and polarizable water models and their properties are reviewed in refs 290 and 291.
AMOEBA Model. The AMOEBA water model uses distributed multipoles and is based on the induced dipole model for polarization. The AMOEBA energy function is the following: In AMOEBA, the short-range valence interactions are described by five terms, namely: bond stretching, angle bending, bond-angle cross term, out-of-plane bending, and torsional rotation. Such terms are complemented by the nonbonded van der Waals interactions, permanent multipolar electrostatics (up to quadrupoles), and explicit point dipole polarization that couples isotropic polarizability to Thole damping. Particularities of AMOEBA for the short-range interactions include the use of bond-angle cross terms, a decomposition of angle bending into  in-plane and out-of-plane components using a Wilson−Decius− Cross form as well as the use of the buffered 14−7 "Halgren" function 254 to deal with van der Waals interactions. Bond stretching, angle bending, and stretching−bending coupling, are identical to those of the MM3 force field. 302 Anharmonicity is taken into account using higher-order deviations from ideal bond lengths and angles. Accuracy of electrostatic interactions is ensured by expansion of the multipolar development beyond point charges and up to quadrupoles, allowing for a better representation of directional effects such as those, for example, found in hydrogen bond interaction networks. These distributed multipoles are extracted form ab initio computations using Stone's Distributed Multipoles Analysis (DMA). 303 The original 2003 AMOEBA model 164 44,309 and other force fields. 310 The advantage of pure solvents is the availability of accurate experimental data on properties including the density (or molecular volume), the heat of vaporization, the isothermal compressibility, and the dielectric permittivity, among others. This offers data that allows for, in particular, Lennard-Jones parameters to be systematically optimized, especially in the context of a hierarchical optimization approach. For example, following optimization of the Drude SWM4-NDP water model, optimization of the alkane parameters was undertaken, from which CH 3 −, −CH 2 −, and −CH− parameters were optimized. 244 Subsequently, alcohol parameters were optimized with the aliphatic parameters initially being transferred from the alkanes, such that optimization focused on the hydroxyl. 311 This was also undertaken for amides and sulfur containing species, 177,312 although in specific cases adjustments of the parameters on the aliphatic carbon covalently linked to heteroatoms was undertaken. A similar hierarchical approach was performed for benzene leading to heterocycles, 313 including nucleic acid bases. 314 A specific advantage for the use of pure solvents with the Drude force field was the ability to systematically scale the polarizabilities from the experimental gas phase values. 217 As the physical justification for this was not clear, as discussed above, an empirical approach was used in which the scaling factor was empirically optimized targeting the pure solvent dielectric constants. This approach leads to scaling factors ranging from 0.6 with some sulfur containing species, 0.7 with amides and alcohols, 0.85 with amides and heterocycles, and 1.0 with aliphatics and halogens. 114 Moving forward, a scaling factor of 0.85 is typically applied in cases where access to experimental dielectric permittivities is not available. During fitting, the APPLE&P force field for alkanes, fluoroalkanes, ethers, carbonates, nitriles, ionic liquids, and battery electrolytes, the polarizabilities of solvents and especially anions were also reduced from the gas-phase values or were fit to QM calculations with a smaller basis set that effectively scales the magnitude of polarizability. 309,315 Visscher et al. recommended to reparametrize not only the Lennard-Jones interactions when considering polarizable forces for alcohols but also the partial charges of the atoms. 310 These were scaled to improve the description of ΔH vap and ΔG hyd . Furthermore, the polarizabilities were determined in gas phase and liquid phase using a QM/MM approach. It turned out that the polarizabilities in the liquid phase of the investigated alcohols are significantly lower than the corresponding values in the gas phase which legitimates the scaling factor of the other approaches mentioned above in some way.
Of course, in the same spirit, various polarizable solvents have been derived for AMOEBA and are available for the general public. 213 3.1.3. Polarizable Ions in Aqueous Solutions. In contrast to simple atomic ions, water interactions with ionic liquids are more complex due to the anisotropic and bulky nature of the molecular ions. 316 Furthermore, multiple possibilities for (bifurcated) hydrogen bonding exist and are in competition with the interionic interactions. Hence, the need for polarizable ion models in aqueous solution has been argued several times. 317−319 Consequently, the preliminary set of Drude and AMOEBA force fields was applied in dilute aqueous solution simulations, typically in the context of estimations of free energies of aqueous solvation (or hydration free energies ΔG hyd ). These efforts include the various neutral species, monoatomic and molecular ions. For Drude, concerning neutral species, the initial hope with the polarizable model was that the pure solvent Lennard-Jones parameters in conjunction with the more sophisticated electrostatic model would yield good agreement with experiment when used to calculate ΔG hyd values. However, this was found not to be the case, leading to the use of atom-pair specific Lennard-Jones parameters for selected solute atoms with water (e.g. use of the NBFIX term in CHARMM that applies specific Lennard-Jones parameters to specific atom pairs rather than assigning them based on combination rules). 320 In retrospect, the need for atom-pair specific Lennard-Jones parameters intuitively makes sense as the variation of the electronic distribution of atoms and molecules associated with the explicit treatment of electron polarizability will be accompanied by variations in the van der Waals features of the system, such that different Lennard-Jones parameters are required for specific environments. Thus, in the absence of a model that allows the van der Waals term to vary as a function of environment the approximation of parameters for specific interacting pairs has been applied. This successfully led to good agreement for the ΔG hyd values for a range of neutral species with experimental data. 320 Moreover, a particular advantage of the Drude polarizable model over other methods to treat polarizability is the use of an explicit particle for the electronic degrees of freedom. This allows for Lennard-Jones parameters to be applied to the electronic degrees of freedom thereby modeling steric effects associated with perturbation of the electron cloud as well as electronic effects in the context of an empirical force field. 82 In the case of Drude halogens, this approach accurately treats both halogen bonds and halogen− hydrogen bond donor interactions as well as reproduces both pure solvent and ΔG hyd experimental data. 120 It has also been applied to more accurately model Mg 2+ , allowing for accurate reproduction of the ΔG hyd and of the energetics of specific water−ion interactions. 121 A central aspect of the Drude ion parameter optimization to yield a set of models that is internally consistent was the adoption of a global free energy scale based on experimental data for neutral salts. One cannot simply rely on absolute experimental values because there is no unambiguous reference scale for charged species; all experiments have to utilize some reference scale. The significance of a global solvation scale is important in force field development because the absolute hydration free energies of the different ions must be internally consistent, for example, to accurately account for the relative binding and ion-pairing affinities. As described in detail elsewhere, handling this issue properly requires that the target ΔG hyd values extracted from different experimental studies be adjusted such that the values are effectively offset to the same counterion. 149,321 Doing this assures that all the ΔG hyd values are constrained to lie on the same consistent scale, such that their relative ΔG hyd values are representative of the experimental regime. Using this approach, in combination with QM calculations, yielded a set of parameters for monoatomic ions including both mono-and divalent cationic species and monoanions. 322 More recently, these efforts have been extended to molecular ions representative of charged moieties in biological macromolecules. 321 An interesting observation from the latter study are differences in the three-dimensional probability of water around the molecule in the Drude versus the CHARMM36 additive model, indicating differences in the nature of the atomic details of the interaction of the ions with the aqueous environment.
A systematic optimization of atomic ions and subsequently molecular ions has also been reported in the context of the AMOEBA induced dipole-based polarizable force field. First, AMOEBA was tested for monovalent ions such as K + and Na + . 323 Absolute solvation free energies were accurately described for such cations as well as for the chloride ions in liquid water and formamide. Such results clearly demonstrated the ability of AMOEBA to capture the thermodynamics and free energies of solvated ions. Extension to divalent cations came later as AMOEBA was extended to Ca 2+ and Mg 2+ . 324 Such an addition required the introduction of a cation specific parametrization of the Thole polarization damping model which was required to be different in cation−water over water− water interactions and adjusted on ab initio polarization energies computed using energy decomposition analysis. This opened the door to hydration free energy for the cations which were found, again, in good agreement with experiment. 183,325 The same ab initio bottom-up strategy coupled to higher-level quantum chemistry was used to enable AMOEBA simulations of tetravalent actinides such as Th 4+ in water. 326 The first polarizable force field estimate of Th(IV) solvation free energy was then predicted. On the basis of these encouraging results, more difficult heavy metal cations were modeled more recently and hydration free energies, structures, and dynamics of openand closed-shell trivalent lanthanide and actinide metal cations were computed. AMOEBA simulations of six cations solvated in bulk water predicted first-and second shell hydration numbers, water residence times, and free energies of hydration are fully consistent with experiment (as illustrated in Figure 18) offering a predictive modeling of f-elements compounds. 327 Such parameters were later used in polarizable simulations of water and an aqueous mixture of 1-ethyl-3-methylimidazolium ethylsulfate [C 2 mim][EtSO 4 ] using the AMOEBA force field and revealed different mechanisms of water exchange processes around the lanthanides Gd 3+ , Dy 3+ , and Ho 3+ . 328 In pure water, the exchange of water molecules in the first solvation shell of the lanthanides can be explained by an associative process. In contrast, in the aqueous mixture of the ionic liquid, the  mechanism is of dissociative nature. Here, the interaction of the lanthanide with the polarizable anion plays an important role.

Experimental Data for Optimizing Pair Specific
Interactions. The studies discussed in the previous section yielded Drude parameters that are in quite good agreement with the experimental ΔG hyd data. However, accurate simulations of heterogeneous systems require that the interactions of the different solutes in the system to be balanced with those with water as pointed out within the AMOEBA results. An efficient approach to address this is using osmotic pressure calculations. 329,330 Such calculations are computational tractable and yield high precision data that can be directly compared to experimental data and may be used to optimize atom-pair specific (NBFIX) parameters between the individual solutes, including ions. In the context of polarizable force fields scaling of the atomic dipole−dipole interactions may also be performed using the through-space Thole scaling term (NBTHOLE). 163 This approach has been used to improve a number of ion−ion parameters in the context of both the additive CHARMM36 and Drude polarizable force fields 121,331 as well as for other force fields. 330 An alternative utilization for dilute aqueous solution and pure solvent data is based on the availability of X-ray or neutron scattering data. 332,333 These experiments yield atomic resolution data on the distribution between molecules and, when isotopic replacement can be performed, between specific atoms in those molecules. The utilization of scattering data for the optimization of water models has a long history and the approach can be used to gain insights into the structure of salt solutions. 334 In the context of the Drude force field neutron scattering data on methanol (MeOH) and tetrahydrofuran (THF), along with water, was used to validate the associated force field parameters. 335 Comparison of simulated and experimental partial structure factors yielded improved agreement with the Drude model for both MeOH and THF over the additive CHARMM36 force field.
The Kirkwood−Buff theory represents another approach to obtain atomistic details from condensed phase simulations that may be directly compared to experimental data. 336 The Kirkwood−Buff theory has been used for a number of years as a guide for force field development 337 and can be used to understand the interactions of solutes with macromolecules. 338 In the context of the Drude force field, the Kirkwood−Buff theory was used to validate the amide and alcohol parameters. Results with the Drude model for amides in aqueous solution led to improved interactions over the additive force field for selected molecular interactions, although deficiencies in both models are evident. 339 Studies on methanol solutions showed some properties being better modeled by the Drude force field (e.g., solution densities and dielectric constants, among others) while better performance for the additive force field was obtained with activity derivative, the excess molar Gibbs energy, and the excess molar enthalpy of mixing. 340 Thus, the Kirkwood−Buff analysis clearly indicates the need for improvement in the Drude force field for the balance of the solute−solute, solute−solvent, and solvent−solvent interactions.
Dielectric spectra of aqueous ionic liquid mixtures can already be reproduced with reasonable agreement using classical nonpolarizable force fields of Canongia-Lopes et al. 18 and TIP3P as the most important feature is the breaking of the cationanion interaction by the interstitial water molecules, accelerating the overall dynamics of the system. However, if one is also interested in conductivities of the various mixtures, nonpolarizable simulations fail for aqueous ionic liquid systems and polarizable MD simulations 140 are mandatory to get closer to the experimental values.
3.1.5. Ionic Solutions at Interfaces. Aqueous solutions of ILs at the solution/water interface were studied by Jungwirth and co-workers using a nonpolarizable and polarizable force field. 341 The water was modeled by TIP4P/2005 and POL3. 296 At low IL concentration, the nonpolarizable and polarizable force fields of the electrolyte yielded reasonable agreement with the experimental surface tension. However, at higher ion concentrations, only the polarizable force field was able to reproduce the experimental increase in surface tension. Voth and co-workers also detected an enhanced concentration of C 2 mim at the vacuum interface, going along with a reduced charge density at the surface to the vacuum. 342 Due to the reduced repulsion, the average distance between the ions is less compared to the nonpolarizable force field, which also favors anions at the vacuum interface. The surface tension of the polarizable force field was lower compared to the nonpolarizable and hence closer to experiment. Wick 343 The potential of mean force for both gases showed a preferred position at the air/liquid interface close to the cations pointing their tails toward the air phase. However, transitioning into the bulk liquid phase, the oscillatory behavior of the potential of mean force indicated that CO 2 and SO 2 interact stronger with the anions.
Interesting application studies on dilute aqueous simulations with the Drude ion parameters in contact with DNA have been reported: The osmotic pressure calculations and QM interactions of the ions with model compounds representative of the phosphate backbone and nucleic acid bases improve the balance of the ion−macromolecule vs. ion−water interactions via the use of atom-pair specific parameters, 344 analogous to the use of the osmotic pressure calculations discussed above. This optimization in conjunction with the optimized Drude DNA and ion parameters was shown to yield a model that gave improved agreement with counterion condensation theory 345 with respect to the neutralization of the DNA charge by the ions. 344 Significantly improved agreement over the additive CHARMM36 force field was also obtained on the competition of counterions (e.g., Na + competition with Li + , K + , or Rb + ) for the ion environment around DNA, yielding improved agreement with buffer exchange-atomic emission spectroscopy experiments. 346 Results from that study also reinforced the need for large simulation systems (e.g. solvation for 25 Å beyond the DNA) for proper sampling of the ion environment of the DNA. Another interesting observation was that the presence of different counterions in DNA simulations leads to changes in the calculated scattering spectra of the DNA with the Drude force field. 347 This effect, which is associated with cooperative basewater−ion hydrogen bond interactions in the grooves of the DNA, is not present with the additive CHARMM36 force field, showing that the DNA conformational properties are sensitive to ion type in the polarizable model, a property that is not present in additive force fields. Similarly to these studies, ion interactions with nucleic acids were very recently studied in a paper published during the completion of this review and concerning the extension of AMOEBA to nucleic acids that now enable full polarizable AMOEBA simulations of complex DNA and RNA systems with various ions. 348

Modeling of Battery Electrolytes
Rapid development of batteries for portable electronic and automotive applications highlighted the need for fundamental understanding of transport mechanisms in battery electrolytes at molecular scale. Seven major classes of electrolyte systems are attracting the attention of the modeling community: (1) Liquid aprotic electrolytes that are widely used in current lithium ion batteries. These simulations are mostly focused on the traditional baseline chemistries comprised of aprotic solvents such as ethylene carbonate (EC) and dimethyl carbonate (DMC). EC has the ability to efficiently dissociate Li-salts as well as its reduction products passivate and stabilize graphite electrode surfaces, while DMC or other linear carbonate or ester solvates decrease electrolyte viscosity and improve ion transport at low temperatures. (2) Liquid aprotic electrolytes for the "beyond Li" chemistry with an emphasis on Mg, Na, and Zn. (3) Aqueous electrolytes due to their intrinsic nonflammability, ease of processing, fast bulk, and interfacial ion transport. (4) Solid polymer electrolytes that are investigated with the aim to eliminate volatile organic components, to increase mechanical stability and flexibility, and to decrease the dendrite growth, therefore allowing usage of Li metal anodes. Poly(ethylene oxide) and other polyethers doped with lithium salts such as Li [TFSI] are considered as baseline systems for this electrolyte class. (5) IL electrolytes doped with the Li, Na, Mg, or Zn salts are investigated due to their low volatility that leads to safety advantages such as delayed thermal runaway and a large variation of available cation and anion combinations that provide an opportunity to tailor electrolyte properties. (6) Solid state conductors that form interphases at the electrodes due to electrolyte reduction and oxidation. (7) Hybrid electrolytes combining multiple classes of electrolytes.

Solvent Polarization and Cation Coordination.
Accurate prediction of the structure and transport in battery electrolytes requires accurate representation of the binding energy of small cations such as Li + , Mg 2+ , or Zn 2+ with solvents and anions. Because of the small size of the Li + cation, about 30% of its binding energy with ether, carbonate, or water comes from the induced polarization interactions, indicating a strong need to include these interactions either in a mean-field sense through the two-body terms or explicitly through the atom dipole polarization or Drude model. 245,246,349 The analysis of the distribution of dipole moments extracted from DFT calculations clearly revealed that molecular dipole moments for typical battery solvents such as EC or propylene carbonate (PC) near Li + are about 50% larger than the gas-phase values and 20% larger than the average dipoles in neat liquid as shown in Figure  19. An increased dipole moment of solvent molecules near Li + observed in DFT calculations is in stark contrast with the scaledcharge mean-field approach that yields smaller dipole moments due to charge reduction. Yet, counterintuitively the charge scaling was empirically found to improve the ion transport often at the expense of predicting a larger coordination number around metal cations.
Another interesting observation made by Pollard et al. 67 is that unlike the increased dipole moment of EC and PC molecules coordinating Li + ion, no increase in the dipole moment of water molecules coordinating halide anions was observed. 349,350 This is likely due to the larger size of anions compared to Li + cation, weaker anion−solvent binding energy and the difference in the hydrogen bonding network formed near anions. The larger size of anions that are of interest for battery electrolytes results in even weaker anion−solvent interactions compared to halide−water. Thus, the solvent dipole is significantly enhanced near small cations such as Li + and unchanged near anions. This consideration is often overlooked during parametrization of the scaled charge (or other mean-field polarizability) models. For the large size cations, e.g., the Drude Mg 2+ model and the SWM4-NDP water, 121 the inner shell waters around Mg 2+ have a dipole of 2.94 ± 0.11 D on average versus 2.48 ± 0.18 D for the water model in bulk solution. In the development of the Drude Mg 2+ ion parameters, the Lennard-Jones potential was applied to the SWM4-NDP water Drude particle−Mg 2+ atom pair interaction that effectively buffered the polarization response of water when coordinated with the ion. This was shown to be essential for the accurate treatment of both the hydration free energy of Mg 2+ and the kinetics of water−ion binding.
Scaled charge force field molecular dynamics (FFMD) simulations were reported for the dilute and low concentration electrolytes based upon LiPF 6 in EC and PC solvents by Chaudhari et al. 351 and are compared with the results from ab initio molecular dynamics (AIMD) simulations. A comparison of the full charge FFMD with AIMD revealed that the Li + first coordination shell is significantly more structured in the FFMD as illustrated in Figure 20  solvent molecules in the Li + coordination shell in carbonatebased solvents, where four to five molecules participate in cation coordination, significantly contributes to the Li + transport mechanism. 354 As the partial charges of atoms in the solvent molecules decrease, i.e., reducing the molecular dipole moment, the first solvation shell of Li + becomes less structured and shifts to larger separations. This leads to an improved agreement of the magnitude of the g Li-O (r) first peak from FFMD with the AIMD results. While with the reduction of partial charges the magnitude of the g Li-O (r) peak decreased, the peak becomes broader and the Li + coordination by EC and PC is effectively increasing as the Li + first solvation shell becomes more diffuse. Thus, the resulting coordination numbers from FFMD with reduced solvent dipoles are larger than AIMD results. 351 Authors recommended a scaling factor of 0.8. Using Bader charge analysis that indicated that the Li + cation charge is reduced by 0.1e, one can partially justify the scaling of charges 351 in these electrolytes, but that is inconsistent with the increased dipole from the condensed phase DFT shown in Figure 19. A similar scaling factor of 0.75 was obtained for water using the ECC 64−66 where scaling the charges of all ions by the inverse of the square root of the electronic part of solvent dielectric constant 1/ ϵ ∞ was applied. It is clear from Figure 20 that it is necessary to decrease the ion−solvent repulsion parameters to avoid the unrealistically high solvent−cation coordination numbers. This correction was made for aqueous electrolytes 64−66 but not in simulations of LiPF 6 in EC and PC by Chaudhari et al. 351 The AIMD simulations of more concentrated EC:DMC electrolytes with LiPF 6 at solvent:Li ratios of 10:1 reported a slightly higher magnitude (35) for the g Li-O (r) of EC peak. 355 MD simulations employing many-body polarizable force fields such as APPLE&P and using full charges for ions (and no rescaling for solvents) reported the first peak of the g Li−O (r) around 30 for EC:LiPF 6 and EC: Li[TFSI] and EC:LiPF 6 electrolytes, which is in good agreement with AIMD results. 260,309 Importantly, the electrolyte conductivity and the extent of ion aggregation were also accurately predicted using the APPLE&P force field for a variety of electrolytes comprised of carbonates, 309,354 glymes, 246,354,356 sulfones (SL), 257 acetonitrile (AN), 357−360 and water 258,259 doped with LiPF 6 , LiFSI, LiTFSI, NaTFSI, LiDFOB, and NaOTf salts. 259 3.2.2. Ion Transport. In dilute and moderately concentrated electrolytes such as traditional 1 M, the Li + often moves primarily with its solvation shell with minimal contribution from exchange of solvent molecules in the shell. Thus, an overestimation of the Li-solvent residence time leads only to moderate errors. As the salt concentration increases and the exchange of solvents and anions in Li + coordination becomes increasingly important for Li + transport, 360,361 the full charge nonpolarizable force fields become increasingly inadequate. For example, Takeuchi et al. 362 reported good predictions of structural properties for PC doped with one molar LiBF 4 or LiPF 6 using nonpolarizable force field while the predicted conductivity was about one order of magnitude lower than in experiment. Another recent MD simulation study from the LBNL group by Rajput et al. 43 clearly demonstrated that ion and solvent diffusion coefficients predicted from MD simulations increasingly deviated from experiments as the LiTFSI concentration increases as shown in Figure 21a. In fact, the predicted Li + and TFSI diffusion coefficients were more than  Figure 21b. 246 MD simulations employing a modified polarizable APPLE&P force field also predicted self-diffusion coefficients for the Li + and Na + cations, TFSI anion, and the same DME solvent in excellent agreement with experiments as reported by Liyana-Arachchi et al. 356 and shown in Figure 22. While the polarizable force fields usually do not require charge rescaling, a slight charge scaling by a factor of 0.94 for the Li + and anion oxygen charges was reported to be important in order to obtain excellent agreement with experimental conductivities at low temperatures over a wide concentration range in H 2 O:LiTFSI and sulfolane(SL):LiFSI electrolytes as shown in Figure 23a,b. This scaling also improved the ability of polarizable force fields to predict ion and water self-diffusion coefficients over a wide concentration range for H 2 O:LiTFSI as shown in Figure 23c. Without scaling, the ion dynamics was up to two times slower compared to experiment at the highest concentration but quite similar at moderate concentrations. 363 On the other hand, the conductivity of concentrated acetonitrile(AN):LiTFSI (AN:Li = 3) predicted from simulations using the same force field agreed well with experiments without charge scaling for a wide temperature range as shown in Figure 23b.
The scaling of the partial charges is rather motivated by mimicking the polarizability than an actual charge transfer. 45,[61][62][63]145 Consequently, one should use a scaling factor f for the partial charges during trajectory production but full charges for analysis. Although the dynamics increases very strongly with decreasing scaling factor and hence an agreement between experimental and computational conductivity is also possible when using scaled charges for its evaluation, the actual values of these scaling factors f̃are quite low (f̃≪ f), much lower than values (≅f) obtained by quantum-mechanical calculations of an ion pair. In addition, these very low scaling factors f̃also have to be applied to the computation of the dipole moment. As a result, the corresponding dipoles vanish and their correlation functions cannot contribute to the dielectric spectrum. On the

Chemical Reviews
Review other hand, if full charges are used for the analysis, dielectric spectra of charge-scaled using the factor f and polarizable MD simulations almost coincide. 45 A more rigorous strategy for designing liquid electrolytes and solid electrolyte force fields was implemented by Jorn et al., 364 who used force matching to fit a nonpolarizable force field to PBE-based AIMD results on a smaller cell. A very good agreement between the two-body and APPLE&P polarizable force field was reported for both Li + solvation and diffusion coefficients in EC-LiPF 6 electrolytes, therefore indicating a great promise of such numerical force matched force fields for bulk electrolytes. Despite the great ability of this effective two-body force field to predict bulk properties of EC-LiPF 6 electrolyte and dilithium ethylene dicarbonate (Li 2 EDC) model solid electrolyte interphase in a good agreement with predictions from polarizable force field-based simulations, 365 the interfacial kinetics of Li + transfer between electrolyte and SEI phases is quite different in the nonpolarizable and polarizable force fields. This indicates a potential limitation of the effective two-body polarization treatment when it is extended to simulations of interfacial properties while being parametrized using bulk properties. 366 MD simulations also investigated ILs doped with Li + , Na + , and Mg 2+ ion-based salts as alternatives for battery electrolytes. Polarizable force fields generally yielded accurate predictions of the transport and structural properties with the representative temperature dependent conductivity and ion self-diffusion predictions shown in Figure 24. Importantly, self-diffusion coefficients of organic and metal salts were accurately predicted. Special attention had to be paid to the inclusion or exclusion of the intramolecular polarization. For example, inclusion of full induced dipole-induced dipole interaction between oxygen atoms on the same TFSI anion (see Figure 24c for structure) penalizes the bidendate orientation of TSFI near the cation in which Li + complexes with two oxygens from the same TFSI. As a result, simulations predict a lower fraction of the Li + −bidentate complexes. 36,367 Decreasing oxygen polarization or exclusion of the induced dipole-induced dipole interactions between oxygens on the same molecule stabilizes the bidentate complex and improves agreement with experimental data.
The charge scaling approach was extensively discussed on the basis of MD simulations of ILs doped with Li + and other salts. The Maginn group 371 noted that "some authors have proposed the use of scaled charge models, in which the charges are uniformly scaled by a factor of 0.8 to represent the charge transfer and the polarization. The use of this kind of model provides better results to dynamical properties such as self-diff usion coef f icients, viscosity, and conductivity without the high cost of the polarizable force f ields. On the other hand, some works have shown that the use of f ull charge models (total ion charge 1.0) provides better results for structure and density, sometimes comparable to polarizable force f ields." A compromise between structural and dynamic properties is typically found on a case by case basis 371 with a typical range of scaling factors from 0.6 to 0.8. In another example, in the MD simulation study of Li + containing a dual-cation ionomer, the charge scaling decreased the ion−ion interaction distance and increased the size of ion aggregates together with speeding up dynamics. 371 Another attempt to incorporate polarization involves an additional short-range two-body function that sharply decays to zero beyond the first Li + solvation shell. 69,70 It was supposed to account for the increased dipole in the first coordination shell but not beyond it. This effective two-body approach also overestimated the size of the Li + coordination shell and predicted a slower Li + −ligand exchange and transport, 70 while MD simulations using polarizable force fields accurately predicted transport and ion aggregation in poly(ethylene oxide)-based electrolytes and molten salts. 261,372

Ionic Liquids
Numerous MD simulation studies of room temperature ionic liquids comprised of various cation/anion combinations have been investigated using polarizable and nonpolarizable force fields and reported thermodynamic, structural, and dynamic

Chemical Reviews
Review properties. These data provide a good opportunity for assessing the importance of inclusion of polarization in systems with chemically diverse sets of cations and anions. Despite an extensive amount of literature on modeling ILs, the direct comparison of simulation predictions using polarizable and nonpolarizable has to be made with caution. The existing nonpolarizable and polarizable force fields usually have different legacy and history of parametrization and empirical adjustments. 29,44,48,192,233,262,373 This leads to variations not only in parameters related to electrostatic interactions 48 but also the parameters for valence and van der Waals interactions. 24,29 For example, when the parameters for the latter are fitted to match ab initio/DFT binding energies or are empirically adjusted to match certain experimental data (e.g., density ρ, heat of vaporization ΔH vap , etc.) they will most likely be quite different depending on whether a polarizable or nonpolarizable model was assumed for the force field. The dispersion term (∼1/r 6 ) describing the van der Waals interactions in a nonpolarizable model might already include some effective (mean field) approximation of the induced polarization effects. 60 Therefore, it is hard to find a consistent set of polarizable and nonpolarizable models for comparison. Even when simulations do report the results for the same system using both polarizable and nonpolarizable force fields, the latter is usually reduced to simulations with a polarizable force field but the polarizability being "turned off" without any adjustment of the nonpolarizable model to effectively match the polarizable one. In such cases, while one can get a good measure of contribution of induced polarization interactions in defining a particular property of interest, it is not necessarily a fair comparison of reliability and predictability for the nonpolarizable model. Nevertheless, a couple of studies reported data where polarizable and nonpolarizable models had the same origin and extra efforts were made to make the two types of models as consistent as possible with each other. For example, Bedrov et al. 52 used the force matching approach to derive a nonpolarizable version of the force field that provides the best description of atomic forces acting on atoms during simulations of ILs using a corresponding polarizable model. In this approach, first, the simulations of IL were conducted using a polarizable model. Then a nonpolarizable model that preserves parameters for van der Waals and valence interactions and partial atomic charges, but has additional two-body interaction terms that are supposed to effectively capture all induced polarization interactions, was fitted to match as close as possible the atomic forces from polarizable simulations. Taking into account that additional terms which effectively approximate polarization interactions did not have any constraints on their functional form (i.e., they were fitted as numerical functions for each possible type of atom−atom interaction) and the fitting of forces was done for the environments corresponding to the condensed state of interest, this nonpolarizable model can be considered as consistent as possible with the polarizable model it was derived from. The comparison of IL properties predicted from this

Chemical Reviews
Review consistent set of polarizable and nonpolarizable force fields facilitates the head-to-head comparison of two types of models as discussed below. However, polarizable parameters are usually transferable, making them superior to fixed charge models as the latter assume a particular environment during the parametrization process. The resulting parameters will not always function well if those assumptions are violated. 323 Also, we look at more generic trends in correlations and properties of ILs predicted from simulations employing various polarizable and nonpolarizable force fields.
3.3.1. Thermodynamics. Chaban and Prezhdo 374 analyzed the impact of temperature on the ionic charges of LiCl, NaCl, and KCl clusters of 10 ionic pairs. Depending on the charge partitioning scheme, the partial charges of the atomic ions varied from ±0.34e to ±0.55e (Hirshfeld method) and from ±0.75e to ±0.87e (ESP charges). Furthermore, the temperature effect on the atomic charges was only visible for those derived by the Hirshfeld method. However, impact of the cations on the respective charges of chloride was more pronounced. The authors pointed out that the Hirshfeld charges measure the electron density localization within a certain radius and is therefore adequate to provide information on nonpolar, polar, ionic, and other types of chemical bonding between the atoms. In contrast, ESP charges reflect the electrostatic potential at the surface of the molecules and are consequently used for the intermolecular interactions. Chaban and Prezhdo also argued for the nonadditivity of the electronic interactions, leading to electronic polarization and charge transfer. The same authors also showed from ab initio MD simulations that in pyridiniumbased ionic liquids the influence of the temperature on the electron delocalization was minor. However, the molecular dipole moments increased as a function of temperature due to growing thermal fluctuations. 375 While the rescaling of ionic charges in nonpolarizable classical MD simulations improved the prediction of cohesive energy density, heat capacity, 33,60,376 and dynamic properties, 45,60 it can hardly be considered as a reliable approximation of induced polarization effects: (1) In MD simulations, usually no explicit hydrogen bonding potential is used. Instead, hydrogen bonding is realized via the Coulomb interaction between the positively charged hydrogen and the negatively charged, electronegative atom. In case of charge scaling, this interaction is drastically weakened, resulting in wiping out all hydrogen bonds as depicted in Figure 25.
(2) For some ionic systems, the reduction of charges still did not provide a satisfactory description as visible for the molecular dipole moments in Figure 26. Although the distribution of dipoles is broadened a little bit due to the charge scaling (compare black dashed and dotted lines in Figure 26), polarizable simulations lead to much broader distributions similar to the QM results. 45 However, shifts to higher dipole moments for the cations and slightly lower dipole moments for the anions are already visible from the charge scaling approach.   The impact of the induced dipoles on the molecular dipole is even more pronounced in the gas phase. 52 In the case of N-methyl-N-propyl-pyrrolidinium bis-(fluorosulfonyl)imide [Pyr 13 ][FSI], the cationic dipole moment in the nonpolarizable force field does not change very much between liquid and gas phase. The changes for the anion were a little bit more pronounced. However, in the case of the polarizable force field, the molecular dipole moments increased for both cations and anions. In addition, the discrepancy between liquid and gas phase values were much more pronounced, e.g., the dipole moment of FSI increased from 1.34 D in the liquid phase to 3.12 D in the gas phase.
(3) These drastic changes for the gas phase compared to the liquid phase naturally have an impact on the heat of vaporization ΔH vap as the difference of the averaged interaction energy in liquid and gas phases is a major component. However, this effect is usually hidden by the fact that ΔH vap is part of the force field parametrization process.  52 turned off the polarizable interactions in their simulations. In contrast to observations in neutral liquids, ΔH vap increased by 20−30%. It is important to keep in mind that ΔH vap for ILs is calculated relative to ion pairs in the gas phase, where polarization energy per ion is larger than in the bulk. Hence, turning off polarization reduces the cation−anion gas-phase energy more than the bulk energy per ion, therefore resulting in the increased ΔH vap , which is the opposite trend from neutral molecules. (4) The modeling of thermal conductivity λ T of molten salts at high temperatures (1200−1300 K) with a polarizable force field 378 is superior compared to switching off the induced forces or other standard nonpolarizable force fields (Fumi−Tosi potentials) as depicted in Table 5. The nonpolarizable force fields for all three molten salts tend to overestimate the thermal conductivity. Another interesting fact is that the contribution from permanent charges is counteracted by the induced contribution. As we will see, this is also true for other properties.

Structural
Correlations. The comparison between polarizable and nonpolarizable models for the description of IL structure has been discussed in several works. Quite generally, charge scaling reduces the density ρ of the ionic liquid, 60 whereas polarizability increases ρ. 44,52 In particular, the polarizability of the cations seems to be important for the density. Turning off the anionic polarizability does not change ρ so much. 44 Bedrov et al. 52 compared the center-of-mass radial distribution functions g(r) for bulk [ ILs at 393 K as obtained from simulations using a polarizable force field (POL), a nonpolarizable force field where polarization was simply turned off (NP), and two versions of the force-matched nonpolarizable force field (NP-FMp and NP-FMe) where polarization interactions were effectively approximated by two-body interactions. Figure 27 shows the comparison of cation−cation g ++ (r), cation−anion g +-(r), and anion−anion g --(r) radial distribution functions. In all three systems, the largest deviation between POL and NP force fields is observed for anion−anion correlation. Moreover, in case of cation−cation and cation−anion, the mean field adjustment to the nonpolarizable model (NP-FMp and NP-FMe models) leads to almost a complete recovery of g(r)s predicted from the POL model. In contrast, for the anion−anion correlation this provides only minor improvement. The anion−anion g --(r)s with any nonpolarizable model showed stronger ordering compared to predictions using the POL model.
Because of the screening of the electrostatic interactions by the induced dipoles, the repulsion between like charges and the attraction between ions of opposite charges is reduced. 45 Consequently, the structural order is less. 122,135 At short distances, the energy barriers to escape the ion cage are lower. 48,138 At longer distances, the oscillations in the radial distribution functions and hence the charge ordering function Q(r) are flattened, 122,135 which can be fitted beyond 5 Å by and is depicted in Figure 28. Negative values of Q(r) indicate ions of opposite charge, whereas positive values depict an accumulation of like charge ions. As visible in the inset, the flattening parameter 1/σ increases linearly with the amount of molecular polarizability. The same trend was observed by Yan et al., 63 who compared ion−ion correlations in [C 2 mim][NO 3 ] and also found significantly more ordered and longer-range anion−anion correlations predicted by a nonpolarizable model compared to the polarizable one. A very recent study by McDaniel and Yethiraj 79 further confirmed that this phenomenon is quite generic by studying bulk [C 4 mim][BF 4 ] and demonstrating that even charge rescaling does not eliminate the overestimated longrange ordering in the anion−anion correlation. These authors also provided a mechanistic explanation why this effect has to be generic to all ionic liquids with asymmetric cation/anion pairs. McDaniel and Yethiraj suggested that the observed structural deviations are the consequence of a qualitatively different electronic dielectric response at infinite frequency by polarizable and nonpolarizable models. The average reciprocal space Coulomb interactions are plotted in Figure 29 as a function of wavevectors k. In the long-range limit (small wavevector) the Coulomb interactions are approaching unity for nonpolarizable

Chemical Reviews
Review models (independently whether the ionic charges were scaled or not), indicating an infinite frequency dielectric response of ϵ ∞ ≈ 1.0. For the polarizable model, the Coulomb interactions are approaching 0.5 value, which corresponds to ϵ ∞ ≈ 2.0. The difference in screening conditions is the primary reason for enhanced ion structuring at 6−7 Å distances, which is independent of the strength of ion−ion interactions (as demonstrated by results from scaled and nonscaled simulations with nonpolarizable force fields). Therefore, McDaniel and Yethiraj suggest that any simulations using pairwise-additive potentials (i.e., nonpolarizable electrostatic interactions) will contain some artifacts in the predicted liquid structure.

Ion Dynamics.
It is well documented that for all ionic liquids, polarizable models consistently predict faster relaxations (including translational and rotational motions), higher selfdiffusion and ionic conductivity, and lower shear viscosity than the nonpolarizable models. 45,60,77,122,379−381 Bedrov et al. 52 showed that even if a nonpolarizable model is fitted (using the force matching approach, NP-FMe and NP-FMp models

Chemical Reviews
Review discussed above) to reproduce (as far as possible) the forces acting on atoms in the simulations using a polarizable model, the self-diffusion coefficients are still a factor of two to three lower than what was predicted by simulations with the polarizable force field. The authors associated such behavior with the inability of the fitted nonbonded force field to capture directionality of instantaneous atomic forces and as a consequence the correct distribution/fluctuations of those forces. Yan et al. 382 suggested that the faster dynamics observed in simulations of ILs with polarizable models is the consequence of attenuated long-range electrostatic interactions caused by enhanced screening from the induced dipoles, implying that simulations using polarizable models are analogous to simulations with the nonpolarizable model but at higher temperatures. The influence of polarization on the ion transport increases with decreasing temperature indicating that a nonpolarizable force field will have a higher activation energy for transport properties. 44 Interestingly, a temperature effect on the charge scaling factor derived by ESP charges could not be observed for simple molten salts 374 and pyridinium-based ionic liquids. 375 Yan et al. also suggested that polarization affects the vibrational spectrum and hence can influence the hydrogen bond dynamics in ILs. 382 McDaniel and Yethiraj suggested that enhanced dynamics in simulations with polarizable force fields is due to the influence of long-range screening conditions on modulated ion structuring. 79 The enhanced dynamics in ILs also could be obtained by reducing ion charges (through uniform rescaling). Chaban et al. 60 reported that scaling factors ranging between 0.64 and 0.78 (depending on IL and property of interest) can reproduce experimental dynamic and thermodynamic properties. Effectively, by scaling ionic charges, the effective ion−ion interactions are weakened resulting in a higher mobility. The damping of electrostatic interactions can also be explained by a continuum model depicted in Figure 30. Here, the charged particles are immersed in an "inner solvent" of the induced dipoles, represented by a dielectric continuum with a dielectric constant ϵ ∞ . Choosing the correct boundary conditions for the electrostatic potential inside the spheres Φ sp and in the dielectric continuum Φ con , In eq 3.3.6, the radius a of the sphere does not play a role as the distance r between the two charges is much larger than a.
(Although this assumption might not be valid in dense liquids and solid stets where r can be comparable to a). The screening of the electrostatic interaction by a factor of 1/ϵ ∞ in eq 3.3.6 corresponds to a screening factor of f 1/ = ϵ ∞ for each partial charge, 45,62 which is basically the reciprocal value of the refractive index and is in the range of 0.7−0.8, which was already obtained by quantum-mechanical calculations of molecular charges of ion pairs. 6,48,61 The corresponding calculation for the interaction inside the cavity yields The acceleration of the ionic diffusion coefficients due to the damping of the electrostatic interactions can be fitted to 60 The constant exponent c 2 is similar for cations and anions 45,60 in charge scaled simulations. If one correlates the scaling factor f with a molecular polarizability α on the basis of the effective Coulomb energy as reported in ref 45, eq 3.3.8 also holds true for the molecular polarizability. 45,135 Again, c 2 values for cations and anions are similar but less compared to the respective scaling factor f. So far, the molecular polarizabilities were scaled in a uniform way.
Chaban et al. used a pseudofluctuating charge model to describe the charge transfer in NaCl and KCl using classical MD simulations of 1000 ion pairs and switching a certain number of atomic charges to zero or ±1e for neighboring pairs. 383 At a fraction of roughly 90% of fully ionic species, the computational diffusion coefficients agree with experimental data. As average scaling factors of 0.80 or even lower for each ion are needed to reproduce experimental diffusion coefficients, the nonlinearity of the scaling factor becomes evident. This further argues for using the factor f scaled charges only for trajectory production to mimic the average effect of polarizable forces and the full charges should be used for the computation of conductivity as much lower charge scaling factors f̃≪ f would be required to reproduce conductivity if the scaled charges f̃q iβ are also applied for the conductivity computation.
However, the cationic polarizability seems to have a bigger impact on the diffusion coefficients than the anionic α. Turning off the respective induced dipoles for cations resulted in much lower diffusion coefficients, whereas nonpolarizable anions but polarizable cations yielded diffusion coefficients comparable to the full polarizable system. To some extent, this can be explained Depending on the ionic liquid, polarizability effects on the rotation are comparable 45 or less pronounced. 52 The rotational relaxation of the dipole moment may change drastically for small changes of the charge scaling factor as visible in Figure 31. For [C 2 mim][OTf], this turning point appears around a scaling factor of 0.85, which is most common for the scaling. Consequently, one has to be very careful for the evaluation of a meaningful scaling factor as several properties may be quite insensitive and others very sensitive for the variation of the scaling factor. As a consequence, an analogous equation to (3.3.8) does not exist for the scaling factor f but for the molecular polarizability α. However, the corresponding c 2 -value of the cations is higher than for the anions.
3.3.4. Collective Dynamics. Because many dynamic properties scale with the viscosity of the system, the reproduction of this physicochemical property is highly desirable. Both the nonpolarizable GAFF-force field using a charge scaling factor f = 0.8 376 as well as the polarizable APPLE&P force field 44 yield reasonable agreement with experiment for many ILs. The latter force field is also capable of reproducing the corresponding conductivities. Only for some NTf 2 -based ILs the polarizable force field has to be optimized to show less discrepancy with the experimental values.
Polarizable force fields for ionic liquids are also suitable to compute frequency-dependent spectra, e.g., the conductivity spectrum, 382 OKE spectrum, 382 and dielectric spectrum. 45,136,138 Interestingly, charge scaled simulations are also capable of reproducing the frequency-dependent dielectric spectrum over a broad frequency regime. 45 Generally, the dielectric spectrum can be decomposed into a dielectric permittivity ϵ(ω) and a dielectric conductivity ( ) 0 ω ϑ . The latter is only present in the case of charged species. Ionic liquids are an interesting solvent class as the molecular ions possess a net charge and dipole moment. 138 Consequently, these ions contribute to both ϵ(ω) describing the relaxation of the collective dipole rotations and ( ) 0 ω ϑ characterizing the mutual motions of the ions. In simple electrolyte solutions, the neutral solvent, e.g., water, is responsible for ϵ(ω) and the atomic ions, e.g., NaCl, yield a small ( ) 0 ω ϑ usually at higher frequencies ω. In ionic liquids ϵ(ω) and ( ) 0 ω ϑ are not separated in frequency space but overlap making an interpretation on the sole basis of the experimental spectrum difficult. In MD simulations, the collective dipole moment M ⃗ tot (t) = M ⃗ D perm (t) + M ⃗ D ind (t) + M ⃗ J (t) can be split into several contributions: Because the ionic liquid ions are charged species with a molecular charge q i , a molecular dipole moment is not-welldefined. Consequently, one has to choose a reference site. The center of mass r⃗ i (t) of that molecule i turned out to be useful, as it is also the reference site for rotation and translation. 138,316 Autocorrelation functions of the collective rotational dipole moment M ⃗ D perm (t) and the collective induced dipole moment M ⃗ D ind (t) contribute to the dielectric permittivity ϵ(ω) 3.3.12) and the autocorrelation function of the time derivative of the collective translational dipole moment M ⃗ J (t) to the dielectric conductivity ( ) 0 ω ϑ . Interestingly, the relaxation constants for the latter autocorrelation function are very similar to those obtained from the autocorrelation function of neighborhood residency. 138 In other words, the stability of the ion cage and the exchange of its members directly influences the jittering of the  central ion and consequently the high-frequency regime of ( ) 0 ω ϑ . The summations in eqs 3.3.9−3.3.11 can also be split into contributions from the cations and anions. A corresponding decomposition is given in Table 6. The static permittivity ϵ(0) is governed by the relaxation of the permanent dipoles of the anions. Pure induced contributions are negligible for the static value. The cross-correlations between permanent anionic and induced cationic dipoles enhances ϵ(0), whereas the interaction of permanent anionic dipoles with induced anionic dipole moments decreases ϵ(0). The influence of the induced dipoles is more of subtle nature and becomes visible in the frequencydependent dielectric spectrum as it shifts ϵ(ω) contributions of the cations and anions to higher frequencies (as expected from the increased dynamics of polarizable systems). Furthermore, at high frequencies an additional small peak arises from the autocorrelation of the induced cationic dipole moments.
Optical Kerr effect (OKE) spectroscopy yields complementary information to dielectric spectra. First, computational OKE spectra were calculated by Margulis and co-workers 384 using molecular polarizabilities due to the immense computational effort. They showed that the ionic cage is important for the longtime decay of the OKE spectrum in agreement with experiment. 385 Ishida et al. 386 also used molecular polarizabilities to compare vibrational density of state spectra to OKE spectra for [C 4 mim][XF 6 ] with X = (P,As,Sb). They also reported that the OKE spectra are dominated by the molecular reorientation of the cations and anions indicating weak interionic interactions.
In principle, the collective polarizability Π(t) is given by The nonlinear response of the liquid probed by OKE is composed of a zero-time electronic response (containing no information on molecular dynamics) and the nuclear response characterized by the anisotropic polarizability response function using the off-diagonal Cartesian components a ≠ b ∈ (x,y,z) and the Heaviside function θ(t). Its Laplace transform yields the OKE response function χ(ω) depicted in Figure 32 for 3 ]. 382 In the first approximation, χ(ω) can be calculated from nonpolarizable simulations using the time-dependent dipole− dipole tensor T(r⃗ ij (t)) and molecular polarizabilities from a QM frequency calculation. The initial response in Figure 32 at 10 fs is attributed to intramolecular vibrational motion. The following underdamped oscillations with a period of 0.037 ps (∼900 cm −1 ) on top of the damped collective response are assigned to bending motions. The collective response itself is mainly due to librational dynamics. As intramolecular bonding and angle vibrations are similar in the polarizable and nonpolarizable simulation, the corresponding R (3) (t) looks alike. However, at longer times the polarizable R (3) (t) of a stronger signal is visible, indicating faster relaxation of rotational dynamics which is already expected because of the overall faster dynamics in polarizable systems. The OKE response function χ(ω) shows a broad peak up to 300 cm −1 . Here, the discrepancy between polarizable and nonpolarizable response is most prominent and be attributed to diffusive reorientation.
Another THz spectroscopy is solvation dynamics spectroscopy 387 probing the solvent response due to an instantaneous change of the local electric field. This solvent response is coupled to the dielectric spectrum, 136 which is expected because dielectric spectroscopy applies an external electric field, whereas solvation dynamics spectroscopy creates a strong local electric field and both probe the solvent response.  4 ] improved the agreement with experimental data significantly, in particular in the initial time regime and the long-time limit. 137 Only the induced dipoles of first-shell solvent molecules around the excited solute are affected by the solute dipole change, although the electrostatic interactions are visible up to 40 Å. Interestingly, the pure induced contribution to the solvent response function is weak, but the cross-correlation between permanent charges and induced dipoles is strong and counteracting the contribution from the permanent charges only.
Furthermore The number in brackets represent the contributions in a nonpolarizable simulation of the otherwise same force field.

Chemical Reviews
Review shift. 388,140 Using a corresponding nonpolarizable force field nonlinear effects (slowing down the Stokes shift relaxation function) could not be detected. 388,389 Consequently, this fact emphasizes again that induced dipoles do not only affect the dynamics of the liquid system but may also change the fundamental response mechanism to a nonequilibrium event changing the local electric field.

Electrolyte−Electrode Interfaces
Development of electrochemical devices such as batteries and supercapacitors critically depends on understanding of the charge transport, charge storage, and interfacial chemistry at the electrolyte−electrode interfaces. 390−393 Atomistic MD simulations have been extensively applied to study various correlations of electrolytes at charged surfaces, including the structure of electric double layers (EDL) of pure ILs and ionic solutions on flat electrodes, 394−398 dependence of EDL capacitance on electrode voltage and temperature, 50,399,400,401 as well as electrode surface topography, curvature, and porosity. 396,402−410 A more comprehensive discussion of these simulations can be found in several recent review papers. 411−415 In such simulations, the influence of polarization effects is even more complex and can be divided into two issues: (1) Electrode polarization; (2) Influence of electrolyte polarization on EDL structure and properties. In this section, we briefly discuss these specific issues. 3.4.1. Electrode Polarization. The first issue is related to how the charged surface is represented in MD simulation methodology and is coupled to the electrolyte. The most straightforward approach is to confine the electrolyte between two solid surfaces and then apply fixed charges +Q and −Q on the electrodes. The total electrode charge Q is usually homogeneously distributed between N atoms on the electrode surface, with each atom being assigned a charge value of q = Q/ N. As the simulation progresses, the electrolyte will rearrange near each charged surface forming the corresponding EDLs. The electrode potentials as well as the potential difference between electrodes are not known a priori and have to be determined by post analysis of simulation average charge density profiles established in the electrolyte in the direction perpendicular to the electrode surfaces. Because of its simplicity, this approach has been employed in the overwhelming majority of simulations of various ILs and other electrolytes on charged surfaces. However, there is an important drawback in the underlying physical assumption of this approach. In real charge storage devices, the control parameter is the electrostatic potential difference between electrodes, while the charge magnitude and distribution on the electrode surface must reflect the electrode polarization associated with the restructuring of EDL near the electrode surface.
A more natural simulation set up is to confine the electrolyte between two electrodes with a specified applied potential difference between them. This method does not assume a priori any charge distribution on the electrode surface. Instead, the potential difference between two electrodes is constrained while the charges on electrode surfaces are allowed to equilibrate to minimize the electrostatic energy of the system. In this constant potential technique, the electrode charges are assigned either by imposing a fixed electrostatic potential on the electrode surface 416−418 or by constraining the total electrode charge to a desired value and computing the electrode chemical potential as a Lagrange multiplier associated with such a constraint. 419 To reflect a semiconducting character of the electrode, one can introduce an additional energy term approximating the density of states dependence on electrode charge. However, in most simulations using a constant potential approach, the electrodes are treated as conductors and the distribution of charges on such electrodes can be effectively captured using Gaussian smeared charges assigned on electrode atoms. 417,420 The distribution of each charge is controlled by a single parameter, the width of the Gaussian, that in turn can be tuned to quantitatively reproduce the expected behavior of a conductor upon the approach of an external charge. 418 It was shown that the width of the Gaussian distribution of about 0.5 Å is optimal for most cases. 418,421 The electrode charges can also be assigned on a finer grid than the locations of atoms representing the electrode, resulting in a better stability and convergence of the method. 422 The choice of the width of Gaussian distribution is often ambiguous, therefore raising an important question regarding the sensitivity of EDL properties predicted from simulations on the choice of this simulation parameter. Moreover, most MD simulation codes treat electrostatic charges as point distributed and hence utilization of Gaussian distributed charges requires modification of the energy, force, and stress functions. However, Vatamanu et al. 423 recently showed that constant electrode potential simulations that employ Gaussian distributed charges can be modified by adding an energy term proportional to the electrode charge squared while keeping other electrostatic interactions similar to those of point charges (i.e., excluding calculations of Gaussian cross terms of the electrode−electrode and electrode− electrolyte electrostatic interactions). The scheme can be straightforwardly implemented in most MD simulation codes, as it requires no modifications of the standard energy and force evaluation routines and has demonstrated accurate prediction of EDL structure, electrode charge density, and differential capacitance.
Holm and co-workers proposed other efficient methods to handle electrostatic interactions in systems with 2D periodicity and dielectric interfaces. They demonstrated that linear scaling with number of charges can be obtained for electrostatic interactions that include image charges (ICMMM2D and ELCIC methods). 424,425 Subsequently, the induced charge computation (ICC) approach 426 was introduced and which can treat arbitrary curved surfaces and calculate the induced polarization charge located on point charges (instead of Gaussian charges discussed above). It can be applied for systems with arbitrary dielectric contrasts, including metallic electrodes, and any point charge Poisson solver can be used for efficient calculation of Coulomb interactions. The method can be straightforwardly applied to constrain a desired potential difference between electrodes 427 and has been employed in recent studies of charging/discharging kinetics of IL electrolyte in charged nanopores. 428,429 In the constant electrode potential approach, the charges on the surface can (1) evolve in response to EDL restructuring, i.e., due to electrode electronic polarization from electrolyte rearrangement near the surface, and (2) have a heterogeneous distribution in accordance with the structural heterogeneities (surface roughness or defects) inherent to a given electrode surface structure. It is clear that in this case, the electrode polarization and electrolyte restructuring in the EDL are coupled, which is physically a more realistic representation. The question is, how

Chemical Reviews
Review important is the inclusion of this additional physics into simulations for prediction of different properties? Figure 33 compares charge distributions on the basal and prismatic plane graphite surfaces as obtained from MD simulations using a constant applied potential approach. 394 In the constant charge method, the charge distribution would correspond to the delta function. In the constant potential method, a relatively broad distribution of charges for the electrode surface atoms is observed. Note that on atomically flat basal plane graphite, the distributions are relatively symmetric and have a single peak, therefore approximating such distribution with a single average value (as in the constant charge approach) might be a reasonable approximation. However, for an atomically more corrugated prismatic graphite electrode, multiple peaks in the charge distribution are visible, indicating that approximating such distributions with a single charge value cannot be adequate. Taking into account that in many energy storage devices, various chemical modifications of electrode surfaces are considered as one of the routes to design more efficient electrode materials, the ability of MD simulations to capture the details of charge distribution on electrode surfaces is crucial for the accurate prediction of EDL properties.
There have been several works that compared properties predicted from simulations using constant charge and constant potential approaches. One of the most comprehensive discussions and comparisons were presented by Haskins and Lawson, 430 where they systematically investigated properties of [C 2 mim][BF 4 ] at charged surfaces using different methods. Figure 34 shows the differential capacitance obtained for this IL at the graphite electrode. Both, constant charge and constant potential approaches result in very similar average capacitance of ∼4.8 μF/cm 2 . However, the dependence of differential capacitance on the electrode potential is different. The constant charge method produces a well-defined camel-shape dependence, while constant potential approach results in more shallow dependence with less pronounced peaks. At larger magnitudes of electrode potential (>|1|V), the constant electrode potential predicts higher values of the capacitance and weaker dependence on the electrode potential. An analysis of rearrangements in the EDL structure with changes of the electrode potential allowed the authors to make the following conclusions: In the constant electrode potential method, more cations and anions are allowed to pack in the surface ion layer compared to a constant charge method. Thus, the energetic favorability to have both ions in the surface layer makes their separation more difficult and hence higher electrode potentials are required to expel the co-ion from the surface layer. This leads to higher values of differential capacitance for the constant charge potential at low voltages and lower values at higher voltages compared to the constant potential method. Because the average capacitance in both methods is the same, Figure 34 illustrates that in the constant electrode method the EDL restructuring responsible for the charge storage is spread out over a larger potential window compared to the constant charge method. Interestingly, in coarse-grained simulation of IL on a smooth electrode surface or atomistically detailed graphene surface, the choice of electrode treatment (constant charge or constant potential approach) did not show a significant effect on electric double layer structure or differential capacitance. 431 In another work, Wang et al. 432 showed that if the electrode charge distribution is restricted to the locations of atoms representing the electrode surface, then it does not provide a uniformly constant electrostatic potential on the electrode surface leading to both spatial and temporal deviations from the target value of the potential. While in the constant potential method these deviations are relatively minor, for the constant charge method they can be up to 20%, particularly at higher voltages. Finally, Merlet et al. 433 and Vatamanu et al. 50 showed that utilization of the constant charge approach can significantly change the relaxation times in the EDL and lead to orders of magnitude faster charging/discharging kinetics compared to the constant potential method.  The inclusion of electrode polarization becomes even more important for the accurate modeling of aqueous electrolytes with small ions such as Li + that have a compact solvation shell and strongly polarize the electrode upon their close approach. 434 Figure 35a compares the distribution of the surface charges from constant potential simulations of the water-in-salt electrolyte (WiSE) with the surface charges for the electrodes in contact with traditional ILs. A much stronger polarization of the surface by small Li + compared to polarization by the bulky C 2 mim, C 4 mim, and NTf 2 ions of IL leads to a much broader charge distribution with higher probability of more negative charges due to close approach of the Li + to the negative electrode. A positive electrode has a narrower charge distribution because the inner Helmholtz layer primarily consists of the bulky NTf 2 anions, leading to weaker electrolyte polarization similar to the IL double layers as seen in Figure 35b. Thus, the constant electrode approaches are expected to become even less accurate for the WiSE electrolytes than for ILs.
3.4.2. Electrolyte Polarization. Another question is, how important is the electrolyte polarizability near charged surfaces? Besides the already discussed influence of induced polarization on bulk properties of ILs and electrolytes, polarizability of electrolytes can also affect the EDL properties. Haskins and Lawson 430 compared the average magnitude of induced dipoles on [C 2 mim][BF 4 ] ions as a function of their separation from the electrode surface (Figure 36a). They found that surface cations have induced dipoles directed away from the surface, while surface anions have their induced dipoles toward the surface. This alters the interaction of these ions with the corresponding surfaces and between ions in the surface layer(s) and hence creates additional energetic barriers to EDL formation. To quantify how much this affects the differential capacitance, Haskins and Lawson 430 conducted two simulations of the same electrolyte at charged surfaces using a constant potential approach: in the first case, they used a fully polarizable force field to represent [C 2 mim][BF 4 ], while in the second, they turned off all polarization effects in the electrolyte (i.e., set atomic polarizabilities to zero). Figure 36b shows the resulting difference of differential capacitance (ΔC D = C D (pol) − C D (nonpol)). They reported that at lower voltages, electrolyte polarizability tends to reduce the differential capacitance values, while at higher voltages it leads to a slight increase. This behavior can be understood by realizing that induced dipoles oppose the electric field imposed by the electrode surface and hence soften the interactions between ions near the surface.

CONCLUSIONS
Nowadays, polarizable MD simulations of electrolytes, ionic liquids, their mixtures with cosolvents like water as well as at interfaces with gases or solids start to become routine as the corresponding computational power is available. For many applications, in particular for structural investigations, nonpolarizable simulations may yield satisfactory agreement with experimental results. Here, system sizes and simulation trajectory length can be extended due to the cheaper computational costs compared to polarizable simulations, leading to statistically more reliable computational results. On the other hand, local interactions in a fluctuating environment

Chemical Reviews
Review cannot be easily modeled by static electrostatic parameters, i.e., fixed charges, and induced dipoles may be inevitable to make simulation predictions close to experimental data. However, there is an unjustified belief that reliability problems disappear once the induced dipoles are switched on.
The dilemma is depicted in Figure 37. Nonpolarizable MD simulations are usually three to ten times faster than corresponding polarizable simulations depending on the exact computational conditions, e.g., the MD program, its algorithm to compute the induced dipoles at a given configuration, the time step, etc. Nonpolarizable MD simulations have proven to reproduce structural features for a plethora of ionic liquids. However, they usually fail to capture the dynamics of various systems by a significant factor. Possible byways in nonpolarizable simulations, e.g., scaling charges or reparametrizing Lennard-Jones interactions, come at a cost as they only model an averaged effect of the induced dipoles. The collective dynamics and sometimes even the molecular dynamics may be reproduced by these methods, but local interactions like hydrogen bonds suffer from downscaling the charges. Spontaneous, nonequilibrium events, e.g., turning on an electric field or exciting a molecule by a laser beam require an instantaneous response of the ions and solvent molecules which may only be given by their induced dipoles. Furthermore, the local environment may also influence the interactions of the ions. For example, the repulsion between ionic liquid cations is reduced at the interface between liquid and air or at a negative electrode. Because many building blocks of the ionic liquid cations and anions are rather rigid (with the exception of the alkyl tails), fixed partial charge distributions may not be able to capture this effect.
Interestingly, the physicochemical properties, which are generally more time-consuming to compute, are usually those which are better described with polarizable force fields as visible in Figure 37. MD simulations of mixtures are more computationally demanding compared to bulk simulations of pure liquids. Here, the interactions of the ions with cosolvents or solutes profit from a polarizable description which is especially profound for the concentrated electrolytes containing small cations such as Li + , Mg 2+ , or Zn 2+ that are of interest to battery applications. As described above, interfaces are even more complicated and also often require induced dipoles. In the case of dynamics, the very large number of nonequilibrium simulations, e.g., in solvation dynamics spectroscopy, to get reliable statistics is often a computational issue. As a result, the total simulation period of these nonequilibrium simulations is longer compared to a single equilibrium run. Nonequilibrium events may necessitate induced dipoles for nonlinear effects. 142,388 The most common algorithms to model induced dipoles are Drude oscillators and induced point dipoles, which are physical and mathematical approximations of induced polarization, result in similar predictions of structure and dynamics for ionic liquids and electrolytes. The computational cost of these two approaches is comparable and significant progress has been made to reduce it in the past decades. Atomic polarizabilities are available from various sources, but only a few complete polarizable force fields (e.g., APPLE&P) of ionic liquids exists. Using atomic polarizabilities in combination with existing nonpolarizable force fields, e.g., the ionic liquids force field from Canongia-Lopes and Padua, requires a reparametrization of the Lennard-Jones interactions which can be done individually for the atoms or by scaling the Lennard-Jones ϵ as a function of the atomic polarizabilities. However, polarizable force fields for ionic liquids from various groups are currently under development and will be available soon to the computational community, enabling large scale simulations with high performance packages such as NAMD 116 for Drude approaches or Tinker-HP 154 for AMOEBA and other point dipole induced models up to Gaussian-based electrostatics ones.  Liquids" and was awarded the venia docendi in "Theoretical Chemistry". From 2013 to 2017, he was appointed as working group leader of European COST action CM1206 responsible for "Physicochemical Properties of Ionic Liquids and Their Modelling". Christian Schroder's current research interests are based on computational spectroscopy including NMR, dielectric, and THz spectroscopy to study solvation behavior with a particular focus on neoteric solvents such as ionic liquids or their solutions.