Excited states with internally contracted multireference coupled-cluster linear response theory.

Excited states with internally contracted multireference coupled-cluster linear response theory Pradipta Kumar Samanta, Debashis Mukherjee, Matthias Hanauer, and Andreas Köhn Citation: The Journal of Chemical Physics 140, 134108 (2014); doi: 10.1063/1.4869719 View online: http://dx.doi.org/10.1063/1.4869719 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/13?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Additional global internal contraction in variations of multireference equation of motion coupled cluster theory J. Chem. Phys. 138, 134108 (2013); 10.1063/1.4796523 Linear-response theory for Mukherjee's multireference coupled-cluster method: Excitation energies J. Chem. Phys. 137, 044116 (2012); 10.1063/1.4734309 Dyson orbitals for ionization from the ground and electronically excited states within equation-of-motion coupledcluster formalism: Theory, implementation, and examples J. Chem. Phys. 127, 234106 (2007); 10.1063/1.2805393 Coupled-cluster response theory with linear- r 12 corrections: The CC2-R12 model for excitation energies J. Chem. Phys. 124, 044112 (2006); 10.1063/1.2161183 Calculation of excited-state properties using general coupled-cluster and configuration-interaction models J. Chem. Phys. 121, 9257 (2004); 10.1063/1.1805494

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 205.155.154.122 On: Sat, 05 Apr 2014 11:12:08

THE JOURNAL OF CHEMICAL PHYSICS 140, 134108 (2014)

Excited states with internally contracted multireference coupled-cluster linear response theory Pradipta Kumar Samanta,1,a) Debashis Mukherjee,1,b) Matthias Hanauer,2,c) and Andreas Köhn2,d) 1

Raman Center for Atomic, Molecular and Optical Sciences, Indian Association for the Cultivation of Science, Kolkata 700 032, India 2 Institut für Physikalische Chemie, Universität Mainz, D-55099 Mainz, Germany

(Received 19 November 2013; accepted 17 March 2014; published online 4 April 2014) In this paper, the linear response (LR) theory for the variant of internally contracted multireference coupled cluster (ic-MRCC) theory described by Hanauer and Köhn [J. Chem. Phys. 134, 204211 (2011)] has been formulated and implemented for the computation of the excitation energies relative to a ground state of pronounced multireference character. We find that straightforward application of the linear-response formalism to the time-averaged ic-MRCC Lagrangian leads to unphysical second-order poles. However, the coupling matrix elements that cause this behavior are shown to be negligible whenever the internally contracted approximation as such is justified. Hence, for the numerical implementation of the method, we adopt a Tamm-Dancoff-type approximation and neglect these couplings. This approximation is also consistent with an equation-of-motion based derivation, which neglects these couplings right from the start. We have implemented the linear-response approach in the ic-MRCC singles-and-doubles framework and applied our method to calculate excitation energies for a number of molecules ranging from CH2 to p-benzyne and conjugated polyenes (up to octatetraene). The computed excitation energies are found to be very accurate, even for the notoriously difficult case of doubly excited states. The ic-MRCC-LR theory is also applicable to systems with open-shell ground-state wavefunctions and is by construction not biased towards a particular reference determinant. We have also compared the linear-response approach to the computation of energy differences by direct state-specific ic-MRCC calculations. We finally compare to Mk-MRCC-LR theory for which spurious roots have been reported [T.-C. Jagau and J. Gauss, J. Chem. Phys. 137, 044116 (2012)], being due to the use of sufficiency conditions to solve the MkMRCC equations. No such problem is present in ic-MRCC-LR theory. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4869719] I. INTRODUCTION

The accurate computation of excited states remains a challenge in electronic structure theory. In many cases, pronounced multireference character is observed even for the ground state, and it is not always possible to express the dominant electronic configurations in the excited state as single excitations with respect to one prominent ground-state configuration. The existing approaches to excitation energies in such situations roughly fall into two categories: (1) those that provide absolute excited state energies of all the considered states, and (2) those that directly target the difference between the excited states and the ground state and thus directly yield excitation energies. Some representatives of the first group are multireference configuration interaction (MRCI), multireference perturbation theory (CASPT2), and related methods, see Refs. 1 and 2 for reviews. They can properly deal with the nondynamical correlations stemming from the multireference character encountered in various states, but the treatment of dynamic eleca) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected]

0021-9606/2014/140(13)/134108/14/$30.00

tron correlation is often not comparable to the high standard that can be reached in single-reference coupled-cluster (CC) theory for ground states.3, 4 The corresponding multireference treatment of excited states at the CC level may be achieved by traditional spinorbital-based SUMRCC5–8 as also by the recently developed spin-free UGA-SUMRCC,9 which is related to pioneering work of Li and Paldus.10 Both of these methods, however, suffer from the problem of intruders.8, 11, 12 The methods belonging to the second group take advantage of the situation that certain excited configurations, such as the core-to-virtual double excitations, have very similar importance in the ground and excited states. They can thus be kept the same while describing the two states. An attractive possibility is to express the excited state with respect to the ground state, so that common correlation terms drop out in the energy differences. This can be realized in two different ways. One may, e.g., generate a linear response (LR) of the ground state to a harmonic time-dependent perturbation and seek the poles of the Fourier transformed linear response function. The other approach would be via an ansatz for the excited state which contains the information of the wave operator of the ground state. Here, one posits a factorized ansatz ˆ ex ˆ of the excited states such as ˆ g for the wave operator ˆ ˆ ˆ or ex g where g is the wave operator for the ground state

140, 134108-1

© 2014 AIP Publishing LLC


134108-2

Samanta et al.

ˆ ex is the wave operator bringing in valence correlations, and differential correlations, etc., for the excited states. Such a strategy has been followed in the valence-universal MRCC scheme13–16 which is based on Fock space and is best suited for spectroscopic energy differences like excitation energies, ionization potentials, and electron affinities. Recently, a factorized ansatz based on the aforementioned UGA-SUMRCC has been employed which gives energy differences directly but without the calculation of all intermediate valence sectors between the states of interest unlike VUMRCC. This has been called UGA-QFMRCC.17 In principle, one can also use ˆ ex , leading to a multireference a linear parameterization for analogue of the so-called equation-of-motion (EOM) method (vide infra). Somewhat intermediate in structure with respect to the group 1 and 2 in a CI context is the difference dedicated CI (DD-CI)18, 19 method, which leaves out the core-to-virtual double excitations in both the ground and excited states. This offers both a more efficient and more balanced description of the excitation energy, but it is not as accurate as the corresponding VUMRCC or UGA-QFMRCC methods. It should be noted here that in the context of singlereference CC theory, the direct methods can be derived by either linear response theory20–24 (CC-LR) or equation-ofmotion theory25–27 (EOM-CC). A very similar formalism called symmetry adapted cluster configuration interaction (SAC-CI) was proposed by Nakatsuji.28 The EOM methodology has also been generalized in order to describe multiconfigurational states via multi-ionization29 or spin-flip,30, 31 starting from a related single-configurational state. As far as excitation energies are concerned, CC-LR and EOMCC methods are perfectly equivalent at the single-reference level.32 However, response theory is a more general framework that offers a route to a multitude of other properties such as nonlinear optical and magnetic properties.24 The singles-doubles approximation of EOM-CC (EOMCCSD), however, cannot provide an accurate description of excited states dominated by doubles and in situations where the ground state itself has quasi-degenerate character. An extension of EOM-CCSD, inspired by the active-space concept, has been introduced by Kowalski and Piecuch33, 34 and is known as EOM-CCSDt. This method has its roots in Adamowicz’s35 state-specific MRCC approach and uses cluster operators and excitation operators of higher than twobody rank, defined through active orbitals. These higher-body operators are of internal and semi-internal type. A different method known as left-eigenstate completely renormalized equation-of-motion coupled cluster approach with singles, doubles, and non-iterative triples [CR-EOMCCSD(2,3)] has been developed by Piecuch and co-workers36 later. In this method, an energy term that has the effect of the triply excited configurations via the triexcited moment of the EOMCCSD is added as a correction to the EOM-CCSD energy. Both these methods and the size-intensive modification δ-CREOMCC(2,3)37 have been used to study excited states dominated by double excitations and potential energy surfaces of excited states.36, 37 These methods, however, are based on a single-reference framework and they cannot be expected to be equally reliable for systems with varying degree of multireference character.

J. Chem. Phys. 140, 134108 (2014)

A promising strategy for extending the power of the linear response approach to systems with multireference character is to apply it to a highly reliable state-specific multireference coupled-cluster38–40 (MRCC) method. Most notably, this has been achieved for the well-known state-specific MRCC method of Mukherjee et al.41–46 (SSMRCC or Mk-MRCC), in the work by Chattopadhyay et al.,47 and the recent work of Jagau and Gauss.48 The latter have also extended the theory to the dynamic linear response function.49 However, these studies47–49 have revealed a problem in the response approach applied to the Mk-MRCC method, when states are targeted which lie outside the model space (i.e., the zeroth-order space of the underlying ground-state calculation). The excitation manifold, and hence the response vector, contains redundancies at linear order, viz., an individual excited determinant may be reached by several excitations starting from different reference determinants (also known as the multiple parentage problem).47 As a consequence, the theory predicts additional unphysical excited states. Mention should also be made in this context of the multireference equation-of-motion (MREOM) methods of Nooijen and co-workers.50–53 These methods employ a two-step procedure, in which the Hamiltonian is first similarity transformed within an internally contracted approach, and then diagonalized in an uncontracted CI space, thereby yielding the energies of ground and several excited states at once. While in principle, a strictly state-specific MRCC method may allow one to also directly target an excited state,54 this strategy would be limited to states that can be described by the model space and could not be readily extended to dynamical response properties. In the present work, we discuss a linear response theory for excited states within what can be termed as an internally contracted MRCC (ic-MRCC) method. An internally contracted MRCC method makes use of a single wave operator which acts on an unperturbed function 0 , which itself is a linear combination of model space functions. Banerjee ˆ = eTˆ and Simons55, 56 had considered an exponential ansatz of the wave operator where Tˆ was confined only to excitations from core and active orbitals and omitted excitations into the active orbitals. This makes the components of Tˆ commuting, which led to the MRCC equations terminating at the quartic power of the cluster amplitudes. Consequently, this theory missed the internal and the semi-internal correlations. Another general method57, 58 using a suitably normal ˆ was suggested by Mukherordered exponential ansatz59 for jee, in which all possible components of Tˆ were included. The method of Mukherjee60 has unfortunately not yet been applied hitherto, but is currently studied.60, 61 On the other hand, the general formulation using an ordinary exponential parameterization has been developed by Evangelista and Gauss62 and also by Hanauer and Köhn63, 64 which led to very promising results. All these variants of icMRCC57, 58, 62, 63 use only a set of linearly independent excitations in Tˆ . It is hence not expected that a linear response of the formulation will lead to spurious excitation energies. The focus of the present work lies on the development and implementation of a linear response theory for the icMRCCSD variant of Hanauer and Köhn.63 We note that we


134108-3

Samanta et al.

J. Chem. Phys. 140, 134108 (2014)

are not exploiting in this paper some of the allied developments within ic-MRCC theory, such as the use of a sequential ansatz in order to include higher commutators with respect to single excitations,65 or the rigorously size-extensive formalism based on a generalized normal ordering.57, 60, 61 Efforts to develop a best practice strategy within ic-MRCC theory are underway and will be the subject of separate publications. During the present study, it has turned out that a Lagrangian-based derivation following a general response framework leads to slightly different equations than an EOMlike derivation of the equations for excitation energies. This study will mainly focus on the former approach and leave the latter one for future studies. Section II of this article falls into a summary of the icMRCC method and the derivation of the equations for excitation energies through the response formalism, as well as a discussion of the pole structure of the dynamic response function. The variant that follows from an EOM-like strategy is also briefly sketched. After providing details on the implementation of the method, we present a number of applications. The CH2 molecule serves as a test system for studying the method’s accuracy, in particular for doubly excited states. Singlet states of para-benzyne are used for a comparison with Mk-MRCC linear response theory, and to re-evaluate the severity of its problem with spurious poles. For selected states of CH2 and p-benzyne, we then compare the linear response approach of ic-MRCC theory with a direct state-specific treatment. We finally apply our ic-MRCC linear response theory to the challenging case of excited states in polyenes.

|H¯ 0 |c = Ec ,

(4)

0 |τˆ † H¯ 0 |0 = 0 .

(5)

ˆ ˆ H¯ 0 = e−T Hˆ 0 eT is the similarity transformed electronic Hamiltonian of the unperturbed system, pq rs Hˆ 0 = 0|Hˆ 0 |0 + aˆ pq , fqp aˆ pq + 14 grs (6) pq

pqrs

which is normal ordered with respect to the core determinant |0 that comprises all doubly occupied orbitals in | 0 . Equation (4) is an eigenvalue equation and reads in a more explicit notation like μ |H¯ 0 |ν cν = Ecμ . (7) ν

The prime in τˆ in Eq. (5) indicates that the projection is restricted to a linearly independent excitation manifold {τˆρ |0 }. It is defined by a linear transformation of the excitation operators, τˆ = τˆ X,

II. THEORY

In the following, we briefly recapitulate the ic-MRCC variant of Hanauer and Köhn, as described in Refs. 63 and 64. The ic-MRCCSD ansatz for the wavefunction of a specific (initial) electronic state is Tˆ

| = e |0 ,

(1)

where the reference function |0 = |μ cμ = |c

(2)

μ

is a linear combination of determinants from a complete active space (CAS), and the cluster operator I J AB Tˆ = aˆ I J = tAI aˆ IA + 14 tAB τˆρ tρ = τˆ t (3) I J AB

(8)

where the transformation matrix X is chosen such that the excited functions become normalized and mutually orthogonal,

A. The ic-MRCCSD method

IA

determinants and excitation operators, respectively, and c and t are column vectors of the corresponding wavefunction parameters. The equations for the latter follow from insertion of Eq. (1) into the time-independent Schrödinger equation, mulˆ tiplication from the left with e−T , and projection onto both the set of reference determinants and the excited states that can be reached from the reference function by the excitation operators:

0 |τˆρ† τˆσ |0 = δρσ .

(9)

For an investigation of different choices of X, we refer to the various procedures outlined in Refs. 58, 60, 61, 63, 66, and 67. The default choice of the present work is the sequential orthogonalization based on excitation operators normal ordered with respect to a closed-shell core determinant.63 In this scheme, a canonical orthogonalization of the singly excited functions is succeeded by a projection of the doubly excited functions onto a space orthogonal to the singly excited ones, and a canonical orthogonalization of the resulting doubly excited functions.

ρ

performs single and double excitations from the combined set of (correlated) doubly occupied spinorbitals and active spinorbitals (labeled by I, J) to the combined set of active and virtual spinorbitals (A, B). From now on, we will refer to the spinorbitals simply as orbitals. Since we treat the coefficients cμ in | 0 as additional wavefunction parameters, we consider cluster amplitudes with indices labeling active orbitals only (u, v, . . .) as redundant and set them to zero, i.e., tvu = 0, uv = 0. Equations (2) and (3) introduce a shorthand notation, twx in which | and τˆ are row vectors comprising the reference

B. ic-MRCC excitation energies

In the present article, we follow the response approach as described by Christiansen et al.24 Within this theory, response functions and response equations are derived from a time-averaged quasi-energy Lagrangian, and excitation energies are obtained as poles of the linear response function. Jagau and Gauss48, 49 have recently applied this approach to response functions and excitation energies in Mk-MRCC theory, and derived equations for excitation energies identical to those of a previous work by Mukherjee and co-workers47


Samanta et al.

134108-4

J. Chem. Phys. 140, 134108 (2014)

who used a slightly different formalism. We note that, for icMRCC theories alternative approaches to excitation energies do not necessarily result in the same working equations. A more detailed discussion will be reserved for future publications. It is tempting to skip the derivation of the time-averaged quasi-energy ic-MRCC Lagrangian and directly postulate the excitation energies as eigenvalues of the ic-MRCC Jacobian. However, we will present the derivation in detail, as the resulting eigensystem turns out to have a nontrivial right-hand side. We start from the time-dependent Schrödinger equation ∂ ˆ | = 0, (10) H −i ∂t

ˆ ˆ H¯ = e−T Hˆ eT , we get

where Hˆ = Hˆ (t) = Hˆ 0 + Vˆ (t) and | = |(t) depend on time. Here,

while the amplitude equations follow from projection onto ˆ 0 |τˆ † e−T eiQt

ˆ −iωt + εX (−ω)Xe ˆ iωt Vˆ (t) = εX (ω)Xe

ˆ

ˆ

e

−iQt

|H¯ |0 − Q|0 − i

−e−iQt

ˆ ∂eT Tˆ ∂0 |0 + ie Qe |0 + i = 0, ∂t ∂t Tˆ

(13) where we have omitted explicit time-dependence of Hˆ , Tˆ , and | 0 for brevity. An expression for the quasi-energy follows from projecˆ tion onto 0 |e−T eiQt , Tˆ

ˆ ˆ ∂e ˆ |0 0 |e−T Hˆ eT |0 − Q − i0 |e−T ∂t ∂0 −i 0 = 0. ∂t

(15)

(16)

Tˆ

(17)

In order to reformulate the second term in Eq. (17), let us first consider e−T

ˆ

ˆ

∂eT = 1 − Tˆ + 12 Tˆ Tˆ − . . . ∂t

× T˙ + 12 T˙ Tˆ + 12 Tˆ T˙ + 16 T˙ Tˆ Tˆ + 16 Tˆ T˙ Tˆ + 16 Tˆ Tˆ T˙ + . . .

= T˙ + 12 [T˙ , Tˆ ] + 16 [[T˙ , Tˆ ], Tˆ ] + · · · . (18) Upon defining the matrix

St = 0 |τˆ † τˆ + 12 [τˆ , Tˆ ] + 16 [[τˆ , Tˆ ], Tˆ ] + . . . |0 , (19) Eq. (17) takes the form 0 |τˆ † H¯ |0 − iSt

Tˆ

Hˆ e |0

∂c = 0, ∂t

ˆ ∂e |0 = 0. 0 |τˆ † H¯ |0 − i0 |τˆ † e−T ∂t

(12)

where eT (t) | 0 (t) is a Floquet state periodic in ω and Q is the quasi-energy, which is time-independent for periodic perturbations.68 The cluster amplitudes in Tˆ and the parameters cμ in the reference function are periodic in time with . However, we follow an orbital-unrelaxed apperiodicity 2π ω proach, in which the orbitals in | 0 remain unperturbed. Insertion of (12) into (10) yields

∂c . ∂t

Ensuring that Eq. (15) remains time-independent requires additional conditions, which can always be fulfilled for periodic perturbations. It is more convenient, however, not to consider these conditions explicitly, but to work with time-averaged quantities,68 as we will do below, see Eq. (24). We obtain the equation for the reference function by projecting Eq. (13) onto ˆ |e−T eiQt ,

(11)

is a time-dependent perturbation. For simplicity, we restrict the discussion to a single periodic perturbation, a generalization to multiple perturbations is straightforward.24 The field strength parameters εX (ω) and εX (−ω) are associated with a time-independent perturbation operator Xˆ and a frequency ω or −ω, respectively. The Hermiticity of Vˆ (t) is ensured by the relationship εX (ω)∗ = εX (−ω). The time-dependent icMRCC ansatz for the time-dependent wavefunction, (t), can be written in the Floquet representation as |(t) = e−iQt eT (t) |0 (t),

Q = 0 |H¯ |0 − ic†

∂t = 0. ∂t

(20)

Here, St is a non-unit matrix. The multi-commutator terms involving τˆ and Tˆ in the expression of St are nonzero because of the non-commutativity of the excitations in the cluster operator. For later reference, we subsume this whole multi-commutator expression as a modified excitation operator τ˜ . Let us now assemble the time-dependent ic-MRCCSD Lagrangian L by adding to the quasi-energy Q the conditions (16) and (20), equipped with time-dependent Lagrange multipliers, to be absorbed into c¯μ μ | = c¯ | (21) 0 | = μ

(14)

The second term in Eq. (14) arises as | 0 is normalized to unity. The third term of the same equation, on the other hand, ˆ vanishes because ∂eT /∂t always creates at least one hole or particle in the inactive orbital spaces when acting on | 0 , and therefore yields a function orthogonal to | 0 . Upon defining

and ˆ =

λρ τˆρ† = λ τˆ † .

(22)

ρ

This gives a form similar to the time-independent ic-MRCC Lagrangian,63 however with an additional term containing the


134108-5

Samanta et al.

J. Chem. Phys. 140, 134108 (2014) TABLE I. Selected tensors appearing in ic-MRCC response theory.

time derivatives of the wavefunction parameters, ˆ H¯ |0 − Q(0 |0 − 1) L = 0 |H¯ |0 + 0 | ∂c 1 0

∂t − i c¯ λ . (23) ∂t 0 St ∂t The Lagrangian in Eq. (23) is the starting point for obtaining response functions and response equations. In essence, we have to expand all operators in a Fourier series and then form the time-averaged Lagrangian {L}T according to 1 t0 +T L(t)dt, (24) {L}T = T t0 which now depends on frequency and not on time. In Eq. (24), T is the period of the perturbation. From this expression, a linear response function can be defined as d 2 {L}T 1 , X; Y ω = Cˆ ±ω (25) 2 dεX (−ω)dεY (ω) ε=0 where the symmetrization operator Cˆ ±ω acts according to Cˆ ±ω f (ω1 , ω2 ) = f (ω1 , ω2 ) + f ∗ (−ω1 , −ω2 ) in order to ensure the correct symmetry of the response function, X; Y ω = X; Y ∗−ω .24 A more detailed study of the ic-MRCC linear response function will be published elsewhere. Here, we will briefly discuss its pole structure, as the poles of X; Yω can be identified as the ic-MRCC excitation energies. The response equations for the first-order parameters cX (ω) and tX (ω) follow from differentiation of the response function, Eq. (25), with respect to c¯ Y (−ω) and λY (−ω), respectively. It turns out that the equations for cX (ω) and tX (ω) are coupled to the equations for cX (−ω)∗ and tX (−ω)∗ due to the appearance of 0 | in the residual equations [Eq. (17)]: ⎞ ⎛ ⎞⎤ ⎡⎛ 1 0 0 0 Act¯ 0 0 Acc ¯ ⎜0 ⎥ ⎢⎜A 0 ⎟ 0 0 ⎟ λt ⎟ ⎜ ⎟⎥ ⎢⎜ λc Aλt Bλc − ω ⎟ ⎜ ⎟ ⎥ ⎢⎜ ⎝0 0 −1 ⎣⎝ 0 0 A∗cc 0 ⎠⎦ A∗ct¯ ⎠ ¯ B∗λc ⎛

0 cX (ω)

A∗λc ⎞

A∗λt ⎛

0 ξX c¯

⎞

⎜ t (ω) ⎟ ⎜ ξX ⎟ ⎜ ⎜ λ ⎟ ⎟ ×⎜ X ⎟ = −⎜ X∗ ⎟. ∗ ⎝ cX (−ω) ⎠ ⎝ ξ c¯ ⎠ tX (−ω)∗

0

0

− ∗λt

(26)

ξ X∗ λ

The definitions for the tensors from Eq. (26) are collected in Table I. This equation diverges if the frequency ω makes the matrix on the left-hand side singular, leading to the well-known eigenvalue problem from which the poles of the (linear) response function—and hence the excitation energies of the system—can be determined. This eigensystem’s structure is similar to the one in time-dependent Hartree-Fock (TDHF) theory,69 and one can show that pairs of solutions exist with corresponding eigenvalues ωi and −ωi , respectively. Since the excitation operators τ are chosen to be linearly independent, all elements of the response vector are non-redundant and it is evident that the ic-MRCC linear response function has the right number of poles, in contrast to the Mk-MRCC linear response theory.47, 49 However, the off-diagonal coupling

Abbreviation

Expression

(Acc ¯ )μν

μ |(H¯ 0 − E)|ν

(Act¯ )μρ

μ | ∂ H(0)0 |0

(Aλc )ρμ (Aλt )ρσ

¯

(0)

∂tρ (0) † 0 |τˆρ H¯ 0 |μ (0) † ¯ (0) 0 |τˆρ ∂ H(0)0 |0 ∂t σ

(Bλc )ρμ

† (0) μ |τˆρ H¯ 0 |0

( λt )ρσ

0 |τˆρ e−T

(ξ X c¯ )μ (ξ X λ )ρ

ˆ (0) ∂eTˆ (0) (0) (0) |0 ∂tσ (0) (0) ˆ ˆ Tˆ | (0) μ |e−T Xe 0 (0) † −Tˆ (0) ˆ Tˆ (0) (0) 0 |τˆρ e Xe |0 (0)

†

through Bλc causes singularities in cX (ω) and tX (ω) at both ω = ωi and ω = −ωi . Terms in the linear response function that are quadratic in these parameters may then give rise to unphysical second-order poles (cf. Refs. 24, 70, and 71). We note that the latter problem does not arise in Mk-MRCC linear response theory. The coupling term Bλc arises due to the occurrence of the (time-dependent) reference function in the projection manifold. It contains projections of H¯ 0 |0(0) onto excited functions μ |τˆρ† . Since projections onto functions from the projection manifold vanish due to fulfillment of the ic-MRCC equations for the reference state, the Bλc term will only contain contributions from an extended ground state residual arising from projections onto functions excited by operators of higher rank, † (0) (0) † (0) μ τˆρ τˆσ 0 0 τˆσ H¯ 0 0 . (Bλc )ρμ = σ rank(σ ) > rank(ρ)

(27) The first factor in Eq. (27) naturally vanishes in the full configuration interaction (FCI) limit, and it vanishes individually for each excitation class (number of inactive annihilation and creation operators in τˆρ ), for which the maximum operator rank is large enough to cover the complete space of functions within this excitation class. Like the internally contracted MRCI method,72 the ic-MRCC method relies on the observation that one can truncate the excitations in each class at a global rank, i.e., rank 2 in ic-MRCCSD, without deteriorating the accuracy of the method. The residual terms in Eq. (27) are thus expected to be small in general. We view the Bλc term as an artifact stemming from our way to account for the relaxation of the reference function through the ic-MRCC Lagrangian via the use of the 0 | which is the adjoint bra of | 0 . Though one can envision a replacement of the left-hand reference function by a biorthogonal complement which would have eliminated the second-order contributions to the poles, this would require a simultaneous solution of the equations for wavefunction parameters and Lagrange multipliers.63 Naturally, the Bλc term is also not present when one chooses to fix the reference function and instead uses internal excitations in the cluster operator. Moreover, one may think of a non-Lagrangian formulation as an alternative


Samanta et al.

134108-6

J. Chem. Phys. 140, 134108 (2014)

EOM-inspired derivation, which will not lead to this coupling term either (vide infra, Eq. (31) and the text following it). Therefore, in this work we choose to simplify the equations in the spirit of the Tamm-Dancoff approximation and neglect the term Bλc , i.e., we remove the coupling between solutions for positive and negative frequencies. The resulting eigenvalue equation for the ic-MRCC transition energies then has only half of the dimension of Eq. (26): 1 0 rc Act¯ Acc ¯ −ω = 0. (28) 0 λt Aλc Aλt rt It is worth mentioning that the terms Act¯ and Aλt can be rewritten by expanding the derivative of the similarity transformed Hamiltonian as −Tˆ

Tˆ (0)

∂ H¯ 0 ∂e ∂e ˆ (0) ˆ (0) = (0) Hˆ 0 eT + e−T Hˆ 0 (0) (0) ∂t ∂t ∂t

Tˆ ˆ (0) ∂e ¯ ˜ , = [H¯ 0 , e−T (0) ] = H0 , τ ∂t (0)

(29)

Eq. (8). The core step in this procedure is the diagonalization of overlap matrices of the excitation operators, and requires one to introduce a numerical threshold parameter η below which eigenvalues of the overlap matrices are discarded. Information about the choice of this threshold as well as on the basis sets and frozen core approximations used is given in Sec. IV. The expressions for the product of the Jacobian of Eq. (28) and the response vector are obtained by differentiation of the ic-MRCC Lagrangian [the time-independent part of Eq. (23)], whereas the product of the metric matrix and the response vector is defined separately, following Eq. (19). The eigenvalue equations for the excitation energies take the following form: A r − ω r = 0,

where both the matrices A and and the vector r are written in the orthogonalized excitation basis. Using Eq. (8), this can be reformulated as a transformed equation in the original excitation basis, X† (Ar − ωr) = 0,

where τ˜ is defined, as also mentioned earlier, according to τ˜ = e−T

ˆ (0)

= τˆ +

Tˆ (0)

∂e ∂ t (0)

(30)

1 ˆ (0) [τˆ , T ] 2

+

1 [[τˆ , Tˆ (0) ], Tˆ (0) ] 6

+ ··· .

The new definition of Rˆ = τ˜ rt = τ r˜ t

III. DETAILS OF IMPLEMENTATION

We have realized the proposed method, to be called icMRCC linear response theory (ic-MRCC-LR) in the following, as an extension of the current implementation of icMRCC theory in the General Contraction Code (G E CC O). We refer to Ref. 64 for a more detailed description. G E CC O reads in the Hamiltonian integrals from CASSCF calculations performed by the programs DALTON73 or G AMESS,74 and provides a set of tools for the automated derivation and evaluation of second quantized expressions as well as subspace methods for the solution of systems of equations. The computation of excited states is preceded by a ground state calculation. Unless otherwise noted, we truncate only the term ˆ operator in the Lagrangian [Eq. (23)] at the containing the twofold commutator level (quadratic power in the cluster operator) and use a normal ordering with respect to a closedshell core determinant and the sequential way to remove redundancies from Tˆ1 and Tˆ2 as described in Ref. 63. The latter affects the choice of the transformation matrix X from

(33)

where the expressions for A and are straightforwardly obtained from those of A and by substituting τˆ for τˆ . In each iteration, we compute the products Ar and r and then evaluate the left-hand side of Eq. (33). A new subspace vector is then generated according to bx[n+1] = −

(31)

gives us the freedom to follow another way to solve Eq. (28) by choosing r˜ t as the variable and truncating it in a similar way as the cluster operator. The resulting equation will then have an EOM-like structure containing the commutator ˆ The same Rˆ will appear in the ω containing term and [H¯ 0 , R]. turn λt into a unit matrix according to its definition in Table I. This alternative approach will be discussed in a more rigorous way in a forthcoming publication.

(32)

[X† (Ar[n] − ω[n] r[n] )]x . dx − ω[n]

As a preconditioner

d=

dc , dt

(34)

(35)

we use dt,ρ = 0 |τˆρ† [Hˆ D , τˆρ ]|0 ,

(36)

with Dyall’s zeroth-order Hamiltonian Hˆ D ,75 and dc,μ = μ |(Hˆ 0 − E)|μ .

(37)

After transforming this vector back to the original excitation basis, b = Xb , it is orthogonalized to the previous subspace vectors to fulfill (b[n+1] )† b[k] = δn+1,k ∀ k ≤ n. Following Davidson’s algorithm,76 the eigenvalue problem is then solved in the subspace to give the solution r[n+1] . Starting from a suitable initial guess r[1] , this procedure is iterated until convergence. In the special case that the initial and final states have the same spatial and spin symmetry, the response vector c (38) r0 = 0 will always be a valid solution with the eigenvalue ω = 0, since the eigenvalue equations then reduce to the ic-MRCC equations for the initial state. Although the response vectors for the excited states may be nonorthogonal to r0 , r = r⊥ + xr0 ,

with x = r0 · r,

(39)


134108-7

Samanta et al.

J. Chem. Phys. 140, 134108 (2014)

the component parallel to r0 does not contribute to the eigenvalue due to Ar0 = 0 (A − ω)r = (A − ω)r⊥ − xωr0 .

(40)

Therefore, we orthogonalize both Ar and r to r0 in each iteration in order to avoid the trivial solution ω = 0. This is equivalent to solving the equation †

X† (1 − r0 r0 )(A − ω)r⊥ = 0

(41)

and will give r⊥ instead of r as eigenvectors, but the correct excitation energies. We note that a similar issue occurs in MkMRCC response theory and is there dealt with in the same manner.49 IV. RESULTS AND DISCUSSION A. Methylene

In the following, we first investigate a variety of vertical excitation energies of singlet methylene, with special focus on the performance of the ic-MRCCSD-LR method for single- and double-excitation dominated states. Being sufficiently small, the CH2 molecule allows to compare our results to full CI as a reference. This molecule is of particular interest because its low-lying excited states comprise both single- and double-replacement dominated excitations with respect to the dominant configuration of the lowest singlet state. Within the present study, the excitations from 1 1 A1 to 2 1 A1 , 2 1 B1 , 3 1 B1 , 2 1 B2 , and 2 1 A2 are double-replacement dominated while the remaining excitations are single-replacement dominated.77 As reference, we use the full CI and single-reference coupled-cluster results of Koch et al.77 The calculations for CH2 are done assuming C2v symmetry with the coordinates C(0,0,0), H(0,±1.644403,1.32213) in atomic units and using the cc-pVDZ basis set of Dunning78 with spherical Gaussians and augmented with diffuse functions. The augmentation consists of one s function with exponent 0.015 for C and one s function with exponent 0.025 for H. All orbitals are included in the correlated methods. We applied the ic-MRCC-LR method in connection with two active spaces, CAS(2,2) and CAS(6,6). While the minimal active space comprises the orbitals 3a1 and 1b1 , the full-valence active space contains in addition the orbitals 2a1 , 4a1 , 1b2 , and 2b2 . In order to demonstrate the impact of the commutator approximation [i.e., the truncation of the residual equations, Eqs. (4) and (17), at the commutator level Ncom ], we report icMRCC results for both Ncom = 2 and Ncom = 3. In case of the larger CAS(6,6), it is also necessary to discard small metric eigenvalues, and we report the results for two choices of the threshold η. The results are collected in Table II. The singlereplacement dominated states are accurately described by both single-reference CCSD and CC3, with errors below 0.02 eV. The latter method includes triply excited clusters to the leading order in perturbation theory.79 The ic-MRCC results using Ncom = 2 and η = 10−5 are somewhat less accurate for some of these states. In particular for the larger active space, errors up to 0.1 eV are found

(for the 4 A1 state). The main reason for this seems to be the truncation of the commutator approximation, as inclusion of threefold commutators (Ncom = 3) significantly reduces these errors to 0.03 eV. Additionally, the excitation energies for the CAS(6,6) are also sensitive to the threshold η for discarding small metric eigenvalues. Notably, an increase of the threshold to η = 10−4 brings down the errors of single-replacement dominated transitions close to those obtained with CCSD and CC3. This is also supported by the root mean square errors computed for these states. Hald et al.80 have analyzed the accuracy of singlereference coupled-cluster methods for single- and doubleexcitation dominated excitations. While CC3 is correct through third-order in perturbation theory for singleexcitation dominated states, an equally good description of the double-excitation dominated states requires to go beyond CCSDT. This is reflected by the much worse performance of CC3 for the doubly excited states of CH2 as shown in Table II. In ic-MRCC-LR, however, excitations involving active orbitals are effectively described by lower excitation ranks or by a response of the reference function. From this, we can expect an increased accuracy of ic-MRCCSD over CC3 for double-excitation dominated states. As obvious from Table II, this is indeed the case. While for CCSD, the errors are in the range 1.5–2.4 eV and in the range 0.5–1.2 eV for CC3, the ic-MRCC-LR method reduces the errors to 0.2–1.2 eV for CAS(2,2) and to 0.05–0.7 for CAS(6,6). In particular for the low-lying doubly excited state 2 1 A1 , the error with CAS(6,6) is as small as for the singly exited states, giving a balanced description of all states up to about 7 eV! For energetically higher-lying doubly excited states, the performance worsens somewhat, which is correlated with the observation of a lower weight from the internal part of the excitation vector rc for these states. This indicates that these states require a larger active space for their improved description. Following a suggestion by a reviewer, we also compare to results from CR-EOMCC methods36, 37 in Table S1 of the supplementary material.81 The ic-MRCCSD-LR and CREOMCC(2,3) methods perform rather similarly for the lowlying excitations, for higher-lying double excitations (where ic-MRCCSD-LR deteriorates as the excitations lie outside the model space) the CR-EOMCC(2,3) results are better. Unfortunately, this accuracy does not transfer to larger systems, as results for oligoenes show (see the supplementary material81 for Sec. IV D). A similar comparison of the excitation energy for three low-lying states of CH2 is provided in Table S2 of the supplementary material.81 Unlike the previous one, we have compared here the adiabatic excitation energy and have used an augmented TZ2P basis set for which FCI results82 are available. Due to the choice of a bigger basis, the difference of results as obtained with our ic-MRCC-LR method and of the other ones is likely to be more indicative of the relative performance. For the double excitation dominated state 21 A1 , the ic-MRCC-LR method gives an error of 0.025 eV using CAS(2,2) as the active space. This error is smaller than that of the single-reference based methods EOM-CCSD (1.1 eV) and CR-EOMCC(2,3)83 (0.1 eV).


134108-8

Samanta et al.

J. Chem. Phys. 140, 134108 (2014)

TABLE II. Vertical excitation energies for CH2 , given in eV. The structure [C(0,0,0), H(0,±1.644403,1.32213)] and the basis set (cc-pVDZ augmented with one s function with exponent 0.015 for C and one s function with exponent 0.025 for each H) are taken from Ref. 77. ic-MRCCSD CAS(6,6) Ncom = 2

CAS(2,2)

4.718 6.505 8.457 1.798 9.360 11.100 7.708 8.141 5.857 10.142

0.065 − 0.065 − 0.104 0.008 0.426 0.502 − 0.057 0.120 0.017 0.653

0.069 − 0.029 − 0.045 0.004 0.456 0.548 − 0.016 0.126 0.004 0.720

0.049 − 0.011 − 0.032 0.007 0.314 0.402 0.003 0.108 0.014 0.671

0.062 − 0.009 − 0.022 0.005 0.454 0.548 0.005 0.124 0.003 0.732

0.061 0.420

0.025 0.457

0.017 0.381

0.011 0.461

Ncom = 3

4.656 6.514 8.479 1.793 8.906 10.553 7.704 8.016 5.853 9.410

6.112 6.509 8.460 1.780 10.705 12.376 7.714 9.619 5.859 11.826

5.127 6.509 8.474 1.788 9.478 11.168 7.720 8.533 5.859 10.583

4.835 6.485 8.438 1.799 9.607 11.327 7.688 8.342 5.839 10.635

4.834 6.496 8.450 1.797 9.600 11.321 7.699 8.349 5.838 10.634

4.721 6.449 8.375 1.803 9.332 11.055 7.647 8.136 5.870 10.063

Difference with respect to FCI 0.180 0.178 − 0.030 − 0.019 − 0.041 − 0.029 0.006 0.004 0.701 0.694 0.774 0.769 − 0.016 − 0.005 0.325 0.333 − 0.014 − 0.015 1.225 1.225

2 1 A1 3 1 A1 4 1 A1 1 1 B1 2 1 B1 3 1 B1 1 1 B2 2 1 B2 1 1 A2 2 1 A2

d s s s d d s d s d

1.456 − 0.005 − 0.019 − 0.014 1.799 1.823 0.011 1.603 0.005 2.416

0.471 − 0.005 − 0.005 − 0.005 0.571 0.615 0.016 0.517 0.005 1.173

RMSEc RMSEc

s d

0.012 1.849

0.009 0.717

c

4.705 6.503 8.447 1.800 9.220 10.954 7.706 8.125 5.867 10.081

Ncom = 2

d s s s d d s d s d

b

4.725 6.486 8.434 1.798 9.362 11.100 7.687 8.142 5.857 10.129

CC3b

2 1 A1 3 1 A1 4 1 A1 1 1 B1 2 1 B1 3 1 B1 1 1 B2 2 1 B2 1 1 A2 2 1 A2

a

η = 10−4

CCSDb

State

0.025 0.739

0.017 0.737

Ncom = 3 10−5

η=

10−5

FCIb

Typea

η=

10−4

η=

Character of excitation with respect to the Hartree-Fock determinant: s = single-excitation dominated, d = double-excitation dominated. Reference 77. Root mean square error.

B. Excited singlet states of p-benzyne

Among all three isomers of benzyne, p-benzyne has the most prominent multireference character in its ground state. For this reason, p-benzyne has been studied extensively using multireference theories, mainly with a focus on its singlettriplet splitting.42, 84–93 A recent study by Jagau and Gauss48 has reported vertical excitation energies for several excited states obtained with Mk-MRCCSD-LR. These results have shown that spurious roots occur in Mk-MRCCSD-LR due to its overcomplete excitation manifold. In our study, we have used the structure from Ref. 48, which has been obtained by optimizing the ground state at the Mk-MRCCSD/cc-pCVTZ level of theory. We have used the same orbital basis, cc-pCVTZ, and correlated all electrons to compare our results with Mk-MRCCSD-LR. The active orbitals used in the study are 6ag and 5b3u (with the molecule oriented such that the z axis is perpendicular to the molecular plane and the x axis runs through the radical centers). The CI coefficients of the configurations within the active space, φ 1 = |6ag α 6ag β and φ 2 = |5b3u α 5b3u β, are c1 = −0.53 and c2 = 0.85, according to ic-MRCCSD theory. Table III compares the results from the methods CCSD-LR, Mk-MRCC-LR, and ic-MRCC-LR. We note that the elements of the ic-MRCC re-

sponse vectors refer to the original excitation basis (not the orthogonalized basis). Within each spatial symmetry, the states are ordered according to their ic-MRCC excitation energies. The excitation energies from CCSD and Mk-MRCC are assigned to those states based on a comparison of the response vectors. Note that this assignment leads to deviations from an energetic ordering of the Mk-MRCC states in the symmetries B2u , B1u , and Au . For excited states in which the dominant configuration does not come from an excitation either from or to an active orbital (1 1 B2u , 2 1 B3u , 4 1 B1u , 3 1 Au ), the excitation energies lie within 0.1 eV for all three methods, CCSD, Mk-MRCCSD, and ic-MRCCSD. Among those roots for which the eigenvectors have dominant contributions from core-to-active or active-to-virtual excitations, the results from Mk-MRCC-LR and ic-MRCC-LR differ to a significant extent. Whenever the dominant excitations are defined with respect to the first reference function (i.e., the one with the smaller weight) in MkMRCCSD (2 1 B2u , 1 1 B1u , 2 1 B1u , 3 1 B2g , 1 1 Au , 2 1 Au ), the excitation energies obtained from ic-MRCC-LR are about 0.4–0.9 eV lower than the corresponding Mk-MRCCSD excitation energies. For states with dominant excitations with respect to the second reference function, we do not observe a particular trend. While for most of these states, the


134108-9

Samanta et al.

J. Chem. Phys. 140, 134108 (2014)

TABLE III. Vertical singlet excitation energies of p-benzyne, given relative to the 1 Ag ground state in eV, using the structure from Ref. 48 and the cc-pCVTZ basis set. The CAS(2,2) in the MRCC calculations comprises the active orbitals 6ag and 5b3u . The Hartree-Fock reference function in CCSD corresponds to the dominant configuration in the active space, φ 2 = |5b3u α 5b3u β. Mk-MRCCSD-LRa State

CCSDa

E

2 1 Ag

6.082

5.089

3 1 Ag

8.566

9.009

1 1 B2u

5.355

5.417

2 1 B2u

...

8.188

3 1 B2u

7.192

7.525

4 1 B2u 1 1 B3u

8.582 4.046

9.321 3.827

...

4.831

2 1 B3u

6.710

6.711

1 1 B1g 2 1 B1g 1 1 B1u

6.514 8.180 ...

7.150 8.864 4.836

2 1 B1u 3 1 B1u

... 5.728

7.641 6.365

4 1 B1u

8.604

8.616

1 1 B3g 2 1 B3g 1 1 B2g 2 1 B2g 1B e 2g

3.133 4.370 2.941 5.401 ...

3.754 5.068 3.556 6.094 6.578

3 1 B2g 1 1 Au 2 1 Au 1A e u 3 1 Au

... 8.413 ... ... 9.058d

... 4.917 7.034 9.142 9.162d

1B

3u

e

ic-MRCCSD-LR

LR vectors dominated by

E

rc (φ 1 ) rc (φ 2 ) 1b3g /5b3u → 6ag /1au (φ 2 )c 1b3g /5b3u → 1au /6ag (φ 2 )c 1b2g → 1au (φ 2 ) 1b3g → 2b1u (φ 2 ) 1b2g → 1au (φ 1 ) 1b3g → 2b1u (φ 1 ) 3b1g → 5b3u (φ 1 ) 3b2u /6ag → 5b3u /5b3u (φ 1 )d

0.834b 0.521b 0.473 0.226 0.466 −0.322 −0.272 0.206 0.657 0.130d

4.489

1b3g → 2b1u (φ 2 ) 4b2u → 6ag (φ 2 ) 1b2g → 1au (φ 2 ) 1b3g → 2b1u (φ 1 ) 3b2u → 6ag (φ 2 ) 5b3u → 6ag (φ 2 ) 6ag → 5b3u (φ 1 ) 1b2g → 2b1u (φ 2 )

0.406 0.306 0.288 −0.248 0.619 0.487 −0.272 −0.256

7.608

6ag → 5b3u (φ 1 ) 5b3u → 6ag (φ 2 ) 5ag → 5b3u (φ 1 ) 1b3g → 1au (φ 2 ) 1b3g → 1au (φ 1 ) 1b2g → 2b1u (φ 2 ) 3b1g → 6ag (φ 2 ) 5b3u → 5b2u (φ 2 ) 1b2g → 5b3u (φ 1 )

0.522 0.309 0.211 0.500 −0.312 0.278 0.659 0.635 0.666

...

6ag → 2b1u (φ 1 ) 1b1u → 6ag (φ 2 ) 1b2g /5b3u → 6ag /6ag (φ 2 ) 3b1g → 1au (φ 2 ) 3b1g → 1au (φ 1 ) 1b3g → 6ag (φ 2 ) 5b3u → 1au (φ 2 ) 1b2g → 6ag (φ 2 ) 5b3u → 2b1u (φ 2 ) 1b2g /6ag → 5b3u /5b3u (φ 1 ) 1b1u → 5b3u (φ 1 )

0.642 0.595 0.304 0.559 −0.340 0.669 0.665 0.671 0.664 0.498 0.456

1b3g → 5b3u (φ 1 ) 6ag → 1au (φ 1 ) 1b3g /5b3u → 6ag /6ag (φ 2 ) 1b2g → 5b2u (φ 2 )d 1b2g → 5b2u (φ 1 )d

0.674 0.655 0.609 0.562d −0.327d

7.882 5.368

LR vectors dominated by rc (φ 1 ) rc (φ 2 ) 1b3g /5b3u → 6ag /1au c 1b3g /5b3u → 1au /6ag c 1b2g → 1au 1b3g → 2b1u

0.712b 0.449b 0.418 0.206 0.535 0.390

3b1g → 5b3u 3b2u → 6ag 1b3g → 2b1u 1b2g → 1au 3b1g → 5b3u 1b3g → 2b1u 4b2u → 6ag 1b2g → 1au 3b2u → 6ag rc (|6ag α 5b3u β) rc (|6ag β 5b3u α) 5ag → 5b3u 1b2g → 2b1u ...

0.855 −0.402 0.252 −0.171 0.640 −0.385 0.344 0.270 0.633 0.592 −0.592 0.359 −0.149

6.661

1b3g → 1au 1b2g → 2b1u

0.604 −0.295

6.754 8.596 4.120

8.688

3b1g → 6ag 5b3u → 5b2u 1b2g → 5b3u 1b1u → 6ag 6ag → 2b1u 1b1u → 6ag 1b2g → 5b3u 3b1g → 1au

0.744 0.762 1.107 0.323 1.062 0.689 −0.488 0.640

3.496 4.900 3.234 5.849 ...

1b3g → 6ag 5b3u → 1au 1b2g → 6ag 5b3u → 2b1u ...

0.784 0.776 0.767 0.764

8.338 4.516 6.502 ... 9.138

1b1u → 5b3u 1b3g → 5b3u 6ag → 1au

1.138 1.223 1.180

1b2g → 5b2u

0.642

7.339

9.126 4.159

6.776 7.451

a

Reference 48. Note that the true response vector rc may be shifted along the direction of the ground state vector c. c The corresponding distribution of spin functions is α/β → α/β. d Reference 110. e Artificial root arising in Mk-MRCCSD due to the problem of overcompleteness of excited space. b


134108-10

Samanta et al.

J. Chem. Phys. 140, 134108 (2014)

TABLE IV. Comparison of vertical ic-MRCCSD-LR excitation energies (in eV) with those obtained using direct (state-specific) ic-MRCC and selected conventional methods. Those states of CH2 and p-benzyne are considered, which can be described within the same CAS(2,2) as the reference state in the linear response treatment, i.e., 1 1 A1 and 1 1 Ag , respectively. Structures and basis sets are the same as in Tables II and III. Values given in parentheses are based on the CASSCF orbitals of the reference state and thus do not contain orbital relaxation effects. ic-MRCC / CAS(2,2) State

CASSCF (CASCI)

CCSD

CC3

2 1 A1 1 1 B1 1 3 B1

5.320 (14.677) 2.026 (8.815) −0.091 (3.021)

6.112a 1.780a

5.127a 1.788a

8.130 (11.788) 7.366 (11.595) 0.082 (0.123)

6.082b 4.046b

4.245c 3.622c

FCI methylene 4.656a 1.793a −0.082

SD-LR

Direct SD

Direct SD(T)

4.835 1.799 −0.068

4.670 (4.745) 1.805 (1.783) −0.077 (−0.097)

4.646 (4.646) 1.795 (1.757) −0.081 (−0.115)

4.949 (5.027) 4.593 (4.730) 0.237 (0.241)

4.105 (3.932) 3.843 (3.724) 0.285 (0.284)

p-benzyne 2 1 Ag 1 1 B3u 1 3 B3u

4.489 4.159 0.261

a

Reference 77. Reference 48. c Computed with the program package C FOUR (Ref. 111). b

methods agree within 0.5 eV, two notable exceptions exist: For the states 3 1 Ag and 3 1 B1u , the excitation energies differ by more than 1 eV. No additional root has appeared in our ic-MRCC-LR study of p-benzyne. This is also expected, because the set of excited functions in ic-MRCC theory does not contain any redundancy. Our results also agree on the presence of some additional roots in Mk-MRCC-LR as pointed out by Jagau and Gauss,48 but not on all of them. Those states which we recognize as additional root in Mk-MRCC-LR are denoted by leaving a blank space in the beginning of a row in Table III. One Mk-MRCC root of each of the symmetries B3u , B2g , and Au is assigned as additional root. As these roots come in pairs (due to the choice of the active space), we consider the one closer in energy to the ic-MRCC root as “proper” and the other one as “spurious.” The pair of corresponding roots with B3u symmetry is dominated by a single excitation from one active orbital to the other and vice versa. The reason for the appearance of a spurious root for Mk-MRCCLR in this case is that the model space response is described by internal excitation operators rather than by a response of the reference function. This is done in order to prevent a symmetry-conditioned decoupling between closed-shell and open-shell determinants through the eigenvalue equations.48 In ic-MRCC-LR, however, there is no obstacle to describing this state through a response of the CI coefficients in the reference function. The additional roots of B2g and Au symmetry, on the other hand, get their main contributions from double excitations that are linearly dependent with single excitations acting on the other reference function. The large admixture of double excitations in the description of these additional roots is likely the reason for the much higher energies of these roots compared to their counterparts. In Ref. 48, further Mk-MRCC-LR roots are reported as spurious, namely, the ones assigned by us as 2 1 B2u , 1 1 B1u , 2 1 B1u , and 2 1 Au . However, all of these roots have ic-MRCCLR counterparts, and a comparison of the response vectors with those of other roots of the same symmetry suggests that they indeed describe distinct physical states. These roots have been mainly considered as spurious because of the absence

of corresponding roots in the CCSD-LR calculations.48 However, we deem it likely that the CCSD counterparts of these roots occur at a much higher energy, since their main contributing excitations are single excitations with respect to the first reference function and thus correspond to double excitations with respect to the Hartree-Fock determinant. In conclusion, ic-MRCC-LR avoids the difficulties due to spurious roots as encountered in Mk-MRCC-LR and is certainly also more accurate than the corresponding singlereference treatment of biradicalic systems such as p-benzyne.

C. Direct ic-MRCC vs. linear response approach

Any state that can be described within the model space may either be directly targeted by solving the ic-MRCC equations, or it may be obtained by the linear response formalism starting from a different initial state. In the following, we compare both approaches, focusing on those low-lying states of methylene and p-benzyne that can be described within the same CAS(2,2) as the initial states 1 1 A1 and 1 1 Ag , respectively. Table IV compares the ic-MRCCSD-LR method with the direct (state-specific) methods ic-MRCCSD and ic-MRCCSD(T)94 and also reports results from CASSCF and the traditional EOM-CC approaches CCSD and CC3. We have employed the same structures and basis sets as in Secs. IV A and IV B, for methylene and p-benzyne, respectively. For methylene, however, a similar study with a triplezeta basis set is available in Table S2 of the supplementary material.81 In case of methylene, Table IV reveals very small deviations (

Communication: automatic code generation enables nuclear gradient computations for fully internally contracted multireference theory.

Linear response theory for the density matrix renormalization group: efficient algorithms for strongly correlated excited states.

Valence Excited States in Large Molecules via Local Multireference Singles and Doubles Configuration Interaction.

First-order derivative couplings between excited states from adiabatic TDDFT response theory.

A new size extensive multireference perturbation theory.

Free electrons and ionic liquids: study of excited states by means of electron-energy loss spectroscopy and the density functional theory multireference configuration interaction method.

Relativistic GVVPT2 multireference perturbation theory description of the electronic states of Y2 and Tc2.

Configuration interaction singles natural orbitals: an orbital basis for an efficient and size intensive multireference description of electronic excited states.

Modelling excited states of weakly bound complexes with density functional theory.

Identifying the Hamiltonian structure in linear response theory.

Wavefunction in density functional theory embedding for excited states: which wavefunctions, which densities?

Constricted Variational Density Functional Theory Approach to the Description of Excited States.

Charge-transfer excited states in aqueous DNA: Insights from many-body Green's function theory.

Spectroscopy of inorganic excited states.

Description of spin-orbit coupling in excited states with two-component methods based on approximate coupled-cluster theory.

Excited-State Vibrations of Solvated Molecules: Going Beyond the Linear-Response Polarizable Continuum Model.

Determining polarizable force fields with electrostatic potentials from quantum mechanical linear response theory.

Polarizable embedding with a multiconfiguration short-range density functional theory linear response method.

Excited quartet states in DNA photolyase.

Excited states of positronic lithium and beryllium.

RNA bases.

Electronically excited states of PANH anions.

Predictable valence excited states of anions.

Charge-displacement analysis for excited states.