Assessment of Ab Initio and Density Functional Theory Methods for the Excitations of Donor-Acceptor Complexes: The Case of the Benzene-Tetracyanoethylene Model

The understanding of the excited-state properties of electron donors, acceptors and their interfaces in organic optoelectronic devices is a fundamental issue for their performance optimization. In order to obtain a balanced description of the different excitation types for electron-donor-acceptor systems, including the singlet charge transfer (CT), local excitations, and triplet excited states, several ab initio and density functional theory (DFT) methods for excited-state calculations were evaluated based upon the selected model system of benzene-tetracyanoethylene (B-TCNE) complexes. On the basis of benchmark calculations of the equation-of-motion coupled-cluster with single and double excitations method, the arithmetic mean of the absolute errors and standard errors of the electronic excitation energies for the different computational methods suggest that the M11 functional in DFT is superior to the other tested DFT functionals, and time-dependent DFT (TDDFT) with the Tamm–Dancoff approximation improves the accuracy of the calculated excitation energies relative to that of the full TDDFT. The performance of the M11 functional underlines the importance of kinetic energy density, spin-density gradient, and range separation in the development of novel DFT functionals. According to the TDDFT results, the performances of the different TDDFT methods on the CT properties of the B-TCNE complexes were also analyzed.


Introduction
In terms of device architecture and materials [1], promising organic devices (organic photovoltaics, organic light-emitting diodes, organic photon detectors, etc.) usually contain electron donor and acceptor materials constructing an active layer [2,3]. Charge transports in donor/acceptor materials and charge transfer (CT) at the interfaces between electron donors and acceptors have great importance in the improvement of device performance [4][5][6][7][8][9][10][11][12]. Usually, the charge transfer in organic optoelectronic devices involves molecular excitation behavior [13,14]. Molecular excitation properties, such as energy level landscape and local and CT transition characters, determine various photophysical and photochemical properties-including florescence and phosphorescence radiation, internal conversion, intersystem thermochemistry and kinetics, etc. For LRCRSH functionals, the Coulomb potential operator is divided into short-range (SR) and long-range (LR) parts by using error function 1 r = er f (ωr) r + er f c(ωr) r (1) where ω is the range-separation parameter [36]. Correspondingly, the SR exchange is usually represented by local or semi-local potential, whereas the LR exchange is adopted as exact HF exchange.
It should be mentioned that the dynamic exchange-correlation kernel was presented for effective SR electron-electron interaction, and its accuracy can compete with the actual static exchange-correlation kernels [38,39]. Due to the correct LR asymptotic behavior, LRCRSH functionals generate more accurate results for LR CT. Recently, it was reported that the ω can be tuned based upon Koopmans' theorem to obtain certain expected physical quantities, such as ionization potential and fundamental gap [40][41][42][43]. Also, it was pointed out that the optimally tuned ω depends on the system under study [44]. It is therefore more accurate to say that the optimally tuned ω depends not only on the system and state, but also on the properties (CT or others). The LRCRSH functionals with the optimally tuned ω can generate sufficiently accurate results for many systems [45][46][47][48][49]. But the important issue that should be noted is that if we adopt the LRCRSH functionals with the optimally tuned ω for different systems, the calculated results will be rooted in different exchange potentials in DFT. Undoubtedly, this reduces the comparability of computational studies for different systems, which is important in the design of novel materials. The performance of some DFT functionals in CT of donor/acceptor complexes, such as tetrathiafulvalene-tetracyanoquinodimethane (TTF-TCNQ) [50][51][52] and pentacene-C60 [45,53,54], have been evaluated. Several DFT functionals, including M11 [55], SOGGA11X [56], N12SX/ MN12SX [57], HISSbPBE [58], and APF [59] had been presented and implemented in a quantum chemistry package recently. But their performances in molecular singlet/triplet excitations (CT or local excitations) were not systematically evaluated. Meanwhile, regarding the full TDDFT and TDA-TDDFT, which one is the better choice for the description of singlet and triplet excitations for molecular systems? Furthermore, in order to generate a balanced description of singlet/triplet CT and local excitations, is it possible to find a general method without further tuning the empirical parameters (such as the optimally tuned range-separation parameter ω) and do so at an affordable computational cost? Therefore, in this work, in order to highlight the effects of the theoretical methods and to avoid experimental uncertainties (temperature, vibronic coupling, and zero-point effects etc. [60]), the selected benchmark is the EOM-CCSD method because of its higher accuracy for molecular excitation [61,62]. The robust computation of excitation energies by the EOM-CCSD method has been further approved in [63]. Based on the benchmark calculations, we assessed several wavefunction-based methods and the TDDFT methods (both the full TDDFT and the TDA-TDDFT) with different functionals for the excitations of the benzene-tetracyanoethylene (B-TCNE) complexes, which was selected as the reference donor-acceptor model system.

The Exchange-Correlation Energy in the Tested DFT Functionals
In order to understand the different performances of the excited-state calculations, it is necessary to analyze the exchange-correlation energy functionals. The tested DFT functionals in this work can be classified into four types. The first is global hybrid functionals, including the SOGGA11X and APF functionals. The second type of functional is range-separated hybrid functionals, which include LC-ωPBE, ωB97XD, N12SX/MN12SX, and HISSbPBE. Hybrid meta-GGA M06-2X [64] is the third type of functional. The fourth type of functional is the range-separated hybrid meta-GGA functional M11. In order to clearly present the differences, the exchange-correlation energy functionals are listed in Table 1, where the superscripts and subscripts in the formulas stand for method and

Singlet Excitations
The EOM-CCSD/SVP benchmark results for the B-TCNE complexes are listed in Table 2, including three singlet and triplet excitation energies (eV), wavelength (nm), and main transition configurations with coefficients larger than 10%. Figure 1 shows the molecular orbitals (MOs) involved in these transitions. In terms of the MOs in Figure 1 and transition configurations in Table 2, the first singlet excited state for the B-TCNE complexes can be assigned as the CT excitation (labeled as S-CT 1 ), while the main character of the third singlet excited state can be assigned as the locally excited (LE) state (labeled as S-LE 1 ), though it exhibits hybridization of the LE state and CT. The calculated vertical electronic excitation energies (eV) of the S-CT 1 and S-LE 1 for the B-TCNE complexes are listed in Table 3. The EOM-CCSD/SVP benchmark results of the S-CT 1 and S-LE 1 are 4.19 and 5.34 eV, respectively. In order to give the computational methods accuracy, the calculated excitation energy errors (defined as the excitation energy difference between the giving method and the benchmark) of S-CT 1 and S-LE 1 are presented in Figures 2 and 3.     The excitation energies of S-CT 1 and S-LE 1 calculated with EOM-CCSD/LanL2DZ are about 0.28 and 0.03 eV higher than those of the EOM-CCSD/SVP results, underlining that the basis sets incompleteness generate more significant effects on the CT excitations than those on the LE state. The excitation energies of S-CT 1 and S-LE 1 calculated by LC-ωPBE/SVP are about 0.11 eV higher than those of LC-ωPBE/aug-cc-pVTZ. The tendencies indicated by these data are that (i) the larger basis sets generate smaller singlet CT and LE energies and (ii) the dependence of ab initio methods are more sensitive on basis sets than those of DFT methods. The TDA results of LC-ωPBE/aug-cc-pVTZ and LC-ωPBE/SVP also support these tendencies.
For the results of the wavefunction-based methods, the error of S-CT 1  The good performance of the LC-ωPBE functional can be attributed to its range-separation parameter, which is relatively suitable for the B-TCNE complexes. The ωB97XD and ωB97XD(0.23) results indicate that the LRCRSH functionals with the optimally tuned range-separation parameter ω improve computational accuracy for singlet CT and LE excitations through the variation of HF exchange in E XC For the hybrid meta-GGA functional, the S-CT 1 and S-LE 1 excitation energies of M06-2X are about 1.18 and 0.75 eV lower than those of benchmark results, respectively. The accuracy of the M06-2X functional for singlet excitations is similar to those of the ωB97XD and ωB97XD(0.23) functionals due to its high nonlocality functional with double the amount of non-local exchange [35]. For the range-separated hybrid meta-GGA M11 functional, the S-CT 1 and S-LE 1 excitation energies are about 0.62 and 0.69 eV lower than those of benchmark results, respectively, generating a balanced description of the LE state and the CT excitation with fewer excitation energy errors. According to these descriptions and the data presented in Figures 2 and 3, for most of the tested DFT functionals, the excitation energy errors are not systematic, and the CT excitation energy errors are larger than those of the LE state because the delocalization errors in E XC generate more significant influence on the CT excitations. The exchange-correlation kernel in the tested DFT functionals contains local/semi-local potential and non-local HF exchange to different extents, meaning that the delocalization errors cannot be completely dismissed and can be partially remedied by increasing HF exchange in the DFT functionals. Sini et al. reported that, for TTF-TCNQ complexes, the more content of HF exchange in the DFT functionals, the higher the CT excitation energy [65]. In this case, the smaller content of HF exchange in the DFT functionals, such as the APF/MN12SX/N12SX/HISSbPBE functionals, the larger the excitation energy errors (for both the CT and the LE state).
For the TDA-TDDFT results with aug-cc-pVTZ basis sets (see Table 3), the excitation energies of S-CT 1 and S-LE 1 range from 2.03~3.91 eV and 4.58~5.08 eV, respectively. The excitation energy errors of S-CT 1 and S-LE 1 , which are presented in Figures 2 and 3, range from 0.27~2.16 eV and 0.26~0.76 eV, respectively. The functional effects in the TDA-TDDFT are similar to those in the full TDDFT. Compared with the results of the full TDDFT, the TDA-TDDFT significantly improves the accuracy of LE (about 0.3 eV), underlining the rationality of ignoring the de-excitation contributions in the TDA-TDDFT. But the accuracy improvement of the CT is tiny, suggesting that the de-excitation contributions for the LE state are more remarkable than those for CT. Another possible reason for the improvement of the TDA-TDDFT might be the compensation effects between the exchange-correlation DFT kernel and the TDA.

Triplet Excitations
In terms of triplet excitations data (listed in Table 2) and related MOs presented in Figure 1, the first triplet excited state (labeled as T 1 ) is LE in TCNE with π→π* transition, and the EOM-CCSD/SVP benchmark excitation energy is 2.71 eV. The T 1 excitation energy of the B-TCNE complexes calculated by the selected ab initio and DFT methods are also listed in Table 3. The calculated excitation energy errors are shown in Figure 4. The T 1 excitation energy of the SAC-CI result is about 0.19 eV smaller than that of the benchmark. The excitation errors indicate that the SAC-CI method is superior to the CIS and CIS(D) methods for triplet excitations. For the full TDDFT results, the calculated T 1 excitation energy range from 1.50~2.44 eV, and the corresponding excitation energy errors are about 0.27~1.21 eV. The M06-2X functional generates the smallest excitation error (~0.27 eV). The excitation energy errors of the LC-ωPBE and HISSbPBE functionals are larger than those of other functionals, indicating that these functionals are not suitable for triplet excitations. The other functionals have similar accuracy. The common feature of the LC-ωPBE and HISSbPBE functionals are their PBE correlation. In spin-flip transition, the importance of electronic correlation is increasing because correlation exists in electrons with opposite spin. For the TDA-TDDFT results, the T 1 excitation energies range from 2.01~2.72 eV, and the accuracy of the T 1 excitation energies is significantly improved. It was found that applying TDA upshifts the triplet excitation energies and reduces the triplet excitation errors for the selected organic molecules [66]. The T 1 excitation energy difference between the TDA-TDDFT results and the full TDDFT results also emphasizes the remarkable de-excitation contributions in the full TDDFT and again underlines the rationality of ignoring the de-excitation contributions in the TDA-TDDFT for triplet excitation. Meanwhile, the performances of the M06-2X and M11 functionals in the TDA-TDDFT are superior to that of other functionals, and the excitation energy errors of the APF/MN12SX/N12SX/HISSbPBE functionals, which contain small contents of HF exchange are larger than those of other functionals. In addition, it was reported that the adiabatic approximation error for the S 0 →T 1 spin-flip transition is larger than that for the singlet transition [67]. Therefore, adiabatic approximation of TDDFT is another reason for excitation energy error.

Assessment of Different Methods for Excited States
In order to quantitatively evaluate the performance of the different computational methods for the excitations of the donor-acceptor B-TCNE complexes, we defined the arithmetic mean of the absolute errors (AMAE) and standard errors (SE) of the vertical electronic excitation energies for the B-TCNE complexes as respectively. In these expressions, the i presents the lowest singlet local/CT excitation or the first triplet excited state, N = 3 means the average of the abovementioned three kinds of excited states, the E i is the corresponding benchmark excitation energies calculated with EOM-CCSD/SVP, and the E m j is the corresponding excitation energies calculated with other methods. In terms of the EOM-CCSD/SVP benchmark results, Table 4 lists the calculated AMAE and SE of the vertical excitation energies for the B-TCNE complexes. Compared with the CIS and CIS(D) results, the smaller AMAE and SE of the SAC-CI method mean it performs better in excitation calculations. According to the full TDDFT results, the smallest SE of the M11 functional indicates that it performs best in showing a balanced description of singlet LE/CT and triplet excitations for the B-TCNE complexes. For the TDA-TDDFT, the M11 functional has the smallest SE except for the LC-ωPBE functional. The AMAE and SE calculated using LC-ωPBE/aug-cc-pVTZ and (TDA) LC-ωPBE/aug-cc-pVTZ supports the finding that the smallest AMAE and SE of (TDA) LC-ωPBE/aug-cc-pVTZ results from TDA, rather than a DFT functional. Furthermore, the AMAE and SE for a given functional in the TDA-TDDFT are smaller than those in the full TDDFT, underling the effectiveness of TDA in excitation calculations. It should be noted that the benchmark EOM-CCSD data raises uncertainty because of the relatively small basis sets of SVP. For the B-TCNE complexes in this work, according to the basis sets effects on the EOM-CCSD data for the tested molecules [62,68] and the low-lying excited states that were calculated in this work, the estimated average uncertainty of the EOM-CCSD data for the S-LE 1 , S-CT 1, and T 1 excitations should be less than 0.1 eV if the basis sets can be extended from SVP to the larger. However, this uncertainty of the EOM-CCSD data does not affect the assessment of the DFT functionals because the AMAE for the excitation energies calculated using the TDDFT methods with different functionals are remarkably larger than the uncertainty of the EOM-CCSD benchmark data.

Transferred Charges from TDDFT
In order to further investigate the computational method effects on the CT properties, the calculated transferred charges q CT (in e) and the CT excitation length D CT (in Å) for the S-CT 1 state of the B-TCNE complexes are listed Table 5. Because the electron densities of EOM-CCSD, CIS/CIS(D), SAC-CI, and LC-ωPBE are not available with the density keyword in the Gaussian package, the CT properties were not investigated using these methods. For a given DFT functional, the q CT and D CT results of the full TDDFT are slightly larger than those of the TDA-TDDFT (about 0.001 e and 0.001 Å, respectively), suggesting that the full TDDFT and the TDA-TDDFT generate similar electron density variation induced by excitation. The q CT and the CT excitation length D CT range from 1.22 to 1.49 e, and from 1.78 to 2.10 Å, respectively. The remarkable q CT confirms the CT character of this excited state. The calculated D CT is smaller than the distance between molecular centers because of polarization effects during CT. It also can be found that the larger q CT corresponds to the shorter D CT . For instance, the SOGGA11X functional in the full TDDFT generates the largest q CT (~1.49 e), while the D CT is the shortest (~1.78 Å) generated. Considering the accuracy of the M11 functional for excitation energy, the differences of transferred charges (∆q, in e) between the TDDFT results with different functionals and the full TDDFT results with the M11 functional are presented in Figure 5 in order to compare the effects of the methods on the CT properties. Apparently, the smallest ∆q of the M06-2X functional suggests that the M11 and M06-2X functionals generate similar CT properties for the B-TCNE complexes.

Computational Method
The isolated monomers of the benzene and TCNE molecules were fully optimized at the level of the B3LYP functional and aug-cc-pVTZ basis sets [69] with a tight self-consistent field convergence threshold. The optimized benzene and TCNE molecules are planar structures with point groups D 6h and D 2h , respectively. The monomers of benzene and TCNE molecular geometries were frozen to construct the B-TCNE complexes because the intra-molecular monomer geometrical parameters are a little perturbed during the formation of donor-acceptor complexes [44,70]. Hence, the relative orientation of the π-stacked B-TCNE complexes was constructed in parallel cofacial configuration with the fixed distance between the molecular centers along the z-axis (see Figure 6, z-axis is also labeled). The typical donor-acceptor inter-planar distances were about 3.3-4.0 Å [46,[50][51][52]65,71]. In this work, the distance between the benzene and TCNE molecular centers is fixed at 3.5 Å, which is close to the value of the TTF-TCNQ [65] and terthiophene-TCNQ [71] complexes. The wavefunction-based methods EOM-CCSD [17], CIS/CIS(D) [18], and SAC-CI [16] were applied to calculate the excited states of the B-TCNE complexes. To compare the effects of the basis sets, the LanL2DZ and SVP were adopted for the EOM-CCSD calculations, while the SVP basis sets were applied in the CIS, CIS(D), and SAC-CI calculations. The EOM-CCSD, CIS, CIS(D), and SAC-CI calculations with aug-cc-pVTZ basis sets are beyond the capability of our computing cluster. In the full TDDFT and TDA-TDDFT calculations, the tested functionals included LC-ωPBE [37], M11 [55], ωB97XD [72], M06-2X [64], SOGGA11X [56], N12SX/MN12SX [57], HISSbPBE [58], and APF [59]. The LC-ωPBE, ωB97XD, and M06-2X functionals were selected for comparison of calculating singlet LE/CT and triplet excitations. The aug-cc-pVTZ basis sets were applied in the TDDFT calculations. To further check the effects of the basis sets, the full TDDFT calculations with the LC-ωPBE functional were also performed at the level of the SVP basis sets. In order to further check the LRCRSH functional performance in the CT and local excitations, the ωB97XD functional with the optimally tuned ω was also tested in the full TDDFT and TDA-TDDFT calculations. The optimally tuned ω value for the B-TCNE complex is 0.23 a 0 −1 , which was obtained according to the reported tuning procedure [20]. The geometrical optimization and excitation calculations were performed with the Gaussian09 package [73]. The transferred charges were analyzed using the Multiwfn program [74], and the details of the method for analyzing the transferred charges were reported in our previous work [75,76].

Conclusions
In this work, in order to get a balanced description of different excitation types for electron-donoracceptor systems, including singlet CT and LE, as well as triplet excited states, B-TCNE complexes were selected as the model system, and then several ab initio methods and the full/TDA-TDDFT were evaluated with different functionals. On the basis of the benchmark calculations of the EOM-CCSD method, the arithmetic mean of the absolute errors and standard errors of the electronic excitation energies for the different computational methods suggest that the accuracy of the SAC-CI method is better than that of the CIS/CIS(D) method, the M11 functional in DFT is superior to the other tested functionals for a balanced description of LE, CT, and triplet excitations, and the TDA-TDDFT improves the accuracy of the calculated excitation energies relative to full TDDFT. The TDA-TDDFT and full TDDFT results give similar transferred charges and CT excitation lengths. Meanwhile, the calculated errors are not systematic and depend on the characters of the excited states. Furthermore, in terms of the exchange-correlation energy formula for the M11 functional, the performance of the M11 functional underlines the importance of kinetic energy density, spin-density gradient, and range separation in the development of novel DFT functionals.