2019 D´avidMester EﬃcientExcited-StateQuantumChemicalApproachesforLargeMolecules

(1)

Faculty of Chemical Technology and Biotechnology George Olah Doctoral School

Efficient Excited-State Quantum Chemical Approaches for Large Molecules

Ph.D. dissertation

D´ avid Mester

Advisor: Prof. Mih´ aly K´ allay

Department of Physical Chemistry and Materials Science MTA-BME Lend¨ ulet Quantum Chemistry Research Group

2019

(2)

Abbreviations

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

HF Hartree–Fock

I/O operation input/output operation

KS Kohn–Sham

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

(5)

Symbols

µ, ν, . . . atomic orbitals

i, j, . . . (quasi-)canonical occupied orbitals i⁰, j⁰, . . . localized occupied orbitals

a, b, . . . (quasi-)canonical virtual orbitals p, q, . . . general (quasi-)canonical orbitals P, Q, . . . auxiliary functions

n_AO, n_aux number of atomic orbitals and auxiliary functions n_occ, n_virt number of occupied and virtual orbitals

C general transformation matrix

c^a_i single excitation CIS coefficients c^ab_ij double excitation CIS(D) coefficients

T_iâ, T_ijâb CC2 ground state single and double amplitudes tâb_ij MP2 ground state amplitudes

r_i^a, r^ab_ij single and double excitation coefficients T₁, T₂ single and double cluster operator R₁, R₂ single and double excitation operator

H,Hˆ normal and similarity-transformed Hamiltonian F,Fˆ normal and similarity-transformed Fockian A,A^eff,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

T_LDF, T_LMO, T_PAO thresholds for the LDF, LMO, and PAO truncation α_X, α_C scaling factors for the exchange- and correlation terms α^OS_C , α^SS_C scaling factors for opposite- and same-spin terms

P local domain

(6)

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 methods 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 moments. 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

(7)

are several actively researched phenomena related to the excited electronic states of molecular systems. For instance, excited states play an important role for photochromic materials, 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 expensive approaches. These results are detailed in Chapter 6.

(8)

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 analytical TDDFT gradients also enable the efficient calculation of excited-state equilibrium structures and other molecular properties [6–9]. The time-dependent Hartree–Fock (TDHF) method [10] 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) [11] 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) [21], 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 approximations 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

(9)

required for the calculations. The simplest approximate TDDFT method is the Tamm–

Dancoff approximation (TDA), which was proposed by Hirata and Head-Gordon [24].

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. [25]

and later enhanced by Whitten [26] and Dunlap and co-workers [27]. 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 [33], it has also been successfully employed for (TDA-)TDDFT [34]. 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 [35] and Mart´ınez and co-workers [36], 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 various methods and accuracy. Such approaches include, for example, the fragment-based TDDFT [37–42] and CIS [43] methods, the frozen-density embedding approximation developed by Neugebauer et al. [44–47], and the density-fragment interaction scheme proposed by Fujimoto and Yang [48]. Many approaches use localized molecular orbitals (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 approaches [56]. 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 [61]. 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 [62] and the auxiliary-density matrix method (ADMM) [63] were developed by Neese and Kjærgaard and co-workers.

(10)

Both schemes reduce the computation time required to calculate the rate-determining exact 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 [67] functional. In ADMM, the exchange contribution is calculated using a smaller auxiliary basis set, and the contribution is extrapolated to the larger basis. 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 [5]. Besides the aforementioned problems, in the well-established 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 [5]. Maitra and co-workers proposed the dressed TDDFT approach, where one double excitation is included at a time [72]. 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 [73]. An even more rigorous polarization propagator approach to dressed TDDFT has been presented by Casida and Huix-Rotllant [74], 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 approximation [75], as well as by Romaniello et al. [76], Sangalli and co-workers [77], and Rebolini and Toulouse [78] 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 [84] following the related multi-coefficient correlation methods of Zhao and Truhlar [85–87]. In Grimme’s B2PLYP (Becke’s exchange with perturbative correction and Lee–Yang–Parr correlation) approach a hybrid KS calculation is performed, and

(11)

the energy is augmented with a second-order perturbation theory (PT2) correction evaluated on the KS orbitals. Several variants of the B2PLYP functional were later proposed with different parametrizations [88–90]. Nonempirical DH functionals, which can be derived 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) [97]

or scaled-opposite-spin (SOS) [98] second-order Møller–Plesset (MP2) correction. The first method in the former category was developed by Chai and Head-Gordon [99], while the combination of SOS-MP2 and DFT was proposed by Goerigk and Grimme [100]. 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 [106], 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 [128].

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 [129]. 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 [130], PBE0-DH, and PBE0-2 [131] DHs, and recently Schwabe and Goerigk successfully combined time-dependent DH DFT theory with spin scaling [132].

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 approach introduced by ˇC´ıˇzek a half a century ago [139]. 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

(12)

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 moments 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. [157], which can also be regarded as an approximate CC2 method [158].

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 relatively easy to implement. The ADC scheme was first derived by Schirmer [159] employing 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 extensive benchmark calculations [166] 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 demanding 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 [170], static polarizability [171], core-valence excitations [172], and excited-state dynamics [173]

were also developed. The ADC method was also combined with the spin-flip [174], the scaled-opposite-spin [169], and the frozen density embedding [175] 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 expenses of CC2 and ADC(2) approaches, which provide results consistently better than

(13)

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 [31], 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- [181] and excited-state [182] 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 contributions. 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 matrix, 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 importance 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 [22] and Craw- ford et al. [210], who generalized the ground-state local CC singles and doubles (CCSD) method of Werner and co-workers [211] to EOM-CCSD. In the local EOM-CCSD method developed by the former authors, the doubles amplitudes were restricted using the information by inspecting the configuration interaction singles (CIS) wave function [22],

(14)

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 [212] 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 theories [218–221]. The chain of spheres exchange [222] 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 transition orbital framework (CorNFLEx) [19] 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 [227].

(15)

Chapter 3 Theoretical background

In this chapter the theoretical background of methods and approximations used in this thesis are outlined. We present briefly the most important equations and considerations, 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 energy, 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)Φ₀(r) = E₀Φ₀(r) (3.2)

form, where F denotes the Fockian, Φ₀ 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.

(16)

The CIS wave function, Ψ^CIS, is constructed as a linear combination of the HF and singly excited determinants:

Ψ^CIS(r) = Φ₀(r) +X

ia

c^a_iΦ^a_i(r), (3.3)

where câ_i is the CIS coefficient and Φâ_i is the excited determinant which is formed as Φâ_i =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

X

ia

hΦ^b_j|H|Φ^a_ii=E_CISc^b_j (3.4) form, with elements of

hΦ^b_j|H|Φ^a_ii= (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

i(r₁)a(r₁) 1

r₁₂j(r₂)b(r₂)dr₁dr₂ , (3.6) where r₁ is the coordinate of electron 1, and r₁₂ = |r₁ −r₂|. 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 ω=E_CIS−E_HF. The final equation using matrix notation can be recast as

Ac=ωc, (3.7)

with the matrix elements

A_ia,jb= (f_ab−f_ij)δ_ijδ_ab+ 2(ia|jb)−(ij|ab). (3.8) The excitation energies can be obtained as the solution of the secular equation by diagonalization of Hermitian matrix A, whose dimension is n_occn_virt, wheren_occ (n_virt) 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.

(17)

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

σ_ia =X

b

f_abc^b_i −X

j

f_ijc^a_j + 2X

jb

(ia|jb)c^b_j−X

jb

(ij|ab)c^b_j . (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⁻|Φ₀i = 0 0 c 0

!

. (3.10)

The oscillator strength (f) in length gauge can be defined by f = 2ω

3 |hΨ^CIS|dˆ|Φ0i|² = 2ω 3

X

α

|hΨ^CIS|dˆ^α|Φ0i|², (3.11) where the sum of squares on the right-hand side is the dipole strength, and ˆd= ( ˆd^x,dˆ^y,dˆ^z) is the dipole operator and its components. The x component of the transition electric dipole moment is obtained as

hΨ^CIS|dˆ^x|Φ₀i=√ 2X

pq

d^x_pqρ_pq =√ 2X

ai

d^x_aiρ_ai, (3.12) where d^x_pq 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 [68] the final equation

(18)

obtained can be written in matrix notation in the form of

A B

−B −A

! c

c⁰

!

=ω c c⁰

!

, (3.16)

where c⁰ is the singles de-excitation vector, and B is the so-called coupling matrix with elements

B_ia,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 thec⁰ 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+Bc⁰ =ωc and −Bc−Ac⁰ =ωc⁰ :

(A−B)(c−c⁰) = ω(c+c⁰), (3.18) and

(A+B)(c+c⁰) =ω(c−c⁰). (3.19) If the second expression is recast forc+c⁰ and this result is inserted into the first equation, the final equation can be read as

(A+B)(A−B)(c−c⁰) = ω²(c−c⁰). (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

b

f_abX

c

f_bc˜c^c_i −X

j

f_ijX

k

f_kjc^a_k+ 4X

jb

(ai|bj)˜c^b_j

+ X

jb

(ab|ij)X

kc

(bc|jk)˜c^c_k−X

jb

(aj|bi)X

kc

(bk|cj)˜c^c_k (3.21) expression, where the ˜c=c−c⁰ shorthand notation is introduced. The transition density matrix used for the calculation of transition properties is obtained as

ρ_pq = 0 c⁰ c 0

!

, (3.22)

(19)

where the excitation and de-excitation vectors can be expressed from the solution of the iterative procedure as

c= ˜c

2 +(A−B)˜c

2ω (3.23)

and

c⁰ = (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[ρ] =T_s[ρ] +J[ρ] +E_XC[ρ] +E_eN[ρ], (3.25) where T_s[ρ] is the exact kinetic energy of a non-interacting system, J[ρ] stands for the classical Coulomb energy, E_XC[ρ] is the exchange-correlation energy, and E_eN[ρ] denotes the potential energy of the electron-nucleus interaction. All contributions in Eq. (3.25), except E_XC[ρ], can be calculated exactly. The performance of DFT methods depends on the accuracy of the XC contribution usually including empirical parameters. In global hybrid density functional theory [231], the XC energy of a ground-state system is expressed as

E_XC^H = (1−α_X)E_X^DFT+α_XE_X^HF+E_C^DFT, (3.26) where E_X^DFT and E_C^DFT are the semilocal exchange and correlation energies, respectively, E_X^HF 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

B_ia,jb^H = 2(ia|jb)−α_X(ib|ja) + (1−α_X)(ia|f_X|jb) + (ia|f_C|jb), (3.28)

(20)

where the DFT exchange term is written in the following form:

(ia|fX|jb) = Z Z

i(r1)a(r1) δ²E_X

δρ(r₁)δρ(r₂)j(r2)b(r2)dr1dr2 , (3.29) where E_X denotes the exchange functional. Similar expression holds for the DFT correlation. 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 systems. 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 A^Hmatrix, 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 B^H 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 [84], E_C^PT2, as E_XC^DH = (1−α_X)E_X^DFT+α_XE_X^HF+ (1−α_C)E_C^DFT+α_CE_C^PT2 (3.30) withα_C as a scaling parameter for the correlation contributions. In the more complicated spin-scaled DHs [103] the opposite-spin (OS) and same-spin (SS) contributions to the PT2 correlation energy, E_C^OS−PT2 and E_C^SS−PT2, are scaled separately by factors α_C^OS and α^SS_C ,

(21)

respectively, as

E_XC^DH = (1−α_X)E_X^DFT+α_XE_X^HF+ (1−α_C)E_C^DFT+α^OS_C E_C^OS−PT2+α^SS_C E_C^SS−PT2 . (3.31) E_XC^DH 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 [129] defined by Eq. (3.30) TDDFT is combined with the CIS(D) method [157]. Eqs. (3.16) or (3.7) are solved with modified A^DH and B^DH 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 [132] 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 [232]. 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 operator acting on the HF ground-state determinant:

Ψ^CC = e^TΦ₀ , (3.32)

where T is the cluster operator, with T =T₁+T₂+T₃+. . ., where T₁ =P

aiT_i^aa⁺i⁻ is the cluster operator of single excitations. For convenience, the T₂ operator together with the higher excitation cluster operators will be denoted in a general form as

T_n = 1 (n!)²

X

µn

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 e^T = e^T¹^+T² operator can be expanded in a Taylor-series in the

e^T = 1 +T +1

2T²+· · ·= 1 +T₁ +T₂ +1

2T₁²+T₁T₂+1

2T₂²+. . . (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.

(22)

In this case, the CC2 correlation energy for a restricted Hartree–Fock reference can be expressed as

∆E_CC2=X

aibj

[2(ai|bj)−(aj|bi)](T_ijâb+T_iâT_j^b). (3.35) The equations for theT_iâsingle andT_ijâbdouble 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 [150]. The resulting equations read as

Ω_µ₁ =hµ₁|Hˆ + [ ˆH, T₂]|Φ₀i= 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^−T¹He^T¹ . (3.38)

The T1-transformed two-electron MO integrals can be expressed in a closed form as (pqˆ|rs) = X

tuxy

(1−t^T₁)tp(1+t1)uq(1−t^T₁)xr(1+t1)ys(tu|xy), (3.39) where t₁ is an n_b×n_b matrix with n_b as the size of the basis. The elements of matrixt₁ are zero except for its virtual-occupied block, where they are equal to theT_µ₁ 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

T_ij^ab = (aiˆ|bj) εi+εj−εa−εb

= (aiˆ|bj)

D^ab_ij , (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 [31]. In practice, the equations for the single excitation amplitudes,

Ωai =X

bjc

(jbˆ|ac) ˆT_ij^cb−X

bjk

(jbˆ|ik) ˆT_kj^ab+X

bj

FˆjbTˆ_ijâb+ ˆFai+ (εa−εi)T_iâ= 0 , (3.41) are iterated, where ˆT_ijâb = 2T_ijâb−T_jiâb, and the matrix elements of the similarity-transformed Fock-operator can be calculated as

Fˆ_jb =X

ai

[2(bj|ai)−(aj|bi)]T_i^a, (3.42)

(23)

Fˆ_ai=X

bj

[2(ai|jb)−(ab|ji)]T_j^b−X

bjk

[2(bj|ik)−(bk|ij)]T_k^aT_j^b+

+X

bjc

[2(bj|ca)−(ba|cj)]T_j^bT_i^c−X

bjck

[2(ck|bj)−(bk|cj)]T_k^cT_i^bT_j^a. (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

!

, (3.44)

where εµ2 = −D^ab_ij 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,

A^eff_µ₁_,ν₁(ω_CC2) =A_µ₁_,ν₁ −X

γ2

A_µ₁_,γ₂A_γ₂_,ν₁

ε_γ₂ −ω_CC2 , (3.45)

is introduced, where ω_CC2 is the CC2 excitation energy [31]. The pseudo-eigenvalue equation of matrix A^eff,

σ(ωCC2,r) =A^eff(ωCC2)r=ωCC2r , (3.46) is only solved in the space spanned by the single excitations, whose coefficients are included in vector r. The working equations for the σ vector can be derived by straight- forward differentiation of Eq. (3.41) with respect to T_µ₁. 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_µ₁ amplitudes results inn terms of similar structure upon differentiation. An important feature of matrix A^eff 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 [158].

The ground-state ADC(2) correlation energy is simply obtained from MP2 theory as

∆E_MP2 =X

ijab

(ia|jb)(2t^ab_ij −t^ab_ji), (3.47)

(24)

where the first-order amplitudes, t^ab_ij, are given in the canonical HF basis as t^ab_ij = (ia|jb)

ε_i+ε_j−ε_a−ε_b = (ia|jb)

D^ab_ij . (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= (R₁+R₂)|Ψ^MP1i, (3.50)

where the spin-coupled single and double excitation operators, R₁ and R₂, 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:

A^ADC(2) = A_µ₁_,ν₁ hµ₁|[H, τ_ν₂]|Φ₀i hµ₂|[H, τ_ν₁]|Φ₀i hµ₂|[F, τ_ν₂]|Φ₀i

!

, (3.51)

with elements A_µ₁_,ν₁ =A^CIS_µ₁_,ν₁ +A^[2]µ1,ν1 in the singles-singles block, where

A^CIS_µ₁_,ν₁ =hµ₁|[H, τ_ν₁]|Φ₀i, (3.52) and

A^[2]_µ₁_,ν₁ = 1

2(hµ₁|[[H, T₂], τ_ν₁]|Φ₀i+hν₁|[[H, T₂], τ_µ₁]|Φ₀i). (3.53) Similarly to the case of LR-CC2 [31], in practice, the

σ =A^ADC(2)r =ω_ADC(2)r (3.54)

eigenvalue problem is recast as a non-linear eigenvalue equation:

σ(ω_ADC(2),r₁) = A^eff(ω_ADC(2))r₁ =ω_ADC(2)r₁, (3.55) where ω_ADC(2) is the ADC(2) excitation energy, and r₁ (r) is a vector composed of the r_µ₁ (r_µ₁ andr_µ₂) coefficients. The benefit is that the resulting equation with the effective Jacobian matrixA^eff(ω_ADC(2)) has to be solved only for ther_µ₁ amplitudes corresponding to single excitations. The elements of the effective Jacobian read explicitly as

A^eff_µ

1,ν1(ω_ADC(2)) = A_µ,₁_ν₁ −X

γ2

A_µ₁_,γ₂A_γ₂_,ν₁

ε_γ₂ −ω_ADC(2) , (3.56)

with εγ2 =−D_ij^ab.

In the following we briefly collect the working equations required for the implementation 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

(25)

expressions corresponding to Eq. (3.55) we arrive at the σ_ia=X

jb

[2(ia|jb)−(ij|ab)]r^b_j+ (ε_a−ε_i)r_i^a (3.57) + 1

2 X

kjbc

[2(kc|jb)−(jc|kb)]r^c_k(2t^ab_ij −t^ab_ji) + 1

2 X

kjbc

[2(ia|jb)−(ja|ib)]r^c_k(2t^cb_kj −t^cb_jk) + 1 2

X

kjbc

(ib|ck)(2t^bc_jk −t^bc_kj)r^a_j + 1

2 X

kjbc

(jb|ck)(2t^bc_ik−t^bc_ki)r^a_j −1 2

X

kjbc

(jb|ck)(2t^ac_jk −t^ac_kj)r_i^b

− 1 2

X

kjbc

(ja|ck)(2t^bc_jk−t^bc_kj)r^b_i +X

bkc

(ab|ck)ˆr_ik^bc−X

cjk

(ij|ck)ˆr^ac_jk

sigma vector elements for singlet excitations, while for the triplet case the sigma vector reads as

σ_ia =−X

jb

(ij|ab)r_j^b+ (ε_a−ε_i)r^a_i +1 2

X

kjbc

(jc|kb)r_k^ct^ab_ji + 1 2

X

kjbc

(ja|ib)r^c_kt^cb_jk (3.58) +1

2 X

kjbc

(ib|ck)(2t^bc_jk−t^bc_kj)r^a_j +1 2

X

kjbc

(jb|ck)(2t^bc_ik−t^bc_ki)r_j^a− 1 2

X

kjbc

(jb|ck)(2t^ac_jk −t^ac_kj)r_i^b

−1 2

X

kjbc

(ja|ck)(2t^bc_jk −t^bc_kj)r_i^b+X

bkc

(ab|ck)ˆr^bc_ik−X

cjk

(ij|ck)ˆr_jk^ac.

The required ˆr_µ₂ intermediates, having different expressions for the two different kinds of spin multiplicities, are obtained using the

r_µ₂ =−X

ν1

A_µ₂_,ν₁r_ν₁

ε_µ₂ −ω_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

c= s

X

ai

(r_i^a)²+X

aibj

(r^ab_ij)²− 1 2

X

aibj

r^ab_ijr_ji^ab. (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Φ₀|(1 +T₂^†)p⁺q⁻(R₁+R₂)(1 +T₂)|Φ₀i

= hΦ₀|p⁺q⁻R₁|Φ₀i+hΦ₀|T₂^†p⁺q⁻R₁|Φ₀i+hΦ₀|T₂^†p⁺q⁻R₁T₂|Φ₀i

+ hΦ0|T₂^†p⁺q⁻R2|Φ0i. (3.61)

2019 D´avidMester EﬃcientExcited-StateQuantumChemicalApproachesforLargeMolecules

Faculty of Chemical Technology and Biotechnology George Olah Doctoral School

Efficient Excited-State Quantum Chemical Approaches for Large Molecules

Ph.D. dissertation

D´ avid Mester

Advisor: Prof. Mih´ aly K´ allay

Department of Physical Chemistry and Materials Science MTA-BME Lend¨ ulet Quantum Chemistry Research Group

2019

Contents

Abbreviations

Symbols

Chapter 1 Introduction

Chapter 2 Overview

Chapter 3

Theoretical background

3.1 CIS and TDHF theories

3.2 TDDFT and TDA-TDDFT theories

3.3 CC2 and ADC(2) theories