Dynamical properties across different coarse-grained models for ionic liquids

The present work makes use of AA simulation trajectories of [C₄mim][PF₆] 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.

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 ${\left[{\mathrm{P}\mathrm{F}}_{6}\right]}^{-}$ 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.

Zoom In Zoom Out Reset image size

2.3.1. Methods

Iterative Boltzmann inversion (IBI): 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:

$$\begin{equation}{U}^{\left(i+1\right)}\left(\xi \right)={U}^{\left(i\right)}\left(\xi \right)-\alpha {k}_{\mathrm{B}}T\enspace \mathrm{ln}\enspace \frac{{P}^{\left(i\right)}\left(\xi \right)}{{P}^{\left(\mathrm{r}\mathrm{e}\mathrm{f}\right)}\left(\xi \right)}\enspace ,\end{equation} \tag{ 1 }$$

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:

$$\begin{equation}\boldsymbol{b}=\boldsymbol{G}\boldsymbol{\phi }.\end{equation} \tag{ 2 }$$

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]:

$$\begin{equation}{\boldsymbol{b}}^{\text{AA}}=\boldsymbol{G}\left({\boldsymbol{\phi }}^{{\ast}}\right){\boldsymbol{\phi }}^{{\ast}}.\end{equation} \tag{ 3 }$$

This approach has been previously denoted as an iter-gYBG method [89–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

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].

In this work, we investigate a spectrum of particle-based CG models for [C₄mim][PF₆], 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⁻³, compared with 1.357 g cm⁻³ for the thermo model, 1.389 g cm⁻³ for the AA model, and 1.369 g cm⁻³ from experiments. Additionally, the surface tension of the thermo* model at 300 K was determined to be 38.1 mN m⁻¹, compared with 39.0 mN m⁻¹ for the thermo model and 42.5 mN m⁻¹ 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).

Zoom In Zoom Out Reset image size

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⁻¹ for the I1–PF and I2–PF interactions and ≈1.2 kJ mol⁻¹ for the I3–CT interaction. The thermo* model (yellow curve) was constructed by adjusting the hard core of each short-range 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–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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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),

$$\begin{equation}\langle {r}^{2}\left(t\right)\rangle =\left\langle \frac{1}{N}\sum\limits _{i}^{N}{\left[{\boldsymbol{r}}_{i}\left(t+{t}_{0}\right)-{\boldsymbol{r}}_{i}\left({t}_{0}\right)\right]}^{2}\right\rangle ,\end{equation} \tag{ 4 }$$

where the pointed brackets denote an average over various time origins t₀. 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²(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 struct-anabond 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.

Zoom In Zoom Out Reset image size

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: D_s = D^CG/D^AA. 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}_{\text{s}}^{\text{ani}}/{D}_{\text{s}}^{\text{cat}}$, 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}_{\text{s}}^{\text{ani}}/{D}_{\text{s}}^{\text{cat}}$, 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 structure-based 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.

Zoom In Zoom Out Reset image size

To ascertain the origin of the observed differences between AA and CG dynamics, we relate long-time and short-time 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–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}_{\text{p}}^{2}\equiv \langle {r}^{2}\left(t=1\enspace \mathrm{p}\mathrm{s}\right)\rangle$. It is apparent that the very different dynamical behaviors of the variety of models characterized above collapse to a similar $D\left({r}_{\text{p}}^{2}\right)$ 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}_{\text{p}}^{2}$ (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}_{\text{p}}^{2}$ is preserved under coarse-graining. This provides a tractable method for CG parameterizations, since ${r}_{\text{p}}^{2}$ can be efficiently determined from (relatively) short AA reference simulations.

Zoom In Zoom Out Reset image size