FULL PAPER

WWW.C-CHEM.ORG

Statistical Analysis of Electronic Excitation Processes: Spatial Location, Compactness, Charge Transfer, and Electron-Hole Correlation* Felix Plasser,*[a,b] Benjamin Thomitzni,[a] Stefanie A. B€appler,[a] Jan Wenzel,[a] Dirk R. Rehn,[a] Michael Wormit,[a] and Andreas Dreuw[a] We report the development of a set of excited-state analysis tools that are based on the construction of an effective exciton wavefunction and its statistical analysis in terms of spatial multipole moments. This construction does not only enable the quantification of the spatial location and compactness of the individual hole and electron densities but also correlation phenomena can be analyzed, which makes this procedure particularly useful when excitonic or charge-resonance effects are of interest. The methods are first applied to bianthryl with a focus on elucidating charge–resonance interactions. It is

shown how these derive from anticorrelations between the electron and hole quasiparticles, and it is discussed how the resulting variations in state characters affect the excited-state absorption spectrum. As a second example, cytosine is chosen. It is illustrated how the various descriptors vary for valence, Rydberg, and core-excited states, and the possibility of using this information for an automatic characterization of state C 2015 Wiley Periodicals, Inc. characters is discussed. V

Introduction

ground-state orbitals. Specifically, the natural orbitals (NO),[4] the natural transition orbitals (NTO),[5–7] and the attachment/ detachment densities[8] follow this path. The third problem is more subtle and not even its existence is widely acknowledged. However, it is obvious that for a complete description of an excited state all orbital pairings and their relative phases matter. This situation is particularly well seen in the case of symmetric dimers where only a sign change in the configurations differentiates between excitonic and charge-resonance (CR) states.[9–12] In a more general sense, the problem can be traced back to correlation effects between the electron and hole quasi-particles.[12,13] A number of strategies to overcome this problem have been introduced by different authors but their application is still quite limited.[11–22] Aside from a qualitative pictorial representation of the excited state, it is often of interest to proceed to a quantitative analysis. The most obvious possibility for obtaining such descriptors is the computation of physical observables like the oscillator strength or the excited-state dipole moment. The

A plethora of methods for the computation of molecular excited states has been developed in recent years, ranging from methods of highest accuracy to very efficient methods, which allow the treatment of large systems.[1–3] In most cases, the center of attention lies on excitation energies and the character of the wavefunctions is often ignored or analyzed on a quite superficial level, mostly by qualitatively assigning pp , np or Rydberg character based on molecular orbitals. Huge amounts of computational effort are reduced to one number, the excitation energy, while a host of additional information is neglected even though it could be quite easily exploited. This is unfortunate as this information could further our understanding of excitation processes as well as help in pinpointing potential shortcomings of a chosen computational protocol. When excited-state wavefunctions are to be characterized, it is most common to regard a response vector and visualize the corresponding orbitals. While this is a straight-forward approach, there are three significant downsides to it. First, the representation depends on the ground-state orbitals, which may or may not be well suited to describe the state. Second, the results depend on the precise wavefunction expansion chosen and a quantitative comparison between different methods is not possible. Third, whenever more than one pair of orbitals are present, the analysis does not only become cumbersome but also nontrivial phase information may become decisive. A number of strategies are available for tackling the first two problems: typically, the dependence on the response vector is removed by the computation of reduced density matrices (DM), and a subsequent orbital transformation is performed to eliminate any dependence on the

DOI: 10.1002/jcc.23975

[a] F. Plasser, B. Thomitzni, S. A. B€ appler, J. Wenzel, D. R. Rehn, M. Wormit, A. Dreuw Interdisciplinary Center for Scientific Computing, Ruprecht-Karls-University, Im Neuenheimer Feld 368, 69120, Heidelberg, Germany [b] F. Plasser Institute for Theoretical Chemistry, University of Vienna, W€ ahringerstr. 17, 1090, Wien, Austria E-mail: [email protected] *Dedicated to Michael Wormit in grateful memory. Contract grant sponsor: Alexander von Humboldt-Foundation (F.P.), Vienna Scientific Cluster (VSC) School (F.P.), Heidelberg Graduate School of Mathematical and Computational Methods for the Sciences (J.W., D.R.R., and S.A.B.), and Landesgraduiertenf€ orderung Baden-W€ urttemberg (S.A.B.) C 2015 Wiley Periodicals, Inc. V

Journal of Computational Chemistry 2015, 36, 1609–1620

1609

FULL PAPER

WWW.C-CHEM.ORG

advantage of this approach is that physically well-defined properties are computed. However, the information gained is often rather limited as many properties of interest are not directly linked to a physical observable. Another interesting approach for obtaining quantitative information about an excited state derives from various quantities emerging in the aforementioned transformation formalisms. The NO eigenvalues are used to compute an effective number of unpaired electrons[23–26] while the NTO eigenvalues are used to give the NTO participation ratio or collectivity index,[16] which measures excited-state correlation.[11,27] The electronic promotion number arising in the attachment/detachment analysis[8] can be used to quantify orbital relaxation.[28] And, finally, the squared norm of the 1TDM, which arises as normalization factor in the aforementioned electron-hole correlation analysis, serves as a universal measure of single-excitation character.[13,21,29] Furthermore, the different densities discussed above can be subjected to population analyses to divide them into different atomic contributions.[13,30] This provides specific information about the involvement of the different atoms with the downside of being dependent on the specific population analysis scheme chosen. Finally, a number of measures for quantifying charge transfer have been suggested.[31–36] While these can provide specific information of interest, the interpretation of the results is hampered by the lack of a universal underlying physical picture. In this work, a new type of quantitative excited-state analysis is introduced using a formalism based on the physical picture of an exciton wavefunction, i.e., a two-body function describing the correlated motions of the electron and hole quasi-particles.[12] While this formalism originates in the description of periodic solids, we show that the same concepts are also advantageous for excited states in molecules. The square of the exciton wavefunction signifies the joint probability distribution of the electron-hole pair, and we suggest analyzing this function by standard statistical procedures. The quantities obtained allow to answer such questions as to where the charge distributions are located, how diffuse or compact they are, and whether or not excitonic correlation effects play a role. Special focus is laid on charge separation where two distinct definitions of the electron-hole separation and their physical interpretations are given. Two molecules are used to exemplify different aspects of our methods (Fig. 1). First, bianthryl is chosen as a molecule with a particularly pronounced spectroscopic signature of CR effects.[37,38] In the following text, we do not only quantify these and their modulations with changing interring torsion, we also provide a phenomenological interpretation: CR effects are identified with negative correlations of the electron and hole quasiparticles, i.e., the electron and hole avoid each other dynamically. Furthermore, the excited-state absorption spectrum of bianthryl is computed and interpreted in terms of the various excitonic and CR states. Second, cytosine is chosen as a prototypical heteroaromatic molecule. Here, we show how the different descriptors vary for valence, Rydberg, and core excitations and discuss how the results could be used for an 1610

Journal of Computational Chemistry 2015, 36, 1609–1620

Figure 1. Molecules studied in this work: (a) bianthryl and (b) cytosine.

automatic assignment of these state characters without visual inspection of orbital shapes.

Quantitative Analysis of Electronic Excitation Processes In this section, the different analysis methods are introduced and their physical and mathematical background is discussed. The physical picture underlying our approach is the analysis of the two-body exciton wavefunction vexc ðrh ; re Þ, which describes the correlated motion of the electron and hole quasi-particles. First, this function is constructed from the 1TDM and its one-body and two-body multipole moments are computed. Second, using this foundation, a statistical analysis of the exciton probability distribution is carried out. Furthermore, the results are interpreted in terms of an expansion in NTOs, and finally, analogous quantities for the difference density matrix (1DDM) are given. The exciton wavefunction and its multipole moments Within many-body Green’s function theory, the exciton is described by the electron-hole amplitudes, which are given by the one-particle transition density matrix (1TDM) c0I between the ground state and the excited state I of interest.[12,39] Using the same idea, we compute an exciton wavefunction vexc by identifying it with the 1TDM of a quantum chemical excited state calculation, i.e., vexc ðrh ; re Þ5c0I ðrh ; re Þ:

(1)

The 1TDM is in turn defined as ð c0I ðrh ; re Þ5N U0 ðrh ; r2 ; :::; rN Þ3UI ðre ; r2 ; :::; rN Þdr2 ; :::; drN ; (2) where U0 and UI are the ground and excited-state wavefunctions. ri denotes the spatial and spin coordinates of the i-th electron: ri 5ð~ x i ; si Þ with ~ x 5ðx; y; zÞT being the coordinates in standard three-dimensional Euclidian space. The matrix representation of the 1TDM in an orbital basis set fvl g is given as †

0 ^l a^m jUI i; D0I lm 5hU ja †

(3)

where a^l and a^m are the standard one-particle creation and annihilation operators. This expression can be used to expand the 1TDM in the following form: WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

c0I ðrh ; re Þ5

X

D0I lm vl ðrh Þvm ðre Þ:

(4)

By virtue of this construction, the exciton wavefunction is a readily available quantity, which can be analyzed by standard quantum-chemical methods. It is not only amenable to population analysis,[13] but also expectation values may be computed[12] according to ^ ^ hvexc jOjvexc i : hOi5 hvexc jvexc i

(5)

The denominator of Eq. (5) contains the squared norm X of the exciton wavefunction, a quantity with the physical interpretation of measuring the single-excitation character of an excitation.[13,21,28,29] It is explicitly given as X5hvexc jvexc i5

ðð

c0I ðrh ; re Þ2 drh dre :

(6)

ð0Þ

and specifically Mn 5S is the overlap matrix in the AO basis. While these multipole moments do not possess a clear intuitive meaning themselves, they can be combined to such quantities as will be shown in the next section. Statistical description of the exciton wavefunction The square of the exciton wavefunction signifies the probability distribution of the exciton, i.e., vexc ðrh ; re Þ2 gives the probability that the hole is located at position rh while the electron is at re. From a mathematical point of view, this is a bivariate probability distribution, which can be characterized by standard statistical methods. From a physical point of view, the resulting quantities give information about the spatial location and compactness of the charge distributions, as well as about charge transfer and correlation effects. In this section, we will define various such quantities, explain their construction using the above defined multipole moments, and discuss their meaning. The discussion is started with the mean separation vector

Insertion of Eq. (4) into this expression yields the matrix representation[11,13] of X with respect to an orbital basis set X5trðDI0 SD0I SÞ;

(7)

where S is the overlap matrix between the basis functions and DI0 is equal to the transpose of D0I. For the following analysis, various statistical moments will be computed, and as will be seen below, the common building blocks for these are expectation values of polynomials of the electron and hole coordinates. The operator O^ is of the form ^ O5

X

akl ðxhk xel 1yhk yel 1zhk zel Þ5

k;l

X X

akl nkh nle

(8)

~h!e 5h~ d x e 2~ x h i5h~ x e i2h~ x h i;

~h!e j, which give the distance and its absolute value dh!e 5jd between the centroids of the hole and electron densities. This measure can be used to quantify charge separation in simple charge-transfer states. However, it vanishes by construction for centrosymmetric systems (cf. Refs. [32, 40, 41]) and dynamical charge-separation effects are neglected. Moving to second-order multipole moments, the first property we introduce is the hole size rh defined as the rootmean-square (rms) deviation of the hole position operator

n2fx;y;zg k;l

rh 5 where the akl are arbitrary coefficients and the symbol n is used to denote a generic Cartesian coordinate (x, y, or z). The expectation value of this operator may be split according to ^ hOi5h

X X

X X

akl nkh nle i5

n2fx;y;zg k;l

akl hnkh nle i

(9)

n2fx;y;zg k;l

yielding individual terms of the form hnkh nle i5

1 X

ðð

nkh nle c0I ðrh ; re Þ2 drh dre

hnkh nle i5

1  I0 ðkÞ 0I ðlÞ  tr D Mn D Mn ; X

(11)

ðkÞ

where Mn refers to the kth-order multipole matrix for coordinate n, a quantity which is routinely available in quantum chemical calculations. Its components are given as ð ðkÞ Mn;lm 5 vl ðrÞnk vm ðrÞdr

(12)

(14)

P

2 2 n2fx;y;zg hnh i2hnh i 5

P

  ð2Þ tr DI0 Mn D0I S

n2fx;y;zg

which are interpreted as multipole moments of the exciton probability distribution v2exc . As explained in Ref. [12], these multipole moments can be evaluated in analogy to Eq. (7)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h~ x 2h i2h~ x h i2 :

This value is small when the excitation process is initiated from one compact orbital while it increases if this orbital is more diffuse or if several spatially distributed orbitals are involved. Using Eq. (11), the explicit form of rh is readily obtained: r2h 5

(10)

(13)

X

12 0  ð1Þ tr DI0 Mn D0I S A : 2@ X

(15)

At this point, it is worth noting that r2h corresponds to the trace of the covariance matrix of the vector ~ x h. A more complete picture of the hole distribution could be obtained by computing the complete matrix. However, in this work, we restrict ourselves to scalar quantities for simplicity. In analogy to Eq. (14), the size of the electron density is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi re 5 h~ x 2e i2h~ x e i2 :

(16)

An interesting application of re is measuring the diffuseness of the acceptor orbital in a Rydberg state, an option which will Journal of Computational Chemistry 2015, 36, 1609–1620

1611

FULL PAPER

WWW.C-CHEM.ORG

be explored below. The three quantities dh!e , rh, and re are depicted in Figure 2. Moving to mixed multipole moments reveals information about correlated electron-hole motions. These correlations give rise to two central phenomena: CR interactions and excitonic effects, which are both notoriously difficult to understand in the standard molecular-orbital picture. We have previously defined the exciton size or dynamic chargeseparation distance qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dexc 5 hj~ x e 2~ x h j2 i

(17)

as the rms separation of the electron and hole positions. This measure does not only provide the expected results in the case of standard charge-transfer states, a case in which also dh!e can be applied, but proved to be very helpful in the characterization of CR and excitonic interactions as well.[12] To analyze electron-hole correlation, we use two quantities from bivariate statistics: the covariance and correlation coefficient. Although these quantities are better known in the context of linear regression analysis, they can be used effectively in the present context as well. The covariance between the hole and electron position vectors is defined as COVð~ x h; ~ x e Þ5hð~ x h 2h~ x h iÞ  ð~ x e 2h~ x e iÞi5h~ xh  ~ x e i2h~ x h i  h~ x e i: (18) Here, again a scalar is computed rather than the full crosscovariance matrix between the two vectors. The covariance itself does not possess an immediate intuitive meaning and is therefore usually normalized against the standard deviations, yielding the Pearson correlation coefficient. In analogy, we define the correlation coefficient of the electron-hole pair as Reh 5

COVð~ x h; ~ x eÞ : rh re

(19)

The values of Reh range from 21 to 1 (a fact which can be seen for example by application of the Cauchy–Schwarz inequality). Positive correlation, as indicated by Reh > 0 points to a concerted motion of the electron and hole quasi-particles, while negative correlation means that they avoid each other dynamically. A value of Reh  0 indicates that the hole and electron behave independently. This simple statistical transformation offers a completely new perspective onto electron-hole correlation as will be outlined in more detail below. However, a restriction of the present formalism is that Reh only allows the quantification of linear correlation effects. A more elaborate analysis scheme including higher multipole moments would be required to quantify more involved types of correlation, e.g., the interplay between negative short-range and positive long-range correlations postulated for excitons in conjugated polymers.[42] In the above discussion, two measures of the electron-hole separation were given, the mean dh!e and the root-meansquare dexc electron-hole separation distances. The former marks the distance between the charge centroids (Fig. 2), 1612

Journal of Computational Chemistry 2015, 36, 1609–1620

Figure 2. Graphical representation of the hole (red) and electron (blue) densities together with the hole size rh, the electron size re, and the distance between their centroids dh!e .

whereas the latter corresponds to a dynamical charge separation, which is identified with the concept of an exciton size (see e.g., Ref. [43]), and we will now shed light on their numerical relation. The two quantities are connected by the following equation 2 2 dexc 5dh!e 1r2h 1r2e 223COVð~ x h; ~ x e Þ;

(20)

which can be verified by inserting the respective definitions. When compared to the distance between the charge centroids, the actual dynamic charge separation is enhanced by the sizes of the electron and hole distributions (a fact also discussed in Ref. [32]), as well as by correlation effects represented by the covariance. Inserting the definition of Reh, Eq. (19), and considering that this value ranges from 21 to 1 leads to the inequality 2 2 2 dh!e 1ðrh 2re Þ2  dexc  dh!e 1ðrh 1re Þ2

(21)

which shows the range of dexc values that can be obtained for 2 a given static charge distribution as characterized by dh!e , rh, and re. Specifically, it follows that dh!e  dexc ;

(22)

i.e., the dynamic charge-separation distance can never be smaller than the distance between the charge centroids “Another way to understand this relation is by realizing that the quadratic mean of a distribution is always greater or equal to the arithmetic mean.” Interpretation in terms of NTOs To gain intuitive insight into the above-defined descriptors, it is helpful to discuss the equations in terms of molecular orbitals. For this purpose, we use the NTOs, which are defined by a singular value decomposition of the 1TDM[5–7] pffiffiffiffiffi pffiffiffiffiffi D0I 5U diagð k1 ; k2 ; . . .ÞVT :

(23)

I0 Here, the pairs of hole/particle NTOs w0I i /wi are contained in the respective columns of the U and V matrices, whereas the weights of the configurations are represented by the ki with

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

P X5 i ki .[11] Using the NTOs, the 1TDM can be expanded in the diagonal form vexc ðrh ; re Þ5

X pffiffiffiffi I0 ki w0I i ðrh Þwi ðre Þ:

(24) I0 vexc ðrh ; re Þ5w0I 1 ðrh Þw1 ðre Þ:

i

Then, the general Eq. (11) for the multipole moments becomes hnkh nle i5

1 X pffiffiffiffiffiffiffiffi 0I ki kj hwi ðrh Þjnkh jw0I j ðrh Þi3 X i;j 1 X pffiffiffiffiffiffiffiffi 0I;k I0;l l I0 ki kj Mn;ij Mn;ij hwI0 i ðre Þjne jwj ðre Þi5 X i;j

(25)

1X 1 X pffiffiffiffiffiffiffiffi 0I;k I0;l I0;l ki M0I;k ki kj Mn;ij Mn;ij : n;ii Mn;ii 1 X i X i6¼j

(26)

The first term refers to the shapes of the NTOs of the individual NTO pairs and can be qualitatively understood by visual inspection. By contrast, the second term connects the different NTOs. It is therefore sensitive to their relative phases and fundamentally difficult to understand intuitively. This second term becomes important whenever the off-diagonal elements of the NTO multipole matrices do not vanish. This is the case if the orbitals occupy similar regions of space, which is usually the case in highly symmetric or delocalized systems. Under these conditions, the 1TDM-based approach is clearly advantageous as alternative strategies to obtain analogous information[9,40] are quite involved. However, if the orbitals are localized, the problem disappears and a visual analysis becomes more straightforward. At this point, it is interesting to compare the formalism presented above to theories that only include diagonal NTO contributions.[32,35] These can lead to the same results in cases of independent NTO pairs, but cannot account for cases where the off-diagonal terms are important. Furthermore, a formal shortcoming of leaving out the off-diagonal terms lies in the fact that the results depend on the resolution of degenerate NTOs leading to potential discontinuities with respect to changes in the molecular geometry. We also want to point out that the problem of cross terms exists only for mixed multipole moments. If for example, only the hole coordinate is considered the expectation value takes up the simple form hnkh i5

1X ki M0I;k n;ii ; X i

(28)

In this case, the multipole moments are simply given as a product of the respective integrals of the hole and electron NTOs k 0I I0 l I0 hnkh nle i5hw0I 1 ðrh Þjnh jw1 ðrh Þihw1 ðre Þjne jw1 ðre Þi:

(29)

As a consequence, the covariance and therefore also Reh vanish:

I0;l where M0I;k n;ij and Mn;ij denote the multipole integrals over the hole and particle NTOs. To analyze this equation in more detail, it is helpful to separate the diagonal and off-diagonal parts:

hnkh nle i5

An interesting situation is obtained in the singleconfigurational case, i.e., if only one of the NTO amplitudes is not vanishing (see also Ref. [27])

(27)

because the fwI0 i g form an orthonormal set. The multipole moment of the hole is simply obtained as the weighted sum of the multipole moments of the individual hole NTOs.

I0 I0 x e Þ5hw0I xh  ~ x e jw0I COVð~ x h; ~ 1 ðrh Þw1 ðre Þj~ 1 ðrh Þw1 ðre Þi2 I0 x h jw0I x e jwI0 hw0I 1 ðrh Þj~ 1 ðrh Þi  hw1 ðre Þj~ 1 ðre Þi50 (30)

This shows that a nonvanishing correlation coefficient Reh can only be obtained for a multiconfigurational excited-state wavefunction (given a closed shell ground state). For measuring the multiconfigurational nature of the excited state and thus the electron-hole correlation, the NTO participation ratio was introduced previously[11] P ð ki Þ2 X2 PRNTO 5 Pi 2 5 P 2 : i ki i ki

(31)

If PRNTO is 1, then it strictly follows that Reh vanishes, while higher values of PRNTO are in general associated with correlated motions between the electron and hole quasi-particles.

Extension to the attachment/detachment analysis of the difference density matrix As an alternative to the 1TDM, it is possible to analyze the one-particle difference density matrix (1DDM), which is constructed by subtracting the density matrix of the ground state D00 from the one of the excited state DII. Subsequently, the attachment/detachment analysis of Head-Gordon et al.[8] is performed. This procedure starts with a diagonalization of the difference density matrix WT ðDII 2D00 ÞW5diagðj1 ; j2 ; . . .Þ

(32)

where the ji are the eigenvalues of the 1DDM, and W is the matrix containing the eigenvectors, termed natural difference orbitals (NDOs).[13] The NDOs can be used to visualize the excitation process similarly to the NTOs[44] while also accounting for many-body and orbital-relaxation effects.[28] The detachment density matrix DD is obtained by only considering the negative eigenvalues of the 1DDM di 5minðji ; 0Þ

(33)

and a back-transformation to the initial orbital basis Journal of Computational Chemistry 2015, 36, 1609–1620

1613

FULL PAPER

WWW.C-CHEM.ORG

DD 5Wdiag ðd1 ; d2 ; . . .ÞWT :

(34)

The attachment density DA is constructed from the positive eigenvalues in an analogous fashion ai 5maxðji ; 0Þ

(35)

DA 5Wdiag ða1 ; a2 ; . . .ÞWT :

(36)

The promotion number p52

X

di 5

X

i

ai

(37)

i

corresponds to the spatial integral over the attachment or detachment densities. It counts the total number of electrons rearranged and is enhanced by double-excitation character and orbital relaxation.[28] In the present context, p serves as a normalization factor in analogy to X. In contrast to the 1TDM, it is not possible to construct a correlated two-body function from the 1DDM.[13] It is therefore not possible to translate the whole formalism into the 1DDM picture and only a reduced set of descriptors is available. First, the distance vector between the charge centers is defined as

~D!A 5 d

  ð1Þ tr DA Mn p

2

  ð1Þ tr DD Mn p

:

(38)

Furthermore, the variance of the attachment and detachment densities are given as

r2A=D 5

  ð2Þ X tr DA=D Mn n2fx;y;zg

p

12 0  ð1Þ tr DA=D Mn A : 2@ p

(39)

For the configuration-interaction singles method, it holds that ~h!e ; rD 5rh , and rA 5 re. However, for correlated ~D!A 5d d methods, these relations do not apply anymore as the difference density contains additional contributions from manybody and orbital-relaxation effects.[13]

Computational Details Geometry optimizations of bianthryl were performed at the MP2/SV(P) level using the resolution-of-the-identity approximation (RI).[45,46] A potential curve was computed by fixing the interring torsion (defined as the average of the four equivalent CCCC interring torsion angles) and optimizing the remaining coordinates in the ground state. Bianthryl excitation energies were computed at the RI-ADC(2)/SV(P) level.[47–51] For these calculations, the D2 symmetry group was used. Oscillator strengths for excited-state absorption were computed with respect to the lowest excited state (11 B2 ). For this purpose, the transition dipole moments between excited states were obtained via the intermediate state representation (ISR)[52] using the implementation described in Refs. [50, 53]. The cytosine molecule was optimized at the MP2/SV(P) level and excitation energies were computed at the ADC(2)/aug-cc1614

Journal of Computational Chemistry 2015, 36, 1609–1620

pVDZ level.[54] Core excitations of cytosine were calculated at the ADC(2)-x level in combination with the core-valence separation (CVS) approximation[55–59] using the Cartesian 6D/10F version of the 6-31111G** basis set,[61,62] a combination which has been shown to provide accurate results when compared to experimental data.[58–60] Geometry optimizations at the RI-MP2 level were generally carried out using the Turbomole 6.3 program package.[63] All excited-state ADC calculations were performed with a developmental version of the Q-Chem 4.2 program package,[64,65] and the methods implemented here are scheduled to be released within Q-Chem 4.3. Post processing of the data was carried out using the TheoDORE 1.0 program.[66] Aside from the quantities defined above, also a previously defined measure for the total charge separation was computed. First, the charge-transfer numbers XAB were constructed (see also Ref. [16]). For two molecular fragments A and B, these are given as[13] XAB 5

1 X X 0I 0I ½ðD SÞlm ðSD0I Þlm 1D0I lm ðSD SÞlm  2 l2A m2B

(40)

where the sums run over all atomic basis functions l,m on the respective fragments. Second, to compute the total amount of charge separation,[11] a summation over all the off-diagonal elements of the X-matrix is performed, yielding xCT 5

1X XAB : X B6¼A

(41)

This value ranges from 0 to 1 where xCT 5 0 points to a locally excited state or Frenkel exciton while xCT 5 1 indicates complete charge separation.

Results and Discussion CR interactions in bianthryl Over the last decades, the bianthryl molecule (Fig. 1a) has been subjected to a number of experimental[37,38,67–69] and computational[70–72] studies, and in particular the formation of CR states has been examined. A CR state is formed through a superposition of charge-transfer states of opposite directionality, which leads to a state without any net charge shifts or dipole moments in spite of its intrinsic charge-separated character.[9,11,73] This phenomenon is highly challenging to grasp in the MO picture, as a meticulous analysis of orbital phases is required,[9,40] but can be readily understood using the tools described above. First, the ground-state (S0) minimum is optimized and a potential energy curve is computed along the interring torsion coordinate, which has been characterized as the decisive reaction coordinate previously.[67,71] The unconstrained optimization of the ground state leads to a torsion angle of 90 , i.e., perpendicular anthryl units. This result can be rationalized by considering that this geometry minimizes the steric repulsion, which outweighs the loss of p-conjugation energy. The situation is different for the excited states for which the effect of WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

Figure 3. Relaxed scan around the central torsion coordinate of bianthryl computed at the RI-ADC(2)/SV(P) level of theory: (a) energies (eV) relative to the ground-state minimum, (b) charge-transfer measures, and (c) exciton sizes. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

p-conjugation is obviously enhanced, as presented in Figure 3. The four excited states shown exhibit rather flat potentials up to higher torsion angles and the first excited state (11 B2 ) even shows a clear minimum, stabilized by 0.10 eV, present at a torsion angle of 120 . Following Refs. [37, 38], we next address the CR character of the different excited states and its modulation with respect to interring torsion. Due to the symmetry of the system, the canonical MOs as well as the NTOs are all evenly distributed between the two anthryl units and no charge transfer can be discerned by simply looking at them. The information of interest lies in the off-diagonal resonance contributions between the different NTOs as given by the second term in Eq. (26). To recover these terms, we use the charge-transfer measure [xCT, Eq. (41)] and the exciton size [dexc, Eq. (17)]. The xCT values

FULL PAPER

are presented in Figure 3b. At the ground-state minimum (90 torsion), there is a clear separation between the excitonic states 11 B2 and 21 A, which are characterized by a xCT value close to zero (xCT 5 0.03) and the CR states 31 A and 21 B2 with xCT around 0.94. Once the rings are twisted, these states interact and the charge-transfer measures converge toward a value of 0.5, i.e., an even mixture between CR and excitonic character. The exciton size (Fig. 3c) dexc follows a very similar trend. When the rings are perpendicular, a clear separation ˚ and the CR between the local states with dexc around 4.0 A ˚ states with dexc around 5.6 A is apparent. And again these quantities converge toward an intermediate value with torsion. To elucidate the differences between the excitonic and CR states, two additional properties are discussed, the electronhole correlation coefficient [Reh, Eq. (19)] and the NTO participation ratio [PRNTO, Eq. (31)]. The Reh values are plotted in Figure 4a. These show that the excitonic states are characterized by a positive correlation between the electron and hole while the CR states show a negative value. At the perpendicular geometry, these are Reh 5 0.36 and Reh 5–0.33. The phenomenological interpretation is that in the case of excitonic states, the electron and hole are either both on the first or both on the second fragment, while they are always on opposing fragments in the case of CR states. The concerted motion in the first case leads to positive correlation, while the dynamical repulsion in the second case leads to negative correlation. The absolute values of Reh at this geometry are about a third of the possible maximum of 1. This shows that there is of course no perfect linear correlation as can be understood by the fact that the chromophores are rather large and positioned close together. With increasing torsion, the Reh values move toward zero in agreement with the previous results. As discussed above, correlation effects are associated with the multiconfigurational character of the excited state, a feature, which can be monitored by the NTO participation ratio (PRNTO) as computed from the NTO singular values [Eq. (31)].[11] The results are presented in Figure 4b. At 90 torsion, the PRNTO value is slightly above two for all four states considered. This corresponds to the idealized situation of coupling between two independent chromophores,[74] where the delocalized state is constructed as a superposition of two independent configurations (see also Refs. [11, 12]). The torsion results in a lowering of PRNTO and at 135 this value is smaller than 1.1 for all four states considered. As shown in Eq. (30), this is expected to coincide with a lowering of Reh, which is in fact the case. However, an important difference between the two quantities is that only Reh allows to judge whether positive or negative correlations exist while PRNTO only allows to gauge the general multiconfigurational character. Finally, we discuss the physical implications of the rather abstract statistical descriptors. For this purpose, the excitedstate absorption (ESA) spectrum of bianthryl out of its first excited state (11 B2 ) is computed and the results are compared to time-resolved experiments. In Figure 5, the ESA energies DDE and oscillator strengths of the four lowest states with nonvanishing oscillator strength to 11 B2 are plotted. These amount to the 21 A and 31 A states discussed above, as well as, Journal of Computational Chemistry 2015, 36, 1609–1620

1615

FULL PAPER

WWW.C-CHEM.ORG

At 90 the 11 B2 ! 31 A, transition proceeds from a locally excited into a completely charge-separated state, whereas after the torsion both possess partial CT character, which allows for a direct transition between them. At higher torsion angles also the 11 B2 ! 21 A transition is activated, which is related to mixing between the 21 A and 31 A states. In summary, the following picture emerges: at the perpendicular geometry, the S1 state is of locally excited character and nonvanishing ESA is only possible with respect to the higher-lying 21 B1 and 21 B3 states. With the torsion, the local and CT states mix, which leads to enhanced CR character for the S1 state. This in turn facilitates ESA into the states of A symmetry, which are reached by lower transition energies. The flat S1 potential along the torsion (Fig. 3a) suggests that the two structures can be readily interconverted explaining why they have the same transient lifetime. This qualitative picture nicely agrees with the general spectroscopic signature reported by Asami et al.[38] A more quantitative comparison of

Figure 4. Relaxed scan around the central torsion coordinate of bianthryl computed at the RI-ADC(2)/SV(P) level of theory: (a) electron-hole correlation coefficient and (b) NTO participation ratio. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

the higher-lying 21 B3 and 21 B1 states. At the ground-stateminimum geometry, there are two strong transitions into the 21 B3 and 21 B1 states with an ESA energy of about 1.432 eV (11,500 cm–1). Note that due to the enhanced symmetry at this geometry these states are degenerate and transform as an E irreducible representation of the D2d symmetry group. As these transitions are polarized perpendicular to the inter-ring CAC bond, they can safely be identified with the locally excited band of Ref. [38] located at 9800 cm–1. As the torsion angle is increased, the absorption strengths into these higherlying states are steadily lowered. However, with the torsion, the transitions to the A states located at lower energy, which are polarized along the inter-ring CAC bond, come into play. First, the 11 B2 ! 31 A transition is activated, peaking at a torsion angle of 125 with an oscillator strength of 0.077, and an excited-state absorption energy of 0.339 eV (2700 cm–1). This energy is lower than the one given in Ref. [38] for the partial CT band (6700 cm–1). However, it is certainly reasonable that it would increase in energy if the S1 geometry were relaxed and if solvation effects were taken into account. The modulation in transition strengths can be understood by revisiting Figure 3: 1616

Journal of Computational Chemistry 2015, 36, 1609–1620

Figure 5. Excited-state absorption of bianthryl for different torsion angles: (a) relative energies and (b) oscillator strengths, both with respect to the 11 B2 state. B2 ! A transitions are polarized along the inter-ring C-C bond, while B2 ! B1 and B2 ! B3 transition are polarized perpendicular to it. The four lowest energy states with non-vanishing oscillator strength are considered. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 1. Analysis of the first nine singlet excited states of cytosine at the ADC(2)/aug-cc-pVDZ level of theory: excitation energies, oscillator strengths, various descriptors computed from the transition density matrix (see text), and state character. State S1 S2 S3 S4 S5 S6 S7 S8 S9

DE (eV)

f

X

PRNTO

dh!e (A˚)

4.565 4.787 5.194 5.325 5.534 5.655 5.754 5.890 5.944

0.055 0.002 0.001 0.006 0.148 0.000 0.030 0.034 0.009

0.85 0.83 0.83 0.88 0.85 0.86 0.86 0.88 0.88

1.095 1.010 1.025 1.024 1.212 1.016 1.080 1.392 1.458

0.96 1.90 1.76 2.07 1.09 1.56 2.83 2.38 1.86

rh[a] (A˚) 2.00 1.73 1.64 2.01 1.94 1.31 1.61 2.10 2.18

(0.72) (0.40) (0.41) (0.74) (0.72) (0.37) (0.40) (0.71) (0.74)

re[a] (A˚) 1.99 1.98 2.00 3.57 2.12 2.15 3.71 3.90 4.08

(0.84) (0.83) (0.83) (1.86) (0.91) (0.87) (1.89) (1.72) (1.73)

dexc (A˚)

˚ 2) COV (A

Reh

Character[b]

2.92 3.21 3.08 4.51 3.05 2.91 4.78 4.74 4.64

0.18 0.13 0.14 0.40 0.07 0.15 0.77 1.43 1.71

0.046 0.037 0.042 0.056 0.017 0.053 0.129 0.174 0.192

p2p nmix 2p nmix 2p p2Ryd. p2p nO 2p nO 2Ryd. p2Ryd. p2Ryd.

[a] z-component given in parentheses. [b] nmix–mixed n-orbital, nO–n-orbital located on oxygen, Ryd.–diffuse Rydberg orbital.

the excitation energies would require the inclusion of structural relaxation, solvation effects, and possibly dynamics and is out of the scope of this work.

Valence, Rydberg, and core-excited states in cytosine As a second example, the cytosine molecule is discussed as a prototypical small heteroaromatic molecule possessing a wide variety of excited states of different character.[75–78] In the following text, we consider valence states of pp and np character, compare them to Rydberg states and later proceed to core-excited states. Specific signatures of these diverse types of excited states are highlighted. Furthermore, numerical examples for the different statistical moments and their relations are given. In addition, the difference density matrix is analyzed.

Figure 6. The two primary NTO pairs of the S8 Rydberg state of cytosine. The number on the arrow indicates the weight ki.

Geometry optimization of cytosine at the RI-MP(2) level leads to a not perfectly planar structure exhibiting some pyramidalization at the amino group and no symmetry is present for this molecule. However, for the following analysis, it is favorable to realize that this molecule is almost planar and we define the z-coordinate as the one pointing out of this plane. The first nine excited states of cytosine are computed at the ADC(2)/aug-cc-pVDZ level. Excitation energies, oscillator strengths, and descriptors of the 1TDM are presented in Table 1. It is observed that the computed descriptors allow for a well-defined separation between states of different character as given in the right column. First, the Rydberg states are clearly separated from the rest by large exciton sizes (dexc > ˚ ) and in particular electron sizes (re > 3:5 A ˚ ). Second, the 4.5 A states involving n-orbitals are characterized by comparably ˚ ) where in particular the out-ofsmall hole sizes (rh < 1:75 A ˚ plane z-component  0.4 A as shown in parenthesis in Table 1. Finally, the states not fulfilling these conditions can be characterized as pp states. The correlation coefficients Reh range from 0.017 for S5 (pp ) to 0.192 for S9 (Rydberg). These values are rather close to zero, which is in agreement with the general expectation that electron-hole correlation, as defined above, does not play an important role in small molecules and supports the general notion of ignoring these effects in such systems. It should, however, be stressed that this analysis does not pertain to dynamic correlation effects, which are certainly important in cytosine for obtaining the correct energetics.[75,76] Furthermore, it is observed that the Reh values are all above zero, which means that there is in general a positive correlation between the electron and hole indicating a concerted motion deriving from their mutual attraction. As pointed out above, a loose connection between Reh and PRNTO exists. The present investigation confirms this but also shows that this is not a strict monotonous relation as Reh depends on the orbital shapes as well. The largest Reh values (0.174, 0.192) are present for the higher-energy Rydberg states, which also possess higher PRNTO values (1.392, 1.458). The enhanced correlation in this case probably derives from the presence of more polarizable orbitals. In Figure 6, the two primary NTO pairs for the S8 state (Reh 5 0.174) are shown as an Journal of Computational Chemistry 2015, 36, 1609–1620

1617

FULL PAPER

WWW.C-CHEM.ORG

Table 2. Analysis of the first nine excited states of cytosine at the ADC(2)/aug-cc-pVDZ level of theory: descriptors computed from the difference density matrix. State S1 S2 S3 S4 S5 S6 S7 S8 S9

p

dD!A (A˚)

rD (A˚)

rA (A˚)

1.39 1.67 1.64 1.56 1.37 1.50 1.79 1.58 1.51

0.56 0.86 0.79 1.26 0.61 0.84 1.32 1.38 1.18

2.06 1.99 1.93 2.33 2.05 1.72 2.11 2.41 2.47

2.12 2.14 2.12 3.31 2.17 2.19 3.27 3.56 3.70

example. These hint at a weak correlation in the sense that for the first transition the amino group plays an important role while for the second transition, both, the hole and particle orbitals are somewhat pushed away to the lower right. But as discussed in the context of Eq. (26), it should be remembered that also cross terms, which cannot be understood from the visual representation, contribute to the mixed multipole moments, and thus Reh. Table 1 also serves as a good example to illustrate the numerical relations between the different 1TDM descriptors. In agreement with Eq. (22), it is seen that the distance between the charge centroids dh!e is always smaller than the dynamic exciton size dexc (by about 2 A˚) and a more precise consideration shows that they are connected by Eq. (20). Finally, the X values [Eq. (6)], which measure the single-excitation character,[13,29] are analyzed. These are all close to 0.85 with little variation indicating predominant single-excitation character of the states. For comparison, also an analysis of the difference density matrix (1DDM) is performed (Table 2). First, the promotion number [p, Eq. (37)] is computed. This quantity counts all effects that go beyond the simple picture of a one-electron transition, which could either derive from multiple excitation character or from orbital-relaxation effects.[28] As discussed previously,[28,79] the origin of these terms lies in the correlated nature of the ADC method and the underlying intermediate state representation.[52] Note that the construction of relaxed

Table 3. Analysis of the first five C 1s and O 1s core-excited states of cytosine at the CVS-ADC(2)-x/6-31111G** level of theory: excitation energies, oscillator strengths, various descriptors computed from the transition density matrix (see text), and state character. DE (eV)

1618

f

X

284.98 286.57 286.86 287.33 287.70

0.017 0.050 0.005 0.014 0.053

0.71 0.73 0.75 0.71 0.75

531.04 532.31 532.63 533.59 533.64

0.027 0.001 0.000 0.000 0.000

0.74 0.70 0.76 0.76 0.76

rh (A˚) C 1s 0.166 0.166 0.166 0.167 0.167 O 1s 0.122 0.122 0.122 0.122 0.122

re (A˚)

˚) dexc (A

Character

1.72 1.60 2.95 2.14 1.70

1.93 1.70 3.31 2.48 1.77

1s-p 1s-p 1s-Ryd. 1s-p 1s-p

1.75 2.00 3.10 4.27 4.17

2.22 3.66 3.96 5.41 5.77

1s-p 1s-p 1s-Ryd. 1s-Ryd. 1s-Ryd.

Journal of Computational Chemistry 2015, 36, 1609–1620

densities through the “Z-vector” method[80] would yield additional relaxation effects (see e.g., Ref. [81]), which are, however, not considered here. The p values vary from 1.37 to 1.79 in the present case. In light of the predominant single-excitation character, as indicated by the X values, this indicates the presence of relaxation effects. The dD!A values, i.e., the distances between the centers of the attachment and detachment densities, give a clue about the underlying physical processes. Since this value is always smaller than the corresponding dh!e value, it indicates a relaxation in the direction of reduced charge separation. By contrast, a comparison of the rD and rA values with the corresponding rh and re values only shows that the former sizes are somewhat blurred, i.e., the smaller values become larger and the larger values become smaller, which relates to the fact that orbital-relaxation effects come into play throughout the whole molecule. The final topic discussed here relates to states of entirely different character, the C 1s and O 1s core-excited states of cytosine (see Table 3). First the O 1s spectrum is discussed. Cytosine possesses only one oxygen atom and hence, only one very localized O 1s orbital from which the electron can be excited into various p and Rydberg orbitals. Compared to the valence-excited states discussed above (Table 1), the major difference is the hole size (rh), which has a value of exactly 0.122 ˚ for the five O 1s states shown here. The small size and the A little variations can be understood by the fact that the excitation starts from the very compact 1s orbital on the oxygen atom. On the other hand, the electron sizes vary and possess similar values as the ones of the valence-excited states showing that the virtual orbitals are more or less unaffected by the core excitation. Similar observations are made for the C 1s excited states of cytosine. In this case, there are four carbon atoms but due to the fact that these are not chemically identical the 1s orbitals are each localized on only one of them and therefore all coreexited states discussed here are dominated by only one carbon 1s orbital leading to very localized excitations. This is ˚ . Note that reflected by the small hole sizes of about 0.166 A this value is slightly larger when compared to the O 1s states, which follows from the decreased nuclear charge. All investigated core-excited states of cytosine exhibit X values between 70% and 76% indicating a dominant singleexcitation character. However, the double-excitation character is stronger than for the valence-excited states. This is in accordance with previous investigations that demonstrated typical amounts of doubly excited amplitudes between 20% and 30% for core-excited states.[58–60] The Reh and COV values are not presented in Table 3 as these are always 0.000. This total lack of correlation can be understood by the fact that the excitation proceeds from a very compact 1s orbital. It is also reflected by the fact that all PRNTO values are 1.000.

Conclusions In this work, a new formalism for the analysis of excited states is introduced and its implementation within the Q-Chem program package is reported. The method is based on the WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

physical picture of an exciton wavefunction, which is constructed from the 1TDM.[12] A subsequent statistical analysis of this two-body function in terms of spatial multipole moments provides a simple and well-defined route to a number of properties with clear physical meaning. It is not only possible to quantify the location and compactness of the independent distributions of the electron and hole quasi-particles, but also correlation effects can be quantified, which are at the heart of such phenomena as CR and excitonic effects. Due to the general definition of our formalism, the analysis is available for any quantum-chemical method that produces the required one-particle density and multipole matrices. Here, we use the ADC(2) method to compute valence and Rydberg-excited states, as well as the CVS-ADC(2)-x method for core-excited states; however, the formalism can be readily extended to ADC(3) and multireference configuration interaction as applied previously.[28,79] Furthermore, an extension of this toolbox to time-dependent density functional theory and coupled-cluster methods is currently in progress. The formalism was explained in some detail in terms of the underlying physical picture and a practical route for implementing the equations was given. The evaluation of various statistical moments and their physical meaning was discussed. The equations were recast into the picture of NTOs to allow for a different viewpoint on the relevant terms. While most of the work is focused on the analysis of the 1TDM, we have also shown how to apply the formalism in the 1DDM case. The analysis methods were first applied to the bianthryl molecule, which is known for its spectroscopic signature of CR states.[37,38] It was first shown how the CR character of the different states is modulated by the interring torsion. A new perspective on the CR phenomenon was then given by identifying it with an anticorrelation between the electron and hole quasiparticles. Moreover, this correlation effect was explained through the multiconfigurational nature of the state as measured by the number of nonvanishing NTOs eigenvalues (PRNTO ). Finally, the effect of the modulations in wavefunction character on the excited-state absorption strengths was studied and the data were compared to transient experiments.[38] As a second example, cytosine was chosen as a prototypical heteroaromatic molecule possessing a variety of excited states of different character. First, a connection between the computed descriptors and the qualitative state character was drawn. A clear separation between the states was found in the sense that (i) np states were characterized by small hole sizes rh, (ii) Rydberg states by large exciton sizes dexc and electron sizes re, and (iii) pp states by intermediate values of these. This finding suggests that our analysis tools can be used for an automatic assignment of excited-state characters without the need of visually inspecting orbitals. Second, the numerical relations between the different descriptors were illustrated. Third, an analysis of the 1DDM was performed to estimate the effects of orbital relaxation. Finally, core-excited states were analyzed, which possess the defining feature of very small hole sizes rh. This work illustrates the potential of our analysis methods in diverse situations, and a number of cases were identified where this approach may be particularly advantageous. The methods

open a route for an automatic assignment of excited-state characters without the need of visualizing orbitals. In addition, the new descriptors provide a way to quantitatively compare wavefunctions produced by different computational protocols and to evaluate possible shortcomings of various methods. More importantly, our analysis strategies provide a basis for understanding correlation effects, which underlie properties such as excitonic and CR effects, which are notoriously difficult to understand in the molecular-orbital picture. The methods discussed represent only the tip of the exciton analysis iceberg and there are a number of possible extensions to this approach, including the computation of higher multipole moments and covariance matrices, the inclusion of different operators such as the electron repulsion or kinetic energy, and the analysis of state DMs rather than 1TDMs. These can be readily integrated along the lines described here in case they are required for future applications.

Acknowledgments We dedicate this paper to Michael Wormit in grateful memory. He was a great colleague and an excellent scientist, this project would not have been possible without him. Keywords: excited states  wavefunction analysis  excitons  correlation

How to cite this article: F. Plasser, B. Thomitzni, S. A. B€appler, J. Wenzel, D. R Rehn, M. Wormit, A. Dreuw. J. Comput. Chem. 2015, 36, 1609–1620. DOI: 10.1002/jcc.23975

[1] A. Dreuw and M. Head-Gordon, Chem. Rev. 2005, 105, 4009. [2] F. Plasser, M. Barbatti, A. J. A. Aquino, and H. Lischka, Theor. Chim. Acta 2012a, 131, 1073. [3] L. Gonzalez, D. Escudero, and L. Serrano-Andres, ChemPhysChem 2012, 13, 28. [4] P.-O. Lowdin, Phys. Rev. 1955, 97, 1474. [5] A. V. Luzanov, A. A. Sukhorukov, and V. E. Umanskii, Theor. Exp. Chem. 1976, 10, 354. [6] R. L. Martin, J. Chem. Phys. 2003, 118, 4775. [7] I. Mayer, Chem. Phys. Lett. 2007a, 437, 284. [8] M. Head-Gordon, A. Grana, D. Maurice, and C. White, J. Phys. Chem. 1995, 99, 14261. [9] A. L. L. East and E. C. Lim, J. Chem. Phys. 2000, 113, 8981. [10] C. R. Kozak, K. A. Kistler, Z. Lu, and S. Matsika, J. Phys. Chem. B 2010, 114, 1674. [11] F. Plasser and H. Lischka, J. Chem. Theory Comput. 2012, 8, 2777. [12] S. A. B€appler, F. Plasser, M. Wormit, and A. Dreuw, Phys. Rev. A 2014, 90, 052521. [13] F. Plasser, M. Wormit, and A. Dreuw, J. Chem. Phys. 2014a, 141, 024106. [14] S. Tretiak and S. Mukamel, Chem. Rev. 2002, 102, 3171. [15] M. Sun, P. Kjellberg, W. J. Beenken, and T. Pullerits, Chem. Phys. 2006, 327, 474. [16] A. V. Luzanov and O. A. Zhikol, Int. J. Quantum Chem. 2010, 110, 902. [17] B. M. Wong and T. H. Hsieh, J. Chem. Theory Comput. 2010, 6, 3704. [18] F. Plasser, A. J. A. Aquino, W. L. Hase, and H. Lischka, J. Phys. Chem. A 2012b, 116, 11151. [19] S. Sharifzadeh, P. Darancet, L. Kronik, and J. B. Neaton, J. Phys. Chem. Lett. 2013, 4, 2197. [20] A. N. Panda, F. Plasser, A. J. A. Aquino, I. Burghardt, and H. Lischka, J. Phys. Chem. A 2013, 117, 2181. [21] L. Blancafort and A. A. Voityuk, J. Chem. Phys. 2014, 140, 095102.

Journal of Computational Chemistry 2015, 36, 1609–1620

1619

FULL PAPER [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64]

1620

WWW.C-CHEM.ORG

J. Huang, L. Du, D. Hu, and Z. Lan, J. Comp. Chem. 2015, 36, 151. K. Takatsuka, T. Fueno, and K. Yamaguchi, Theor. Chim. Acta 1978, 48, 175. M. Head-Gordon, Chem. Phys. Lett. 2003, 372, 508. F. Plasser, H. Pasalic, M. H. Gerzabek, F. Libisch, R. Reiter, J. Burgd€ orfer, T. M€ uller, R. Shepard, and H. Lischka, Angew. Chem. Int. Ed. 2013, 52, 2581. A. V. Luzanov, J. Struct. Chem. 2014, 55, 799. P. R. Surjan, Chem. Phys. Lett. 2007, 439, 393. F. Plasser, S. B€appler, M. Wormit, and A. Dreuw, J. Chem. Phys. 2014b, 141, 024107. S. Matsika, X. Feng, A. V. Luzanov, and A. I. Krylov, J. Phys. Chem. A 2014, 118, 11943. R. Richard and J. Herbert, J. Chem. Theor. Comput. 2011, 7, 1296. M. J. G. Peach, P. Benfield, T. Helgaker, and D. J. Tozer, J. Chem. Phys. 2008, 128, 44118. C. A. Guido, P. Cortona, and C. Adamo, J. Chem. Phys. 2014, 140, 104101. T. Etienne, X. Assfeld, and A. Monari, J. Chem. Theory Comput. 2014, 10, 3896. E. Ronca, M. Pastore, L. Belpassi, F. De Angelis, C. Angeli, R. Cimiraglia, and F. Tarantelli, J. Chem. Phys. 2014a, 140, 054110. H. Ma, T. Qin, and A. Troisi, J. Chem. Theory Comput. 2014, 10, 1272. T. Etienne, J. Chem. Theory Comput. 2015, 11, 1692. T. Takaya, S. Saha, H. Hamaguchi, M. Sarkar, A. Samanta, and K. Iwata, J. Phys. Chem. A 2006, 110, 4291. N. Asami, T. Takaya, S. Yabumoto, S. Shigeto, H.-o. Hamaguchi, and K. Iwata, J. Phys. Chem. A 2010, 114, 6351. M. Rohlfing and S. G. Louie, Phys. Rev. B 2000, 62, 4927. W. Hieringer and A. G€ orling, Chem. Phys. Lett. 2006, 419, 557. P. Petelenz and B. Pac, J. Am. Chem. Soc. 2013, 135, 17379. S. Brazovskii and N. Kirova, Chem. Soc. Rev. 2010, 39, 2453. M. Knupfer, J. Fink, E. Zojer, G. Leising, and D. Fichou, Chem. Phys. Lett. 2000, 318, 585. D. Kannar and P. G. Szalay, J. Mol. Model. 2014, 20, 2503. F. Weigend and M. Haser, Theor. Chem. Acc. 1997, 97, 331. A. Schafer, H. Horn, and R. Ahlrichs, J. Chem. Phys. 1992, 97, 2571. J. Schirmer, Phys. Rev. A 1982, 26, 2395. A. B. Trofimov and J. Schirmer, J. Phys. B 1995, 28, 2299. E. Epifanovsky, D. Zuev, X. Feng, K. Khistyaev, Y. Shao, and A. I. Krylov, J. Chem. Phys. 2013, 139, 134105. M. Wormit, D. R. Rehn, P. H. Harbach, J. Wenzel, C. M. Krauter, E. Epifanovsky, and A. Dreuw, Mol. Phys. 2014, 112, 774. A. Dreuw and M. Wormit, WIREs Comput. Mol. Sci. 2015, 5, 82. J. Schirmer and A. B. Trofimov, J. Chem. Phys. 2004, 120, 11449. S. Knippenberg, D. R. Rehn, M. Wormit, J. H. Starcke, I. L. Rusakova, A. B. Trofimov, and A. Dreuw, J. Chem. Phys. 2012, 136, 064107. R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 1992, 96, 6796. L. S. Cederbaum, W. Domcke, and J. Schirmer, Phys. Rev. A 1980, 22, 206. A. Barth and L. S. Cederbaum, Phys. Rev. A 1981, 23, 1038. A. Barth and J. Schirmer, J. Phys. B: Atom. Mol. Phys. 1985, 18, 867. J. Wenzel, M. Wormit, and A. Dreuw, J. Comp. Chem. 2014a, 35, 1900. J. Wenzel, M. Wormit, and A. Dreuw, J. Chem. Theory Comput. 2014b, 10, 4583. J. Wenzel, A. Holzer, M. Wormit, and A. Dreuw, J. Chem. Phys. 2015, 142, 214104. R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, J. Chem. Phys. 1980, 72, 650. T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. V. R. Schleyer, J. Comp. Chem. 1983, 4, 294. R. Ahlrichs, M. Bar, M. Haser, H. Horn, and C. Kolmel, Chem. Phys. Lett. 1989, 162, 165. A. I. Krylov and P. M. W. Gill, WIREs: Comp. Mol. Sci. 2013, 3, 317.

Journal of Computational Chemistry 2015, 36, 1609–1620

[65] Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kus, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. DiStasio, H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. D. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, E. Neuscamman, C. M. Oana, R. Olivares-Amaya, D. P. O’Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. St€ uck, Y.-C. Su, A. J. W. Thom, T. Tsuchimochi, V. Vanovschi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, J. Yang, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K. L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Z. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill, M. HeadGordon, Mol. Phys. 2015, 113, 184. [66] F. Plasser, TheoDORE: A Package for Theoretical Density, Orbital Relaxation, and Exciton Analysis, Available at http://theodore-qc. sourceforge.net/. [67] K. Yamasaki, K. Arita, O. Kajimoto, and K. Hara, Chem. Phys. Lett. 1986, 123, 277. [68] N. Mataga, H. Yao, T. Okada, and W. Rettig, J. Phys. Chem. 1989, 93, 3383. [69] S. A. Kovalenko, J. L. P erez Lustres, N. P. Ernsting, and W. Rettig, J. Phys. Chem. A 2003, 107, 10228. [70] T. J. Kang, W. Jarzeba, P. F. Barbara, and T. Fonseca, Chem. Phys. 1990, 149, 81. [71] F. C. Grozema, M. Swart, R. W. J. Zijlstra, J. J. Piet, L. D. A. Siebbeles, and P. T. van Duijnen, J. Am. Chem. Soc. 2005, 127, 11019. [72] G. -J. Zhao, Y. -H. Liu, K. -L. Han, and Y. Dou, Chem. Phys. Lett. 2008, 453, 29. [73] A. Dreuw and M. Head-Gordon, Chem. Phys. Lett. 2006, 426, 231. [74] I. Mayer, Chem. Phys. Lett. 2007b, 443, 420. [75] N. Ismail, L. Blancafort, M. Olivucci, B. Kohler, and M. A. Robb, J. Am. Chem. Soc. 2002, 124, 6818. [76] M. Merchan and L. Serrano-Andres, J. Am. Chem. Soc. 2003, 125, 8108. [77] H. R. Hudock and T. J. Martinez, ChemPhysChem 2008, 9, 2486. [78] M. Barbatti, A. J. A. Aquino, J. J. Szymczak, D. Nachtigallova, P. Hobza, and H. Lischka, Proc. Natl. Acad. Sci., USA 2010, 107, 21453. [79] F. Plasser and A. Dreuw, J. Phys. Chem. A 2015, 119, 1023. [80] N. C. Handy and H. F. Schaefer, J. Chem. Phys. 1984, 81, 5031. [81] E. Ronca, C. Angeli, L. Belpassi, F. D. Angelis, F. Tarantelli, and M. Pastore, J. Chem. Theory Comput. 2014b, 10, 4014.

Received: 22 April 2015 Revised: 26 May 2015 Accepted: 28 May 2015 Published online on 29 June 2015

WWW.CHEMISTRYVIEWS.COM