Faculty of Chemical Technology and Biotechnology George Olah Doctoral School
Efficient Excited-State Quantum Chemical Approaches for Large Molecules
D´ avid Mester
Advisor: Prof. Mih´ aly K´ allay
Department of Physical Chemistry and Materials Science MTA-BME Lend¨ ulet Quantum Chemistry Research Group
1 Introduction 1
2 Overview 3
3 Theoretical background 10
3.1 CIS and TDHF theories . . . 10
3.2 TDDFT and TDA-TDDFT theories . . . 14
3.3 CC2 and ADC(2) theories . . . 16
3.4 Related reduced-cost approximations . . . 21
3.4.1 Density fitting and natural auxiliary function approximations . . 21
3.4.2 Natural orbital approximation . . . 22
3.5 Localization of molecular orbitals and its benefits . . . 22
4 Reduced-cost approximations for excited states 25 4.1 CC2 singlet excitation energies. . . 25
4.1.1 Theory and implementation . . . 25
4.1.2 Results . . . 32
4.2 ADC(2) excitation energies and transition moments . . . 37
4.2.1 Theory and implementation . . . 37
4.2.2 Results . . . 43
5 Reduced-scaling approximations for excited states 51 5.1 Local density fitting approach for TDDFT and CIS methods . . . 51
5.1.1 Theory and implementation . . . 51
5.1.2 Results . . . 56
5.2 Local domain-based approach for correlation methods . . . 64
5.2.1 Theory and implementation . . . 65
5.2.2 Results . . . 71
6 Improved double hybrid TDDFT approach 85 6.1 Theory and implementation . . . 85
6.2 Results . . . 87
7 Summary 100
A Appendix: Molecules and their sizes 105
B Appendix: Supporting calculations 109
ADC(2) second-order algebraic-diagrammatic construction
AO atomic orbital
BP atom list Boughton–Pulay atom list
CC2 approximate second-order coupled cluster singles and doubles
CIS configuration interaction singles
CIS(D) CIS with perturbative second-oder correction double excitations CT excitation charge transfer excitation
DIIS algorithm direct inversion in the iterative subspace algorithm
DF density fitting
DH double hybrid
EOM-CC theory equation-of-motion coupled-cluster theory
ERI electron repulsion integral
I/O operation input/output operation
LDF local density fitting
LMO local molecular orbital
LR-CC theory linear-response coupled-cluster theory
LT Laplace transform
MAE mean absolute error
MAX maximum absolute error
ME mean error
MO molecular orbital
MP2 second-order Møller–Plesset
NAF natural auxiliary function
CS-NAF complete MO space NAF
RS-NAF restricted MO space NAF
NO, VNO natural orbital and virtual natural orbital
PAO projected atomic orbital
RMS deviation root mean square deviation
SCF self-consistent field
SCS correction spin-component-scaled correction SOS correction scaled-opposite-spin correction
SS or OS contribution same spin or opposite spin contribution SG training set training set of Schwabe and Goerigk
SVD singular value decomposition
TDA Tamm–Dancoff approximation
TDDFT time-dependent density functional theory TDHF theory time-dependent Hartree–Fock theory
µ, ν, . . . atomic orbitals
i, j, . . . (quasi-)canonical occupied orbitals i0, j0, . . . localized occupied orbitals
a, b, . . . (quasi-)canonical virtual orbitals p, q, . . . general (quasi-)canonical orbitals P, Q, . . . auxiliary functions
nAO, naux number of atomic orbitals and auxiliary functions nocc, nvirt number of occupied and virtual orbitals
C general transformation matrix
cai single excitation CIS coefficients cabij double excitation CIS(D) coefficients
Tia, Tijab CC2 ground state single and double amplitudes tabij MP2 ground state amplitudes
ria, rabij single and double excitation coefficients T1, T2 single and double cluster operator R1, R2 single and double excitation operator
H,Hˆ normal and similarity-transformed Hamiltonian F,Fˆ normal and similarity-transformed Fockian A,Aeff,A˜ normal, effective, and modified effective Jacobian
σia sigma vector
ω excitation energy
f oscillator strength
dˆ dipole operator
J three-center integrals in the (quasi-)canonical basis
ˆJ,J,˜ J three-center integrals in the CS-NAF, NO, and RS-NAF basis εVNO, εNAF thresholds for the virtual NO and NAF truncation
TLDF, TLMO, TPAO thresholds for the LDF, LMO, and PAO truncation αX, αC scaling factors for the exchange- and correlation terms αOSC , αSSC scaling factors for opposite- and same-spin terms
P local domain
Chapter 1 Introduction
Our efforts on the development of quantum chemical methods focus on the efficient and accurate solution of the Schr¨odinger equation. The first contributions to this field dates back to the 1930s, while the improvements were most accelerated by advancements in algorithm design and the performance of the available computers, which began in 1980 and it still lasts today. Since then, the results achieved in the field of quantum chemistry have been awarded with two Nobel Prizes in Chemistry. Over the past decades, one significant avenue of the developments has been the elaboration of more accurate meth- ods and the implementation of the corresponding complicated equations in the program packages, for example, via automatic code generation techniques. Many high-accuracy quantum chemical methods were successfully developed, which allow the theoretical pre- diction or understanding of experimentally measurable physical and chemical properties, such as equilibrium structures, thermodynamic data, absorption spectra, or dipole mo- ments. The accuracy of the obtained results is comparable to or more and more often even better than that of the experiments.
However, the solution of the corresponding systems of equations containing billions of unknown variables is extremely demanding. For more elaborate methods, the number of the required operations scales as the 7th or 8th power of the system size. It is easy to see that in the latter case, by doubling the size of the studied molecule, the time required for the calculations increases by a factor of 256. Therefore, the upper limit of the applicability of the most accurate methods is only a few atoms. Accordingly, the substantial aim of the method development in recent years has been to reduce the calculation times of our existing, highly reliable methods, or to elaborate new schemes by modifying or combining the theories.
External perturbations, such as electromagnetic radiation, of certain strength can lead to the rearrangement of the electronic structure and to electronic excitation. There
are several actively researched phenomena related to the excited electronic states of molec- ular systems. For instance, excited states play an important role for photochromic ma- terials, photo-initialized chemical processes, and energy transfer and storage. Quantum- chemical methods have now become routine tools in the investigation of excited-state properties and processes in the fields of spectroscopy, analytical- and biochemistry, where often extensive molecular systems are studied. Fortunately, the error of the calculated ground- and excited-state energy cancel each other in general, thus it is permissible to use lower scaling methods for fairly accurate excited-state properties. However, as it can be seen in Fig. 1.1, the upper limit of the applicability of such methods is still easily accessible, despite their relatively low scaling. Our goal is to expand the current limit of
Figure 1.1: The approximate time required for a fourth- (blue) and fifth-power-scaling (orange) calculation as a function of the system size.
accurate excited-state methods and to be able to model chemically relevant systems that can not be treated currently due to the extensive computational requirements.
The most important aim of this thesis is to introduce approximations which can effectively decrease the number of variables. In some cases, these are just numerical mathematical transformations, however, in the other group of the approximations we carry out transformations controlled by specific physical or chemical considerations. In this case, drastic reduction in the computational demand can be gained. The final goal is to find and combine the most appropriate transformations and minimize the number of variables, while the error of the approximations introduced by us is on average still an order of magnitude smaller than the intrinsic error of the original method. These results are presented in Chapters 4and 5.
In addition, we have also developed an effective and relatively cheap method to calculate accurate excited-state properties. The proposed fourth-power-scaling method outperforms not only the existing similar methods but also the considerably more expen- sive approaches. These results are detailed in Chapter 6.
Chapter 2 Overview
Nowadays, time-dependent density functional theory (TDDFT), which is derived from density functional theory (DFT) through the linear-response formalism, is the most common choice to investigate time-dependent properties of molecular systems, such as excitation energies, polarizabilities, and chiroptical properties [1–5]. The developed an- alytical TDDFT gradients also enable the efficient calculation of excited-state equilib- rium structures and other molecular properties [6–9]. The time-dependent Hartree–Fock (TDHF) method  is an analogue of TDDFT where one chooses the Hartree–Fock (HF) solution as the ground-state reference. Another related HF-based excited-state method is the configuration interaction singles (CIS)  approach. It is well-known that TDDFT is generally more accurate despite the similar computational requirements, at least for valence excitations [12–16], due to the consideration of dynamic electron correlation. In most practical applications for bigger molecules TDDFT is used, and the simpler CIS method is rarely employed directly due to its relatively large inaccuracy. Nevertheless, CIS has growing importance in the emerging reduced-cost and reduced-scaling higher- order correlated excited-state techniques, which utilize the information obtained from the CIS wave function. Such approaches have been developed for the second-order coupled- cluster (CC) [17–20], the second-order algebraic-diagrammatic construction (ADC) , and higher-order CC [22,23] methods. In the corresponding calculations for extended systems, the time spent on the initial canonical CIS calculation, due to the effective ap- proximations for the correlated model, is significant compared to the overall computation time.
The computational expenses of the CIS, TDHF, and TDDFT methods, as long as the DFT functional contains HF exchange contribution, scale with the fourth power of the system size. Accordingly, the upper limit for the applicability of these theories, if no further approximations are introduced, is roughly 150 atoms. Since the size of the chemically relevant systems in many cases exceeds this region, a number of attempts have been made over the last two decades to reduce the computation time and the resources
required for the calculations. The simplest approximate TDDFT method is the Tamm–
Dancoff approximation (TDA), which was proposed by Hirata and Head-Gordon .
The TDA-TDDFT approach gives similarly accurate results for hybrid functionals as the parent TDDFT method, while computation times are approximately halved. If TDA is applied to TDHF, the CIS method is retrieved. A commonly used and implemented scheme is the density fitting (DF) approach, which was introduced by Shavitt et al. 
and later enhanced by Whitten  and Dunlap and co-workers . In the DF approach the four-center electron repulsion integrals are expressed as the products of two- and three-center integrals, therefore the operation count, the memory requirement, and the number of the input-output (I/O) operations can be greatly decreased [28–32]. While the DF technique was originally developed for ground-state HF and Kohn–Sham (KS) calculations , it has also been successfully employed for (TDA-)TDDFT . Besides the speedup of the calculation of the corresponding matrix elements, the strategies for the solution of the TDDFT equations have also been improved. Noteworthy is the smart trick recently introduced by Furche  and Mart´ınez and co-workers , where the equations are solved in a non-orthonormal vector space, thus the prescreening of the integrals can be performed more efficiently.
Another widely used technique to reduce the computational expenses is to split up the system studied into smaller parts and treat the individual subsystems with var- ious methods and accuracy. Such approaches include, for example, the fragment-based TDDFT [37–42] and CIS  methods, the frozen-density embedding approximation developed by Neugebauer et al. [44–47], and the density-fragment interaction scheme proposed by Fujimoto and Yang . Many approaches use localized molecular or- bitals (MOs), such as the “divide-and-conquer” TDDFT introduced by Nakai and co- workers [49,50], and it is also worthwhile mentioning the absolutely localized MO- [51–53], in situ optimized localized MO- [54,55], and non-orthogonal localized MO-based ap- proaches . Further related methods are the simplified TDA-TDDFT [57,58] approach of Grimme, the state-specific TDDFT [59,60], and the restricted virtual space TDDFT theory . The approximations outlined above have many advantages, but some defi- ciencies are also present. For instance, some of them cannot be applied in a black box manner, which means that the user has to define the subsystems of the molecule. Par- ticular methods cannot be used with hybrid functionals, whereas other schemes can only be applied to periodic systems. Furthermore, quantitative results cannot be obtained by most of the approaches due to the large error of the approximations, or, compared to genuine TDDFT calculations, the size of the investigated systems cannot be significantly increased. To accelerate TDDFT calculations utilizing hybrid functionals with a tolera- ble error, the chain of spheres exchange (COSX) approach  and the auxiliary-density matrix method (ADMM)  were developed by Neese and Kjærgaard and co-workers.
Both schemes reduce the computation time required to calculate the rate-determining ex- act exchange contributions. In the case of COSX, a seminumerical integration technique related to the pseudospectral approach [64–66] is used, with which the average speedup factor is about 2−4 for the excitation energy calculations for systems including 50−150 atoms. The mean (maximum) error originating from the approximation is 0.02 (0.07) eV using the BHLYP  functional. In ADMM, the exchange contribution is calculated using a smaller auxiliary basis set, and the contribution is extrapolated to the larger ba- sis. This approach brings similar speedups as COSX, but the error of the approximation for TDDFT excitation energies may exceed 0.1 eV. For the calculation of the exchange contribution, an asymptotic linear scaling can be achieved with both methods.
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 incorrect [68–71]. These problems can be mitigated by the more advanced range-separated hybrid functionals . Besides the aforementioned problems, in the well-established adiabatic approximation of TDDFT only states dominated by one-electron excitations can be modeled. However, the ex- plicit 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 . Maitra and co-workers proposed the dressed TDDFT approach, where one double excitation is included at a time . A more gen- eral propagator equation, which allows for any number of doubly excited states, has been developed by Casida utilizing the equation-of-motion superoperator approach . An even more rigorous polarization propagator approach to dressed TDDFT has been pre- sented by Casida and Huix-Rotllant , who have also derived explicit formulas for a second-order ADC-like method based upon a Kohn–Sham (KS) zeroth-order Hamilto- nian. Improved frequency-dependent exchange-correlation (XC) kernels have also been derived by Gritsenko and Baerends utilizing the common energy denominator approxima- tion , as well as by Romaniello et al. , Sangalli and co-workers , and Rebolini and Toulouse  relying on the Bethe–Salpeter equation. The performance of the dressed TDDFT methods have been assessed in several publications, and significant improvement in excitation energies and state characters were observed for states of double excitation character [79–83].
The performance of density functional approximations for ground-state calculations 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 pio- neered by Grimme  following the related multi-coefficient correlation methods of Zhao and Truhlar [85–87]. In Grimme’s B2PLYP (Becke’s exchange with perturbative correc- tion and Lee–Yang–Parr correlation) approach a hybrid KS calculation is performed, and
the energy is augmented with a second-order perturbation theory (PT2) correction eval- uated on the KS orbitals. Several variants of the B2PLYP functional were later proposed with different parametrizations [88–90]. Nonempirical DH functionals, which can be de- rived from the adiabatic connection formalism, were developed by Toulouseet al. [91,92], Adamo and co-workers [93,94], and Chaiet al.[95,96]. Spin-scaled DH variations were also proposed, where the PT2 correction is replaced by the spin-component-scaled (SCS) 
or scaled-opposite-spin (SOS)  second-order Møller–Plesset (MP2) correction. The first method in the former category was developed by Chai and Head-Gordon , while the combination of SOS-MP2 and DFT was proposed by Goerigk and Grimme . Sub- sequently, various spin-scaled DHs were tested in several groups [101–105]. Another pro- totype of DHs was put forward by Zhang and co-workers , 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 [100,105,107,108].
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 [109–111], the CC and DFT blends [112–115], or the multi- configurational DFT methods [116–127]. For an exhaustive survey on combined wave function and DFT methods see the recent review by Truhlar and co-workers .
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 . In their approach, which has been followed in all 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 , PBE0-DH, and PBE0-2  DHs, and recently Schwabe and Goerigk successfully combined time-dependent DH DFT theory with spin scaling .
The performance of DHs for excited states has been benchmarked in several studies [133–
136]. 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, the spurious low-lying states induced by the self interaction error, and also give accurate excitation energies for extendedπ-electron systems [137,138].
One of the most accurate electronic structure theories is the coupled-cluster ap- proach introduced by ˇC´ıˇzek a half a century ago . The hierarchical CC methods based on the exponential parametrization of the wave function, that is, the CC singles and doubles (CCSD), the CC singles, doubles, and triples (CCSDT), . . . , enable the
consideration of electron correlation with arbitrary accuracy. The CC methods can also be generalized to excited states, but it is not a trivial task. One option is to invoke linear-response (LR) theory, which was first extended to CC theory by Monkhorst and co-workers [140,141], and later by Koch et al.[142,143] An alternative approach, which is equivalent to LR-CC theory for excitation energies, is the equation-of-motion (EOM) CC method developed by Bartlett and co-workers [144,145]. The theory, implementation, and performance of the various LR-CC and EOM-CC methods is reviewed in several publications in the literature [146–149]. Concerning bigger molecules even the relatively cheap CCSD approach is too expensive, and further approximations are required. A simple approximate CCSD approach is the second-order CC (CC2) method proposed and first implemented by Christiansen et al. [150–152] and later perfected by H¨attig and co-workers [31,153,154]. The CC2 method supplies excitation energies and transition mo- ments with a moderate error, at least for valence states, with respect to higher-order CC methods [155,156]. It is also worth mentioning the CIS(D) approach introduced by Head- Gordon et al. , which can also be regarded as an approximate CC2 method .
The CIS(D) approach improves the CIS excitation energy with a perturbative correction for the double excitations and scales as the fifth power of the systems size as the CC2 method, but, in practice, it is considerably less expensive than CC2.
Among the theories suitable for excited-state property calculations ADC is one of the most promising approaches. It is a Hermitian and size-consistent method, and it is rel- atively easy to implement. The ADC scheme was first derived by Schirmer  employ- ing a diagrammatic perturbation expansion of the polarization propagator, utilizing the Møller–Plesset partitioning of the Hamiltonian. A similar result was later obtained with the so-called intermediate state representation (ISR) approach developed by Schirmer et al. [160–162] While the initial implementations of the theory were limited to its second- order variant [ADC(2)], later it was extended to the third-order [ADC(3)] [163,164]. A more efficient implementation [165–168] of the ADC(2) and ADC(3) methods and ex- tensive benchmark calculations  were reported by Dreuw et al. In these studies the performance of ADC(2) was also compared to that of the closely-related but more de- manding CC2 approach, and it has been proven that the ADC(2) method is practically as accurate as CC2 [166,169]. Furthermore, tools to compute two-photon absorption , static polarizability , core-valence excitations , and excited-state dynamics 
were also developed. The ADC method was also combined with the spin-flip , the scaled-opposite-spin , and the frozen density embedding  approaches.
Despite the excellent numerical results the applicability of CC and ADC methods is limited due to the steep scaling of their costs with the system size. An alternative solution instead of the TDDFT methods could be the reduction of computational ex- penses of CC2 and ADC(2) approaches, which provide results consistently better than
common TDDFT methods [134,176,177]. In the case of most reduced-cost approaches the bottleneck related to the processing of four-center electron repulsion integrals (ERIs) is removed exploiting DF or Cholesky-decomposition (CD) techniques, which can also be combined with Laplace transform (LT) approximations. The DF approximation was adapted for CC2 by H¨attig and co-workers , who also developed analytic gradients for CC2 employing DF [153,154]. In the tensor hypercontraction (THC) scheme of Sherrill, Mart´ınez, and their co-workers, which is a generalization of the DF technique, an even lower-order representation of the ERI tensor is used [178–180]. The method was also successfully applied to ground-  and excited-state  CC2 calculations. In the CD approach proposed by Koch et al.[183,184] the four-center ERI tensor is decomposed as a product of triangular matrices neglecting the columns/rows that give negligible contri- butions. The efficiency of the approximation was also demonstrated for CC2 [185–187].
The LT approximation developed by Alml¨of and H¨aser [188–190] eliminates the orbital energy denominators appearing in many-body methods like CC2. Together with the DF and further approximations, it can be efficiently used for the scaling reduction of the CC2 method [191,192]. Another simple technique for reducing the costs of correlated excited- state methods is the restricted virtual space approach, where the high-lying canonical virtual molecular orbitals are neglected. This approach was also tested at the ADC(2) and CC2 levels [193–197].
Significantly more efficient methods can be developed by determining the MOs that play an important role in the excitation. One of the most popular schemes is the natural orbital (NO) approximation [198–200], with which the MO space where the equations are solved can be effectively reduced. In the NO approach, a one-particle density ma- trix, which is formed using a lower-level wave function, is diagonalized, and the orbitals with significant importance are selected from the resulting NOs. The approach is widely utilized for ground-state calculations [201–204], and after a few early attempts, its im- portance for excited-state theories started to increase recently [19,20,205–207]. The developed approaches are not only suitable for relatively cheap methods, such as ADC(2) and CC2 but could also extend the applicability of higher-order ab initio methods to medium-sized molecules.
Further computational savings can be achieved if one takes advantage of the local- ity of the MOs [208,209]. In this case, not only the time required for the calculations is decreased, but at the same time the scaling of the methods is also reduced. The first excited-state local approaches were presented by Korona and Werner  and Craw- ford et al. , who generalized the ground-state local CC singles and doubles (CCSD) method of Werner and co-workers  to EOM-CCSD. In the local EOM-CCSD method developed by the former authors, the doubles amplitudes were restricted using the in- formation by inspecting the configuration interaction singles (CIS) wave function ,
which idea has been taken over in several subsequent studies. Thereafter, Korona, Sch¨utz, Kats, and their co-workers developed various excited-state CC methods using local approaches [17,192,212–215]. In later publications the development of local CC2 and ADC(2) methods was reported [17,21], which were also extended to the calculation of molecular properties  and improved with Laplace transform techniques [192,213–215].
Parallel to those efforts, further papers were published by Russ and Crawford about the calculation of excited-state properties [216,217]. Promising results were also obtained by H¨attig et al. extending the pair natural orbital (PNO) approach to excited-state the- ories [218–221]. The chain of spheres exchange  and the back transformed PNO based [223–225] approaches developed by Izs´ak et al. The recent local framework for calculating excitation energies (LoFEx) [18,23,226] and the correlated natural transi- tion orbital framework (CorNFLEx)  approaches of Baudin and Kristensen introduce somewhat different strategies. The latter is an encouraging combination of the NO and the local approaches, where the reduced domains of the MOs are constructed by ana- lyzing an approximate second-order density matrix and considering distance criteria for the orbitals. A comprehensive study was recently published on the topic of reduced-cost approximations by Crawford, Kumar, and co-workers .
In this chapter the theoretical background of methods and approximations used in this thesis are outlined. We present briefly the most important equations and consid- erations, firstly, for the simplest methods, such as CIS, TDHF, and TDDFT, thereafter for the more elaborate ones, for example, the CC2 and ADC(2) approaches. Finally, the theoretical background of the approximations utilized in the method development is assessed in the last two subsections.
3.1 CIS and TDHF theories
The exact non-relativistic energy of a stationary system can be obtained by the solution of the time-independent Schr¨odinger equation:
H(r)Ψ(r) = EΨ(r), (3.1)
whereH is the Hamiltonian, Ψ stands for the exact wave function,Edenotes the total en- ergy, andris a vector for electron coordinates. The analytical solution of the Schr¨odinger equation is not possible for systems containing more than two particles. In such cases, the wave function describing the ground state system can only be constructed with certain approximations, for example, the Born–Oppenheimer approximation, where the move- ment of nuclei and electrons is decoupled. In addition, the exact wave function of the system is approximated by determinants which are constructed as the products of single- particle functions. For the starting point, the solution of the Hartree–Fock equations is used most widely, which can be written in the
F(r)Φ0(r) = E0Φ0(r) (3.2)
form, where F denotes the Fockian, Φ0 is the single Hartree–Fock determinant, which corresponds to the best single Slater determinant describing the electronic ground state of the system, and E0 is the sum of occupied orbital energies.
The CIS wave function, ΨCIS, is constructed as a linear combination of the HF and singly excited determinants:
ΨCIS(r) = Φ0(r) +X
where cai is the CIS coefficient and Φai is the excited determinant which is formed as Φai =a+i−Φ0, wherea+ andi− are the creation and annihilation operators corresponding to MOs a and i, respectively. In other words, the ith occupied orbital is replaced by the ath virtual orbital in the determinant. The notation follows the convention that i, j, . . ., a, b, . . ., and p, q, . . . denote the occupied, virtual, and general MO indices, respectively.
Due to Brillouin’s theorem, the excited state wave functions are orthogonal to the ground state. Consequently, inserting Eq. (3.3) into Eq. (3.1) and projecting both sides onto the space of singly excited determinants, the resulting equation can be written in the
hΦbj|H|Φaii=ECIScbj (3.4) form, with elements of
hΦbj|H|Φaii= (EHF+fab−fij)δijδab+ 2(ia|jb)−(ij|ab) (3.5) in spatial orbital basis, where ECIS (EHF) is the CIS (HF) total energy, fpq stands for the corresponding Fock-matrix element, and the two-electron integral in the Mulliken notation, (ia|jb), is defined as
(ia|jb) = Z Z
r12j(r2)b(r2)dr1dr2 , (3.6) where r1 is the coordinate of electron 1, and r12 = |r1 −r2|. Since the CIS approach is a variational method, all excited-state total energies are upper bounds to their exact values. The excitation energy is simply obtained as ω=ECIS−EHF. The final equation using matrix notation can be recast as
with the matrix elements
Aia,jb= (fab−fij)δijδab+ 2(ia|jb)−(ij|ab). (3.8) The excitation energies can be obtained as the solution of the secular equation by diag- onalization of Hermitian matrix A, whose dimension is noccnvirt, wherenocc (nvirt) is the number of occupied (virtual) orbitals. In general, excitation energies obtained with the CIS method are usually overestimated even with 2 eV. This is can be explained by the fact that virtual orbital energies are calculated for the (N+ 1)-electron system, whereN stands for the number of electors. Therefore, the orbital energy difference, which is the leading term in the above equation, is hard to be interpreted as an excitation.
Since only the few lowest solutions of the eigenvalue problem are important in practical applications, and the exact diagonalization cannot be performed for extended systems, the eigenvalue equation is solved iteratively for the lowest roots [228–230]. The rate-determining step of the iterative process is the construction of the so-called sigma vector, whose elements can be written as
fijcaj + 2X
(ij|ab)cbj . (3.9) At the end of the iterative procedure the converged CIS wave function is utilized for the evaluation of transition moments. The transition density matrix needed for the ground to excited state transition moments can be simply expressed as
ρpq =hΨCIS|p+q−|Φ0i = 0 0 c 0
The oscillator strength (f) in length gauge can be defined by f = 2ω
3 |hΨCIS|dˆ|Φ0i|2 = 2ω 3
|hΨCIS|dˆα|Φ0i|2, (3.11) where the sum of squares on the right-hand side is the dipole strength, and ˆd= ( ˆdx,dˆy,dˆz) is the dipole operator and its components. The x component of the transition electric dipole moment is obtained as
dxpqρpq =√ 2X
dxaiρai, (3.12) where dxpq stands for the x component of the dipole moment integrals in the MO basis.
The general expression to describe a physical system evolving with time is the time- dependent Schr¨odinger equation:
H(r, t)Ψ(r, t) = i∂
∂tΨ(r, t), (3.13)
where the time-dependent Hamiltonian, H(r, t), can be split up as
H(r, t) =H(r) +V(r, t), (3.14) whereV(r, t) is an arbitrary time-dependent external potential. Analogously to the above discussion, the exact wave function can also be approximated with a single reference determinant in this case. The resulting TDHF equations read as
[F(r) +V(r, t)]Φ(r, t) = i∂
∂tΦ(r, t). (3.15)
At t = 0 the system is stationary and can be described with the determinant Φ0. If a weak external time-dependent perturbation is turned on, the system responds to the changes. This response function is treated up to first order in linear response theory.
After lengthy algebraic modifications and Fourier transformation  the final equation
obtained can be written in matrix notation in the form of
=ω c c0
where c0 is the singles de-excitation vector, and B is the so-called coupling matrix with elements
Bia,jb= 2(ia|jb)−(ib|ja). (3.17)
As it can be seen, TDHF contains not only singly excited determinants but also singly de-excited determinants. Accordingly, it can be considered as an extension of the CIS method. However, the de-excitations are, of course, nonphysical since one cannot de- excite from the HF reference state. The magnitude of thec0 amplitudes can be regarded as a measure of the ground-state correlation, which, as a consequence, should be small with respect to the c amplitudes.
Inspecting the dimension of Eq. (3.16) we can observe that it is twice larger than that of Eq. (3.7). However, with some algebraic modifications the same dimension can be reached. For this purpose, we take the combinations of two separate equations, Ac+Bc0 =ωc and −Bc−Ac0 =ωc0 :
(A−B)(c−c0) = ω(c+c0), (3.18) and
(A+B)(c+c0) =ω(c−c0). (3.19) If the second expression is recast forc+c0 and this result is inserted into the first equation, the final equation can be read as
(A+B)(A−B)(c−c0) = ω2(c−c0). (3.20) Similar to the CIS method, the above equation is solved iteratively, and the rate-determining step is also the sigma vector calculation. The element of the vector is defined by the
σia = X
(bk|cj)˜cck (3.21) expression, where the ˜c=c−c0 shorthand notation is introduced. The transition density matrix used for the calculation of transition properties is obtained as
ρpq = 0 c0 c 0
where the excitation and de-excitation vectors can be expressed from the solution of the iterative procedure as
c0 = (A−B)˜c
ω −c. (3.24)
It is important to note that, if matrix Bis set to zero, the CIS method is recovered. This approximation is known as the Tamm–Dancoff approximation and will be utilized in the following section.
3.2 TDDFT and TDA-TDDFT theories
The density functional theory relies on the first Hohenberg–Kohn theorem, which states that one-to-one correspondence is present between the exact energy, the wave function, and the electron density, hereafter denoted by ρ(r). The total energy can be split up to the contribution of various terms:
E[ρ] =Ts[ρ] +J[ρ] +EXC[ρ] +EeN[ρ], (3.25) where Ts[ρ] is the exact kinetic energy of a non-interacting system, J[ρ] stands for the classical Coulomb energy, EXC[ρ] is the exchange-correlation energy, and EeN[ρ] denotes the potential energy of the electron-nucleus interaction. All contributions in Eq. (3.25), except EXC[ρ], can be calculated exactly. The performance of DFT methods depends on the accuracy of the XC contribution usually including empirical parameters. In global hy- brid density functional theory , the XC energy of a ground-state system is expressed as
EXCH = (1−αX)EXDFT+αXEXHF+ECDFT, (3.26) where EXDFT and ECDFT 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. The equations of the time-dependent density functional theory can be obtained through the linear response formalism with the similar considerations as in the TDHF theory if one chooses the Kohn–Sham solution as the ground state reference. The hybrid TDDFT excitation energies are also obtained through the solution of Eq. (3.16), in which the elements of matrices A and B can be expressed as
AHia,jb= (fab−fij)δijδab+ 2(ia|jb)−αX(ij|ab) + (1−αX)(ia|fX|jb) + (ia|fC|jb) (3.27) and
Bia,jbH = 2(ia|jb)−αX(ib|ja) + (1−αX)(ia|fX|jb) + (ia|fC|jb), (3.28)
where the DFT exchange term is written in the following form:
(ia|fX|jb) = Z Z
δρ(r1)δρ(r2)j(r2)b(r2)dr1dr2 , (3.29) where EX denotes the exchange functional. Similar expression holds for the DFT corre- lation. The algorithmic considerations are very similar to the TDHF case, only the DFT contributions should be added to the sigma vector. These cubic-scaling terms do not significantly affect the time required for calculations inspecting sufficiently large systems.
Furthermore, the overhead can be reduced by prescreening techniques for extended sys- tems. However, the prefactor is quite large, thus these steps can be rate-determining for medium-sized molecules.
In general, TDDFT methods provide more accurate excitation energies than the CIS or TDHF approaches. This can, at least partly, be explained by the fact that the difference of the Kohn–Sham orbital energies, which are the leading term of the diagonal elements of the AHmatrix, are usually better approximations for excitation energies. This is because the virtual KS orbital energies are calculated for the N-electron system contrary to the HF theory, where the virtual orbital energies are obtained for the (N+1)-electron system.
However, it is now well-known that TDDFT has severe problems even with qualitative description of Rydberg states, transitions of extended π-systems dominated by double excitation, and CT states. The problems with Rydberg states and extended π-systems can be attributed to the wrong long-range behavior of the standard XC functionals. In the case of excited CT states, the excitation energies are significantly underestimated, and potential energy curves do not exhibit the correct asymptotic behavior with respect to a distance coordinate between the occupied and virtual orbitals involved in the CT excitation.
The expression for the transition density is also identical to Eq. (3.22). The Tamm–
Dancoff approximation, when matrix BH is set to zero, can be applied as well. In this case the eigenvalue problem will be Hermitian, and the TDA-TDDFT excitation energies are obtained. Obviously, the transition properties could be calculated with the density defined in Eq. (3.10).
In some cases it may be required to use the more sophisticated global double hybrid functionals for accurate properties. In the simplest global double hybrid functionals the semilocal correlation energy is combined with an MP2-like PT2 correction , ECPT2, as EXCDH = (1−αX)EXDFT+αXEXHF+ (1−αC)ECDFT+αCECPT2 (3.30) withαC as a scaling parameter for the correlation contributions. In the more complicated spin-scaled DHs  the opposite-spin (OS) and same-spin (SS) contributions to the PT2 correlation energy, ECOS−PT2 and ECSS−PT2, are scaled separately by factors αCOS and αSSC ,
EXCDH = (1−αX)EXDFT+αXEXHF+ (1−αC)ECDFT+αOSC ECOS−PT2+αSSC ECSS−PT2 . (3.31) EXCDH 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. In the excited-state extension of the genuine DHs  defined by Eq. (3.30) TDDFT is combined with the CIS(D) method . Eqs. (3.16) or (3.7) are solved with modified ADH and BDH matrices, where the last term is scaled by (1 − αC). The CIS(D) excitation energy correction is computed with excitation vector c, multiplied by αC, and added to the TDDFT excitation energy. In the existing spin-scaled TD DH variants  the same- and opposite-spin contributions are scaled by different parameters in the “direct” and
“indirect” terms of the CIS(D) correction following the SCS-CIS(D) approach of Rhee and Head-Gordon . In both cases no correlation correction is added to the transition moments, but they are calculated by the standard TDDFT formulas [Eqs. (3.22) or (3.10)]. The only improvement at the calculation of transition probabilities, such as oscillator or rotator strengths, is that the corrected excitation energy is used in the corresponding formulas.
3.3 CC2 and ADC(2) theories
In the coupled cluster theory, the wave function is obtained as an exponential op- erator acting on the HF ground-state determinant:
ΨCC = eTΦ0 , (3.32)
where T is the cluster operator, with T =T1+T2+T3+. . ., where T1 =P
aiTiaa+i− is the cluster operator of single excitations. For convenience, the T2 operator together with the higher excitation cluster operators will be denoted in a general form as
Tn = 1 (n!)2
Tµnτµn , (3.33)
whereTµn is the cluster amplitude associated with theτµn excitation operator. One of the most common ground-state CC methods is the CC singles and doubles approach, where the excitations are treated up to doubles. In this case the eT = eT1+T2 operator can be expanded in a Taylor-series in the
eT = 1 +T +1
2T2+· · ·= 1 +T1 +T2 +1
2T22+. . . (3.34) form. Neglecting higher order terms in the CC equations, the approximate second order CC method can be gained, which is a very efficient approach for excited-state calculations.
In this case, the CC2 correlation energy for a restricted Hartree–Fock reference can be expressed as
[2(ai|bj)−(aj|bi)](Tijab+TiaTjb). (3.35) The equations for theTiasingle andTijabdouble excitation amplitudes were derived retain- ing the CCSD singles equations, while approximating the doubles equations to be correct through first order with assuming the singles to be zeroth-order parameters . The resulting equations read as
Ωµ1 =hµ1|Hˆ + [ ˆH, T2]|Φ0i= 0 , (3.36) Ωµ2 =hµ2|Hˆ + [F, T2]|Φ0i= 0 , (3.37) where Ω is the CC2 residual and |µni stands for n-fold excited determinants. ˆH is the similarity-transformed Hamiltonian, which is obtained from the original Hamiltonian H as
Hˆ = e−T1HeT1 . (3.38)
The T1-transformed two-electron MO integrals can be expressed in a closed form as (pqˆ|rs) = X
(1−tT1)tp(1+t1)uq(1−tT1)xr(1+t1)ys(tu|xy), (3.39) where t1 is an nb×nb matrix with nb as the size of the basis. The elements of matrixt1 are zero except for its virtual-occupied block, where they are equal to theTµ1 amplitudes.
The simple equations for the doubles amplitudes, Eq. (3.37), can be regarded as the doubles equations of the MP2 method with an effective Hamiltonian, and the doubles amplitudes can simply be expressed as
Tijab = (aiˆ|bj) εi+εj−εa−εb
Dabij , (3.40)
where εp is the corresponding orbital energy. Substituting these into Eq. (3.36) it is obvious that the doubles amplitudes can be calculated on-the-fly, and their storage can be avoided . In practice, the equations for the single excitation amplitudes,
FˆjbTˆijab+ ˆFai+ (εa−εi)Tia= 0 , (3.41) are iterated, where ˆTijab = 2Tijab−Tjiab, and the matrix elements of the similarity-transformed Fock-operator can be calculated as
[2(ck|bj)−(bk|cj)]TkcTibTja. (3.43) Concerning excited states, the LR-CC theory calculates the excitation energies as the eigenvalues of the Jacobi-matrix, which is defined by the derivative of the residual with respect to the cluster amplitudes as
Aµi,νj = ∂Ωµi
∂Tνj = hµ1|[ ˆH, τν1] + [[ ˆH, τν1], T2]|Φ0i hµ1|[ ˆH, τν2]|Φ0i hµ2|[ ˆH, τν1]|Φ0i δµ2ν2εµ2
where εµ2 = −Dabij if τµ2 =a+b+i−j−. Similar to the ground-state problem, the storage of the doubles amplitudes can also be avoided at the solution of the eigenvalue equation of matrix A if an effective Jacobian,
Aeffµ1,ν1(ωCC2) =Aµ1,ν1 −X
εγ2 −ωCC2 , (3.45)
is introduced, where ωCC2 is the CC2 excitation energy . The pseudo-eigenvalue equation of matrix Aeff,
σ(ωCC2,r) =Aeff(ωCC2)r=ωCC2r , (3.46) is only solved in the space spanned by the single excitations, whose coefficients are in- cluded in vector r. The working equations for the σ vector can be derived by straight- forward differentiation of Eq. (3.41) with respect to Tµ1. The structure of the resulting expressions is very similar to the that of the corresponding ground-state terms, thus, the excited-state equations can be solved defining intermediates of the same type. Only the number of intermediates will be higher since each term in Eq. (3.41) including n pieces of Tµ1 amplitudes results inn terms of similar structure upon differentiation. An impor- tant feature of matrix Aeff is that it also depends on the excitation energy of the given excited state, hence the eigenvalue equation cannot be solved simultaneously for all the considered states using the conventional Davidson-type diagonalization techniques.
Another prominent method for excited-state calculations is the ADC(2) approach, which is often called the Hermitian CC2 or iterative MP2 theory for excited states. The derivation of the ADC theory can be approached from the propagator theory [159,163]
or from the so-called intermediate state representation [160,161]. However, we omit the first steps and focus on the analogy between the ADC(2) and CC2 methods solely .
The ground-state ADC(2) correlation energy is simply obtained from MP2 theory as
(ia|jb)(2tabij −tabji), (3.47)
where the first-order amplitudes, tabij, are given in the canonical HF basis as tabij = (ia|jb)
εi+εj−εa−εb = (ia|jb)
Dabij . (3.48)
Utilizing this the first-order Møller–Plesset (MP1) wave function reads as
|ΨMP1i= (1 +T2)|Φ0i . (3.49) The ADC(2) ansatz for the wave function of the excited states is given in the form of
|ΨADC(2)i= (R1+R2)|ΨMP1i, (3.50)
where the spin-coupled single and double excitation operators, R1 and R2, respectively, can be defined similar to Eq. (3.33) withrµ1 andrµ2 as the corresponding coefficients. The excitation energy, being correct up to second-order, is obtained via the diagonalization of the following Hermitian Jacobian:
AADC(2) = Aµ1,ν1 hµ1|[H, τν2]|Φ0i hµ2|[H, τν1]|Φ0i hµ2|[F, τν2]|Φ0i
with elements Aµ1,ν1 =ACISµ1,ν1 +Aµ1,ν1 in the singles-singles block, where
ACISµ1,ν1 =hµ1|[H, τν1]|Φ0i, (3.52) and
Aµ1,ν1 = 1
2(hµ1|[[H, T2], τν1]|Φ0i+hν1|[[H, T2], τµ1]|Φ0i). (3.53) Similarly to the case of LR-CC2 , in practice, the
σ =AADC(2)r =ωADC(2)r (3.54)
eigenvalue problem is recast as a non-linear eigenvalue equation:
σ(ωADC(2),r1) = Aeff(ωADC(2))r1 =ωADC(2)r1, (3.55) where ωADC(2) is the ADC(2) excitation energy, and r1 (r) is a vector composed of the rµ1 (rµ1 andrµ2) coefficients. The benefit is that the resulting equation with the effective Jacobian matrixAeff(ωADC(2)) has to be solved only for therµ1 amplitudes corresponding to single excitations. The elements of the effective Jacobian read explicitly as
1,ν1(ωADC(2)) = Aµ,1ν1 −X
εγ2 −ωADC(2) , (3.56)
with εγ2 =−Dijab.
In the following we briefly collect the working equations required for the implemen- tation of the ADC(2) method in spatial MO basis for a restricted HF reference, because, to the best of our knowledge, they are not published in the literature. Deriving the
expressions corresponding to Eq. (3.55) we arrive at the σia=X
[2(ia|jb)−(ij|ab)]rbj+ (εa−εi)ria (3.57) + 1
[2(kc|jb)−(jc|kb)]rck(2tabij −tabji) + 1
[2(ia|jb)−(ja|ib)]rck(2tcbkj −tcbjk) + 1 2
(ib|ck)(2tbcjk −tbckj)raj + 1
(jb|ck)(2tbcik−tbcki)raj −1 2
− 1 2
sigma vector elements for singlet excitations, while for the triplet case the sigma vector reads as
(ij|ab)rjb+ (εa−εi)rai +1 2
(jc|kb)rkctabji + 1 2
(ja|ib)rcktcbjk (3.58) +1
(ib|ck)(2tbcjk−tbckj)raj +1 2
(jb|ck)(2tbcik−tbcki)rja− 1 2
The required ˆrµ2 intermediates, having different expressions for the two different kinds of spin multiplicities, are obtained using the
εµ2 −ωADC(2) (3.59)
amplitudes. At the end of the iterative procedure the converged ADC(2) solution vector is normalized, which is necessary for the evaluation of transition moments. To achieve this in spatial orbital basis, the amplitudes obtained are divided by the normalization constant
(rabij)2− 1 2
rabijrjiab. (3.60) Then the transition density matrix required for the ground to excited state transition moments can be obtained as
ρpq = hΨMP1|p+q−|ΨADC(2)i=hΦ0|(1 +T2†)p+q−(R1+R2)(1 +T2)|Φ0i
+ hΦ0|T2†p+q−R2|Φ0i. (3.61)