Targeting excited states in all-trans polyenes with electron-pair states
TTargeting excited states in all-trans polyenes with electron-pair states
Katharina Boguslawski ∗ Institute of Physics, Faculty of Physics,Astronomy and Informatics,Nicolaus Copernicus University in Torun,Grudziadzka 5, 87-100 Torun, PolandandFaculty of Chemistry, Nicolaus Copernicus University in Torun,Gagarina 7, 87-100 Torun, Poland (Dated: May 1, 2019)Wavefunctions restricted to electron pair states are promising models for strongly-correlated sys-tems. Specifically, the pair Coupled Cluster Doubles (pCCD) ansatz allows us to accurately de-scribe bond dissociation processes and heavy-element containing compounds with multiple quasi-degenerate single-particle states. Here, we extend the pCCD method to model excited states usingthe equation of motion (EOM) formalism. As the cluster operator of pCCD is restricted to electron-pair excitations, EOM-pCCD allows us to target excited electron-pair states only. To model singlyexcited states within EOM-pCCD, we modify the configuration interaction ansatz of EOM-pCCDto contain also single excitations. Our proposed model represents a simple and cost-effective alter-native to conventional EOM-CC methods to study singly excited electronic states. The performanceof the excited state models is assessed against the lowest-lying excited states of the uranyl cationand the two lowest-lying excited states of all-trans polyenes. Our numerical results suggest thatEOM-pCCD including single excitations is a good starting point to target singly excited states.
I. INTRODUCTION
Quantum chemical modeling of electronic structures ofexcited states is usually more challenging than those ofthe ground state [1]. Specifically, we have to be ableto describe various types of excited states simultane-ously, each governed by different amounts of dynamicand nondynamic/static electron correlation effects [2–4].Furthermore, excited states usually require larger one-electron basis functions that allow us to describe all cor-relation effects on an equal footing. The presence of min-ima, transition states, and surface crossings poses addi-tional challenges when modeling their potential energysurfaces [ ? ]. Therefore, an accurate and reliable de-scription of electronically excited states still remains anopen problem in quantum chemistry.In the past two decades, many promising quantumchemistry approaches suitable for excited states haveemerged (see, for instance, recent reviews [5, 33 ? ? ?? ]). In wavefunction-based methods, the most popu-lar approaches are the multi-reference configuration in-teraction, complete active space self-consistent field the-ory [6, 7], including dynamic energy corrections [8–12],and coupled cluster ans¨atze. Specifically, in coupled clus-ter theory, considerable development was done in theFock-Space coupled cluster (FSCC) formulation [13, 14]and the equation of motion formalism [15–19]. In thestandard formulation of the FSCC approach, one startsfrom a closed-shell electronic configuration described bya standard coupled cluster approach and calculates the ∗ k.boguslawski@fizyka.umk.pl electronically excited states using the information fromionization potentials and electron affinities. Alterna-tively, electronic spectra can be obtained from the high-spin reference wavefunction. Despite its sound theoret-ical basis, the FSCC approach suffers form convergencedifficulties. To remedy this problem, intermediate Hamil-tonian techniques have been developed [20–24], but con-vergence difficulties still remain for large and complexmolecular systems [25]. In the equation of motion cou-pled cluster (EOM-CC) method, one obtains electronicspectra by introducing a linear excitation operator on topof the CC ground state wavefunction. This theory hasbeen successfully applied to various problems in chem-istry and physics [26–31]. However, the standard EOM-CC method does not produce accurate excitation ener-gies for systems with significant multi-reference charac-ter. This deficiency led to the development of its modifiedvariants, such as, the spin-flip EOM-CC [32, 33], com-pletely renormalized EOM-CC [34–36] and active spaceEOM-CC methods [37–42]. Large-scale modeling of ex-cited states using the EOM-CC formalism has been avail-able with the NWChem software package [43, 44]. Sim-plified EOM-CC models represent alternative approachesto describe excited states in large molecular system [45].Another example is the recently formulated equation ofmotion linear coupled cluster ansatz [46].In this work we derive and implement another simpli-fied version of the EOM method to describe electroni-cally excited states in large molecular systems. Specif-ically, we use the pair coupled cluster doubles (pCCD)ansatz [47, 48] to model the initial ground-state electronicstructure. Recent applications of this model to ground-state properties, including bond-breaking processes, havebeen very encouraging [47–60]. The pCCD wavefunction a r X i v : . [ c ond - m a t . s t r- e l ] A p r is a simplified CC-type wavefunction, where the clusteroperator is restricted to electron pair excitations [61, 62], | pCCD (cid:105) = exp (cid:32) P (cid:88) i =1 K (cid:88) a = P +1 c ai a † a ↑ a † a ↓ a i ↓ a i ↑ (cid:33) | (cid:105) , (1)where | (cid:105) is some reference determinant. Indices i and a correspond to occupied and virtual orbitals with respectto | (cid:105) , P and K denote the number of electron pairs andorbitals, respectively. The pCCD ansatz is also known asthe antisymmetric product of 1-reference orbital geminal(AP1roG) ansatz, where the pair-coupled-cluster ampli-tudes correspond to the geminal coefficients. This wave-function ansatz is, by construction, size-extensive and hasmean-field scaling if the geminal coefficients/amplitudesare optimized using the projected Schr¨odinger equationapproach. Note that | (cid:105) is usually optimized as well andhence differs from the Hartree–Fock determinant.This work is organized as follows. In section II, webriefly summarize the equation of motion formalism witha pCCD reference function as well as a simple and cost-effective extension to account for single excitations. Nu-merical examples are presented in section V, where weinvestigated the lowest lying excited states of the UO molecule as a test case and the two lowest-lying excitedstates of all-trans polyenes. Finally, we conclude in sec-tion VI. II. TARGETING EXCITED STATES IN PCCD
In order to model excited states in pCCD, we will usethe equation of motion (EOM) formalism [1, 17, 63]. InEOM-CC, excited states are parametrized by a linear CI-type ansatz [19], ˆ R = (cid:88) µ c µ ˆ τ µ , (2)where the summation is over all excitation operatorspresent in the cluster operator ˆ T as well as the identityoperator ˆ τ . The operator ˆ R is then used to generate thetarget state from the initial CC state, | Ψ (cid:105) = ˆ R exp( ˆ T ) | (cid:105) = (cid:88) µ c µ ˆ τ µ exp( ˆ T ) | (cid:105) , (3)with | (cid:105) being the CC reference determinant.To arrive at the EOM-CC working equations, it is con-venient to use the normal-product form of the Hamilto-nian, ˆ H N = ˆ H −(cid:104) | ˆ H | (cid:105) , written as a sum of the Fock op-erator ˆ F N and the electron repulsion term ˆ W N (in physi-cists’ notation),ˆ H N = (cid:88) pq f pq { a † p a q } + 12 (cid:88) pqrs (cid:104) pq | rs (cid:105){ a † p a † q a s a r } . (4)Furthermore, we will disregard any excitation properties,like dipole moments, and focus on excitation energies in-stead. In that case, we have to solve for the ˆ R amplitudes only. Our target-state Schr¨odinger equation then readsˆ H N ˆ R exp ( ˆ T ) | (cid:105) = ∆ E ˆ R exp ( ˆ T ) | (cid:105) , (5)where ∆ E is the energy difference with respect to theFermi vacuum expectation value | (cid:105) . Introducing thesimilarity transformed Hamiltonian in normal-productform ˆ H N = exp ( − ˆ T ) ˆ H N exp ( ˆ T ) and subtracting theequation for the CC ground state, we obtain the EOM-CC equations for the ˆ R amplitudes,[ ˆ H N , ˆ R ] | (cid:105) = ω ˆ R | (cid:105) , (6)where ω are the excitation energies with respect to theCC ground state, exp( ˆ T ) | (cid:105) . The excitation energies arethus the eigenvalues of a non-Hermitian matrix, (cid:20) (cid:104) | ˆ H N | µ (cid:105) (cid:104) ν | [ ˆ H N , τ µ ] | (cid:105) (cid:21) , (7)where the first row is associated with the CC referencestate and the subsequent rows correspond to the excitedconfigurations ν >
0. The EOM-CC working equationsmay be solved using, for instance, non-Hermitian exten-sions of the Davidson algorithm to determine the lowest-lying excited electronic states.
A. Electron-pair excitations
In this work, we are considering an pCCD refer-ence function | pCCD (cid:105) as a special CC state confined toelectron-pair states. As indicated in eq. (1), the pCCDcluster operator contains only electron-pair excitations,ˆ T = ˆ T p = (cid:80) ia t ai a † a a † ¯ a a ¯ i a i . In the corresponding EOMmodel, the ˆ R operator is thus restricted to the identityoperator ˆ τ as well as all pair excitations present in thecluster operator ˆ T p ,ˆ R p = c ˆ τ + (cid:88) ia c a ¯ ai ¯ i ˆ τ a ¯ ai ¯ i , (8)where ˆ τ a ¯ ai ¯ i = a † a a † ¯ a a ¯ i a i creates an electron pair in thevirtual orbital a . To obtain the target-state Schr¨odingerequation of EOM-pCCD, we have to substitute the gen-eral cluster operator ˆ T by the pair-excitation operator ˆ T p in eq. (5),ˆ H N ˆ R p exp ( ˆ T p ) | (cid:105) = ∆ E ˆ R p exp ( ˆ T p ) | (cid:105) . (9)The ˆ R p amplitudes are determined from the EOM-pCCDequations (restricting ˆ R to ˆ R p and ˆ T to ˆ T p in eq. (6)),[ ˆ H ( p ) N , ˆ R p ] | (cid:105) = ω p ˆ R p | (cid:105) , (10)where ˆ H ( p ) N indicates the similarity transformed Hamil-tonian of pCCD, ˆ H ( p ) N = exp ( − ˆ T P ) ˆ H N exp ( ˆ T P ), and ω p are the electron-pair excitation energies. Thus, EOM-pCCD allows us to model electron-pair excited statesonly. In order to target singly excited or general dou-bly excited states, we have to extend the pCCD clus-ter operator to include excitations beyond electron pairs.This can be done either by changing to a frozen-pairCCSD [64] formalism or to a linearized CCSD correctionwith an pCCD reference function [56]. In this work how-ever, we will consider a different, cost-effective approachto account for single excitations in the pCCD model thatdoes not scale as O ( n ) as conventional EOM-CCSDmethods. Note that EOM-pCCD scales as O ( o v ),where o is the number of occupied orbitals (equivalentto the number of electron pairs) and v is the number ofvirtual orbitals. B. Accounting for single excitations
As proposed by Forseman et al. [65], configuration in-teraction with only single substitutions (CIS) representsan accurate model to investigate (singly) excited elec-tronic states, even for large systems. In the CIS method,the reference is a single Slater determinant obtained froman SCF procedure, while the CIS wavefunction is ex-panded as | CIS (cid:105) = c | (cid:105) + (cid:88) ia c ai | ai (cid:105) , (11)with | ai (cid:105) being a singly excited determinant where the(occupied) orbital i of | (cid:105) has been substituted by the(virtual) orbital a . Similar to CIS, we will include singleexcitations in EOM-pCCD by extending the ˆ R p operatorof eq. (8). In addition to the identity operator ˆ τ and allpair excitations ˆ T p , ˆ R p also contains a summation overall single excitations,ˆ R ps = c ˆ τ + (cid:88) ia c ai ˆ τ ai + (cid:88) ia c a ¯ ai ¯ i ˆ τ a ¯ ai ¯ i , (12)where ˆ τ ai is a singlet excitation operator ˆ τ ai = a † a a i + a † ¯ a a ¯ i that creates a singly excited electronic state with respectto | (cid:105) , | ai (cid:105) = ˆ τ ia | (cid:105) . The ˆ R ps amplitudes are determinedfrom solving [ ˆ H ( p ) N , ˆ R ps ] | (cid:105) = ω ps ˆ R ps | (cid:105) , (13)where we still have the similarity transformed Hamilto-nian of pCCD, ˆ H ( p ) N , while ω ps are the excitation ener-gies of both singly excited and pair excited states. Wewill label this simplified model as EOM-pCCD+S to in-dicate that single excitations are included a posteriori in the ˆ R p operator. In contrast to CIS that uses a singleSlater determinant as reference, EOM-pCCD+S employsthe pCCD wavefunction as reference state. Furthermore,electron correlation effects are included through the ˆ T p operator in the similarity transformed Hamiltonian.Similar to EOM-pCCD, the excitation energies are ob-tained by diagonalizing a non-Hermitian matrix of the form (cid:104) | ˆ H ( p ) N | bj (cid:105) (cid:104) | ˆ H ( p ) N | b ¯ bj ¯ j (cid:105)(cid:104) ai | ˆ H ( p ) N | (cid:105) (cid:104) ai | ˆ H ( p ) N | bj (cid:105) (cid:104) ai | ˆ H ( p ) N | b ¯ bj ¯ j (cid:105) (cid:104) a ¯ ai ¯ i | ˆ H ( p ) N | bj (cid:105) (cid:104) a ¯ ai ¯ i | ˆ H ( p ) N | b ¯ bj ¯ j (cid:105) . (14)Note that in contrast to conventional EOM-CC methods,the first column does not equal zero because single excita-tions are not included in the cluster operator of pCCD.Thus, terms like (cid:104) ai | ˆ H ( p ) N | (cid:105) do not vanish as they arenot incorporated in the ground-state CC amplitude equa-tions. Although we can account for single excitations ina rather straightforward way, we loose size-intensivity inthe EOM model. For the molecular systems investigatedin this work, the error introduced by extending only theˆ R p operator is approximately three orders of magnitudesmaller than the actual excitation energies, while thecomputational cost increases insignificantly compared toEOM-pCCD (the Hamiltonian still contains terms thatscale as O ( o v ), but with a larger pre-factor). Thus,EOM-pCCD+S represents a cost-effective starting pointto study singly excited electronic states in the pCCDmodel. III. COMPUTATIONAL DETAILSA. Molecular structures
In this work, we investigate all-trans polyenes contain-ing 2 to 7 π -bonds (the corresponding Lewis structuresare shown in Figure 1). Specifically, we study five dif-ferent geometries, three for the lowest lying verticallyexcited states and two for lowest lying adiabatically ex-cited states. The first set of structures was optimized foreach ground state using the ADF software package [66–68] with a DZ basis set [69] and the B3LYP exchange–correlation functional [70]. The remaining sets of geome-tries were taken from ref. [71] and were optimized us-ing the density matrix renormalization group (DMRG).Those include molecular geometries for the verticallyand adiabatically excited states of C H , C H , andC H . Furthermore, for both the vertically and adia-batically excited states, we study geometries optimizedusing two different active spaces as discussed in ref. [71]:an active space containing only π -electrons (abbreviatesas v π and a π ) and an active space comprising both σ -and π -electrons (denoted by v and a , respectively).For the linear UO molecule, we calculated the fourlowest-lying vertically and adiabatically excited states.The position of the minimum of each potential energysurface was obtained by simultaneously stretching theU–O bond between 1.610 and 1.870 ˚A. H CH CH CH CH H ( a ) C H H CH CH CH CH CH CH CH CH H ( c ) C H H CH CH CH CH CH CH CH CH CH CH CH CH H ( e ) C H H CH CH CH CH CH CH H ( b ) C H H CH CH CH CH CH CH CH CH CH CH H ( d ) C H H CH CH CH CH CH CH CH CH CH CH CH CH CH CH H ( f ) C H FIG. 1. Lewis structures of all investigated all-trans polyenes.
B. Calculation of ground and excited states
The pCCD, EOM-pCCD, and EOM-pCCD+S cal-culations were performed with our in-house quantum-chemistry code employing restricted canonical Hartree–Fock orbitals as molecular orbital basis. The lowest-lyingexcitation energies were determined employing a mod-ified Davidson algorithm. We should emphasize thatthe molecular orbital basis was not optimized withinthe pCCD method. Orbital optimization in pCCD typi-cally results in localized, symmetry-broken orbitals thatpreclude us from identifying the symmetry of excitedstates. Furthermore, orbital optimization is only per-formed for the pCCD ground state wavefunction and re-sults in molecular orbitals that are biased toward theground state. This bias can be reduced by introducinga state-average protocol yielding a set of compromise or-bitals (as usually done in multi-configuration SCF meth-ods). Since we need to identify the symmetry of the(lowest-lying) excited states, we chose to not optimizethe orbital basis and to use orbitals that transform ac-cording to the irreducible representations of the molecu-lar point group (numerical examples for excitation ener-gies within pCCD obtained in the optimized orbital basiscan be found in the Supporting Information). Moreover,recent numerical studies confirm that pCCD using re-stricted Hartree–Fock orbitals captures, at least qualita-tively, the most important correlation effects around theequilibrium geometry [60].For the investigated all-trans polyenes, we used the6-31G basis set as well as the cc-pVDZ basis set of Dun-ning [72] for direct comparison with either CIS(D) (6-31G) or MRMP/DMRG (cc-pVDZ) reference data, re-spectively.For UO , scalar relativistic effects were incorporatedthrough relativistic effective core potentials (RECP). Inall calculations, we have used a small core (SC) RECP(60 electrons in the core) with the following contractionscheme (12s11p10d8f) → [8s7p6d4f] [73]. For the lighterelements (O), the cc-pVDZ basis set of Dunning [72] wasemployed, (10s5p1d) → [4s3p1d]. The CIS and EOM-CCSD calculations have been performed with the Mol- pro software package [74]. No frozen core was used inall calculations.
IV. THE UO
MOLECULE AS A TEST CASE
Before scrutinizing the performance of EOM-pCCDand EOM-pCCD+S in modeling the lowest-lying excitedstates in all-trans polyenes, we will assess the accuracy ofEOM-pCCD+S against molecular systems whose lowest-lying excited states are purely singly-excited states. Forthat purpose, we chose the uranyl cation (UO ) as itcontains a large number of strongly-correlated electronsdistributed among the 5 f -, 6 d -, and 7 s -orbitals and itselectronic spectrum is well understood [75–80].The vertical and adiabatic excitation energies for thefour lowest-lying excited states in EOM-pCCD+S andCIS are summarized in Table I. Note that for the UO molecule the two lowest-lying excited states can be ac-curately described within the EOM-CCSD model, whichyields excitations energies that are similar to completelyrenormalized EOM-CCSD(T) reference values [30]. Forthe higher-lying excited states, however, dynamic corre-lation effects become important and a triples correctionhas to be included to accurately model those states. TheEOM-CCSD results in Table I can thus be considered asupper bounds of the excitation energies for the π u → δ u state.In general, EOM-pCCD+S overestimates vertical ex-citations energies of the two lowest-lying excited states( σ u → φ u and σ u → δ u ) by approximately 0.5 eV, whilethe corresponding adiabatic excitation energies are over-estimated by about 0.7 eV. Note that for the uranylcation, CIS outperforms EOM-pCCD+S and deviatesfrom EOM-CCSD reference data by approximately 0.3(vertical excitations) to 0.6 eV (adiabatic excitations).As expected, neither EOM-pCCD+S nor CIS are able toaccurately predict the excitation energies for the π u → δ u states as these excited states are dominated by dynamiccorrelation effects which are not included in CIS and onlymarginally accounted for in pCCD (see also [60]). Weshould emphasize that EOM-pCCD+S provides equilib-rium U–O bond lengths of excited states that deviateless from the EOM-CCSD reference values (∆ r e ≈ . r e ≈ .
08 ˚A in CIS).To conclude, the EOM-pCCD+S model predicts thecorrect order of the lowest-lying excited states in theuranyl cation and provides excitation energies of decentaccuracy with errors of about 0.5 eV with respect toEOM-CCSD reference values. However, errors in exci-tation energies are slightly worse than in the simple CISmodel.
V. EXCITED STATES IN ALL-TRANSPOLYENES
All-trans polyenes are model systems for carotenoidsand polyene chromophores that play an important rolein photoprocesses. The proper description of the twolowest-lying excited states poses a challenge to both ex-periment and quantum chemistry approaches [81, 85–97],especially because doubly excited configuration are re-quired to accurately model ground and excited states oflonger polyenes. The excitation energies for the two-lowest vertical excited states of C H to C H usingthe DFT optimized structures and two different basis sets(cc-pVDZ and 6-31G) are summarized in Table II. Theexcitation energies obtained in the 6-31G basis are givenin parenthesis next to the corresponding results calcu-lated for the cc-pVDZ basis set. Note that the 6-31Gbasis set was used in CIS(D) [81], while a cc-pVDZ ba-sis set was employed in MRMP calculations [82]. SinceEOM-pCCD can only model pair excited states, only theexcitation energies of the 2 A − g state are shown in theTable.For all investigated polyenes, EOM-pCCD overesti-mates the excitation energies of the first excited state2 A − g by approximately 4 eV compared to MRMP and1.5 eV compared to CIS(D) data. This error can be par-tially attributed to missing single excitations that haveto be included in the model to describe the character ofthe 2 A − g state accurately. Including singles on top ofEOM-pCCD a posteriori reduces the error by approxi-mately 2 to 1 eV, depending on the number of π bonds.While EOM-pCCD+S predicts excitation energies thatare lower than the corresponding excitation energies ofCIS(D) (differences amount up to 1.7 eV), it overesti-mates the excitation energies of the 2 A − g state by about2 eV compared to MRMP data. Furthermore, the simpleEOM-pCCD+S model does not predict the right order ofstates where the first dark state 2 A − g should lie belowthe first bright state 1 B + u . The excitation energies ofthe first bright state 1 B + u , however, deviate only by ap-proximately 0.3 eV from MRMP data. We should notethat EOM-pCCD+S outperforms CIS(D) in predictingthe excitation energies of the 1 B + u state. Despite itssimplicity, EOM-pCCD+S is thus a good starting pointto model singly excited states.We can extrapolate the excitation energies of the two lowest-lying excited states for longer polyenes in the limitof infinite number of π bonds and predict the ordering ofstates for longer polyenes (see, for instance, refs. [81, 82,86, 88]). To do so, the excitation energies can be fittedto a particle-in-a-box model [81] of the form E ( n π ) = aL ( n π ) − b + c, (15)where L ( n π ) is the length of the polyene that dependson the number of π bonds n π , while a, b, c are the fittingparameters. The function for the length of the polyeneis obtained by linear regression of the terminal C–C dis-tance of the DFT optimized structures and defines thedistance between the terminal carbon atoms as a func-tion of the number of double bonds, L ( n π ) = 2 . n π − . . (16)The excitation energy at infinite number of π bonds isobtained from the fitting parameter c . For EOM-pCCD,the extrapolated excitation energy of the 2 A − g state is5.61 eV and reduces to 4.28 eV for EOM-pCCD+S. Forthe 1 B + u state, EOM-pCCD+S results in an extrapo-lated excitation energy of 3.06 eV. The resulting fits aswell as excitation energies are shown in Figure 2. Notethat EOM-pCCD+S model overestimates the extrapo-lated excitation energy of the 1 B + u state by approx-imately 1.30 eV compared to experiment [97] and byabout 0.94 eV compared to MRMP reference data (2.12eV) [82]. E xc i t a t i on ene r g y [ e V ] Number of π bonds A g− (EOM−pCCD+S)B u+ (EOM−pCCD+S)A g− (EOM−pCCD) FIG. 2. Extrapolation of excitation energies calculated forthe EOM-pCCD and EOM-pCCD+s models for the first darkstate 2 A − g (orange line) and the first bright state 1 B + u (blueline; EOM-pCCD+S only). Finally, we will discuss how the molecular structure af-fects excitation energies and their character in the EOM-pCCD+S method. Note, however, that the labels “ v ”and “ a ” in Table III correspond to the vertically andadiabatically excited state in DMRG calculations. Theydo not represent the true adiabatic and vertical excita-tion energies in EOM-pCCD+S as the DMRG-optimized TABLE I. Vertical and adiabatic excitation energies and equilibrium U–O bond lengths of the four lowest-lying excited statesin the UO molecule calculated for EOM-pCCD+S and different quantum chemistry methods.method vertical adiabatic ω [eV] ω [eV] r e [˚A] σ u → φ u σ u → δ u π u → δ u π u → δ u σ u → φ u σ u → δ u π u → δ u π u → δ u σ u → φ u σ u → δ u π u → δ u π u → δ u EOM-pCCD+S 4 .
45 4 .
83 7 .
79 7 .
93 4 .
43 4 .
82 7 .
63 7 .
80 1 .
711 1 .
709 1 .
746 1 . .
33 4 .
79 7 .
84 8 .
00 4 .
18 4 .
66 7 .
47 7 .
62 1 .
694 1 .
690 1 .
719 1 . .
02 4 .
36 5 .
28 5 .
36 3 .
70 4 .
08 4 .
63 4 .
72 1 .
772 1 .
768 1 .
805 1 . H to C H calculated forEOM-pCCD, EOM-pCCD+S, and different quantum chemistry methods. Note that the 6-31G basis set was used in CIS(D),while the cc-pVDZ basis set was utilized in MRMP. CASSCF was performed in a double-zeta basis set (see correspondingreferences). The excitation energies of EOM-pCCD and EOM-pCCD+S are determined for the cc-pVDZ basis set, while thecorresponding results for the 6-31G basis set are given in parenthesis.2 A − g B + u C −− C EOM-pCCD EOM-pCCD+S CIS(D) a MRMP b CASSCF c EOM-pCCD+S CIS(D) a MRMP b CASSCF c .
56 (10 .
49) 7 .
45 (7 .
37) 9 .
01 6 .
31 6 .
67 7 .
20 (7 .
44) 8 .
09 6 .
21 7 .
733 9 .
11 (9 .
04) 6 .
79 (6 .
75) 7 .
81 5 .
10 5 .
64 5 .
98 (6 .
16) 6 .
78 5 .
25 7 .
064 8 .
11 (8 .
02) 6 .
15 (6 .
09) 6 .
78 4 .
26 5 .
16 5 .
19 (5 .
34) 5 .
95 4 .
57 6 .
625 7 .
42 (7 .
32) 5 .
69 (5 .
63) 6 .
12 3 .
68 4 .
32 4 .
62 (4 .
75) 5 .
43 4 .
17 6 .
376 6 .
93 (6 .
82) 5 .
37 (5 .
29) 5 .
55 3 .
19 – 4 .
20 (4 .
31) 5 .
00 3 .
87 –7 6 .
58 (6 .
46) 5 .
13 (5 .
05) 5 .
14 2 .
80 – 3 .
87 (3 .
97) 4 .
70 3 .
60 – a Taken from Ref. [81] b Taken from Ref. [82] c Taken from Ref. [83] structures are used. To simplify the labeling of the ex-cited states and to facilitate direct comparison to DMRG,we will use the same labeling (“ v ” and “ a ”) in DMRGand EOM-pCCD+S for all investigated excited states.Figure 3 and Table III show the two lowest-lying excita-tion energies of C H , C H , and C H for differ-ent molecular geometries optimized by DMRG. For the2 A − g state, the excitation energies predicted by EOM-pCCD+S reduce by approximately 1 eV if the molec-ular structure is allowed to relaxed (adiabatic excita-tions). The choice of the active space in the DMRG struc-ture optimization, however, does not significantly affectthe accuracy of EOM-pCCD+S for longer polyene chainlengths (differences are less than 0.1 eV between v ( a )and v π ( a π )), while DMRG is more sensitive to the activespace (differences in DMRG excitation energies amountto 1 eV). Similar observation can be made for the 1 B + u state. For relaxed molecular structures, the excitationenergies slightly improve by 0.2 eV compared to the ver-tical excitation energies. Furthermore, the active spaceused in DMRG calculations does not change the excita-tion energies considerably and similar excitation energiesare obtained for the vertically and adiabatically excitedstates, respectively, with EOM-pCCD+S. We should em-phasize that the excitation energies of the 1 B + u statepredicted by EOM-pCCD+S are closer to experimentalvalues than the corresponding DMRG excitation ener-gies. Differences between the theoretically determined values and experiment can be partly attributed to thebasis set used in calculations as well as missing electroncorrelation effects not included in the theoretical model.Figure 4 shows the contributions (in percent) of thethree dominant configurations of the 2 A − g excited statefor different molecular geometries. These are the dou-bly excited b g → a ∗ u (HOMO → LUMO ), the singlyexcited b g → b ∗ g (HOMO → LUMO+1), and the singly ex-cited a u → a ∗ u (HOMO − → LUMO) [86, 88]. For theDFT optimized structures shown in Figure 4(a) (verti-cal excitations), the contribution from the doubly-excitedconfiguration HOMO → LUMO increases to approxi-mately 50-60%, while both singly excited configurationsconstitute 15 to 25%. For longer polyene chain lenghts,the contribution from the doubly-excited configurationdecreases.The character of the 2 A − g state changes consider-ably if the molecular structure is modified. Specifi-cally, the weights of the configurations of the verticallyexcited 2 A − g state strongly depend on the number of π -bonds if DMRG optimized structures are used (seesolid lines in Figure 4(b)). While the dark state inC H is dominated by the HOMO → LUMO configu-ration for the shorter polyene chain lengths, the weight ofthe HOMO − → LUMO configuration gradually increasesand constitutes more than 60 % for C H . Extendingthe polyene chain to 7 double bonds further decreasesthe weight of the HOMO → LUMO configuration to TABLE III. Vertical and adiabatic excitation energies of the two lowest-lying excited states in all-trans polyenes C H to C H calculated with EOM-pCCD, EOM-pCCD+S, and DMRG for different DMRG-optimized geometries. Note thatdifferent active spaces are used in DMRG calculations (see computational details and ref. [71]), while all orbitals are active inEOM-pCCD and EOM-pCCD+S. The DMRG reference data is taken from ref. [71]. Experimental data is taken from ref. [84].2 A − g v a v π a π C −− C EOM-pCCD+S DMRG EOM-pCCD+S DMRG EOM-pCCD+S DMRG EOM-pCCD+S DMRG Exp.5 6 .
28 5 .
43 5 .
18 4 .
01 6 .
07 4 .
51 4 .
85 3 .
36 3 .
036 5 .
99 4 .
76 4 .
61 3 .
41 5 .
84 4 .
15 4 .
60 2 .
99 2 .
697 5 .
76 4 .
64 4 .
44 3 .
22 5 .
63 3 .
91 4 .
45 2 .
73 2 . B + u v a v π a π EOM-pCCD+S DMRG EOM-pCCD+S DMRG EOM-pCCD+S DMRG EOM-pCCD+S DMRG Exp.5 4 .
91 5 .
35 5 .
18 4 .
98 4 .
79 5 .
77 4 .
62 5 .
49 3 .
576 4 .
50 4 .
98 4 .
30 4 .
60 4 .
41 5 .
41 4 .
26 5 .
13 3 .
317 4 .
25 4 .
66 4 .
04 4 .
29 4 .
14 5 .
16 3 .
99 4 .
87 3 . E xc i t a t i on ene r g y [ e V ] Number of π bonds DMRG A g− (v)DMRG B u+ (v)DMRG A g− (a)DMRG B u+ (a)DMRG A g− (v π )DMRG B u+ (v π )DMRG A g− (a π )DMRG B u+ (a π )Experiment A g− Experiment B u+ ( a ) DMRG reference E xc i t a t i on ene r g y [ e V ] Number of π bonds EOM−pCCD+S A g− (v)EOM−pCCD+S B u+ (v)EOM−pCCD+S A g− (a)EOM−pCCD+S B u+ (a)EOM−pCCD+S A g− (v π )EOM−pCCD+S B u+ (v π )EOM−pCCD+S A g− (a π )EOM−pCCD+S B u+ (a π )Experiment A g− Experiment B u+ ( b ) EOM-pCCD+S FIG. 3. Excitation energies in eV of the two lowest-lying excited states for C H , C H , and C H calculated withEOM-pCCD+S. The molecular structures used in each calculation are given in parentheses and explained in the text (seecomputational details). The DMRG reference values are taken from ref. [71]. less than 10 %, while the contribution of the singly ex-cited HOMO → LUMO+1 configuration slightly increases.Thus the vertically excited 2 A − g state of C H is dom-inated by only one configuration.The picture changes completely if the molecu-lar structure of the excited state is allowed to re-lax (see dashed lines in Figure 4(b)). For grow-ing polyene chains, the weight of the doubly excitedHOMO → LUMO configuration gradually increases (upto 70 % for C H ), while both singly excited config-urations (HOMO → LUMO+1) and HOMO − → LUMO)approach weights of approximately 15 %. AlthoughEOM-pCCD+S predicts a dominant contribution of thedoubly-excited HOMO → LUMO configuration in the2 A − g excited state, it underestimates the overall weight of all doubly-excited configurations ( <
70 % in EOM-pCCD+S) that may increase to approximately 80 % forlonger polyene chain lenghts [81]. This underestimationcan be attributed to broken-pair configurations that arenot included in the pCCD model and have to be ac-counted for a posteriori [56, 64]. To conclude, EOM-pCCD+S can appropriately describe the three most im-portant configurations present in the dark state (see also[81, 86, 88]) if molecular structures are relaxed. The qual-itatively good description of the character of the firstdark state is also evident in the excitation energies. Theadiabatic excitation energies are approximately 1 to 1.5eV lower in energy than the corresponding vertical ones(compare Tables II and III). E xc i t a t i on c on t r i bu t i on [ % ] Number of π bonds HOMO → LUMO HOMO−1 → LUMOHOMO → LUMO+1 ( a ) DFT optimized structures E xc i t a t i on c on t r i bu t i on [ % ] Number of π bonds HOMO → LUMO (v)HOMO−1 → LUMO (v)HOMO → LUMO+1 (v)HOMO → LUMO (a)HOMO−1 → LUMO (a)HOMO → LUMO+1 (a) ( b ) DMRG optimized structures FIG. 4. Contributions (in percent) of the three dominant con-figurations of the 2 A − g state for the (a) DFT optimized struc-tures (vertical excitations) and (b) DMRG optimized struc-tures (vertical and adiabatic excitations). VI. CONCLUSIONS
In this work, we have extended the pCCD modelto excited states using the equation of motion formal-ism. Since the cluster operator of pCCD is restrictedto electron-pair states, only electron-pair excitations canbe modeled with EOM-pCCD. To augment the excited-state model with single and general double excitations,we have to modify the cluster operator to include singleand general double excitations as well which consider-ably increases the computational cost compared to thepCCD reference calculation. To arrive at a cost-effectivemodel similar to pCCD, we have included single exci-tation a posteriori in the CI-type ansatz. Our method(EOM-pCCD+S) is motivated by the CIS approach thatallows us to efficiently model singly excited states withan SCF reference function. Since only the CI-type ansatzis modified, the computational scaling remains similarto EOM-pCCD (the Hamiltonian contains terms thatscale as O ( o v )), but with a larger pre-factor. Despitethe simplicity of the model, EOM-pCCD+S breaks size-intensivity, which can be restored by adjusting the clusteroperator.Although pCCD is commonly combined with anorbital-optimization protocol, we did not optimize themolecular orbital basis in the pCCD reference calcula-tion. Orbital optimization in pCCD usually results in lo-calized, symmetry-broken orbitals that prevent us fromidentifying the symmetry of the targeted excited states.Furthermore, if the orbitals are optimized within pCCD,the optimal set of orbitals is biased toward the groundstate and the corresponding excitation energies are gen-erally overestimated. This bias can be reduced by opti-mizing the orbitals simultaneously for both ground andexcited states, which increases the computational cost ofEOM-pCCD and EOM-pCCD+S significantly. To ben- efit from the O ( o v ) scaling of both excitation models,we have thus skipped the orbital optimization step andused restricted Hartree–Fock orbitals in all calculations.Both excitation models have been assessed againstexcitation energies of the uranyl cation, whose lowest-lying excited states contain purely singly-excited states,and all-trans polyenes containing 2 to 7 π -bonds. Forthe UO molecule, EOM-pCCD+S provides excita-tion energies that are close to CIS results, overesti-mates, however, the excitation energies of the two lowest-lying excited states by about 0.5 eV compared to EOM-CCSD reference calculations. For all-trans polyenes,EOM-pCCD+S fails to predict the correct order ofthe first dark and first bright state, while the exci-tation energies of the 1 B + u state are in good agree-ment with MRMP results and closer to experimentthan DMRG reference data. Furthermore, the excita-tion energies and the character of the excited statesstrongly depend on the molecular structures used incalculations. Specifically, the character of the 2 A − g state can only be properly predicted if molecular struc-tures are allowed to relax, resulting in a dominantdoubly-excited HOMO → LUMO configuration andtwo singly-excited configurations (HOMO → LUMO+1)and HOMO − → LUMO). In order to accurately modelthe first dark state in all-trans polyenes, double excita-tions beyond electron-pair excitations as well as higherexcitations (triples, etc.) might be important to cap-ture the missing correlation effects in the targeted ex-cited states that cannot be described within pCCD. Pos-sible extension of EOM-pCCD using CC-type correctionson top of the pCCD reference function are presently be-ing developed in our laboratory. To conclude, EOM-pCCD+S represents a good and cost-effective startingpoint to investigate singly-excited states in the pCCDmodel. Additional numerical studies are currently underinvestigation in our laboratory.
SUPPLEMENTARY MATERIAL
See Supplementary Material for excitation energies ofselected all-trans polyenes obtained by EOM-pCCD andEOM-pCCD+S in the optimized pCCD orbital basis.
ACKNOWLEDGEMENTS
Financial support from a SONATA BISgrant of the National Science Centre, Poland(no. 2015/18/E/ST4/00584) and a Marie-Sk(cid:32)lodowska-Curie Individual Fellowship project no. 702635–PCCDXis gratefully acknowledged.Calculations have been carried out using resources pro-vided by Wroclaw Centre for Networking and Supercom-puting (http://wcss.pl), grant No. 412.We had many helpful discussions with Leszek Meissner,Piotr Piecuch, and Karol Kowalski. [1] R. J. Bartlett and M. Musia(cid:32)l, Rev. Mod. Phys. , 291(2007).[2] T. J. Lee and P. R. Taylor, Int. J. Quantum Chem. ,199 (1989).[3] R. J. Bartlett and J. F. Stanton, Rev. Comput. Chem. , 65 (1994).[4] K. Boguslawski, P. Tecmer, O. Legeza, and M. Reiher,J. Phys. Chem. Lett. , 3129 (2012).[5] P. Piecuch, K. Kowalski, I. S. Pimienta, and M. J.Mcguire, Inter. Rev. Phys. Chem. , 527 (2002).[6] B. Roos and P. R. Taylor, Chem. Phys. , 157 (1980).[7] P. E. M. Siegbahn, J. Alml¨of, A. Heiberg, and B. O.Roos, J. Chem. Phys. , 2384 (1981).[8] J. Finley, P.-A. Malmqvist, B. O. Roos, , and L. Serrano-Andr´es, Chem. Phys. Lett. , 299 (1998).[9] K. Andersson, P.-A. Malmqvist, B. O. Roos, A. J. Sadlej,and K. Woli´nski, J. Phys. Chem. , 5483 (1990).[10] C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger,and J.-P. Malrieu, J. Chem. Phys. , 10252 (2001).[11] P. Pulay, Int. J. Quantum Chem. , 3273 (2011).[12] P.-A. Malmqvist, K. Pierloot, A. R. M. Shahi, C. J.Cramer, and L.Gagliardi, J. Chem. Phys. , 204109(2008).[13] L. Meissner and M. Musia(cid:32)l, in Recent Progress in CoupledCluster Methods (Springer, 2010) p. 395.[14] D. I. Lyakh, M. Musia(cid:32)l, V. F. Lotrich, and . J. Bartlett,Chem. Rev. , 182 (2012).[15] D. J. ROWE, Rev. Mod. Phys. , 153 (1968).[16] T. H. Dunning and V. McKoy, J. Chem. Phys. , 1735(1967).[17] J. Jankowski, K. Kowalski, and P. Jankowski,Chem. Phys. Lett. , 608 (1994).[18] I. Shavitt and R. J. Bartlett, Many-body methods inchemistry and physics (Cambridge University Press,2009).[19] T. Helgaker, P. Jørgensen, and J. Olsen,
MolecularElectronic-Structure Theory (Wiley, 2000).[20] A. Landau and E. Eliav, J. Chem. Phys. , 6862(2000).[21] A. Landau, E. Eliav, Y. Ishikawa, and U. Kaldor,J. Chem. Phys. , 6862 (2001).[22] L. Visscher, E. Eliav, and U. Kaldor, J. Chem. Phys. , 9720 (2001).[23] L. Meissner, Int. J. Quantum Chem. , 2199 (2007).[24] M. Musia, J. Chem. Phys. , 134111 (2012).[25] P. Tecmer, H. van Lingen, A. S. P. Gomes, and L. Viss-cher, J. Chem. Phys. , 084308 (2012).[26] K. Kowalski, S. Hirata, M. Wloch, P. Piecuch, and T. L.Windus, J. Chem. Phys. , 074319 (2005).[27] K. Kowalski and P. Piecuch, J. Chem. Phys. , 643(2001).[28] K. Kowalski, R. M. Olson, S. Krishnamoorthy, V. Tip-paraju, and E. Apra, J. Chem. Theory Comput. , 2200(2011).[29] K. Lopata, R. Reslan, M. Kowaska, D. Ne/uhauser,N. Govind, and K. Kowalski, J. Chem. Theory Com-put. , 3686 (2011).[30] P. Tecmer, N. Govind, K. Kowalski, W. A. de Jong., andL. Visscher, J. Chem. Phys , 034301 (2013).[31] E. Epifanovsky, K. Klein, S. Stopkowicz, J. Gauss, andA. I. Krylov, J. Chem. Phys. , 064102 (2015). [32] S. V. Levchenko and A. I. Krylov, J. Chem. Phys. ,175 (2004).[33] A. I. Krylov, Acc. Chem. Res. , 83 (2006).[34] K. Kowalski and P. Piecuch, J. Chem. Phys. , 1715(2004).[35] M. Wloch, J. R. Gour, K. Kowalski, and P. Piecuch,J. Chem. Phys. , 214107 (2005).[36] K. Kowalski, S. Krishnamoorthy, O. Villa, J. R. Ham-mond, and N. Govind, J. Chem. Phys. , 154103(2010).[37] K. Kowalski and P. Piecuch, J. Chem. Phys. , 8490(2000).[38] K. Kowalski, S. Hirata, M. W(cid:32)loch, P. Piecuch, and T. L.Windus, J. Chem. Phys. , 074319 (2005).[39] P. Piecuch, Mol. Phys. , 2987 (2010).[40] K. Kowalski, S. Krishnamoorthy, O. Villa, J. R. Ham-mond, and N. Govind, J. Chem. Phys. , 154103(2010).[41] J. Shen and P. Piecuch, J. Chem. Phys. , 194102(2013).[42] J. Shen and P. Piecuch, J. Chem. Theory Comput. ,4968 (2012).[43] K. Kowalski, S. Krishnamoorthy, R. M. Olson, V. Tippa-raju, and E. Apra, in (IEEE, 2011) pp. 1–10.[44] H. J. van Dam, W. A. De Jong, E. Bylaska, N. Govind,K. Kowalski, T. Straatsma, and M. Valiev, WIRES Com-put. Mol. Sci. , 888 (2011).[45] S. R. Gwaltney, M. Nooijen, and R. J. Bartlett,Chem. Phys. Lett. , 189 (1996).[46] J. N. Byrd, V. Rishi, A. Perera, and R. J. Bartlett,J. Chem. Phys. , 164103 (2015).[47] P. A. Limacher, P. W. Ayers, P. A. Johnson, S. DeBaerdemacker, D. Van Neck, and P. Bultinck, J. Chem.Theory Comput. , 1394 (2013).[48] K. Boguslawski, P. Tecmer, P. W. Ayers, P. Bultinck,S. De Baerdemacker, and D. Van Neck, Phys. Rev. B , 201106(R) (2014).[49] P. Tecmer, K. Boguslawski, P. A. Limacher, P. A. John-son, M. Chan, T. Verstraelen, and P. W. Ayers, J. Phys.Chem. A , 9058 (2014).[50] P. A. Limacher, T. D. Kim, P. W. Ayers, P. A. Johnson,S. De Baerdemacker, D. Van Neck, and P. Bultinck, Mol.Phys. , 853 (2014).[51] K. Boguslawski, P. Tecmer, P. A. Limacher, P. A. John-son, P. W. Ayers, P. Bultinck, S. De Baerdemacker, andD. Van Neck, J. Chem. Phys. , 214114 (2014).[52] K. Boguslawski, P. Tecmer, P. W. Ayers, P. Bultinck,S. De Baerdemacker, and D. Van Neck, J. Chem. TheoryComput. , 4873 (2014).[53] P. Tecmer, K. Boguslawski, and P. W. Ayers, Phys.Chem. Chem. Phys. , 14427 (2015).[54] P. A. Limacher, J. Chem. Theory Comput. , 3629(2015).[55] P. Limacher, P. Ayers, P. Johnson, S. De Baerdemacker,D. Van Neck, and P. Bultinck, Phys. Chem. Chem. Phys , 5061 (2014).[56] K. Boguslawski and P. W. Ayers, J. Chem. Theory Com-put. , 5252 (2015).[57] A. J. Garza, A. G. S. Alencar, and G. E. Scuseria, J. Chem. Phys. , 244106 (2015).[58] T. M. Henderson, I. W. Bulik, and G. E. Scuseria,J. Chem. Phys. , 214116 (2015).[59] J. A. Gomez, T. M. Henderson, and G. E. Scuseria,J. Chem. Phys. , 244117 (2016).[60] K. Boguslawski, P. Tecmer, and ¨O. Legeza, Phys. Rev.B , 155126 (2016).[61] T. M. Henderson, G. E. Scuseria, J. Dukelsky, A. Signo-racci, and T. Duguet, Phys. Rev. C , 054305 (2014).[62] T. Stein, T. M. Henderson, and G. E. Scuseria, J. Chem.Phys. , 214113 (2014).[63] J. F. Stanton and R. J. Bartlett, J. Chem. Phys. , 7029(1993).[64] T. M. Henderson, I. W. Bulik, T. Stein, and G. E. Scuse-ria, J. Chem. Phys. , 244104 (2014).[65] J. B. Foresman, M. Head-Gordon, J. A. Pople, and M. J.Frisch, J. Phys. Chem. , 135 (1992).[66] (2014), ADF2014.01, SCM, Theoretical Chemistry,Vrije Universiteit, Amsterdam, The Netherlands, .[67] G. te Velde, F. M. Bickelhaupt, S. J. A. van Gisber-gen, C. F. Guerra, E. J. Baerends, J. G. Snijders, andT. Ziegler, J. Comput. Chem. , 931 (2001).[68] C. F. Guerra, J. G. Snijders, G. te Velde, and E. J.Baerends, Theor. Chem. Acc. , 391 (1998).[69] E. van Lenthe and E. J. Baerends, J. Comput. Chem. ,1142 (2003).[70] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, andM. J. Frisch, J. Phys. Chem. , 11623 (1994).[71] W. Hu and G. K.-L. Chan, J. Chem. Theory Comput. , 3000 (2015).[72] T. H. Dunning, J. Chem. Phys. , 1007 (1989).[73] W. K¨uchle, M. Dolg, H. Stoll, and H. Preuss, Mol. Phys. , 1245 (1991).[74] H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, P. C.M. Sch¨utz, T. Korona, A. Mitrushenkov, G. Rauhut,T. B. Adler, R. D. Amos, A. Bernhardsson, A. Bern-ing, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn,F. Eckert, E. Goll, C. Hampel, G. Hetzer, T. Hrenar,G. Knizia, C. K¨oppl, Y. Liu, A. W. Lloyd, R. A. Mata,A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura,A. Nicklass, P. Palmieri, K. Pfl¨uger, R. Pitzer, M. Rei-her, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni,T. Thorsteinsson, M. Wang, and A. Wolf, “Molpro, ver-sion 2012.1, a package of Ab Initio , 6880(1999).[76] K. Pierloot and E. van Besien, J. Phys. Chem. ,204309 (2005).[77] R. G. Denning, J. Phys. Chem. A , 4125 (2007).[78] F. R´eal, V. Vallet, C. Marian, and U. Wahlgren, J. Phys.Chem. , 214302 (2007).[79] P. Tecmer, A. S. P. Gomes, U. Ekstr¨om, and L. Visscher,Phys. Chem. Chem. Phys. , 6249 (2011).[80] P. Tecmer, R. Bast, K. Ruud, and L. Visscher,J. Phys. Chem. A , 7397 (2012).[81] J. H. Starcke, M. Wormit, J. Schirmer, and A. Dreuw,Chem. Phys. , 39 (2006).[82] Y. Kurashige, H. Nakano, Y. Nakao, and K. Hirao,Chem. Phys. Lett. , 425 (2004).[83] K. Nakayama, H. Nakano, and K. Hirao, Int. J. Quan-tum Chem. , 157 (1998).[84] B. E. Kohler, J. Chem. Phys. , 2788 (1988).[85] B. Hudson and B. Kohler, Chem. Phys. Lett. , 299(1972).[86] K. Schulten, I. Ohmine, and M. Karplus, J. Chem. Phys. , 4422 (1976).[87] R. J. Cave and E. R. Davidson, J. Chem. Phys. , 4481(1987).[88] P. Tavan and K. Schulten, Phys. Rev. B , 4337 (1987).[89] R. J. Cave and E. R. Davidson, J. Chem. Phys. , 614(1988).[90] W. J. Buma, B. E. Kohler, and K. Song, J. Chem. Phys. , 6367 (1991).[91] G. Orlandi, F. Zerbetto, and M. Z. Zgierski, Chem. Rev. , 867 (1991).[92] L. Serrano-Andres, M. Merchan, I. Nebot-Gil, R. Lindh,and B. O. Roos, J. Chem. Phys. , 3151 (1993).[93] L. Serrano-Andres, R. Lindh, B. O. Roos, and M. Mer-chan, J. Chem. Phys. , 9360 (1993).[94] R. P. Krawczyk, K. Malsch, G. Hohlneicher, R. C. Gillen,and W. Domcke, Chem. Phys. Lett. , 535 (2000).[95] C.-P. Hsu, S. Hirata, and M. Head-Gordon,J. Phys. Chem. A , 451 (2001).[96] K. B. Wiberg, A. E. de Oliveira, , and G. Trucks,J. Phys. Chem. A , 4192 (2002).[97] R. L. Christensen, M. G. I. Galinato, E. F. Chu,J. N. Howard, R. D. Broene, and H. A. Frank,J. Phys. Chem. A112