Targeting excited states in all-trans polyenes with electron-pair states

TL;DR

This paper extends the pCCD method to model excited states in strongly-correlated systems, offering a simpler and cost-effective alternative to traditional methods, with promising results on polyenes and uranyl cation.

Contribution

The authors develop an EOM-pCCD method with single excitations to efficiently target singly excited states in strongly-correlated molecules.

Findings

EOM-pCCD with single excitations accurately describes excited states.

The method performs well on polyenes and uranyl cation.

It offers a computationally cheaper alternative to EOM-CC methods.

Abstract

Wavefunctions restricted to electron pair states are promising models for strongly-correlated systems. Specifically, the pair Coupled Cluster Doubles (pCCD) ansatz allows us to accurately describe 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 using the 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 singly excited states within EOM-pCCD, we modify the configuration interaction ansatz of EOM-pCCD to contain also single excitations. Our proposed model represents a simple and cost-effective alternative to conventional EOM-CC methods to study singly excited electronic states. The performance of the excited state models is assessed against the lowest-lying excited states of the uranyl cation and the two lowest-lying excited states of all-trans polyenes. Our numerical results suggest that EOM-pCCD including single excitations is a good starting point to target singly excited states.

Figures

Tables

Table 1. Table 1: Vertical and adiabatic excitation energies and equilibrium U–O bond lengths of the four lowest-lying excited states in the \ce UO2^2+ molecule calculated for EOM-pCCD+S and different quantum chemistry methods.

method	vertical				adiabatic
	$\omega$ [eV]				$\omega$ [eV]				$r_e$ [Å]
	${\sigma }_u\to {\varphi }_u$	${\sigma }_u\to {\delta }_u$	${\pi }_u\to {\delta }_u$	${\pi }_u\to {\delta }_u$	${\sigma }_u\to {\varphi }_u$	${\sigma }_u\to {\delta }_u$	${\pi }_u\to {\delta }_u$	${\pi }_u\to {\delta }_u$	${\sigma }_u\to {\varphi }_u$	${\sigma }_u\to {\delta }_u$	${\pi }_u\to {\delta }_u$	${\pi }_u\to {\delta }_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.744$
CIS	$4.33$	$4.79$	$7.84$	$8.00$	$4.18$	$4.66$	$7.47$	$7.62$	$1.694$	$1.690$	$1.719$	$1.720$
EOM-CCSD	$4.02$	$4.36$	$5.28$	$5.36$	$3.70$	$4.08$	$4.63$	$4.72$	$1.772$	$1.768$	$1.805$	$1.805$

Table 2. Table 2: Vertical excitation energies of the two lowest-lying excited states in all-trans polyenes \ce C4H6 to \ce C14H16 calculated for EOM-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 corresponding references). The excitation energies of EOM-pCCD and EOM-pCCD+S are determined for the cc-pVDZ basis set, while the corresponding results for the 6-31G basis set are given in parenthesis.

	$2^1{A}_g^{-}$					$1^1{B}_u^{+}$
\ceC=C	EOM-pCCD	EOM-pCCD+S	CIS(D)111Taken from Ref. Starcke et al. (2006)	MRMP222Taken from Ref. Kurashige et al. (2004)	CASSCF333Taken from Ref. Nakayama et al. (1998)	EOM-pCCD+S	CIS(D)111Taken from Ref. Starcke et al. (2006)	MRMP222Taken from Ref. Kurashige et al. (2004)	CASSCF333Taken from Ref. Nakayama et al. (1998)
2	$10.56\left(10.49\right)$	$7.45\left(7.37\right)$	$9.01$	$6.31$	$6.67$	$7.20\left(7.44\right)$	$8.09$	$6.21$	$7.73$
3	$9.11\left(9.04\right)$	$6.79\left(6.75\right)$	$7.81$	$5.10$	$5.64$	$5.98\left(6.16\right)$	$6.78$	$5.25$	$7.06$
4	$8.11\left(8.02\right)$	$6.15\left(6.09\right)$	$6.78$	$4.26$	$5.16$	$5.19\left(5.34\right)$	$5.95$	$4.57$	$6.62$
5	$7.42\left(7.32\right)$	$5.69\left(5.63\right)$	$6.12$	$3.68$	$4.32$	$4.62\left(4.75\right)$	$5.43$	$4.17$	$6.37$
6	$6.93\left(6.82\right)$	$5.37\left(5.29\right)$	$5.55$	$3.19$	–	$4.20\left(4.31\right)$	$5.00$	$3.87$	–
7	$6.58\left(6.46\right)$	$5.13\left(5.05\right)$	$5.14$	$2.80$	–	$3.87\left(3.97\right)$	$4.70$	$3.60$	–

Table 3. Table 3: Vertical and adiabatic excitation energies of the two lowest-lying excited states in all-trans polyenes \ce C10H12 to \ce C14H16 calculated with EOM-pCCD, EOM-pCCD+S, and DMRG for different DMRG-optimized geometries. Note that different active spaces are used in DMRG calculations (see computational details and ref.  Hu and Chan ( 2015 ) ), while all orbitals are active in EOM-pCCD and EOM-pCCD+S. The DMRG reference data is taken from ref.  Hu and Chan ( 2015 ) . Experimental data is taken from ref.  Kohler ( 1988 ) .

	$2^1{A}_g^{-}$
	$v$		$a$		$v_{\pi }$		$a_{\pi }$
\ceC=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.03$
6	$5.99$	$4.76$	$4.61$	$3.41$	$5.84$	$4.15$	$4.60$	$2.99$	$2.69$
7	$5.76$	$4.64$	$4.44$	$3.22$	$5.63$	$3.91$	$4.45$	$2.73$	$2.44$
	$1^1{B}_u^{+}$
	$v$		$a$		$v_{\pi }$		$a_{\pi }$
	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.57$
6	$4.50$	$4.98$	$4.30$	$4.60$	$4.41$	$5.41$	$4.26$	$5.13$	$3.31$
7	$4.25$	$4.66$	$4.04$	$4.29$	$4.14$	$5.16$	$3.99$	$4.87$	$3.12$

Equations

$$|\mathrm{pCCD}⟩=\exp \left(\sum _{i=1}^P\sum _{a=P+1}^K{c}_i^a{a}_{a\uparrow }^{†}{a}_{a\downarrow }^{†}{a}_{i\downarrow }{a}_{i\uparrow }\right)|0⟩,$$

$$\hat{R}=\sum _{\mu }{c}_{\mu }{\hat{\tau }}_{\mu },$$

$$|\Psi ⟩=\hat{R}\exp (\hat{T})|0⟩=\sum _{\mu }{c}_{\mu }{\hat{\tau }}_{\mu }\exp (\hat{T})|0⟩,$$

$${\hat{H}}_N=\sum _{pq}{f}_{pq}\{a_p^{†}{a}_q\}+\frac{1}{2}\sum _{pqrs}⟨pq|rs⟩\{a_p^{†}{a}_q^{†}{a}_s{a}_r\}.$$

$${\hat{H}}_N\hat{R}\exp (\hat{T})|0⟩=\Delta E\hat{R}\exp (\hat{T})|0⟩,$$

$$[{\hat{\mathcal{H}}}_N,\hat{R}]|0⟩=\omega \hat{R}|0⟩,$$

$$\left[\begin{array}{cc}0& ⟨0|{\hat{\mathcal{H}}}_N|\mu ⟩\\ 0& ⟨\nu |[{\hat{\mathcal{H}}}_N,{\tau }_{\mu }]|0⟩\end{array}\right],$$

$${\hat{R}}_p=c_0{\hat{\tau }}_0+\sum _{ia}{c}_{i\bar{i}}^{a\bar{a}}{\hat{\tau }}_{a\bar{a}i\bar{i}},$$

$${\hat{H}}_N{\hat{R}}_p\exp ({\hat{T}}_p)|0⟩=\Delta E{\hat{R}}_p\exp ({\hat{T}}_p)|0⟩.$$

$$[{\hat{\mathcal{H}}}_N^{(p)},{\hat{R}}_p]|0⟩={\omega }_p{\hat{R}}_p|0⟩,$$

$$|\mathrm{CIS}⟩=c_0|0⟩+\sum _{ia}{c}_i^a{|}_i^a⟩,$$

$${\hat{R}}_{ps}=c_0{\hat{\tau }}_0+\sum _{ia}{c}_i^a{\hat{\tau }}_{ai}+\sum _{ia}{c}_{i\bar{i}}^{a\bar{a}}{\hat{\tau }}_{a\bar{a}i\bar{i}},$$

$$[{\hat{\mathcal{H}}}_N^{(p)},{\hat{R}}_{ps}]|0⟩={\omega }_{ps}{\hat{R}}_{ps}|0⟩,$$

$$\left[\begin{array}{ccc}0& ⟨0|{\hat{\mathcal{H}}}_N^{(p)}{|}_j^b⟩& ⟨0|{\hat{\mathcal{H}}}_N^{(p)}{|}_{j\bar{j}}^{b\bar{b}}⟩\\ {⟨}_i^a|{\hat{\mathcal{H}}}_N^{(p)}|0⟩& {⟨}_i^a|{\hat{\mathcal{H}}}_N^{(p)}{|}_j^b⟩& {⟨}_i^a|{\hat{\mathcal{H}}}_N^{(p)}{|}_{j\bar{j}}^{b\bar{b}}⟩\\ 0& {⟨}_{i\bar{i}}^{a\bar{a}}|{\hat{\mathcal{H}}}_N^{(p)}{|}_j^b⟩& {⟨}_{i\bar{i}}^{a\bar{a}}|{\hat{\mathcal{H}}}_N^{(p)}{|}_{j\bar{j}}^{b\bar{b}}⟩\end{array}\right].$$

$$E(n_{\pi })=\frac{a}{L(n_{\pi })-b}+c,$$

$$L(n_{\pi })=2.4709n_{\pi }-\mathrm{1.2398.}$$

Full text

I Introduction

Quantum chemical modeling of electronic structures of excited states is usually more challenging than those of the ground state Bartlett and Musiał (2007). Specifically, we have to be able to describe various types of excited states simultaneously, each governed by different amounts of dynamic and nondynamic/static electron correlation effects Lee and Taylor (1989); Bartlett and Stanton (1994); Boguslawski et al. (2012). Furthermore, excited states usually require larger one-electron basis functions that allow us to describe all correlation effects on an equal footing. The presence of minima, transition states, and surface crossings poses additional challenges when modeling their potential energy surfaces SerranoAndrés200599. Therefore, an accurate and reliable description of electronically excited states still remains an open problem in quantum chemistry.

In the past two decades, many promising quantum chemistry approaches suitable for excited states have emerged (see, for instance, recent reviews Piecuch et al. (2002); SerranoAndrés200599; excitations-review; Krylov (2006); Roos_rev_2008; nishikawa2012preface). In wavefunction-based methods, the most popular approaches are the multi-reference configuration interaction, complete active space self-consistent field theory Roos and Taylor (1980); Siegbahn et al. (1981), including dynamic energy corrections Finley et al. (1998); Andersson et al. (1990); Angeli et al. (2001); Pulay (2011); Malmqvist et al. (2008), and coupled cluster ansätze. Specifically, in coupled cluster theory, considerable development was done in the Fock-Space coupled cluster (FSCC) formulation Meissner and Musiał (2010); Lyakh et al. (2012) and the equation of motion formalism ROWE (1968); Dunning and McKoy (1967); Jankowski et al. (1994); Shavitt and Bartlett (2009); Helgaker et al. (2000). In the standard formulation of the FSCC approach, one starts from a closed-shell electronic configuration described by a standard coupled cluster approach and calculates the electronically excited states using the information from ionization potentials and electron affinities. Alternatively, electronic spectra can be obtained from the high-spin reference wavefunction. Despite its sound theoretical basis, the FSCC approach suffers form convergence difficulties. To remedy this problem, intermediate Hamiltonian techniques have been developed Landau and Eliav (2000); Landau et al. (2001); Visscher et al. (2001); Meissner (2007); Musiał (2012), but convergence difficulties still remain for large and complex molecular systems Tecmer et al. (2012a). In the equation of motion coupled cluster (EOM-CC) method, one obtains electronic spectra by introducing a linear excitation operator on top of the CC ground state wavefunction. This theory has been successfully applied to various problems in chemistry and physics Kowalski et al. (2005a); Kowalski and Piecuch (2001); Kowalski et al. (2011a); Lopata et al. (2011); Tecmer et al. (2013); Epifanovsky et al. (2015). However, the standard EOM-CC method does not produce accurate excitation energies for systems with significant multi-reference character. This deficiency led to the development of its modified variants, such as, the spin-flip EOM-CC Levchenko and Krylov (2004); Krylov (2006), completely renormalized EOM-CC Kowalski and Piecuch (2004); Wloch et al. (2005); Kowalski et al. (2010a) and active space EOM-CC methods Kowalski and Piecuch (2000); Kowalski et al. (2005b); Piecuch (2010); Kowalski et al. (2010b); Shen and Piecuch (2013, 2012). Large-scale modeling of excited states using the EOM-CC formalism has been available with the NWChem software package Kowalski et al. (2011b); van Dam et al. (2011). Simplified EOM-CC models represent alternative approaches to describe excited states in large molecular system Gwaltney et al. (1996). Another example is the recently formulated equation of motion linear coupled cluster ansatz Byrd et al. (2015).

In this work we derive and implement another simplified version of the EOM method to describe electronically excited states in large molecular systems. Specifically, we use the pair coupled cluster doubles (pCCD) ansatz Limacher et al. (2013); Boguslawski et al. (2014a) to model the initial ground-state electronic structure. Recent applications of this model to ground-state properties, including bond-breaking processes, have been very encouraging  Limacher et al. (2013); Boguslawski et al. (2014a); Tecmer et al. (2014); Limacher et al. (2014a); Boguslawski et al. (2014b, c); Tecmer et al. (2015); Limacher (2015); Limacher et al. (2014b); Boguslawski and Ayers (2015); Garza et al. (2015); Henderson et al. (2015); Gomez et al. (2016); Boguslawski et al. (2016). The pCCD wavefunction is a simplified CC-type wavefunction, where the cluster operator is restricted to electron pair excitations Henderson et al. (2014a); Stein et al. (2014),

[TABLE]

where $|0\rangle$ is some reference determinant. Indices $i$ and $a$ correspond to occupied and virtual orbitals with respect to $|0\rangle$, $P$ and $K$ denote the number of electron pairs and orbitals, respectively. The pCCD ansatz is also known as the antisymmetric product of 1-reference orbital geminal (AP1roG) ansatz, where the pair-coupled-cluster amplitudes correspond to the geminal coefficients. This wavefunction ansatz is, by construction, size-extensive and has mean-field scaling if the geminal coefficients/amplitudes are optimized using the projected Schrödinger equation approach. Note that $|0\rangle$ is usually optimized as well and hence differs from the Hartree–Fock determinant.

This work is organized as follows. In section II, we briefly summarize the equation of motion formalism with a pCCD reference function as well as a simple and cost-effective extension to account for single excitations. Numerical examples are presented in section V, where we investigated the lowest lying excited states of the \ceUO2^2+ molecule as a test case and the two lowest-lying excited states of all-trans polyenes. Finally, we conclude in section VI.

II Targeting excited states in pCCD

In order to model excited states in pCCD, we will use the equation of motion (EOM) formalism Stanton and Bartlett (1993); Jankowski et al. (1994); Bartlett and Musiał (2007). In EOM-CC, excited states are parametrized by a linear CI-type ansatz Helgaker et al. (2000),

[TABLE]

where the summation is over all excitation operators present in the cluster operator $\hat{T}$ as well as the identity operator $\hat{\tau}_{0}$. The operator $\hat{R}$ is then used to generate the target state from the initial CC state,

[TABLE]

with $|0\rangle$ being the CC reference determinant.

To arrive at the EOM-CC working equations, it is convenient to use the normal-product form of the Hamiltonian, $\hat{H}_{N}=\hat{H}-\langle 0|\hat{H}|0\rangle$, written as a sum of the Fock operator $\hat{F}_{N}$ and the electron repulsion term $\hat{W}_{N}$ (in physicists’ notation),

[TABLE]

Furthermore, we will disregard any excitation properties, like dipole moments, and focus on excitation energies instead. In that case, we have to solve for the $\hat{R}$ amplitudes only. Our target-state Schrödinger equation then reads

[TABLE]

where $\Delta E$ is the energy difference with respect to the Fermi vacuum expectation value $|0\rangle$. Introducing the similarity transformed Hamiltonian in normal-product form $\hat{\mathcal{H}}_{N}=\exp{(-\hat{T})}\hat{H}_{N}\exp{(\hat{T})}$ and subtracting the equation for the CC ground state, we obtain the EOM-CC equations for the $\hat{R}$ amplitudes,

[TABLE]

where $\omega$ are the excitation energies with respect to the CC ground state, $\exp(\hat{T})|0\rangle$. The excitation energies are thus the eigenvalues of a non-Hermitian matrix,

[TABLE]

where the first row is associated with the CC reference state and the subsequent rows correspond to the excited configurations $\nu>0$. The EOM-CC working equations may be solved using, for instance, non-Hermitian extensions of the Davidson algorithm to determine the lowest-lying excited electronic states.

II.1 Electron-pair excitations

In this work, we are considering an pCCD reference function $|\mathrm{pCCD}\rangle$ as a special CC state confined to electron-pair states. As indicated in eq. (1), the pCCD cluster operator contains only electron-pair excitations, $\hat{T}=\hat{T}_{p}=\sum_{ia}t_{i}^{a}a^{\dagger}_{a}a^{\dagger}_{\bar{a}}a_{\bar{i}}a_{i}$. In the corresponding EOM model, the $\hat{R}$ operator is thus restricted to the identity operator $\hat{\tau}_{0}$ as well as all pair excitations present in the cluster operator $\hat{T}_{p}$,

[TABLE]

where $\hat{\tau}_{a\bar{a}i\bar{i}}=a^{\dagger}_{a}a^{\dagger}_{\bar{a}}a_{\bar{i}}a_{i}$ creates an electron pair in the virtual orbital $a$. To obtain the target-state Schrödinger equation of EOM-pCCD, we have to substitute the general cluster operator $\hat{T}$ by the pair-excitation operator $\hat{T}_{p}$ in eq. (5),

[TABLE]

The $\hat{R}_{p}$ amplitudes are determined from the EOM-pCCD equations (restricting $\hat{R}$ to $\hat{R}_{p}$ and $\hat{T}$ to $\hat{T}_{p}$ in eq. (6)),

[TABLE]

where $\mathcal{\hat{H}}^{(p)}_{N}$ indicates the similarity transformed Hamiltonian of pCCD, $\mathcal{\hat{H}}^{(p)}_{N}=\exp{(-\hat{T}_{P})}\hat{H}_{N}\exp{(\hat{T}_{P})}$, and $\omega_{p}$ are the electron-pair excitation energies. Thus, EOM-pCCD allows us to model electron-pair excited states only. In order to target singly excited or general doubly excited states, we have to extend the pCCD cluster operator to include excitations beyond electron pairs. This can be done either by changing to a frozen-pair CCSD Henderson et al. (2014b) formalism or to a linearized CCSD correction with an pCCD reference function Boguslawski and Ayers (2015). In this work however, we will consider a different, cost-effective approach to account for single excitations in the pCCD model that does not scale as $\mathcal{O}(n^{6})$ as conventional EOM-CCSD methods. Note that EOM-pCCD scales as $\mathcal{O}(o^{2}v^{2})$, where $o$ is the number of occupied orbitals (equivalent to the number of electron pairs) and $v$ is the number of virtual orbitals.

II.2 Accounting for single excitations

As proposed by Forseman et al. Foresman et al. (1992), configuration interaction with only single substitutions (CIS) represents an accurate model to investigate (singly) excited electronic states, even for large systems. In the CIS method, the reference is a single Slater determinant obtained from an SCF procedure, while the CIS wavefunction is expanded as

[TABLE]

with $|_{i}^{a}\rangle$ being a singly excited determinant where the (occupied) orbital $i$ of $|0\rangle$ has been substituted by the (virtual) orbital $a$. Similar to CIS, we will include single excitations in EOM-pCCD by extending the $\hat{R}_{p}$ operator of eq. (8). In addition to the identity operator $\hat{\tau}_{0}$ and all pair excitations $\hat{T}_{p}$, $\hat{R}_{p}$ also contains a summation over all single excitations,

[TABLE]

where $\hat{\tau}_{ai}$ is a singlet excitation operator $\hat{\tau}_{ai}=a^{\dagger}_{a}a_{i}+a^{\dagger}_{\bar{a}}a_{\bar{i}}$ that creates a singly excited electronic state with respect to $|0\rangle$, $|_{i}^{a}\rangle=\hat{\tau}_{ia}|0\rangle$. The $\hat{R}_{ps}$ amplitudes are determined from solving

[TABLE]

where we still have the similarity transformed Hamiltonian of pCCD, $\mathcal{\hat{H}}^{(p)}_{N}$, while $\omega_{ps}$ are the excitation energies of both singly excited and pair excited states. We will label this simplified model as EOM-pCCD+S to indicate that single excitations are included a posteriori in the $\hat{R}_{p}$ operator. In contrast to CIS that uses a single Slater determinant as reference, EOM-pCCD+S employs the pCCD wavefunction as reference state. Furthermore, electron correlation effects are included through the $\hat{T}_{p}$ operator in the similarity transformed Hamiltonian.

Similar to EOM-pCCD, the excitation energies are obtained by diagonalizing a non-Hermitian matrix of the form

[TABLE]

Note that in contrast to conventional EOM-CC methods, the first column does not equal zero because single excitations are not included in the cluster operator of pCCD. Thus, terms like $\langle_{i}^{a}|\hat{\mathcal{H}}_{N}^{(p)}|0\rangle$ do not vanish as they are not incorporated in the ground-state CC amplitude equations. Although we can account for single excitations in a rather straightforward way, we loose size-intensivity in the EOM model. For the molecular systems investigated in this work, the error introduced by extending only the $\hat{R}_{p}$ operator is approximately three orders of magnitude smaller than the actual excitation energies, while the computational cost increases insignificantly compared to EOM-pCCD (the Hamiltonian still contains terms that scale as $\mathcal{O}(o^{2}v^{2})$, but with a larger pre-factor). Thus, EOM-pCCD+S represents a cost-effective starting point to study singly excited electronic states in the pCCD model.

III Computational details

III.1 Molecular structures

In this work, we investigate all-trans polyenes containing 2 to 7 $\pi$-bonds (the corresponding Lewis structures are shown in Figure 1). Specifically, we study five different geometries, three for the lowest lying vertically excited states and two for lowest lying adiabatically excited states. The first set of structures was optimized for each ground state using the ADF software package adf (2014); te Velde et al. (2001); Guerra et al. (1998) with a DZ basis set van Lenthe and Baerends (2003) and the B3LYP exchange–correlation functional Stephens et al. (1994). The remaining sets of geometries were taken from ref. Hu and Chan (2015) and were optimized using the density matrix renormalization group (DMRG). Those include molecular geometries for the vertically and adiabatically excited states of \ceC10H12, \ceC12H14, and \ceC14H16. Furthermore, for both the vertically and adiabatically excited states, we study geometries optimized using two different active spaces as discussed in ref. Hu and Chan (2015): an active space containing only $\pi$-electrons (abbreviates as $v_{\pi}$ and $a_{\pi}$) and an active space comprising both $\sigma$- and $\pi$-electrons (denoted by $v$ and $a$, respectively).

For the linear \ceUO_2^2+ molecule, we calculated the four lowest-lying vertically and adiabatically excited states. The position of the minimum of each potential energy surface was obtained by simultaneously stretching the U–O bond between 1.610 and 1.870 Å.

III.2 Calculation of ground and excited states

The pCCD, EOM-pCCD, and EOM-pCCD+S calculations were performed with our in-house quantum-chemistry code employing restricted canonical Hartree–Fock orbitals as molecular orbital basis. The lowest-lying excitation energies were determined employing a modified Davidson algorithm. We should emphasize that the molecular orbital basis was not optimized within the pCCD method. Orbital optimization in pCCD typically results in localized, symmetry-broken orbitals that preclude us from identifying the symmetry of excited states. Furthermore, orbital optimization is only performed for the pCCD ground state wavefunction and results in molecular orbitals that are biased toward the ground state. This bias can be reduced by introducing a state-average protocol yielding a set of compromise orbitals (as usually done in multi-configuration SCF methods). Since we need to identify the symmetry of the (lowest-lying) excited states, we chose to not optimize the orbital basis and to use orbitals that transform according to the irreducible representations of the molecular point group (numerical examples for excitation energies within pCCD obtained in the optimized orbital basis can be found in the Supporting Information). Moreover, recent numerical studies confirm that pCCD using restricted Hartree–Fock orbitals captures, at least qualitatively, the most important correlation effects around the equilibrium geometry Boguslawski et al. (2016).

For the investigated all-trans polyenes, we used the 6-31G basis set as well as the cc-pVDZ basis set of Dunning Dunning (1989) for direct comparison with either CIS(D) (6-31G) or MRMP/DMRG (cc-pVDZ) reference data, respectively.

For \ceUO_2^2+, scalar relativistic effects were incorporated through relativistic effective core potentials (RECP). In all calculations, we have used a small core (SC) RECP (60 electrons in the core) with the following contraction scheme (12s11p10d8f) $\rightarrow$ [8s7p6d4f] Küchle et al. (1991). For the lighter elements (O), the cc-pVDZ basis set of Dunning Dunning (1989) was employed, (10s5p1d) $\rightarrow$ [4s3p1d]. The CIS and EOM-CCSD calculations have been performed with the Molpro software package Werner et al. (2012). No frozen core was used in all calculations.

IV The UO${}_{2}^{2+}$ molecule as a test case

Before scrutinizing the performance of EOM-pCCD and EOM-pCCD+S in modeling the lowest-lying excited states in all-trans polyenes, we will assess the accuracy of EOM-pCCD+S against molecular systems whose lowest-lying excited states are purely singly-excited states. For that purpose, we chose the uranyl cation (\ceUO_2^2+) as it contains a large number of strongly-correlated electrons distributed among the 5$f$-, 6$d$-, and 7$s$-orbitals and its electronic spectrum is well understood Zhang and Pitzer (1999); Pierloot and van Besien (2005); Denning (2007); Réal et al. (2007); Tecmer et al. (2011, 2012b).

The vertical and adiabatic excitation energies for the four lowest-lying excited states in EOM-pCCD+S and CIS are summarized in Table 1. Note that for the \ceUO_2^2+ molecule the two lowest-lying excited states can be accurately described within the EOM-CCSD model, which yields excitations energies that are similar to completely renormalized EOM-CCSD(T) reference values Tecmer et al. (2013). For the higher-lying excited states, however, dynamic correlation effects become important and a triples correction has to be included to accurately model those states. The EOM-CCSD results in Table 1 can thus be considered as upper bounds of the excitation energies for the $\pi_{u}\rightarrow\delta_{u}$ state.

In general, EOM-pCCD+S overestimates vertical excitations energies of the two lowest-lying excited states ($\sigma_{u}\rightarrow\phi_{u}$ and $\sigma_{u}\rightarrow\delta_{u}$) by approximately 0.5 eV, while the corresponding adiabatic excitation energies are overestimated by about 0.7 eV. Note that for the uranyl cation, CIS outperforms EOM-pCCD+S and deviates from 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 to accurately predict the excitation energies for the $\pi_{u}\rightarrow\delta_{u}$ states as these excited states are dominated by dynamic correlation effects which are not included in CIS and only marginally accounted for in pCCD (see also Boguslawski et al. (2016)). We should emphasize that EOM-pCCD+S provides equilibrium U–O bond lengths of excited states that deviate less from the EOM-CCSD reference values ($\Delta r_{e}\approx 0.06$ Å compared to $\Delta r_{e}\approx 0.08$ Å in CIS).

To conclude, the EOM-pCCD+S model predicts the correct order of the lowest-lying excited states in the uranyl cation and provides excitation energies of decent accuracy with errors of about 0.5 eV with respect to EOM-CCSD reference values. However, errors in excitation energies are slightly worse than in the simple CIS model.

V Excited states in all-trans polyenes

All-trans polyenes are model systems for carotenoids and polyene chromophores that play an important role in photoprocesses. The proper description of the two lowest-lying excited states poses a challenge to both experiment and quantum chemistry approaches Hudson and Kohler (1972); Schulten et al. (1976); Cave and Davidson (1987); Tavan and Schulten (1987); Cave and Davidson (1988); Buma et al. (1991); Orlandi et al. (1991); Serrano-Andres et al. (1993a, b); Krawczyk et al. (2000); Hsu et al. (2001); Wiberg et al. (2002); Starcke et al. (2006); Christensen et al. (2008), especially because doubly excited configuration are required to accurately model ground and excited states of longer polyenes. The excitation energies for the two-lowest vertical excited states of \ceC4H6 to \ceC14H16 using the DFT optimized structures and two different basis sets (cc-pVDZ and 6-31G) are summarized in Table 2. The excitation energies obtained in the 6-31G basis are given in parenthesis next to the corresponding results calculated for the cc-pVDZ basis set. Note that the 6-31G basis set was used in CIS(D) Starcke et al. (2006), while a cc-pVDZ basis set was employed in MRMP calculations  Kurashige et al. (2004). Since EOM-pCCD can only model pair excited states, only the excitation energies of the $2^{1}A_{g}^{-}$ state are shown in the Table.

For all investigated polyenes, EOM-pCCD overestimates the excitation energies of the first excited state $2^{1}A_{g}^{-}$ by approximately 4 eV compared to MRMP and 1.5 eV compared to CIS(D) data. This error can be partially attributed to missing single excitations that have to be included in the model to describe the character of the $2^{1}A_{g}^{-}$ state accurately. Including singles on top of EOM-pCCD a posteriori reduces the error by approximately 2 to 1 eV, depending on the number of $\pi$ bonds. While EOM-pCCD+S predicts excitation energies that are lower than the corresponding excitation energies of CIS(D) (differences amount up to 1.7 eV), it overestimates the excitation energies of the $2^{1}A_{g}^{-}$ state by about 2 eV compared to MRMP data. Furthermore, the simple EOM-pCCD+S model does not predict the right order of states where the first dark state $2^{1}A_{g}^{-}$ should lie below the first bright state $1^{1}B_{u}^{+}$. The excitation energies of the first bright state $1^{1}B_{u}^{+}$, however, deviate only by approximately 0.3 eV from MRMP data. We should note that EOM-pCCD+S outperforms CIS(D) in predicting the excitation energies of the $1^{1}B_{u}^{+}$ state. Despite its simplicity, EOM-pCCD+S is thus a good starting point to model singly excited states.

We can extrapolate the excitation energies of the two lowest-lying excited states for longer polyenes in the limit of infinite number of $\pi$ bonds and predict the ordering of states for longer polyenes (see, for instance, refs. Schulten et al. (1976); Tavan and Schulten (1987); Kurashige et al. (2004); Starcke et al. (2006)). To do so, the excitation energies can be fitted to a particle-in-a-box model Starcke et al. (2006) of the form

[TABLE]

where $L(n_{\pi})$ is the length of the polyene that depends on the number of $\pi$ bonds $n_{\pi}$, while $a,b,c$ are the fitting parameters. The function for the length of the polyene is obtained by linear regression of the terminal C–C distance of the DFT optimized structures and defines the distance between the terminal carbon atoms as a function of the number of double bonds,

[TABLE]

The excitation energy at infinite number of $\pi$ bonds is obtained from the fitting parameter $c$. For EOM-pCCD, the extrapolated excitation energy of the $2^{1}A_{g}^{-}$ state is 5.61 eV and reduces to 4.28 eV for EOM-pCCD+S. For the $1^{1}B_{u}^{+}$ state, EOM-pCCD+S results in an extrapolated excitation energy of 3.06 eV. The resulting fits as well as excitation energies are shown in Figure 2. Note that EOM-pCCD+S model overestimates the extrapolated excitation energy of the $1^{1}B_{u}^{+}$ state by approximately 1.30 eV compared to experiment Christensen et al. (2008) and by about 0.94 eV compared to MRMP reference data (2.12 eV) Kurashige et al. (2004).

Finally, we will discuss how the molecular structure affects excitation energies and their character in the EOM-pCCD+S method. Note, however, that the labels “$v$” and “$a$” in Table 3 correspond to the vertically and adiabatically excited state in DMRG calculations. They do not represent the true adiabatic and vertical excitation energies in EOM-pCCD+S as the DMRG-optimized structures are used. To simplify the labeling of the excited states and to facilitate direct comparison to DMRG, we will use the same labeling (“$v$” and “$a$”) in DMRG and EOM-pCCD+S for all investigated excited states. Figure 3 and Table 3 show the two lowest-lying excitation energies of \ceC10H12, \ceC12H14, and \ceC14H16 for different molecular geometries optimized by DMRG. For the $2^{1}A_{g}^{-}$ state, the excitation energies predicted by EOM-pCCD+S reduce by approximately 1 eV if the molecular structure is allowed to relaxed (adiabatic excitations). The choice of the active space in the DMRG structure optimization, however, does not significantly affect the accuracy of EOM-pCCD+S for longer polyene chain lengths (differences are less than 0.1 eV between $v$ ($a$) and $v_{\pi}$ ($a_{\pi}$)), while DMRG is more sensitive to the active space (differences in DMRG excitation energies amount to 1 eV). Similar observation can be made for the $1^{1}B_{u}^{+}$ state. For relaxed molecular structures, the excitation energies slightly improve by 0.2 eV compared to the vertical excitation energies. Furthermore, the active space used in DMRG calculations does not change the excitation energies considerably and similar excitation energies are obtained for the vertically and adiabatically excited states, respectively, with EOM-pCCD+S. We should emphasize that the excitation energies of the $1^{1}B_{u}^{+}$ state predicted by EOM-pCCD+S are closer to experimental values than the corresponding DMRG excitation energies. Differences between the theoretically determined values and experiment can be partly attributed to the basis set used in calculations as well as missing electron correlation effects not included in the theoretical model.

Figure 4 shows the contributions (in percent) of the three dominant configurations of the $2^{1}A_{g}^{-}$ excited state for different molecular geometries. These are the doubly excited $b_{g}\rightarrow a_{u}^{*}$ (HOMO${}^{2}\rightarrow$LUMO2), the singly excited $b_{g}\rightarrow b_{g}^{*}$ (HOMO$\rightarrow$LUMO+1), and the singly excited $a_{u}\rightarrow a_{u}^{*}$ (HOMO$-$1$\rightarrow$LUMO) Schulten et al. (1976); Tavan and Schulten (1987). For the DFT optimized structures shown in Figure 4(a) (vertical excitations), the contribution from the doubly-excited configuration HOMO2 $\rightarrow$ LUMO2 increases to approximately 50-60%, while both singly excited configurations constitute 15 to 25%. For longer polyene chain lenghts, the contribution from the doubly-excited configuration decreases.

The character of the $2^{1}A_{g}^{-}$ state changes considerably if the molecular structure is modified. Specifically, the weights of the configurations of the vertically excited $2^{1}A_{g}^{-}$ state strongly depend on the number of $\pi$-bonds if DMRG optimized structures are used (see solid lines in Figure 4(b)). While the dark state in \ceC10H12 is dominated by the HOMO${}^{2}\rightarrow$LUMO2 configuration for the shorter polyene chain lengths, the weight of the HOMO$-$1$\rightarrow$LUMO configuration gradually increases and constitutes more than 60 % for \ceC14H16. Extending the polyene chain to 7 double bonds further decreases the weight of the HOMO${}^{2}\rightarrow$LUMO2 configuration to less than 10 %, while the contribution of the singly excited HOMO$\rightarrow$LUMO+1 configuration slightly increases. Thus the vertically excited $2^{1}A_{g}^{-}$ state of \ceC14H16 is dominated by only one configuration.

The picture changes completely if the molecular structure of the excited state is allowed to relax (see dashed lines in Figure 4(b)). For growing polyene chains, the weight of the doubly excited HOMO${}^{2}\rightarrow$LUMO2 configuration gradually increases (up to 70 % for \ceC14H16), while both singly excited configurations (HOMO$\rightarrow$LUMO+1) and HOMO$-$1$\rightarrow$LUMO) approach weights of approximately 15 %. Although EOM-pCCD+S predicts a dominant contribution of the doubly-excited HOMO${}^{2}\rightarrow$LUMO2 configuration in the $2^{1}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 % for longer polyene chain lenghts Starcke et al. (2006). This underestimation can be attributed to broken-pair configurations that are not included in the pCCD model and have to be accounted for a posteriori Henderson et al. (2014b); Boguslawski and Ayers (2015). To conclude, EOM-pCCD+S can appropriately describe the three most important configurations present in the dark state (see also  Schulten et al. (1976); Tavan and Schulten (1987); Starcke et al. (2006)) if molecular structures are relaxed. The qualitatively good description of the character of the first dark state is also evident in the excitation energies. The adiabatic excitation energies are approximately 1 to 1.5 eV lower in energy than the corresponding vertical ones (compare Tables 2 and 3).

VI Conclusions

In this work, we have extended the pCCD model to excited states using the equation of motion formalism. Since the cluster operator of pCCD is restricted to electron-pair states, only electron-pair excitations can be 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 single and general double excitations as well which considerably increases the computational cost compared to the pCCD reference calculation. To arrive at a cost-effective model similar to pCCD, we have included single excitation a posteriori in the CI-type ansatz. Our method (EOM-pCCD+S) is motivated by the CIS approach that allows us to efficiently model singly excited states with an SCF reference function. Since only the CI-type ansatz is modified, the computational scaling remains similar to EOM-pCCD (the Hamiltonian contains terms that scale as $\mathcal{O}(o^{2}v^{2})$), but with a larger pre-factor. Despite the simplicity of the model, EOM-pCCD+S breaks size-intensivity, which can be restored by adjusting the cluster operator.

Although pCCD is commonly combined with an orbital-optimization protocol, we did not optimize the molecular orbital basis in the pCCD reference calculation. Orbital optimization in pCCD usually results in localized, symmetry-broken orbitals that prevent us from identifying the symmetry of the targeted excited states. Furthermore, if the orbitals are optimized within pCCD, the optimal set of orbitals is biased toward the ground state and the corresponding excitation energies are generally overestimated. This bias can be reduced by optimizing the orbitals simultaneously for both ground and excited states, which increases the computational cost of EOM-pCCD and EOM-pCCD+S significantly. To benefit from the $\mathcal{O}(o^{2}v^{2})$ scaling of both excitation models, we have thus skipped the orbital optimization step and used restricted Hartree–Fock orbitals in all calculations.

Both excitation models have been assessed against excitation energies of the uranyl cation, whose lowest-lying excited states contain purely singly-excited states, and all-trans polyenes containing 2 to 7 $\pi$-bonds. For the \ceUO2^2+ molecule, EOM-pCCD+S provides excitation energies that are close to CIS results, overestimates, 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 of the first dark and first bright state, while the excitation energies of the $1^{1}B_{u}^{+}$ state are in good agreement with MRMP results and closer to experiment than DMRG reference data. Furthermore, the excitation energies and the character of the excited states strongly depend on the molecular structures used in calculations. Specifically, the character of the $2^{1}A_{g}^{-}$ state can only be properly predicted if molecular structures are allowed to relax, resulting in a dominant doubly-excited HOMO${}^{2}\rightarrow$LUMO2 configuration and two singly-excited configurations (HOMO$\rightarrow$LUMO+1) and HOMO$-$1$\rightarrow$LUMO). In order to accurately model the first dark state in all-trans polyenes, double excitations beyond electron-pair excitations as well as higher excitations (triples, etc.) might be important to capture the missing correlation effects in the targeted excited states that cannot be described within pCCD. Possible extension of EOM-pCCD using CC-type corrections on top of the pCCD reference function are presently being developed in our laboratory. To conclude, EOM-pCCD+S represents a good and cost-effective starting point to investigate singly-excited states in the pCCD model. Additional numerical studies are currently under investigation in our laboratory.

Supplementary Material

See Supplementary Material for excitation energies of selected all-trans polyenes obtained by EOM-pCCD and EOM-pCCD+S in the optimized pCCD orbital basis.

Acknowledgements

Financial support from a SONATA BIS grant of the National Science Centre, Poland (no. 2015/18/E/ST4/00584) and a Marie-Skłodowska-Curie Individual Fellowship project no. 702635–PCCDX is gratefully acknowledged.

Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), grant No. 412.

We had many helpful discussions with Leszek Meissner, Piotr Piecuch, and Karol Kowalski.