Comparison of the Performances of Different Computational Methods to Calculate the Electrochemical Stability of Selected Ionic Liquids

The electrochemical stability windows (ESW) of selected ionic liquids have been calculated by comparing different computational approaches previously suggested in the literature. The molecular systems under study are based on di-alkyl imidazolium and tetra-alkyl ammonium cations coupled with two different imide anions (namely, bis-fluorosulfonyl imide and bis-trifluoromethyl sulfonyl imide), for which an experimental investigation of the ESW is available. Thermodynamic oxidation and reduction potentials have here been estimated by different models based on calculations either on single ions or on ionic couples. Various Density Functional Theory (DFT) functionals (MP2, B3LYP, B3LYP including a polarizable medium and empirical dispersion forces) were exploited. Both vertical and adiabatic transitions between the starting states and the oxidized or reduced states were considered. The approach based on calculations on ionic couples is not able to reproduce the experimental data, whatever the used DFT functional. The best quantitative agreement is obtained by calculations on single ions when the MP2 functional in vacuum is considered and the transitions between differently charged states are vertical (purely electronic without the relaxation of the structure). The B3LYP functional underestimates the ESW. The inclusion of a polar medium excessively widens the ESW, while a large shrinkage of the ESW is obtained by adopting an adiabatic transition scheme instead of a vertical transition one.


Introduction
Electrochemical devices, such as lithium batteries, fuel cells and solar cells, require suitable electrolytes able to allow a facile migration of ions while being inert at the electrodes. Electrolytes must fulfill many requirements: high ionic conductivity, high thermal and electrochemical stability and wide liquid ranges, among the others. Under these aspects, many ionic liquids (ILs) are apparently suitable to compete and replace standard organic solvents, thanks to their specific physico-chemical properties [1,2]. With the development of high voltage lithium batteries, the need for innovative electrolytes able to sustain highly oxidizing electric potentials is becoming more apparent and poses new challenges.
Experimentally, many ionic liquids have been investigated to establish their electrochemical stability window [3,4]. Overall, the anions of the ionic liquids play a central role, as they inevitably determine the highest voltage at which the devices can operate without the degradation of the electrolyte (anodic limit) [4]. However, from an experimental point of view, the determination of the electrochemical stability of ionic liquids poses some challenges and the interpretation of data is not always straightforward. Indeed, it is well known that the chemical nature of the electrodes (e.g., Pt, Ni, stainless steel, aluminum, carbon composites, etc.) can increase or decrease the electrochemical stability window; moreover, the acceptable current level for the definition of non-degrading electrolyte varies from study to study [4], as well as between data analysis algorithms. Finally, different reference electrodes are often used for different applications, thus shifting the anodic and cathodic ranges. In general, however, it has been shown that the ionic liquids with the highest anodic stability are highly fluorinated.
The electron affinity and ionization of anions and cations composing ionic liquids have also been investigated computationally by many authors, and some correlation between computational results and the electrochemical stability window (ESW) has been proposed. A simple approach suggests a linear correlation between the energy of the highest unoccupied molecular orbital (HOMO) of anions and their oxidation potential [5,6]. On the other hand, the computation of the energy difference for a vertical transition from the charged ion and the neutral state is easily performed, either in the gas phase or in the presence of a suitable polarizable medium (PM). This second approach allows us to mimic ionic liquids' oxidation and reduction in liquid media [7,8]. Overall, the inclusion of the PM strongly overestimated the electrochemical stability while, in the gas phase, DFT calculations with the B3LYP or VSXC functional provided results in quantitative agreement with experimental data [7]. Computational studies have also been carried out on cations, and a detailed investigation of the changes of their reduction potential in respect to substitutions in their alkyl chain has been reported, exploiting hybrid DFT with the "classic" B3LYP functional [9]. Furthermore, the linear correlation between electron affinities estimated through vertical and adiabatic transitions between charged and neutral states hasbeen sketched [9]. A combination of DFT calculations and molecular dynamics was used by Ong et al. [10] to estimate the ESW of neutral ionic pairs between common cations (1-butyl-3-methylimidazolium (BMIM) and N,N-propylmethylpyrrolidinium (P 13 )) and simple anions, PF 6 , BF 4 , and bis(trifluoromethylsulfonyl)imide (TFSI). Additionally, ab-initio molecular dynamics have been proposed for the investigation of ESW [11]. Various DFT functionals have been applied to establish a correlation between HOMO-LUMO energy differences, ionization potentials and electron affinities for both adiabatic or vertical transitions of single ions or ion pairs in various chemical systems [12][13][14]. Cheng [15] and Kazemiabnavi [16] exploited full thermodynamic cycles to calculate the reduction and oxidation potentials of ILs, following a strategy suggested by Winget et al. and Marenich et al. for simple systems [17,18]. Cheng and Kemiabnavi modelled the thermodynamics of oxidation and reduction processes, also considering ion salvation [15,16]. Similar approaches have also been used for a high throughput screening of hundreds of possible anions and cations obtained by substitutions with defined functional groups in common anions and cations [15,19], using hybrid-DFT calculations and the B3LYP functional.
We recently reported an experimental physico-chemical investigation of novel liquids based on imidazolium or tetra-alkyl-ammonium cations, coupled with bis(perfluroalkylsulfonyl)imide anions (i.e., EMIFSI, EMITFSI, N1114FSI, N1114TFSI and N122(2O1)TFSI, respectively). We demonstrated their use as solvents in electrolytes for lithium-ion batteries working up to 5 V [20,21]. The anodic stability of these ionic liquids was measured on carbon composite electrodes over aluminum using a 0.8:0.2 molar mixture with the LiTFSI salt [20] and, in all cases, they overcame 4.5 V vs. Li. Here, we tackle the challenge of a systematic computational investigation of the thermodynamic stability widows of these ionic liquids by comparing the experimental data to various computational approaches proposed in the literature. Our aim is to evaluate the reliable computational approximation for a precise estimate of the ESWs of ionic liquids. We carried out this study on the five different ionic liquids studied by us previously [20]; the structure of the ILs ionic constituents is shown in Figure 1.

Materials and Methods
All computational investigations of the single ions or of the ionic couples composing the ionic liquids were performed by means of Spartan software [22]. For each ion or ionic couple, possible conformers were investigated by means of a systematic rotation of flexible bonds, exploiting built-in-routines. Unrelaxed and relaxed structures were computed at MP2 [23] or DFT levels of theory [24]: 0K electronic energy as well as Gibbs energy at finite temperatures considering thermal effects were calculated. Various computational methods were used: MP2 [25,26], B3LYP [27], B3LYP with a polar solvent (dimethylformamide, ε r = 37.22) using the C-PCM algorithm [28] and B3LYP with the same polar solvent and dispersion forces. For all calculations, a 6-31G** basis set was employed.
Oxidation and reduction reactions were mimicked using different approaches, starting from a structurally relaxed minimum energy molecule/ion (initial state): 1. Vertical transition of single ions, where the final state preserves the structure of the initial state. 2. Adiabatic transition of single ions, where the final state is energetically relaxed to its minimum. 3. Vertical or adiabatic transition of neutral ionic couples.
In the case of vertical transitions, ESW were estimated, starting from oxidation/reduction reaction energies calculated by subtracting the computed total energy, i.e., E tot , of the final state (product) to that of the initial one (reagent). For an oxidation reaction where the reagent (initial state) is the chemical specie A with net charge n:

4.
For a reduction reaction where the reagent (initial state) is the chemical specie B with net charge m: For anions, i.e., A − , the oxidation reactions involve the transition from the anion to a neutral state, while the reduction reaction is due to the transition from the anion to a dianion state. The first reactions give the anodic limit, while the second is responsible for the cathodic one.
For cations, the oxidation transforms the cation into a di-cation, and it is responsible for anodic stability, while the reduction transforms the cation into a neutral species and is linked to the cathodic limit.
In the case of adiabatic transitions, either starting from ions or ionic couples, the energy differences can be calculated by considering explicitly zero-point energies (ZPE) for both reagents and products, both of which have vibration real frequency values. Thus, one may obtain the ∆G o 0K,anodic and the ∆G o 0K,cathodic starting from the Gibbs energies at 0K for the initial and final states (namely G tot,0K = E tot + ZPE).
The representation of ESW using the usual relative electrochemical potentials has been obtained in the case of vertical transitions by the following relations: where F is the Faraday constant and ∆E is the oxidation or reduction electronic reaction energy, and the term−1.46 V is due to the necessity of referring these limits to the standard Li + /Li 0 electrode, as discussed in ref. [12]. In the case of adiabatic transitions, the ∆E is substituted by ∆G o 0K .

Vertical Transitions Approximation Starting from Ions
Tables 1-4 report the electronic energy of the ions and neutral species for the two conformers of the FSI and TFSI anions, and for the lowest energy conformer of the cations; the corresponding anodic and cathodic limits, calculated according to Equations (7) and (8), are also reported. Each table comprises a specific computational method: 3. B3LYP in a simulated polar medium 4. B3LYP in a simulated polar medium with the inclusion of empirical dispersion forces.
In all cases, the structure of the initial anion or cation was optimized at the corresponding level of theory.
In all cases, as expected, the anodic and cathodic limit of the anion is lower than that of the cation; therefore, the limits of stability of whole ionic liquid are the anodic limit of the anion and the cathodic limit of the cation.
The MP2 method (see Table 1) gives anodic limits of the conformers for the two bis(perfluroalkylsulfonyl)imide anions values between 4.8 and 5.0 V and estimates in the range of −1.5 to 0.5 V for the cathodic limit of the cations. These values are in good agreement with our experimental determinations [20]. Table 1. Electronic energy (in atomic units) of anions and cations and their unrelaxed neutral or double charged forms calculated at the MP2/6-31G** level of theory and derived anodic and cathodic limit (V).  Table 3. Electronic energy (in atomic units) of anions and cations and their unrelaxed neutral or double charged forms calculated at the B3LYP/6-31G** level of theory in a polar medium and derived anodic and cathodic limit (V).  The electrochemical stability window shrinks slightly while using the B3LYP functional: the anodic limit decreases to 3.9/4.5 V (Table 2) while the cathodic limit increases around −0.6/+1.2 V. These values slightly underestimate the oxidation limit while overestimating the reduction stability. The introduction of a polar medium widens the electrochemical stability window (−2.4/+6.4 V) (see Table 3), overcoming the experimental determinations, while the addition of empirical dispersion forces has a marginal effect (see Table 4).

Adiabatic Transition Approximation Starting from Ions
Tables 5-7 report the values obtained at the levels of MP2, B3LYP in vacuum and B3LYP in a polar medium. Due to the increased computational time needed to calculate the vibrational properties of ions, we restricted calculations to the anodic limit of the anions and the cathodic limit of the cation; these are the relevant values to derive the electrochemical stability of the ionic liquids. Compared to the figures obtained for vertical transitions, the adiabatic transitions values are lower for the anodic stability limit and higher for the cathodic limit. This effect is expected since the structural relaxation of reduction/oxidation products leads to a stabilization of the final state and thus to smaller reaction energies. Overall, the EWSs range between 2.1/2.3 V and 3.5/4.1 V for the MP2 model (see Table 5), 2.2/2.9 V and 3.2/3.6 V for the B3LYP functional (see Table 6) and between−0.1/+1.1 V and 5.1/5.6 V for B3LYP in a polar medium (see Table 7). These values are clearly not consistent with the experimental ones [20]. Moreover, in the case of relaxed neutral EMI at the MP2 level and of N122(2O1) at the B3LYP level, no convergence of the structure could be obtained, even after more than 1000 iterations, due to the structural instability of both doublet final states. One should also note that after the adiabatic transitions, starting from ammonium ions, one of the alkyl chains tended to detach from the N atom, suggesting a tendency for cleavage by the molecule. It must be noted that the variation of ZPE between the initial and final states of the transition is small (<30 kJ/mol) in comparison to the electronic energy difference; therefore, ZPE practically cancels out upon subtraction.
To help to visualize the performances of various functionals in describing the ESW of TFSI and FSI anion, Figure 2 reports a graphical representation of the anodic limits calculated at different theory levels compared to the experimental values. Tables 8 and 9 report the electronic energy of the neutral, positive and negative ionic couple, as well as the derived anodic and cathodic limits calculated at the B3LYP level of theory in a vacuum (Table 8) or in a polar medium (Table 9), adopting either a vertical or an adiabatic approximation. Calculations based on the B3LYP functional in a vacuum clearly give an extremely high anodic limit and an extremely low cathodic limit, either considering or not considering any structural relaxation. Additionally, values obtained by the B3LYP model in a polar medium from vertical transitions display the same behavior. The estimated B3LYP in a polar medium considering the structural relaxation of the ionic couples (adiabatic approximation) is apparently closer to that in the experimental determination. To confirm the reliability of this approach (ionic couples in polar media under adiabatic transition approximation) beyond the five ionic couples, we calculated the anodic and cathodic limits of a quite different ionic liquid: Butyl methyl imidazolium chloride (BMIMCl). The ESW of BMIMCl was experimentally investigated in Ref. [29] using an Ag/Ag+ reference electrode. Taking into account the shift between the Li and Ag reference electrode, BMIMCl is expected to have anodic and cathodic limits at 3.38 V and 1.12 V, respectively, which are far from our predictions for ionic couples. Therefore, even considering the partial error compensation provided by the simultaneous adiabatic approximation and the adoption of a continuous polar medium, it is unlikely that any of the methods based on calculations on ionic couples can provide reliable values of anodic or cathodic limits of ionic liquids.  Table 8. Electronic energy (in atomic units) of the neutral ionic couples and of the ionic couple with the addition or removal of an electron, without or with the possible relaxation of the geometry and derived cathodic and anodic limits (in V). All calculations were performed at the B3LYP/6-31G** level of theory in gas phase.  Table 9. Electronic energy (in atomic units) of the neutral ionic couples and of the ionic couple with the addition or removal of an electron, without or with the possible relaxation of the geometry and derived cathodic and anodic limits (in V). All calculations were performed at the B3LYP/6-31G** level of theory in a polar medium.

Discussion
From the previously reported data, it is quite evident that the adoption of a specific approximation, either in terms of computational conditions or thermodynamic approximation, has a remarkable impact on the accuracy and precision of ESW prediction.
The accuracy of the adopted methodology can be evaluated for the anodic limits (oxidation potentials) for all five ionic liquids by comparing the computational data with the experimental determination reported by us previously [20]. Two statistic quantities can be evaluated:

MAE: mean absolute error
2. MRE: mean relative error where AL computational and AL experimental are the anodic limits calculated or measured experimentally, respectively, whereas N = 5 is the number of ionic liquids considered in this investigation. In Table 10, the summary of the MEA and MRE are reported for the eleven different computational/thermodynamic models.
It can be noted from Figure 2 that, for the estimate of the oxidation potentials of ionic liquids, the adoption of the MP2 method in vacuum mimicking a vertical transition from the anion largely overcomes all the other computation/methodological approaches in terms of accuracy. Lowering the level of theory from MP2 to B3LYP leads to large errors in the estimate of the oxidation potentials. Of course, the experimental anodic stability of different ionic liquids slightly depends on the cation; however, the spread of these values is extremely small when compared to the large differences introduced by the calculations with different computational methods, especially in the case of TFSI-based ILs (see Figure 2). Overall, the accuracy of the MP2 predictions is surprising since the "vertical transition approximation from ion in vacuum" is crude compared to the complexity of a real case. From a physical point of view at the atomic scale, an irreversible oxidation of ionic liquid is a vertical ionization of an adsorbed anion at the electrolyte/electrode interface. Thus, the "vertical transition approximation from ion in vacuum" neglects two major real physical effects: (a) the partial solvation provided by the molecules/ions surrounding the reagent/products in proximity to the electrode, and (b) the ZPE of reagents and products. One may speculate that the partial compensation of solvation/ZPE from reagent/products results in precise predictions beyond expectations.
On the other hand, once the level of theory decreases from MP2 to B3LYP, accuracy drops. Furthermore, all the other adopted approximations fail in their attempts to compensate the worse electronic structure description modeled by the hybrid-DFT functional B3LYP compared to the MP2 perturbational method, although there is better mimicking of solvation and ZPE.
The partial solvation effect cannot be modeled by (a) adopting a continuous model solvent (polar media), (b) computing the entire ionic couple or (c) combining continuous model solvent and the ionic couple. In all approaches, reagents are too stabilized, thus leading to unrealistic large oxidation potentials.
Turning to the impact of the ZPE variations during the transition, one may consider that it is computationally feasible but with a weak physical meaning. In fact, whereas the ZPE of the reagent can be very easily and accurately computed, the neutral product of the vertical transition is not a minimum structure and has a large majority of imaginary vibrational frequencies. One possible way to circumvent this drawback is to mimic oxidation through an adiabatic transition. In fact, after structural relaxation, the oxidized product can be a vibrational minimum with all positive vibrational frequencies. However, this approach leads to an unrealistic stabilization of the electron energy of the oxidation product, thus leading to unrealistic small oxidation potentials.
One may note that the simultaneous occurrence of these two computational artifacts (i.e., overstabilization of the reagents by the continuous solvation model and the overstabilization of the products by following the adiabatic transition approach) can partially compensate each other. In fact, both the estimates of the oxidation potentials using adiabatic transitions in a continuous solvent from ions show errors of 9% or 0.45 V, which is smaller in comparison to the predictions from vertical transition from ions in vacuum at the same level of theory (B3LYP, i.e., −14%, −0.71 V).

Conclusions
The anodic and cathodic limits of selected ionic liquids were investigated computationally by means of various methods to compare them with the available experimental results. In previous literature, different techniques based on calculations on single ions, as well as on ionic couples, were proposed. In the present investigation, methods involving calculations on ionic couples either overestimate electrochemical stability or give values of the anodic and cathodic limits which tend to be independent of the specific couple. Techniques based on calculations on single ions seem to be more reliable, especially those based on vertical electronic transition computed with the MP2 functional. DFT methods relying on the B3LYP functional either underestimate the ESW, when calculations are performed in vacuum, or largely overestimate the ESW, when ions are placed in a polarizable medium. Overall, the calculations considering vertical electronic transitions in vacuum of the single ions using the MP2 theory are in the best agreement with the previous experiments. This fact opens new perspectives for future validation of the model for the electrochemical stability of other liquids already experimentally investigated, as well as for the screening of anions and cations that have not already been studied from an experimental point of view.
Author Contributions: Conceptualization, A.P. and S.B.; methodology, A.P. and S.B.; validation, A.P.; investigation, A.P.; writing-original draft preparation, A.P. and S.B.; funding acquisition, A.P. and S.B. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors would like to acknowledge the financial support from the European Union Horizon 2020 research and innovation programme within the Si-DRIVE project; grant agreement No. 814464.