Which ab initio calculation methods describe the charge transfer in ionic liquid vapors? Case study on EMIM-OTF reveals the frontier orbitals and the ionization threshold

Abstract

Time-of-flight mass-spectrometry (TOF-MS) measurements of 1-ethyl-3-methylimidazolium trifluoromethanesulfonate ([EMIM][OTF]) vapors have been conducted in order to gain insight into the VUV (vacuum ultraviolet) absorption thresholds. The experimental data were interpreted with hybrid density-functional theory (DFT), time-dependent DFT (TD-DFT), and several flavours of many-body perturbation theory (MBPT): second-order MBPT, Hedin's GW and Bethe–Salpeter equation (BSE) calculations. To our knowledge, this is the first comprehensive comparison of DFT functionals for an ionic liquid (IL) vapor. More than 60 exchange–correlation approximations, including many range-separated hybrids (RSHs) and optimally tuned variants, were assessed. The theoretical methods predicted a wide range of energies; only the experimental TOF-MS data allowed the truly predictive approaches to be identified. The experimental spectra reveal charge transfer (CT) states at an optical gap of 6.75 eV and a maximum ion yield signal at 8.3 eV. Combining the experiment and theory enabled the re-examination of the ion-pair geometry and the reconstruction of the frontier orbitals of [EMIM][OTF]. The HOMO and LUMO energies were found to be −9.05 eV and −2.30 eV, respectively; energies for HOMO−1 and LUMO+1 were also determined. Although a ΔSCF calculation with the non-hybrid M11-L functional gives a reasonable ionization threshold, ≥50% exact Hartree–Fock (HF) exchange is required to correctly describe both the optical gap and the ultraviolet photoemission spectrum (UPS). The RCAM-B3LYP functional and the Bartlett's QTP family of functionals showed good performance in predicting the CT excitations from TD-DFT calculations. No method was able to predict the proposed LUMO energy, though the M08-SO functional was closest overall. Within the limitations of the present work, the global hybrid PBE-2X functional demonstrates the most consistent overall performance. Its range-separated hybrid (RSH) analogue CAM-PBE50* also performed well for the calculation of CT excitations.

---

1. Introduction

Ionic liquids (ILs) are molten organic salts that melt below 100 °C. Their negligible vapor pressure, broad liquid range, high thermal and electrochemical stability, and strong solvating power have promoted their use in electrodeposition, heterocyclic synthesis, and diverse electrochemical processes. Accurate knowledge of their frontier-orbital energies and valence-band composition is essential for both fundamental insight and device design.^1–6

Yet the gas-phase electronic structure of many ILs, including [EMIM][OTF], remains poorly understood. For example, HOMO–LUMO gaps reported for gas-phase [EMIM][OTF] span almost 6 eV: in 2024 Wang et al. published a 4.5 eV gap using B3LYP and a 10.2 eV gap value calculated with the ωB97X-D functional.⁷ Some authors use the MP2 method,⁸ while other groups recommend the ωB97X-D and M06-2X functionals.⁹ However, just very recently Fogarty et al. continued to use the old B3LYP functional for IL calculations.¹⁰

This persistent inconsistency primarily arises from two factors. First, the lack of reliable experimental gas-phase data has led to the simultaneous publication of these widely differing gap values of 4.5 and 10.2 eV. Even when experimental data are available, their interpretation is not always straightforward. For example, there are claims of adducts in the vapor phase photoemission spectra of [EMIM][BF₄].² Unexplained low-energy tails extending down to ∼2 eV have been observed in the UPS spectra of liquid-phase [EMIM][TFSI].^11–13 Second, although the 4.5 eV gap predicted for gas-phase [EMIM][OTF] is almost certainly an underestimate, it is still widely used because it conveniently matches the electrochemical stability window, which is 3.9–4.1 eV experimentally.⁷ This superficial agreement sustains the uncertainty surrounding the true HOMO and LUMO levels and is unlikely to be resolved until experimental data and high-level ab initio calculations converge.

Beyond the gap values themselves, several other challenges complicate the ab initio modeling of ILs. The localization of the HOMO state in IL ion-pairs is non-trivial and depends strongly on the specific cation–anion combination. It may reside on the cation or the anion, or be delocalized across both.^1,14,15 Moreover, the identification of the lowest-energy conformer is highly method-dependent, further complicating ab initio calculations.¹⁶

Thirdly, while hybrid DFT calculations provided very good qualitative agreement with most experimental UPS spectra a systematic shift was still necessary. There is no good way of a priori knowing what the required shift is for a particular hybrid DFT functional. Although using shifted energies can lead to good agreement with experiment,^10,14 a calculation method that requires no systematic shift would be much more preferable. To achieve this, a comparison between experimental data and different calculation methods is needed.

Fourthly, in some ILs, such as [EMIM][BF₄]¹⁷ and [EMIM][DCA],¹⁸ valence photoemission involves electronic relaxation, charge transfer, or shake-up processes. These effects, possibly due to the dual ionic and aromatic character of these systems, complicate the interpretation experimental spectra and quasiparticle energies.

In contrast, the electronic structure of [EMIM][OTF] was captured very well in previous studies and showed no evidence of any charge transfer during photoemission.^2,14 This supports the validity of the ion-pair approximation for the vapor phase and confirms that hybrid DFT provides a reliable description of its ground-state electronic structure.

The next logical step is to investigate excited states. Regrettably, excited states and charge transfer energies in ionic liquids remain virtually unexplored. To address this gap, we measured energy-dependent TOF-MS mass spectra of the [EMIM][OTF] vapor under VUV photoexcitation, thereby gaining insight into its excited and unoccupied electronic states. This approach builds on our previous investigations of TIY (total ion yield) and PIY (partial ion yield) for [EMIM][BF₄] and other ILs.^19,20

The paper begins with an analysis of the electronic structure and charge transfer processes in [EMIM][OTF], followed by an overview of the computational methods employed. Ab initio description of electronic excitations generally relies on TD-DFT, which is known to struggle with charge transfer type excitations. Moreover, the intrinsically charged nature of IL ion-pairs poses an additional challenge for ground state calculations, as many computational methods suffer from excessive charge delocalization or over-localization.

The primary goal is to identify ab initio methods that best reproduce experimental data, thereby providing guidance for future studies of ionic liquids and related systems. Although drawing general conclusions from a single case is difficult, no comparative data for other ionic liquids currently exist.

Despite its apparent simplicity, [EMIM][OTF] will be shown to surprisingly resist quantitative description by most current ab initio methods!

2. Experimental

2.1 TOF-MS measurements

The [EMIM][OTF] under study was purchased from Iolitec GmbH. It had a stated purity of 98% and was transparent. The TOF-MS measurements were carried out using the ion TOF spectrometer of the gas-phase end station²¹ at the FinEstBeAMS beamline of the MAX-IV laboratory (Lund, Sweden). The beamline is equipped with a collimating grazing incidence plane grating monochromator and covers the excitation photon energy range from 4.5 eV to about 1300 eV.^22,23

The energy resolution of the beamline was set to about 2 meV during the experiment. A LiF optical window was used to filter out any higher diffraction orders and scattered high-energy photons while offering transparency up to 10.8 eV photon energies.

After inserting the IL into the experimental chamber, heating of IL at around 120 °C for several hours was performed to remove residual water from the chemical. Despite the drying procedure, water was still present in the sample and it evaporated to the measurement chamber. The temperature of the effusion cell and the quartz crucible holding the IL was increased to 180 °C in order to evaporate the IL. The pressure in the experimental chamber during the measurements was 3 × 10⁻⁶ mbar while the experimental chamber base pressure was 2 × 10⁻⁷ mbar or better.

Unfortunately, due to technical issues with the TOF-MS spectrometer the background is elevated and some prominent gases also appear in the spectra, most notably nitrogen and water. The signal to noise ratio is worsened due to this background. However, because the background is mostly independent of the photon energy, it does not significantly influence the main results.

The TOF-MS spectra were recorded for a fixed acquisition time at each excitation energy. The resulting spectra were then smoothed under a moving-average window of 20 ns. The PIY of the EMIM cation was obtained by subtracting a linearly sloped background, defined by the peak edges at TOF 17 625 and 18 650 ns, and integrating the background-corrected signal over the same TOF range. All these integrated counts were then subsequently normalized to the beamline mirror current to account for the varying photon flux at different excitation energies. A different PIY, representing an alternative fragmentation pathway for the higher mass fragments, is shown in the graphical abstract. An example of a recorded TOF-MS spectrum is shown in Fig. 2.

2.2 Computational details

Ab initio ground state and TD-DFT calculations were performed using the MolGW (version 3.3)²⁴ and NWChem (versions 7.0.0 and 7.2.3)²⁵ software packages. Meta-GGA functionals, such as M06-2X and SCAN, were exclusively evaluated using NWChem. Both software packages utilized the LIBXC library (version 7.0.0) to access a wide range of DFT functionals.²⁶ The corresponding LIBXC functional identifiers are listed in Table 3.

Full matrix TD-DFT calculations, excluding the Tamm–Dancoff approximation, were performed using MolGW. These calculations included only singlet excited states and employed the linear response formalism to electromagnetic perturbation to compute the spatially averaged photoabsorption cross section.²⁴

MolGW was also employed for MBPT calculations within the GW approximation as well as for solving the BSE. All presented calculations—ground state, TD-DFT, GW, and BSE—were performed using the same molecular geometry, shown in the SI. The screened Coulomb interaction was computed within the random phase approximation (RPA). For evaluating the self-energy, both the contour deformation technique and the spectral decomposition method were tested, with the latter ultimately being used due to comparable results and improved computational efficiency. The BSE kernel was constructed using the RPA and the static screening approximation, while vertex corrections were neglected to reduce computational cost.

To simulate the spectra shown in Fig. 6 and 7, transitions from both BSE and TD-DFT calculations were broadened using a Gaussian function with a full width at half maximum (FWHM) of 0.2 eV. A broader 0.25 eV FWHM Gaussian was applied for spectra associated with the ‘3 + 2 model’ (described below), as well as for the density of states (DOS)-type spectra shown in Fig. 3. This broadening provides a simplified approximation of real spectral broadening effects, such as thermal and vibronic.

The molecular structure of the conformer was not relaxed under the functional of interest before performing the ΔSCF calculation. In most cases, the ΔSCF energy remained largely insensitive to geometry optimization, as the energetic contributions from structure relaxation tended to cancel out when computing the energy difference between the ground state and the ionized cation. However, in some GGA functionals (HCTH407+, HCTH-A, HLE16) there was a substantial ΔSCF energy change after the structure relaxation. This is due to the fact that these functionals predict different (wrong) bond lengths and equilibrium structure compared to most hybrid functionals.

All ground-state calculations were carried out using both Dunning's aug-cc-pVDZ and aug-cc-pVTZ basis sets. The HOMO energy was found to be largely insensitive to the choice of basis set. The third column of Table 3 provides an evaluation of whether each ΔSCF energy is converged with respect to both geometry (relaxed vs. unrelaxed) and basis set (aug-cc-pVDZ vs. aug-cc-pVTZ). A convergence threshold of 0.2 eV was used. Except for the BSE profiles shown in Fig. 6, all reported computational data were obtained using the aug-cc-pVDZ basis set, unless stated otherwise. This choice was made due to the significantly higher computational cost of TD-DFT calculations with the aug-cc-pVTZ basis set, which becomes prohibitive for larger ILs.

To accelerate the computation of Coulomb and exchange integrals, resolution of identity (RI) approximations, charge density fitting, or auxiliary basis sets were employed where available. These approximations introduced an average change of only 23 meV to the HOMO energy, with a maximum absolute deviation of approximately 60 meV.

Although an extensive conformer search was performed in our previous work, a revised [EMIM][OTF] ion-pair conformer has been found and its structure was optimized. Initial MP2 calculations with the aug-cc-pVDZ basis set predicted the revised conformer to be 35 meV lower in total energy than the previously considered conformer. In contrast, calculations using various hybrid DFT functionals with the same aug-cc-pVDZ basis set indicated that the previous conformer was 200–320 meV lower in energy.

This energetic ordering reversed upon employing the larger aug-cc-pVTZ basis set. Optimization of both conformers with the aug-cc-pVTZ basis set using nine different DFT functionals showed the revised conformer to be, on average, 45 meV lower in energy. This energy difference is relatively small and comparable to the thermal energy (kT) at the evaporation temperature. The HOMO energy of the revised conformer is more negative, indicating a more strongly bound orbital, consistently across both aug-cc-pVDZ and aug-cc-pVTZ basis sets. While the conformer search remains somewhat incomplete, the main focus is on comparing different calculation methods, and the conformer selection is not expected to influence the presented results.

No counterpoise correction was applied to account for basis set superposition error (BSSE) in the MP2 calculations. Additionally, zero-point energy and vibronic effects were not included in any of the calculations, though the former is not expected to significantly impact the conclusions presented.²⁷

Aside from spectral broadening and rounding, all reported data are presented without empirical adjustment, except the photoabsorption profiles calculated by the BSE@HF shown on Fig. 6 and the BSE@LC-BOP shown on Fig. 7, which have been shifted to align the main experimental PIY peak at 8.3 eV. This should make the visual comparison between the experimental spectra and the simulated BSE photoabsorption profiles easier.

The rounded ab initio calculation results without any other adjustments are shown in Table 3.

3. General background of EMIM-OTF

A revised [EMIM][OTF] ion-pair conformer was selected and its structure was optimized. This structure along with its HOMO and LUMO is depicted in Fig. 1. The most significant structural difference compared to the conformer used in our previous work¹⁴ appears to be the proximity of the triflate anion's fluorine atoms to the imidazolium ring, whereas they were previously oriented away from it. The HOMO remains predominantly composed of oxygen 2p orbitals from the triflate anion (see Fig. 1).

Fig. 1 The calculated HOMO and LUMO orbitals of the [EMIM][OTF] ion-pair (left) and the “3 + 2” orbital scheme with the transition energies (right). The transitions obtained from this model are shown in Fig. 6 and 7.

Both the revised conformer and the conformer used in our previous work reproduce the experimental UPS spectrum well, but can differ more in their photoabsorption profile. This is likely because VUV excitation involves an anion-to-cation charge transfer transition, making the anion orientation more critical in this excitation regime. Some functionals agree better with the experimental PIY when their TD-DFT calculations are based on the old conformer (CAM-QTP-00, QTP17, CAMY-BLYP) while others compellingly prefer the new: CAMH-B3LYP, wb97X-D, CAM-PBE50*, PBE-2X. As previously noted, the total energy difference between the two conformers is very small indeed (45 meV). These findings collectively underscore the difficulty in predicting lowest energy IL conformer structures.

The connection between the experimental fragmentation spectrum measured by TOF-MS and the ab initio photoabsorption calculation may need some elucidation. The underlying assumption is that the intact cation peak in the TOF-MS spectra (see Fig. 2) also provides information about excited states below the ionization threshold. This is due to the fact that the cation–anion binding energy is smaller than the CT energy. The cation–anion binding energy, defined as the energy required to separate the ion-pair into an isolated EMIM cation and a triflate anion, is estimated by our ab initio calculations to be around 3.4–3.9 eV. This value is smaller than the first observed CT feature in the PIY at 6.75 eV. [EMIM][OTF] is assumed to behave similarly to [EMIM][BF₄], for which the appearance energy of the EMIM cation was shown to be lower than the ionization energy of the ion-pair.²⁰ However, the reader has to keep in mind that the experimental PIY spectrum is not directly be proportional to the photoabsorption cross-section, as the fragmentation probability and branching ratios between different channels may vary non-linearly with excitation energy.

Fig. 2 The experimental TOF-MS curve at the excitation energy of 8.1 eV.

Further insight into ILs themselves and the VUV excitation process can be obtained by considering the special limiting case of a donor–acceptor complex. In an ideal donor–acceptor system, direct chemical interactions between the donor and acceptor entities are considered negligible.

In the context of ionic liquids, the donor–acceptor model becomes applicable when considering a neutral IL cation positioned adjacent to a neutral IL anion. This conceptual arrangement represents the state immediately preceding charge transfer—the point at which an electron transitions from the cation to the anion, thereby forming the charged ion-pair or “charging” the ion-pair. The experimentally observed CT process can then be viewed as the reverse of this event, where an electron is transferred back from the anion to the cation. Following the latter CT, the resulting EMIM-OTF complex can readily dissociate.

However, this arrangement is only a rough approximation of a true donor–acceptor pair. In addition to electrostatic interactions, there are also significant chemical interactions between the cation and anion due to their close proximity. Despite this limitation, the Mulliken formula provides the necessary components for estimating the energy required for electron transfer from the donor to the acceptor:

In this formula we have the ionization energy of the donor (cation) and the electron affinity of the acceptor (anion). This energy release E_ioniz − E_aff is partially compensated by the coulombic attraction 1/R occurring after the CT. The experimental CT process is observed as two distinct transitions: the HOMO–LUMO transition at 6.75 eV and the 3× HOMO−1 to LUMO transition at 7.9 eV.

A ‘back of the envelope’ calculation can be performed to get a rough estimate of the terms in the Mulliken formula. Using the top-performing PBE-2X DFT functional via a ΔSCF approach, the energy difference E_ioniz − E_aff is estimated to be 9.5 eV. However, this value is likely an overestimate, serving as an upper bound, since it neglects the chemical interactions between the cation and anion. By subtracting the cation–anion binding energy to account for these interactions, a CT energy of 6.0 eV is obtained. This value instead underestimates the CT energy as we have overestimated the cation–anion binding energy by using the binding energy of the charged species instead of the neutral radicals. The same 6.0 eV value can be further supported by a similar ΔSCF calculation. Namely, the total energy gain from an electron transfer from the cation to the anion is 8.3 − 2.3 = 6.0 eV. Here the two terms represent the energy gains of the isolated anion and isolated cation radicals, respectively.

Using a cation anion distance of 4.1 Å, which corresponds to the distance between the carbon of the OTF anion and a nitrogen of the EMIM cation, the electrostatic energy of the ion-pair is roughly estimated to be about 7.75 eV. This corresponds closely to the weighted average of the experimental CT energy, which is at 7.81 eV for the proposed (CT) LUMO-only transitions.

The close agreement between the ab initio estimated minimum Mulliken CT energy of 6.0 eV and the experimental minimum of 6.75 eV; the experimental mean CT energy of 7.81 eV, the donor–acceptor model mean CT energy of 7.75 eV and the estimated electrostatic energy of 7.75 eV supports the conclusion that the CT process is primarily electrostatic in nature. Therefore, ILs represent an interesting case intermediate between a fully coupled electronic system and an ideal donor–acceptor complex.

The HF exchange is mostly responsible for the E_ioniz − E_aff energy as it promotes charge localization and allows for the buildup of the electrostatic energy between the cation and anion. The CT energy and degree of charge separation are dependent on the HF exchange content in the calculation method. Standard DFT functionals on the contrary tend to delocalize the cation and anion charge due to the so-called self-interaction error, which will be discussed later in more detail. This explains why HF exchange is crucial for accurate calculations of ILs and why results strongly depend on its content (e.g., through parameters α and µ), as will be shown later.

The experimental Fermi energy of [EMIM][OTF] vapor was estimated to be about 8.4–8.5 eV (see Fig. 3 and¹⁴). It is important to note that this so-called “Fermi energy” refers to the minimum photoelectron binding energy. It will be shown that this is not the vertical ionization energy nor the true HOMO energy of the IL. Our previous work estimated the HOMO energy at −(8.0–8.2) eV using methods that are now outdated.¹⁴ The G₀W₀@PBE0 calculations by Kahk et al. yielded a (true) HOMO energy value of −9.06 eV.² Kanai et al. reported an energy gap of 8.1 eV for the similar liquid-phase [BMIM][OTF] IL.²⁸ Seymore determined the liquid-phase ionization energy of [BMIM][OTF] to be 9.4 ± 0.4 eV.⁴ Fogarty et al. performed valence band XPS measurements on [BMIM][OTF] among other ILs¹ and found that B3LYP-D3(BJ) calculations using a lone ion-pair solvation model, shifted by approximately 1.7 eV, could reproduce a wide range of non-resonant valence photoemission spectra.¹⁰ They estimated the liquid phase [EMIM/BMIM/OMIM][OTF] ionization energy to be about 7.3 eV.

Fig. 3 Comparison between the experimental UPS spectrum (black curve) and the DOS type spectra of several DFT functionals. The UPS spectrum is reproduced from ref. 2 and 14. The DOS type spectra were obtained by uniformly broadening the calculated ground state MOs. Four DFT functionals with low absolute energy errors are shown: LC-BLYP-EA (topmost red curve), CAM-LDA0 (blue curve), PBE-2X (green curve) and M08-SO (bottom brown curve).

These five estimations were partly or fully derived from experimental data and although the liquid and gas phase properties are different, the predicted HOMO energy range is −(7.3–9.4) eV. In comparison, HOMO energies calculated in this study span an even wider range: −5.4 to −13.2 eV.

A proposed simplified energy level scheme for the [EMIM][OTF] ion-pair, derived by combining experimental data and ab initio calculation results, is shown in Fig. 1. The figure shows 3 occupied states (HOMO, HOMO−1, HOMO−2) labelled A, B and C; and 2 unoccupied states (LUMO and LUMO+1). Furthermore, the HOMO−1 is triply degenerate, while HOMO−2 and the LUMO+1 are doubly degenerate. The scheme also shows some transitions from the occupied states to the LUMO(s) and to the continuum. This energy level scheme and its derived transitions are collectively labeled the ‘3 + 2 model’.

The optical gap is defined by a neutral excitation that corresponds to the difference between the energies of the lowest dipole-allowed excited state and the ground state.²⁹ The fundamental gap is defined by charged excitations as the difference between the first ionization potential and the first electron affinity.²⁹ The charge transfer excitation is where photon absorption induces the transfer of an electron from a donor unit to an acceptor unit. Such transitions are identified by their dominant character of excitation from a donor-localized orbital to an acceptor-localized orbital.

Based on this information, we can draw the following conclusions about the electronic structure of the [EMIM][OTF] vapor. The optical gap of [EMIM][OTF] is the lowest energy excitation at 6.75 eV (peak A₁ on Fig. 6 and 7). Because the HOMO of [EMIM][OTF] is located on the anion and the LUMO levels are on the cation (see Fig. 1) then the transition responsible for this optical gap is by definition a CT transition. This CT transition in [EMIM][OTF] should, in principle, be similar to the proposed CT state of [EMIM][DCA],¹⁸ i.e., from the anion to the (non-bonding) π* orbitals of the EMIM cation (associated with the C=C double bond region). Indeed, the analysis of the natural transition orbitals (NTO), obtained via singular value decomposition of the TD-DFT transition density matrix, indicates that the final state of the CT excitations mostly corresponds to the LUMO. For example, when using the ωB97X-V functional, the dominant NTO pair describes 93.4% of the CT excitation at 6.7 eV, with the corresponding particle NTO exhibiting 99.8% LUMO (depicted in Fig. 1) character.

The experimental UPS spectrum of [EMIM][OTF] vapor is presented alongside the experimental PIY spectrum and the ab initio calculated photoabsorption profiles (Fig. 7). The experimental PIY spectrum shows transitions to the LUMO(s), which start from the gap and peak around 8.3 eV.

The UPS photoionization spectrum shows transitions to the continuum states. Its first peak is located at around 9.1 eV and the spectrum reaches maximum intensity around 10.55 eV. The energy difference between the PIY and UPS maxima is about −2.25 eV and this corresponds well with the proposed LUMO energy of −2.3 eV.

The energy scale of the experimental UPS spectrum is reversed compared to the conventional plotting style and furthermore shifted by 2.3 eV from the excitation energy x-axis below. This shift visually aligns peaks corresponding to the same occupied electronic states involved in photoemission and photoabsorption (labeled A–C in Fig. 1), thereby enabling direct visual comparison of initial state features. For example, this aligns the UPS photoionization peak I_a with the peak A₁ of PIY and the UPS I_b, I_c doublet with the corresponding PIY peaks B₂, C₁ (see Fig. 7).

4. Theoretical methods

4.1 TD-DFT

Good overviews of TD-DFT can be found in ref. 30 and 31. Certain types of optical excitations, notably CT excitations typically defy proper description using conventional approximations employed in TD-DFT calculations.²⁹ Most functionals produce spurious, low-energy CT states, which often contaminate TD-DFT calculations in large molecules.^32,33 The TD-DFT description of Rydberg states or diffuse valence states can also be rather poor due to the incorrect asymptotic behavior of the available functionals.^33,34

The TD-DFT predicted excitation energies are shown in Table 3. Note that, indeed, most functionals have at least one spurious CT excitation and only the energy of the first spurious CT transition is shown in Table 3. These spurious CT excitations are determined by aligning the calculated CT transitions corresponding to the first excitation (A₁₂) and main excitations (B₁₂, C₁₂) and identifying the leftover low-energy excitations. The overall intensity of the predicted spurious CT transitions varies between close to zero (CAM-PBE50*) up to about to the first CT intensity predicted from the BSE@HF method (HCTH407+, for example). Even the RSH functionals that do not produce spurious CT excitations tend to overestimate the intensity of the CT excitation at 6.75 eV. For example, QTP17, CAM-PBE50*, ωb97X, PBE50, PBE-2X, the CAM-QTP family, all have this deficiency (see also Fig. 6 and 7 and the ‘3 + 2 model’). Only functionals that produce a spurious (first) CT excitation have an A_1,2 intensity ratio closer to the expected 1 : 2 value. For example, LDA0, PBE0 and B3LYP show roughly a (1 : 1.1)–(1 : 2.5) ratio for the A_1,2 peaks (not shown).

4.2 MBPT

Dyson equations of MBPT are based on the one-particle propagator and it would be highly valuable if they could be solved. This is because the Dyson equations provide all quasiparticle energies, which among other excitations have all the correct ionization potentials:³⁵Here T is the kinetic energy operator, v is the sum of the external (nuclear attraction) and Hartree potentials, ε_i is the quasiparticle energy, φ_i is the Dyson orbital and Σ_i is the non-local electron self-energy that captures relaxation, correlation, self-interaction corrections i.e. it contains all the many-body interactions experienced by the single particle. It is this Σ_i term, that is very difficult to calculate.

The Dysons orbitals form an overdetermined or overcomplete non-orthogonal set and contain too much ‘excess baggage’.³⁶ The single particle Green's function is the next best tool and it can be derived from the Dyson orbitals:³⁷

Here, only the electron addition (ψ_α) and removal Dyson orbitals (ψ_β) contribute. The Green's function is therefore highly important in MBPT³⁸ and it's poles correspond to the electron addition and removal energies (eqn (3)).

Under conditions where there is a dominant solution of the Dyson eqn (2), one with the largest residue, the Green's function can be approximated by simply just taking this pole into eqn (3) and neglecting satellites, for example. The Green's function constructed this way is referred to as the quasi-particle Green's function or just the Green's function as opposed to the more general single particle Green's function.³⁷

If Ẑ is the nuclear potential operator of all the nuclei in the system then the many-body N-particle ground state density and energy can be expressed through the Dyson orbitals:

The ground state electronic density is a functional only of the ‘occupied’ Dyson orbitals {ψ_β}.

The Kohn–Sham theorem proves the existence of (single determinant) orbitals for the calculation of the ground state density. However, more interestingly, the DFT ground state energy expression is somewhat analogous to the eigenvalue expression of eqn (2) if the exchange–correlation potential v_xc would approximate the self-energy Σ:³⁷

Now it is possible to sketch the approximation ladder for the COT principles. Starting from the most general Dysons orbitals (eqn (2)), then simplifying further to the single-particle Green's function (eqn (3)), then to the quasiparticle Green's function and finally reaching the Generalized Kohn–Sham system (GKS) orbitals (eqn (5)). This is the mathematical preamble needed for understanding the COT principles and GKS, which will be discussed further below, under modern DFT.

However, in practice, Hedin's GW method³⁹ has become the de facto standard for the calculation of quasiparticle energies and approximating solutions to the Dyson equations. In the GW approximation the electron self-energy is expressed using the single particle Green's function and the screened Coulomb potential (W):⁴⁰

This Σ can be considered to be a nonlocal and dynamic effective potential. The initial Green's function G₀ is found from the ground state using the (G)KS orbitals instead of the Dysons orbitals of eqn (2) and (3).

Usually the GW@RPA approximation is used for GW calculations. This means that the RPA is used to calculate the inverse dielectric screening function ε⁻¹ that is convoluted with the bare Coulomb potential v to obtain the screened Coulomb potential W₀:⁴¹

In this one-shot G₀W₀ calculation, basically only the screening of the Coulomb interaction is accounted for. It does not account for the second order exchange (SOX) i.e. screening of the exchange interaction. This G₀W₀@RPA is a simple and somewhat crude approximation, but performs remarkably well for HOMO values due to cancellation of higher order terms.^42,43

A recent review confirmed that the G₀W₀ method provides a very good estimation of the ionization energies of molecules.⁴³ Indeed, it also performs well for our IL ion-pair (see Table 3).

In practice, GW means not just one, but a family of methods with various levels of self-consistency⁴⁴ and various correction terms. Very often, the DFT functional, that is used to generate the orbitals for the initial Green's function G₀, must be chosen carefully.

4.3 Bethe–Salpeter equation

The BSE formalism is a MBPT method that provides realistic neutral excitation energies, in a similar manner as TD-DFT does within DFT. The theoretical foundation of the BSE approach is completely different from that of TD-DFT as it is not based on the electron density, but on the one-particle Green's function (eqn (3)). The KS eigenvalues are replaced with the GW quasiparticle energies and the electronic Hessian does not depend on the f_xc kernel but on the BSE kernel, which effectively screens the Coulomb interaction for a given frequency.⁴⁵ However, the BSE working equations can be cast into a very similar form to TD-DFT.

BSE builds on top of a GW calculation by adding electron–hole pair or excitonic effects to the GW HOMO–LUMO gap.⁴¹ This is a third order calculation or analogous to a second order perturbation calculation, because the HOMO–LUMO gap derived from the GW calculation is itself a perturbative evaluation of the DFT gap.

The excitonic effect i.e. the attraction between the excited electron and the hole left behind, always leads to stabilization near the gap. This is why the BSE optical gap is always smaller than the GW gap.⁴¹

In our case, the performance of the BSE is evident for the LC-PBEOP, LC-BLYP, LC-BOP functionals, for example, where it successfully introduces exited CT states below the ionization threshold (see Table 3). In addition, the BSE excitation cross sections have good qualitative agreement with the experimental PIY.

The BSE calculation tends to systematically underestimate the energies of the CT excitations i.e. overestimate the excitonic binding energy. For example, in case of the LC-wPBEh-whs functional the TD-DFT calculation slightly overestimates the CT excitation energies, but BSE instead underestimates (see Table 3). No BSE calculation is able to exactly match the experimental CT energies and similarly to TD-DFT, it also tends to overestimate the A_1,2 peaks.

One reason for this is due to the basis set limitation. Improving the basis set to aug-cc-pVTZ will decrease the G₀W₀ quasiparticle energies and the BSE excitations correspondingly. To investigate this effect further, the BSE calculation with four functionals was performed in the aug-cc-pVTZ basis set. In case of the LC-BOP functional, which has the smallest CT excitation error of about 0.25 eV, the BSE energy scale shift disappears and there is very good agreement with the experimental PIY (Fig. 6). Also, the energies of the BSE photoabsorption profiles calculated by the LC-wPBE08-whs, LC-VV10 and LC-BLYP-EA functionals improve by 0.2–0.3 eV upon a similar basis set upgrade. However, in the latter case the spectrum still does not match the experiment. This is in spite the fact that the G₀W₀ calculation now already underestimates the HOMO energy by about 0.5 eV. Therefore, the LC-BOP, LC-VV10 and LC-wPBE08-whs BSE calculations give the right answer for the wrong reason, their GKS orbital and G₀W₀ self-energy values are underestimated in the aug-cc-pVTZ basis set thereby compensating for the tendency of BSE to underestimate the CT excitation energies.

The reason for the systematic underestimation of the CT energies is probably caused by the overscreening of the HF exchange term, which arises from neglecting higher-order derivatives of the screened exchange.⁴⁵ Future studies that add dynamic excitations, use higher order vertex correction terms or are based on optimally tuned RSH functionals, may be needed.

4.4 Modern (RSH) DFT

Most DFT failures in molecular systems can be traced to the so-called self-interaction error (SIE): an artificial interaction of the electrons with themselves. SIE is also known as the delocalization error, the fractional electron problem, missing derivative discontinuity or just curvature error. SIE can be assessed by doing calculations for fractional occupation numbers. A straight line is expected for total electronic energies when the occupation number is varied smoothly between integer numbers.⁴⁶ However, it is well known that SIE results in too low energies for the intermediate species with fractional occupations.⁴⁷

Systems where the DFT SIE is clearly evident are the X₂⁺ dimer, like H₂⁺ or He₂⁺. Standard DFT will make each fragment incorrectly have a fractional charge at large separation and too low total energy. This can occur in one-electron systems such as H₂⁺, many-electron radicals and also in neutral molecules. (Semi)local exchange functionals do not cancel the spurious coulombic self-interaction, leading to a self-repulsion of localized one-electron densities and thus to delocalization.⁴⁸ For the same reason, SIE also leads to too low transition energies in charge transfer complexes.

The Coulomb and exchange energies of a single electron need to cancel to avoid self-interaction. Thus, there is a deep connection between the SIE and correct exchange.

One way to correct the SIE is to mix in some correct Hartree–Fock exchange to the DFT functional. This works because of the opposite tendency of HF, as it tends to overlocalize electrons and overestimate the bandgaps. By mixing DFT with HF exchange we can get a (partial) cancellation of these errors. In fact, the SIE cannot be corrected in non-hybrid functionals such as meta-GGAs.⁴⁹

More recently, these global hybrids have been replaced by range-separated hybrid (RSH) functionals, which have a fraction of long-range (LR) exact HF exchange that typically smoothly increases to 1. The exact LR HF exchange should balance the Hartree energy and cancel the self-repulsion of the electrons.

E_x = αE^SR_HFx + (α + β)E^LR_HFx + (1 − α)E^SR_DFTx + (1 − [α + β])E^LR_DFTx (7)

The α parameter controls the fraction of short-range (SR) HF exchange (HFx). The DFT exchange (DFTx) is usually some derivate of the homogenous electron gas exchange energy in the GGA or LDA form. Another term that is used with these kinds of functionals is long-range-corrected (LRC or LC), which means that there is 100% HF exchange for long-range interactions (α + β = 1). This full exact LR exchange cancels the spurious coulombic self-interaction but it is mostly possible in finite systems, like our ion-pairs considered here. For periodic systems, such as solids, or partially periodic systems such as surfaces, nanowires or nanotubes this 100% LR exchange treatment can usually not be used.

In the initial LRC hybrid DFT treatments, the contribution of HF exchange was dampened to zero at short range (α = 0). An important development of modern DFT was made by Yanai et al. when they introduced the Coulomb-attenuated B3LYP (CAM-B3LYP) functional in 2004.⁵⁰ Exchange is partitioned into SR and LR components that are scaled by the error function such that there is seamless transition between the two regimes and they can be calculated independently based on the following identity:

Besides the error function erf, the Yukawa functional is a similar way of separating the exchange integrals. The latter approach is represented by the CAMY-B3LYP, LCY-BLYP, CAMY-BLYP and LCY-PBE functionals in Table 3.

The range-separation is determined by the parameter µ (often also denoted by γ), where 1/µ defines a characteristic length scale for the transition between the SR and LR descriptions. In this scheme as µ → ∞ the RSH incorporates full HF exchange, while µ → 0 leads to a usual global hybrid (with α). In this way the RSH is able to smoothly cover pure LDA/GGA exchange and full HF exchange. Fig. 4 illustrates the three main classes of hybrid DFT functionals.

Fig. 4 The fraction of the exact HF exchange used in the calculations as a function of the separation of the orbitals. Three different classes of hybrid functionals are shown. From left to right: global hybrid functionals, RSH functionals with zero SR HF exchange and RSH functionals with nonzero SR HF exchange.

By treating LR interactions with (localizing) full HF exchange and SR interactions by (delocalizing) semi-local DFT exchange, LRC hybrids intrinsically balance the spurious localization error of HF and the spurious delocalization error of semi-local DFT.⁴⁸ The functional behaves “electron gas-like” in the short range, but has the correct long-range and asymptotic behavior of wavefunction based methods.

The LC-PBEOP, LC-BLYP and LC-BOP were one of the first RSH functionals to be developed in 2004. They already provided much better results for optical excitations and CT excitations.⁵¹ For simple molecules the LC-BOP functional provides good ionization potentials and electron affinities from its frontier orbital energies.⁵⁷ It has been claimed that there is clear-cut proof that the correct description of both Rydberg and CT states can be restored using RSHs.²⁷ The other successful RSH functional was the rCAM-B3LYP, which was developed in 2007 to overcome the SIE in many-electron systems.⁴⁶ Several groups have shown that for a well-constructed RSH functional, the curves of the energy as a function of the fractional number of electrons are much more linear than those obtained with conventional functionals.^29,52,53

Bartlett et al. coined the term ab initio DFT,³¹ which is based on the principles of Correlated Orbital Theory (COT). These principles were established already in 2005,⁵⁴ but the first practical COT DFT functional – CAM-QTP-00 – was introduced in 2014. It was inspired by the good performance of the CAM-B3LYP functional.⁵⁵ The CAM-QTP-02 and QTP17 were introduced in 2018.^47,56

The key idea here is that the Kohn–Sham system is fictitious as the KS orbitals themselves can be viewed as mere mathematical devices en route to obtaining the correct ground-state density. By putting additional constraints on theses KS orbitals, we arrive at the generalized KS system. In practice, the GKS contains some amount of exact HF exchange or other non-local component such as the meta-GGA kinetic energy gradients. In principle it can contain other orbital-dependent potentials.

Bartlett's ab initio DFT is GKS scheme where the DFT functionals are not be empirically fitted to the problem at hand (molecular geometries, polarizabilities, excitation energies, thermochemical properties etc.) but rather must reproduce the known properties of the correct orbitals and density.⁴⁷

It is known that the asymptotic decay of the density (r → ∞) should only involve the ionization potential of the many-electron system. The Koopmans theorem of HF theory shows that this is the HOMO energy. We would like the same behavior in ab initio DFT. This fact, that the RSH functionals can be used to obtain the correct long-range asymptotic behavior, was also shown in 2005.⁵⁷

In effect, the DFT version of Koopmans' theorem – Janak's theorem⁵⁸ – must hold. In principle this condition is imposed for both the neutral and the anionic state, where the ionization potential corresponds to the electron affinity of the neutral state:

Sometimes deeper levels will have to be considered as well and a dication calculation needs to be performed:⁵²

E_tot(M_N−2) − E_tot(M_N−1) = −ε_HOMO−1

Besides the correct asymptotic decay of electron density, the goal is that the determinant that emerges from COT, will have the property that all its orbitals will correspond to exact principal ionization energies (eqn (2) and (3)). This is in contrast to the original KS idea where the ground state density was the objective (eqn (5)). The COT believes that the ‘exactness’ of the one-particle energy spectrum is actually far more important than an exact density from a single determinant.³⁶ This is easily confirmed by this work, as we are interested in the energies of the frontier orbitals rather than the exact ground state electron density. Furthermore, there is strong evidence that the SIE is mitigated in functionals where the orbital energies reflect the vertical ionization energies.^47,48,56

The COT MOs reflect observables, like ionization potentials and electron affinities, instead of some nebulous approximations.⁴⁷ This also implies that the need for the detailed treatment of electron–electron repulsion is removed as whatever effects introduced by it are now hidden in the orbital energies themselves.⁴⁷ The Bartlett's ab initio DFT scheme can also be considered to be an extension of the HF theory to include correlation.³⁶ Additionally, this CGS scheme allows us to retain the single particle picture as these GKS orbitals are now quasi-electrons or quasi-holes, i.e. effective particles that contain (“are dressed by”) the relaxation effects of the other electrons.²⁹ As an electron is inserted or ejected (leaving a hole behind), all other electrons in the system respond to the presence of this extra electron or hole. The ionization potential and electron affinity are now the lowest quasi-hole and quasi-electron excitation energies, respectively. The HOMO–LUMO gap should mimic the quasiparticle gap and the energies of the GKS orbitals should be similar to the quasiparticle energies. However, this should be achieved without using a perturbative calculation that includes excited states.

For example, the first order G₀W₀ correction to the GKS orbital energies would be:

So, if the DFT exchange–correlation energy is roughly equal to the G₀W₀ self-energy Σ₀, then the orbital energies would not change and in principle the GKS orbitals would approximate the G₀W₀ quasiparticles. Thus, modern DFT has made full turnaround: the KS orbitals that were initially formally introduced to calculate the kinetic energy functional are now expected to approximate the ionization Dyson orbitals of eqn (3) and (4). The GKS with COT principles will be considered the final (or latest) development of modern DFT.

Since the first implementation of the RSH method, several schemes were proposed to determine the range separation parameter µ. Some of them focused on an empirical parameterization to minimize errors on general or specific energy properties.⁵⁹ They led to values for µ lying typically between 0.2 and 0.5 bohr⁻¹, thus covering characteristic separation lengths between 2 and 5 bohrs.

Starting from very simple molecules such as the diatomic Li₂, S₂, P₂, O₂, CO and HF it was found that the fundamental gap is a very sensitive function of the range parameter µ. For example, for F and O the gap changes by as much as 8 eV when µ changes from 0.3 to 1.²⁹ The optimal µ varies from 0.3 bohr⁻¹ in Li₂ to about 0.7 bohr⁻¹ in F₂.⁶⁰ This phenomenon has been attributed to delocalization trends: as orbitals delocalize with increasing system size, semi-local approximations become more accurate, and the spatial range, 1/µ, at which exact-exchange–corrections are needed, becomes larger.²⁹

It has been demonstrated that the optimal µ value can vary significantly with system size. A similar picture is seen for other LRC hybrids, independently of the form of the underlying semi-local functionals and the amount of HF exchange used in the SR (α).⁴⁸

When considering the optimal value of µ for the LC-ωPBE functional vs. the number of repeat units of four different molecular chains (polyenes, alkanes, oligoacenes, and oligothiophenes) it was found that in all instances, the range-separation parameter decreases significantly with chain length.⁴⁸ This clearly demonstrates that not only the size but also the extent of conjugation plays an important role in determining the optimal µ values for a specific system of interest.⁴⁸

Thus, generally speaking, the standard µ-values in LRC hybrids are not suitable for a reliable description of a ‘new’ system. Furthermore, there is no universal, system-independent value for the range-separation parameter.⁴⁸ It is already well established that it is difficult to optimize an LRC or CAM functional to reproduce both ground and excited state properties.³² Also, simultaneous good performance on small and large molecules is difficult. A rigorous analysis, based on the adiabatic connection theorem, shows that µ is itself a functional of the electron density.^57,60

So what RSH µ-value should be used for ionic liquids? The answer lies in the GKS COT principles we discussed previously. The RSH µ value will have to be found self-consistently to satisfy eqn (9). The ionization energy will be calculated with the ΔSCF method. This idea is captured by the following equation:

E_ioniz[µ] ≈ E_tot(M_N−1) − E_tot(M_N) = ΔSCF = −ε_HOMO[µ] (10)

The value of µ that satisfies eqn (10) is called optimal and the procedure of finding this value is called optimal tuning (OT). The RSH functional with the optimal µ value for the system of interest is called optimally tuned or IP-optimized and it is denoted with an asterisk.

The tuning idea can be generalized, for example the first triplet excited state can be tuned for.⁶¹ Fig. 5 illustrates the [EMIM][OTF] OT procedure for several functionals.

Fig. 5 The ionization energy error of several RSH functionals as a function of the range-separation parameter µ. Optimal values of µ are highlighted with arrows.

This self-consistent OT approach allows to overcome the difficulty of the absence of a universal, system-independent value for a range-separation parameter. Instead it gives a natural, physically motivated procedure for finding a system-dependent range-separation parameter µ. OT does not introduce any empiricism to the underlying functional, as the calculated ΔSCF and HOMO are not matched to any experimental value.⁶¹ Instead, the range-separation parameter µ is optimized to achieve an internal consistency.

The OT approach initially demonstrated its success with long-range through-space CT excitations occurring in molecular dimers and has largely proved its robustness through numerous ground and excited state applications especially prone to one- and many-electron SIEs.⁶² The OT functionals can produce very similar results to GW calculations, for example for the valence band of aromatic molecules.⁵² Optimally tuned RSH also significantly remedies the cation dimer problem, improves the CT and Rydberg excitation energies.⁶³

However, in some cases the charge localization can occur too abruptly and at too short distances. In other words, the OT RSH functionals are still not perfect in the long and intermediate range. This can be explained by the fact that RSH functionals are unable to describe static correlation in the LR⁴⁸ and this error increases as the HF exchange content in the functionals is increased.⁶⁴

Still, the localization-delocalization error at the HOMO level is significantly reduced in the OT RSH functionals compared to HF and all commonly used semi-local, global hybrids and LRC-hybrid functionals. In fact, OT can be interpreted as a minimization of the delocalization error at the level of the HOMO.^48,63

5. Results and discussion

The reader is now referred to Table 3, where the calculation results are presented. Four main classes of DFT functionals were evaluated: non-hybrid (black-colored), global hybrids (brown-colored), range-separated hybrids (blue-colored) and meta-GGA functionals (magenta-colored). The latter can be either non-hybrids, global hybrids or RSHs. They can be distinguished by looking at the second column to see the amount of exact HF exchange they contain, if any.

Most commonly, eqn (9) does not hold and the HOMO energies may not correspond to the vertical ionization energy. The HOMO or the HOMO–LUMO gap may correspond to the optical gap instead. For example, it has been repeatedly observed that HOMO–LUMO gaps derived from the HSE06 are often close to the optical gap, rather than the fundamental one.²⁹ In this study, the only functional that exhibits signs of this behavior is CAMH-B3LYP. This functional yields a HOMO energy of −7.9 eV, which roughly aligns with the main PIY threshold energy. It also produces a HOMO–LUMO gap of 6.9 eV, which approximately matches the CT excitation energy. Finally, its ΔSCF energy of 9.1 eV matches the ionization energy. Further investigation is needed to determine whether this alignment is coincidental or indicative of a more systematic trend.

The HF calculation itself (not shown) does not provide a reliable description of the [EMIM][OTF] electronic structure. To qualitatively match the experimental UPS spectrum, a typical energy scale compression of 15–20% and a shift of approximately 1.0 eV are required. This is typical for HF/MP2.^6,15,65 Moreover, both the HF and MP2 ΔSCF calculations overestimate the vertical ionization energy by approximately 1.6 eV. The RPA TD-HF calculation also overestimates the main CT peak by roughly 0.4 eV, while the BSE@HF calculation overestimates the experimental CT energies by about 0.6 eV. Using the larger aug-cc-pVTZ basis set further increases the energy scale shift to 1.0 eV. However, after correcting for this large energy shift, the BSE-calculated cross sections closely resemble those predicted by the simple ‘3 + 2 model’. Additionally, the BSE@HF is the only method that reproduces the expected A_1,2 intensity ratio 1 : 2 without introducing spurious excitations, and it clearly resolves the B₁ peak (see Fig. 6).

Fig. 6 From top to bottom: the topmost green curve labelled BSE@LC-BOP refers to the Bethe–Salpeter equation calculation from the LC-BOP functional using the extended aug-cc-pVTZ basis set. The experimental PIY spectrum showing the intensity of the EMIM cation is labelled PIY (blue points). The red curve labelled BSE@HF refers to the Bethe–Salpeter response calculation from Hartree–Fock (G₀W₀) orbitals, also in the extended basis set. The arrow indicates the direction and the label next to it indicates the magnitude of the shift (1.0 eV) used to align this calculated BSE spectrum with the experimental PIY spectrum. The black line is calculated by modelling transitions in a simple 3 + 2 model in which 3 occupied states get excited to 2 unoccupied states. The individual transitions with their magnitudes are shown with arrows below.

5.1 GW results

The GW calculation improves the HF HOMO energy by 1.1 eV, resulting in a G₀W₀@HF HOMO value of −9.5 eV.

Our results indicate that G₀W₀ and fully self-consistent GW (scGW) yield very similar outcomes for most hybrid functionals tested, suggesting that GW values are nearly converged with respect to both the number of iterations (G_nW_n) and the method applied. However, in the case of HF, scGW@HF provides better agreement with experiment, improving the HOMO energy to −9.2 eV. For certain functionals, such as LCY-PBE, the G₀W₀ result of −8.8 eV is actually closer to the experimental value than the fully self-consistent scGW result of −8.5 eV. Indeed, there are cases where G₀W₀ offers better agreement with experiment, while in other instances, quasiparticle self-consistent GW (qsGW) performs better.⁶⁶ In principle, it is difficult to predict a priori which GW variant will yield more accurate results.

In addition to differences between one-shot G₀W₀, scGW, and quasiparticle self-consistent methods (such as qsGW₀, which is self-consistent only in the Green's function), basis set dependence introduces further complications.^66,67 All G₀W₀ and BSE values presented in Table 3 were obtained from calculations using the aug-cc-pVDZ basis set. When upgrading to the aug-cc-pVTZ basis, all G₀W₀ HOMO energies decrease by approximately 0.5 eV. Complete basis set (CBS) extrapolation, which is available in MolGW,⁶⁸ reveals that the total decrease in the HOMO energy is about 0.7 eV. This systematic shift arises because both the Green's function (G₀) and the self-energy (Σ₀) depend on the excited states used in the calculation of the screened interaction (W₀) and the susceptibility (χ₀), both of which are sensitive to the choice of basis set. The extra energy originates from the correlation part of the self-energy.

Therefore, in the CBS limit, all G₀W₀ calculations that have the correct −9.0 eV HOMO values in Table 3 would actually underestimate the HOMO energy and have −9.7 eV listed instead. And contrary, for example, in the CBS limit the G₀W₀@PBE0 value of −8.4 eV would indeed improve to the correct −9.1 eV value, in agreement with previous results.² The abovementioned scGW@LCY-PBE value of −8.5 eV would similarly improve to −9.2 eV.

The energy variations resulting from the choice of functional, basis set, and GW approach are too large to allow reliable predictions without extensive benchmarking. Differences between atom-centered and plane-wave basis set implementations further compound these challenges.⁶⁷

Based on the above discussion, although the G₀W₀@aug-cc-pVDZ offers indeed excellent results for the HOMO energy of [EMIM][OTF], it is difficult to predict whether this combination will also work on other ILs.

5.2 Non-hybrid DFT

The HLE16 and HLE17 functionals were modeled specifically for bandgaps and molecular excitation energies. They perform indeed better than the HCTH-A and HCTH-407+ functionals, which were optimized for HOMO energies.

The non-hybrid meta-GGA functionals TPSS, revM06-L, M06-L, M11-L, and SCAN also perform poorly for our [EMIM][OTF] system, similarly underestimating the absorption threshold substantially and exhibiting a large ΔSCF ionization energy vs. HOMO energy discrepancy (≈3 eV).

Notably, the HCTH-407+ seems flawed as it predicts a HOMO–LUMO gap of 3.1 eV and optical TD-DFT excitations at 3.2 eV, 4.3 eV and 5.2 eV. Some researchers consider these predicted spurious CT excitations at low energies (3.2–3.5 eV, in this case) to be a catastrophic problem.^27,69 While GW and BSE calculations improve upon these values, substantial discrepancies remain. Overall, the poor performance of this class of non-hybrid functionals underscores the necessity of incorporating Hartree–Fock exchange for accurate predictions.

5.3 Meta-GGA functionals

The meta-GGA functionals depend on the kinetic energy density and thus on the occupied KS orbitals. The orbitals are implicit, semi-local functionals of the density and therefore, in principle, can improve band-gap prediction when evaluated in a GKS scheme. However, it is hard to produce v_xc derivative discontinuities in the meta-GGA scheme that would increase the energy gaps.⁷⁰ For this reason, the non-hybrid meta-GGA functionals still underestimate the HOMO energy substantially.

The SCAN functional is widely used for band-gap estimations. SCAN stands for strongly constrained and appropriately normed and this functional was constructed by imposing 17 known exact constraints. The SCAN functional shows numerical instability due to unphysical divergence of the exchange–correlation potential. Subsequently, the R2SCAN functional was developed, which results in a functional that not only satisfies all constraints but also is numerically stable.³⁴ The TPSS functional is the other nonempirical meta-GGA functional and it was developed to satisfy the correct exchange-hole properties i.e. to correct for the SIE.⁴⁹ Indeed, these two functionals predict LUMO energies close to the proposed LUMO energy of −2.3 eV. Even more usefully, the ionization energies of the hybrid versions of these meta-GGA functionals SCAN0 and R2SCAN0, estimated using the ΔSCF method, are in agreement with the experimental HOMO energy.

However, for the non-hybrid meta-GGA functionals like (R2)SCAN, M06-L, M11-L, the ionization energies from the ΔSCF and Janak's theorem differ by up to 3 eV, as mentioned before. Furthermore, the same M11-L, SCAN and TPSS provide gaps of about 3.4–3.8 eV, which are totally inaccurate. Therefore, in practice, the band gaps from the non-hybrid meta-GGAs are not quantitatively accurate for ILs.

Increasing the HF fraction to α = 0.5 as in R2SCAN50 (ref. 71) reduces the HOMO error to just 0.2 eV compared to experimental photoemission data. The other meta-GGA functional that already predicts the correct HOMO energy (see Fig. 3) is the M08-SO, which contains the ‘magic’ 56% HF exchange. Heavily parameterized and optimized, M08-SO ranks first and second among 200 tested functionals in benchmarks for non-covalent dimer interactions and barrier heights, respectively.⁷² It also performs well for the LUMO energy (see section below).

In contrast, other Minnesota functionals do not achieve the same level of accuracy. For example, the new global hybrid M15 and the previous RSHs M11 and revM11 do not offer good agreement with experimental values. The M06-2X, which also includes a 54% HF exchange, shows relatively good performance in several studies,²⁷ but still exhibits a ∼0.5 eV deviation in ground-state orbital energies for the [EMIM][OTF] ion pair.

5.4 Global hybrids

Note that, in the HSE06 functional, unlike in the conventional RSH scheme, HF exchange is used in the SR only, with no LR HF exchange. It's range-separation parameters are chosen to maximize computational expedience and optimize metallic properties rather than to correct the SIE, making it more like a GGA functional when it comes to understanding transitions involving localized defect orbitals and conduction/valence band orbitals.⁷³

It is evident that standard global hybrid functionals such as B3LYP, HSE06, LDA0, and PBE0 generally fail to predict accurate values for experimental observables. These functionals tend to produce spurious CT transitions around 4.5 eV and predict band gaps in the range of 5–6 eV—significantly deviating from experimental results. While functionals like PBE1/3 and MPW1K perform slightly better, they still yield unphysical CT excitations.

As shown in Table 3, the inclusion of a substantial fraction of exact HF exchange is crucial for obtaining reliable results. Previous studies have consistently identified 50% exact exchange as a key threshold for improving accuracy. Among global hybrids, the best-performing functionals typically feature HF exchange fractions between 41% and 54%.²⁷ For example, in a study of CT excitations, BHHLYP (which includes 50% HF exchange) was the only functional that delivered reasonable excitation energies.²⁷ Moreover, intermolecular parasitic CT states are significantly suppressed when α ≥ 0.5.²⁷

In other words, the HF exchange fraction is the primary factor in determining accurate gap values. While the elimination of long-range exchange or the semi-empirical tuning of exchange–correlation components can contribute, they play a secondary role.²⁹

Looking specifically at PBE1/3 (33% HF), which still produces spurious CT states, it is clear that 33% HF exchange is insufficient. Even PBE50, with 50% HF, shows a weak CT excitation at 6.53 eV, suggesting that in some cases, even 50% may not be enough.

A study on small oligoacenes found that increasing the HF content in PBE0 (α = 0.25) to α = 0.75 improved gap predictions and hole localization in molecular solids.²⁹ Similar behavior was described by Laurent et al. for other small molecules and other underlying functionals.²⁷ BH&HLYP, which incorporates 50% HF exchange, has been found among the most successful global hybrid functional starting points for G₀W₀.⁷⁴ Later Bruneval et al. showed that for G₀W₀ calculations on the GW100 set the global hybrid PBEh with an increased α = 0.75 showed the best systematic agreement with CCSD(T) values.⁴³ The BSE calculation with 100% HF exchange shows no spurious CT peaks and has the intensities for the A_1,2 CT peaks close to the expected 1 : 2 ratio (see Fig. 6), providing further evidence for the benefits of high HF exchange in IL calculations.

However, it is important to note that increasing the HF exchange fraction to such high levels is known to compromise the accuracy of thermochemical and ground-state electronic structure predictions.²⁹ As Kronik et al. pointed out, pushing α beyond 0.5 may break compatibility with the semi-local correlation component.⁵² Furthermore, the full HF exchange methods, such as MP2 and M06-HF already produce too high binding energies and overbind the ion-pair. Similarly, the aforementioned BSE@HF calculation already substantially overestimates CT excitation energies.

Examining the trend across functionals—PBE0, PBE1/3, MPW1K, PBE50, as well as meta-GGAs like SCAN, R2SCAN0, and R2SCAN50—it appears that an HF exchange fraction in the range of 50–60% is close to optimal. Bartlett also mixed in 50% HF, but called this choice arbitrary.³⁶ Jiang et al. also noticed that a suitable amount of exact exchange, in the range of 44–56%, is in favor of producing neither too delocalized nor too localized electronic structures.⁷⁵

That said, as Holzer et al. emphasized, in typical cases such as valence excitations in organic molecules, where CT states are absent, high HF exchange is not necessarily beneficial. In fact, it can increase mean errors and worsen the performance of functionals when compared to their low or no-HF exchange counterparts.⁷⁴

5.5 RSH functionals

Bartlett's QTP functionals indeed do avoid spurious CT excitations predicted by many previous generation DFT functionals such as LDA0, PBE0, B3LYP, (TUNED-)CAM-B3LYP, HSE06 etc. Among them, CAM-QTP-01 and CAM-QTP-02 show a slight shift of 0.15 eV in the main PIY absorption peak, predicting it at 8.15 eV instead of the experimental 8.3 eV. However, CAM-QTP-00 fails to capture the lower-energy shoulder of the main peak and does not predict the B₁ absorption feature at 7.9 eV. The global hybrid QTP17 is the best among them not only in TD-DFT excitation energies but also it has the smallest overall absolute error in the description of the experimental UPS. However, the QTP17 severely overestimates the intensity of the A_1,2 doublet compared to the expected intensity from the ground state occupations.

Weintraub et al. proposed an alternative improvement to CAM-type functionals by incorporating accurate exchange hole properties.⁷⁶ The Weintraub's LC-ωPBEh-whs functional performs the best out of the tree tested. The “h” indicates the presence of some short-range HF exchange i.e. a nonzero α value and for this reason it predicts neutral excitations in its TD-DFT calculation. One can argue that the Weintraub's LC-ωPBEh-whs performs comparably to the ωb97X family.

RSH functionals with α = 0, like LC-BLYP, LC-BLYP-EA or LC-wPBE are able to describe the ground state and charged excitations, but require a BSE calculation to describe the neutral excitations. Out of the α = 0 RSH functionals the LC-BLYP-EA stands out. The experimental UPS spectrum agrees well with its ground state DOS (Fig. 3) and the TD-DFT calculation provides ionization cross sections (see Fig. 7). Nevertheless, the α = 0 RSH functionals still struggle to simultaneously provide accurate descriptions of both charged and neutral excitations, highlighting the limitations of fully long-range-corrected functionals without SR HF exchange.

Fig. 7 The experimental PIY spectrum from this work showing the intensity of the EMIM cation is labelled PIY (blue points). The photoabsorption profile calculated using the PBE-2X functional is labelled TD-DFT PBE-2X (green curve). The gray curve labelled BSE refers to the Bethe–Salpeter response calculation from the LC-BOP functional. The arrow indicates the direction and the label next to it indicates the magnitude of the shift (0.25 eV) used to align this calculated BSE spectrum with the experimental PIY spectrum. The photoabsorption profile calculated using the optimally tuned CAM-PBE50* functional is shown with a magenta curve below. The black line is calculated by modelling transitions in a simple 3 + 2 model in which 3 occupied states get excited to 2 unoccupied states. The individual transitions with their magnitudes are shown with arrows below. The bottom blue curve is from the TD-DFT calculation for the LC-BLYP-EA functional. It shows charged excitations and ionization (peaks I_a, I_b, I_c). The experimental UPS spectrum (black points, labelled Exp. UPS) is reproduced from ref. 2 and 14 and the energy scale of the experimental UPS spectrum is shown on the top axis.

5.6 Optimal tuning

Interestingly, the LC-BLYP default value for the range separation parameter of µ = 0.33a₀⁻¹ is the exact optimal value for [EMIM][OTF] (Fig. 5). The OT LC-BLYP* functional offers an excellent description of the ground state DOS (Fig. 3). LC-BLYP also performs very well in the fractional occupation test.⁴⁶

The CAM-LDA0 functional was evaluated for OT. Since eqn (10) could not be satisfied for this functional, the default range-separation parameter µ = 0.33a₀⁻¹ was used.⁷⁷ This yielded a TD-DFT excitation spectrum in good agreement with experiment, with an overall energy shift of only 0–0.2 eV. The ground-state orbitals also matched the experimental UPS spectrum exceptionally well, showing no energy shift (see Fig. 3). Notably, CAM-LDA0 is one of the few functionals that does not predict a near-zero or positive-energy LUMO. It is also computationally faster than many GGA-based functionals. However, it does predict a spurious CT excitation at 5.7 eV.

Several functionals using µ = 0.3–0.33a₀⁻¹—such as LC-BLYP, LC-PBEOP, LRC-ωPBE, LC-BLYP-EA, CAM-LDA0, and CAMY-B3LYP—exhibit small energy scale errors of just 0–0.2 eV compared to experiment. As previously discussed, the OT procedure for LC-BLYP also converged to µ = 0.33a₀⁻¹ (see Fig. 5), suggesting that this value is near-optimal for functionals with α = 0 or low HF exchange content, such as CAMY-B3LYP. As the µ value is increased, the energy scale errors of the α = 0 RSH functionals also increase. For example, the µ = 0.4a₀⁻¹ functional LC-wPBE and the µ = 0.45a₀⁻¹ functional LC-VV10 overestimate the orbital binding energies. At even larger values, such as µ = 0.75a₀⁻¹ in LCY-BLYP, the errors become severe, approaching a 4 eV overestimation. In contrast, functionals with higher HF exchange content require smaller µ values for optimal performance (see Fig. 5).

The importance of optimal µ tuning is illustrated by the LRC-PBE0 (also known as CAM-PBE0 (ref. 78) or RSX-PBE0 (ref. 62)) functional. The commonly used default µ values—0.35a₀⁻¹ or 0.47a₀⁻¹ (ref. 69) or 0.35–0.39a₀⁻¹ (ref. 62)—are likely too large and may contribute to its poor performance in thermochemical calculations.⁷⁸ Rohrdanz et al. also suggested that a much smaller value of µ is appropriate for TD-DFT.³² In our system, the OT procedure yields µ = 0.23a₀⁻¹ for LRC-PBE0, which improves the predicted energies by more than 1.5 eV. This demonstrates the significant performance gains achievable through OT. Indeed, the choice of the range-separation parameter µ is a decisive factor in the accuracy of RSH functionals.⁴⁸ Even so, the LRC-PBE0* still has a single weak spurious excitation at 6.36 eV.

Increasing the HF exchange further leads to the CAM-PBE50 functional, also known as RSX-0DH(SCF).⁶² Now the spurious CT has almost disappeared being just off 0.03 eV from the 6.8 eV CT state. When the CAM-PBE50 is tuned according to the second part of eqn (9) i.e. trying to optimize for the electron affinity, then an even smaller value of µ = 0.09 is found to be optimal. However, this still gives a very similar agreement with the UPS and PIY spectra. Interestingly, the CAM-PBE50* functional seems to predict the stronger triply degenerate HOMO−1 as the HOMO. This can be seen from its predicted ionization energy of 10.25 eV. However, in reality this is due to an energy scale shift of 1.0 eV, possibly caused by the high HF exchange content and resulting overbinding. The small µ and reduced long-range HF exchange cause CAM-PBE50* to behave similarly to its global hybrid counterpart, PBE50. Indeed, the TD-DFT photoabsorption profiles of CAM-PBE50* and PBE50 are similar. As previously discussed, further increases in α are expected to worsen overbinding and reduce accuracy.

It's worth noting that while OT α = 0 functionals such as LC-PBEOP, LC-BLYP*, and LC-BLYP-EA perform excellently for [EMIM][OTF], our preliminary results indicate that this performance does not generalize across other ILs. As such, these functionals should always be optimally tuned prior to application in new systems.

The overall recommendations of the evaluated calculation methods are shown in Table 1.

Table 1 Recommendations of studied functionals for use in similar systems

Rank	Ground state and ionization threshold description	TD-DFT (neutral) excitations
1	PBE-2X	PBE-2X
	M08-SO	QTP17
2	Optimally tuned α = 0 RSHs:
	LC-BLYP-EA*	CAM-PBE50*
	LC-PBEOP*	ωB97X
	LC-BLYP*	ωB97X-V
3	CAM-LDA0	LC-ωPBEh-whs

---

5.7 Comparison with literature benchmarks

There are no prior comprehensive comparison studies of IL calculations. In 2013 Zahn et al. evaluated about 25 functionals.⁷⁹ The LC-BOP and M06 were among the top performers.

However, there are numerous vertical ionization benchmark studies. The GW100 set is a popular reference set containing small organic and inorganic molecules. In case of the ionization energies from the GW100 and GW27 datasets, calculated in the G₀W₀ approximation, the CAM-B3LYP and ωB97X-D performed average, while CAM-QTP-00 and CAM-QTP-02 were somewhat worse.⁷⁴ The set also showed that the (OT) RSH functionals were an excellent starting point for the G₀W₀ calculations.^74,80 In the study by Loos et al. it was found that the very popular B3LYP and PBE0 tend to underestimate vertical transition energies with errors on the order of 1.0 eV.⁸¹ Interestingly, this is very close to the error of these two functionals for the [EMIM][OTF] case.

The Caricato set is composed of 30 valence and 39 Rydberg experimental excitation energies of 11 small molecules. The CAM-B3LYP, M06-2X and ωb97X-D showed the best performance in that set.²⁷ Besides the benzoic acid set the CAM-B3LYP and M06-2X came out on top in another small molecule benchmark.²⁷ In a very recent study by Holzer and Franzke ωB97X-D was identified as the top-performing functional for TD-DFT calculations of CT excitations.⁸² Additionally, Loos et al. reported that CAM-B3LYP, M06-2X, and ωB97X-D exhibited negligible mean errors in TD-DFT calculations of CT and strong CT transitions.⁸¹ These three functionals outperformed not only other DFT functionals such as M11 and ωB97X, but also several wavefunction-based approaches including EOM-MP2, CIS(D), (EOM-)CCSD, ADC(2), ADC(3), and even CC2.

Indeed, the CAM-B3LYP and M06-2X seem to predict very similar energies for our [EMIM][OTF] system. Moreover, CAM-B3LYP and ωB97X-D yield comparable TD-DFT photoabsorption profiles. However, in this work, the ωB97X functional performs better in describing CT excitations than ωB97X-D. Notably, none of the three, CAM-B3LYP, ωB97X-D, or M06-2X, offer the best overall agreement with the [EMIM][OTF] experimental data. This is likely due to the more demanding HF exchange requirements of the IL ion-pair system. As a result, functionals with higher exact exchange content, such as QTP17 and PBE-2X, outperform the aforementioned functionals.

In a recent study, Patra et al.³⁴ investigated a range of organic chromophores and identified the ωB97 family, M11, and ωM06-D3 as the top-performing functionals. In contrast, CAM-B3LYP, SCAN0, and M06-2X showed noticeably lower accuracy, while PBE0, TPSSh, SCAN, and M06-L performed significantly worse. The latter finding is consistent with this work.

Our preliminary analysis suggests that, for some functionals the average of the ΔSCF and −HOMO energy provides a more reliable estimate of ionization energy, particularly across different ILs. For instance, the M06-2X functional yields a value of 9.05 eV, which matches perfectly with the [EMIM][OTF] case. CAM-B3LYP and R2SCAN50 also seem to show similar behavior.

5.8 LUMO

It is known that in DFT the fulfillment of Janak's theorem and simultaneously satisfying a similar relation between the LUMO and electron affinity, is very difficult.²⁹ Furthermore, DFT struggles to accurately describe unoccupied (virtual) states without the aid of TD-DFT or BSE calculations. This issue is similar with the HF virtual orbitals, which are positive for neutral systems.³¹ Bartlett proposed that incorporating an appropriate fraction of exact Hartree–Fock (HF) exchange may, in principle, provide a reasonable zeroth-order approximation to the orbital energy spectrum,³¹ consistent with the previously discussed “3 + 2 model”. Within this framework, the calculated LUMO energy should correspond to the energy of the first free valence-band (CT) state, NOT to the electron affinity E_a defined by eqn (11):

E_a(M) = E(M) − E(M⁻) (11)

The experimental electron affinity of [EMIM][OTF] is unknown, but based on similar molecules such as pyrrole, imidazole⁸³ and especially the MMIM cation⁸⁴ it is likely negative or close to zero.

As the LUMO energies are widely used in the literature,^85,86 it remains important to compare the energies of the virtual states between DFT functionals. Thus, based on the previous discussion, the calculated LUMO energy of [EMIM][OTF] is expected to approximate the energy of the CT state. Indeed, for many functionals, a π-symmetry orbital on the EMIM cation (Fig. 1) or similar appears as the LUMO.

Table 2 presents calculated electron affinities, LUMO and higher unoccupied orbital energies for several functionals. Many RSH functionals were omitted, as they predict a close to zero energy LUMO similar to the RSH functionals included in Table 2.

Table 2 Comparison of calculated electron affinities, LUMO energies and the next unoccupied virtual states for selected functionals. The ordering of the methods is subjective and reflects how well each approach reproduces (i) strong LUMO binding and (ii) the expected intensity ratio between the proposed LUMO and the (degenerate) LUMO+1 orbitals, which should follow a 1 : 2 ratio based on the derived frontier orbital model

Method	ΔSCF electron affinity eqn (11) (eV)	LUMO (eV)	Next LUMO in DOS (eV)	Intensity ratio in DOS
M08-HX	0.3	−1.4	−1.0	1 : 2
M08-SO	0.25	−1.3	−0.9	1 : 2
B3LYP	0.2	−1.4	−0.9	1 : 2
LDA0 and PBE0	0.1–0.2	−1.25	−0.82	1 : 2
M11-L	0.4	−2.65	−2.3	1 : 2 : 2
M06	−0.1	−1.8	−1.45	1 : 3
TPSS	0.1	−2.1	−1.7	1 : 1
SCAN and R2SCAN	0.1	−2.0	−1.5	1 : 1
M06-L	−0.2	−1.9	−1.45	1 : 1
M11	0.0	−0.9	−0.5	1 : 2
PBE1/3	0.1	−0.85	−0.45	1 : 2
revM06	−0.15	−0.6	−0.2	1 : 2
LRC-PBE0*	1.8	−1.4	−0.8	1 : 2
M06-2X	−0.05	−0.8	−0.4	1 : 2 : 2
CAM-LDA0	0.7	−1.6	−1.1	1 : 1
TPSS0	NA	−1.05	−0.55	1 : 1
revM11	0.1	−0.1	0.3	1 : 2
LC-BLYP*	−0.1	0.07	0.45	1 : 2
LC-BLYP-EA	−0.1	0.07	0.42	1 : 2
LC-PBEOP	−0.15	0.25	0.65	1 : 2
PBE50	0.0	−0.25	0.15	2 : 3
CAM-PBE50*	0.6	−0.35	0.15	2 : 3
PBE-2X	0.0	−0.25	0.2	2 : 3
ωB97X-V	−0.1	0.3	0.65	1 : 2

---

Table 3 Comparison of different calculation methods. All energies are in electronvolts (eV). Global hybrid functionals have brown labels, blue labels are used for RSH functionals and meta-GGA functionals are colored magenta. The number after the functional in parenthesis indicates the corresponding functional number in the LIBXC library. Green highlighted values show excellent agreement with experimental values. The second column indicates the convergence of the ΔSCF shown in the fourth column. Y means that the ΔSCF is stable under basis set change to aug-cc-pVTZ and/or geometry optimization, N means that the ΔSCF energy changes more than 0.2 eV. A question mark indicates that it was not tested

a RPA TD-HF calculation.

---

Functionals closest to the expected behavior are listed at the top of Table 2. The focus here is on whether these methods can accurately predict both the energy separation and intensity ratio among the LUMOs.

Among the tested functionals, M08 predicts the proposed LUMO : LUMO+1 intensity ratio of 1 : 2 and the expected energy separation of 0.4 eV. The M08 functionals are followed by the widely used B3LYP functional; indeed, B3LYP also performs very well for calculating electron affinities of simple molecules.⁸³

In contrast, high-HF-exchange non-meta-GGA functionals—such as LC-BLYP*, PBE-2X, and members of the ωB97X family, though excellent for describing occupied ground-state orbitals, fail to provide meaningful LUMO energies. Increasing HF exchange causes the virtual orbitals to shift into positive energy, similar to what is observed in pure HF. For instance, while PBE0 and PBE1/3 still produce negative LUMOs energy, PBE50 already predicts them as positive. PBE-2X and the ωB97 family, which are among the most reliable functionals for predicting ground-state orbitals, now appear at the bottom of Table 2.

This deficiency is less severe in meta-GGA functionals, which appear to offer a more accurate quasiparticle-like description of the virtual orbital manifold than conventional hybrids or RSH functionals. In meta-GGAs, increasing HF exchange shifts the occupied orbitals more toward negative energies, while the virtual orbitals do not shift toward positive energies to the same extent.

The fact that almost all ab initio methods overestimate the proposed LUMO energy may explain why, for example, the electrochemical stability windows are ‘calculated’ by combining results from the isolated cations and anions calculations.^3,87

6. Conclusions

It is important to emphasize that the present conclusions are subject to several limitations. Only a single conformer was considered in the calculations. The basis set employed was predominantly aug-cc-pVDZ, with limited comparisons to aug-cc-pVTZ and CBS extrapolations. In addition, the presented experimental PIY data are relatively coarse in resolution. Finally, any extension of these findings to other ionic liquids remains speculative.

The IL conformer present at the evaporation temperature of 180 °C was re-evaluated compared to previous studies, guided by the additional constraints provided by the absorption spectra. Indeed, the prediction of the correct IL conformer structure was shown to be highly nontrivial. The charge transfer process in the [EMIM][OTF] can be understood as an electron transfer from the anion to the cation. The observed CT energy is in good agreement with the estimated electrostatic stabilization energy of the ion-pair.

By combining ab initio calculations with experimental data, a proposed set of quasiparticle frontier orbitals for the [EMIM][OTF] ion pair was derived. These results also enabled critical evaluation of various computational methods. Hybrid functionals proved to be essential for accurately describing the electronic structure of this ion-pair. All calculation results were heavily dependent on the Hartree–Fock exchange content included. In general, functionals with a high percentage of HF exchange effectively removed spurious low-energy CT states that commonly affect TD-DFT calculations in large systems. However, most TD-DFT calculations still overestimated the intensity of the first part of the A_1,2 doublet. Optimally tuned RSH functionals proved to be superior to their standard off-the-shelf counterparts. The TD-DFT calculations from both the CAM-PBE50* and LRC-PBE0* showed good quantitative agreement with the experimental PIY spectrum.

The BSE calculation also helps to reveal the excited states. The BSE calculations based on the LC-BOP and LC-VV10 functionals offered the best overall agreement with the experimental PIY spectrum. Nonetheless, similar to TD-DFT, BSE calculations systematically underestimated CT excitation energies, often introduced spurious CT states and overestimated the intensity of the A_1,2 doublet.

Several functionals successfully reproduced the experimental vertical ionization energy using either the ΔSCF or GW approach. Among them, the ΔSCF calculation with the non-hybrid M11-L offered the most computationally efficient estimate.

In addition to the PBE-2X and the M08-SO functionals, the range-separated hybrid LRC-wPBE is the only functional to correctly predict the HOMO energy of −9.0 or −9.1 eV from ground state calculations. The so-called ‘magic’ of the GW method was also observed: it corrected HOMO energies, which varied from −10.6 eV to −8.1 eV across methods, to the experimentally consistent value of −9.0 or −9.1 eV. However, the GW-derived quasiparticle energies remained sensitive to the choice of basis set.

Indeed, no method succeeded in simultaneously providing the HOMO and the expected LUMO energies. Unfortunately, many studies still rely on (G)KS LUMO energies for predicting chemical properties, which is not recommended. However, if a zeroth-order estimation is necessary, the meta-GGA functionals should be used for the LUMO energy calculation.

In conclusion, the experimental data indicate that the vertical ionization energy is 9.05 eV, the first CT excitation occurs at 6.75 eV photon energy and the main CT transition is at 8.3 eV. According to the ab initio DFT philosophy, the HOMO energy should reflect the negative of the ionization potential and TD-DFT should provide the excited states. Among the 62 computational methods tested, only the PBE-2X functional quantitatively met these requirements. It accurately reproduced all key experimental energies within the measurement and computational uncertainties. The energy of the optical gap is described very well and there are no spurious CT states in the TD-DFT calculation.

Conflicts of interest

The authors declare that they have no known competing financial interests, intellectual property or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

All research data will be provided upon reasonable request.

Supplementary information (SI): Cartesian coordinates of the conformer used and a representative input file for ground state and TD-DFT calculations in MolGW. See DOI: https://doi.org/10.1039/d5ra07324c.

Acknowledgements

This work was supported by the Estonian Ministry of Education and Research (grant number TK210, Centre of Excellence in Sustainable Green Hydrogen and Energy Technologies). The activities of Estonian researchers at the MAX IV Lab were partially supported by Estonian Research Council core facility projects TT20 and TARISTU24-TK27. The project “Developing new research services and research infrastructures at MAX IV synchrotron radiation source” (2014-2020.4.01.20-0278) of the European Regional Development Fund is gratefully acknowledged. This study has received financial support from European Commission project EXANST (contract 101159716). We acknowledge MAX IV Laboratory for time on the FinEstBeAMS beamline under Proposals 20180539 and 20190253. Research conducted at MAX IV, a Swedish national user facility, is supported by the Swedish Research council under contract 2018-07152, the Swedish Governmental Agency for Innovation Systems under contract 2018-04969, and Formas under contract 2019-02496. Fruitful discussions with colleagues dr. Arvo Kikas, prof. Marco Kirm and dr. Harry Alles were helpful in preparing the manuscript. Assistance in running the various calculations was provided by dr. Juhan Matthias Kahk.

References

1. R. M. Fogarty, R. G. Palgrave, R. A. Bourne, K. Handrup, I. J. Villar-Garcia, D. J. Payne, P. A. Hunt and K. R. J. Lovelock, Electron spectroscopy of ionic liquids: experimental identification of atomic orbital contributions to valence electronic structure, Phys. Chem. Chem. Phys., 2019, 21, 18893.
2. J. M. Kahk, I. Kuusik, V. Kisand, K. R. J. Lovelock and J. Lischner, Frontier Orbitals and Quasiparticle Energy Levels in Ionic Liquids, npj Comput. Mater., 2020, 6, 148.
3. S. Rangan, J. Viereck and R. A. Bartynski, Electronic Properties of Cyano Ionic Liquids: a Valence Band Photoemission Study, J. Phys. Chem. B, 2020, 124, 7909–7917.
4. J. M. Seymour, E. Gousseva, A. I. Large, C. J. Clarke, P. Licence, R. M. Fogarty and D. A. Duncan, Experimental measurement and prediction of ionic liquid ionisation energies, Phys. Chem. Chem. Phys., 2021, 20957.
5. I. Kuusik, M. Berholts, J. Kruusma, V. Kisand, A. Tõnisoo, E. Lust and E. Nõmmiste, Valence electronic structure of [EMIM][BF4] ionic liquid: photoemission and DFT+D study, RSC Adv., 2018, 8(53), 30298–30304.
6. I. Kuusik, M. Berholts, J. Kruusma, A. Tõnisoo, E. Lust, E. Nommiste and V. Kisand, Valence electronic structure of [EMIM][B(CN)4]: ion-pair vs. bulk description, RSC Adv., 2019, 9, 33140.
7. J. Wang and Y. Wang, Strategies to Improve the Quantum Computation Accuracy for Electrochemical Windows of Ionic Liquids, J. Phys. Chem. B, 2024, 128, 1943–1952.
8. S. Kazemiabnavi, Z. Zhang, K. Thornton and S. Banerjee, Electrochemical Stability Window of Imidazolium-Based Ionic Liquids as Electrolytes for Lithium Batteries, J. Phys. Chem. B, 2016, 120, 5691–5702.
9. N. Ilawe, J. Fu, S. Ramanathan, B. Wong and J. Wu, Chemical and Radiation Stability of Ionic Liquids: A Computational Screening Study, J. Phys. Chem. C, 2016, 120(49), 27757–27767.
10. R. M. Fogarty, R. P. Matthews, P. A. Hunt and K. R. J. Lovelock, Accurate Prediction of Ionic Liquid Density-of-States from Low-Cost Calculations, Phys. Chem. Chem. Phys., 2025, 27(17), 9068–9075.
11. O. Höfft, S. Bahr, M. Himmerlich, S. Krischok, J. Schaefer and V. Kempter, Electronic structure of the surface of the ionic liquid [EMIM][Tf2N] studied by metastable Impact Electron Spectroscopy (MIES), UPS, and XPS, Langmuir, 2006, 22(17), 7120–7123.
12. M. Reinmöller, A. Ulbrich, T. Ikari, J. Preis, O. Höfft, F. Endres, S. Krischok and W. J. D. Beenken, Theoretical reconstruction and elementwise analysis of photoelectron, Phys. Chem. Chem. Phys., 2011, 13, 19526–19533.
13. S. Krischok, M. Eremtchenko, M. Himmerlich, P. Lorenz, J. Uhlig, A. Neumann, R. Öttking, W. J. D. Beenken, M. Himmerlich, P. Lorenz, O. Höfft, S. Bahr, V. Kempter and J. A. Schaefer, A Comparative Study on the Electronic Structure of the 1-Ethyl-3-Methylimidazolium Bis(trifluoromethylsulfonyl)amide RT-Ionic Liquid by Electron Spectroscopy and First Principles Calculations, Z. Physiol. Chem., 2006, 220, 1407–1416.
14. I. Kuusik, M. Kook, R. Pärna and V. Kisand, Ionic Liquid Vapors in Vacuum: Possibility to Derive Anodic Stabilities from DFT and UPS, ACS Omega, 2021, 6, 5255–5265.
15. I. Kuusik, M. Kook, R. Pärna, A. Kivimäki, T. Käämbre, L. Reisberg, A. Kikas and V. Kisand, The electronic structure of ionic liquids based on the TFSI anion: A gas phase UPS and DFT study, J. Mol. Liq., 2019, 294, 111580.
16. I. Lage-Estebanez, L. del Olmo, R. López and J. García de la Vega, The role of errors related to DFT methods in calculations involving ion pairs of ionic liquids, J. Comput. Chem., 2017, 38(8), 530–540.
17. I. Kuusik, M. Tarkanovskaja, J. Kruusma, V. Kisand, A. Tõnisoo, E. Lust and E. Nõmmiste, Valence band photoelectron spectra of [EMIM][BF4] ionic liquid vapor: Evidences of electronic relaxation, J. Mol. Liq., 2016, 223, 939–942.
18. I. Kuusik, M. Kook, R. Pärna and V. Kisand, Charge transfer and electronic relaxation effects in the photoemission of EMIM-DCA ionic liquid vapor, Chem. Phys., 2023, 111971.
19. M. Kook, I. Kuusik, R. Pärna, T. Käämbre, A. Kikas, A. Tõnisoo and J. M. Kahk, Ion fragmentation study of [EMMIM][TFSI], [EMIM][OTf] and [EMIM][DCA] by vacuum ultraviolet light, Int. J. Mass Spectrom., 2022, 471, 116732.
20. I. Kuusik, M. Tarkanovskaja, J. Kruusma, V. Reedo, R. Välbe, A. Lõhmus, V. Kisand, E. Lust, E. Kukk and E. Nõmmiste, Near threshold photodissociation study of EMIMBF4 vapor, RSC Adv., 2015, 5(9), 6834–6842.
21. K. Kooser, A. Kivimäki, P. Turunen, R. Pärna, L. Reisberg, M. Kirm, M. Valden, M. Huttula and E. Kukk, Gas-phase endstation of electron, ion and coincidence spectroscopies for diluted samples at the FinEstBeAMS beamline of the MAX IV 1.5 GeV storage ring, J. Synchrotron Radiat., 2020, 27(4), 1080–1091.
22. R. Pärna, R. Sankari, E. Kukk, E. Nõmmiste, M. Valden, M. Lastusaari, K. Kooser, K. Kokko, M. Hirsimäki, S. Urpelainen and M. Huttula, FinEstBeaMS A wide-range Finnish-Estonian Beamline for Materials Science at the 1.5 GeV storage ring at the MAX IV Laboratory, Nucl. Instrum. Methods Phys. Res., Sect. A, 2017, 859, 83–89.
23. K. Chernenko, A. Kivimäki, R. Pärna, W. Wang and R. Sankari, Performance and characterization of the FinEstBeAMS beamline at the MAX IV Laboratory, J. Synchrotron Radiat., 2021, 28(5), 1620–1630.
24. F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang and J. B. Neaton, molgw 1: Many-body perturbation theory software for atoms, molecules, and clusters, Comput. Phys. Commun., 2016, 208, 149–161.
25. E. Aprà, E. Bylaska, W. de Jong, N. Govind, K. Kowalski, T. Straatsma, M. Valiev, H. van Dam, Y. Alexeev, J. Anchell, V. Anisimov, F. Aquino, R. Atta-Fynn, J. Autschbach, N. Bauman, J. Becca, D. Bernholdt, K. Bhaskaran-Nair, S. Bogatko, P. Borowski, J. Boschen, J. Brabec, A. Bruner, E. Cauët, Y. Chen, G. N. Chuev, C. J. Cramer, J. Daily, M. J. O. Deegan, T. H. Dunning Jr, M. Dupuis, K. G. Dyall, G. I. Fann, S. A. Fischer, A. Fonari, H. Früchtl, L. Gagliardi, J. Garza, N. Gawande, S. Ghosh, K. Glaesemann, A. W. Götz, J. Hammond, V. Helms, E. D. Hermes, K. Hirao, S. Hirata, M. Jacquelin, L. Jensen, B. G. Johnson, H. Jónsson, R. A. Kendall, M. Klemm, R. Kobayashi, V. Konkov, S. Krishnamoorthy, M. Krishnan, Z. Lin, R. D. Lins, R. J. Littlefield, A. J. Logsdail, K. Lopata, W. Ma, A. V. Marenich, J. Martin Del Campo, D. Mejia-Rodriguez, J. E. Moore, J. M. Mullin, T. Nakajima, D. R. Nascimento, J. A. Nichols, P. J. Nichols, J. Nieplocha, A. Otero-de-la-Roza, B. Palmer, A. Panyala, T. Pirojsirikul, B. Peng, R. Peverati, J. Pittner, L. Pollack, R. M. Richard, P. Sadayappan, G. C. Schatz, W. A. Shelton, D. W. Silverstein, D. M. A. Smith, T. A. Soares, D. Song, M. Swart, H. L. Taylor, G. S. Thomas, V. Tipparaju, D. G. Truhlar, K. Tsemekhman, T. Van Voorhis, Á. Vázquez-Mayagoitia, P. Verma, O. Villa, A. Vishnu, K. D. Vogiatzis, D. Wang, J. H. Weare, M. J. Williamson, T. L. Windus, K. Woliński, A. T. Wong, Q. Wu, C. Yang, Q. Yu, M. Zacharias, Z. Zhang, Y. Zhao and R. J. Harrison, NWChem: Past, present, and future, J. Chem. Phys., 2020, 152, 184102.
26. S. Lehtola, C. Steigemann, M. J. T. Oliveira and M. A. L. Marques, Recent developments in Libxc - A comprehensive library of functionals for density functional theory, SoftwareX, 2018, 7, 1.
27. A. D. Laurent and D. Jacquemin, TD-DFT Benchmarks: A Review, Int. J. Quantum Chem., 2013, 113, 2019–2039.
28. K. Kanai, T. Nishi, T. Iwahashi, Y. Ouchi, K. Seki, Y. Harada and S. Shin, Electronic structures of imidazolium-based ionic liquids, J. Electron Spectrosc. Relat. Phenom., 2009, 174, 110–115.
29. L. Kronik, T. Stein, S. Refaely-Abramson and R. Baer, Functionals, Excitation Gaps of Finite-Sized Systems from Optimally Tuned Range-Separated Hybrid, J. Chem. Theory Comput., 2012, 8, 1515–1531.
30. M. Marques and E. Gross, Time-Dependent Density Functional Theory, Annu. Rev. Phys. Chem., 2004, 55, 427–455.
31. R. J. Bartlett, I. V. Schweigert and V. F. Lotrich, Ab initio DFT: Getting the right answer for the right reason, J. Mol. Struct., 2006, 771, 1–8.
32. M. A. Rohrdanz and J. M. Herbert, Simultaneous benchmarking of ground- and excited-state properties with long-range-corrected DFT, J. Chem. Phys., 2008, 129, 034107.
33. M. A. Rohrdanz, K. M. Martins and J. M. Herbert, A long-range-corrected density functional that performs well for both ground-state properties and time-dependent density functional theory excitation energies, including charge-transfer excited states, J. Chem. Phys., 2009, 130, 054112.
34. A. Patra, G. B. Pipim, A. I. Krylov and S. M. Sharada, Performance of Density Functionals for Excited-State Properties of Isolated Chromophores and Exciplexes: Emission Spectra, Solvatochromic Shifts, and Charge-Transfer Character, J. Chem. Theory Comput., 2024, 20, 2520–2537.
35. N. L. Nguyen, G. Borghi, A. Ferretti, I. Dabo and N. Marzari, First-Principles Photoemission Spectroscopy and Orbital Tomography in Molecules from Koopmans-Compliant Functionals, Phys. Rev. Lett., 2015, 114, 166405.
36. R. J. Bartlett, Towards an exact correlated orbital theory for electrons, Chem. Phys. Lett., 2009, 484, 1–9.
37. C. J. N. Coveney and D. Tew, A regularized second-order correlation method from Green's function theory, J. Chem. Theory Comput., 2023, 19, 3915–3928.
38. D. Golze, M. Dvorak and a. P. Rinke, The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy, Front. Chem., 2019, 7, 377.
39. L. Hedin, New method for calculating the one-particle green's function with application to the electron-gas problem, Phys. Rev. A, 1965, 139, 796.
40. M. Govoni and G. Galli, GW100: Comparison of Methods and Accuracy of Results Obtained with the WEST Code, J. Chem. Theory Comput., 2018, 14, 1895–1909.
41. D. Mejia-Rodriguez, E. Aprà, J. Autschbach, N. P. Bauman, E. J. Bylaska, N. Govind, J. R. Hammond, K. Kowalski, A. Kunitsa, A. Panyala, B. Peng, J. J. Rehr, H. Song, S. Tretiak, M. Valiev and F. D. Vila, NWChem: Recent and Ongoing Developments, J. Chem. Theory Comput., 2023, 19, 7077–7096.
42. A. Förster and F. Bruneval, Why Does the GW Approximation Give Accurate Quasiparticle Energies? The Cancellation of Vertex Corrections Quantified, J. Phys. Chem. Lett., 2024, 15, 12526–12534.
43. F. Bruneval, N. Dattani and M. J. van Setten, The GW Miracle in Many-Body Perturbation Theory for the Ionization Potential of Molecules, Front. Chem., 2021, 9, 749779.
44. S. Jana and J. Herbert, Slater transition methods for core-level electron binding energies, J. Chem. Phys., 2023, 158, 094111.
45. C. Holzer and Y. J. Franzke, A Guide to Molecular Properties from the Bethe–Salpeter Equation, J. Phys. Chem. Lett., 2025, 16, 3980–3990.
46. A. J. Cohen, P. Mori-Sánchez and W. Yang, Development of exchange-correlation functionals with minimal many-electron self-interaction error, J. Chem. Phys., 2007, 126, 191109.
47. R. L. A. Haiduke and R. J. Bartlett, Non-empirical exchange-correlation parameterizations based on exact conditions from correlated orbital theory, J. Chem. Phys., 2018, 148, 184106.
48. T. Körzdörfer and J.-L. Brédas, Organic Electronic Materials: Recent Advances in the DFT Description of the Ground and Excited States Using Tuned Range-Separated Hybrid Functionals, Acc. Chem. Res., 2014, 47, 3284–3291.
49. J. Perdew, A. Ruzsinszky, J. Tao, V. N. Staroverov, G. E. Scuseria and G. I. Csonka, Prescription for the design and selection of density functional approximations: More constraint satisfaction with fewer fits, J. Chem. Phys., 2005, 123, 062201.
50. T. Yanai, T. P. Tew and N. C. Handy, A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP), Chem. Phys. Lett., 2004, 393, 51–57.
51. Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai and K. Hirao, A long-range-corrected time-dependent density functional theory, J. Chem. Phys., 2004, 120, 8425–8433.
52. S. Refaely-Abramson, S. Sharifzadeh, N. Govind, J. Autschbach, J. B. Neaton, R. Baer and L. Kronik, Quasiparticle Spectra from a Nonempirical Optimally Tuned Range-Separated Hybrid Density Functional, Phys. Rev. Lett., 2012, 109, 226405.
53. T. Tsuneda, J.-W. Song, S. Suzuki and K. Hirao, On Koopmans' theorem in density functional theory, J. Chem. Phys., 2010, 133, 174101.
54. R. J. Bartlett, V. F. Lotrich and I. V. Schweigert, Ab initio density functional theory: The best of both worlds?, J. Chem. Phys., 2005, 123, 062205.
55. P. Verma and R. J. Bartlett, Increasing the applicability of density functional theory. IV.Consequences of ionization-potential improved exchange-correlation potentials, J. Chem. Phys., 2014, 140, 18A534.
56. Y. Jin and R. J. Bartlett, Accurate computation of X-ray absorption spectra with ionization potential optimized global hybrid functional, J. Chem. Phys., 2018, 149, 064111.
57. R. Baer and D. Neuhauser, Density Functional Theory with Correct Long-Range Asymptotic Behavior, Phys. Rev. Lett., 2005, 94, 043002.
58. J. F. Janak, Proof that sE/a n; = e; in density-functional theory, Phys. Rev. B: Condens. Matter Mater. Phys., 1978, 18, 7165–7167.
59. E. Brémond, A. J. Pérez-Jiménez, J. C. Sancho-García and C. Adamo, Range-separated hybrid density functionals made simple, J. Chem. Phys., 2019, 150, 201102.
60. E. Livshits and R. Baer, A well-tempered density functional theory of electrons in molecules, Phys. Chem. Chem. Phys., 2007, 9, 2932–2941.
61. Z. Lin and T. Voorhis, Triplet Tuning: A Novel Family of Non-Empirical Exchange–Correlation Functionals, J. Chem. Theory Comput., 2019, 15, 1226–1241.
62. E. Brémond, A. J. Pérez-Jiménez, J. C. Sancho-García and C. Adamo, Range-separated hybrid and double-hybrid density functionals: A quest for the determination of the range-separation parameter, J. Chem. Phys., 2020, 152, 244124.
63. R. Baer, E. Livshits and U. Salzner, Tuned Range-Separated Hybrids in Density Functional Theory, Annu. Rev. Phys. Chem., 2010, 61, 85–109.
64. P. Mori-Sanchez and A. J. Cohen, The derivative discontinuity of the exchange–correlation functional, Phys. Chem. Chem. Phys., 2014, 16, 14378.
65. J. Cornil, S. Vanderdonckt, R. Lazzaroni, D. Santos, G. Thys, H. Geise, L. Yu, M. Szablewski, D. Bloor, M. Lögdlund and W. Salaneck, Valence Electronic Structure of π-Conjugated Materials: Simulation of the Ultraviolet Photoelectron Spectra with Semiempirical Hartree–Fock Approaches, Chem. Mater., 1999, 11, 2436–2443.
66. F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Evers and M. J. van Setten, Quasi-Particle Self-Consistent GW for Molecules, J. Chem. Theory Comput., 2016, 12, 2528–2541.
67. M. J. van Setten, F. Caruso and S. Sharifzadeh, GW100: Benchmarking G0W0 for Molecular Systems, J. Chem. Theory Comput., 2015, 11, 5665–5687.
68. F. Bruneval, I. Maliyov, C. Lapointe and M.-C. Marinica, Extrapolating Unconverged GW Energies up to the Complete Basis Set Limit with Linear Regression, J. Chem. Theory Comput., 2020, 16, 4399–4407.
69. A. W. Lange, M. A. Rohrdanz and J. M. Herbert, Charge-Transfer Excited States in a π-Stacked Adenine Dimer, As Predicted Using Long-Range-Corrected Time-Dependent Density Functional Theory, J. Phys. Chem. B, 2008, 112(20), 6304–6308.
70. T. Lebeda, T. Aschebrock, J. Sun, L. Leppert and S. Kümmel, Right band gaps for the right reason at low computational cost with a meta-GGA, Phys. Rev. Mater., 2023, 7, 093803.
71. M. Bursch, H. Neugebauer, S. Ehlert and S. Grimme, Dispersion Corrected r2SCAN Based Global Hybrid Functionals: r2SCANh, r2SCAN0, and r2SCAN50, J. Chem. Phys., 2022, 156, 134105.
72. N. Mardirossian and M. Head-Gordon, Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals, Mol. Phys., 2017, 115(19), 2315–2372.
73. M. Li, J. R. Reimers, M. J. Ford, R. Kobayashi and R. D. Amos, Accurate prediction of the properties of materials using the CAM-B3LYP density functional, J. Comput. Chem., 2021, 42, 1486–1497.
74. C. Holzer, Y. J. Franzke and M. Kehry, Assessing the Accuracy of Local Hybrid Density Functional Approximations for Molecular Response Properties, J. Chem. Theory Comput., 2021, 17, 2928–2947.
75. Y. Jiang, Z. Hu, B. Zhou, C. Zhong, Z. Sun and H. Sun, Accurate Prediction for Dynamic Hybrid Local and Charge Transfer Excited States from Optimally Tuned Range-Separated Density Functionals, J. Phys. Chem. C, 2019, 123, 5616–5625.
76. E. Weintraub, T. M. Henderson and G. E. Scuseria, Long-Range-Corrected Hybrids Based on a New Model Exchange Hole, J. Chem. Theory Comput., 2009, 5, 754–762.
77. M. A. Mosquera, C. H. Borca, M. A. Ratner and G. C. Schatz, Connection between Hybrid Functionals and Importance of the Local Density Approximation, J. Phys. Chem. A, 2016, 120, 1605–1612.
78. J. R. Gómez-Pérez, F. A. Delesma, P. Calaminici and A. M. Köster, Accuracy of auxiliary density functional theory hybrid calculations for activation and reaction enthalpies of pericyclic reactions, J. Mol. Model., 2018, 24, 223.
79. S. Zahn, D. R. MacFarlane and E. I. Izgorodina, Assessment of Kohn–Sham density functional theory and Møller–Plesset perturbation theory for ionic liquids, Phys. Chem. Chem. Phys., 2013, 15, 13664.
80. C. A. McKeon, S. M. Hamed, F. Bruneval and J. B. Neaton, An optimally tuned range-separated hybrid starting point for ab initio GW plus Bethe–Salpeter equation calculations of molecules, J. Chem. Phys., 2022, 157, 074103.
81. P.-F. Loos, M. Comin, X. Blase and D. Jacquemin, Reference Energies for Intramolecular Charge-Transfer Excitations, J. Chem. Theory Comput., 2021, 17, 3666–3686.
82. C. Holzer and Y. J. Franzke, A General and Transferable Local Hybrid Functional for Electronic Structure Theory and Many-Fermion Approaches, J. Chem. Theory Comput., 2025, 21, 202–217.
83. M. Puiatti, D. M. A. Vera and A. B. Pierini, In search for an optimal methodology to calculate the valence electron affinities of temporary anions, Phys. Chem. Chem. Phys., 2009, 11, 9013–9024.
84. I. Anusiewicz, S. Freza, M. Bobrowski and P. Skurski, Electron attachment to representative cations composing ionic liquids, J. Chem. Phys., 2021, 154, 104302.
85. R. Vijayalakshmi, R. Anantharaj and A. B. Lakshmi, Evaluation of Chemical Reactivity and Stability of Ionic Liquids Using Ab Initio and COSMO-RS model, J. Comput. Chem., 2020, 41, 885–912.
86. M. I. Ofem, C. A. Ayi, H. Louis, T. E. Gber and A. A. Ayi, Influence of anionic species on the molecular structure, nature of bonding,reactivity, and stability of ionic liquids-based on1-butyl-3-methylimidazolium, J. Mol. Liq., 2023, 387, 122657.
87. C. Lian, H. Liu, C. Li and J. Wu, Hunting Ionic Liquids with Large Electrochemical Potential Windows, Thermodynamics and molecular-scale phenomena, AIChE J., 2019, 65(2), 804–810.

---