Targeting doubly-excited states with Equation of Motion Coupled Cluster theory restricted to double excitations

TL;DR

This paper introduces a new EOM coupled cluster method tailored for doubly-excited states, combining pCCD and LCCSD to improve accuracy in challenging electronic excitations.

Contribution

The authors develop EOM-pCCD-LCCSD, a simplified yet accurate approach for describing doubly-excited states, outperforming traditional methods in specific molecular systems.

Findings

Accurately predicts excitation energies in polyenes.

Outperforms conventional EOM-CCSD in state ordering.

Closest to experimental results among tested methods.

Abstract

The accurate description of doubly-excited states using conventional electronic structure methods is remarkably challenging, primarily because such excited states require the inclusion of doubly or higher excited configurations or the application of multi-reference methods. We present a new approach to target electronically excited states that feature a double-electron transfer. Our method uses the equation of motion (EOM) formalism with a pair coupled cluster doubles (pCCD) reference function, where dynamical correlation is accounted for by a linearized coupled cluster correction with singles and doubles (LCCSD). Specifically, our proposed EOM-pCCD-LCCSD model represents a simplification of the conventional EOM-CCSD approach, where the electron-pair amplitudes of CCSD are tailored by pCCD. The performance of EOM-pCCD-LCCSD is assessed for the lowest-lying excited states in CH$^+$ and all-trans polyenes. In contrast to conventional EOM-CC methods with at most double excitations, EOM-pCCD-LCCSD predicts the right order of states in polyenes with excitation energies closest to experiment, outperforming even highly-accurate methods such as the density matrix renormalization group algorithm.

Figures

Tables

Table 1. Table 1: Error measures [eV] for the ground and lowest-lying excited states of \ce CH+ determined for pCCD-LCCSD (with and without orbital optimization) and EOM-CCSD with respect to FCI.  20 The error measures are calculated for the total electronic energy E el subscript 𝐸 el E_{\rm el} of each state and the corresponding excitation energies ω el subscript 𝜔 el \omega_{\rm el} . NPE: non-parallelity error; MAE: maximum absolute error; ME: mean error; RME: root mean error.

		EOM-pCCD-LCCSD(HF)				EOM-pCCD-LCCSD(opt)				EOM-CCSD 21
	State	NPE	MAE	ME	RME	NPE	MAE	ME	RME	NPE	MAE	ME	RME
$E_{\mathrm{el}}$	X ${}_{}{}^1\Sigma _{}^{+}$	0.811	1.019	0.394	0.453	0.113	0.139	0.110	0.115	0.153	0.205	0.099	0.109
	$2^1{\Sigma }^{+}$	0.676	0.884	-0.573	0.620	0.234	0.510	0.442	0.446	0.507	0.751	0.607	0.623
	$3^1{\Sigma }^{+}$	0.244	0.266	0.010	0.189	0.282	0.447	0.221	0.232	0.775	0.883	0.334	0.422
	$1^1\Pi$	0.372	0.530	0.298	0.321	0.396	0.540	0.275	0.301	0.420	0.498	0.226	0.261
	$2^1\Pi$	0.317	0.716	0.637	0.645	0.353	0.735	0.638	0.648	0.322	0.650	0.566	0.575
	$1^1\Delta$	0.406	0.707	0.481	0.506	0.572	1.426	1.114	1.126	0.459	1.408	1.177	1.186
${\omega }_{\mathrm{el}}$	$2^1{\Sigma }^{+}$	0.411	1.229	-0.967	0.974	0.124	0.373	0.332	0.334	0.628	0.668	0.508	0.539
	$3^1{\Sigma }^{+}$	1.282	1.282	-0.383	0.555	0.382	0.421	0.110	0.150	0.622	0.678	0.235	0.317
	$1^1\Pi$	0.460	0.489	-0.096	0.154	0.499	0.514	0.164	0.225	0.267	0.293	0.127	0.153
	$2^1\Pi$	0.400	0.431	0.243	0.317	0.456	0.710	0.528	0.544	0.288	0.563	0.467	0.476
	$1^1\Delta$	0.277	0.312	0.088	0.151	0.675	1.400	1.003	1.022	0.306	1.203	1.079	1.083

Table 2. Table 2: Vertical excitation energies [eV] of the two lowest-lying excited states in all-trans polyenes \ce C4H6 to \ce C14H16 calculated for various flavours of EOM-pCCD 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 and all EOM-CC calculations. CASSCF was performed in a double-zeta basis set (see corresponding references).

\ceC=C	2	3	4	5	6	7
	$2^1{A}_g^{-}$
EOM-pCCD 15, 16	10.56	9.11	8.11	7.42	6.93	6.58
EOM-pCCD+S 15, 16	7.45	6.79	6.15	5.69	5.37	5.13
EOM-pCCD-LCCSD	6.26	5.29	4.58	4.11	3.79	3.57
EOM-LCCSD	7.91	6.98	6.28	5.70	5.09	–
EOM-CCSD	7.59	6.64	5.94	5.43	5.05	4.75
CIS(D) 35	9.01	7.81	6.78	6.12	5.55	5.14
MRMP 37	6.31	5.10	4.26	3.68	3.19	2.80
CASSCF 38	6.67	5.64	5.16	4.32	–	–
	$1^1{B}_u^{+}$
EOM-pCCD+S 15, 16	7.20	5.98	5.19	4.62	4.20	3.87
EOM-pCCD-LCCSD	7.00	5.91	5.23	4.77	4.43	4.17
EOM-LCCSD	7.24	6.27	5.70	5.37	5.28	–
EOM-CCSD	6.86	5.74	5.03	4.53	4.16	3.87
CIS(D) 35	8.09	6.78	5.95	5.43	5.00	4.70
MRMP 37	6.21	5.25	4.57	4.17	3.87	3.60
CASSCF 38	7.73	7.06	6.62	6.37	–	–

Table 3. Table 3: Vertical and adiabatic excitation energies [eV] of the two lowest-lying excited states in all-trans polyenes \ce C10H12 to \ce C14H16 calculated for various flavours of EOM-CC and DMRG. Note that all orbitals are active in all performed EOM-CC calculations.

	Excited state	$2^1{A}_g^{-}$			$1^1{B}_u^{+}$
	\ceC=C	5	6	7	5	6	7
v	EOM-pCCD+S 15, 16	6.28	5.99	5.76	4.91	4.50	4.25
	EOM-pCCD-LCCSD	4.67	4.40	4.32	4.98	4.65	4.44
	EOM-LCCSD	6.11	5.68	5.39	5.42	5.15	4.98
	EOM-CCSD	5.79	5.43	5.22	4.78	4.41	4.19
	DMRG39	5.43	4.76	4.64	5.35	4.98	4.66
a	EOM-pCCD+S 15, 16	5.18	4.61	4.44	5.20	4.30	4.04
	EOM-pCCD-LCCSD	3.47	2.77	1.85	5.08	4.38	4.15
	EOM-LCCSD	4.81	6.47	6.02	5.11	4.98	4.51
	EOM-CCSD	5.16	4.81	4.59	4.85	4.14	3.90
	DMRG39	4.01	3.41	3.22	4.98	4.60	4.29
Exp.40		3.03	2.69	2.44	3.57	3.31	3.12

Equations

$$|\mathrm{pCCD}⟩=\exp \left(\sum _{i=1}^{\mathrm{occ}}\sum _{a=1}^{\mathrm{virt}}{t}_i^a{a}_a^{†}{a}_{\bar{a}}^{†}{a}_{\bar{i}}{a}_i\right)|{\Phi }_0⟩=e^{{\hat{T}}_{\mathrm{p}}}|{\Phi }_0⟩,$$

$$|\Psi ⟩=\exp (\hat{T})|\mathrm{pCCD}⟩.$$

$$⟨{\Phi }_{\nu }|({\hat{H}}_N+[{\hat{H}}_N,{\hat{T}}_{\mathrm{np}}])|\mathrm{pCCD}⟩=0,$$

$$⟨{\Phi }_{\nu }|({\hat{H}}_N+[{\hat{H}}_N,{\hat{T}}_{\mathrm{np}}]+[{\hat{H}}_N,{\hat{T}}_{\mathrm{p}}]+[[{\hat{H}}_N,{\hat{T}}_{\mathrm{np}}],{\hat{T}}_{\mathrm{p}}])|{\Phi }_0⟩=0.$$

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

$$|\Psi ⟩=\hat{R}\exp (\hat{T})|{\Phi }_0⟩=\sum _{\mu }{c}_{\mu }{\hat{\tau }}_{\mu }\exp (\hat{T})|{\Phi }_0⟩.$$

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

Full text

Institute of Physics, Nicolaus Copernicus University in Torun] Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University in Torun, Grudziadzka 5, 87-100 Torun, Poland \alsoaffiliation[Faculty of Chemistry, Nicolaus Copernicus University in Torun] Faculty of Chemistry, Nicolaus Copernicus University in Torun, Gagarina 7, 87-100 Torun, Poland \phone(+48)56 611 3338

\abbreviationsAP1roG, LCCSD, EOM, pCCD

1 Introduction

The efficient and reliable description of electronically excited states of atoms and molecules is gaining in importance in many areas of chemistry, physics, biology, and materials science. This trend promotes the development of new quantum chemistry methods dedicated to specifically model electronically excited states properties in large molecules and complex systems 1, 2, 3, 4, 5, 6, 7. Such methods are, however, usually designed to accurately treat electronic excitations with a dominant transfer of one single electron. Excitation energies featuring a double-electron transfer are more problematic as their description requires the inclusion of higher order excitations or even a multi-reference treatment 8. Both methodologies typically suffer from an unfavorable computational scaling, which limits their application to small model systems. To advance a reliable, but inexpensive description of doubly-excited states in large molecules, we have developed a new approach based on the computationally tractable pair coupled cluster doubles (pCCD) model 9, 10. Specifically, the pCCD approach represents a simplified version of CCD, where only electron-pair excitations are kept in the cluster operator,

[TABLE]

where $a^{\dagger}_{p}$ and $a_{p}$ ($a^{\dagger}_{\bar{p}}$ and $a_{\bar{p}}$) are the electron creation and annihilation operators for $\alpha$ ($\beta$) electrons and $|\Phi_{0}\rangle$ is some reference determinant. By construction, pCCD can only describe correlations restricted to electron pairs. The missing correlation effects that can be attributed to broken-pair states have to be included a posteriori. This can be achieved, for instance, using perturbation theory 11, 12 or CC corrections 13, 14. In the latter approach, the electronic wave function is written using an exponential ansatz with some cluster operator $\hat{T}=\sum_{\nu}t_{\nu}\tau_{\nu}$ ($\hat{\tau}$ is an excitation operator) and the pCCD wave function as reference state,

[TABLE]

$\hat{T}$ may contain single excitations $\hat{T}_{1}$, (non-pair) double excitations $\hat{T}_{2}^{\prime}$, and higher excitations 13. Note that the $\prime$ indicates the exclusion of pair excitations, that is, $\hat{T}_{2}=\hat{T}_{2}^{\prime}+\hat{T}_{\rm p}$. To facilitate the distinction between electron-pair (p) excitations and non-pair (np) excitations, we will label the excitation operators accordingly. The CC correction can be further simplified by considering a linearized CC (LCC) ansatz. 14 The cluster amplitudes $t_{\nu}$ are then determined by solving a linear set of coupled equations

[TABLE]

where the Baker–Campbell–Hausdorff expansion $e^{-\hat{T}_{\rm np}}\hat{H}_{N}e^{\hat{T}_{\rm np}}=\hat{H}_{N}+[\hat{H}_{N},\hat{T}_{\rm np}]+\frac{1}{2}[[\hat{H}_{N},\hat{T}_{\rm np}],\hat{T}_{\rm np}]+\ldots$ has been truncated after the second term, $\hat{H}_{N}=\hat{H}-\langle\Phi_{0}|\hat{H}|\Phi_{0}\rangle$ is the quantum chemical Hamiltonian in its normal product form, and $|\Phi_{\nu}\rangle=\hat{\tau}_{\nu}|\Phi_{0}\rangle$. Using eq. (1), eq. (3) can be brought into the familiar form of single-reference CC theory,

[TABLE]

Thus, the hybrid pCCD-LCCSD (linearized coupled cluster singles and doubles) approach represents a simplification of CCSD and frozen-pair (fp)CCSD 13, where all non-pair amplitudes enter linearly into the CC amplitudes equation eq. (4). The $\hat{T}_{\rm p}$ operator introduces non-linear terms into the single and double amplitudes equations, which explicitly couple the pair amplitudes with all non-pair amplitudes. Most importantly, all proposed CC corrections on top of pCCD represent some tailored CC flavor, where the electron-pair amplitudes of the CC correction are taken from a previously optimized pCCD wave function.

Numerical studies suggest 14, 12 that pCCD-LCCSD is a reliable wave function ansatz to model both static and dynamic electron correlation and allows us to reach chemical accuracy ($\sim 1$ kcal/mol) for many challenging systems. 12 The pCCD-LCCSD method is, however, applicable to ground-state electronic structures only and must be extended to target excited states. Recently, we presented the first extension of pCCD to describe excited-state electronic structures, where single excitations are treated approximately in the excited state model. 15, 16 The proposed EOM-pCCD+S ansatz can reliably describe singly-excited states when dynamic electron correlation effects are similar in both the ground and excited states. If, however, correlation effects differ in the ground and excited states, the corresponding excitation energies will be inaccurate. To improve the EOM-pCCD+S model, we will present an excited state model based on the pCCD-LCCSD reference function. The proposed excited state ansatz will be assessed against the lowest-lying excited states of the \ceCH^+ molecule and all-trans polyenes with up to 7 double bonds.

2 Extending pCCD-LCCSD to model excited states

We will target excited states using the equation-of-motion (EOM) formalism, where excited states are modelled using a linear CI-type ansatz,

[TABLE]

The above sum runs over all excitations present in the cluster operator as well as the identity operator $\hat{\tau}_{0}$. The $\hat{R}$ operator, then, generates the excited state by acting on the CC reference state,

[TABLE]

Since we will focus on excitation energies, we have to solve for the $\hat{R}$ amplitudes only. 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,

[TABLE]

where $\omega=(E-E_{0})$ are the excitation energies with respect to the CC ground state, $\exp(\hat{T})|\Phi_{0}\rangle$. As the pCCD-LCCSD cluster operator contains all single and double excitations, the $\hat{R}$ operator becomes $\hat{R}=\hat{R}_{0}+\hat{R}_{1}+\hat{R}_{2}$. Furthermore, the similarity transformed Hamiltonian of pCCD-LCCSD has a special form 14. Specifically, for all non-pair excitations, we have the linearized approximation $\mathcal{\hat{H}}^{\rm np}_{N}=e^{-\hat{T}_{\rm p}-\hat{T}_{\rm np}}\hat{H}_{N}e^{\hat{T}_{\rm p}+\hat{T}_{\rm np}}\approx\hat{H}_{N}+[\hat{H}_{N},\hat{T}]+[[\hat{H},\hat{T}_{\rm np}],\hat{T}_{\rm p}]$, while the electron-pair amplitudes are optimized in pCCD with $\mathcal{\hat{H}}^{\rm p}_{N}=e^{-\hat{T}_{\rm p}}\hat{H}_{N}e^{\hat{T}_{\rm p}}=\hat{H}_{N}+[\hat{H}_{N},\hat{T}_{\rm p}]+\frac{1}{2}[[\hat{H},\hat{T}_{\rm p}],\hat{T}_{\rm p}]$. Although the similarity transformed Hamiltonian differs for electron-pair and non-pair excitations, we have considered the most general form of the $\hat{R}$ operator ($\hat{R}=\hat{R}_{0}+\hat{R}_{1}+\hat{R}_{2}$) when calculating the excitation energies from eq. (7). The diagrammatic representation of the left-hand-side of the EOM-pCCD-LCCSD equations eq. (7) are summarized in Fig. S1 of the Supporting Information.

Finally, we should note that the orbital-pairing scheme as well as orbital optimization are crucial in pCCD calculations. Specifically, size consistency can be recovered if the orbitals are optimized within pCCD. Numerical studies suggest 12 that dynamic energy corrections on top of pCCD yield similar equilibrium properties in canonical Hartree–Fock and pCCD optimized natural orbitals. Hence, orbital optimization can be omitted if (relative) energies or properties of equilibrium structures or molecular structures close to the equilibrium are sought. For stretched bonds or in the vicinity of dissociation, orbital optimization might become important to reach spectroscopic accuracy due to size-consistency errors. Since, however, the molecular orbitals are typically optimized for the ground state, the conventional pCCD orbital optimization protocol results in molecular orbitals that are biased toward the ground state. Furthermore, (variational) orbital optimization usually results in symmetry-broken (localized) orbitals 17, 10, 18, 19 that prevent us from identifying excited states according to their (point group) symmetry. To provide information about spatial symmetry, excited state calculations can either be performed within the canonical Hartree–Fock orbital basis or using orbitals that have been optimized imposing point-group symmetry. In the following, we will investigate both choices for excited state calculations with a pCCD-LCCSD reference function.

3 Targeting doubly-excited states in \ceCH+ and all-trans polyenes

Our first test system comprises the excited-state potential energy surfaces (PESs) of a benchmark system for which full configuration interaction (FCI) results are available. Specifically, we focus on the lowest-lying excited states of the \ceCH+ molecule, whose valence excited states contain significant biexcited components when the \ceC-H bond is stretched. Theoretical studies showed that the conventional EOM-CCSD approach yields large errors in excitation energies and breaks the degeneracy of states in the dissociation limit. 20, 21 Thus, the \ceCH+ molecule represents an ideal test system to assess the accuracy of the simplified EOM-pCCD-LCCSD formalism.

Fig. 1 illustrates the PESs for the ground and lowest-lying excited states determined for FCI, EOM-CCSD, and EOM-pCCD-LCCSD and the three lowest-lying ${}^{1}\Sigma^{+}$ states (including the ground state), the two lowest-lying ${}^{1}\Pi$ states, and the lowest-lying ${}^{1}\Delta$ state. Note that the first excited ${}^{1}\Sigma^{+}$ and ${}^{1}\Delta$ states feature a doubly-excited character, while the second ${}^{1}\Pi$ state has partial biexcited character. Hence, we will focus our discussion on these particular states. To facilitate our analysis, Table 1 summarizes some selected error measures in excitation and total electronic energies with respect to FCI results. These include the non-parallelity error (NPE $=\max\limits_{r_{\rm CH}}(|\Delta E_{r_{\rm CH}}|)-\min\limits_{r_{\rm CH}}(|\Delta E_{r_{\rm CH}}|)$ and $\Delta E_{r_{\rm CH}}=E^{\rm CC}_{r_{\rm CH}}-E^{\rm FCI}_{r_{\rm CH}}$), the maximum absolute error (MAE $=\max\limits_{r_{\rm CH}}(|\Delta E_{r_{\rm CH}}|)$), the mean error (ME $=\sum\limits_{r_{\rm CH}}\Delta E_{r_{\rm CH}}/N$), and the root mean error (RME $=\sqrt{\sum\limits_{r_{\rm CH}}\Delta E_{r_{\rm CH}}^{2}/N}$) determined for all C–H distances $r_{\rm CH}$ along each PES, where $N$ is the total number of points along the PES. Note that for pCCD-LCCSD, we performed two different calculations: (i) using the canonical Hartree–Fock (HF) basis and (ii) using the optimized pCCD natural orbitals (opt) for the ground state molecule imposing $C_{2h}$ symmetry (to label the excited states according to their symmetry). While all ${}^{1}\Sigma^{+}$ states are better described within the optimized pCCD basis, the ${}^{1}\Pi$ and ${}^{1}\Delta$ states have smaller error measures using canonical HF orbitals. Nonetheless, EOM-pCCD-LCCSD improves the description of all doubly-excited states and results in MAE of about 0.3 ($1^{1}\Delta$) to 0.4 ($2^{1}\Sigma$ and $2^{1}\Pi$) eV for excitation energies, while EOM-CCSD excitation energies yield an MAE of 1.2 ($1^{1}\Delta$) and about 0.6 ($2^{1}\Sigma$ and $2^{1}\Pi$) eV. Furthermore, the corresponding RMEs are reduced from 1.1 eV to 0.15 eV ($1^{1}\Delta$) and from 0.5 to 0.3 eV ($2^{1}\Sigma$ and $2^{1}\Pi$). Similar observation can be made for the total electronic energies. Note, however, that all error measures for $E_{\rm el}$ increase by 0.1–0.4 eV. In general, the excited-state PESs predicted by EOM-pCCD-LCCSD are qualitatively closer to the FCI reference curve and result in smaller NPEs (0.2 to 0.4 eV) compared to EOM-CCSD (0.3 to 0.78 eV). Most importantly, EOM-pCCD-LCCSD(HF) can—at least partly—numerically restore the asymptotic degeneracy of the $1^{1}\Delta$ and $2^{1}\Pi$ state, fails, however, for the $2^{1}\Sigma$ state (those three states dissociate into C(1D) and H+). Specifically, the energy difference between the $1^{1}\Delta$ and $2^{1}\Pi$ state at 5.0 bohr is 0.066 eV for EOM-pCCD-LCCSD(HF) and increases to 0.865 eV for EOM-CCSD, while the $2^{1}\Sigma$ state lies 0.917 (1.163) eV below the $1^{1}\Delta$ state for EOM-pCCD-LCCSD (EOM-CCSD). Orbital optimization within pCCD, however, breaks the asymptotic degeneracy of the $2^{1}\Sigma$, $1^{1}\Delta$, and $2^{1}\Pi$ states and results in similar errors as conventional EOM-CCSD. Nonetheless, the proposed EOM-pCCD-LCCSD formalism significantly improves both excitation energies and PESs for excited states that feature a biexcited character compared to the conventional EOM-CCSD model.

Our second test system covers all-trans polyenes comprising two to seven double bonds. All-trans polyenes are considerably challenging to model theoretically, primarily because doubly-excited configuration are required to accurately describe ground and excited states of longer polyenes 22. Furthermore, the proper description of the two lowest-lying excited states poses a problem to both experiment and quantum chemistry approaches 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36. From a theoretical point of view, the lowest-lying excited states of all-trans polyenes have to be studied using highly-accurate electronic structure methods, like complete-active-space second-order perturbation theory (CASPT2) and the density-matrix renormalization-group algorithm (DMRG). Such calculations predict the $2^{1}A_{g}^{-}$ state to be the lowest excited state and suggest a significant double excitation character with a dominant HOMO${}^{2}\rightarrow$LUMO2 transition.

The excitation energies for the two lowest-lying vertical excited states of \ceC4H6 to \ceC14H16 obtained by various EOM-CC methods using the DFT-optimized structures of ref. 15 are summarized in Table 2. Conventional multi-reference methods like CASSCF and MRMP predict the correct order of states. As expected, EOM-CCSD overestimates the excitation energies of the $2^{1}A_{g}^{-}$ state and suggests that the dark state lies above the bright state by 0.9 eV. Similar results are obtained from its linearized version EOM-LCCSD and the simple pCCD-based excited state model EOM-pCCD+S. We should note that we encountered convergence difficulties in EOM-LCCSD for longer polyene chain lengths. Thus, the corresponding excitation energies are not shown in the Table. In contrast to the conventional EOM-CCSD method and its simplified variants, the EOM-pCCD-LCCSD approach predicts the first dark state to be the lowest-lying excited state. Although being similar in computational cost as conventional EOM-CCSD, EOM-pCCD-LCCSD yields the correct order of states. We should emphasize that all EOM-CC methods perform equally good for the first bright state $1^{1}B_{u}^{+}$, whose state is mostly characterized by the single HOMO${}^{1}\rightarrow$LUMO1 transition. Differences with respect to EOM-CCSD typically amount up to 0.3 eV, while the largest deviations are found for EOM-LCCSD (from 0.4 to 1.1 eV).

The influence of relaxation effects of the molecular structure on the excitation energies are summarized in Table 3. Note that the vertical and adiabatic excitation energies are calculated for the DMRG-optimized molecular geometries and hence do not equal the corresponding adiabatic excitation energies of each EOM-CC method. We should note that two different active spaces have been used in DMRG calculations. Since the EOM-CC excitation energies are similar for both DMRG structures, only one set of excitation energies is presented in the Table (see also Table S3 of the Supporting Information). As observed for the DFT-optimized structures, all studied EOM-CC methods yield similar (vertical and adiabatic) excitation energies for the first bright state. In general, EOM-LCCSD and DMRG yield excitation energies that differ most from experimental reference data (1.5 to 2.0 eV), followed by EOM-pCCD-LCCSD (approximately 1.2 eV) and EOM-CCSD (approximately 1 eV).

In contrast to the first bright state, the accurate description of the first dark state is more challenging. While EOM-pCCD+S, EOM-LCCSD, and EOM-CCSD yield the wrong order of states, DMRG and EOM-pCCD-LCCSD correctly predict the dark state to be the lowest-lying excited state. Compared to experimental results, EOM-CCSD overestimates excitation energies by 2.8 eV (vertical) and 2.1 eV (adiabatic), respectively, while DMRG yields differences between 2-2.4 eV (vertical) and 0.3-1.0 eV (adiabatic), respectively. Most importantly, the EOM-pCCD-LCCSD excitation energies deviate less from experimental data (1.7-1.9 eV for vertical, 0.06-0.5 eV for adiabatic excitations). The large differences for \ceC14H16 might originate from the fact that the corresponding ground and excited state molecular structures have not been optimized for the pCCD-LCCSD wave function or that pCCD-LCCSD is insufficient to describe the electron correlation effects in ground and excited states.

Furthermore, all ground and excited state calculations have been performed using the cc-pVDZ basis set to allow for a direct comparison to DMRG reference calculations. For a better comparison to experiment, larger basis sets (of triple-zeta quality and augmented functions) should be applied. Nonetheless, despite of the small size of the basis set and the errors in excitation energies, EOM-pCCD-LCCSD significantly outperforms conventional EOM-CC methods restricted to at most double excitations and predicts the right order of states in all-trans polyenes.

The good performance of EOM-pCCD-LCCSD compared to the remaining EOM-CC models can be explained by analyzing the contributions of doubly- and singly-excited configurations, where the overall weight of all doubly-excited configurations may increase to approximately 80 % for longer polyene chain lengths. 35 Fig. 2(a) shows the weight of singly-excited and doubly-excited configurations of the $2^{1}A_{g}^{-}$ state with respect to the number of double bonds (DFT-optimized structures) obtained by EOM-CCSD, its linearized version EOM-LCCSD, and the pair-amplitude tailored EOM-pCCD-LCCSD variant. Both EOM-CCSD and EOM-LCCSD predict the first dark state to feature 70 to 80 % singly-excited configurations. For shorter polyene chain lengths, EOM-LCCSD behaves similar to EOM-CCSD. If the number of double bonds exceeds 4, the contribution of singly-excited configurations rapidly decreases, while the weight of doubly-excited configurations significantly increases. For \ceC12H14, both types of configurations have similar weights of approximately 50 %. By contrast, EOM-pCCD-LCCSD correctly predicts the individual weights of singly- and doubly-excited configurations. For an increasing number of $\pi$ bonds, the contribution of doubly-excited configurations gradually accumulates approaching 80 %, while the weight of singly-excited configurations decreases accordingly. To conclude, the character of the first dark state is appropriately described by the proposed EOM-pCCD-LCCSD model and lies energetically below the first bright state. Similar trends can be observed for the DMRG-optimized structures (see Fig. 2(b) and (c)). Specifically, EOM-LCCSD predicts larger weights of singly-excited configurations in the adiabatic $2^{1}A_{g}^{-}$ state, which may explain the poor performance of EOM-LCCSD for adiabatic excitation energies, while the contribution of singly-excited configurations in EOM-pCCD-LCCSD decreases from 20 % to about 10 % if the excited state structure is relaxed, while the weight of the reference state increases to 10-20 %.

We should emphasize that EOM-LCCSD yields an adiabatically excited $2^{1}A_{g}^{-}$ state of \ceC10H12 that features similar excitation contributions as obtained within EOM-pCCD-LCCSD, while the corresponding excitation energy lies between EOM-pCCD-LCCSD and EOM-CCSD results. However, this rather good performance of EOM-LCCSD has been observed only for this particular molecular geometry and might be coincidental. For all investigated all-trans polyenes with an increasing number of double bonds, LCCSD and EOM-LCCSD featured slow convergence in contrast to EOM-pCCD-LCCSD and EOM-CCSD calculations. Furthermore, the excitation contributions of the EOM-LCCSD excited states are more sensitive to the molecular geometry than those obtained in EOM-pCCD-LCCSD and EOM-CCSD calculations.

4 Conclusions

We have presented a new excited state model based on the pCCD wave function using the equation-of-motion formalism. Our method represents a simplified version of EOM-CCSD where the electron-pair amplitudes are tailored by the pCCD cluster amplitudes and non-pair amplitudes enter only linearly into the CC and EOM equations. The proposed EOM-pCCD-LCCSD model has been benchmarked against excitation energies that feature a doubly-excited character. In general, EOM-pCCD-LCCSD allows for an improved description of such excited states including at most double excitations in the CC ansatz. Specifically, errors in excitation energies are reduced to 0.06–0.6 eV with respect to FCI or experimental reference data.

We should emphasize that pCCD calculations are commonly combined with an orbital optimization protocol (for the ground state) to recover size-consistency, while post-pCCD calculations are performed within the pCCD-optimized orbital basis. The corresponding natural orbitals are typically symmetry-broken (localized), which prevents us from identifying excited states according to their spatial symmetry, and biased toward the ground state. Symmetry-breaking can be prevented by imposing (point-group) symmetry in the orbital-optimization protocol.

To further improve both excitation energies and PESs and to eliminate the orbital bias toward the ground state, the EOM-pCCD-LCCSD model can be combined with an orbital optimization protocol for excited states. This is currently being investigated in our laboratory.

Computational Methodology

The molecular structures of all investigate all-trans polyenes containing 2 to 7 $\pi$-bonds were taken from ref. 15 (DFT-optimized structures, vertical excitations) and ref. 39 (DMRG-optimized structures, vertical and adiabatic excitations). For direct comparison with DMRG reference data, we used the cc-pVDZ basis set of Dunning 41 in all calculations. For the \ceCH^+ ion, the basis set was taken from ref. 42. All pCCD, EOM-pCCD-LCCSD, and EOM-LCCSD calculations were performed with our in-house quantum-chemistry code employing restricted canonical Hartree–Fock orbitals (\ceCH^+ and polyenes) and pCCD natural orbitals optimized for the pCCD ground state reference function (\ceCH^+, imposing $C_{2h}$ point group symmetry) as molecular orbital basis. The EOM-CCSD calculations were performed with the Molpro program suite 43, 44.

Acknowledgement

K.B. acknowledges financial support from a Marie-Skłodowska-Curie Individual Fellowship project no. 702635–PCCDX and a scholarship for outstanding young scientists from the Ministry of Science and Higher Education. Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), grant no. 412.

Supporting Information

Diagrammatic representation of EOM-pCCD-LCCSD, EOM-pCCD-LCCSD total electronic energies and excitation energies for the lowest-lying excited states of \ceCH+, and vertical and adiabatic excitation energies for \ceC10H12 to \ceC14H16 using all DMRG-optimized structures.