Performance Analysis and Optimization of Mixed-Reference Spin-Flip

Subscriber access provided by Macquarie University

A: Spectroscopy, Molecular Structure, and Quantum Chemistry

Performance Analysis and Optimization of Mixed-Reference SpinFlip Time-Dependent Density Functional Theory (MRSF-TDDFT) for Vertical Excitation Energies and Singlet-Triplet Energy Gaps Yevhen Horbatenko, Seunghoon Lee, Michael Filatov, and Cheol Ho Choi J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b07556 • Publication Date (Web): 22 Aug 2019 Downloaded from pubs.acs.org on August 29, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Performance Analysis and Optimization of Mixed-Reference Spin-Flip Time-Dependent Density Functional Theory (MRSFTDDFT) for Vertical Excitation Energies and Singlet-Triplet Energy Gaps Yevhen Horbatenko,†,§ Seunghoon Lee,‡,§ Michael Filatov,*,† Cheol Ho Choi*,†

†Department

of Chemistry, Kyungpook National University, Daegu, 41566, Republic of

Korea ‡Department

of Chemistry, Seoul National University, Seoul, 151747, Republic of Korea

Abstract The mixed-reference spin-flip (MRSF) time-dependent density functional theory (TDDFT) method eliminates the notorious spin-contamination of SF-TDDFT; thus, enabling identification of states of proper spin-symmetry for automatic geometry optimization and molecular dynamics simulations. Here, we analyze and optimize the MRSF-TDDFT in the calculations of the vertical excitation energies (VEEs) and the singlet-triplet (ST) gaps. The dependence of the obtained VEEs and ST gaps on the intrinsic parameters of the MRSFTDDFT method is investigated and prescriptions for the proper use of the method are formulated. For VEEs, MRSF-TDDFT displays similar or better accuracy than SF-TDDFT (ca. 0.5 eV), while considerably outperforming the LR-TDDFT for the ST gaps. As a result, a new functional of STG1X (dubbed here) especially for ST gaps is suggested on the basis of splitting between the components of the atomic multiplets.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction The linear response (LR) time-dependent density functional theory (TDDFT) is a popular computational method widely used to study excited states of large molecular systems.1-4 With the use of the currently available approximate density functionals, LR-TDDFT yields sufficiently accurate results for the low-lying excited states.5-7 However, it fails to describe several types of the excited states important for photochemical and photovoltaic (PV) applications. In particular, excitations of the ground state multi-reference systems,8-10 avoided crossings and conical intersections between the ground and excited state potential energy surfaces (PESs),11-13 long range charge transfer excitations14-16 are among the situations poorly described by LR-TDDFT based on the conventional single-determinant Kohn-Sham (KS) reference state. Some of these problems, e.g., description of the conical intersections, can be ameliorated with the use of SF-TDDFT, which employs a high-spin component of the triplet state as a reference.17-19 The use of a single component of the reference triplet state results in a substantial spin-contamination, which often impedes identification of the response excited state as a singlet or a triplet.20 To resolve this issue, a number of spin-adaptation (SA) techniques were developed in the context of SF-TDDFT21-24 and the related SF configuration interaction with single excitations (CIS)20,

25-27

methods. Perhaps, the most advanced SA

approach employs the tensor equations-of-motion (TEOM) formalism,21-24 which enables one to completely eliminate the erroneous spin-contamination of the SF-CIS excited states. However, due to complexity of the TEOM formalism and its incompatibility with approximate DFT functionals,24 implementation of the analytic energy gradients for the TEOM based SA-SF-TDDFT methods is still absent. This considerably impedes application of the TEOM method to practical problems.


Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Recently, some of us have developed a new method, dubbed MRSF-TDDFT,28 where the spin-contamination is ameliorated due to the use of both high-spin components of the triplet state, i.e., MS = +1 and MS = −1, as a reference for obtaining response excited states. Thus, a huge advantage of MRSF-TDDFT (further on denoted as MRSF for brevity) is that it recovers almost all the spin configurations necessary for restoring the proper total spin of the SF excited states. Another considerable advantage of this method is that the analytic energy gradients are straightforwardly available for the response states.29 Combined with the ease of identification of the excited states, e.g., as a singlet or a triplet, the availability of the analytic energy gradient considerably enhances application possibilities of the MRSF method and enables one to carry out automatic geometry optimization of targeted states and to run on-thefly non-adiabatic molecular dynamics simulations. So far, however, a proper benchmarking of the MRSF method against a database of VEEs standardly used to assess the accuracy of new computational techniques was lacking. Such information is needed, e.g., to evaluate applicability of the emerging methods to various types of excited state phenomena. Typically, the new methods, e.g., density functionals, or computational techniques, are evaluated for a broad class of excitations tabulated by Thiel et al.,30-32 which cover predominantly the singlet excited states. For certain applications, e.g., organic light emitting diodes (OLEDs) and PV systems exploiting the singlet exciton fission (SEF) phenomenon, the knowledge of the triplet excited states is also required. Hence, in this work, the MRSF method will be tested in the calculation of both types of excited states, singlets as well as triplets, and prediction of the method will be compared with the available theoretical best estimates (TBE). Special attention will be paid to description of the ST gaps by MRSF, which are important, e.g., for identifying suitable precursor compounds for the SEF PV systems. Towards this goal, the dependence of the MRSF excitation energies on the spin-pairing coupling which is an intrinsic parameter of the method (see below) will be 3 ACS Paragon Plus Environment


Page 4 of 35

studied, as well as the dependence on the choice of approximate density functional. It is believed that the results obtained in this work will assist in establishing the MRSF method as a simple yet accurate computational procedure applicable to a wide class of excited states.

2. MRSF Method In this Section, a brief account of the MRSF method will be given; detailed description can be found somewhere else.28-29 The MRSF method employs the zeroth-order reduced density matrix (RDM) for a mixed reference state, 1 M S 1 ( x, x)   0M S 1 ( x, x) , 0  2

0MR ( x, x) 

(1)

which comprises the RDMs of high-spin components, MS = +1 and MS = −1, of the triplet state. Within Tamm-Dancoff approximation (TDA), pole analysis of the dynamic polarizability with linear response of the density from this RDM results in eigenvalue equations for singlet (k = S) and triplet states (k = T),

A

( k )(0) pq , rs

( k )(0) X rs( k )(0)  ( k )(0) X pq

rs

,

(2)

( k )(0) ( k )(0) from which the excitation amplitudes X rs and X pq , and the excitation energies ( k )(0) ( k )(0) with respect to the reference state are obtained. In the above equation, Apq ,rs is the orbital

Hessian matrix,28 and the 0 superscript indicates that the amplitudes and excitation energies are obtained for the uncoupled components of the reference state. Spatial Kohn-Sham molecular orbitals for the MRSF and SF density are exactly the same.29 However, the excited states of MRSF-TDDFT are described by responses from the multireference density, while those of SF-TDDFT are from the single component (MS= +1) reference density. Thus, there is a clear advantage at response states. The response electronic configurations in Eq. (2) are shown in Figure 1. Compared to the conventional SF-TDDFT 4 ACS Paragon Plus Environment

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


formalism, the MRSF method generates electronic configurations, which are missing in SFTDDFT and which are important for restoring proper expectation value of the total spin operator Sˆ 2 . Intrinsic parameters of MRSF: Although MRSF restores the correct total spin of the response excited states, which considerably simplifies identification of the states of interest (triplet or singlet), the electronic configurations originating from different components of the mixed reference state remain uncoupled.28 To alleviate the situation, an a posteriori coupling between these configurations was introduced,28 much in the spirit of the spin-complete SFCIS method of Sears et al.20 Consequently, the coupling takes the form of: 1 ˆ M S 1 Apq ,rs  cSP  MpSq ,  H  r s 

(3)

M S 1 M S 1 where cSP is a spin-pairing coefficient and  p q and  r s are the configurations

generated from the

MS = +1 and MS = −1 components of the mixed reference state,

respectively. The spin-coupled excited states are then obtained from Eq. (4) by replacing the ( k )(0) MRSF orbital Hessian Apq ,rs by the matrix:

( k ) A(pqk ),rs  A(pqk )(0) , rs  Apq , rs .

(4)

The use of the MRSF formalism helped to mitigate erroneous splitting between the different components, e.g., iPx,y and iPz (i=1,3), of the atomic and molecular multiplet states that was predicted by the usual SF-TDDFT method.28 However, there is no exact prescription as to the choice of the cSP coefficient; in Ref. [28] it was initially set to be equal to the amount of the exact exchange cHF in the approximate density functional used. The cHF is a parameter that defines the fraction of the exact or Hartree-Fock (HF) exchange energy in an approximate global hybrid exchange-correlation (XC) functional, such as B3LYP, PBE0, or BHHLYP via33



E XC  cHF E XHF  (1  cHF ) E XDFT  ECDFT ,

Page 6 of 35

(5)

HF where E X is the HF exchange energy. In the present work, the MRSF method is used in

connection with global hybrid functionals only, where cHF is a constant; hence, there is no ambiguity with identifying cHF. In the future, the method will be extended for the use in connection with local hybrid and range separated functionals, where the contribution of the exact exchange is given by a function of the interelectronic distance, rather than a constant.

3. Computational details All single point calculations were performed using locally modified GAMESS-US program package.34 Two all-electron basis sets of 6-31G*35 and cc-pVTZ36 were employed. Vertical excitation energies: The accuracy of description of the valence excitation energies was assessed in the calculation of the benchmark set of medium-sized molecules as introduced by Thiel and co-workers.30-32 This set contains the most common organic chromophores such as polyenes, polyaromatic hydrocarbons, carbonyl compounds, amides, heterocycles, and nucleobases and is widely used to benchmark existing and new computational methods. In this work, the latest version of the aforementioned set (TBE-2) was used, where the aug-cc-pVTZ basis was employed instead of the TZVP basis.32 The LRTDDFT, SF-TDDFT, and MRSF-TDDFT methods in their collinear formulation were used in the calculations. The hybrid DFT functionals with various amount of exact exchange were employed (the amount of the exact exchange is given in parentheses), i.e., B3LYP (0.21),37-38 PBE0 (0.25),39 BHHLYP (0.5),33 M08-SO (0.5679).40 B3LYP was chosen as it is a generally accepted standard and is believed to have a good average performance over many chemical problems.5,

41-42

The PBE0 functional was included owing to its good performance for the

singlet excitation energies observed in the previous works.43-45 BHHLYP was chosen since it 6 ACS Paragon Plus Environment

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


performs well for calculations of the vibronic structure in electronic spectra of organic molecules.46 In this work, only the VEEs were calculated. The dependence of VEEs on the spin-pairing coefficient cSP was studied by varying it in the range of 0.3–0.9. All the calculated values of the vertical excitation energies were compared against the TBE-2 reference data. The accuracy of the calculations was estimated by means of the mean absolute error (MAE): MAE 

1 E TBE  2  E calc ,  n n

(6)

where n is the number of the singlet exited states, ETBE-2 and Ecalc are the TBE-2 reference and calculated energy, respectively. Degeneracy of atomic multiplet states: To study degeneracy of various multiplet state components, the MRSF calculations for a set of neutral and ionized atoms of the first and second row elements, i.e., Be, B+, C, N+, O, F+, Mg, Al+, Si, P+, S, Cl+, were carried out. For this type of calculations, a hybrid DFT functional based on the Becke exchange47 and the Lee-Yang-Parr correlation48 functionals and employing a varying fraction of the exact exchange B(cHF)HF(1−cHF)-LYP was used (the amount of the exact exchange is given in parentheses): BHHLYP = B(0.5)HF(0.5)-LYP, B(0.2)HF(0.8)-LYP, B(0.15)HF(0.85)-LYP, B(0.1)HF(0.9)-LYP, B(0.0)HF(1.0)-LYP. Simultaneously, a value of the spin-pairing coefficient cSP was varied in the range of 0.1–0.9 to get insight into of its effect on the degeneracy. Singlet-triplet gaps: For the calculations of the ST gaps, the aforementioned atomic set was augmented by two additional molecular sets. The first molecular set was taken from the work of Zimmerman et al.49 and it includes small and medium-sized molecules. This benchmark set, probably, is a good option for testing the ST gaps, as it comprises molecules that have either diradical or diradicaloid electronic characteristic in their electronic ground states. The 7 ACS Paragon Plus Environment


Page 8 of 35

second set was adopted from a recent work of Van Voorhis,50 where various kinds of polycyclic aromatic hydrocarbons,51-55 new organic PV materials,56-59 and thermally activated delayed fluorescence emitters were studied.60-62 The majority of the molecules in the set are large-size organic molecules and are used or potentially can be used in various real-life applications. The calculations with the set of Zimmerman were performed with the B3LYP, PBE0, BHHLYP, and B(0.15)HF(0.85)-LYP density functionals. In addition, the M08-HX (0.5223)40 was employed due to its good performance for the ST gaps of a benchmark set introduced by Truhlar.63 These computations were carried out by using the MRSF and LR method. For the calculations with the set of Van Voorhis, only B(0.15)HF(0.85)-LYP functional was employed. The computations of the ST gaps for this set were carried out using the MRSF formalism only. The results were compared with those from Ref.

50

obtained by

the LR method. For both types of calculations, i.e., with the Zimmerman and van Voorhis sets, the cSP value was set to be equal to the cHF value of the approximate DFT functional. All the calculated values of the energies of the atomic multiplet state components as well as the molecular ST gaps were compared against available experimental data. The accuracy of these calculations was estimated by means of the absolute error (AE): AE  E exp  E calc ,

(7)

where Eexp and Ecalc are experimental and calculated energy, respectively.

4. Results and Discussion 4.1. Valence excited states. When calculating the singlet valence excited states, diffuse functions were not employed; this choice of the basis set is suitable for prediction of the energies of the low-lying excited states, see, e.g., Thiel et al.30 Dependence on spin-pairing coupling: VEEs of the benchmark set of molecules were calculated with the density functionals described in the previous section, see Figures S3 and 8 ACS Paragon Plus Environment

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


S4 of the Supporting Information. To investigate the effect of cSP on the calculated VEEs, the MRSF calculations with cSP varying in the range 0.3 – 0.9 were carried out using the BHHLYP functional. The calculated VEEs display only marginal dependence on the choice of the cSP constant. Figure 2a shows correlation of the calculated VEEs with the TBE-2 values for two limiting cSP values, 0.3 and 0.9. Hence, in practical application of the MRSF method, any cSP value within the studied range can be chosen, e.g., a value equal to the fraction of the exact exchange in the approximate XC functional. Dependence on basis set: Although our main goal is to study the dependence of VEE on cSP, we calculated VEE with two different basis sets (double- and triple-zeta qualities) to check the basis set dependencies. Correlation plot for vertical excitation energies calculated with BHHLYP and two basis sets, 6-31G* and cc-pVTZ, is shown in Figure 2b. The VEEs calculated with cSP= 0.3 and 0.9 are very close to each other and show only a weak dependence on the cSP value. For both values of cSP, the excitation energies somewhat decrease upon enlarging the basis set. The similar effect of decreasing of the excitation energies was reported in the works of Thiel et al.30,

32

Generally, correlation between the

VEEs obtained using the two basis sets is excellent, with a correlation coefficient ca. 0.99, see Figure 2b. For correlation plots obtained using the other functionals, see Figure S1 of the Supporting Information. Somewhat larger deviations of the VEEs obtained with the 6-31G* basis from the cc-pVTZ ones was observed only for the energies above 6.5 eV; this can be attributed to the occurrence of the mixed valence-Rydberg states.30,

32

Hence, the MRSF-

TDDFT/6-31G* method is expected to produce sufficiently reliable excitation energies and can be confidently used for automatic optimization of the geometry of the excited states and the non-adiabatic molecular dynamics simulations. Dependence on the fraction of the exact exchange: Statistical performance of various TDDFT-based methods (MRSF, SF, and LR) for VEEs of the singlet excited states obtained 9 ACS Paragon Plus Environment


using 6-31G* basis set is given in Figure 3 and Table 1 (for the calculations employing the cc-pVTZ basis set, see Figure S2 and Table S1 in the Supporting Information). The MRSFTDDFT VEEs obtained with various functionals (and cSP = cHF) follow the same trend as was observed earlier for density functionals with varying percentage of the HF exchange,42 i.e., the hybrid functionals with a moderate percentage of the exact exchange (22-25%), such as PBE0, provide a better agreement with the reference data, whereas the functionals with a higher percentage of the exact exchange (>50%), such as BHHLYP and M08-SO, tend to overestimate the excitation energies. The obtained MAEs of the MRSF VEEs from the TBE2 reference values is 0.46 eV (PBE0), 0.52 eV (B3LYP), and 0.57 eV (BHHLYP and M08SO) are to be compared to the SF-TDDFT MAEs of 0.52 eV (PBE0), 0.58 eV (B3LYP), 0.57 eV (BHHLYP), and 0.54 eV (M08-SO), and to the LR-TDDFT MAEs of 0.34 eV (PBE0), 0.30 eV (B3LYP), 0.71 eV (BHHLYP), and 0.48 eV (M08-SO). It is noteworthy that both MRSF and SF-TDDFT outperform LR-TDDFT when using the BHHLYP functional, whereas the latter method favors functionals with small percentage of the HF exchange, such as B3LYP and PBE0. Figure 3 shows also variances for each functional and TDDFT method. The variances were calculated as squares of the corresponding standard deviation values. In the case of the MRSF method they are 0.14, 0.15, 0.23, and 0.23 eV2 for PBE0, B3LYP, BHHLYP, and M08-SO, respectively. The SF method delivers somewhat larger variance values, i.e., 0.17, 0.19, 0.21, 0.28 eV2 in the same series of DFT functionals, whereas the smallest variances are found for the LR method, i.e., 0.07, 0.07, 0.13, and 0.11 eV2. Correlation plots between calculated VEEs and TBE-2 values, as well as the respective Pearson correlation coefficient r values, are shown in Figures S3 and S4 for the 6-31G* and cc-pVTZ basis sets, respectively. As seen in Figure S3, MRSF-TDDFT slightly outperforms the conventional SF-TDDFT method; the correlation coefficient averaged over all the DFT functionals is 0.894, for the former, and 0.869, for the latter. Both methods however are 10 ACS Paragon Plus Environment

Page 10 of 35

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


outperformed by the standard LR-TDDFT which yields average value of 0.956. The VEEs produced by LR-TDDFT are narrowly distributed around the ideal correlation line, whereas the MRSF- and SF-TDDFT VEEs show somewhat larger deviation from the TBE-2 values. There are cases where the calculated MRSF excitation energies, especially for highlying valence excited states, deviate significantly from the TBE-2 reference. The most probable reason for these deviations is either valence-Rydberg mixing or contribution from double excitations. This occurs in the case of such molecules as formaldehyde, acetone, pyrrole, pyrazine, pyridine, naphthalene, uracil, thymine, adenine. For example, in formaldehyde and acetone there are    * and    * excitations around 9 eV. It is known that both states have strong valence-Rydberg mixing,64 and a basis set with diffuse functions is needed to reproduce these excitations. Moreover,    * excited state is characterized by a configurational mixing with the ground state as well as with    * double excitations.64 For naphthalene, there is a    * excited state around 6.5 eV that has strong contribution from doubly excited configurations.65 For the MRSF method, the variation of VEEs and MAEs with respect to the cSP coefficient is presented in Table 1. With the increasing cSP value (0.3 – 0.9), the MAE (and VEE) of MRSF increase non-linearly; however, the overall variation is small. Note that, when using cSP=0, the MRSF method recovers the results obtained with SF-TDDFT.28 At an intermediate value of the cSP coefficient, MAE of MRSF passes through a shallow minimum. The cSP value at the minimum is dependent on the functional employed; with the B3LYP, PBE0, and BHHLYP functionals, the minimum occurs within the 0.4-0.5, 0.3-0.5, and 0.00.3 range of the

cSP values, respectively. However, the overall decrease of MAE is

insignificant and no attempt to optimize the cSP value was done.



4.2. Atomic multiplet states. One of the typical problems encountered by TDDFT formalism in general is related to erroneous splitting between the components of the excited state multiplets.66-69 As shown in Ref. [28], MRSF reduces considerably the magnitude of the erroneous splitting, introduced by SF-TDDFT. Nearly perfectly degenerate iPx,y and iPz, (i=1, 3) components of the singlet and triplet atomic states were obtained with the use of non-zero spin-pairing coefficient cSP. Although an insignificant splitting between the multiplet components still remains due to the self-interaction error of the approximate density functionals (the splitting disappears, when the HF exchange-only functional is used), the results of Ref. [28] encourage a more systematic study into the effect of the cSP magnitude on the multiplet splitting. As shown in Figure 4, the MRSF response states originating from both components of the mixed reference state span markedly larger set of electronic configurations than could be obtained from the closed-shell reference state of LR-TDDFT. In some cases, e.g, Be and B+, the MRSF set of response configurations is sufficient to represent the target electronic state accurately; in the other cases, only a few electronic configurations are missing in the MRSF set. By contrast, LR-TDDFT generates much fewer response configurations, in most cases not nearly sufficient to represent some multiplet states, e.g., the 1D and 1S states of C, N+, O, and F+. Dependence of the multiplet splitting on spin-pairing constant: Figure 5a shows the energies of the components of the 3P, 1D, and 1S multiplet states of carbon in dependence on the cSP value. The components of the multiplet states can be separated into two types: For the former, the energies are independent on the cSP values, and, for the latter, there is a pronounced dependence on this value. The first type of multiplet components originates from the OO transitions from the reference state, and the second is due to the OV transitions, see Figures 1 and 4. For the latter response states, the spin-pairing term in Eqs. (3) and (4) couples the 12 ACS Paragon Plus Environment

Page 12 of 35

Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


configurations obtained from the different components, MS = +1 and MS = −1, of the mixed reference; this term includes only the exchange integral (OV|VO) and contribution of the DFT exchange-correlation kernel vanishes.28 To provide for the degeneracy between the multiplet components obtained due to the two types of transitions, the magnitude of the spinpairing term needs to be the same as the magnitude of the coupling between the electronic configurations originating from the OO transitions. The latter is dominated by the exchange integral between the open-shell orbitals of the reference state. Given that, for carbon, only two of the three valence p-orbitals are occupied in the reference state – hence, belong to the OO subspace – and the third p-orbital is vacant, an approximate equalization between the energies of the OO and OV response states belonging to the same multiplet can be achieved by setting cSP ≈ cHF and requiring that the fraction of the HF exchange cHF be sufficiently large. The validity of this conjecture is tested in the atomic calculations reported in the following. Dependence of the multiplet splitting on fraction of the exact exchange: As the full spherical symmetry is not available in the GAMESS-US code, in the following, the atomic states will be classified by the symmetry species of the D2h point group. For carbon, the 1Ag, 1B1g, 1Ag (different 1Ag), 3B1g states originate from the OO transitions and the

1B

2g,

1B

3g,

3B

2g,

3B

3g

states are due to the OV transitions from the mixed reference state. Analysis of the atomic orbital populations in these states lets us make the following assignment: The 3B1g, 3B2g, and 3B

3g

belong to the 3P atomic multiplet, one of the 1Ag states and the 1B1g, 1B2g, 1B3g states

belong to the 1D atomic multiplet, and the remaining 1Ag state represents the 1S multiplet. Some of the multiplet components are missing due to the incompleteness of the MRSF excitation space, see Figures 1 and 4. Ideally, the components of the 3P multiplet and the 1D multiplet states need to be degenerate among each other. Figure 5 shows the dependence of the energies of the singlet 13 ACS Paragon Plus Environment


Page 14 of 35

and triplet response states of carbon on the cSP value for several cHF values. A nearly perfect degeneracy between the components of the 3P multiplet and the 1D multiplet is restored by setting cSP = cHF = 0.8, see Figure 5b. The relative energies of the corresponding states (either singlets or triplets) generated within the OO and OV spaces change with the amount of the exact exchange cHF. To understand this behavior, the dependence of the O and V orbital energies on cHF was investigated, see Figure S5 of the Supporting Information. With the increasing fraction of the exact exchange, the O orbital energies remain almost constant, whereas the energy of the V-type orbital increases. This shifts the center of gravity of the configurations originating from the OV transitions upward in energy. Simultaneously, the splitting between the components of the 3P and 1D multiplets originating in the OO space is increasing with the increasing fraction of the exact exchange. At a sufficiently large cHF value, the two effects offset each other and a nearly perfect degeneracy between the components of the multiplet states is reached. Besides restoring degeneracy between the components of the multiplet states of carbon, the use of cSP = cHF ≈ 0.8 results in very good agreement between the calculated and experimental 3P-1D and 3P-1S multiplet splittings; 1.28 (calc.) and 1.26 (exp.70) eV for the former and 2.68 (calc.) and 2.68 (exp.70) eV for the latter. Hence, it can be recommended to use cSP = cHF ~ 0.8 to obtain both nearly degenerate components of the multiplet states and accurate values of the multiplet splittings of other atoms. Optimizing intrinsic parameters of MRSF: The observation that the atomic multiplet states can be recovered by a suitable choice of the cSP and cHF parameters suggests that these parameters can be optimized for a set of atoms and atomic ions. Using the non-degeneracy error ΔEND (Eq.8) E ND   ( SOO  SOV )2  (TOO  TOV )2  J


J

,

(8)

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where XOO and XOV (X = S, T) are the response states (singlets and triplets) obtained by excitations within the OO subspace and the OV subspace, respectively, as the target function the optimal values of the cSP and cHF parameters were obtained for a set of the first and second row atoms and atomic ions, J = Be, B+, C, N+, O, F+, Mg, Al+, Si, P+, S, Cl+. Setting cSP = cHF, the optimal values of cHF were found for a functional of the form B(1−cHF)HF(cHF)LYP in connection with the 6-31G* and cc-pVTZ basis sets. For both basis sets, ΔEND reaches minimum at very close cHF values, 0.850 and 0.848, respectively (Figure S6, the Supporting Information). Hence, cHF = 0.85 is recommended in this work as a basis set independent value. In the following, the functional B(1−cHF)HF(cHF)-LYP with cHF = 0.85 is referred to as STG1X. The relative energies of the 3P, 1D, and 1S multiplets calculated using the 6-31G* basis set for the set of atoms and atomic ions are given in Table S3 of the Supporting Information. For each multiplet state, the energy was calculated as an average of its components; for instance, the energy of the 3P state of carbon, see Figure 5a, was calculated as an average of the 3B1g and 3B2g components. The magnitudes of the splitting between the atomic multiplet states, i.e., 3P − 1S, 1P − 1S, 1D − 3P, and 1S − 3P were obtained as differences of the respective average values. Figure 6 shows the splitting values calculated with the STG1X and the BHHLYP functionals, which are compared against the experimental data (see Figure S7 for absolute deviations from the experimental data). Except the 3P − 1S splitting (Figure 6a left panel), STG1X clearly outperforms BHHLYP. Overall, the STG1X results are in a very good agreement with the experimental values. In the case of the 1D − 3P splitting, STG1X almost reaches “chemical accuracy” of 0.05 eV according to its MAE of 0.08 eV (Figure 6b). The accuracy improvement by STG1X is remarkable, considering the fact that STG1X was not optimized to predict accurate multiplet splittings, but to reduce the degeneracy errors. This observation underlines an internal consistency of the method.



4.3. Molecular ST gaps. Performance of the STG1X functional was further tested in the MRSF calculations of the singlet-triplet (ST) gaps in a series of molecules. Figure 7a shows the results for a series of small and medium size molecules. Besides the absolute errors in the MRSF/STG1X ST gaps Figure 7a reports the results of the calculations obtained with several other functionals; B3LYP, PBE0, BHHLYP and M08-HX. The adiabatic ST gaps were obtained by taking energies of the singlet and the triplet states at the geometries optimized for the respective state and the vertical ST gaps were calculated at the geometry of the lowest energy state, singlet or triplet. The calculated MRSF values for the ST gaps along with their signs are collected in Table S4 of the Supporting Information. For molecules NH to PH2+, the best agreement with the experimental ST gaps is obtained with the use of the STG1X functional. The use of lower fraction of the HF exchange resulted in greater errors and even in prediction of the wrong ground state of CH2; the singlet 1A

1

state was wrongly predicted to be below the triplet 3B1 state by the B3LYP, PBE0, M08-

HX and BHHLYP functionals. The same occurred for the ground state of NH2+ in the case of B3LYP and PBE0. Similar findings were previously reported for the SF-TDDFT calculations with some of these functionals.71 For molecules of medium size, it was the M08-HX and BHHLYP functionals that showed the best agreement with the experimental ST gaps. However, the STG1X functional displays sufficiently good accuracy as well; there are only a few cases where the deviations from the experimental ST gaps exceed 0.36 eV (a value suggested by Truhlar et al.63 to classify a functional as accurate or not), in particular, o- and m-benzyne. However, the ground singlet states of the latter molecules possess only a little diradical characteristic (as can be judged by the occupation numbers of the natural orbitals shown in Figure S8(c,d)) and their relative stability with respect to the triplets is underestimated by the functionals with large fraction of the HF exchange. Nevertheless, the MAE of the STG1X functional (0.16 16 ACS Paragon Plus Environment

Page 16 of 35

Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


eV) for the full set of molecules reported in Table S4 is lower than that for the M08-HX (0.34 eV) and BHHLYP (0.37 eV) functionals. The B3LYP and PBE0 functionals showed the worst performance with MAE 0.79 and 0.73 eV, respectively. The ST gaps were also calculated using the LR method for the same set of molecules, see Figure 7b for the absolute errors and Table S5 of the Supporting Information for the ST gap values and their signs. For smaller molecules NH to PH2+, the best agreement with the experiment is achieved by the M08-HX functional. For the medium-size molecules all functionals yield deviations from the experiment considerably exceeding 0.36 eV. All the functionals tested wrongly predict the triplet state to be the ground state for c-C4H4 and pbenzyne instead of singlet. The total smallest MAE is observed for B3LYP (0.55 eV) and M08-HX (0.58 eV) followed by PBE0 (0.71 eV) and BHHLYP (0.84 eV). Overall, MAE increases as the amount of the HF exchange increases. Thus, in connection with LR-TDDFT, STG1X is expected to yield large value of MAE since this functional was mainly optimized for the MRSF method only. As the original SF method introduces the spin contamination, which impedes identification of the singlet and triplet states, the SF method was not used here to calculate the ST gaps. For the case of the Be atom, the occurrence of the erroneous spin contamination was already shown in the previous study.28 The molecules reported in Figure 8 and Table S6 (Supporting Information) were selected from the molecular set used by Van Voorhis et al.50 for benchmarking of their triplettuned (TT) functional TT-PBEh; only the molecules, for which the largest deviations from the experimental ST gaps were obtained in Ref. [50], were selected here. Figure 10 shows the absolute errors between the ST gaps calculated with MRSF/STG1X and the experimental values. For comparison, the absolute errors for the same molecules calculated with LRTDDFT/TT-PBEh50 are provided as well. On average, the STG1X functional displays accuracy comparable with the TT-PBEh functional. The MAE of STG1X is 0.38 eV, only 17 ACS Paragon Plus Environment


insignificantly worse than the MAE of the TT functional (0.29 eV). Considering, that the value of the cHF parameter of the STG1X functional was not optimized for the ST gaps (it was obtained to mitigate the non-degeneracy error of the atomic multiplets), the good performance of the latter functional in comparison with an individually tuned functional is gratifying.

5. Conclusions In this work, the accuracy of the vertical excitation energies (VEE) and of the singlet-triplet (ST) energy gaps obtained with the use of the MRSF method was evaluated by comparing with the available theoretical best estimates. The dependence of the calculated VEEs and ST gaps on the intrinsic parameters of the MRSF method (the spin-pairing coefficient, cSP) and the parameters of the approximate exchange-correlation functional (percentage of the HF exchange, cHF) was investigated. For VEEs, the standard benchmark set of Thiel et al.30 was used. A number of density functionals with varying fraction of the HF exchange, B3LYP, PBE0, BHHLYP, and M08SO, were employed in connection with MRSF and two basis sets, 6-31G* and cc-pVTZ. For comparison, the standard LR-TDDFT and SF-TDDFT calculations were performed. The VEEs obtained with MRSF display a weak dependence on the basis set size; similar with the other computational techniques. Overall, MRSF yields more accurate VEEs than SF-TDDFT and is somewhat less accurate that LR-TDDFT; the best MAEs are 0.46 eV (MRSF/PBE0), 0.52 eV (SF/PBE0), and 0.30 eV (LR/B3LYP). Only an insignificant dependence of the VEEs calculated with MRSF on the spin-pairing parameter cSP was found. Hence, cSP can be chosen equal to the amount of the HF exchange cHF in the approximate density functional, which was our initial prescription.


Page 18 of 35

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In contrast to VEE, it was observed that the energies of the singlet and triplet states depend strongly on cSP. MRSF calculations of the singlet and triplet states were carried out for a number of atoms, atomic ions, and molecules of varying size. For atoms and atomic ions, a special emphasis was put on the elimination of the erroneous non-degeneracy of the response states belonging to the same atomic multiplet state. On the basis of these calculations, the latter parameter is recommended to be chosen equal to the HF exchange fraction, i.e., cSP = cHF. Thus, setting cSP = cHF = 0.85 mitigates the non-degeneracy error and considerably improves agreement of the calculated multiplet energies with the experiment. Consequently, the STG1X = B(0.15)HF(0.85)-LYP functional proposed on the basis of the atomic calculations was used in connection with MRSF (cSP = 0.85) to calculate the ST gaps in a series of molecules of varying size. Although for some molecules, especially of medium and large size, better ST gaps were obtained with a functional featuring lower fraction of the HF exchange (cHF ≈ 0.5), the STG1X functional showed sufficiently good overall performance. Thus, the latter functional describes the ST gaps of molecules representing interest for PV and biological applications with an accuracy (MAE = 0.38 eV) comparable with the recently proposed triplet-tuned TT-PBEh functional of Van Voorhis et al.50 (MAE = 0.29 eV). Given that the intrinsic parameters of the MRSF/STG1X method were never adjusted to reproduce the experimental ST gaps, the good overall performance for the ST gaps is gratifying and the MRSF/STG1X method can be suggested for computation of triplet excited states.

Associated Content Supporting Information Figures S1–S8 and Tables S1–S6 (PDF).



Author Information Corresponding Author *E-mail: [email protected] Tel: +82-53-950-5332. *E-mail: [email protected] Tel: +82-53-950-5332. Author Contribution §Y.H.

and S.L. contributed equally to this work.

Notes The authors declare no competing financial interest.

Acknowledgements This work was supported by the Samsung Science and Technology Foundation (SSTFBA1701-12).


Page 20 of 35

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References

1. Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 1984, 52, 997-1000. 2. Burke, K.; Werschnik, J.; Gross, E. K. U. Time-dependent density functional theory: Past, present, and future. J. Chem. Phys. 2005, 123, 062206. 3. Casida, M. E.; Huix-Rotllant, M. Progress in Time-Dependent Density-Functional Theory. Annu. Rev. Phys. Chem. 2012, 63, 287-323. 4. Maitra, N. T. Perspective: Fundamental aspects of time-dependent density functional theory. J. Chem. Phys. 2016, 144, 220901. 5. Matsuzawa, N. N.; Ishitani, A.; Dixon, D. A.; Uda, T. Time-Dependent Density Functional Theory Calculations of Photoabsorption Spectra in the Vacuum Ultraviolet Region. J. Phys. Chem. A 2001, 105, 4953-4962. 6. Zhao, Y.; Truhlar, D. G. Density Functional for Spectroscopy: No Long-Range Self-Interaction Error, Good Performance for Rydberg and Charge-Transfer States, and Better Performance on Average than B3LYP for Ground States. J. Phys. Chem. A 2006, 110, 13126-13130. 7. Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215-241. 8. Cave, R. J.; Zhang, F.; Maitra, N. T.; Burke, K. A dressed TDDFT treatment of the 21Ag states of butadiene and hexatriene. Chem. Phys. Lett. 2004, 389, 39-42. 9. Maitra, N. T.; Zhang, F.; Cave, R. J.; Burke, K. Double excitations within time-dependent density functional theory linear response. J. Chem. Phys. 2004, 120, 5932-5937. 10. Neugebauer, J.; Baerends, E. J.; Nooijen, M. Vibronic coupling and double excitations in linear response time-dependent density functional calculations: Dipole-allowed states of N2. J. Chem. Phys. 2004, 121, 6155-6166. 11. Levine, B. G.; Ko, C.; Quenneville, J.; MartÍnez, T. J. Conical intersections and double excitations in time-dependent density functional theory. Mol. Phys. 2006, 104, 1039-1051. 12. Huix-Rotllant, M.; Filatov, M.; Gozem, S.; Schapiro, I.; Olivucci, M.; Ferré, N. Assessment of Density Functional Theory for Describing the Correlation Effects on the Ground and Excited State Potential Energy Surfaces of a Retinal Chromophore Model. J. Chem. Theory Comput. 2013, 9, 39173932. 13. Gozem, S.; Melaccio, F.; Valentini, A.; Filatov, M.; Huix-Rotllant, M.; Ferré, N.; Frutos, L. M.; Angeli, C.; Krylov, A. I.; Granovsky, A. A.; Lindh, R.; Olivucci, M. Shape of Multireference, Equation-ofMotion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection. J. Chem. Theory Comput. 2014, 10, 3074-3084. 14. Dreuw, A.; Weisman, J. L.; Head-Gordon, M. Long-range charge-transfer excited states in time-dependent density functional theory require non-local exchange. J. Chem. Phys. 2003, 119, 2943-2946. 15. Dreuw, A.; Head-Gordon, M. Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules. Chem. Rev. 2005, 105, 4009-4037. 16. Maitra, N. T. Undoing static correlation: Long-range charge transfer in time-dependent density-functional theory. J. Chem. Phys. 2005, 122, 234104. 17. Shao, Y.; Head-Gordon, M.; Krylov, A. I. The spin–flip approach within time-dependent density functional theory: Theory and applications to diradicals. J. Chem. Phys. 2003, 118, 4807-4818.



18. Rinkevicius, Z.; Vahtras, O.; Ågren, H. Spin-flip time dependent density functional theory applied to excited states with single, double, or mixed electron excitation character. J. Chem. Phys. 2010, 133, 114104. 19. Li, Z.; Liu, W. Theoretical and numerical assessments of spin-flip time-dependent density functional theory. J. Chem. Phys. 2012, 136, 024107. 20. Sears, J. S.; Sherrill, C. D.; Krylov, A. I. A spin-complete version of the spin-flip approach to bond breaking: What is the impact of obtaining spin eigenfunctions? J. Chem. Phys. 2003, 118, 90849094. 21. Li, Z.; Liu, W. Spin-adapted open-shell random phase approximation and time-dependent density functional theory. I. Theory. J. Chem. Phys. 2010, 133, 064106. 22. Li, Z.; Liu, W.; Zhang, Y.; Suo, B. Spin-adapted open-shell time-dependent density functional theory. II. Theory and pilot application. J. Chem. Phys. 2011, 134, 134101. 23. Li, Z.; Liu, W. Spin-adapted open-shell time-dependent density functional theory. III. An even better and simpler formulation. J. Chem. Phys. 2011, 135, 194106. 24. Zhang, X.; Herbert, J. M. Spin-flip, tensor equation-of-motion configuration interaction with a density-functional correction: A spin-complete method for exploring excited-state potential energy surfaces. J. Chem. Phys. 2015, 143, 234107. 25. Casanova, D.; Head-Gordon, M. The spin-flip extended single excitation configuration interaction method. J. Chem. Phys. 2008, 129, 064104. 26. Tsuchimochi, T. Spin-flip configuration interaction singles with exact spin-projection: Theory and applications to strongly correlated systems. J. Chem. Phys. 2015, 143, 144114. 27. Mato, J.; Gordon, M. S. A general spin-complete spin-flip configuration interaction method. Phys. Chem. Chem. Phys. 2018, 20, 2615-2626. 28. Lee, S.; Filatov, M.; Lee, S.; Choi, C. H. Eliminating spin-contamination of spin-flip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix. J. Chem. Phys. 2018, 149, 104101. 29. Lee, S.; Kim, E. E.; Nakata, H.; Lee, S.; Choi, C. H. Efficient implementations of analytic energy gradient for mixed-reference spin-flip time-dependent density functional theory (MRSF-TDDFT). J. Chem. Phys. 2019, 150, 184111. 30. Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3. J. Chem. Phys. 2008, 128, 134110. 31. Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. Benchmarks for electronically excited states: Time-dependent density functional theory and density functional theory based multireference configuration interaction. J. Chem. Phys. 2008, 129, 104103. 32. Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. Benchmarks of electronically excited states: Basis set effects on CASPT2 results. J. Chem. Phys. 2010, 133, 174318. 33. Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. J. Chem. Phys. 1993, 98, 1372-1377. 34. Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery Jr, J. A. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347-1363. 35. Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chim. Acta 1973, 28, 213-222. 36. Jr., T. H. D. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007-1023. 37. Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652. 38. Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623-11627. 22 ACS Paragon Plus Environment

Page 22 of 35

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


39. Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158-6170. 40. Zhao, Y.; Truhlar, D. G. Exploring the Limit of Accuracy of the Global Hybrid Meta Density Functional for Main-Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2008, 4, 1849-1868. 41. Bauernschmitt, R.; Ahlrichs, R. Treatment of electronic excitations within the adiabatic approximation of time dependent density functional theory. Chem. Phys. Lett. 1996, 256, 454-464. 42. Jacquemin, D.; Wathelet, V.; Perpète, E. A.; Adamo, C. Extensive TD-DFT Benchmark: SingletExcited States of Organic Molecules. J. Chem. Theory Comput. 2009, 5, 2420-2435. 43. Adamo, C.; Scuseria, G. E.; Barone, V. Accurate excitation energies from time-dependent density functional theory: Assessing the PBE0 model. J. Chem. Phys. 1999, 111, 2889-2899. 44. List, N. H.; Olsen, J. M.; Rocha-Rinza, T.; Christiansen, O.; Kongsted, J. Performance of popular XC-functionals for the description of excitation energies in GFP-like chromophore models. Int. J. Quantum Chem. 2012, 112, 789-800. 45. Leang, S. S.; Zahariev, F.; Gordon, M. S. Benchmarking the performance of time-dependent density functional methods. J. Chem. Phys. 2012, 136, 104101. 46. Dierksen, M.; Grimme, S. Density functional calculations of the vibronic structure of electronic absorption spectra. J. Chem. Phys. 2004, 120, 3544-3554. 47. Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At. Mol. Opt. Phys. 1988, 38, 3098-3100. 48. Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter 1988, 37, 785-789. 49. Zimmerman, P. M. Singlet–Triplet Gaps through Incremental Full Configuration Interaction. J. Phys. Chem. A 2017, 121, 4712-4720. 50. Lin, Z.; Van Voorhis, T. Triplet Tuning: A Novel Family of Non-Empirical Exchange–Correlation Functionals. J. Chem. Theory Comput. 2019, 15, 1226-1241. 51. Schmidt, W. Photoelectron spectra of polynuclear aromatics. V. Correlations with ultraviolet absorption spectra in the catacondensed series. J. Chem. Phys. 1977, 66, 828-845. 52. McClure, D. S. Excited Triplet States of Some Polyatomic Molecules. I. J. Chem. Phys. 1951, 19, 670-675. 53. Clar, E. Absorption spectra of aromatic hydrocarbons at low temperatures. LV-Aromatic hydrocarbons. Spectrochim. Acta 1950, 4, 116-121. 54. Bolovinos, A.; Tsekeris, P.; Philis, J.; Pantos, E.; Andritsopoulos, G. Absolute vacuum ultraviolet absorption spectra of some gaseous azabenzenes. J. Mol. Spectrosc. 1984, 103, 240-256. 55. Eland, J. H. D. Photoelectron spectra and ionization potentials of aromatic hydrocarbons. Int. J. Mass Spectrom. Ion Phys. 1972, 9, 214-219. 56. Siegert, S.; Vogeler, F.; Marian, C. M.; Weinkauf, R. Throwing light on dark states of αoligothiophenes of chain lengths 2 to 6: radical anion photoelectron spectroscopy and excited-state theory. Phys. Chem. Chem. Phys. 2011, 13, 10350-10363. 57. Huang, D.-L.; Dau, P. D.; Liu, H.-T.; Wang, L.-S. High-resolution photoelectron imaging of cold C60− anions and accurate determination of the electron affinity of C60. J. Chem. Phys. 2014, 140, 224315. 58. Pasinszki, T.; Krebsz, M.; Vass, G. Ground and ionic states of 1,2,5-thiadiazoles: An UVphotoelectron spectroscopic and theoretical study. J. Mol. Struct. 2010, 966, 85-91. 59. Maeyama, T.; Yagi, I.; Fujii, A.; Mikami, N. Photoelectron spectroscopy of microsolvated benzophenone radical anions to reveal the origin of solvatochromic shifts in alcoholic media. Chem. Phys. Lett. 2008, 457, 18-22. 60. Gan, S.; Luo, W.; He, B.; Chen, L.; Nie, H.; Hu, R.; Qin, A.; Zhao, Z.; Tang, B. Z. Integration of aggregation-induced emission and delayed fluorescence into electronic donor–acceptor conjugates. J. Mater. Chem. C 2016, 4, 3705-3708. 23 ACS Paragon Plus Environment


61. Zhang, Q.; Li, B.; Huang, S.; Nomura, H.; Tanaka, H.; Adachi, C. Efficient blue organic lightemitting diodes employing thermally activated delayed fluorescence. Nat. Photonics 2014, 8, 326. 62. Huang, S.; Zhang, Q.; Shiota, Y.; Nakagawa, T.; Kuwabara, K.; Yoshizawa, K.; Adachi, C. Computational Prediction for Singlet- and Triplet-Transition Energies of Charge-Transfer Compounds. J. Chem. Theory Comput. 2013, 9, 3872-3877. 63. Isegawa, M.; Truhlar, D. G. Valence excitation energies of alkenes, carbonyl compounds, and azabenzenes by time-dependent density functional theory: Linear response of the ground state compared to collinear and noncollinear spin-flip TDDFT with the Tamm-Dancoff approximation. J. Chem. Phys. 2013, 138, 134111. 64. Müller, T.; Lischka, H. Simultaneous calculation of Rydberg and valence excited states of formaldehyde. Theor. Chem. Acc. 2001, 106, 369-378. 65. Rubio, M.; Merchán, M.; Ortí, E.; Roos, B. O. A theoretical study of the electronic spectrum of naphthalene. Chem. Phys. 1994, 179, 395-409. 66. Becke, A. D.; Savin, A.; Stoll, H. Extension of the local-spin-density exchange-correlation approximation to multiplet states. Theor. Chim. Acta 1995, 91, 147-156. 67. Wang, F.; Ziegler, T. The performance of time-dependent density functional theory based on a noncollinear exchange-correlation potential in the calculations of excitation energies. J. Chem. Phys. 2005, 122, 074109. 68. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Insights into Current Limitations of Density Functional Theory. Science 2008, 321, 792-794. 69. Grofe, A.; Chen, X.; Liu, W.; Gao, J. Spin-Multiplet Components and Energy Splittings by Multistate Density Functional Theory. J. Phys. Chem. Lett. 2017, 8, 4838-4845. 70. Moore, C. E., Atomic Energy Levels. Circular of the National Bureau of Standards 467: Washington DC, U.S., 1949; Vol. 1. 71. Bernard, Y. A.; Shao, Y.; Krylov, A. I. General formulation of spin-flip time-dependent density functional theory using non-collinear kernels: Theory, implementation, and benchmarks. J. Chem. Phys. 2012, 136, 204103.


Page 24 of 35

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure captions. Figure 1. The upper panel shows the two components of the MRSF-TDDFT reference state denoted by black (MS = +1) and red (MS = −1) arrows. The zeroth-order MR RDM

0MR

which is a combination of RDMs of MS

= +1 and −1 triplet components is used in MRSF-TDDFT, whereas only MS = +1 RDM is used in SF-TDDFT. The lower panel shows the electronic configurations generated by spin-flip linear response from MR-RDM, which are given by blue, black, and red arrows. Closed, open, and virtual orbital sub-spaces are denoted by O, V, and C, respectively. The blue arrows represent the electronic configurations obtained from both components of the reference. They require a symmetrization procedure to eliminate the OO-type spin contamination. The black and red arrows represent the electronic configurations generated from Ms = +1 and −1 components, respectively. The blue and black arrows show configurations generated by the standard SF-TDDFT. Configurations which cannot be recovered with the MRSF-TDDFT method are shown by gray dashed arrows.

Figure 2. Correlation plots for the vertical excitation energies (in eV) obtained at the MRSF-TDDFT level of theory with the BHHLYP DFT functional and the 6-31G* basis set compared (a) against the TBE-2 reference set of excitation energies and (b) against the excitation energies obtained with cc-pVTZ basis set. The plots are shown at two extreme values of the spin-pairing coefficient used for these calculations, i.e., cSP = 0.3 (black circles) and cSP = 0.9 (red circles). Values of the correlation coefficient r at these extreme values of cSP are also given.

Figure 3. Mean absolute errors and variances for the vertical excitation energies between the TBE-2 reference and calculated values obtained at SF-, MRSF-, and LR-TDDFT level of theory with various DFT functionals and 6-31G* basis set. In the case of the MRSF calculations, mean absolute errors are shown at a value of spinpairing coefficient cSP = 0.3. Variances (in eV2) are shown as error bars.

Figure 4. Configurations generated by (a) the MRSF-TDDFT method and (b) the LR-TDDFT method for atoms and atomic ions. For brevity, only the MRSF configurations originating from the MS= +1 component of the mixed reference are shown; configurations originating from the MS= −1 component can be obtained by exchanging alpha and beta spins. In the case of the LR method, the closed-shell MS= 0 reference and configurations obtained from it are sown.



Figure 5. Dependence of the excitation energies of the five lowest singlet and three lowest triplet states on the spin-pairing coefficient cSP for the C atom. The filled circles and squares indicate the energies of the singlet and triplet states, respectively. Symmetries and multiplicities of the states are shown in the legend. The dependence was calculated with various hybrid DFT functionals based on the Becke exchange and the Lee-Yang-Parr correlation, i.e., B(1−cHF)HF(cHF)-LYP: (a) BHHLYP = B(0.5)HF(0.5)-LYP, (b) B(0.2)HF(0.8)-LYP, (c) B(0.1)HF(0.9)-LYP, and (d) B(0.0)HF(1.0)-LYP with 6-31G* basis set. Zero-energy level corresponds to the MS= +1 ROHF reference of MRSF-TDDFT.

Figure 6. Splitting between different atomic multiplet components (a) 3P − 1S (left panel), 1P − 1S (right panel), (b) 1D − 3P, and (c) 1S − 3P obtained by using MRSF-TDDFT with the STG1X functional developed in this work as well as the BHHLYP functional with the 6-31G* basis set. The experimental values for the splitting are also shown. The values of total MAE (in eV) are given in the legend.

Figure 7. Comparison of the absolute errors of the ST gaps for the small- and medium-sized molecules calculated with various density functionals by using (a) the MRSF-TDDFT and (b) LR-TDDFT method with the STG1X functional and 6-31G* basis set. For the majority of the molecules, adiabatic ST gaps were calculated, except c-C4H4, C3H6, i-C7H14, n-C7H14, C2H4, butadiene, and hexatriene for which vertical ST gaps were calculated. ‘c’, ‘i’, and ‘n’ stand for cyclo-, iso-, and normal-, respectively. The values of total MAE (in eV) are given in the legend.

Figure 8. Comparison of the absolute errors of the adiabatic ST gaps for the large-sized molecules calculated by using MRSF-TDDFT with the STG1X functional and 6-31G* basis set and compared against the ST gaps calculated with the TT-PBEh functional by using LR-TDDFT formalism from the work of van Voorhis.50 Only a few cases with a large absolute error in Ref. [50] are considered.


Page 26 of 35

Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figures. Figure 1.



Figure 2.

Figure 3.


Page 28 of 35

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4.



Figure 5.


Page 30 of 35

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6.



Figure 7.


Page 32 of 35

Page 33 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8.



Page 34 of 35

Table captions. Table 1. Mean absolute errors MAE (in eV) for the vertical excitation energies computed at MRSF-, SF, and LR-TDDFT level of theory with various DFT functionals and 6-31G* basis set. In the case of MRSF-TDDFT, the dependence of MAE on the spin-pairing coefficient cSP is also reported.

cSP

MRSF-TDDFT B3LYP

PBE0

BHHLYP

M08-SO

0.3

0.52

0.46

0.57

0.57

0.4

0.51

0.46

0.59

0.59

0.5

0.51

0.46

0.61

0.61

0.6

0.52

0.48

0.63

0.64

0.7

0.54

0.50

0.66

0.67

0.8

0.56

0.51

0.70

0.71

0.9

0.58

0.54

0.74

0.75

0.57

0.54

0.71

0.48

SF-TDDFT (0.0)

0.58

0.52 LR-TDDFT

0.30

0.34


Page 35 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic.


Performance Analysis and Optimization of Mixed-Reference Spin-Flip

Recommend Documents