Physics Medicine
Computer Science-2021.gz:143377
Computer Science
cdeb8a654c8542ece1b1041ad4ea8dacbf271e81
Multireference Algebraic Diagrammatic Construction Theory for Excited States: Extended Second-Order Implementation and Benchmark

We present an implementation and benchmark of new approximations in multireference algebraic diagrammatic construction theory for simulations of neutral electronic excitations and UV/Vis spectra of strongly correlated molecular systems (MR-ADC). Following our work on the first-order MR-ADC approximation [J. Chem. Phys. 2018, 149, 204113], we report the strict and extended second-order MR-ADC methods (MR-ADC(2) and MR-ADC(2)-X) that combine the description of static and dynamic electron correlation in the ground and excited electronic states without relying on state-averaged reference wavefunctions. We present an extensive benchmark of the new MR-ADC methods for excited states in several small molecules, including the carbon dimer, ethylene, and butadiene. Our results demonstrate that for weakly-correlated electronic states the MR-ADC(2) and MR-ADC(2)-X methods outperform the third-order single-reference ADC approximation and are competitive with the results from equation-of-motion coupled cluster theory. For states with multireference character, the performance of the MR-ADC methods is similar to that of an N-electron valence perturbation theory. In contrast to conventional multireference perturbation theories, the MR-ADC methods have a number of attractive features, such as a straightforward and efficient calculation of excited-state properties and a direct access to excitations outside of the frontier (active) orbitals.


Introduction
Ab initio electronic structure methods 1-4 are often used in combination with experiment to gain a deeper understanding of the excited-state properties of molecules and to study the behavior of these molecules as they undergo various photochemical transformations. Among many quantum chemical methods, algebraic diagrammatic construction theory [5][6][7][8][9] (ADC) has become one of the widely used approaches for studying electronically excited states, [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] due to a combination of its low computational cost and systematically improvable accuracy. In its original formulation proposed almost 40 years ago, 5 ADC is based on the conventional (single-reference) perturbation theory, which significantly simplifies its working equations and allows for efficient calculation of excitation energies and excited-state properties. However, due to its single-reference nature, conventional ADC approximations become unreliable in systems with strong electron correlation in the ground or excited electronic states.
In this work, we present the implementation and benchmark of the strict and extended secondorder MR-ADC methods (MR-ADC(2) and MR-ADC(2)-X) for neutral electronic excitations that provide a more balanced description of static and dynamic electron correlation effects. We present an overview of the MR-ADC theory for neutral excitations in Section 2, before outlining the principle components of our implementation in Section 3. We lay out the details of our computations in Section 4 and present our benchmark results for excited states of small molecules (HF, CO, N 2 , F 2 , H 2 O), the carbon dimer (C 2 ), as well as the ethylene (C 2 H 4 ) and butdiene (C 4 H 6 ) molecules in Section 5. Our results demonstrate that the secondorder MR-ADC approximations provide accurate predictions of the excited-state energies for weaklyand strongly-correlated electronic states at equilibrium geometries and near dissociation. We outline our conclusions in Section 6.

MR-ADC for the Polarization Propagator
We begin with a brief overview of the principal components of the MR-ADC theory in the context of neutral (particle-number-conserving) electronic excitations. A more detailed discussion of MR-ADC using the formalism of the effective Liouvillean theory 64 can be found in Ref. 26. Electronic excitations arise from the perturbation of a molecule by an external electric field. The response of the molecule is described by the spectral function 5 where D is a dipole moment operator matrix with elements D pq expressed in a finite single-particle basis set and Π(ω) is a polarization propagator, which is the central object of interest in this work. The polarization propagator Π(ω) is a particular type of the retarded frequency-dependent propaga-tor 3,4 G(ω), which can be generally written as: where G + µν (ω) and G − µν (ω) are the forward and backward components of the propagator, H is the electronic Hamiltonian, |Ψ and E are the eigenfunction and eigenenergy of H, respectively. The polarization propagator Π(ω) is characterized by the choice of q † ν = a † p a q − Ψ| a † p a q |Ψ , expressed in terms of the fermionic annihilation and creation operators, and the negative sign for the second term in Eq. (2). In the case of the polarization propagator, the forward and backward components are related as Π † + (−ω) = −Π − (ω), so in computing one component, the other is trivially obtained. For this reason, in this work we define Π(ω) ≡ Π + (ω) and consider only the forward component of the propagator henceforth.
When written in the equivalent Lehmann 65 (spectral) representation, in which a resolution of identity is introduced over all excited eigenstates of the system (|Ψ k ), the poles of the polarization propagator Π pqrs (ω) correspond to excitation energies (ω k = E k − E 0 ) while the expectation values in the numerator describe the probabilities of the corresponding transitions. The frequency ω ≡ ω + iη is defined in terms of real (ω ) and infinitesimal imaginary (iη) components. Eq. (3) can be written compactly in a matrix form whereΩ is a diagonal matrix of excitation energies (ω k ) andX is a matrix of spectroscopic amplitudes X pqk = Ψ|a † q a p |Ψ k . The expressions for the polarization propagator in Eqs. (3) and (4) are exact if written in the basis of the exact eigenstates, however in practice these states must be approximated.
MR-ADC provides a hierarchy of efficient approximations to the polarization propagator (4) by expressing it in a time-independent multireference perturbative expansion and truncating this series at a low order in perturbation theory. Working equations of the MR-ADC approximations are de- In Eqs. (5) and (6), T is an excitation operator that generates all of the internally-contracted excitations between the core, active, and external orbitals ( Figure 1). Next, the full electronic Hamiltonian H is partitioned into the zeroth-order part H (0) and a perturbation V = H −H (0) , where H (0) is the Dyall Hamiltonian 33-35,72 with an eigenfunction |Ψ 0 . This partitioning of H leads to the perturbative expansions for the wavefunction |Ψ and the polarization propagator Π(ω): Truncating these expansions at the n th order defines the polarization propagator of the MR-ADC(n) approximation. The MR-ADC(n) working equations are derived by writing Π(ω) in the form where M is a non-diagonal effective Hamiltonian matrix that can be used to compute the excitation energies and T is an effective transition mo-ments matrix that contains information about the probabilities of electronic transitions. The matrices M and T are evaluated in an approximate basis of many-electron (internally-contracted) states, which are in general non-orthogonal with an overlap matrix S. In Eq. (9) each matrix is evaluated up to the same (n th ) order in perturbation theory. The MR-ADC(n) excitation energies Ω and complementary eigenvectors Y are obtained by solving the generalized eigenvalue problem Combining the eigenvectors Y with the effective transition moments matrix T allows to compute the spectroscopic amplitudes and the MR-ADC(n) polarization propagator along with the spectral function T (ω) in Eq. (1). The spectroscopic amplitudes X can be also used to compute the oscillator strength for each electronic transition k where Ψ|D|Ψ k is the transition dipole moment matrix element.

Perturbative Structure of the MR-ADC Matrices
In this work, we consider the strict and extended second-order MR-ADC approximations for electronic excitations, denoted as MR-ADC(2) and MR-ADC(2)-X. In the MR-ADC(2) approximation, contributions to Π(ω) are included strictly up to the second order in the MR-ADC perturbation expansion, while the MR-ADC(2)-X method partially incorporates selected third-order terms as discussed below. At each perturbation order n, contributions to the matrix elements of M, T, and S in Eq. (9) are (k = 0, 1) operators (Eqs. (17) and (18)) on the reference state |Ψ 0 . A circle connected with an arrow denotes a single excitation. Black, green, and red lines represent core, active, and external orbitals, respectively. expressed as whereH (l) andq (l) µ are the l th -order operators obtained from the perturbative analysis of the Baker-Campbell-Hausdorff expansions for the effective HamiltonianH = e −A He A and observableq µ = e −A q µ e A operators, respectively. The excitation operators h (m) † ν are used to construct the manyelectron basis states that are necessary for representing the wavefunctions of the excited electronic states. Defining a p q ≡ a † p a q and a pq rs ≡ a † p a † q a s a r , the low-order h (m) † ν used in MR-ADC(2) and MR-ADC(2)-X have the form: h (1) † = {a ab ij , a ax ij , a ab ix , a xy ij , a ab xy , a yz ix , a az xy , a ay ix } Eqs. (14) and (15) define the perturbative structure of the M and T matrices shown in Figure 3, with the effective operatorsH (k) andq (l) expanded to different orders in different sectors defined by the h

Amplitudes of the Effective Hamiltonian
Computation of the M and T matrix elements in Figure 3 requires solving for the amplitudes of the excitation operator T (k) (Eq. (6)). We express T (k) in a general form as where t (k) µ are the kth-order amplitude coefficients and τ † µ denotes the corresponding string of creation and annihilation operators. In MR-ADC(2) and MR-ADC(2)-X only the low-order amplitudes t (1) and t (2) need to be evaluated. These amplitudes are determined by solving a system of projected linear equations For the details of solving these amplitude equations we refer the readers to our earlier publications. 62,63 Both MR-ADC(2) and MR-ADC(2)-X require calculating all first-order amplitudes, specifically: The first three sets of amplitudes correspond to single excitations (t x ) and the next five sets represent external double excitations. The last three sets of t (1) in Eq. (21) describe the semiinternal double excitations (t The MR-ADC(2) and MR-ADC(2)-X equations for the effective Hamiltonian matrix M also depend on the second-order singles and semi-internal doubles amplitudes of the form: We note that in addition to the external singles am- x ), Eq. (22) also includes the internal singles amplitudes t y (2) x (x > y) that we did not encounter in our previous work. These amplitudes are necessary to ensure the Hermiticity of the effective Hamiltonian matrix M as discussed below. Additionally, the MR-ADC(2)-X equations for the T matrix elements (Figure 3b) depend on the external double-excitation amplitudes (2) xy (23) In our previous work on EA/IP-MR-ADC(2) and EA/IP-MR-ADC(2)-X, 62,63 we demonstrated that neglecting all second-order amplitudes except t a(2) i has a negligible effect on the computed charged excitation energies and transition probabilities. However, for neutral electronic excitations, some of the second-order amplitudes need to be computed to prevent asymmetry of the M matrix. To illustrate this, we consider the second-order contributions to one of the off-diagonal blocks of M and its trans-pose that can be expressed and simplified as: The two blocks of the M matrix are equal to each other provided that the second term in Eq. (25) is zero, which is satisfied if the second-order effective HamiltonianH (2) is parametrized with the t b(2) x amplitudes that fulfill the projected amplitude (20)). As another example we consider one of the diagonal blocks of M: Here, the last two terms in Eqs. (26) and (27) are in general different, unless the effective Hamiltonian includes contributions from the internal second-order amplitudes t y(2) x (x > y) ensuring that ix,jy block. The solution of the amplitude equations for the t y (2) x (x > y) amplitudes is described in the Appendix.
Including all second-order singles amplitudes ensures the full Hermiticity of the M matrix in all sectors. On the other hand, neglecting the secondorder doubles amplitudes does not affect the symmetry of M and has a negligible effect on the computed excitation energies and oscillator strengths (see Table S1 of the Supporting Information for details). For this reason, in our implementation of MR-ADC for electronic excitations, we approximate:

Implementation
Our implementation of MR-ADC(2) and MR-ADC(2)-X for electronic excitations follows the steps outlined below: 1. Given an active space, compute the reference (usually, ground-state) CASSCF wavefunction |Ψ 0 and molecular orbitals.
2. Using the reference CASSCF molecular orbitals compute the CASCI wavefunctions |Ψ I (I > 0) for the user-defined number of excited states (N CI ).
3. Calculate the reduced density matrices  20)) to compute the first-and secondorder amplitudes t (1) and t (2) discussed in Section 2.2.2.
As discussed in our previous work, 62,63 the MR-ADC generalized eigenvalue problem is solved in the symmetrically-orthogonalized form (MỸ = Y Ω,M = S −1/2 M S −1/2 ,Ỹ = S 1/2 Y) using the multiroot Davidson algorithm 75,76 that iteratively optimizes the eigenvectorsỸ until convergence by forming the matrix-vector products σ =MỸ. We note that the S −1/2 matrix has a block-diagonal structure that allows to avoid storing the full matrix in memory. The small non-diagonal blocks of S −1/2 are computed very efficiently by diagonalizing small sectors of the overlap matrix S. 62 Equations and computer code for the calculation of the σ vectors were generated using a modified version of the SecondQuantizationAlgebra program (SQA) developed by Neuscamman and co-workers. 77 To optimize computational efficiency, our code generator implements the tensor contractions by automatically separating them into the efficient tensor intermediates, which are reused throughout the calculation.
For a fixed active space, the MR-ADC(2) and MR-ADC(2)-X implementations have the O(M 5 ) and O(M 6 ) computational scaling with the basis set size M , respectively, which is equivalent to the scaling of their single-reference counterparts. 9 Both methods have the formal O(N det N 8 act ) scaling with the number of active orbitals (N act ) and Slater determinants (N det ), similar to that of the conventional second-order multireference perturbation theories. The active-space scaling can be further lowered to O(N det N 6 act ) using a technique described in our previous work, 62,63 which forms efficient intermediates to avoid computation of the four-particle RDMs, although we did not employ this approach in our current implementation.

Computational Details
The MR-ADC(2) and MR-ADC(2)-X methods were implemented in Prism, a standalone program that is being developed in our group. The Prism program interfaces with PySCF 78 to obtain integrals and CASSCF/CASCI reference wavefunctions. The MR-ADC results were benchmarked against excitation energies obtained from the semi-stochastic heat-bath configuration interaction (SHCI) 79-81 calculations extrapolated to the full configuration interaction limit. The SHCI method was implemented in the Dice program. [79][80][81] We also compare the MR-ADC results to those computed using strongly-contracted N -electron valence second-order perturbation theory (sc-NEVPT2), 33,34 single-reference ADC (ADC(n), n = 2, 3), 5,9 and equation-of-motion coupled cluster theory with single and double excitations (EOM-CCSD). [82][83][84][85] We used PySCF to obtain the sc-NEVPT2 energies. The ADC(2), ADC(3), and EOM-CCSD results were computed using Q-Chem. 86 Performance of the MR-ADC(2) and MR-ADC(2)-X methods was tested for a benchmark set of five small molecules (HF, F 2 , CO, N 2 , and H 2 O), a carbon dimer (C 2 ), and two alkenes: ethylene and butadiene. As in our earlier work, 62,63 computations of the five small molecules were performed using two geometries: equilibrium and stretched. The equilibrium geometries were taken from Ref. 13. The stretched geometries of all molecules were obtained by increasing their bond lengths by a factor of two. The C−C bond distance in the carbon dimer was set to 2.4 a 0 , which is close to its equilibrium bond length. 87 The geometries of the alkenes were obtained as described in Ref. 88 and are reported in the Supporting Information. Calculations of HF, F 2 , CO, N 2 , and H 2 O employed the aug-cc-pVDZ basis set. 89 Carbon dimer calculations were performed using cc-pVDZ. 90 For the ethylene and butadiene molecules we used the ANO-L-PVTZ and ANO-L-PVDZ basis sets, respectively. 91 The active spaces used in the MR-ADC and sc-NEVPT2 computations (Section 5.1) were chosen as follows: where the molecular orbital occupations shown correspond to the those in the Hartree-Fock Slater determinant. The active spaces used in the calculations of the alkene molecules ranged from the full valence π-orbital space up to the triple-π space for ethylene and up to double-π for butadiene, as described in the Supporting Information. The sc-NEVPT2 computations were performed using the state-averaged CASSCF reference wavefunctions, where the ground and all excited states of interest were averaged with equal weights. All MR-ADC calculations were performed using the state-specific (ground-state) CASSCF reference wavefunctions. The MR-ADC calculations were converged with the number of CASCI states by including at least 15 excited singlet CASCI states and no less than 10 excited triplet CASCI states. The η s = 10 −5 and η d = 10 −10 truncation parameters were used to eliminate redundant excitations in the solution of the MR-ADC equations. 62,63 As in our previous work, 62,63 we have confirmed that the com- Table 1: Vertical excitation energies (Ω, eV) and oscillator strengths (f ) of molecules with equilibrium geometries. See Section 4 for details of the calculations. Also shown are mean absolute errors (∆ MAE ), standard deviations (∆ STD ), and maximum absolute errors (∆ MAX ) of the energies, relative to SHCI. Only orbitals participating in excitations are listed in the electron configuration of each excited state. puted MR-ADC excitation energies and oscillator strengths are size-consistent (see the Supporting Information for details).

Excited States in HF, CO, N 2 , F 2 , and H 2 O
We begin our discussion of results by benchmarking the MR-ADC(2) and MR-ADC(2)-X excitation energies against the accurate excited-state energies computed using the semi-stochastic heat-bath CI algorithm (SHCI) 79-81 for five small molecules (HF, CO, N 2 , H 2 O, and F 2 ) at equilibrium and stretched geometries. In addition to SHCI, we compare the performance of our MR-ADC methods with the conventional (single-reference) ADC approximations (ADC(n), n = 2, 3), equation-ofmotion coupled cluster theory with single and double excitations (EOM-CCSD), as well as stronglycontracted second-order N-electron valence perturbation theory (sc-NEVPT2). Table 1 compares vertical excitation energies of molecules at their near-equilibrium geometries computed using six approximate methods with the reference results from SHCI. EOM-CCSD shows the best performance for this benchmark set, with the mean absolute error ∆ MAE = 0.12 eV and a Table 2: Vertical excitation energies (Ω, eV) and oscillator strengths (f ) of molecules with stretched geometries. See Section 4 for details of the calculations. Also shown are mean absolute errors (∆ MAE ), standard deviations (∆ STD ), and maximum absolute errors (∆ MAX ) of the energies, relative to SHCI. Only orbitals participating in excitations are listed in the electron configuration of each excited state.     Table 2 shows results for electronic excitations of molecules with stretched geometries where the description of multireference effects is important. Unsurprisingly, the accuracy of results produced by the single-reference methods decreases dramatically. In particular, the errors of the singlereference ADC(n) methods become larger with increasing perturbation order (n) as ∆ MAE roughly doubles from ADC (2)  The performance of all methods for this benchmark set is summarized in Figure 4. Overall, our results demonstrate that MR-ADC(2) and MR-ADC(2)-X produce accurate results for many electronic transitions outperforming the singlereference ADC(2) and ADC(3) methods at nearequilibrium geometries and showing accuracy similar to sc-NEVPT2 near dissociation. In contrast to sc-NEVPT2, the MR-ADC methods allow for a straightforward calculation of transition properties (such as oscillator strengths) and do not require using state-averaged CASSCF reference wavefunctions.

Carbon Dimer
Next, we consider the C 2 molecule, known for the multireference nature of its ground and excited states, which also require accurate description of dynamic correlation. 81,87,[92][93][94][95] In our earlier work, 26 we demonstrated that the first-order MR-ADC approximation (MR-ADC(1)) captures the strongly correlated nature of the C 2 excited states, but produces large errors in vertical excitation energies due to the missing description of the two-electron dynamic correlation effects. Here, we reinvestigate C 2 using MR-ADC(2) and MR-ADC(2)-X that provide a more balanced description of static and dynamic correlation. Table 3 reports the MR-ADC(2), MR-ADC(2)-X, and sc-NEVPT2 vertical excitation energies for C 2 at a near equilibrium bond distance (r = 2.4 a 0 ) computed using the cc-pVDZ basis set and (8e, 8o) active space. The results of these multireference methods are compared to the accurate excitation energies from density matrix renormalization group (DMRG) computed by Wouters et al. 87 Figure 5 summarizes the performance of each method for predicting the vertical excitation energies of C 2 . The best agreement with DMRG is shown by MR-ADC(2)-X that reproduces the reference energies with ∆ MAE and ∆ STD of ∼ 0.05 eV. The sc-NEVPT2 method (∆ MAE and ∆ STD of 0.08 eV) exhibits significantly worse accuracy than MR-ADC(2)-X, but performs slightly better than MR-ADC(2) (∆ MAE = 0.11 eV, ∆ STD = 0.09). The MR-ADC(2) and MR-ADC(2)-X mean absolute errors are reduced tenfold compared to that of MR-ADC(1), 26 indicating that the secondorder contributions to the MR-ADC propagator are absolutely critical for quantitive predictions of the C 2 excited-state energies. As illustrated in Figure 5a, MR-ADC(2)-X has a fairly consistent Table 3: Carbon dimer vertical excitation energies (Ω, eV) and oscillator strengths (f ) computed using the cc-pVDZ basis set with r(C−C) = 2.4 a 0 . The MR-ADC and sc-NEVPT2 calculations were performed using the state-specific and state-averaged CASSCF reference wavefunctions, respectively, with the same (8e, 8o) active space. Also shown are mean absolute errors (∆ MAE ), standard deviations (∆ STD ), and maximum absolute errors (∆ MAX ) in the computed excitation energies, relative to the DMRG results from Ref. 87.
performance for all excited states resulting in a small ∆ STD , while MR-ADC(2) shows significantly larger errors (∼ 0.2 eV) for the 3 Σ + u , 1 ∆ g , and 1 Σ + g states. The errors of sc-NEVPT2 gradually increase reaching ∼ 0.3 eV for the highest-energy 3 ∆ u state. The MR-ADC methods predict only one dipole-allowed transition ( 1 Π u ) with an oscillator strength of 0.003, in agreement with the selection rules of optical spectroscopy.

Alkenes: Ethylene and Butadiene
Finally, we benchmark MR-ADC(2) and MR-ADC(2)-X for the low-lying excited states of ethylene (C 2 H 4 ) and butadiene (C 4 H 6 ). Both molecules have a low-lying triplet state ( 3 B 1u or 3 B u ), as well as a dipole-allowed π − π * transition to the singlet state ( 1 B 1u or 1 B u ) that requires a very accurate description of dynamic correlation between the σ and π electrons. [96][97][98][99][100] In addition, the butadiene molecule features a dipole-forbidden tran- sition to the 1 A g state with a substantial doubleexcitation character, requiring an accurate description of static correlation in the π and π * orbitals. For this reason, excitation energies of the butadiene 1 B u and 1 A g states are very sensitive to the level of electron correlation treatment and have been a subject of many theoretical studies. 42,88,[98][99][100][101][102][103][104][105][106][107][108][109][110][111] In this work, we compute the ethylene and butadiene excited states using MR-ADC(2) and MR-ADC(2)-X and compare their excitation energies to accurate reference results from SHCI obtained by Chien et al. 100 for the same molecular geometries and basis set (see Section 4 for details). Table 4 reports vertical excitation energies for the 3 B u and 1 B u excited states of ethylene computed using MR-ADC(2), MR-ADC(2)-X, sc-NEVPT2, and EOM-CCSD, along with the reference results from SHCI. For the multireference methods we employ three different active spaces ranging from the full-π (2e, 2o) to the triple-π (2e, 6o) active space (see the Supporting Information for details). The EOM-CCSD excitation energies are in a good agreement with SHCI with errors less than 0.15 eV, indicating that a highlevel description of dynamic correlation is sufficient for accurate description of the 3 B u and 1 B u excited states. Using the minimal (2e, 2o) active space, all three multireference methods overesti-  (2)). Increasing the active space lowers the computed 3 B u excitation energy reducing the errors of all three methods to less than 0.1 eV. In contrast, including more π-orbitals in the active space increases the computed energies of the 1 B u state. Out of all multireference methods, sc-NEVPT2 shows the strongest active-space dependence of the 1 B u excitation energy that changes from 7.68 (2e, 2o) to 8.34 eV (2e, 6o) and has the worst agreement with SHCI for the largest active space with 0.29 eV error. MR-ADC(2)-X combined with the (2e, 6o) active space shows the best agreement with SHCI for the 1 B u state with a small 0.02 eV error. The performance of MR-ADC (2) is similar to that of sc-NEVPT2, although the former method shows a weaker active-space dependence. The 1 B u oscillator strengths computed using the two MR-ADC methods do not change significantly with increasing active space and agree well with each other. Vertical excitation energies of butadiene are shown in Table 5. As in our study of ethylene, we employ three active spaces for the multireference methods ranging from the full-π (4e, 4o) to the double-π (4e, 8o) space. Notably, for this system the EOM-CCSD method shows a large 0.7 eV  Figure 6: Vertical excitation energies for the butadiene dipole-allowed 1 B u (plot (a)) and dipole-forbidden 1 A g (plot (b)) electronic states computed using MR-ADC(2), MR-ADC(2)-X, and sc-NEVPT2 as a function of the number of active π-orbitals (n) with 4 active electrons (4e, no). The reference excitation energies computed using SHCI 100 are shown as dashed horizontal lines.
error for the dipole-forbidden 1 A g state, which is known to have a double-excitation character and requires a multireference treatment. The activespace dependence of the multireference methods for the 1 A g and 1 B u excited states is illustrated in Figure 6. When combined with the smallest (4e, 4o) active space, all three multireference methods outperform EOM-CCSD for the 1 A g excitation energy with errors less than 0.4 eV. Using the largest (4e, 8o) active space, the best results are shown by MR-ADC(2)-X and sc-NEVPT2 with errors of 0.12 eV and 0.19 eV, respectively. A stronger activespace dependence is observed for the dipole-allowed 1 B u state with the sc-NEVPT2 and MR-ADC(2)-X excitation energies varying by ∼ 0.6 eV. Employing the (4e, 8o) active space, the best results are shown by MR-ADC(2)-X that outperforms sc-NEVPT2 and EOM-CCSD with an error of 0.06 eV in the 1 B u relative energy. sc-NEVPT2 and MR-ADC(2) overestimate the 1 B u energy by 0.09 and 0.07 eV, respectively. The MR-ADC methods also show a significant active-space dependence in the computed 1 B u oscillator strength that changes by ∼ 15 % from (4e, 4o) to (4e, 8o). A weaker dependence on the active space is observed for the 3 B u state, where the best agreement with SHCI is shown by sc-NEVPT2. Overall, the presented results demonstrate that the MR-ADC methods are competitive in accuracy with sc-NEVPT2 for all three electronic states of butadiene while capable of describing the multireference nature of the 1 A g state that is not captured by EOM-CCSD.

Conclusions
In this work, we presented an implementation and benchmark of the strict and extended second-order multireference algebraic diagrammatic construction theory for neutral electronic excitations (MR-ADC(2) and MR-ADC(2)-X). The MR-ADC(2) approximation incorporates all contributions to the polarization propagator up to the second order in its time-independent multireference perturbative expansion, while the MR-ADC(2)-X method includes additional third-order terms in the description of double excitations outside of the active space. We benchmarked the performance of MR-ADC(2) and MR-ADC(2)-X for the excited states in five small molecules (HF, CO, N 2 , F 2 , and H 2 O) at equilibrium and stretched geometries, a carbon dimer (C 2 ), as well as ethylene (C 2 H 4 ) and butadiene (C 4 H 6 ).
Our results demonstrate that the MR-ADC methods provide accurate predictions of excitedstate energies for the weakly-and stronglycorrelated electronic states, at geometric equilibrium or near dissociation. For the weaklycorrelated systems, we find that the MR-ADC(2) and MR-ADC(2)-X methods outperform the third-order single-reference ADC approximation (ADC(3)) and are often competitive with the results from equation-of-motion coupled cluster theory (EOM-CCSD). For the multireference problems, the performance of MR-ADC(2) and MR-ADC(2)-X is similar to that of strongly-contracted N-electron valence perturbation theory (sc-NEVPT2). The MR-ADC methods have a number of added benefits, such as a straightforward and efficient calculation of excited-state properties, a direct access to excitations outside of the active space, and an ability to calculate multireference excited states without state-averaged reference wavefunctions.
Our current work can be extended in many possible directions. One important direction is an efficient implementation of the MR-ADC(2) and MR-ADC(2)-X methods that will enable calculations of larger molecular systems with more orbitals and electrons in the active space. Another extension is to take advantage of the MR-ADC ability to calculate excitations outside of the active space, which can be used for efficient calculations of the X-ray absorption spectra of multireference systems. Finally, the methods presented in the current work can be extended to account for relativistic effects, such as spin-orbit and spin-spin coupling, which will enable calculations of magnetic systems with heavy elements. Work along these directions is ongoing in our group.
7 Appendix: Equations for the t y(2) x

Amplitudes
As discussed in Section 2.2.2, the MR-ADC(2) and MR-ADC(2)-X second-order effective Hamiltoniañ H (2) incorporates terms that depend on the internal second-order amplitudes t y(2) x (x > y) that ensure Hermiticity of the effective Hamiltonian matrix M. The t y (2) x (x > y) amplitudes parametrize a contribution to the second-order excitation operator T (2) (Eq. (6)) that has the form Following our previous work, 26,62,63 the amplitudes t y(2) x (x > y) are obtained by solving the projected linear equations (20) that can be written as: Eq. (30) can be simplified and converted to a tensor form: where T (2) contains t y (2) x (x > y) and the elements of K and V (2) are defined as follows: K xy,wz = Ψ 0 |(a † x a y − a † y a x )[H (0) , a † z a w − a † w a z ]|Ψ 0 (32) Expression for V (2) can be found in Eq. (38) of Ref. 62. Eq. (31) has the form of the Newton-Raphson equation for the second-order parameters t y (2) x (x > y) that describe relaxation of the active-space orbitals in response to dynamical correlation from outside of the active space. It can be solved by diagonalizing the K matrix according to the generalized eigenvalue problem and inverting the amplitude equation (31) to obtain: In Eqs. (34) and (35) the elements of S are defined as S xy,wz = Ψ 0 |(a † x a y − a † y a x )(a † z a w − a † w a z )|Ψ 0 (36) andZ = S 1/2 Z.

Supporting Information Available
Benchmark of the approximation for the secondorder amplitudes in Eq. (28). Size-intensivity tests for MR-ADC(2) and MR-ADC(2)-X. Cartesian geometries of the ethylene and butadiene molecules. Active space orbitals in the ethylene and butadiene calculations.