Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
pubs.acs.org/JCTC
Combined Density Functional and Algebraic-Diagrammatic Construction Approach for Accurate Excitation Energies and Transition Moments Dávid Mester* and Mihály Kállay* Department of Physical Chemistry and Materials Science, Budapest University of Technology and Economics, P.O. Box 91, H-1521 Budapest, Hungary Downloaded via RUTGERS UNIV on August 8, 2019 at 13:15:23 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
S Supporting Information *
ABSTRACT: A composite of time-dependent density functional theory (TDDFT) and the second-order algebraicdiagrammatic construction [ADC(2)] approach is presented for efficient calculation of spectral properties of molecules. Our method can be regarded as a new excited-state doublehybrid (DH) approach or a dressed TDDFT scheme, but it can also be interpreted as an empirically tuned ADC(2) model. Several combinations of exchange−correlation functionals and spin-scaling schemes are explored. Our bestperforming method includes the Perdew, Burke, and Ernzerhof exchange and Perdew’s 1986 correlation functional and employs the scaled-opposite-spin approximation for the higher-order terms. The computation time of the new method scales as the fourth power of the system size, and an efficient cost-reduction approach is also presented, which further speeds up the calculations. Our benchmark calculations show that the proposed model outperforms not only the existing DH approaches and ADC(2) variants but also the considerably more expensive coupled-cluster methods.
1. INTRODUCTION The spectral properties of molecular systems can be accurately calculated using coupled-cluster (CC) methods,1 such as the CC singles and doubles (CCSD)2 or the approximate CC singles, doubles, and triples (CC3)3 models. However, the computational expenses of these approaches scale as the sixth or higher power of the system size, and to bigger molecules, only more approximate methods are applicable. One of them is the second-order approximate CCSD (CC2) method proposed by Christiansen and co-workers,4 which has a scaling of N5, where N is a measure of the system size. An even more approximate method is the configuration interaction singles (CIS) with perturbative doubles [CIS(D)] approach of HeadGordon et al.,5 where the CIS excitation energy is improved by a second-order correction for double excitations. One of the most successful approaches among the N5-scaling theories is the second-order algebraic-diagrammatic construction [ADC(2)] method.6 ADC(2) was developed by Schirmer exploiting a diagrammatic perturbation expansion of the polarization propagator and Møller−Plesset (MP) partitioning of the Hamiltonian. The major advantage of ADC(2) in comparison to CC2 is its efficiency. The ground-state energy is evaluated just at the second-order MP (MP2) level, and the matrix to be diagonalized for the excited states is Hermitian. Consequently, only one system of equations must be solved for each excited state to compute excitation energies and transition moments, whereas for CC2, the solution of three systems per © XXXX American Chemical Society
excited state is needed in addition to two further systems of equations for the ground state. Over the past decade, the scope of ADC(2) has been significantly extended by Dreuw and coworkers,7−17 Köhn et al.,18,19 and Hättig and co-workers.20,21 The performances of ADC(2), CC2, and other approximate methods were analyzed in extensive benchmark studies,10,22−26 and it has been found that the accuracy of ADC(2) is very close to that of CC2. It means that the computed transition energies for valence excitations are of similar quality as or even better than those obtained with the more expensive CCSD method, while for Rydberg excitations, less satisfactory results can be expected. For the excited states of even bigger molecules, the somewhat less accurate time-dependent density functional theory (TDDFT) with the usual density functional approximations is the method of choice. It can derived from density functional theory (DFT) relying on the linear-response formalism.27−31 TDDFT using pure exchange−correlation (XC) functionals scales as N3, but for semiquantitative accuracy, hybrid functionals are recommended, which scale as N4. While hybrid TDDFT excitation energies and spectral intensities are quite good for valence excitations, those for Rydberg and charge transfer (CT) states or excitations of extended π-electron systems can still be qualitatively Received: April 21, 2019 Published: July 2, 2019 A
DOI: 10.1021/acs.jctc.9b00391 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation incorrect.32−35 These problems can be mitigated by the more advanced range-separated hybrid functionals.31 Besides the aforementioned problems, in the wellestablished adiabatic approximation of TDDFT, only states dominated by one-electron excitations can be modeled. However, the explicit inclusion of double and higher excitations is a must for the adequate description of particular excited states. To remedy this problem several attempts have been made to include excitations beyond singles.31 Maitra and co-workers proposed the dressed TDDFT approach where one double excitation is included at a time.36 A more general propagator equation, which allows for any number of doubly excited states, has been developed by Casida utilizing the equation-of-motion superoperator approach.37 An even more rigorous polarization propagator approach to dressed TDDFT has been presented by Casida and Huix-Rotllant,38 who also derived explicit formulas for an ADC(2)-like method based upon a Kohn−Sham (KS) zeroth-order Hamiltonian. Improved frequency-dependent XC kernels have also been derived by Gritsenko and Baerends utilizing the common energy denominator approximation,39 as well as by Romaniello et al.,40 Sangalli and co-workers, 41 and Rebolini and Toulouse42 relying on the Bethe−Salpeter equation. The performances of the dressed TDDFT methods have been assessed in several publications, and significant improvement in excitation energies and state characters have been observed for states of double-excitation character.43−47 The performance of density functional approximations can also be improved by combining them with wave function methods. Nowadays, the most popular mixed approaches are the double-hybrid (DH) functionals, which were pioneered by Grimme48 following the related multicoefficient correlation methods of Zhao and Truhlar.49−51 In Grimme’s B2PLYP (Becke’s exchange with perturbative correction and Lee− Yang−Parr correlation) approach, a hybrid KS calculation is performed, and the energy is augmented with an MP2-like second-order perturbation theory (PT2) correction evaluated on the KS orbitals. Several variants of the B2PLYP functional were later proposed with different parametrizations.52−54 Nonempirical DH functionals, which can be derived from the adiabatic connection formalism, were developed by Toulouse et al.,55,56 Adamo and co-workers,57,58 and Chai et al.59,60 Spin-scaled DH variations were also proposed, where the PT2 correction is replaced by the spin-component-scaled (SCS)61 or scaled-opposite-spin (SOS)62 MP2 correction. The first method in the former category was developed by Chai and Head-Gordon,63 while the combination of SOS-MP2 and DFT was proposed by Goerigk and Grimme.64 Subsequently, various spin-scaled DHs were tested in several groups.65−69 Another prototype of DHs was put forward by Zhang and coworkers,70 where the KS equations are solved with a functional different from that used at the evaluation of the final energy. The accuracy and efficiency of DH functionals have been demonstrated in numerous studies, and their superiority to conventional DFT methods has been proven.64,69,71,72 Besides the MP2-based DH functionals, numerous other attempts have been made to combine wave function and DFT approaches. We should mention the DHs utilizing the random phase approximation,73−75 the CC and DFT blends,76−79 or the multiconfigurational DFT methods.80−91 For an exhaustive survey on combined wave function and DFT methods, see the recent review by Truhlar and co-workers.92
The application of DHs to excited states is a relatively unexplored field. The first step in this direction was made by Grimme and Neese, who extended the B2PLYP method to excited states.93 In their approach, which has been followed in all of the excited-state DH schemes, a TDDFT calculation is performed omitting the PT2 correction of the functional, and subsequently, the effect of double excitations is added by calculating the (D) correction of the CIS(D) method using the TDDFT excitation amplitudes. The transition moments are evaluated using the conventional TDDFT expressions, and no second-order correction is added. The excited-state DH theory was later extended to the B2GPPLYP,94 PBE0-DH, and PBE0295 DHs, and recently, Schwabe and Goerigk successfully combined time-dependent DH DFT theory with spin scaling.96 The performance of DHs for excited states has been benchmarked in several studies.97−100 The results show that DHs have significantly better error statistics for excitation energies than do common functionals. DHs are also free from the major problem of conventional functionals and the spurious low-lying states induced by the self-interaction error and also give accurate excitation energies for extended πelectron systems.101,102 In this paper, we go one step further and explore various possibilities for the combination of TDDFT with ADC(2) theory. We demonstrate that the composite method improves the efficiency of excitation energy and transition probability calculations with respect to both DH approaches and ADC(2).
2. THEORY 2.1. Double-Hybrid Density Functionals. In global hybrid density functional theory,103 the XC energy of a ground-state system is expressed as H E XC = (1 − αX)E XDFT + αXE XHF + ECDFT
EXDFT
(1)
ECDFT
where and are the semilocal exchange and correlation energies, respectively, EXHF denotes the exact exchange energy, and αX stands for the mixing factor of the semilocal and exact exchange contributions. In the simplest global DH functionals, the semilocal correlation energy is combined with an MP2-like PT2 correction,48 EPT2 C , as DH E XC = (1 − αX)E XDFT + αXE XHF + (1 − αC)ECDFT
+ αCECPT2
(2)
with α C as a scaling parameter for the correlation contributions. In the more complicated spin-scaled DHs,67 the opposite-spin (OS) and same-spin (SS) contributions to the PT2 correlation energy, EOS‑PT2 and ESS‑PT2 , are scaled C C OS SS separately by factors αC and αC , respectively, as DH E XC = (1 − αX)E XDFT + αXE XHF + (1 − αC)ECDFT
+ αCOSECOS‐PT2 + αCSSECSS‐PT2
(3)
EDH XC is frequently augmented with dispersion corrections, which are not considered here. We also note that only DHs based on PT2 correlation are studied here, and we do not consider other DHs utilizing more advanced correlation models. Excitation energies for hybrid DFT approaches can be obtained via time-dependent DFT27−29,104 by solving the ij A B yzij X yz iXy jj zzjj zz = ωjjj zzz k−B −A {k Y { kY{
B
(4) DOI: 10.1021/acs.jctc.9b00391 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation generalized eigenvalue problem with ω as the excitation energy and X and Y as the single excitation and de-excitation vectors, respectively. Using the spin−orbital representation, the elements of matrices A and B are defined by
generates double excitations with amplitudes tab ij . Here, we have introduced the a+ and i− creation and annihilation operators, respectively, for the corresponding spin orbitals. The excited-state wave function is obtained as |Ψ ADC(2)⟩ = (R1̂ + R̂ 2)|ΨMP1⟩
Aai , bj = (εa − εi)δijδab + (ai|bj) − αX(ij|ab) + (1 − αX)(ai|fX |bj) + (ai|fC |bj)
where excitation operators R̂ 1 and R̂ 2 are analogous to T̂ 2 with rai and rab ij as the corresponding amplitudes. The latter and the excitation energy can be computed by diagonalizing the Hermitian ADC(2) Jacobian, AADC(2), with elements21
(5)
and Bai , bj = (ai|bj) − αX(bi|aj) + (1 − αX)(ai|fX |bj) + (ai|fC |bj)
AaiADC(2) = AaiCIS , bj + , bj
(6)
ab ADC(2) Aaibj , ckdl = δacδbdδikδjlDij
(13)
|Ψab ij ⟩
[ 2] Aaĩ ADC(2) = AaiCIS , bj + Aai , bj − , bj
(7)
∑ c
