Dynamical properties across different coarse-grained models for ionic liquids

Room-temperature ionic liquids (RTILs) stand out among molecular liquids for their rich physicochemical characteristics, including structural and dynamic heterogeneity. The significance of electrostatic interactions in RTILs results in long characteristic length- and timescales, and has motivated the development of a number of coarse-grained (CG) simulation models. In this study, we aim to better understand the connection between certain CG parameterization strategies and the dynamical properties and transferability of the resulting models. We systematically compare five CG models: a model largely parameterized from experimental thermodynamic observables; a refinement of this model to increase its structural accuracy; and three models that reproduce a given set of structural distribution functions by construction, with varying intramolecular parameterizations and reference temperatures. All five CG models display limited structural transferability over temperature, and also result in various effective dynamical speedup factors, relative to a reference atomistic model. On the other hand, the structure-based CG models tend to result in more consistent cation–anion relative diffusion than the thermodynamic-based models, for a single thermodynamic state point. By linking short- and long-timescale dynamical behaviors, we demonstrate that the varying dynamical properties of the different CG models can be largely collapsed onto a single curve, which provides evidence for a route to constructing dynamically-consistent CG models of RTILs.


Introduction
Of the broad variety of molecular liquids, ionic liquids (ILs) stand out for their rich physicochemical characteristics [1,2]. ILs are salts, with a melting point or glass-transition temperature that can reach low temperatures-notably, 'roomtemperature' ionic liquids (RTILs) are in the liquid state at ambient conditions. RTILs are commonly composed of an organic cation and an inorganic anion. ILs play an important role as a solvent in sustainable chemistry, with applications including biomaterials and catalysis [3][4][5][6]. Being conductive, ILs are also strong candidates in electrochemical applications.
The organic cations in ILs often consist of a polar ring group along with nonpolar side chains. The amphiphilic nature of these cations facilitates the formation of nanoscale segregation. Nanoheterogeneous structures have been investigated using both x-ray and neutron scattering [7][8][9][10][11]. Analysis of the structure factor revealed structural inhomogeneity, together with complex heterogeneous dynamics, observed via dielectric spectroscopy [12,13]. Nanoscale segregation increases when the temperature is decreased toward the glass transition, evidenced by shifts in the static structure factor [10]. Low temperatures can also yield specific dynamical effects, such as a breakdown of the Stokes-Einstein-Debye relation [14]. Despite a wealth of studies characterizing the properties of ILs, a clear link between the structural and dynamical properties of these systems, especially close to the glassy regime, remains elusive.
Going up the multiscale-modeling ladder facilitates the study of phenomena occurring at longer length-and timescales. The need for large systems and the associated extensive timescales, due in no small part to the strong electrostatic interactions in ILs, lend themselves to a coarsegrained (CG) description of the system. By lumping together several atoms into super-particles or beads, CG models not only decrease the number of particles to be simulated, but also effectively smooth the underlying free-energy landscape [32][33][34].
A variety of CG approaches have been previously applied to study ILs. Wang and Voth used the force-matching-based multiscale coarse-graining method to coarse-grain 1-n-alkyl-3methylimidazolium ([C n mim] + ) nitrate, for both n = 2 (ethyl) and n = 4 (butyl) [35]. The mapping involved a single bead for the imidazolium ring, one bead per methylene group, as well as a separate bead for the terminal methyl group. Explicit Coulomb interactions were used, where the partial charge of each bead was determined as the sum of the partial charges of the underlying atoms. The study highlighted the interplay between long-ranged electrostatics and the collective shortranged interactions between the side chain groups, along with their role in forming and stabilizing spatially heterogeneous domains.
Bhargava et al proposed a CG model for the family of [C n mim] + cations together with the anion hexafluorophosphate ([PF 6 ] − ) [36]. Three CG beads were devoted to the methylimidazolium atoms to ensure the planarity of the ring, while using one to three CG beads to represent the alkyl chain, depending on its length. Transferability between butyl, heptyl, and decyl side chains was obtained by means of two bead types: terminal and interior beads. A single bead was used to represent the anionic PF 6 . Both intra-and intermolecular interactions were represented by simple functional forms: harmonic potentials for the bonds and bending angles; 9-6type Lennard-Jones potentials for the short-ranged nonbonded interactions; and explicit Coulomb electrostatics. The parameters of the model were optimized to reproduce several properties of [C 4 mim][PF 6 ] under ambient conditions: (i) the mean of reference all-atom (AA) distributions along each order parameter that governs an intramolecular interaction in the CG potential, (ii) the density, and (iii) the surface tension. The authors' analysis highlighted the role of the alkyl chains in the morphology of the liquid, leading to nanoscale ordering, in good agreement with x-ray scattering experiments. The study remains, to date, an excellent landmark for CG models in their capability to reproduce various features of RTILs.
Karimi-Varzaneh et al constructed a model that focuses on the same chemistries, i.e., [C n mim][PF 6 ] with n = {4, 7, 10}, and additionally considered two mappings [37]. While the first mapping resembled the one used by Bhargava et al, the second one notably used a single bead to represent the imidazolium ring. The CG model was parameterized using iterative Boltzmann inversion (IBI), fitting all radial distribution functions (RDFs) while employing a single short-ranged pairwise potential per pair type, i.e., without employing explicit electrostatics. While only targeting structural features, excellent agreement was also found for the surface tension, when compared against both the CG model from Bhargava et al and experiments. Aspects of temperature transferability were probed, with a focus on the change in the major peak positions of the x-ray scattering structure factors. The choice of mapping had a direct impact on the quality of the peak positions, highlighting that the aromatic ring and the alkyl chains impact different regions of the structure factor. Finally, the dynamics of the CG models were characterized. Previous studies employing both pulsed field gradient NMR [38] and AA simulations have demonstrated that the diffusion coefficient (D) of the cation is higher than that of the anion for n = 4. Karimi-Varzaneh et al reported that the ratio D cation /D anion depends on the mapping scheme, the value of n, and temperature. Moreover, a better performance of the second mapping in terms of reproducing the non-Gaussian particle displacement statistics of the AA model highlighted the role of a more detailed representation of the alkyl chains.
More recently, Deichmann and van der Vegt revisited the bottom-up parameterization of [C 4 mim] + with three different anions: [PF 6 ] − , tetrafluoroborate, and chloride [39]. The CG model employed a two-bead representation for the cation: one for the imidazolium ring and one for the alkyl chain. They relied on the conditional reversible work method, which determines the CG parameters from (biased) AA simulations by invoking a thermodynamic cycle. They used Morsetype potentials for the short-ranged nonbonded interactions, in addition to explicit Coulomb electrostatics. The resulting thermodynamic properties of the model-both the liquid-vapor surface tension and mass density-agreed well with AA simulations after an a posteriori correction of the force field. Overall, good agreement was found in terms of the RDFs, although a few outliers were observed, in particular those involving both the imidazolium ring and the [PF 6 ] − anion. Regarding dynamics, it was found that the diffusive speedup with respect to the AA model differs for the cation and anions. We also note the development of other CG models for RTILs using alternative methodologies: Newton inversion [40] (also known as inverse Monte Carlo), relative entropy [41], and graph neural networks [42]. While each of these studies features different parameterization schemes and a variety of resulting properties, they share the common difficulty of representability-simultaneous reproduction of structural, thermodynamic, and dynamical properties.
The dynamical properties of CG RTIL models are of particular interest, since the link between structural and dynamical heterogeneity in these systems remains unclear. Unfortunately, the interpretation of CG dynamics is inherently difficult, as the removal of degrees of freedom from the system results in both a loss of friction and a 'smoothing' of the underlying free-energy landscape [43]. Although generalized Langevin dynamics can be applied to correct for these effects, this approach remains extremely challenging for complex soft matter systems, and also (partially) removes the beneficial speedup provided by the CG model [44]. Instead, many researchers have taken a much simpler approach by determining an effective dynamical speedup factor, e.g., by calibrating against a long-timescale dynamical observable. For specific systems, and in particular when a clear timescale separation exists between characteristic processes, this time-rescaling approach can be highly effective [45,46]. In general, though, the presence of multiple coupled kinetic processes, which may be independently accelerated by the coarse-graining, leads to inconsistent dynamics [47,48]. A similar dynamical inconsistency typically arises in CG models of multicomponent systems with distinct molecular species, e.g., cations and anions in RTILs [37,39]. For bottom-up models, improving the description of the many-body potential of mean force (i.e., the theoretically-ideal CG potential) offers a systematic route toward one important aspect of dynamic (in particular, kinetic) consistency: barrier-crossing dynamics [49]. This link further justifies efforts to improve the structural accuracy of CG models via many-body [42,[50][51][52][53][54][55][56][57][58][59][60] or environmental-dependent [49,59,[61][62][63][64][65][66][67][68] interactions.
Recently, Pal and Vogel investigated the relationship between various dynamical modes and the relevance of spatially heterogeneous dynamics for AA and CG models of [C 4 mim][PF 6 ] [69]. They found that, despite the nontrivial transformation of dynamical processes upon coarse-graining, several relationships between dynamical modes on very different timescales are preserved for both ionic species. In particular, time scales of maximum non-Gaussian and spatially heterogeneous dynamical behaviors collapse onto the same curve when plotted as a function of the structural relaxation time, independent of the value of the ionic charges and, thus, of the electrostatic interactions [70], and consistent with results for various types of viscous liquids [71]. Somewhat related, Douglas and coworkers have put forth an energyrenormalization method for coarse-graining that makes use of a relationship between the structural relaxation time of a liquid and the temperature-dependent activation free energy, which is linked to the configurational entropy of the system [72]. This approach determines a temperature-dependent adjustment to the cohesive energy of the molecular interactions, which has been demonstrated to result in consistent CG dynamics over temperature and frequency (for small-amplitude oscillatory shear molecular dynamics) for liquid ortho-terphenyl [73] and several polymer melts [72,74,75]. These studies demonstrate that by preserving certain short-timescale quantities under coarse-graining of these representative systems, the long-timescale dynamics can be predicted.
In this study, we aim to better understand the connection between certain CG parameterization strategies and the dynamical properties and temperature transferability of the resulting models. We systematically compare five CG models: (a) 'thermo': the model from Bhargava et al [36], largely parameterized from experimental thermodynamic observables; (b) 'thermo * ': a refinement of the thermo model to increase its structural accuracy; (c) 'struct-anabond': a structure-based model focused on reproducing intermolecular distributions; (d) 'struct': a structure-based model that targets both intraand intermolecular distributions; (e) 'struct-260 K': a structure-based model that targets both intra-and intermolecular distributions, parameterized at a lower reference temperature than the struct model.
All five CG models display limited structural transferability over temperature, and also result in various effective dynamical speedup factors relative to a reference AA model. For each model, we quantify the speedup, analyze the relative diffusivity of cations and anions, and study the relationship between vibrational and diffusive motions. In this way, we provide further evidence of a route to constructing dynamically-consistent CG simulation models, in particular for RTILs.

All-atom simulations
The present work makes use of AA simulation trajectories of [C 4 mim][PF 6 ] generated and analyzed in previous works [31,69,70]. Briefly, the force field from Bhargava and Balasubramanian [76] was used, which features cations and anions with partial charges of +0.8e and −0.8e, respectively, harmonic potentials for bond stretching and bond-angle bending as well as dihedral interactions. The simulations were performed using the GROMACS software package [77] for N = 256 ion pairs and a time step of 1 fs. Periodic boundary conditions were applied and the particle mesh Ewald (PME) method [78] was utilized to calculate the Coulomb interactions. At all set temperatures T, the system was first equilibrated at a pressure P = 1 bar to adjust the density, employing the Nosé-Hoover thermostat [79,80] and Parrinello-Rahman barostat [81]. The production runs were carried out in the isochoric-isothermal (NVT) ensemble. We note that, as can be seen in figure 5, the simulations at the lowest temperatures (i.e., below 300 K) may not be fully equilibrated, as the particles do not completely reach the diffusive regime. For this reason, we do not include the data from these temperatures in the dynamical property analysis. Nonetheless, we have used the lowest temperature simulation to construct one of the CG models (see below), based on the structural properties, which presumably are much less sensitive to the full equilibration of these simulations. Additionally, we use the structural properties from these lower temperatures to compare to the properties of the CG models.

Coarse-grained representations and interactions
We employ the CG mapping proposed by Bhargava et al [36] and used previously [69], which represents each imidazolium cation with 4 CG sites and each [PF 6 ] − anion with a single CG site. The imidazolium ring is represented by 3 sites, I1, I2 and I3, mapped to the center of mass of a corresponding group of atoms, as illustrated with the large transparent spheres in figure 1. Note that the I1 and I2 sites overlap, sharing contributions from the 2-carbon of the five-membered ring (i.e., the carbon flanked by two nitrogens). The butyl chain is represented by an additional site, denoted CT, while the anion site is denoted PF. In the following, this mapping is applied to 5 different CG models, each employing the same set of interactions, e.g., nonbonded pairwise interactions between all unique pair types, but with variations in the functional forms and parameters of these interactions. 4 bonded interactions are employed to retain the imidazolium connectivity: I1-I2, I1-I3, I2-I3, I2-CT. Accordingly, each model employs 5 bond-angle interactions: I1-I2-I3, I1-I3-I2, I3-I1-I2, I1-I2-CT, I3-I2-CT. There are no dihedral angle interactions in these models. Finally, 15 pairwise nonbonded interactions, corresponding to all unique pair combinations of the 5 site types, are employed, while excluding intramolecular imidazolium pairs. There are two distinct contributions to each nonbonded interaction: (i) a short-ranged van der Waals-like interaction that is parameterized separately for each model, and (ii) a long-ranged Coulomb interaction that is kept fixed for all models. The Coulomb forces are calculated by mapping the partial charges of the AA model to the CG representation: q = +0.356e, +0.292e, +0.152e, 0e, and −0.8e for the I1, I2, I3, CT and PF sites, respectively.  IBI is a bottom-up method that determines the CG potential that will reproduce a given set of 1D distribution functions through an iterative refinement that assumes independence of the CG potential parameters [82,83]. Given an initial model and a set of reference distributions, each interaction potential, U, is updated at each iteration according to: where k B is the Boltzmann constant, ξ is a scalar order parameter that governs the interaction, and P(ξ) is the Jacobian-transformed distribution function along ξ. For pairwise nonbonded interactions, P(ξ) corresponds to the RDF, g(r), where r is the distance between two particles. α ∈ [0, 1] is a damping factor used to increase the stability of the procedure. Iterative generalized Yvon-Born-Green (iter-gYBG): we also considered an alternative iterative parameterization scheme, using a generalization of the Yvon-Born-Green (YBG) integral equation framework [84,85]. This approach relates a set of structural correlation functions, b, to the parameters of the CG interaction potentials, φ, via a matrix equation, which quantifies the cross-correlations, G, between each CG degree of freedom: For nonbonded, pairwise interactions represented by a set of spline functions, b is directly related to the corresponding RDF, and G characterizes the average cosine of the angle between triplets of particles [86,87]. b can also be expressed in terms of force correlation functions through a connection to the multiscale coarse-graining method [85,88]. This method employs force and structural correlation functions that are determined from a set of AA reference simulations, b AA and G AA , i.e., calculated from a set of AA simulation trajectories mapped to the CG representation, to determine an optimal set of CG parameters φ. Note that if the model derived from this method fails to reproduce the target vector of the equations, i.e., b AA , it implies that the cross-correlation matrix generated by the higher resolution model does not accurately represent the correlations that would be generated by the resulting CG model. This indicates a fundamental limitation of the model representation and interaction set. Nonetheless, equation (2) can be solved self-consistently to determine the interaction parameters φ * that reproduce the target correlations b AA [89]: This approach has been previously denoted as an iter-gYBG method [89][90][91]. When the multiscale coarse-graining model is employed as the initial set of parameters, this procedure has been demonstrated to converge very quickly (e.g., in less than 10 iterations), although it may be less robust than IBI [89,92].

2.3.2.
Models. thermo: as already mentioned in the introduction, the model proposed by Bhargava et al [36] was parameterized in a top-down fashion to reproduce experimental thermodynamic properties, and will thus be denoted as 'thermo' in the following. This model employs analytic functional forms for the CG interaction potentials: harmonic potentials for intramolecular interactions and 9-6 Lennard-Jones potentials for the nonbonded interactions. Initial equilibrium distances and force constants for the intramolecular potentials were determined by fitting the Boltzmann-inverted potentials from AA simulations to the harmonic functional form. These parameters were then refined to reproduce the mean of the 1D distributions along each order parameter governing an intramolecular interaction in the CG model. The Lennard-Jones energies and particle radii were then tuned to reproduce the density and surface tension from experimental measurements at 300 K. It is also reported that the site-site RDFs from AA simulations were taken into account when determining these parameters, although it is not clear to what extent or how this was carried out. thermo * : we performed a refinement of the thermo model to improve its structural accuracy, while attempting to minimally perturb the thermodynamic properties of the model. To achieve this, we applied the iter-gYBG approach to each of the short-ranged nonbonded interactions, but restricted our calculation to relatively short distances between CG sites. More specifically, we calculated the change in the force coefficients up to the distance corresponding to the minimum of each of the potentials [see supporting information (https://stacks.iop.org/JPCM/33/224001/mmedia) and manuscript repository [93] for details]. We also calculated the change in the intramolecular interactions, but ignored these when updating the force field at each step. Previous investigations have shown that taking into account the coupling between intra-and intermolecular interactions can have a substantial impact on the resulting force functions from the g-YBG approach (unpublished). We also performed a linear interpolation between the old and new forces over a range of 0.05 nm preceding the cutoff distance used for the force field calculation, followed by a smoothing of the resulting force functions to avoid unphysical kinks. The procedure was terminated after 6 iterations, since further iterations resulted in numerical artifacts, perhaps due to the strong restriction of the tunable force field parameters. Further details of the iter-gYBG calculations can be found in the description of the parameterization of the struct model below.
struct-anabond: the struct-anabond model was parameterized using IBI to adjust the short-ranged nonbonded interactions to reproduce the RDFs from the reference AA simulation at 300 K, while keeping the intramolecular interactions and long-ranged electrostatics fixed to those used in the thermo model. The procedure was terminated after 160 iterations (see manuscript repository [93] for details). IBI calculations were performed using the VOTCA package [94].
struct: the struct model was parameterized with the iter-gYBG method, using the reference AA simulations at 300 K and the BOCS package [95]. All bond and bending-angle force functions were represented with linear spline basis functions with grid spacings of 0.001 nm and 1 deg, respectively. The short-ranged nonbonded force functions were represented with B-spline basis functions with a grid spacing of 0.01 nm. The initial model was determined by solving equation (3) while using the cross correlations from the AA simulation to compute G (i.e., the initial model was obtained via the multiscale coarse-graining method). The iter-gYBG model was then calculated following the iter-gYBG framework, by solving equation (3) iteratively, using cross correlations determined from the CG simulations at the previous step, until a preset accuracy threshold was achieved with respect to the 1D distribution functions along each order parameter governing a CG interaction (8 iterations in this case). The first 3 iterations scaled the calculated parameter update by a factor of 0.25, while further iterations applied a factor of 0.5. For these calculations, the fixed long-ranged Coulomb interactions were incorporated via the reference potential method [96], which subtracts the contributions to the force correlations, b, from a given set of fixed terms in the force field. To increase numerical stability, we also used the direct Boltzmann inverted forces for each intramolecular interaction as reference, even though we still calculate the optimal forces for these interactions (i.e., we calculate the change in the force parameters, relative to the reference forces). Solution of each set of linear equations and post-processing of the potentials followed previous work [89,92].
struct-260 K: the struct-260 K model was parameterized in the same way as the struct model, but using the reference AA simulations at 260 K. In this case, 7 iterations were required to achieve the accuracy threshold.

Coarse-grained simulations
All CG simulations were performed using the GROMACS simulation package [97]. Each simulation, with the exception of those associated with the thermo model (see further details below), was performed at the same density as the corresponding AA simulation of the same temperature (see table  S2 for the specific densities and comparisons to experiments).
Starting from an equilibrated AA configuration mapped to the CG representation, each CG model was applied to energy minimize the configuration, followed by a 10 ns simulation in the NVT ensemble at the respective temperature. All structural analysis was performed using these simulations. For the dynamical analysis, some models/temperatures required longer simulations. In particular, for the thermo * model, we performed additional 25 ns simulations at 400, 350, 320, and 300 K and additional 100 ns simulations at 280, 270, and 260 K. Similarly, for the struct-anabond, struct and struct-260 K models, we performed additional 25 ns simulations at 280, 270, and 260 K. All CG simulations used the stochastic dynamics integrator with a temperature coupling constant of 2 ps, a 1 fs time step, and periodic boundary conditions. Electrostatic interactions were employed using the PME method [78] with a Fourier grid spacing of 0.10 nm. A cutoff of 1.5 nm was used for both the short-range nonbonded interactions and for the real-space contribution to electrostatic interactions. Configurations were sampled every 0.1 ps during the simulations, and then used for subsequent dynamical analysis. The structural analysis was performed using trajectories parsed to yield configurations every 0.4 ps.
The dynamical properties of the thermo model were taken from previously published simulations [69]. These simulations were carried out at slightly different densities than the AA model (see table S2), corresponding to the equilibrium density of the model at 1 bar. To assess the impact of this change in density on the model properties, we ran 10 ns simulations of the thermo model at each temperature and corresponding AA density, following the protocol described above. All structural analysis presented in this work was taken from these simulations for consistent comparisons with the other models. We note that the slight change in density has negligible impact on the RDFs for nearly all of the simulations. At the lowest temperature (260 K) there are noticeable, albeit very small, deviations in a few of the RDFs. (The full set of RDFs are available in the manuscript repository [93] for comparison). Similarly, we expect only a relatively small impact on the dynamics, although we did not explicitly verify this.
Throughout this manuscript we will describe the dynamical timescales of each CG model in terms of physical time units, in order to achieve a consistent comparison with the AA model. However, as described in the Introduction, the process of coarse-graining results in a lost connection between the CG and AA dynamics. Thus, the reported CG timescales are not meaningful without a proper reference or calibration. However, comparison of relative timescales can assist in assessing if a CG model exhibits consistent dynamics compared with the AA model [47].

CG models
In this work, we investigate a spectrum of particle-based CG models for [C 4 mim][PF 6 ], each parameterized to target specific structural or thermodynamic observables of the underlying system, characterized either with respect to experiments or using AA simulations. We consider five models-denoted thermo, thermo * , struct-anabond, struct and struct-260 K-as described in detail in the Methods section. In brief, the thermo model was previously parameterized to reproduce the experimental density and surface tension [36]. The thermo * model is a refinement of the thermo model, which more accurately reproduces the anion-cation RDFs by adjusting the (very) short-ranged region of the nonbonded interactions, while keeping the remaining interactions fixed. By focusing on adjusting the force functions at short distances, the thermo * model largely retains the thermodynamic accuracy of the original (thermo model) parameterization. In particular, the equilibrium density of the thermo * model at 300 K was determined to be 1.301 g cm −3 , compared with 1.357 g cm −3 for the thermo model, 1.389 g cm −3 for the AA model, and 1.369 g cm −3 from experiments. Additionally, the surface tension of the thermo * model at 300 K was determined to be 38.1 mN m −1 , compared with 39.0 mN m −1 for the thermo model and 42.5 mN m −1 from experiments (see supporting information for calculation details). (Note that, for consistent comparisons in the following analysis, the CG models were simulated at the corresponding AA equilibrium density. However, the dynamical analysis of the thermo model, taken from previously published work [69], employed slightly different densities. See the methods section for details.) The struct-anabond model uses the analytic (harmonic) functional forms and parameters for the intramolecular interactions from the thermo model, while employing tabulated potentials for the short-ranged nonbonded interactions to reproduce the RDFs from reference AA simulations at 300 K. The struct model was also parameterized to reproduce these same RDFs, but additionally employs tabulated intramolecular interactions to reproduce the corresponding 1D distributions along each CG (intra)molecular degree of freedom. Finally, the struct-260 K model is equivalent to the struct model, but parameterized using the reference distributions at 260 K. Figure 2 compares the force functions for each model, for a representative set of interactions. Row (a) presents the only three intramolecular interactions whose corresponding 1-D reference distributions (i.e., calculated from the AA simulations after mapping each configuration to the CG representation) displayed multimodal behavior. It is worth noting that the I1-I2-CT and I3-I2-CT interactions are significantly coupled to one another. As mentioned above, the thermo, thermo * and struct-anabond models (red, yellow and green curves, respectively) employ identical intramolecular force functions, in this case a linear function corresponding to the harmonic interaction potential form. The struct and struct-260 K models (blue and purple curves, respectively), which employ tabulated force functions to match the corresponding 1-D distribution functions, demonstrate overall weaker forces, within the central interaction regime, than their analytic counterparts. The linear extrapolation of these forces to large values at the ends of the interaction range is largely arbitrary and does not impact the accuracy of the resulting distribution (within some reasonable range).
Row (b) presents the three short-range nonbonded interactions corresponding to the RDFs with the largest errors in the thermo model (described further in the following section). The thermo model (red curve) employs a 9-6 Lennard-Jones functional form with potential minima of ≈0.35 kJ mol −1 for the I1-PF and I2-PF interactions and ≈1.2 kJ mol −1 for the I3-CT interaction. The thermo * model (yellow curve) was constructed by adjusting the hard core of each shortrange nonbonded interaction (specifically, for distances shorter than the location of the potential minimum), in an attempt to increase the structural accuracy of the model while minimally perturbing the thermodynamic properties. For example, the I1-PF and I2-PF interactions were refined for distances shorter than 0.584 nm while the I3-CT interaction was refined for distances shorter than 0.435 nm. These adjustments had the largest impact on the I1-PF and I2-PF force functions, as can be seen more clearly in the insets of panels (bi) and (bii). The struct-anabond, struct, and struct-260 K models (green, blue and purple curves, respectively) each employ tabulated force functions to reproduce a given set of RDFs. The resulting forces tend to have significantly larger magnitude than the forces of the thermo model, along with more complicated features, as expected. The struct and struct-260 K models tend to be much more similar than the struct-anabond model, as is clear from figure 2(b). This is almost certainly because the former two models are both determined from just a few iterations (i.e., less than 10) starting from the force-matching-based multiscale coarse-graining model at each temperature, while the struct-anabond model is determined from over 100 iterations starting from the thermo model (see methods section for details). All three of these models demonstrate large pressures in the NVT simulations that prohibit the characterization of the typical thermodynamics properties, e.g., density and surface tension, as is common for structure-based CG models of liquids [98][99][100][101]. Figure 3 presents the 1-D distributions corresponding to the interactions presented in figure 2. Panel (ai) shows that the harmonic potential of the thermo model leads to an accurate description of even the most complicated bond distribution in this molecule, the I2-CT bond, although the fine features (e.g., shoulders) of the distribution are not reproduced. The struct model nearly quantitatively reproduces the distribution, by construction. The small remaining discrepancies are due to numerical difficulties of the iter-gYBG procedure which can occur at the tails of the distribution. This issue could likely be resolved by employing a regularization scheme, but we do not pursue this here, as we do not expect these discrepancies to affect our results. The struct-260 K model matches this distribution for most of its range, but completely neglects the small shoulder at short distances, due to this same numerical issue. Panels (aii) and (aiii) demonstrate that the angular distributions involving the alkyl chain of the cation are quite complex, and cannot be accurately described with the harmonic potential employed by the thermo model. The struct and struct-260 K models reproduce these distributions by construction, albeit with some small discrepancies at the tails of the distributions. Panel (aiv) presents the average mean squared error (mse) for each model over all intramolecular distributions along the CG degrees of freedom that govern an interaction in the CG model (see figure S4 for explicit comparisons of each distribution). For each intramolecular distribution, the mse was calculated over a range encompassing all sampled instances of the corresponding order parameter, e.g., bond distance.

Structural properties
Row (b) presents three sets of RDFs, demonstrating the improvement of the thermo * model, relative to the thermo model, in terms of the description of the anion-cation packing. The remainder of the RDFs show similar behavior between the two models, analogous to panel (biii). The struct-anabond, struct, and struct-260 K models reproduce all RDFs by construction. Panel (biv) presents the average mse for each model over all 15 RDFs (see figures S5 and S6 for explicit comparisons of each RDF). For each RDF, the mse was calculated from 0 to 1.2 nm.
We also analyzed the structural properties of the models as a function of temperature. Figure 4 characterizes the structural accuracy of each model, with respect to the reference AA model, in terms of the average mse of the RDFs (panel (a)) and the full width at half maximum (FWHM) of the first peak in the I1-PF RDF (panel (b)). Although it is somewhat difficult to see on the logarithmic scale, panel (a) demonstrates a decrease in structural accuracy for the thermo and thermo * models upon cooling, despite these models being parameterized from reference data at 300 K. On the other hand, the structure-based models (struct-anabond, struct, struct-260 K), with much lower errors overall, show divergence of the error away from the reference temperature of parameterization (300 K for struct-anabond and struct; 260 K for struct-260 K). Panel (b) clearly demonstrates an improved accuracy of the thermo * model's description of the I1-PF RDF over the entire temperature range, relative to the original thermo model. The struct-anabond and struct models very accurately reproduce the AA FWHM for temperatures near the reference temperature of parameterization, but demonstrate increasing error for the much higher temperatures. The struct-260 K model provides a similarly accurate description of the I1-PF FWHM, with slightly larger errors at the highest temperatures, further indicating a limited range of temperature transferability for these structure-based models.

Dynamical properties
To investigate the molecular dynamics of the AA and CG models, we calculate the mean-square displacement (MSD) of the cations and anions based on the time-dependent bead coordinates r i (t), where the pointed brackets denote an average over various time origins t 0 . From the long-time diffusive regime of the MSD, we determine the self-diffusion coefficients D of the cations and anions by fits to r 2 (t) = 6Dt. We note that the AA model generates diffusion constants in good agreement with experiments (see table S3 in the supporting information). Figure 5 presents the MSD of the AA model and two representative CG models (the thermo and struct models) at various temperatures. Typical of viscous liquids, a regime of vibrational motion at short times and a regime of diffusive motion at long times are separated by a plateau, which becomes longer upon cooling [102]. This plateau regime occurs when the ions are temporarily trapped in cages formed by their neighbors and, as a consequence, display only local rattling motions. While the MSD is qualitatively similar for the three models, we observe quantitative differences. In particular, the diffusive motion has a higher temperature dependence in the AA model than in the CG models. This is consistent with previous work [72], and is due to the inherent difference in energy-enthalpy compensation in the CG models [101]. Panel (a) demonstrates for the AA model that the cations (dashed curves) show higher displacements than the anions (solid curves) in the diffusive regime, consistent with experimental results [38]. On the other hand, the thermo model (panel (b)) exhibits relatively similar diffusivity of the ions, with slightly larger cation displacement at low temperatures but a switch-over to lower cation displacement at high temperatures. The thermo * model demonstrates qualitatively similar behavior (see figure S25). The struct model (panel (c)) provides a slightly more consistent description of the relative ion diffusivity, albeit with an underestimate of the cation versus anion displacements. The structanabond and struct-260 K models demonstrate qualitatively similar behavior (see figure S25). We note that the diffusive regime is not fully reached for the AA model at the lower temperatures of the studied range, interfering with a determination of reliable diffusion coefficients D. Therefore, we restrict the following quantitative analysis to temperatures T 300 K.
To quantitatively compare the diffusive behavior of the different models, we calculated the diffusion coefficients for the anions and cations, using the MSD curves presented in figure 5. Figure 6 presents the temperature-dependent ratio of these coefficients for each CG model, relative to the AA model: These ratios, presented for the cations and anions separately in panels (a) and (b), respectively, represent an effective speedup factor for connecting the CG and AA dynamics, at a particular thermodynamic state point and for a specific molecular species. As one may expect based on the smoothening of the energy landscapes upon coarse-graining, the speedup factors D s are significantly larger than unity. All CG models also demonstrate an increase in D s upon cooling, consistent with the above observation of the difference in temperature dependence of the MSDs between the AA and CG models. By refining the structural accuracy of the thermo model, the thermo * model results in slower dynamics overall, and displays the smallest speedup factors of all the models. Similarly, by refining the intramolecular structure relative to the struct-anabond model (which displays the largest D s values), the struct model also exhibits slower dynamics. On the other hand, the change in the reference state of parameterization between the struct and struct-260 K models appears to have minimal impact on the effective dynamical speedup. Although the thermo * model yields the best absolute reproduction of the AA dynamics, according to this metric, the more relevant quantity for assessing the dynamical consistency of each CG model is the relative mobility of the distinct ionic species. Panel (c) of figure 6 presents the ratio of speedup factors between the anions and cations, D ani s /D cat s , at each thermodynamic state point. In this plot, a value of unity indicates consistent CG dynamics at a single temperature, since the anions and cations experience the same speedup upon coarse-graining. The structure-based models (struct-anabond, struct, and struct-260 K) demonstrate slightly more consistent dynamics, i.e., lower D ani s /D cat s , over the entire range of temperatures. Interestingly, the struct-anabond model appears to show the most consistent dynamics, despite having the largest absolute speedup factors, although the scattering of the speedup ratio is slightly larger than for the other structurebased models. The thermo and thermo * models display larger speedup ratios, which increase slightly upon cooling. By refining the structural accuracy relative to the thermo model, the thermo * model does result in slightly more consistent dynamics for all but the highest temperature.
To ascertain the origin of the observed differences between AA and CG dynamics, we relate long-time and shorttime motions for both levels of resolution. We anticipate that smoother energy landscapes and higher temperatures result in not only faster diffusive motion but also larger amplitudes of cage rattling motion in the plateau regime. Elastic [103] and elasticity-based [104] theories of glass-forming liquids propose relations between the structural relaxation time and the cage rattling amplitude, which are consistent with simulation results [105][106][107]. Figure 7 presents the diffusion coefficients, D, versus the amplitude of the cage rattling motion, which we characterize by the MSD in the plateau regime, r 2 p ≡ r 2 (t = 1 ps) . It is apparent that the very different dynamical behaviors of the variety of models characterized above collapse to a similar D(r 2 p ) dependency, where the curves for the cations and anions are somewhat shifted relative to each other. Overall, the structure-based models (struct-anabond, struct, struct-260 K) demonstrate a tighter collapse than the thermo and thermo * models, relative to the AA behavior, although these models also show deviations at small r 2 p (i.e., lower T). The collapse of CG dynamics implies that a smoothening of the energy landscape leads to a facilitation of the diffusive motion via a softening of the local cages. This may provide a promising route for constructing dynamically-consistent CG models, as recently suggested in a study on ILs [69] and put into practice for liquid ortho-terphenyl [73] and several polymer melts [72,75]. More specifically, the latter studies argue that a quantitative prediction of long-time dynamics is possible when r 2 p is preserved under coarse-graining. This provides a tractable method for CG parameterizations, since r 2 p can be efficiently determined from (relatively) short AA reference simulations.

Conclusions
This work employed five distinct structure-and/or thermodynamic-based coarse-graining parameterization strategies for modeling the room temperature IL [C 4 mim][PF 6 ]. All five models displayed limited transferability in accurately representing structural properties over a range of temperatures, although the absolute error of the structure-based models was much lower overall, by construction. These structural discrepancies might be systematically addressed via recently proposed approaches for estimating the entropic contribution to the many-body potential of mean force (i.e., the theoretically-ideal CG potential) [108][109][110][111].
The focus of our study was to better understand the link between the various parameterization strategies and the dynamical properties of the resulting models. We quantified the dynamical speedup for each model and for each ionic species via the ratio of diffusion constants, relative to the AA model, D s = D CG /D AA . All five models showed a systematic increase in D s , for both cations and anions, with decreasing temperature, indicating dynamical inconsistency relative to the AA model (i.e., a lack of transferability). At a given temperature, the structure-based models tended to provide a more consistent description of the relative anion-to-cation behavior. Interestingly, the accuracy with which the model described the intramolecular distributions of the cation and also the reference parameterization temperature for the structure-based models appeared to play a limited role in the resulting diffusive properties of the model. Moreover, our results demonstrate that the absolute speedup is independent of the dynamical consistency of the model: in this case the model with the largest speedup provided the most consistent description of anion-cation relative diffusion.
We also investigated the relationship between vibrational and diffusive motions for each model. Following previous work on theories of glass-forming liquids, we assessed the relationship between the structural relaxation time and the cage rattling amplitude via the diffusion constant and the value of the mean-squared displacement in the plateau regime, respectively. Despite their range of parameterization strategies, we found a collapse of CG model behavior in terms of this relationship (figure 7). The structure-based models demonstrate a slightly tighter collapse, relative to the AA behavior, although clear discrepancies still arise for the lowest temperatures considered. While these models reproduce a set of low-dimensional structural distribution functions by construction, the need for iterative refinement relative to the initial (force-matching-based multiscale coarse-graining) model indicates fundamental limitations of the chosen CG interaction functions to accurately describe higher-order structural correlations [89,92,112]. This then motivates investigations that probe the relationship between improved representations of the many-body potential of mean force and the resulting dynamic consistency of the CG model [49]. Overall, the results in this work provide evidence of a clear route to constructing dynamically consistent CG simulation models of RTILs via a temperature-dependent matching of AA vibra-tional amplitudes, as recently carried out for model glass formers [72,73,75].

Data availability
Additional methodological details and results can be found in the supporting information of this work. We have also compiled a repository [93] containing all the models (in GRO-MACS format) developed in this work (i.e., the thermo * , struct-anabond, struct, and struct-260 K models), some parameter files for the calculations performed with the VOTCA and BOCS packages, and also some of the raw data used to construct the manuscript figures. Additional data can be obtained from the corresponding author upon request.

Author contributions
JFR, KK, TB, and MV designed the work. JFR, SK, SW, and TP performed calculations for the work. JFR, SK, KK, TB, and MV analyzed the results. JFR, SK, KK, TB, and MV wrote the manuscript.