Graphene-ionic liquid interfacial potential drop from DFT-MD simulations
Heigo Ers Meeri Lembinen Maksim Mišin Ari P. Seitsonen Maxim V. Fedorov Vladislav B. Ivaništšev
cond-mat.mtrl-sci












§ ABSTRACT


Ionic liquids (IL) are promising electrolytes for electrochemical applications due to their remarkable stability and high charge density.
Molecular dynamics simulations are essential for better understanding the complex phenomena occurring at the electrode–IL interface. 
In this work, we have studied the interface between graphene and 1-ethyl-3-methyl-imidazolium tetrafluoroborate IL, using density functional theory-based molecular dynamics simulations at variable surface charge densities. 
We have disassembled the electrical double layer potential drop into two main components: one involving atomic charges and the other – dipoles.
The latter component arises due to the electronic polarisation of the surface and is related to concepts hotly debated in the literature, such as the Thomas–Fermi screening length, effective surface charge plane, and quantum capacitance.







§ INTRODUCTION


The ever-growing, world energy consumption has created a demand for higher capacity energy storage devices.<cit.> 
Supercapacitors were developed to meet that demand by storing energy in the electrical double layer (EDL) – at the interface between an electrode and an electrolyte.<cit.>
Processability and a low cost are essential characteristics when choosing electrode material for supercapacitors. 
That is why one of the most widely used electrode material is graphitic carbon.<cit.> 
To allow for the accumulation of more charge, the specific surface area of graphitic carbon is increased via chemical or thermal activation,<cit.> resulting in capacitance as high as 120, 300, and 245 F/g in organic, aqueous, and ionic liquid (IL) electrolytes, respectively.<cit.> 
In addition to activated carbon,  the potential usage of nanostructured carbon and graphene (Gr) are gaining broad interest.<cit.>
However, all carbon-based materials have a limitation due to their quantum properties. 
Upon an increase of a specific surface area, the electrode thickness reduces to the size of a space-charge region. 
That, in turn, reduces the number of mobile charge carriers, which decreases the energy storage capacity.<cit.> 
In general, the capacitance of a Gr electrode–electrolyte interface is limited by the capacitance of Gr.<cit.> 
In particular, the electronic structure of Gr directly relates to the characteristic capacitance minimum near the potential of zero charge.<cit.>
Still, by choosing a suitable electrolyte, one can enhance the energy storage capacity. 
For example, by widening the electrochemical stability window with the use of an IL electrolyte, the storage capacity of the supercapacitors can be improved. 
That is why thermally stable and moisture resistant room temperature ILs are promising electrolytes for electrochemical application.<cit.>

Because macroscopic quantities such as capacitance arise from the microscopical process taking place in the EDL, an understanding of the structure and dynamics of interfaces is crucial for the development of supercapacitors. 
The properties of ILs near surfaces have been actively studied since 2007.<cit.>
Numerous computational studies have focused on the interfacial structure, dynamics, and properties of ILs, using classical molecular dynamics (CMD) <cit.>, density functional theory (DFT) <cit.>, and DFT-based molecular dynamics (DFT-MD) <cit.>. 

For the electrode model, different carbon-based materials were examined; for the electrolyte model, various imidazolium-based ILs were investigated. Specific attention was given to the interface between a Gr sheet and 1-ethyl-3-methylimidazolium tetrafluoroborate (EMImBF_4) because of the high specific surface area of Gr coupled with the large electrochemical stability window of EMImBF_4.<cit.> 
Previously, several capacitor models consisting of carbon-based electrode and EMImBF_4 were studied with CMD by Shim et al.<cit.>, Merlet et al.,<cit.> and Paek et al.<cit.>
Shim and co-workers noticed that interfacial capacitance is higher for the positively charged electrode due to the size differences between the ions.
Merlet et al. reproduced the results of the earlier atomistic simulation using a coarse-grained model of EMImBF_4.<cit.>
Paek et al. conducted, in addition to CMD, a 1.5 ps-long DFT-MD simulation using a model made up of 15 EMImBF_4 ion pairs and a Gr sheet of 60 atoms. By comparing the potential profiles, Paek et al. concluded that the polarisation effects should be included when evaluating the interfacial capacitance.<cit.>

In this study, we have examined the Gr–EMImBF_4 interface with large-scale DFT-MD simulations focusing on its microscopic (electrode charge screening and polarisation) and macroscopic (surface charge density and potential drop) characteristics. 
We show how the analysis of the obtained electronic and geometric structures leads to the estimation of the central quantity – the potential drop across the EDL.







§ COMPUTATIONAL METHODS




 §.§ Classical molecular dynamics simulations


To run CMD simulations of the Gr–EMImBF_4 interface, we used a model with two rigid, parallel Gr sheets with an area of 3.408×3.4433 nm^2.
The 10.472 nm space between the sheets was filled with 450 EMImBF_4 ion pairs using Packmol<cit.>.
All the simulations were run using Gromacs  5.1.4<cit.> software and NaRIBaS scripting framework<cit.>. 
OPLS-AA force field was used with an effective dielectric constant of 1.6.<cit.> 
Each system was pre-equilibrated for 0.1 ns, annealed at 1000, 900, 800 K for 3 ns to produce three replicas, and simulated in the NVT ensemble for 10 ns at 450 K, which was controlled by the v-rescale thermostat <cit.>. 
In between the pre-equilibration and annealing, the electrodes were equally but oppositely charged within 2 ns.
The surface charge densities (σ) of 0 and ±0.5 e/nm^2 were set as point charges of equal magnitude on the surface atoms. 
All other parameters were set the same as in Ref. .



 §.§ DFT-based molecular dynamics simulations


All model systems used in the DFT-MD simulations were cut out from the last step of the CMD trajectory.
Figure <ref> illustrates the DFT-MD simulation cell, where one side of the cell is the Gr sheet (448 carbon atoms) in contact with IL, which consists of 400 ions. To have a similar σ as in the CMD simulations, an excess of EMIm^+ cations or BF_4^- anions was introduced to cause a natural charge redistribution between IL and Gr. Periodic boundary conditions were applied in x,y-directions, parallel to the surface. A vacuum layer was added to extend the z-axis along the surface normal of Gr from 5 to 8 nm. In this manner, two replicas were created for two systems with the charged electrode and three replicas for the system with the neutral electrode. 
The exact number of ions in each system are given in Figure <ref>. 
The DFT-MD simulations at 353 K, in the NVT ensemble, were performed with CP2k software package versions 2.6 and 5.0 <cit.>, implementing the Gaussian plane wave method <cit.>. 
The computations were carried out using the energy cut-off of 600 Ry for the plane wave grid, Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional <cit.>, optimised double-zeta basis set, with corresponding Goedecker–Teter–Hutter (GTH) pseudopotentials,<cit.> and the D3 dispersion correction by Grimme et al. <cit.>.
It was previously shown in Refs.  that the PBE+D3 combination gives accurate interaction energies for IL ion pairs. 
The duration of the simulations was from 2.3 to 8.1 ps with a time step of 0.5 fs. Within the 0.5–1 ps the bond lengths and angles relaxed from the values given by the CMD force field to the ones consistent with the PBE+D3/GTH/DZVP-MOLOPT theory level.
For further analysis 32–37 random snapshots were chosen after 2 ps of simulation for charge analysis with Chargemol 3.5 software implementing the sixth version of density derived electrostatic and chemical (DDEC6) partitioning approach that shows accurate reproduction of the electrostatic potential for a variety of periodic and non-periodic systems <cit.>. 
The averages over these snapshots were used for the system characterisation. 





 §.§ Density of States calculation

Unlike in the CMD simulations, DFT calculations explicitly account for the electronic density distribution, which can be used for the calculation of the projected density of states (pDOS). Based on the pDOS of Gr, one can estimate σ as

    σ_DOS=n(e^-)/S-∫_-∞^∞D(E)P(E)dE,

where D(E) is the projected density of states for Gr at the energy level E, P(E) the occupancy of the level, n(e^-) the number of carbon atoms' valence electrons in case of neutral Gr in vacuum, and S the surface area of Gr in contact with IL.
The additional potential resulting from electronic density redistribution in Gr caused by charging can be estimated from the DOS,<cit.>
 
 
    Δϕ_DOS=E^0_F-E^σ_F/e,


where E^σ_F is the Fermi energy of the Gr sheet with a surface charge density σ, E^0_F the energy at the minimum in the DOS shown in Figure <ref>, and e the elementary charge.

Using the pDOS of Gr, it is also possible to evaluate the Thomas–Fermi screening length of Gr using the equation <cit.>:


    l_TF=ϵ_0ħ/4e^2√(2E_F/π m_e n),


where the concentration of mobile charge carriers (n) in Gr can be estimated as


    n=∫_-∞^E_FD(E)[1-f(E-E_F,T)]   dE  + ∫_E_F^∞D(E)f(E-E_F,T)]   dE.


In Eq. <ref> f(E-E_F,T) stands for the Fermi–Dirac distribution function [exp(E-E_F/k_BT)+1]^-1.




 §.§ Overscreening and potential drop


To quantify the overscreening phenomenon, we used the overscreening factor (β):

    β(z)=-1/σ∫_0^zρ(z')dz',

where ρ_ion(z) is the ionic charge density at distance z from the electrode, and σ is the surface charge density.
The surface charge is screened when β=1 and overscreened when β>1.

To calculate the potential drop, we used the charge density obtained by the means of DDEC6 analysis and Poisson's equation. 

When studying the potential drop in a direction perpendicular to the Gr surface, integration of Poisson's equation yields the expression for electrostatic potential at a point z from the electrode:

    ϕ(z)=-1/ϵ_0∫_0^z (z-z^')ρ(z^')dz^'


The potential drop (ϕ_ion) between the electrode and the electrolyte was calculated as


    Δϕ_ion=ϕ_Gr-ϕ_electrolyte.


where ϕ_Gr is electrostatic potential at the geometric centre of the Gr sheet (z=0) and ϕ_electrolyte is an average of electrostatic potential in IL 2.5–4.0 nm from the electrode. 

When estimating the potential drops, instead of mean average, we chose the weighted average approach to reduce the influence of noise on the resulting values. 
First, we calculated the potential drop for each replica by finding the mean average of the potential drops over its snapshots. 
Then, taking the inverse square of the standard deviation of the potential drops of a replica as a weighting factor, we estimated the potential drop of a system with σ as the weighted average over the replicas.







§ RESULTS AND DISCUSSION




 §.§ Overscreening phenomena


One of the main motivations of the presented study was to provide a DFT-level examination of the overscreening phenomena.
While overscreening naturally appears in the CMD, previously it was not clear how it is affected by the polarisation inherent at the DFT level. 
For instance, calculations of Valencia et al. indicated a possibility of charge transfer between the contact ions and the surface,<cit.> which in principle should diminish the overscreening.

DDEC6 analysis revealed no charge transfer between the electrode surface and the ions in the first layer – the charges of anions and cations of -0.82e and +0.82e, respectively, were the same for ions in contact with the surface and 2.5–4 nm away from it.
Also, within the length of the obtained trajectories, we did not notice any statistically significant changes in the potential profiles (shown in Figure <ref>a) that would indicate the disappearance of the alternate layering.

As an example of density fluctuations, Figure <ref>d illustrates how the layering arises due to geometric constraint on ions caused by the neutral Gr surface.<cit.> The first layer of the IL, which is considered to start at the van der Waals radius for carbon (0.17 nm from the surface), is 0.40 nm broad and the second layer is 0.42 nm wide. Herewith, near the neutral electrode, the layers contain almost the same number of anions and cations so that the densities of the layers equal the bulk density of 1.24 g/cm^3.<cit.>

The densities of the layers become different from the bulk one when the electrode charging causes the overscreening and segregation of anions and cations into distinct layers. Changes in the overscreening parameter β, shown in Figure <ref>b and e, reflect the alternation of ionic layers.
As can be seen in Figures <ref>c and f, β values are larger than one for both negative and positive electrodes implying dense packing of counter-ions and overscreening.
Due to the smaller size of BF_4^-, a denser (2.74 ions/nm^2) and more ordered structure appears at the positive electrode than at the negative one (1.75 ions/nm^2).
Herewith, in all cases the non-uniform ion distribution correlates with the non-uniform surface charge distribution on the Gr electrode as indicated in Figures <ref>b, <ref>c and <ref>f.

Overall, Figures <ref>a–f illustrate that the overscreening phenomena is not an artifact of the CMD simulations as it persists in the more realistic DFT-MD simulations with implicitly included polarisation.



 §.§ Interfacial potential drop


Figure <ref>a shows the weighted average potential (ϕ(z)) profiles in a perpendicular direction to the Gr surface.
As follows from Figure <ref>a, the potential drop occurs mostly within 0.3 nm of the Gr surface. 
Then the potential continuously changes until 2.5 nm with pronounced fluctuations. 
At larger distances, the fluctuations dampen and the potential reaches a plateau. 



Table <ref>  presents values for two corrections to the potential drop that can be evaluated from the DFT data. 
The difference between levels defined in Eq. <ref> and shown in Figure <ref>b describes the additional work done to bring an electron into Gr.
In addition to DOS, the potential change due to electronic density redistribution can also be calculated from the DDEC6 analysis providing dipoles and quadruples in addition to point charges.
The dipoles characterise well the uneven distribution of charge shown in Figure <ref>a. 
Thus, at first approximation, to incorporate the electrode polarisation, the correction term (Δϕ_dip) can be calculated from the electrode dipole moment in z-direction and added to Eq. <ref>:

 
    Δϕ_dip=p_z/ϵ_0 S,


where p_z is the average of dipole moments of the Gr atoms in z-direction and S is area of the electrode. 
The comparison between the correction terms in Table <ref> reveals a similarity in their values caused by the interrelation between the filling of electron levels and the dipole moment formation. 
Small differences between the correction terms arise due to the constraints of the DDEC6 and DOS analyses.
Furthermore, the potential correction terms can be recalculated into the position of the effective surface charge plane (l) along the z-axis relative to the center of the electrode: 
 
    l=ϵ_0Δϕ/σ.

Theoretically, the surface charge plane should be located one-half of an interplanar spacing. Due low concentration of mobile charge carriers, in Gr the plane is shifted by the Thomas–Fermi screening length towards the surface <cit.>. Comparison of the theoretical a/2-l_TF and obtained l_dip/DOS reveals significant difference, yet shows similar σ-dependence.


To the best of our knowledge, the dipole correction term is defined for the first time in the context of the electrode–IL interface and computer simulations. 
Earlier Kornyshev and Schmickler suggested estimating the effect of electric field on the electrode by evaluating the effective charge density plane from electronic density distribution at the DFT level of theory.<cit.> 
A similar approach was used to correct the MD simulations results by Ruzanov et al.<cit.> 
Using the DOS, Paek et al. accounted for the potential drop within the electrode to express the interfacial capacitance as 1/C_total = 1/C_Gr + 1/C_IL.<cit.>
Dufils et al. also used the DOS to evaluate Thomas–Fermi screening length, which they used to tune the metallicity of electrodes in MD simulations <cit.>.
The suggested dipole correction provides an alternative way to evaluate the potential drop within the electrode.
Similarly to the mentioned works , the dipole correction can be expressed as an average quantity for a surface as a whole.
In that case, a separate calculation is required to obtain a value of dipole correction for each field, potential, or surface charge value.
Alternatively, a correlation between the dipole moment and charge on individual surface atoms (shown in Figure <ref>c) can be used.
Atomic charges can be already self-consistently calculated using the constant potential method, for example developed or used in Refs. <cit.>.
Applying the dipole correction over the atomic charges opens a door for the low-cost augmentation of the CMD simulations to account for the electrode polarisation.









§ CONCLUSIONS

In this work, we used DFT-based MD to explore the Gr–1-ethyl-3-methyl-imidazolium tetrafluoroborate interface to examine the charge distribution and polarisation on the Gr surface accounting for the dynamics of ions near the surface. The results indicate that overscreening is preserved upon the inclusion of quantum effects. No indication of the charge transfer between the IL and the electrode was seen and the layered structure in the IL near the Gr surface remained stable during up-to 8 ps long DFT-based MD simulation.

The analysis showed that the atomic charges on the electrode are distributed non-uniformly due to electronic lateral polarisation of the electrode. Polarisation and redistribution of Gr electron density in the direction perpendicular to the surface induces a surface dipole that significantly affects the interfacial potential drop.

To reconcile the differences between classical and DFT-based MD simulations, we introduced a correction term, which quantifies the redistribution of the Gr electron density using atomic dipoles that correlate with the atomic charges. That correction term accounts for the effect of polarisation on the interfacial potential drop and can be used in classical molecular dynamics simulations.



§ ACKNOWLEDGEMENTS

This research was supported by the Estonian Research Council (grants IUT20-13, PSG249, and PUT1107), the EU through the European Regional Development Fund (TK141 “Advanced materials and high-technology devices for energy recuperation systems”), and the Estonian–French cooperation program Parrot, funded by the Estonian Research Council and Campus France. For providing us with the computational resources, we would like to acknowledge the Partnership for Advanced Computing in Europe (PRACE), the Distributed European Computing Initiative (DECI), the HPC Center of the University of Tartu, and Skoltech Pardus HPC cluster.













