Coarse-grained model of a nanoscale-segregated ionic liquid for simulations of low-temperature structure and dynamics

As AA model of [Omim][TFSI], we used the interaction potential proposed by Doherty et al [45]. It is based on the OPLS-AA force field and has reduced ionic charges of ±0.8e to account for polarization effects. The AA model has flexible bonds, expect those of the hydrogen atoms, which are constrained by the LINCS algorithm [46]. Moreover, it allows bond bending and torsional motion.

The simulated AA system contained N = 255 ion pairs, resulting in box lengths L of ∼5 nm. A temperature interval from 400 to 250 K was covered and the lengths of the trajectories ranged from 100 ns at the highest temperature to 3000 ns at the lowest one. The time step of the leap-frog integrator was set to 2 fs. The temperature T was controlled using the velocity-rescale thermostat [47] with a time constant of 1 ps. The pressure time constant for the barostat was set to 2 ps. The long-range part of electrostatic and Lennard-Jones interactions were calculated utilizing the particle mesh Ewald (PME) method [48] with a direct space cutoff of 1.3 nm, a Fourier grid spacing of 0.15 nm⁻¹, and a PME interpolation order of four. Verlet neighbor list updates were performed at least every 10 time steps.

The IBI method was used to derive a CG model from the results of the above AA model at 300 K. For this purpose, the VOTCA 1.5 package [49] was employed and seven beads (A–G) were defined according to the mapping scheme in figure 1, where the position of each bead is given by the center of mass (COM) of the contained atoms.

From an AA simulation at 300 K, we calculated target radial distribution functions (RDFs) for all 28 bead-bead combinations. Moreover, we determined the distributions for all bonds and all angles between neighboring beads. In the CG model, there were no charges on the beads and all bonds were constrained to their mean lengths using the LINCS algorithm because flexible bonds are a major impediment toward the use of larger time steps. We ran 50 IBI iterations with 28 non-bonded potentials and six angle potentials (A–B–D, B–D–D, C–B–D, D–D–D, D–D–E, G–F–G). In each iteration, a system of 255 ion pairs was simulated for 800 ps with a time step of 8 fs. Configurations were stored every 80 fs. The temperature was controlled at 300 K using a Langevin thermostat with a time constant of 1 ps. The start configuration was a mapped AA configuration with a density very close to the average density during the NpT run. The nonbonded and angle distributions were calculated from the trajectory after an equilibration time of 80 ps. Because the iterative method was found to be unstable when using the full IBI update, we scaled it by a factor of 1/4. After each iteration, the Wang–Junghans–Kremer (WJK) correction was applied to the nonbonded potentials. This correction adds a ramp to retain the AA pressure of p^ref = 1 bar [50].

The WJK correction has been described for one component systems and has also been applied to aqueous solutions of urea, but only to update the water–water interactions [51]. Therefore, we shortly describe how to use this correction in systems with multiple different beads. The correction to the interaction potential of beads i and j is a linear ramp

$$\begin{equation}{\Delta}{u}_{\text{ramp},ij}\left(r\right)={A}_{ij}\left(1-\frac{r}{{r}_{\text{cut}}}\right),\end{equation} \tag{ 1 }$$

where A_ij is a parameter and r_cut is the cutoff distance which is 1 nm for all CG interactions. The pressure change that results from applying the ramp to all interactions is approximately

$$\begin{align}\hfill & {\Delta}{p}_{\text{ramp}}\hfill \\ \hfill & \quad \approx \sum\limits _{i=1}^{n}\sum\limits _{j{\geqslant}i}^{n}\left[-\frac{2\pi {N}_{i}{N}_{j}}{3{V}^{2}}{\int }_{0}^{{r}_{\text{cut}}}\mathrm{d}r\enspace {g}_{ij}\left(r\right)\frac{\mathrm{d}{\Delta}{u}_{\text{ramp},ij}}{\mathrm{d}r}{r}^{3}\right]\hfill \\ \hfill & \quad =\sum\limits _{i=1}^{n}\sum\limits _{j{\geqslant}i}^{n}{A}_{ij}\left[\frac{2\pi {N}_{i}{N}_{j}}{3{V}^{2}{r}_{\text{cut}}}{\int }_{0}^{{r}_{\text{cut}}}\mathrm{d}r\enspace {g}_{ij}\left(r\right){r}^{3}\right],\hfill \end{align} \tag{ 2 }$$

where n is the number of bead types, N_i the number of beads of type i, V the box volume and g_ij the RDF of beads i and j. This approximation builds on the assumption that the ramp does not change the RDF significantly. Equation (2) has a unique solution for A_ij only for a single bead type. In the general case, it is necessary to make a choice on how the pressure update is distributed among the interactions. We choose to define A_ij such that each interaction is contributing proportionally to the virial update:

$$\begin{equation}{A}_{ij}=\frac{\frac{{N}_{i}{N}_{j}}{{\sum }_{i=1}^{n}{\sum }_{j{\geqslant}i}^{n}{N}_{i}{N}_{j}}}{\frac{2\pi {N}_{i}{N}_{j}}{3{V}^{2}{r}_{\text{cut}}}{\int }_{0}^{{r}_{\text{cut}}}\mathrm{d}r\enspace {g}_{ij}\left(r\right){r}^{3}}\left(p-{p}^{\text{ref}}\right),\end{equation} \tag{ 3 }$$

where p is the current pressure. The interaction potentials resulting from this IBI approach are used in their tabulated forms.

For the CG model, we simulated a smaller system comprising N = 255 ion pairs in a temperature range from 400 to 170 K and a larger system with N = 1020 ion pairs from 250 to 170 K. The simulation times for the smaller system range from 40 ns at the highest temperature to 4000 ns at the lowest one. For the larger system, the simulation times were chosen accordingly. To integrate the equations of motion, we used the leap-frog stochastic dynamics integrator from GROMACS and a time step of 8 fs. The time constant for the thermostat was set to 2 ps and that for the barostat to 10 ps. The cutoff for the nonbonded interactions is at 1 nm. Neighbor list updates were performed every 10 time steps.

We note that the experimental glass transition temperature of [Omim][TFSI] is between 183 and 189 K [52]. However, we found no indications of glass formation in the studied temperature ranges, reflecting the speedup of dynamics in CG models, see below.

In figure 2, we show the convergence measures Δ_ij for the interactions of bead types i and j, which are given by

$$\begin{equation}{{\Delta}}_{ij}=\sqrt{\frac{1}{{r}_{\text{cut}}}{\int }_{0}^{{r}_{\text{cut}}}\mathrm{d}r{\left[{g}_{ij}\left(r\right)-{g}_{ij}^{\text{ref}}\left(r\right)\right]}^{2}},\end{equation} \tag{ 4 }$$

where ${g}_{ij}^{\text{ref}}\left(r\right)$ is the AA reference RDF and g_ij (r) is the RDF at the current iteration. Furthermore, $\bar{{\Delta}}$ is the arithmetic mean of all Δ_ij . We find only little systematic improvement in the convergence after the twentieth iteration. However, $\bar{{\Delta}}$ continues to fluctuate strongly. The E–E interaction between the tail ends, which is highlighted in figure 2, is mainly responsible for a not ideal match in the structure, and a mismatch in the A–E interaction is correlated with it. This implies that the CG system reproduces the spacial correlations of the AA system to a large extent, but fails to replicate the distribution of the tail ends in some configurations. Below, we will show that tail-end aggregation is key to the nanoscale segregation of the AA and CG models and, hence, particularly large fluctuations may be expected for the E–E interactions.

Zoom In Zoom Out Reset image size

The pressure amounts to 35 bar after 50 iterations. Hence, it is not fully converged to the target value of 1 bar but is substantially improved compared to a pressure of 1400 bar, which is obtained when running IBI without the WJK pressure correction.

We begin by comparing various RDFs of the AA and CG [Omim][TFSI] models at T = 300 K where the IBI was performed. In order to enable a direct comparison, this and the following analyses for the AA model are also performed on the level of beads, i.e., based on the coordinates of the COM of all atoms belonging to a considered bead according to the mapping scheme. In figure 3, we see that the g(r) of both models are in good agreement for most bead types. Some deviations occur for the E–E RDF, which shows a higher first neighbor peak for the CG model than for the AA model, implying that the former is slightly overstructured relating to the arrangement of neighboring tail ends, which will become more evident in the following analysis of nanoscale segregation. Altogether, these findings for the RDFs are consistent with the above results for the IBI convergence.

Zoom In Zoom Out Reset image size

A prominent difference between the AA and CG models is found for the conformation of the alkyl tail. To illustrate the effect, the length distribution of the alkyl tail is shown in figure 4. We see that a D1–E distance of ∼0.73 nm, which corresponds to all but one torsion angles being in gauche conformation, prevails in the CG model, while the alkyl tail conformations are more broadly distributed in the AA model. These findings are a result of the parameterization of the CG force field, where all bond lengths are constrained to allow for a longer time step. On top of the definition of CG beads, this further reduction in the intramolecular degrees of freedom leads to less flexible molecules and, thus, to sharper distributions of structure parameters. Nonetheless, the averages of both distributions agree within the statistical error with mean alkyl tail lengths of 0.70 nm for both the AA and the CG model.

Zoom In Zoom Out Reset image size

For additional comparisons of the AA and CG structures, we further increase the level of chemical abstraction in the analyses. Specifically, we determine the COM positions of the beads (A, B, C) in the imidazolium ring (cation ring), of the beads (D1, D2, D3, E) in the alkyl tails (cation tails), and of the beads (F, G) in the anions. Moreover, we analyze the spatial distributions of the polar groups, comprising the cation rings and the anions, and of the nonpolar groups, i.e., the cation tails, see figure 1, based on the respective COM positions. As the RDFs of all these points are no target functions of the IBI approach, they provide us with an estimate for the agreement of the general structural properties of both models. In figure 5, we compare g(r) for the cation rings, anions, polar groups, and nonpolar groups. We find that the CG model reproduces the RDFs of the AA one reasonably but not fully. Most of the features can be seen at both levels of description, but most first and second neighbor peaks are sharper in the CG than in the AA model. The latter effect cannot be explained by a mismatch of the bead–bead distribution functions in the IBI run. Rather, it shows that a match at the CG resolution is not a guarantee for a match of structure motives on larger scales.

Zoom In Zoom Out Reset image size

To investigate a possible nanoscale segregation of the AA and CG [Omim][TFSI] models, we first inspect some snapshots of the studied systems at T = 180 K and T = 250 K. In figure 6, we see that the polar groups (red spheres) and the nonpolar groups (blue spheres) are not homogeneously distributed in the simulation box, but tend to segregate in both the AA and the CG models. Furthermore, the visual inspection of the configurations suggests that the polar and nonpolar regions grow when the temperature is decreased. Notably, we observe for the larger system comprising 1020 ion pairs that the length scale of the nanoscale segregation at 180 K is comparable to the size of the smaller system with 255 ion pairs, which, thus, may be subject to finite-size effects, see below. Therefore, studies of nanoscale segregation at cryogenic temperatures require not only long simulations to ensure an appropriate equilibration at sluggish particle dynamics but also large systems to avoid size-dependent effects. These conditions cannot be met in AA studies in a reasonable amount of time, demonstrating the need of CG models for such investigations.

Zoom In Zoom Out Reset image size

To ascertain the morphology of the nonpolar regions, we remove all but the E beads from the snapshot of the CG system with N = 1020 ion pairs at T = 180 K and use an isosurface plot of a Gaussian density map to illustrate their connectivity. In figure 6(f), we observe that the E beads and, hence, the tail ends tend to aggregate into worm-like structures with diameters of 1–2 nm. This finding implies that the nonpolar groups do not form isolated clusters in a network of the polar groups [19, 37, 38], but rather that the nanoscale segregation of the studied CG model is characterized by a bicontinuous morphology comprising intertwined nonpolar and polar regions [36].

To develop a first idea of the arrangement of the alkyl tails inside the nonpolar regions, the isosurface plot of the E-rich regions is supplemented by the positions of the other beads of the alkyl tail in figure 7. Explicitly, the D1–D2, D2–D3 and D3–E interbead vectors are included as thin tubes. We see that the end beads (E) constitute the centers of the worm-like nonpolar regions, while the other beads (D1, D2, D3) form a surrounding shell, in which the alkyl tails preferentially point outward, allowing for an integration of the cation rings into the polar regions.

Zoom In Zoom Out Reset image size

To quantify the nanoscale segregation of the models, we compute the static structure factor S(q) of the nonpolar groups according to

$$\begin{equation}S\left(q\right)=1+\rho {\int }_{V}\mathrm{d}r\enspace \mathrm{exp}\left(-\mathrm{i}qr\right)\cdot \left[g\left(r\right)-1\right].\end{equation} \tag{ 5 }$$

Here, g(r) is calculated using the COM coordinates of all nonpolar groups up to distances of half the box length L, limiting the accessible q range to q ⩾ 4π/L. Figure 8 shows the results for the AA and CG models with 255 ion pairs at various temperatures. For both models, S(q) shows a peak in the range q = 12–15 nm⁻¹, corresponding to the next-neighbor distance of the nonpolar beads (r = 0.42–0.52 nm). Of higher interest in the present context is the peak at small values of q, which is associated with the nanoscale segregation [6–10]. When the temperature is decreased, this prepeak strongly grows and shifts to lower q, consistent with the results in a small-wide angle x-ray scattering study on an imidazolium-based IL [7]. While the position of the prepeak is similar in both models, the height shows a weaker temperature dependence in the AA than in the CG description. However, we refrain from a more detailed comparison because, at lower temperatures, the prepeak occurs at the lower end of the accessible q range for the smaller systems so that finite-size effects play a role and the dynamics of the AA model is not fast enough to ensure an ergodic average over the structural heterogeneity, see below. Because of these size and sampling effects, a comparison of the AA and CG results does not inform about the relative performance of both models.

Zoom In Zoom Out Reset image size

Therefore, we focus on the static structure factor S(q) of the CG model to ascertain the temperature-dependent nanoscale segregation. In figure 9, we see that the shift of the prepeak can be followed for the larger CG system with 1020 ion pairs down to 170 K. Thus, we can determine a characteristic length scale of the nanoscale segregation from the prepeak position q_c according to l_c = 2π/q_c. The results for all studied systems are compiled in figure 10. We see that the smaller and larger CG systems show a consistent temperature dependence of l_c. Specifically, the increase of l_c of the small system at high temperatures is continued by that from the large system at low temperatures, where care has to be taken to avoid finite-size effects. In the covered temperature range, the characteristic length scale nearly doubles, resulting in values of ∼3 nm below 200 K. Near room temperature, the results for the CG model are in harmony with l_c = 1.9 nm reported in x-ray scattering work on [Omim][TFSI] [53]. Moreover, the results for the CG model are in reasonable agreement with those from the AA counterpart, although the former yields slightly higher l_c values than the latter. The minor discrepancy may be a result of the reduced flexibility of the alkyl tails in the CG description, which facilitates the alignment of neighboring chains, and of the aforementioned facts that the AA data suffer from finite-size and, due to incomplete sampling of the phase space, imperfect-averaging effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To ascertain the microscopic interpretation of the length scale l_c, we estimate the extension of the nonpolar regions directly from the bead positions. For this purpose, we randomly distribute a large number of points in the simulation box and calculate, for each of the points, the distance to the nearest polar bead. Thus, the largest of these distances in a given configuration yields the radius of the largest sphere that can be placed inside a nonpolar aggregate without enclosing the position of any polar bead. In other words, twice this value is an upper bound for the diameter of the worm-like nonpolar regions. Averaging over many configurations, this method yields characteristic diameters d_c, which increase from ∼1.5 nm near 350 K to a largely temperature-independent value of ∼1.9 nm below 200 K. These diameters d_c are smaller than the corresponding l_c, in particular, at lower temperatures. This confirms that the tails do not form spherical but rather worm-like aggregates so that the value of l_c depends not only on the diameter but also on the length of these structures. Thus, the detailed interpretation of the scattering data depends on the morphology of the nanoscale segregation. Our finding that d_c saturates at a value of ∼1.9 nm at low temperatures is a first hint that the tail arrangement in the nonpolar regions resembles that in cylindrical micelles, where the tail ends meet near the cylinder axis and the tail orientations point radially outward. Specifically, the value d_c ≈ 1.9 nm agrees with the diameter expected for such configuration in our case, which is the sum of two alkyl chain lengths (∼0.7 nm each) and the distance between neighboring E beads (∼0.5 nm).

To explore the arrangement of the alkyl tails inside the nonpolar aggregates in more detail, we analyze their relative orientations. For this purpose, the orientation of the alkyl tails is defined as the unit vector pointing from the D1 bead to the E bead. First, we analyze the relative orientations of alkyl tails in their mutual next-neighbor shells. In doing so, two alkyl tails are considered as neighbors when their COM distance is smaller than 0.75 nm, which is the position of the first minimum in the corresponding RDF for the CG model at 190 K. Moreover, the angle φ between the unit vectors of neighboring alkyl tails is specified based on the value of cos φ, which is readily available from the scalar product. Thus, a probability distribution P(cos φ) ≡ 0.5 is expected in the case of a random orientation of the alkyl tails.

Figure 11 presents the distributions P(cos φ) for the AA and CG models at various temperatures. Independent of the level of description, high values at cos φ = ±1 indicate that neighboring alkyl chains prefer parallel or antiparallel orientations. At 300 K, the probability distributions of both models are in reasonable quantitative agreement. However, P(cos φ) shows a weaker temperature dependence for the AA than for the CG model, where a clear preference for a parallel alignment develops upon cooling.

Zoom In Zoom Out Reset image size

To further investigate the low-temperature behavior, we compare the probability distributions P(cos φ) for the smaller and larger CG systems in figure 12. The relative orientation of neighboring alkyl tails hardly depends on the system size. The high values at cos φ ≈ 1 indicate that parallel alignments clearly dominate the arrangement of the alkyl tails in the nonpolar domains of the CG model at low temperatures.

Zoom In Zoom Out Reset image size

In figure 13, we specify the relative orientation of two alkyl tails in the CG model as a function of their COM distance r. Specifically, we depict the mean value ⟨cos φ⟩ together with the RDF for the COM of the alkyl tails. We see that the parallelism measure ⟨cos φ⟩ and the RDF show maxima at the same distances, indicating that a tendency for a parallel tail alignment is preserved for at least three neighbors along the contour length of the worm-like aggregates. On the other hand, we propose that the intermediary minima are related to the relative orientation of the alkyl tails in the cross-sectional areas of these regions. To discuss this aspect, we again resort to the situation in cylindrical micelles, where chains on opposing sides have antiparallel arrangements and those separated by a quarter turn have perpendicular orientations, considering the alkyl length of ∼0.7 nm and the E–E distance of ∼0.5 nm, these configurations are expected at alkyl distances of, respectively, r ≈ 1.2 nm and r ≈ 0.8 nm, where the minima of ⟨cos φ⟩ are located. Thus, the analysis suggests that the alkyl tails tend to organize in configurations known from cylindrical micelles, though with significantly lower overall order.

Zoom In Zoom Out Reset image size

Finally, we study the MD of the AA and CG models. For this purpose, the mean-square displacement (MSD) of the cations and anions is calculated based on the time-dependent COM coordinates ${ \overrightarrow {r}}_{i}\left(t\right)$:

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

Here, the pointed brackets denote an average over various time origins t₀ along the trajectories. Figure 14 displays the results for the AA and CG models with 255 ion pairs at various temperatures. As expected for viscous liquids, the MSD shows an extended plateau between the regimes of vibrational motion at short times and diffusive motion at long times for sufficiently low temperatures. This indicates that the ions become temporarily trapped in cages formed by their neighbors when the temperature is decreased. While the qualitative behavior is similar for both models, there are distinct quantitative differences. Most notably, the temperature dependence is much higher for the AA model than for its CG counterpart. Moreover, the difference between the cation and anion displacements is smaller for the former than the latter model at lower temperatures. In passing, we note that we found no evidence for a system-size dependence of the MSD in the CG model. Moreover, we point out that the diffusive regime is not fully reached at low temperatures, in particular, for the AA model. Therefore, the above results for the structural properties and the following analysis of the dynamical features at these temperatures should be treated with some reservations.

Zoom In Zoom Out Reset image size

To quantify the differences between the AA and CG dynamics, we determine the self-diffusion coefficients D from fits of the MSD in the diffusive regime with ⟨Δr²(t)⟩ = 6Dt. The results are presented in figure 15. For both models, we observe prominent deviations from Arrhenius behavior, as expected for molecular glass-forming liquids. However, the temperature dependence is much higher for the AA than for the CG model. This effect is further illustrated in the inset, which reveals that the speedup factor D_CG/D_AA increases from ∼10 at T = 400 K to values of more than 1000 near T = 250 K. Similar discrepancies between AA and CG models of glass-forming ILs were reported in previous simulation studies [31, 41]. In the literature, various strategies, e.g., generalized Langevin equation approaches, were proposed to cope with such speedup of dynamics when removing degrees of freedom and smoothing the energy landscapes upon coarse graining [54]. Since the structural properties of ILs are the main focus of the present study, we do not pursue this promising but still challenging route. Nevertheless, our findings that the self-diffusion coefficients of the cations and anions are similar for the AA model in the whole accessible temperature range, whereas the temperature dependence of D is higher for the cations than the anions in the CG model suggest that the prominent nanoscale segregation in the CG model at low temperatures has stronger effects on cation than anion diffusion.

Zoom In Zoom Out Reset image size