## 2 Simulation System and Methods  
The molecular-level structures of [EMIM]\({}^{+}\)/[BF\({}_{4}\)]\({}^{-}\) ion pairs, vertically oriented graphene and the graphene/ionic liquid interface are illustrated in Figure 1, where the atoms on the ions are labeled by the numerical and element symbols in Figure 1a. The non-rough vertically oriented graphene electrode is modeled by 12 successive graphene nano-ribbon layers which are stacked in A-A sequence [30]. All the graphene layers orient toward the ions with their edge planes, forming the electrode surface. According to previous studies [30], as shown in Figure 1c, the layer space of the multi-graphene layers in a non-rough electrode is set as 0.34 nm. Compared with the multi-graphene layers in the non-rough electrode, the rough electrode has the same layer space and different rough surface, possessing grooves with a width of 1.01 nm and depth of 0.85 nm, as depicted in Figure 1b. The non-rough and rough electrodes both have lateral sizes of 4.26 \(\times\) 3.94 nm\({}^{2}\). In the electrolyte near the electrode surface, 300 [EMIM]\({}^{+}\)/[BF\({}_{4}\)]\({}^{-}\) ion pairs are randomly packed with a thickness of about 4.0 nm. Outside the electrolyte region, another 4.0 nm vacuum space is added to eliminate the interaction between the particles and their images outside of the unit cell in the \(z\)-direction.  
Figure 1: The detailed illustration of the simulation system: (**a**) ion structure of [EMIM]\({}^{+}\) and [BF\({}_{4}\)]\({}^{-}\), (**b**) stack structure of graphene electrode, (**c**) interface models of graphene/ion liquid with non-rough and rough graphene surface.  
The simulations are conducted in the canonical ensemble (NVT) using a universal force field (UFF) [32]. For each [EMIM]\({}^{+}\) cation and [BF\({}_{4}\)]\({}^{-}\) anion, the partial charge \(q\) on each atom are modified and assigned according to the previous studies [33] (as shown in Table 1) and remain unchanged during the simulations. Periodic boundary conditions are imposed in all three directions. Under these boundary conditions, particles leave the domain at one side and reenter the domain at the opposite side with the same velocity. During the simulation, the electrode atoms are fixed in their initial positions in all cases. A Nose thermostat [34] is used to maintain the system temperature at 300 K. Starting from a random distribution of [EMIM]\({}^{+}\)/[BF\({}_{4}\)]\({}^{-}\) ion pairs near the electrode surfaces, the systems above are firstly relaxed to eliminate the undesirable interactions between ions and electrode. In addition, the SMART method [35] is used in this step to minimize the system energies. Then, the equilibrium and production processes are continued for 2 \(\times\) 10\({}^{6}\) steps with a time step of 1 fs, and the atom positions and micro-structure of the electrode/electrolyte interface are saved every 20 ps. The Ewald summation method [36] is used for Coulombic interactions with an accuracy of 10\({}^{-3}\) kcal/mol, and the atom-based summation method [37] is used for van der Waals interactions with a cutoff of 1.2 nm. To simulate the charge/discharge process, different constant charge densities \(\sigma\) (the ratio between total charge and atom numbers at the electrode surface) are set on the electrode surface as \(\sigma\) = 0, \(\pm\)0.02, \(\pm\)0.05 and \(\pm\)0.08 e/atom. According to the results of the ion-concentration distributions normal to the graphene surface under different charge densities, the differential capacitances are finally calculated by solving the Poisson equation.