Journal Information
Journal ID (publisher-id): chemical
Title: Journal of the Korean Chemical Society
Translated Title (ko): 대한화학회지
ISSN (print): 1017-2548
ISSN (electronic): 2234-8530
Publisher: Korean Chemical Society대한화학회
INTRODUCTION▼
The electronic properties of organic mixed-valence (MV) systems have been extensively studied due to their critical role in charge transport, molecular electronics, and redoxactive materials.1-14 These systems, characterized by the coexistence of two distinct formal oxidation states within the same molecular framework, exhibit unique phenomena such as intervalence charge transfer (IVCT) and the separation of redox potentials for the two respective redox centers. These behaviors are influenced by the extent of electronic coupling (H12) between two diabatic states, each representing a fully charge-localized state.4,7 The Robin-Day classification2 provides a framework for categorizing MV systems based on the magnitude of H12. In class I, where H12 is negligible, the system exhibits fully charge-localized electronic states. Class II arises when H12 is smaller than the Marcus reorganization energy (λ), leading to partially delocalized electronic states. In contrast, class III occurs when H12 exceeds λ, resulting in fully delocalized states. Figure 1 depicts the perturbation of diabatic free energy surfaces as H12 increases, illustrating the transition from localized states (class I) to partially delocalized states (class II), and finally to fully delocalized states (class III).
The enduring interest in MV systems arises from their potential to exemplify the fundamental principles of electron transfer reactions.15-28 Materials employed in diverse applications, including charge transport, molecular electronics, and redox-active systems, frequently require properties characteristic of either class II or class III MV systems. The ability to precisely modulate these properties would significantly expand their practical applicability. As a result, MV systems near the class II/III borderline have attracted considerable attention in recent decades. Leveraging computational chemistry to predict the properties of MV systems with high accuracy during the design stage could save substantial experimental time and effort.
There have been numerous studies on MV systems that combine experimental and computational approaches, but research aimed at establishing appropriate computational methodologies remains scarce. The general considerations regarding computational methodologies can be summarized as follows:11,14,30 (i) Post-Hartree-Fock (HF) methods that include Møller–Plesset perturbation theory (MPn), coupled-cluster approaches (CCSD), and multi-configuration interaction (CASSCF), are generally too computationally demanding for the relatively large molecules commonly encountered in MV systems. (ii) HF theory, due to its lack of electron correlation, tends to overly localize or polarize the electronic properties of MV systems. (iii) As a result, the most practical approach is the use of density functional theory (DFT). However, the vast array of exchange-correlation functionals developed to date presents a challenge in selecting the most appropriate one for MV systems.
It is generally accepted that hybrid functionals incorporating exact Hartree-Fock exchange (EXexact) are superior to local density approximation (LDA), generalized gradient approximation (GGA), and meta-generalized gradient approximation (mGGA) functionals in terms of describing MV system properties.11,29-30 However, the optimal amount of EXexact required remains a subject of debate. Local and semilocal functionals often over-delocalize electronic states, while range-separated hybrid functionals may over-localize charge distributions. Kaupp proposed that global hybrid functionals with 35% EXexact (BLYP35) best describe experimental IVCT excitation energies, thermal electron-transfer barriers, and ESR hyperfine couplings.14,29 Nevertheless, in our recent study of the QPNPQ system (Chart 1),31 DFT calculations using the BMK functional, while successfully reproducing experimental IVCT excitation energies, revealed notable discrepancies in predicting electrochemical behavior. This underscores the importance of evaluating the performance of exchange-correlation functionals not only in predicting IVCT excitation energies but also in accurately capturing electrochemical behavior.
To address this, we systematically assessed the influence of different exchange-correlation functionals on the geometries and electronic properties of three prototypical organic MV systems: tetracene-1,4,7,10-tetraone (TBQ),15,17 hexahelicene-9,12,13,16-tetraone (HBQ),16,18 and 2,2'',5,5''-tetramethoxy-4,4''-dimethyl-1,1':4',1''-terphenyl (DphD).19 The experimental information regarding these compounds is summarized in Table 1.
Table1.
Experimental Data of the Title Systems
aThe Robin-Day class interpreted by author’s in the original paper.15-19 bCalculated (TBQ and HBQ) and X-ray crystallographic (DphD) results of centroid-to-centroic distances between two redox centers. cThe first reduction/oxidation potential referenced to SCE. dThe second reduction/oxidation potential referenced to SCE. eΔΦ = |Φ1-Φ2|. fAbsorption wavelength of the IVCT band maximum. gExtinction coefficient of IVCT band maximum. hFull width at half maximum of IVCT band.
TBQ features two benzoquinones connected by a naphthalene bridge within a planar tetracene framework, resulting in a delocalized π-electronic structure. Due to the complete coplanarity of the two redox centers and the bridging units, the Robin-Day classification of this compound appears to be class III. However, the authors interpret its classification as class II based on experimental data. HBQ, on the other hand, consists of two benzoquinones linked by a benzo[c] phenanthrene unit. Significant steric hindrance between the two terminal benzoquinone moieties causes them to adopt a helical structure due to distortion. The strong π-stacking interaction and the close centroid-to-centroid distance (3.7 Å) between the redox centers result in a large separation between the first and second reduction potentials (~380 mV), which supports its classification as a class III system. Additionally, the temperature independence of the electronic spin resonance (ESR) splitting further supports this classification.
In contrast to the two aforementioned systems, DphD is a cation radical MV system that comprises two dimethoxytoluene redox centers connected via a phenyl bridge. X-ray crystallographic data reveal that charge distribution parameter (ζ ) values of two redox centers (D1 and D2 for more positive and less positive moiety, respectively) are 0.8 and 0.2, respectively, which indicate borderline class II/III behavior.19 The H12 values determined through Mulliken-Hush analysis (~760 cm-1), redox potentials (~1600 cm-1) and crystal structure analysis (~2400 cm-1) highlight the intricate interchange between charge delocalization and localization in this compound.
These compounds were selected for several compelling reasons. First, all three systems exhibit borderline behavior between class II and class III rather than falling neatly into one category, making them excellent test cases for evaluating computational methods. Second, their structural diversity enhances the robustness of the analysis: TBQ features a coplanar arrangement, HBQ demonstrates strong π-stacking interactions, and DphD showcases a linearly connected bridge with variable dihedral angles between its constituent moieties. Third, experimentally accessible data that include redox potentials and IVCT excitation energies are available for these systems, allowing for direct comparisons with computational results. Finally, these compounds include bisquinone and bisdimethoxytoluene redox centers, which remain relatively unexplored compared to extensively studied dinitro and bistriphenylamine systems.
In this study, we evaluate 26 exchange-correlation functionals. By systematically analyzing their performance, we aim to elucidate the interplay between functional choice and electronic properties, ultimately offering guidelines for computational studies of MV systems.
COMPUTATIONAL MODELING▼
All calculations were conducted using the Gaussian 16 program package (A.03 version).32 Geometries corresponding to the stationary points on the ground state potential energy surface of all structures were fully optimized using DFT. The tested functionals included mGGA functionals (M11L,33 τHCTH34), hybrid GGA (H-GGA) functionals (B3LYP,35-37 B3PW91,37-39 mPW3PBE,36,40-41 B972,42 HSEH1PBE,43 mPW1PW91,41,44 N12SX,45 OHSE2PBE,43 LC-ωPBE,46-48 ωB97X,49 ωB97Xd,50 CAM-B3LYP,51), and hybrid mGGA (H-mGGA) functionals (τHCTHhyb,52 M06,53 BMK,54 M062X53). Additionally, eight of these functionals (B3LYP-D3, B3PW91-D3, LC-ωPBE-D3, ωB97Xd,55 CAM-B3LYP-D3, BMK-D3, M062X-D3) incorporated Grimme’s empirical dispersion corrections.56
Since TBQ●- and HBP●- are anionic MV system, all calculations were performed using the correlation-consistent basis sets with minimally augmented diffuse functions at the double-ζ level, originally developed by Dunning group and modified by Truhlar group (maug-cc-pVDZ).57-60 Charged species were treated as open-shell systems: monoanionic and monocationic species adopted doublet spin states, while dianionic and dicationic species adopted triplet spin states.
Geometry optimizations for all compounds were performed within a dielectric continuum using the self-consistent reaction field (SCRF) approach61-63 by the conductor-like polarizable continuum model (CPCM).64-65 The dielectric constants used corresponded to experimental conditions: dimethylformamide (εs = 37.781) for TBQ and HBQ, and dichloromethane (εs = 8.93) for DphD.
Vibrational frequency calculations were then followed under the same solvent environment in order to confirm the energy minima of the optimized geometries and to determine thermodynamic parameters, which were used to calculate the free energy (ΔG) of redox reactions.
The oxidation (Φox) and reduction (Φred) potentials were derived from the reaction free energies in solution using the equation:
where n is the number of electrons involved in the redox reaction, and F is the Faraday constant, 23.061 kcal mol-1V-1. For anionic MV systems, ΔG values for the first and second reduction steps were calculated as follows:
Similarly, for cationic MV systems, ΔG values for the first and second oxidation steps were derived from:
Experimental redox potentials, originally reported against the saturated calomel reference electrode (SCE), were converted to an absolute vacuum scale by using conventional reference potentials: SCE is 0.244 V vs. SHE, and SHE is 4.44 V vs. the vacuum level.
In general, solvation energies of redox couples were calculated using the Born-Haber cycle,66-67 as shown in Scheme 1, allowing computational cost savings by utilizing optimized gas-phase structures and single-point calculations in solution. However, imaginary frequencies often emerge when frequency calculations are performed in solution with gas-phase-optimized structures. Moreover, the gas-phase-optimized structures influence the MV classification (e.g., class III to class II shifts in borderline cases) since the solvent effects were particularly crucial for the geometries in solutions as polar solvents typically stabilize localized, valence-trapped states over delocalized charge distributions. To circumvent such errors, geometry optimizations, frequency calculations, and TD-DFT calculations were all conducted in a dielectric continuum to ensure consistency and avoid artifacts arising from gas-phase approximations in this study.14,68 The IVCT excitation energies of monocationic and monoanionic MV systems were obtained at their ground-state optimized geometries using TD-DFT calculations.
RESULTS AND DISCUSSION▼
In this study, we compare the experimental and computational results for key physical properties of MV systems, specifically focusing on geometries, first and second redox potentials, gaps between these potentials, and the IVCT excitation energies. The experimental data for these properties are summarized in Table 1.
As previously described, the computational results were obtained using 26 DFT functionals. The selection of these functionals was informed by a recent study by Araújo and coworkers,66 which evaluated the impact of exchange-correlation functionals on the description of electronic properties of non-fullerene acceptors in organic photovoltaics. The BMK functional was particularly included, which was identified as highly effective for evaluating MV systems.29-30
Among the range-separated hybrid functionals, CAM-B3LYP, LC-ωPBE, and ωB97X(d) were employed, using their default range-separation parameters (ω) of 0.33, 0.4, and 0.3, respectively. These parameters correspond to short-range/long-range transition distances of 3.0 Å, 2.5 Å, and 3.3 Å. Given that the shortest centroid-to-centroid distance in the studied systems (HBQ) is 3.7 Å, the long-range behavior of these functionals is particularly relevant for all compounds considered. Consequently, the amount of EXexact for these functionals were taken as 65% for CAM-B3LYP and 100% for LC-ωPBE and ωB97X(d) EXexact proportion category.
Two mGGA functionals, which lack EXexact, were also included in this study to establish reference points for evaluating the influence of EXexact across the full range from 0% to 100%.
Geometry
Evaluating the accuracy of molecular geometry optimized by computational chemistry is inherently challenging. Molecular geometry is characterized by key parameters that include bond lengths, bond angles, and dihedral angles. For molecules with the typical size of MV systems, these parameters can number in the dozens or even exceed 100. It is often a laborious task to compare such a large number of parameters with their respective references.
X-ray crystallographic data can be used as a reference for comparison; however, this approach has its own limitations. Crystal structures reflect molecular geometries in the solid state, which correspond to the geometries at very low temperatures, and are influenced by crystal packing effects, including intermolecular interactions or interactions with counterions. These interactions can differ significantly from the conditions under which electron-transfer behavior in MV systems is experimentally observed. For instance, a MV system may exhibit a localized structure in the crystalline state but adopt a delocalized structure in solution at ambient temperature. Additionally, the geometry of MV systems in solution is highly sensitive to solvent polarity. In polar solvents, charge localization is generally favored, whereas in nonpolar solvents, charge delocalization tends to dominate.
To evaluate the results of geometry optimization in this study, we employed charge distribution parameters (eq. 6) derived from quinoidal distortion at the redox centers, as depicted in Scheme 2.27-28,69-70
Scheme2.
Left: quinoidal distortions of 1,4-benzoquinone (top) and 1,4-dimethoybenzene (bottom). right: bond specification scheme.
Here, d0 and d1 denote the bond lengths of the redox center in the neutral and one-electron-reduced (or oxidized) states of Q (D), respectively, while di corresponds to the relevant bond length in the MV state. The values for d0 and d1 were derived from the optimized geometries of the Q (or D) moieties in their neutral and dianionic (or dicationic) states, respectively, whereas di was obtained from the optimized geometry of the MV system.
For anion radical MV systems such as TBQ●- and HBQ●-, we calculated the average of four parameters (α to δ), while for DphD●+ cation radical MV systems, we averaged five parameters (α to ε). The averaged parameter was simply denoted as ζ, while the individual ζ values for α to ε are represented as ζi. The ζ values were used to assess the degree of charge localization or delocalization in the MV systems.
It should be noted that each redox center has its own ζ value. For anion radical MV systems, two redox centers are designated as Q1 and Q2, with Q1 having a larger degree of negative charge. Similarly, two redox centers are denoted as D1 and D2, with D1 having a larger degree of positive charge for cation radical MV systems. The ζ value for Q1 (D1) should range from 0.5 to 1, and the sum of the ζ values for Q1 (D1) and Q2 (D2) should equal 1, representing the ideal case where charge separation occurs between the two redox moieties. However, in our calculations, there were instances where the sum of the ζ values for Q1 (D1) and Q2 (D2) exceeded 1. Using these ζ values in the calculation of Eq. 11 (vide infra) led to unintentional charge localization. To address this, all ζ values were normalized such that ζ1 + ζ2 = 1, where ζ1 and ζ2 correspond to the ζ values of Q1 (D1) and Q2 (D2), respectively. These normalized ζ values are summarized in Table 2, while the original values are provided in Table S1. The only instance where ζ exceeded 1 was observed for HBQ●- calculated with the BMK functional. In this case, the calculation failed to reproduce the geometry of the MV system and was therefore excluded from the table.
Table2.
The charge distribution parameters (ζ) of TBQ●-, HBQ●-, and DphD●+ calculated with different DFT functionals
aCalculated according to eq. 3. ζ is the averaged value of ζi (i = α, β, γ, and δ for anion radical systems and α, β, γ, δ, and ε for the cation radical system) bInclusion of Grimme’s D3 empirical dispersion correction. cInclusion of Grimme’s D2 empirical dispersion correction. d% Exact-exchange in the long range in the range-separated hybrid functional. eRemoved due to erroneous calculations.
The ζ values of D1/D2 for DphD●+ have been reported as 0.8/0.2 based on X-ray crystallography. If the DFT geometry optimization is accurate, these values should be reproduced.
However, as discussed earlier, the molecular structure in solution is likely more delocalized, implying that DFT calculations performed under the assumption of a dielectric continuum could yield a more symmetric value. Furthermore, given the general trend of increased charge localization with higher proportions of EXexact, it can be anticipated that the ζ values of D1/D2 will rise correspondingly. Nevertheless, the computational results deviated from this expectation. Most functionals predicted delocalized charge distribution across the range of functionals. Only the BMK, BMK-D3, M062X, and LC-ωPBE functionals predicted an asymmetric structure.
In fact, a ζ value of 0.8/0.2 does not clearly indicate a class II system, as this degree of charge separation is relatively high. Consequently, the authors have previously interpreted this compound as being on the borderline between class II and class III.19
Although the X-ray results suggest a partially chargelocalized structure, the possibility of delocalization in solution should also be taken into account, and computational predictions of delocalization cannot be entirely ruled out.
Nevertheless, it is difficult to consider the M062X and LC-ωPBE functionals superior based on their prediction of charge localization, as the partial charge distributions predicted by these functionals are excessively localized compared to the X-ray crystallographic results. In this context, the M062X and LC-ωPBE functionals can instead be regarded as yielding suboptimal results.
BMK and BMK-D3 predict charge distributions of 0.86/0.1 and 0.87/0.13, respectively, closely matching the experimentally determined values. Furthermore, these two functionals were the only ones to predict charge localization pattern with ζ values of approximately 0.86/0.14 for TBQ●-, which had been interpreted as class II.15,17 However, for HBQ●-, BMK yielded a physically unreasonable ζ value exceeding 1 for Q1, while BMK-D3 produced an extremely localized value of approximately 0.97/0.03. These results deviate significantly from the interpretation of this system as class III.16,18
The geometry optimization results for HBQ●- are intriguing. Despite being experimentally reported as class III, results of charge localized geometries were observed with several functionals that include B3LYP/B3LYP-D3 (20%), mPW1PW91 (25%), M06/M06-D3 (27%), BMK/BMK-D3 (42%), M062X-D3 (54%), and LC-ωPBE (100%).
The primary reason for including dispersion corrections was to assess their potential to improve the description of π-stacked systems. However, the results were inconsistent, making it difficult to draw definitive conclusions. Among the eight tested functionals, only M062X and LC-ωPBE showed a dependence of charge localization/delocalization tendencies on dispersion corrections. Specifically, dispersion corrections guided M062X toward charge localization, whereas LC-ωPBE exhibited the opposite trend, leaving the effect of dispersion corrections ambiguous. Given the classification of this system as class III,16,18 it can be concluded that B3LYP, mPW1PW91, M06, BMK, M062X-D3, and LC-ωPBE are unsuitable for accurately reproducing the structures of π-stacked systems, irrespective of whether dispersion corrections are applied. Therefore, BMK, when combined with dispersion corrections, demonstrates excellent performance for planar and linear systems. However, its performance for π-stacked systems remains inconclusive, underscoring the need for developing functionals capable of reliably modeling such systems.
Redox Potentials and Potential Gap
The first and second reduction potentials (Φred,1, Φred,2) of TBQ and HBQ and the first and second oxidation potentials (Φox,1, Φox,2) of DphD, were calculated according to Eq. 1, with ΔG values obtained from the corresponding half-wave reduction reactions (eqs. 2-5). The calculated Φ values, along with their corresponding ΔΦ values, are listed in Table 3. The difference between the calculated and experimental Φ values (ε(Φ) in V) of are illustrated in Fig. 2.
Table3.
The electrochemical redox potentials (Φ) of TBQ, HBQ, and DphD calculated with different DFT functionals.a
aAbsolute energy values in units of V. See table footnote of Table 1 for column 1 and 2. bFor the EXexact percentage, see Table 2. cThe first reduction potential of TBQ (HBQ). Calculated by the relation Φred,1 = -ΔG1/nF (F = Faraday constant, 23.061 kcal mol-1V-1) where ΔG1 is the free energy of TBQ (HBQ)+ e- → TBQ●- (HBQ●-) reaction. dThe second reduction potential of TBQ (HBQ). Calculated by the relation Φred,2 = -ΔG2/nF where ΔG2 is the free energy of TBQ●- (HBQ●-) + e- → TBQ2- (HBQ2-) reaction. eΔΦred = Φred,1 - Φred,2. fThe first oxidation potential of DphD. Calculated by the relation Φox,1 = -ΔG1/nF where ΔG1 is the free energy of DphD●+ + e- → DphD reaction. gThe second oxidation potential of DphD. Calculated by the relation Φox,2 = -ΔG2/nF where ΔG2 is the free energy of DphD2+ + e- → DphD●+ reaction. hΔΦox = Φox,2 – Φox,1. iMeasured in DMF, 0.1 M nBu4NBF4, SCE reference electrode. The values were converted to vs. vacuum level. jMeasured in CH2Cl2, 0.1 M nBu4NBF4, SCE reference electrode. The values were converted to vs vacuum level.
Figure2.
The difference between the calculated and experimental Φ values (ε(Φ) in V) of TBQ, HBQ, and DphD.
The ε(Φ) values exhibited varying trends depending on the functional and compound under investigation. Even within the same compound, error patterns diverged based on the amount of EXexact in the functional. For TBQ, most global hybrid functionals with EXexact proportions ranging from 15% (τHCTHhyb) to 27% (M06-D3) accurately predicted Φred,1, with the exception of OHSE2PBE, which overestimated by 0.27 V. The mGGA functionals M11L and τHCTH showed overestimations of approximately 0.2 V. Beyond BMK (≥ 42%), an increasing underestimation was observed, culminating in a significant -0.55 V underestimation by ωB97X. Regarding Φred,2, functionals with EXexact of 42% or less exhibited underestimations around -0.4 V. Starting from M062X (54%), the errors decreased to approximately -0.1 V. However, LC-ωPBE and LC-ωPBE-D3 showed exceptions, exhibiting an overestimation of about 0.13 V.
For DphD, the trends were inversed relative to TBQ. Functionals with EXexact less than 42% exhibited substantial underestimations of Φox,1 (~-0.4 V), while those with EXexact larger than 54% produced relatively accurate predictions, except for LC-ωPBE, which overestimated by 0.35 V. Predictions Φox,2 were generally precise across most global hybrid functionals, with the exceptions of MN12SX and OHSE2PBE, which overestimated by 0.17 V and 0.23 V, respectively. Conversely, range-separated hybrid functionals commonly underestimated by -0.33 to -0.44 V, except for LC-ωPBE, which provided accurate predictions. These opposing trends for TBQ (reduction) and DphD (oxidation) arise from the underlying energetic landscapes. The Φred values are governed by the free energy difference between the highest-energy neutral state and the intermediate monoanionic state (first reduction) and between the intermediate monoanionic state and the lowest-energy dianionic state (second reduction). The Φox values on the other hand, reflect the free energy differences between the lowest-energy neutral state and the intermediate monocationic state (first oxidation) and between the intermediate monocationic state and the highest-energy dicationic state (second oxidation). Consequently, the apparent inverse trends stem from a shared tendency to underestimate energy differences between the lowest and intermediate states in these systems. For HBQ, global hybrid functionals with EXexact less than 42% provided accurate predictions for both Φred,1 and Φred,2 values. Exceptions included OHSE2PBE, which overestimated both potentials by over 0.2 V, and M06, which severely underestimated the Φred,2 by -0.62 V. Functionals with EXexact larger than 54% displayed ~-0.4 V underestimations for the Φred,1 and ~0.4 V overestimations for the Φred,2. Notably, dispersion corrections significantly reduced errors in M06 and M062X predictions.
While this study analyzed errors in individual redox potentials to evaluate the general performance of exchange-correlation functionals, the principal aim in predicting the redox behavior of MV systems lies in accurately estimating the potential gap (ΔΦ) between consecutive redox reactions. This gap is particularly significant in MV systems, as it directly correlates with the resonance stabilization energy (ΔGr) via the equation:19
Under the class II and class III limits, ΔGr can be related as following two equations:7
These relationships enable the evaluation of H12 values, which is central to classifying MV systems. In class II, λ corresponds to the maximum energy of the IVCT band (λ = νmax) and can be readily obtained, whereas under the class III limit, λ must be independently determined, adding complexity.4,7
The H12 values derived from ΔΦ can be compared to those obtained through spectroscopic methods, such as the Mulliken-Hush relation (eq. 10) for class II systems or the H12 = νmax/2 relation for class III systems. Such comparisons provide a robust basis for accurately classifying the MV system and the discussion regarding H12 values will be addressed in Section IV.
The error trends in ΔΦ (ε(ΔΦ)) showed a clear dependence on the amount of EXexact (Fig. 3). As the amount of EXexact increased, ε(ΔΦ) generally decreased, transitioning from overestimations at lower amount of EXexact to underestimations at higher amount of EXexact. For TBQ, overestimations persisted up to EXexact less than 42%, while underestimations emerged from EXexact larger than 54%. Similar trend was observed for DphD with the exception of M062X-D3 that showed 0.21V of overestimation. For HBQ, most global hybrid functionals predicted relatively accurate ΔΦ values with absolute errors below 0.1 V. The notable exception was M06, which exhibited a pronounced overestimation of 0.56 V. BMK showed marginal -0.2V of underestimation. The mGGA functionals (M11L and τHCTH) showed overestimations of 0.23 V and 0.15 V, respectively.
Figure3.
The difference between the calculated and experimental ΔΦ values (ε(ΔΦ) in V) of TBQ, HBQ, and DphD along with their mean absolute error (MAE).
It is important to highlight that the pronounced underestimations of ΔΦ values observed for functionals with EXexact larger than 54% across all three MV systems arise from the potential inversions predicted by these functionals.68,71 Notable exceptions to this trend include the results calculated for HBQ using M062X-D3 (ΔΦ = 0.30V) and for DphD using LC-ωPBE (ΔΦ = 0.26V). In all other cases with EXexact larger than 54%, ΔΦ values were negative.
A negative ΔΦ value implies that the Frost diagram for the two sequential reduction or oxidation reactions is concave, indicating that an MV system with an intermediate oxidation state would immediately disproportionate into the neutral and fully reduced or oxidized states.28,70 Consequently, the experimental observation of such MV systems becomes impossible. These predictions starkly contradict experimental evidence confirming the existence of these three MV systems, rendering the results evidently erroneous.
The consistent prediction of potential inversion not only by range-separated hybrid functionals but also by functionals with large amount of EXexact underscores their unsuitability for accurately modeling the redox behavior of MV systems.
The mean absolute error (MAE) of ε(ΔΦ) across the three compounds serves as a reliable benchmark for evaluating functional performance. As amount of EXexact increased, MAE decreased but showed an upturn beyond M062X. The functionals M06-D3, BMK, BMK-D3, and M062X-D3 (27% ≤ EXexact ≤ 54%) emerged as the most effective, with MAE values of 0.260, 0.201, 0.172, and 0.266 V, respectively. Notably, when dispersion corrections were omitted from M06-D3 and M062X-D3, errors were significantly increased. Nevertheless, the fact that the MAE values for even the most proficient functionals exceed 0.17 V underscore a fundamental limitation of DFT in predicting ΔΦ with sufficient accuracy. Given that a ΔΦ difference as small as 0.1 V can drastically alter the Robin–Day classification and properties of MV systems, further fine-tuning of exchange-correlation functionals is imperative for reliable predictions of redox properties in these systems. However, the geometry optimization results presented in Table 2 reveal that range-separated hybrid functionals optimized the structures of MV systems as delocalized geometries except LC-ωPBE. This observation indicates that the overestimation of IVCT excitation energy by range-separated hybrid functionals is not solely attributable to their tendency toward charge localization. Even in fully delocalized geometries, these functionals overestimate IVCT excitation energy. Further in-depth investigations are required to elucidate the underlying causes of this behavior.
IVCT Excitation Energy
The calculated IVCT excitation energies of three MV systems are compiled in Table 4 and the differences between calculated and experimental Eex values (ε(Eex)) are displayed in Fig. 4. The trends in IVCT excitation energies obtained from TD-DFT calculations showed consistent patterns across the three MV systems. Functionals with low amount of EXexact (0% ≤ EXexact ≤ 27%) exhibited relatively small errors, while those with intermediate amount of EXexact (42% ≤ EXexact ≤ 65%) displayed moderate errors. Functionals with 100% of EXexact consistently produced substantial overestimations. Specifically, the error ranges for TBQ●−, HBQ●−, and DphD●+ were 0.11–0.23, 0.12–0.28, and 0.09–0.20 eV, respectively, for low amount of EXexact; 0.21–0.58, 0.30–0.75, and 0.17–0.39 eV, respectively, for intermediate amount of EXexact; and 0.60–0.63, 1.00–1.28, and 0.65–1.13 eV, respectively, for 100% of EXexact.
Figure4.
The difference between the calculated and experimental Eex values (ε(Eex) in V) of TBQ, HBQ, and DphD along with their mean absolute error (MAE).
Table4.
The IVCT excitation energy (Eex) of TBQ●-, HBQ●-, and DphD●+ calculated via TD-DFT calculation.a
At low amount of EXexact, TBQ●− showed overestimations, while HBQ●− and DphD●+ showed underestimations. Intermediate amount of EXexact yielded mixed trends, with alternating over- and underestimations. Functionals with 100% of EXexact invariably produced large overestimations, reflecting an excessively localized treatment of MV systems. This behavior is consistent with previous reports, such as those by Kaupp et al., which noted that higher amount of EXexact lead to over-localized geometries and exaggerated IVCT excitation energies.14,72-76
To assess the predictive performance for excitation energies, the mean absolute error (MAE) was employed as the evaluation metric, as in previous sections. Functionals with low amount of EXexact values achieved relatively low MAEs, ranging from 0.14 to 0.20 eV. By contrast, intermediate amount of EXexact values resulted in increased MAEs, ranging from 0.35 to 0.42 eV, while functionals with 100% of EXexact displayed MAEs exceeding 0.8 eV.
Kaupp and coworkers had recommended BLYP35 (35%) as the optimal functional for predicting IVCT excitation energy and suggested BMK (42%) as a viable alternative. In the present study, a correlation between the amount of EXexact and MAE was observed, particularly for global hybrid functionals within the 15% to 27% range of EXexact, where an increase in EXexact led to a reduction in MAE. However, this trend reversed sharply at 42% EXexact, suggesting that the minimum MAE may occur around 35% EXexact or within the 25% to 35% range.
Despite Kaupp's recommendation, BMK (42%) exhibited substantial MAEs exceeding 0.35 eV across all three MV systems, rendering it unsuitable for IVCT excitation energy predictions. In contrast, M06 (27%) demonstrated a low MAE of 0.13 eV, indicating its viability for predicting IVCT excitation energies, particularly when used with dispersion corrections, as previously emphasized in Sections I and II. The performance of mPW1PW91 (25%) was similarly reliable. These results emphasize the crucial importance of functional selection in accurately modeling IVCT excitation energies and underscore the necessity of identifying optimal EXexact content for reliable predictions.
Electronic coupling
We evaluated H12 using four different methods. First, the Mulliken-Hush (MH) relation (Eq. 10) was applied with the νmax and oscillator strength (f) derived from TD-DFT calculations (H12,MH). For the centroid-to-centroid distance between two redox centers (R12), fixed values of 7.35 Å, 3.7 Å, and 8.6 Å were used for TBQ●−, HBQ●−, and DphD●+, respectively. Second, H12 was determined using Eq. 8 under the assumption that the MV system corresponds to the class II limit (H12,II). Since the reorganization energy (λ) is equivalent to Eex in class II, this value can be readily viable from TD-DFT calculation. Meanwhile, ΔGr was calculated using Eq. 7 using ΔΦ values. Third, assuming the class III limit, H12 was calculated using the relationship Eex = H12/2 (H12,III ). Lastly, H12 was determined via Eq. 11 using ζ, which was obtained from geometric parameters analyzed in section I, while ΔE represents Eex from TD-DFT calculation (H12,ζ). The theoretical background of Eq. 11 has been extensively discussed in precedent studies, including our previous work.4,7,27-28 The four types of H12 values are summarized in Table 5.
Table5.
The ζ based H12 of TBQ●-, HBQ●-, and DphD●+.a
aSee table footnote of Table 2 for column 1. bEvaluated according to eq. 10. cEvaluated according to eq. 8 in class II limit. dEvaluated according to eq. 9 in class III limit. eEvalated according to eq. 11. fOscillator strength is 0. gDue to potential inversion.
An examination of Table 5 reveals that for DphD●+, assigned as a class II/III borderline system in the original paper,70 the four H12 values differ significantly. This suggests that the H12 results can vary depending on how the system is classified and which experimental parameters are used for the calculations. In contrast, for TBQ●− and HBQ●−, the H12,II and H12,III values are similar regardless of classification, with H12,MH being notably smaller than the other two. Furthermore, since four types of H12 were derived using parameters already assessed in this work, evaluating their deviations from experimentally determined values would be redundant. As a result, we present the values in Table 5 as reference data without further performance evaluations.
CONCLUSION▼
This study provides a systematic evaluation of the performance of 26 exchange-correlation functionals in predicting the electronic properties of three prototypical organic MV systems: TBQ, HBQ, and DphD. By assessing geometry optimization, redox potentials, IVCT excitation energies, and electronic coupling values, the interplay between functional choice and computational accuracy was thoroughly investigated.
Key findings include the following: (1) For geometry optimization, BMK-D3 demonstrates excellent performance for planar and linear systems but fails to reproduce π-stacked systems, underscoring the need for developing functionals capable of reliably modeling such systems. (2) Redox potential predictions demonstrated that global hybrid functionals with moderate amount of EXexact (27–54%) generally performed better, with M06-D3 emerging as the most reliable functional for accurately capturing ΔΦ. Functionals with EXexact ≥ 54% often predicted potential inversions, rendering them unsuitable for modeling the redox behavior of MV systems. (3) For IVCT excitation energies, low-to-moderate EXexact (15–35%) functionals provided the most accurate predictions. In particular, M06-D3 and mPW1PW91 achieved low mean absolute errors, highlighting their utility for modeling optical properties.
These results indicate that the optimal functional choice for MV systems depends on the property of interest. While M06-D3 consistently delivered robust performance across all categories, no single functional was universally superior. Instead, tailoring the functional selection to the specific system and property remains crucial. This study provides valuable reference data for such tailoring. The insights gained from this study highlight the limitations of current density functional theory approaches and underscore the need for further development of exchange-correlation functionals tailored for mixed-valence systems. Such advancements would significantly enhance the predictive power of computational methodologies, facilitating the rational design of redox-active materials and molecular electronic systems.