p–p stacking of imidazolium cations enhances molecular layering of room temperature ionic liquids at their interfaces†

The interfacial structure of room temperature ionic liquids (RTILs) controls many of the unique properties of RTILs, such as the high capacitance of RTILs and the efficiency of charge transport between RTILs and electrodes. RTILs have been experimentally shown to exhibit interfacial molecular layering structures over a 10 Å length scale. However, the driving force behind the formation of these layered structures has not been resolved. Here, we report ab initio molecular dynamics simulations of imidazolium RTIL/air and RTIL/graphene interfaces along with force field molecular dynamics simulations. We find that the p–p interaction of imidazolium cations enhances the layering structure of RTILs, despite the electrostatic repulsion. The length scales of the molecular layering at the RTIL/air and RTIL/graphene interfaces are very similar, manifesting the limited effect of the substrate on the interfacial organization of RTILs.


Introduction
Room temperature ionic liquids (RTILs) constitute a class of salts, which are in the liquid state at room temperature. Intermolecular ionic interactions between their cations and anions differentiate RTILs from non-ionic liquids, and render RTILs unique solvents: 1,2 RTILs are environmentally friendly solvents thanks to their low vapor pressure. Moreover, owing to their ionic nature, RTILs have potential applications in supercapacitors, as battery electrolytes, and as catalysts for generating graphene sheets from graphite. [3][4][5] The interfacial properties of RTILs govern their performance in supercapacitors and in the graphene exfoliation process. Unveiling and understanding the local structure of RTILs at their interfaces is therefore crucial for designing and optimizing the physical and chemical properties of RTILs.
Since computational approaches based on the force field/ ab initio molecular dynamics (MD) simulation allow us to visualize the molecular structures of RTILs, they have been frequently used for connecting both the bulk 22-37 and interfacial structures [38][39][40][41][42][43][44] of RTILs with their physical and chemical properties. The non-polarizable force field MD (NP-FFMD) simulation was the first to be applied to RTILs. NP-FFMD can reproduce the structural properties of RTILs such as the radial distribution function and density, whereas dynamical properties such as the diffusion constant and viscosity were poorly reproduced, because of the slow RTIL dynamics in the NP-FFMD simulation. [22][23][24][25] To improve the description of the dynamical properties of RTILs, polarizable force field MD has been applied for RTILs. 26,[33][34][35][36][37] The molecular polarizability of RTILs in the polarizable force field MD (P-FFMD) simulation screens the electrostatic interactions, accelerating the RTIL dynamics. 26 Furthermore, the interfacial structure of RTILs at the RTIL/air interface was examined using P-FFMD. 42 The comparison of P-FFMD and NP-FFMD shows that the molecular polarizability randomizes the specific orientation of the imidazolium cation. 42 However, both FFMD simulations are known to have serious drawbacks for balancing various types of intermolecular interactions such as the electrostatic interactions, p + -p + interactions, and hydrogen bond interactions. 27 Both the inclusion of dispersion corrections 45 and the interplay between p + -p + interactions and hydrogen bonding 46 are of importance to stabilize the cation-cation arrangements in the RTIL bulk system. However, the role of p + -p + interactions at the interface of RTILs is still not clear.
Here, by using the ab initio MD (AIMD) simulation with the van der Waals corrections and comparing AIMD data with NP-and P-FFMD data, we explore the RTIL layering at the RTIL/air and RTIL/graphene interfaces. Our simulation indicates that the RTIL interfacial structure results from the competing effects of electrostatic forces (causing layering), the molecular polarizability (disordering the interface) and p + -p + interactions (giving rise to additional layering and ordering). Furthermore, the very similar length scale of molecular layering at the RTIL/air and RTIL/graphene interfaces reveals that the molecular layering is governed by the intrinsic RTIL interactions and is not strongly influenced by the substrate.

Ab initio MD simulations
We employed the Becke-Lee-Yang-Parr (BLYP) 47,48 and Perdew-Burke-Ernzerhof (PBE) 49 exchange-correlation functionals for the AIMD simulation at the 1,3-dimethylimidazolium chloride [MMIM + ][Cl À ]/air interface and the BLYP exchange-correlation functional for the AIMD simulation at the [MMIM + ][Cl À ]/ graphene interface. We used the triple-zeta valence plus two polarization (TZV2P) basis sets and the real-space density cutoff of 400 Ry. The core electrons were described by the Goedecker-Teter-Hutter pseudopotential. 50 The van der Waals correction was included via the Grimme's D3 method. 51 Fig. 1(a). The time step for integrating the equation of motion was set to 0.6 fs. The AIMD simulation employing the QUICKSTEP 52 method was performed by using the CP2K code. 53 We started the BLYP + D3-AIMD simulation at 450 K using four different initial conformations for the RTIL/air interface, while we ran the BLYP + D3-AIMD simulation at 450 K with one initial conformation for the RTIL/graphene interface. Note that since the melting temperature of [MMIM + ][Cl À ] is 398.15 K, 54 our system is in liquid phase. For these initial configurations, we conducted 5 ps AIMD simulations in the NVT ensemble for equilibrating the systems. We used the canonical sampling through a velocity rescaling (CSVR) thermostat 55 to control the temperature. Sequential 425 ps AIMD runs were performed for sampling the AIMD trajectories, and finally we obtained total 123.4 ps AIMD trajectories for the RTIL/air interface and 28.6 ps AIMD trajectory for the RTIL/graphene interface. The snapshots of the simulated RTIL/air interface and RTIL/graphene interface are shown in Fig. 1(b) and (c), respectively.
For the PBE + D3-AIMD simulation at the RTIL/air interface, we used the same MD simulation conditions, and the same initial configurations as the BLYP + D3-AIMD simulation. We used the four different initial configurations for the PBE + D3-AIMD simulation. These initial configurations were the same as used in the BLYP + D3-AIMD simulation. For the simulation at 450 K, we performed 5 ps PBE + D3-AIMD runs for equilibrating the systems and then performed over 27 ps AIMD runs for obtaining the MD trajectories. The total 109 ps AIMD trajectories were used for analyzing the data. More detailed information on the simulation is given in the ESI. †

Force field MD simulations
For the NP-FFMD simulation of the [d-MMIM + ][Cl À ]/air interface, we employed the OPLS-AA-based force field model developed by Canongia Lopes, Deschamps and Pádua. 56 The C-D (C-H) bond lengths were fixed by using the SHAKE algorithm. A time step of 1 fs was used for integrating the equation of motion. We used the same simulation cell size and the same number of ion pairs in the NP-FFMD run. Periodic boundary conditions were used. We ran a 5 ns NP-FFMD simulation in the NVT ensemble to equilibrate the systems at 450 K. The temperature was controlled by the Andersen temperature coupling scheme in the production run. 57 Subsequently, we performed a 38 ns MD run at 450 K to obtain the trajectories, which were used for analyzing the data. The NP-FFMD simulation was performed by using the PMEMD module of the Amber14 package. 58 For the P-FFMD simulation of the [d-MMIM + ][Cl À ]/air interface, we used the SANDER module of the Amber14 software package. 58 The parameters for the intermolecular interactions, point charges, and Lennard-Jones parameters were the same as those used in the NP-FFMD simulation. In the P-FFMD simulation, we added the atomic polarizability models to the NP-FFMD simulation. One model was the Amber polarizability model which is designed to apply for any atom species. 59 The other was the polarizability model optimized by the imidazolium room temperature ionic liquids (the Voth model 26 ). The electrostatic interactions were calculated within the Applequist model. 60 The target temperature was set to 450 K. The other simulation conditions were the same as the NP-FFMD simulation. We ran 5 ns P-FFMD runs to equilibrate the systems. Subsequently, we performed a 20 ns P-FFMD run for the Amber polarizability model and a 19 ns P-FFMD run for the Voth polarizability model at 450 K to obtain the P-FFMD trajectories, which were used for analyzing the data.

Geometry-definitions of stacked, displaced and T-shaped conformations
We analyzed the p + -p + configurations of cations from the MD trajectories using distinct geometry-based definitions. These geometry-based definitions are based on the relative configuration of a pair of [d-MMIM + ], for which we calculated the three vectors, n 1 ! , n 2 ! , and R 1 R 2 ! ; n 1 ! and n 2 ! are the normal of the planes formed by the three carbon atoms of the imidazolium rings of two cations, while R 1 R 2 ! is the vector pointing from R 1 to R 2 where R 1 and R 2 are the centers of mass of three carbon atoms of the imidazolium rings. These are schematically shown in Fig. 2.
We evaluated the p + -p + interaction configuration pairs of [d-MMIM + ] using the following conditions: for the stacked p + -p + configuration (1), we set the criteria of j o 301, The positions of these pairs are represented by the midpoints of R 1 and R 2 . By using these geometry-based definitions, we quantified the propensity of the different cation-cation conformations through the p + -p + interactions.

Layering structure of RTILs at the RTIL/air interface
We calculated the number densities of [MMIM + ] and [Cl À ] (denoted by r + and r À , respectively) at the RTIL/air interface. The axial profiles of the number densities, r + , r À , and the difference of the densities, (r + À r À ), are plotted in Fig. 3(a-d).
These indicate that near the RTIL/air interface, the anion-rich and cation-rich layers appear repeatedly. The (r + À r À ) data also show that the alternating anion-rich and cation-rich layers decay towards the bulk, consistent with the previous NP-FFMD simulation studies. 43,44 To quantify the decay of (r + À r À ) between AIMD, P-FFMD and NP-FFMD simulations, we plot the normalized area of the positive and negative (r + À r À ), that is, we integrated the function of (r + À r À ) between adjacent zero crossings, and took the absolute value of the integrals. This results in a discrete function DS, which represents the successive areas filled with red and blue colors in Fig. 3(a-d). DS is plotted as a function of the peak maxima in Fig. 3(e). This figure shows that DS decays more quickly in the P-FFMD simulation than in the NP-FFMD simulation. This indicates that the molecular polarizability, which is accounted for in P-FFMD but not in NP-FFMD, substantially randomizes the RTIL layering structure near the RTIL/air interface. On the other hand, the AIMD simulation indicates a much slower decay of DS than the P-FFMD simulation, although both the AIMD and P-FFMD simulations include the effects of the molecular polarizability. Surprisingly, (r + À r À ) is similar for the AIMD and NP-FFMD simulations. This raises the question why the AIMD simulation predicts more extensive RTIL layering than the P-FFMD simulation. DS can be described well by a phenomenological function Df (z) = f 0 (1 À exp(À(z À c)/l))exp(Àz/x), where the term exp(Àz/x) represents the decay of (r + À r À ) with a length scale of x (Lorentzian model 10 ). The term (1 À exp(À(z À c)/l)) represents the rapid density increase at the interface and different surface activity of the cation and anion; c and l denote the coordinate shift from the Gibbs dividing surface and the decay constant of the surface modulation (l o x), respectively. f 0 is the amplitude of Df (z). We obtained l from global fitting for all the NP-FFMD, P-FFMD, and AIMD simulations, while x, c, and f 0 were allowed to vary between the different simulations. The normalized fit curves Df (z) are also plotted in Fig. 3(e). The decay constant of the density oscillation in the AIMD simulation amounts to x = 10.4 Å, about twice as large as that inferred from the P-FFMD simulation (x = 6.2 Å for Amber polarizability model and x = 4.2 Å for the Voth polarizability model), but remarkably similar to the x = 9.6 Å of the NP-FFMD simulation.

p + -p + conformations at the RTIL/air interface
To elucidate the mechanism for enhanced RTIL layering in the AIMD simulation against the structure randomization due to the molecular polarizability, we consider the effect of p + -interactions on the RTIL interfacial structure. In fact, the ion pairs interact substantially through the p + -anion interactions in the AIMD simulation, while they do not interact in the FFMD simulation. [29][30][31] However, quantifying the p + -p + interactions of [MMIM + ] is not straightforward, since these p + -p + interactions mediated by [Cl À ] can have several different conformations -stacked, displaced, and T-shaped conformations. [61][62][63] The axial profiles of the p + -p + conformations calculated from the AIMD trajectory, P-FFMD trajectory using the Amber polarizability model, and NP-FFMD trajectory are displayed in Fig. 4(a)-(c). The AIMD simulation shows that at the RTIL/air interface, the stacked conformation is dominant and substantially more prevalent at the interface than in the bulk. In contrast, the displaced and T-shaped conformations are not prominent at the interface. This trend differs from that inferred from the NP-FFMD Fig. 3 (a-d) Axial profiles of the RTIL number densities at the RTIL/air interface. We obtained the data by averaging the number density at the two RTIL/air interfaces in the slab models. The origin point (z = 0 Å) is set to the Gibbs dividing surface. The areas of (r + À r À ), DS(z), are filled with red (blue) color, when the peak is positive (negative). (e) Comparison of DS obtained from different models. The fit function is given by Df (z) = f 0 (1 À exp(À(z À c)/l))exp(Àz/x). The data were normalized to the maximum value of the fits. The fit parameters are summarized in Table 1. and P-FFMD simulations; these three conformations show similar distributions and no preferable stacked conformation at the interface. This highlights the importance of the p + -p + interaction of imidazolium cations through the stacked conformation at the RTIL/air interface, despite the cations being electrostatically repulsive.

Layering structure of RTILs at the RTIL/graphene interface
We have discussed the molecular layering of RTILs at the RTIL/air interface. To demonstrate that the B10 Å molecular layering distance is intrinsic to the nature of RTILs, we also investigated the RTIL/graphene interface. The number density profiles at the RTIL/graphene interface are displayed in Fig. 5(a). The amplitude of the (r + À r À ) peak area is much larger at the RTIL/graphene interface than that at the RTIL/air interface, indicating that the graphene sheet has strong templating effects on the RTILs. This is consistent with the previous NP-FFMD simulation. [64][65][66][67][68] The peak area of (r + À r À ), DS, at the RTIL/graphene interface is plotted in Fig. 5(b) together with DS at the RTIL/air interface. The normalized fit curves of Df (z) are also displayed in Fig. 5(b), while the fit parameters are summarized in Table 1. These plots show very similar decay lengths for both the RTIL/air and RTIL/ graphene interfaces, illustrating that the range of the molecular layering is not affected by the RTIL-substrate interactions but is governed by the nature of the RTIL.

p + -p + conformations with PBE XC functional
Above, we showed the comparison of interfacial RTIL structures predicted by the AIMD and FFMD simulations and found that the p + -p + stacked conformation is essential for the layering structure of RTILs at the RTIL interfaces. In this section, by comparing the RTIL layering structures predicted using the BLYP + D3 and PBE + D3 methods, we again demonstrate that the rich p + -p + stacked conformation at the RTIL/air interface is very sensitive to the p + -p + interaction energy. Here, a key property is that the PBE + D3 level of theory significantly lowers the binding energy of the stacked conformation relative to the BLYP + D3 level of theory. 51 Therefore, one can expect that the PBE + D3-AIMD simulation provides much less stacked conformations compared with the BLYP + D3-AIMD simulation. The axial profiles of the stacked conformations calculated using the PBE + D3-AIMD simulation are displayed in Fig. 6. Indeed, the comparison of the PBE + D3 and BLYP + D3 levels of theory in Fig. 6 shows that more stacked conformations are observed in the BLYP + D3-AIMD simulation. This illustrates that, as expected, the p + -p + stacked conformations are very sensitive to the p + -p + interaction energy.

Consistency between simulation and experiment
Finally, we discuss the consistency of the current simulation and experimental data. We show that the p + -p + interaction of imidazolium cations stabilizes the RTIL layering structure, against both electrostatic repulsion and the molecular polarizability, which both serve to reduce the molecular layering near the RTIL interface. As a result, the AIMD simulation coincidentally shows a very similar length scale of molecular layering as the NP-FFMD simulation. A good agreement between the NP-FFMD and experimental data has been reported; both X-ray scattering and NP-FFMD show 10 Å length of molecular layering Fig. 5 (a) Axial profiles of the RTIL number densities at the RTIL/graphene interface in the AIMD simulation. The number density of graphene (shown in black) was divided by 50. The origin point (z = 0 Å) was set to the position of the first peak position of r + . (b) Normalized DS(z), of the RTIL/graphene and RTIL/air interfaces in the AIMD simulation. The fitting function is given by Df (z) = f 0 (1 À exp(À(z À c)/l))exp(Àz/x), and the fit parameters are given in Table 1. Table 1 The fit parameters of Df (z) = f 0 (1 À exp(À(z À c)/l))exp(Àz/x) for the AIMD (BLYP + D3), NP-FFMD, and P-FFMD simulations at the RTIL/air interface ( Fig. 3(a-d)) and the AIMD simulation at the RTIL/graphene interface (Fig. 5(a)). l was obtained from global fits to all the data, yielding l = 4.0 Å. The standard deviations of S 0 , c, and x for the RTIL/air interface are at most 5%, while those for the RTIL/graphene interface are at most 15%  interfaces. The NP-FFMD simulation gets the correct length scale for the wrong reasons: our finding indicates that a coincidental cancellation of the lack of the p + -p + interaction and molecular polarizability results in a good agreement of the NP-FFMD simulation and the experiment. Furthermore, a similar length scale of RTIL layering at the mica 7 and sapphire 8 interfaces also suggests that the molecular layering structure is less influenced by the substrate, consistent with our findings.

Conclusion
We studied the length scale of the molecular layering for imidazolium RTILs at the RTIL/air and RTIL/graphene interfaces, by using the AIMD simulation. Comparison of the AIMD, NP-FFMD, and P-FFMD simulation data indicates the competing effects of the randomization of RTILs due to the molecular polarizability and the enhanced molecular layering due to the p + -p + interactions of the imidazolium rings. A very similar length scale of RTIL layering at the RTIL/air and RTIL/graphene interfaces shows that the length scale of the molecular layering structure is dictated by the RTIL-RTIL interactions and is not affected by the substrates. Our study highlights the importance of p + -p + interactions of the cation imidazolium rings in stabilizing the surface structure of the imidazolium RTILs, which is essential, for example, for the charge transport properties in the RTILs and between the RTILs and the electrodes.