Time-dependent N-electron valence perturbation theory with matrix product state reference wavefunctions for large active spaces and basis sets: Applications to the chromium dimer and all-trans polyenes

Sokolov, Alexander Yu.; Guo, Sheng; Ronca, Enrico; Chan, Garnet Kin-Lic

doi:10.1063/1.4986975

In earlier work [A. Y. Sokolov and G. K.-L. Chan, J. Chem. Phys. 144, 064102 (2016)], we introduced a time-dependent formulation of the second-order N-electron valence perturbation theory (t-NEVPT2) which (i) had a lower computational scaling than the usual internally contracted perturbation formulation and (ii) yielded the fully uncontracted NEVPT2 energy. Here, we present a combination of t-NEVPT2 with a matrix product state (MPS) reference wavefunction (t-MPS-NEVPT2) that allows us to compute uncontracted dynamic correlation energies for large active spaces and basis sets, using the time-dependent density matrix renormalization group algorithm. In addition, we report a low-scaling MPS-based implementation of strongly contracted NEVPT2 (sc-MPS-NEVPT2) that avoids computation of the four-particle reduced density matrix. We use these new methods to compute the dissociation energy of the chromium dimer and to study the low-lying excited states in all-trans polyenes (⁠ $C_{4} H_{6}$ to $C_{24} H_{26}$ ⁠), incorporating dynamic correlation for reference wavefunctions with up to 24 active electrons and orbitals.

I. INTRODUCTION

Recent advances in molecular electron correlation methods have made it possible to describe dynamic correlation in weakly correlated systems with several hundreds of atoms by taking advantage of accurate local correlation approximations.^1–4 Similarly, efficient methods for static correlation^5–11 now exist that allow us to simulate systems with a large number of strongly correlated open-shell orbitals.^12–15 Yet, despite these efforts, the challenge still remains to efficiently describe dynamic correlation in the presence of significant static correlation. This is a scenario one faces when treating some complicated chemical systems, such as transition metal compounds with multiple metals.^12,15 Technically, the challenges are the high computational cost of the existing canonical algorithms, algebraic complexity of the underlying equations,¹⁶ and numerical instabilities in the simulations.^17–20

The standard approach to electron correlation in multi-reference (strongly correlated) systems is to first carry out a high-level description of static correlation in a small subset of near-degenerate (active) orbitals, followed by a lower-level description of dynamic correlation in the remaining large set of core and external orbitals. Here, the main challenge is to combine the high-level and low-level descriptions without sacrificing their accuracy, while retaining a low computational cost. For this purpose, most multi-reference dynamic correlation theories use the internal contraction approximation,^21–23 which represents the complicated active-space wavefunction in terms of simpler quantities, such as the reduced density matrices (RDMs). While this approximation has been enormously useful in quantum chemistry,^21–36 internally contracted methods require the computation of high-order (three- and four-particle) RDMs. These become prohibitively expensive for larger active spaces.

We have recently proposed a time-dependent formulation of multi-reference perturbation theory³⁷ that efficiently describes dynamic correlation by representing the active-space wavefunction in terms of compact time-dependent quantities (active-space imaginary-time Green’s functions (GFs)). This does not introduce any additional approximations nor do high-order RDMs appear in the equations. The resulting time-dependent theory is equivalent to the fully uncontracted perturbation theory but in fact displays a lower computational scaling than the internally contracted approximations, particularly with respect to the number of active orbitals.

Our previous work described the implementation of the time-dependent second-order N-electron valence perturbation theory (t-NEVPT2) for complete active-space self-consistent field (CASSCF) reference wavefunctions. Here, we present a new implementation combining t-NEVPT2 with matrix product state (MPS) reference wavefunctions, thus allowing us to describe dynamic correlation in multi-reference systems with much larger active spaces. As we will demonstrate, the resulting t-MPS-NEVPT2 approach requires computation of up to two-particle imaginary-time Green’s functions that can be evaluated using the time-dependent density matrix renormalization group (td-DMRG) algorithm^5,8–11 in a highly parallel fashion. In addition, for comparison, we present a low-scaling MPS-based implementation of the strongly contracted NEVPT2 (sc-MPS-NEVPT2),^22,23,35 which does not require computation of the four-particle density matrix. To demonstrate the capabilities of these new methods, we apply them to compute the dissociation energy of the chromium dimer and to study the low-lying excited states in the all-trans polyenes.

This paper is organized as follows. We begin with a brief overview of t-NEVPT2 (Sec. II A) and matrix product state wavefunctions (Sec. II B). We then describe the details of our t-MPS-NEVPT2 implementation (Sec. II C) and outline the computational details (Sec. III). In Sec. IV, we use t-MPS-NEVPT2 to compute energies along the $N_{2}$ dissociation curve (Sec. IV A), the dissociation energy of the chromium dimer (Sec. IV B), and the vertical excitation energies in all-trans polyenes (Sec. IV C). In Sec. V, we present the conclusions of our work. The details of the efficient sc-MPS-NEVPT2 implementation are described in Appendix B.

II. THEORY

A. Time-dependent formulation of N-electron valence perturbation theory (t-NEVPT2)

We begin with a brief overview of multi-reference perturbation theory in its time-dependent form. Our starting point is a zeroth-order electronic wavefunction $| Ψ_{0} (⟩$ computed in a complete active space (CAS) of molecular orbitals. We require $| Ψ_{0} (⟩$ to be an eigenfunction of a zeroth-order Hamiltonian ${\hat{H}}^{(0)}$ ⁠. Following our previous work,³⁷ we assume that ${\hat{H}}^{(0)}$ is the Dyall Hamiltonian,³⁸ ${\hat{H}}^{(0)} = {\hat{H}}_{D}$ ⁠, defined as

{\hat{H}}_{D} = \sum_{i} 𝜀_{i} a_{i}^{†} a_{i} + \sum_{a} 𝜀_{a} a_{a}^{†} a_{a} + {\hat{H}}_{a c t},

(1)

where we partition the spin-orbitals into three sets: (i) core (doubly occupied) with indices i, j, k, l; (ii) active with indices u, v, w, x, y, z; and (iii) external (unoccupied) with indices a, b, c, d. The orbital energies $𝜀_{i}$ and $𝜀_{a}$ are defined as the core and external eigenvalues of the generalized Fock operator,

f_{p}^{q} = h_{p}^{q} + \sum_{r s} v_{p r}^{q s} γ_{s}^{r},

(2)

where $h_{p}^{q}$ and $v_{p q}^{r s}$ are the usual one- and antisymmetrized two-electron integrals, and $γ_{q}^{p} = ⟨ Ψ_{0} | a_{p}^{†} a_{q} | Ψ_{0} ⟩$ is the one-particle density matrix of $| Ψ_{0} (⟩$ ⁠. The ${\hat{H}}_{a c t}$ operator describes the two-electron interaction in the active space

{\hat{H}}_{a c t} = \sum_{x y} (h_{x}^{y} + \sum_{i} v_{x i}^{y i}) a_{x}^{†} a_{y} + \frac{1}{4} \sum_{x y w z} v_{x y}^{z w} a_{x}^{†} a_{y}^{†} a_{w} a_{z} .

(3)

Having specified $| Ψ_{0} (⟩$ and ${\hat{H}}^{(0)}$ ⁠, we now consider an expansion of the energy with respect to the perturbation $λ \hat{V} = λ (Ĥ - {\hat{H}}_{D})$ ⁠. Rather than expressing the energy in a time-independent perturbation series, as is commonly done in the Rayleigh-Schrödinger perturbation theory, we consider a time-dependent expansion³⁷ with respect to the perturbation operator $\hat{V} (τ) = e^{({\hat{H}}_{D} - E_{D}) τ} \hat{V} e^{- ({\hat{H}}_{D} - E_{D}) τ}$ ⁠. The second-order energy contribution can be written as

\begin{array}{rcl} E^{(2)} & = & - \int_{0}^{\infty} d τ ⟨Ψ_{0} | {\hat{V}}^{' †} (τ) {\hat{V}}^{'} (0) | Ψ_{0}⟩ \\ = & - \int_{0}^{\infty} d τ ⟨ Ψ_{0} | {\hat{V}}^{' †} e^{- ({\hat{H}}_{D} - E_{D}) τ} {\hat{V}}^{'} | Ψ_{0} ⟩, \end{array}

(4)

where $τ$ is the imaginary time and ${\hat{V}}^{'}$ denotes the part of $\hat{V}$ that contributes to the first-order wavefunction $| Ψ^{(1)} (⟩$ (⁠ ${\hat{V}}^{'} = \hat{Q} \hat{V}$ ⁠, $\hat{Q} = 1 - | Ψ_{0} (⟩ ⟨) Ψ_{0} |$ ⁠). Equation (4) is the Laplace transform of the second-order energy expression from the Rayleigh-Schrödinger perturbation theory, which yields the exact (uncontracted) energy of the second-order N-electron valence perturbation theory (NEVPT2).^22,23 However, unlike the standard uncontracted theory, the time-dependent energy expression (4) does not require the costly representation of the first-order wavefunction in a large space of determinants and can be evaluated very efficiently. Expanding the perturbation operator $\hat{V}$ into classes of excitations, the second-order energy can be expressed as a sum of 8 terms related to the numbers of holes and particles created in the core and external spaces (see Refs. 22 and 39 for a complete definition)

\begin{array}{rcl} E^{(2)} & = & E^{[0]} + E^{[+ 1]} + E^{[- 1]} + E^{[+ 2]} + E^{[- 2]} \\ + E^{[+ 1^{'}]} + E^{[- 1^{'}]} + E^{[0^{'}]} . \end{array}

(5)

The central quantities to compute in Eq. (5) are the one- and two-particle imaginary-time Green’s functions (1- and 2-GF) in the active space. Specifically, the E^[−1] and E^[+1] terms require computation of 1-GFs (e.g., $G (τ) = ⟨) Ψ_{0} | a_{x}^{†} (τ) a_{y} | Ψ_{0} ⟩)$ ⁠), while 2-GFs are necessary to compute the E^[−2], E^[+2], and $E^{[0^{'}]}$ contributions (e.g., $G (τ) = ⟨) Ψ_{0} | a_{w}^{†} (τ) a_{x}^{†} (τ) a_{y} a_{z} | Ψ_{0} ⟩)$ ⁠). The $E^{[- 1^{'}]}$ and $E^{[+ 1^{'}]}$ components formally involve the three-particle active-space Green’s function but can be efficiently evaluated by forming intermediates. Defining

{\hat{h}}_{p}^{†} = \sum_{x} {\tilde{h}}_{x}^{p} a_{x}^{†} + \frac{1}{2} \sum_{x y z} v_{x y}^{p z} a_{x}^{†} a_{y}^{†} a_{z},

(6)

{\tilde{h}}_{p}^{q} = h_{p}^{q} + \sum_{i} v_{p i}^{q i},

(7)

the $E^{[- 1^{'}]}$ and $E^{[+ 1^{'}]}$ contributions can be expressed as expectation values of the single-index operators

E^{[+ 1^{'}]} = - \int_{0}^{\infty} \sum_{i} e^{𝜀_{i} τ} ⟨ Ψ_{0} | {\hat{h}}_{i} (τ) {\hat{h}}_{i}^{†} | Ψ_{0} ⟩ d τ,

(8)

E^{[- 1^{'}]} = - \int_{0}^{\infty} \sum_{a} e^{- 𝜀_{a} τ} ⟨ Ψ_{0} | {\hat{h}}_{a}^{†} (τ) {\hat{h}}_{a} | Ψ_{0} ⟩ d τ .

(9)

The bottleneck of the time-dependent NEVPT2 (t-NEVPT2) algorithm, when $| Ψ_{0} (⟩$ is expanded in determinants, is the evaluation of the 2-GF that has $O (N_{τ} \times N_{d e t} \times N_{a c t}^{6})$ computational cost, where N_act is the number of active orbitals, N_det is the dimension of the CAS Hilbert space, and $N_{τ}$ is the number of time steps (⁠ $N_{τ}$ $\sim$ 10-20). As a result, t-NEVPT2 has a significantly lower scaling with N_act than the internally contracted NEVPT2 approaches,^22,23 which require computation of up to the four-particle reduced density matrix (4-RDM) with $O (N_{d e t} \times N_{a c t}^{8})$ cost. However, the computational cost of these two types of NEVPT2 formulations still scales exponentially with N_act, since the number of determinants in the zeroth-order wavefunction $N_{d e t} \sim O (e^{N_{a c t}})$ ⁠, and solving for $| Ψ_{0} (⟩$ becomes a bottleneck for large active spaces. In Secs. II B and II C, we will describe how this limitation can be ameliorated by expressing $| Ψ_{0} (⟩$ in the matrix product state (MPS) representation.

B. Matrix product state (MPS) wavefunctions and the density matrix renormalization group (DMRG) algorithm

In this section, we very briefly introduce matrix product state (MPS) wavefunctions, which will be used extensively in this work. Further details about MPS algorithms can be found in references such as Refs. 40 and 41. The MPS is a nonlinear wavefunction composed of a product of tensors where each tensor corresponds to one or more orbitals in the basis. The most common form of the MPS wavefunction is the one-site MPS written in the following form:

| Ψ (⟩ = \sum_{n_{1} \dots n_{k}} A^{n_{1}} A^{n_{2}} \dots A^{n_{p}} \dots A^{n_{k}} | n_{1} n_{2} \dots n_{k} (⟩,

(10)

where $| n_{1} n_{2} \dots n_{k} (⟩$ is a Slater determinant, $A^{n_{p}}$ is a tensor that corresponds to only one orbital p with occupancy $n_{p} \in {| (⟩, | ↑ (⟩, | ↓ (⟩, | ↑ ↓ (⟩}$ (also referred as the site p), and the total number of tensors $A^{n_{p}}$ equals the number of orbitals k in the basis (for a CAS wavefunction, k = N_act). For each p in the range ${2, \dots, k - 1}$ ⁠, $A^{n_{p}}$ is a matrix with fixed dimensions $M \times M$ ⁠. The dimensions of $A^{n_{1}}$ and $A^{n_{k}}$ are $1 \times M$ and $M \times 1$ ⁠, respectively. As a result, for a specified set of occupation numbers $n_{1} n_{2} \dots n_{k}$ ⁠, the product of tensors in Eq. (10) yields a scalar corresponding to the coefficient of the determinant $| n_{1} n_{2} \dots n_{k} (⟩$ in the expansion of the wavefunction $| Ψ (⟩$ ⁠. The parameter M, referred to as the bond dimension, controls the flexibility of the MPS, which increases as M increases.

The most common algorithm when working with MPS is the sweep algorithm, where linear algebra operations are carried out on a single tensor at a time, sweeping successively through the orbitals $1 \dots k$ ⁠. The density matrix renormalization group (DMRG) is the prototypical version of such an algorithm and corresponds to a variational energy minimization using a sweep through the tensors. At a given site in the sweep algorithm, the tensor that is being operated on yields the linear expansion coefficients of the wavefunction in a many-body renormalized basis. For example, at site p, we define the left and right renormalized bases as

| l_{α_{p}} (⟩ = \sum_{n_{1} \dots n_{p - 1}} {[A^{n_{1}} A^{n_{2}} \dots A^{n_{p - 1}}]}_{α_{p}} | n_{1} n_{2} \dots n_{p - 1} (⟩,

(11)

| r_{α_{p + 1}} (⟩ = \sum_{n_{p + 1} \dots n_{k}} {[A^{n_{p + 1}} A^{n_{p + 2}} \dots A^{n_{k}}]}_{α_{p + 1}} | n_{p + 1} n_{p + 2} \dots n_{k} (⟩,

(12)

and the total wavefunction becomes

| Ψ (⟩ = \sum_{l_{α_{p}} n_{p} r_{α_{p + 1}}} {[A^{n_{p}}]}_{α_{p} α_{p + 1}} | l_{α_{p}} n_{p} r_{α_{p + 1}} (⟩ .

(13)

In this context, we can think of $A^{n_{p}}$ as representing a “site wavefunction,” and operators projected into the renormalized space ${| l_{α_{p}} n_{p} r_{α_{p + 1}} (⟩}$ act on this site wavefunction as “site operators.” General numerical algorithms for wavefunctions can be converted into sweep algorithms for MPS by identifying wavefunctions with site wavefunctions and operators with site operators. This idea is used in the time-dependent DMRG (td-DMRG) algorithm that we employ below.

C. t-NEVPT2 with MPS reference wavefunctions (t-MPS-NEVPT2)

1. Overview of the algorithm

We now describe the implementation of t-NEVPT2 for the MPS reference wavefunctions, which we denote as t-MPS-NEVPT2. To aid the discussion, we first rewrite the t-NEVPT2 energy in a compact form

E^{(2)} = E_{e x t}^{(2)} + \sum_{[i]} E_{a c t}^{[i]},

(14)

where we combine terms involving the core and external orbitals

E_{e x t}^{(2)} = \frac{1}{4} \sum_{i j a b} \frac{v_{i j}^{a b} v_{a b}^{i j}}{𝜀_{i} + 𝜀_{j} - 𝜀_{a} - 𝜀_{b}} + \sum_{i a} \frac{{\tilde{h}}_{i}^{a} {\tilde{h}}_{a}^{i}}{𝜀_{i} - 𝜀_{a}} + 2 \sum_{i x y a} \frac{{\tilde{h}}_{a}^{i} v_{i x}^{a y} γ_{y}^{x}}{𝜀_{i} - 𝜀_{a}}

(15)

and introduce a shorthand notation $E_{a c t}^{[i]}$ for the active-space components of the energy contributions E^[i] in Eq. (5) (⁠ $i \in \{+ 1, - 1, + 2, - 2, + 1^{'}, - 1^{'}, 0^{'}\}$ ⁠). Each component $E_{a c t}^{[i]}$ can be expressed in the general form

E_{a c t}^{[i]} = - \int_{0}^{\infty} 𝜖^{[i]} (τ) G^{[i]} (τ) d τ,

(16)

where $G^{[i]} (τ)$ is the Green’s function tensor and $𝜖^{[i]} (τ)$ is the integral prefactor tensor that can be computed using the one- and two-electron integrals and orbital energies (⁠ $𝜀_{i}$ and $𝜀_{a}$ ⁠). Explicit equations for the elements of tensors $𝜖^{[i]} (τ)$ and $G^{[i]} (τ)$ are given in Table I.

TABLE I.

Equations for the elements of tensors $𝜖^{[i]} (τ)$ and $G^{[i]} (τ)$ used in the evaluation of the t-NEVPT2 active-space energy contributions $E_{a c t}^{[i]}$ in Eq. (16).

[i]	$𝜖^{[i]} (τ)$	$G^{[i]} (τ)$
[+1]	$𝜖_{x}^{y} (τ) = \frac{1}{2} \sum_{i j a} v_{i j}^{a y} v_{a x}^{i j} e^{(𝜀_{i} + 𝜀_{j} - 𝜀_{a}) τ}$	$G_{x}^{y} (τ) = ⟨ Ψ_{0} \| a_{x} (τ) a_{y}^{†} \| Ψ_{0} ⟩$
[−1]	$𝜖_{y}^{x} (τ) = \frac{1}{2} \sum_{i a b} v_{i y}^{a b} v_{a b}^{i x} e^{(𝜀_{i} - 𝜀_{a} - 𝜀_{b}) τ}$	$G_{y}^{x} (τ) = ⟨ Ψ_{0} \| a_{x}^{†} (τ) a_{y} \| Ψ_{0} ⟩$
[+2]	$𝜖_{x y}^{z w} (τ) = \frac{1}{8} \sum_{i j} v_{i j}^{z w} v_{x y}^{i j} e^{(𝜀_{i} + 𝜀_{j}) τ}$	$G_{x y}^{z w} (τ) = ⟨ Ψ_{0} \| a_{x} (τ) a_{y} (τ) a_{w}^{†} a_{z}^{†} \| Ψ_{0} ⟩$
[−2]	$𝜖_{z w}^{x y} (τ) = \frac{1}{8} \sum_{a b} v_{z w}^{a b} v_{a b}^{x y} e^{- (𝜀_{a} + 𝜀_{b}) τ}$	$G_{z w}^{x y} (τ) = ⟨ Ψ_{0} \| a_{x}^{†} (τ) a_{y}^{†} (τ) a_{w} a_{z} \| Ψ_{0} ⟩$
$[0^{'}]$	$𝜖_{y z}^{x w} (τ) = \sum_{i a} v_{i z}^{a w} v_{a y}^{i x} e^{(𝜀_{i} - 𝜀_{a}) τ}$	$G_{y z}^{x w} (τ) = ⟨ Ψ_{0} \| a_{x}^{†} (τ) a_{y} (τ) a_{w}^{†} a_{z} \| Ψ_{0} ⟩$
$[+ 1^{'}]$	$𝜖_{i}^{i} (τ) = e^{𝜀_{i} τ}$	$G_{i}^{i} (τ) = ⟨ Ψ_{0} \| {\hat{h}}_{i} (τ) {\hat{h}}_{i}^{†} \| Ψ_{0} ⟩$
$[- 1^{'}]$	$𝜖_{a}^{a} (τ) = e^{- 𝜀_{a} τ}$	$G_{a}^{a} (τ) = ⟨ Ψ_{0} \| {\hat{h}}_{a}^{†} (τ) {\hat{h}}_{a} \| Ψ_{0} ⟩$

The general t-NEVPT2 algorithm consists of three steps: (i) Green’s function evaluation, (ii) computation of the energy as a function of imaginary time (⁠ $τ$ ⁠), and (iii) integration in imaginary time. In step (i), the elements of $G^{[i]} (τ)$ are computed for a set of $τ$ values (time steps, $τ_{k}$ ⁠). In step (ii), the prefactors $𝜖^{[i]} (τ)$ are evaluated using the equations in Table I, and the energy contributions at each time step are computed as $E_{a c t}^{[i]} (τ_{k}) = 𝜖^{[i]} (τ_{k}) G^{[i]} (τ_{k})$ ⁠. For example, $E_{a c t}^{[- 2]} (τ_{k})$ is computed as the dot product of two tensors: $𝜖^{[- 2]} (τ_{k}) = 𝜖_{z w}^{x y} (τ_{k}) = \frac{1}{8} \sum_{a b} v_{z w}^{a b} v_{a b}^{x y} e^{- (𝜀_{a} + 𝜀_{b}) τ_{k}}$ and $G^{[- 2]} (τ_{k}) = G_{z w}^{x y} (τ_{k}) = ⟨ Ψ_{0} | a_{x}^{†} (τ_{k}) a_{y}^{†} (τ_{k}) a_{w} a_{z} | Ψ_{0} ⟩$ ⁠. Finally, in step (iii), each energy contribution $E_{a c t}^{[i]}$ is evaluated by fitting the computed values $E_{a c t}^{[i]} (τ_{k})$ to an exponential function $\sum_{i} E (τ) = a_{i} e^{- b_{i} τ}$ and integrating the obtained result. We note that the above-outlined algorithm can be used to implement t-NEVPT2 for reference wavefunctions in any representation (e.g., determinant-based³⁷ or MPS-based). Only the details of step (i) depend on the explicit form of $| Ψ_{0} (⟩$ ⁠. In Sec. II C 2, we will describe how imaginary-time Green’s functions $G^{[i]} (τ)$ are evaluated in the MPS representation.

2. Green’s function evaluation

To organize our MPS implementation of the imaginary-time Green’s functions, we express the elements of $G^{[i]} (τ)$ (Table I) in the general form $G_{μ ν}^{[i]} (τ) = ⟨ Φ_{μ}^{[i]} (τ) | Φ_{ν}^{[i]} ⟩$ ⁠, where $| Φ_{ν}^{[i]} (⟩$ is obtained by applying all creation and annihilation operators at $τ = 0$ on $| Ψ_{0} (⟩$ (e.g., $a_{x} | Ψ_{0} (⟩$ ⁠, $a_{x}^{†} a_{y} | Ψ_{0} (⟩$ ⁠), while the application of operators at $τ = τ^{'}$ yields $| Φ_{μ}^{[i]} (τ^{'}) (⟩$ (e.g., $a_{x} (τ^{'}) | Ψ_{0} (⟩$ ⁠, $a_{x}^{†} (τ^{'}) a_{y} (τ^{'}) | Ψ_{0} (⟩$ ⁠). In our notation, indices $μ$ and $ν$ run over the total number of states $| Φ_{μ}^{[i]} (⟩$ (i.e., N_act for $a_{x} | Ψ_{0} (⟩$ or $N_{a c t}^{2}$ for $a_{x}^{†} a_{y} | Ψ_{0} (⟩$ ⁠).

We begin the computation of $G^{[i]} (τ)$ by optimizing the zeroth-order MPS wavefunction $| Ψ_{0} (⟩$ using the DMRG algorithm with bond dimension M₀. To evaluate the elements of $G^{[i]} (τ)$ ⁠, we employ the algorithm below.

Compute initial wavefunctions $| Φ_{μ}^{[i]} (⟩$ for every $μ$ ⁠. Each wavefunction $| Φ_{μ}^{[i]} (⟩$ is computed as an individual MPS and stored on disk. Evaluation of $| Φ_{μ}^{[i]} (⟩$ ([i] $\in$ {[+1], [−1], [+2], [−2], $[0^{'}]$ }) requires applying the operator $a_{x}^{†}$ or a_x on $| Ψ_{0} (⟩$ one or more times. The resulting MPS is of exactly the same bond dimension as the original MPS. For $[i] \in {[+ 1^{'}], [- 1^{'}]}$ ⁠, the wavefunctions $| Φ_{μ}^{[+ 1^{'}]} (⟩$ and $| Φ_{μ}^{[- 1^{'}]} (⟩$ are computed by compressing the MPS form of ${\hat{h}}_{i}^{†} | Ψ_{0} (⟩$ and ${\hat{h}}_{a} | Ψ_{0} (⟩$ ⁠, where ${\hat{h}}_{i}^{†}$ and ${\hat{h}}_{a}$ are defined as in Eq. (6). Since applying $ĥ_{i}^{†}$ or $ĥ_{a}$ involves a summation over the active-space labels x, the resulting MPS is of larger bond dimension than $| Ψ_{0} (⟩$ ⁠. We use variational compression to obtain $| Φ_{μ}^{[+ 1^{'}]} (⟩$ and $| Φ_{μ}^{[- 1^{'}]} (⟩$ MPS by minimizing $| | | Φ_{i}^{[+ 1^{'}]} (⟩ - {\hat{h}}_{i}^{†} | Ψ_{0} (⟩ | |$ and $| | | Φ_{a}^{[- 1^{'}]} (⟩ - {\hat{h}}_{a} | Ψ_{0} (⟩ | |$ ⁠, where $| Φ_{i}^{[+ 1^{'}]} (⟩$ and $| Φ_{a}^{[- 1^{'}]} (⟩$ have bond dimension $M_{1} > M_{0}$ ⁠. To maximize the efficiency of computing the overlaps appearing in the compression (e.g., $⟨) Φ_{i}^{[+ 1^{'}]} | {\hat{h}}_{i}^{†} | Ψ_{0} (⟩$ ⁠), we use the standard DMRG formalism of “normal” and “complementary” operators, which requires building $O (k)$ such operators from the left and right blocks. In practice, we find that $M_{1} \approx 2 M_{0}$ usually provides a sufficiently accurate compression. Overall, this step of the t-MPS-NEVPT2 algorithm has $O (N_{a c t}^{2} M_{0}^{3}) + O (N_{e x t} N_{a c t}^{2} M_{1}^{2} M_{0})$ scaling for computing all $| Φ_{μ}^{[i]} (⟩$ wavefunctions with $[i] \in {[+ 2], [- 2], [0^{'}]}$ and $[i] \in {[+ 1^{'}], [- 1^{'}]}$ ⁠, respectively (N_ext is the number of external orbitals).
Propagate wavefunctions $| Φ_{μ}^{[i]} (⟩$ in imaginary time $τ$ according to the time-dependent Schrödinger equation $| Φ_{μ}^{[i]} (τ) (⟩ = e^{- ({\hat{H}}_{D} - E_{D}) τ} | Φ_{μ}^{[i]} (⟩$ ⁠. In this step, the algorithm we use is based on the time step targeting time-dependent DMRG algorithm of Feiguin and White described in Ref. 42, with some small modifications. Specifically, we use the embedded Runge-Kutta [ERK (4,5)] time step algorithm,⁴³ which allows us to automatically control the time step $δ τ$ based on the error estimate of the fifth-order approximation of the propagator $e^{- ({\hat{H}}_{D} - E_{D}) τ}$ ⁠. In the ERK (4,5) method, the wavefunction at time $τ^{'} = τ + δ τ$ is approximated to $O ({(δ τ)}^{5})$ as
$| Φ_{μ}^{[i]} (τ^{'}) (⟩ \approx | Φ_{μ}^{[i]} (τ) (⟩ + \sum_{n = 1}^{6} c_{n} | k_{n} (δ τ) (⟩,$
(17)
where states $| k_{n} (δ τ) (⟩$ are obtained by successively applying the zeroth-order Hamiltonian six times
$| k_{n} (δ τ) (⟩ = - δ τ ({\hat{H}}_{D} - E_{D}) [| Φ_{μ}^{[i]} (τ) (⟩ + \sum_{m = 1}^{n - 1} b_{n m} | k_{m} (δ τ) (⟩] .$
(18)
The coefficients b_nm and c_n are given in Ref. 43. In addition to the fifth-order approximation, Eq. (17) can be used to compute the fourth-order estimate of $| Φ_{μ}^{[i]} (τ^{'}) (⟩$ ⁠, which we denote as $| Φ_{μ}^{[i]} (τ^{'}) [4] (⟩$ ⁠. This can be done by setting c_n = c_n[4], where the fourth-order coefficients c_n[4] can be found in Ref. 43. This gives an estimate of the time step error, $Δ = \sqrt{|) |) | Φ_{μ}^{[i]} (τ^{'}) (⟩ - | Φ_{μ}^{[i]} (τ^{'}) [4] (⟩ |) |)}$ ⁠.
As mentioned in Sec. II B, we can adapt the above general wavefunction propagation to the propagation of a MPS within a sweep algorithm. This is the idea behind the time step targeting td-DMRG. The quantities in the above equations are then to be interpreted as applying to each site, i.e., the wavefunction $| Φ_{μ}^{[i]} (⟩$ is now the site wavefunction, the states $| k_{n} (⟩$ are vectors in the renormalized site basis, and the Hamiltonian ${\hat{H}}_{D}$ is projected into the renormalized site basis. The site error in the propagator approximation is estimated from the site wavefunction as $Δ_{p}$ (where p is the site index). As in Ref. 42, the site wavefunctions $| Φ_{μ}^{[i]} (τ) (⟩$ ⁠, $| Φ_{μ}^{[i]} (τ + δ τ ∕ 3) (⟩$ ⁠, $| Φ_{μ}^{[i]} (τ + 2 δ τ ∕ 3) (⟩$ ⁠, and $| Φ_{μ}^{[i]} (τ + δ τ) (⟩$ are averaged in the density matrix to construct density matrix eigenvectors to update the site wavefunction, before proceeding to the next step in the sweep.
We begin the time propagation with a small value of $δ τ \sim 1 0^{- 3}$ and determine the new time step after each propagation as $δ τ^{'} = min (2 \times δ τ, δ τ_{e m b})$ ⁠, where $δ τ_{e m b} = δ τ {| \frac{Δ E_{c o n v}}{{max}_{p} ({Δ_{p}})} |}^{1 ∕ 5}$ and $Δ E_{c o n v}$ is the specified accuracy threshold. Note that if ${max}_{p} ({Δ_{p}}) > Δ E_{c o n v}$ ⁠, the time step is decreased. In this case, we repeat the time propagation using the smaller time step. The computational cost of a single time step for all states $| Φ_{μ}^{[i]} (⟩$ has $O (N_{a c t}^{5} M_{0}^{3}) + O (N_{a c t}^{6} M_{0}^{2}) + O (N_{e x t} N_{a c t}^{3} M_{1}^{3}) + O (N_{e x t} N_{a c t}^{4} M_{1}^{2})$ scaling.
Compute Green’s function elements at every time step $τ_{k}$ as the overlap between two MPSs $G_{μ ν}^{[i]} (τ_{k}) = ⟨ Φ_{μ}^{[i]} (τ_{k}) | Φ_{ν}^{[i]} ⟩$ ⁠. The computational cost of computing all $G_{μ ν}^{[i]} (τ_{k})$ elements has $O (N_{a c t}^{5} M_{0}^{3})$ scaling.

The cost of the t-MPS-NEVPT2 algorithm outlined above is dominated by the computation of the 2-GF, which requires time-propagation of $O (N_{a c t}^{2})$ MPS wavefunctions $| Φ_{μ}^{[i]} (⟩$ and evaluation of $O (N_{a c t}^{4})$ matrix elements $G_{μ ν}^{[i]} (τ_{k})$ ⁠, leading to a total $O (N_{τ} N_{a c t}^{5} M_{0}^{3}) + O (N_{τ} N_{a c t}^{6} M_{0}^{2})$ computational scaling, where $N_{τ}$ $\sim$ 10-20 is the number of time steps. While this is less than the cost of computing the 4-RDM in DMRG, computing the 3-RDM, which has $O (N_{a c t}^{4} M_{0}^{3}) + O (N_{a c t}^{6} M_{0}^{2})$ scaling, is of lower cost.³⁵ The higher scaling of our 2-GF implementation is because each MPS $| Φ_{μ}^{[i]} (⟩$ is propagated in imaginary time independently, while in the efficient 3-RDM implementation, a common set of renormalized operators is reused to compute all elements of the density matrix together. However, the higher scaling of this particular implementation is not an intrinsic property of the time-dependent theory. In the determinant basis, the computation of the 2-GF and 3-RDM has the same computational scaling. For MPS, it is possible to similarly devise an algorithm to compute the 2-GF with the same cost as the 3-RDM by propagating multiple $| Φ_{μ}^{[i]} (⟩$ wavefunctions using the same renormalized basis (similar to performing a state-averaged DMRG optimization, and related to the algorithm used in Ref. 44). An advantage of our current t-MPS-NEVPT2 implementation, however, is that it is easily parallelized by separating the correlation energy contributions into $O (N_{a c t}^{2})$ components,

E_{a c t}^{[i]} = - \sum_{μ} \int_{0}^{\infty} \sum_{ν} 𝜖_{μ ν}^{[i]} (τ) G_{μ ν}^{[i]} (τ) d τ = \sum_{μ} E_{μ}^{[i]},

(19)

where each component $E_{μ}^{[i]}$ can be evaluated independently with $O (N_{τ} N_{a c t}^{3} M_{0}^{3}) + O (N_{τ} N_{a c t}^{4} M_{0}^{2})$ cost. In addition, as we will demonstrate in Sec. IV, the energy terms that depend on the 2-GF converge very quickly with increasing $M_{0}$ and do not require a large bond dimension.

3. Spin-adaptation

In Sec. II C 2, we discussed the evaluation of the $G_{μ ν}^{[i]} (τ)$ matrix elements in terms of creation and annihilation operators $a_{p}^{†} (τ)$ and $a_{p} (τ)$ ⁠. In a non-spin-adapted DMRG algorithm, $a^{†}$ and a are the spin-orbital operators, and the index p corresponds to a spatial orbital with a particular spin. For example, for [i] = [−2], the spin-orbital 2-GF has the following form: $G_{z_{λ} w_{κ}}^{x_{ρ} y_{σ}} (τ) = ⟨ Ψ_{0} | a_{x_{ρ}}^{†} (τ) a_{y_{σ}}^{†} (τ) a_{w_{κ}} a_{a_{λ}} | Ψ_{0} ⟩$ ⁠, where we use the Roman characters x, y, w, and z to denote spatial orbitals and Greek characters $ρ, σ, κ, λ$ to denote the spin of these orbitals. Here, $a_{w_{κ}}$ and $a_{z_{λ}}$ are simple operators, and application of a pair of operators $a_{w_{κ}} a_{z_{λ}} | Ψ_{0} (⟩$ results in a single MPS wavefunction $| Φ_{w_{κ} z_{λ}}^{[- 2]} (⟩$ ⁠.

In a spin-adapted DMRG implementation,⁴⁵ $a^{†}$ and a refer to spin tensor creation and annihilation operators, and strings of these operators generate multiple eigenstates of the ${\hat{S}}^{2}$ operator when acting on a reference state. As a result, expectation values of spin tensor operators depend on the spin quantum numbers, e.g., $⟨ Ψ_{0} | {{[a_{x}^{†} (τ) a_{y}^{†} (τ)]}_{M_{1}}^{S_{1}} {[a_{w} a_{z}]}_{M_{2}}^{S_{2}}}_{M_{3}}^{S_{3}} | Ψ_{0} ⟩$ ⁠, where we use brackets to denote the coupling of spins S₁, S₂, and S₃ for different tensor products. As an example, we consider the evaluation of $G_{μ ν}^{[- 2]} (τ)$ for a reference state $| Ψ_{0} (⟩$ with S = 0. (Generalization to reference states with $S \neq 0$ is straightforward and is simplified using singlet embedding.)^45,46 Applying a pair of tensor operators on the reference state ${[a_{w} a_{z}]}_{M_{1}}^{S_{1}} | Ψ_{0} (⟩$ results in two sets of eigenstates $| {Φ_{w z}^{[- 2]}}_{M_{1}}^{S_{1}} (⟩$ with $S_{1} = 0,1$ ⁠. The wavefunctions $| {Φ_{w z}^{[- 2]}}_{M_{1}}^{S_{1}} (⟩$ are propagated in imaginary time for each value of S₁, and the matrix elements are computed using Clebsch-Gordan coefficients $c_{M_{1}, M_{2}, M_{3}}^{S_{1}, S_{2}, S_{3}}$ as ${(G_{z w}^{x y})}^{S_{1}} (τ) = \sum_{M_{1}} {(- 1)}^{S_{1} + M_{1}} c_{M_{1}, - M_{1}, 0}^{S_{1}, S_{1}, 0} ⟨ Ψ_{0} | {[a_{x}^{†} (τ) a_{y}^{†} (τ)]}_{M_{1}}^{S_{1}} {[a_{w} a_{z}]}_{M_{1}}^{S_{1}} | Ψ_{0} ⟩$ ⁠, where we used the fact that S₁ = S₂ and S₃ = 0 for a closed-shell $| Ψ_{0} (⟩$ ⁠. Expanding the spin tensor operators ${[a_{w} a_{z}]}^{S_{1}}$ as a combination of spin-orbital operators $a_{w_{ρ}}$ and $a_{z_{σ}}$ leads to a system of linear equations that can be solved to obtain the spin-orbital 2-GF elements $G_{z_{λ} w_{κ}}^{x_{ρ} y_{σ}} (τ)$ from the matrix elements ${(G_{z w}^{x y})}^{S_{1}} (τ)$ (see Appendix A for explicit equations).

III. COMPUTATIONAL DETAILS

All DMRG computations for the reference wavefunctions, reduced density matrices, and Green’s functions were performed using the Block program.⁴⁵ The t-MPS-NEVPT2 and sc-MPS-NEVPT2 correlation energies were computed using Pyscf,⁴⁷ through its interface with Block. For all NEVPT2 computations, the DMRG reference wavefunctions were fully optimized with respect to active-external orbital rotations using the DMRG-SCF algorithm implemented in Pyscf. We denote active spaces used in DMRG-SCF as (ne, mo), where n is the number of active electrons and m is the number of orbitals. We used tight convergence parameters for the energy in the DMRG sweeps (10⁻⁹ $E_{h}$ ⁠) and orbital optimization iterations (10⁻⁶ $E_{h}$ ⁠). In t-MPS-NEVPT2, the accuracy of time propagation was controlled by setting $Δ E_{c o n v} = 1 0^{- 4}$ $E_{h}$ for all systems but the N₂ molecule, where the tighter threshold $Δ E_{c o n v} = 1 0^{- 6}$ $E_{h}$ was used (see Sec. II C 2 for details). Specific details pertaining to our computations of the all-trans polyenes are described in Sec. IV C.

As described in Sec. II C, the bond dimensions necessary to represent the reference MPS (⁠ $M_{0}$ ⁠) and the compressed MPS (⁠ $M_{1}$ ⁠) are usually different (⁠ $M_{0} < M_{1}$ ⁠). However, to simplify the discussion, in our computations, we set M₀ and M₁ to the same value, and whenever we refer to bond dimension M, we imply that $M_{0} = M_{1} = M$ ⁠.

IV. RESULTS

A. Analysis of numerical accuracy: $N_{2}$ dissociation

We first investigate the numerical accuracy of our algorithm by computing the errors of t-MPS-NEVPT2 as a function of the bond dimension relative to the uncontracted NEVPT2 energies computed using our t-NEVPT2/CASSCF implementation.³⁷ For this purpose, we evaluate the t-MPS-NEVPT2 energies for a range of geometries along the ground-state potential energy curve (PEC) of the $N_{2}$ molecule. At each geometry, we optimize the molecular orbitals using the CASSCF method in the (10e, 10o) active space and perform a single DMRG computation using these orbitals to obtain a MPS wavefunction with a large bond dimension $M^{'} = 1000$ ⁠. We subsequently compress this reference wavefunction down to an MPS with a smaller bond dimension (M = 50, 100, 200) to be used in the NEVPT2 computation. In this section, we use Dunning’s cc-pVQZ basis.⁴⁸

Figure 1 plots the t-MPS-NEVPT2 correlation energy errors for different values of the bond dimension M, as well as the error in sc-NEVPT2/CASSCF due to the strong contraction approximation. Notably, even with a small M = 50, the error from finite bond dimension in t-MPS-NEVPT2 is significantly smaller than the contraction error of sc-NEVPT2 for all geometries along the $N_{2}$ potential energy curve (PEC). As the bond dimension M increases, the t-MPS-NEVPT2 errors decrease, becoming smaller than 0.1 $m E_{h}$ at M = 200 for all energy points. As can be seen from Fig. 1, the convergence of t-MPS-NEVPT2 energies is not monotonic along the dissociation curve. In particular, in the range of short bond distances (⁠ $r < 2.1$ Å), the correlation energies converge significantly faster to the exact (uncontracted) NEVPT2 energies than the energies at the dissociation limit (r $\geq$ 2.1 Å). This result is consistent with the fact that both the zeroth-order and first-order wavefunctions at longer bond distances contain contributions from a larger number of determinants than at the shorter bond distances. To illustrate this, we computed the t-MPS-NEVPT2 energies using a large bond dimension M = 1000 (Fig. 1). For all geometries with r $\leq$ 2.1 Å, our t-MPS-NEVPT2 algorithm achieves an accuracy of better than 10⁻⁶ $E_{h}$ ⁠, while errors in the range of $\sim 1 0^{- 5}$ $E_{h}$ still remain in the dissociation limit. The slow convergence of correlation energies with bond dimension near the dissociation limit has also been observed in Sharma and Chan’s MPS-PT2 study of $N_{2}$ dissociation using the cc-pVDZ basis set and the (10e, 8o) active space.⁴⁹

FIG. 1.

View large Download slide

Error in the t-MPS-NEVPT2 correlation energy relative to that of the uncontracted NEVPT2 theory for $N_{2}$ as a function of the N–N bond distance and the DMRG bond dimension M. Computations employed the (10e, 10o) active space and the cc-pVQZ basis set. Also shown are the strong contraction errors computed as the difference between sc-NEVPT2 and t-NEVPT2 based on the CASSCF reference.

We now analyze the errors in the t-MPS-NEVPT2 correlation energy from two different contributions: (i) the sum of active-space energy components that depend on the 1- and 2-GF $\sum_{[i]} E_{a c t}^{[i]}$ (⁠ $[i] \in {[+ 1], [- 1], [+ 2], [- 2], [0^{'}]}$ ⁠) and (ii) the sum of energy terms that depend on the contracted 3-GF (⁠ $E_{a c t}^{[+ 1^{'}]} + E_{a c t}^{[- 1^{'}]}$ ⁠). These errors are plotted in Figs. 2(a) and 2(b), for different values of r and M. Between the two energy contributions, the contribution (i) from the 1- and 2-GF exhibits significantly faster convergence with increasing M than the contribution (ii) from the contracted 3-GF. At M = 50, contribution (i) shows errors of $\sim 1 0^{- 4}$ $E_{h}$ for a large range of geometries, while the errors in contribution (ii) are an order of magnitude larger. Comparing Figs. 2(a) and 2(b), we observe that the energy components (i) and (ii) achieve similar accuracy when computed with bond dimensions of M and 2M, respectively. In fact, for a specified M, the errors in the total t-MPS-NEVPT2 correlation energy are dominated by the errors in contribution (ii) [cf. Figs. 1 and 2(b)]. As discussed in Sec. II C 2, the observed slower convergence of contribution (ii) is the result of the MPS compression used to evaluate $E_{a c t}^{[+ 1^{'}]}$ and $E_{a c t}^{[- 1^{'}]}$ ⁠, which requires a larger bond dimension than the one used in the optimization of the reference MPS for an accurate representation. Since the two contributions can be evaluated separately, we recommend to compute contribution (i) using a small bond dimension M and to use a larger bond dimension $\sim 2 M$ to obtain contribution (ii).

FIG. 2.

View large Download slide

Error in the t-MPS-NEVPT2 correlation energy contributions for $N_{2}$ as a function of the N–N bond distance and the DMRG bond dimension M. Plot (a) shows the errors for $E_{a c t}^{[i]}$ (⁠ $[i] \in {[+ 1], [- 1], [+ 2], [- 2], [0^{'}]}$ ⁠) components that depend on the 1- and 2-GF, while plot (b) presents the errors for $E_{a c t}^{[i]}$ (⁠ $[i] \in {[+ 1^{'}], [- 1^{'}]}$ ⁠) that depend on the contracted 3-GF. Computations employed the (10e, 10o) active space and the cc-pVQZ basis set.

B. Chromium dimer

The chromium dimer (⁠ ${Cr}_{2}$ ⁠) is a notoriously challenging system for multi-reference methods. Computational studies^{29,35,37,50–60} have demonstrated that including the dynamic correlation in combination with the multi-configurational description of static correlation and large basis sets is crucial to properly describe the ground-state potential energy curve (PEC) of this molecule. In 2011, Kurashige and Yanai showed that the ${Cr}_{2}$ dissociation curve can be accurately described by combining the CASPT2 variant of the internally contracted perturbation theory with DMRG in a large (12e, 28o) active space,²⁹ which included the 3d, 4d, 4s, and 4p orbitals of the Cr atoms. Very recently (2016), Guo et al. computed ${Cr}_{2}$ PEC using a combination of sc-NEVPT2 and DMRG (sc-DMRG-NEVPT2) with the (12e, 22o) active space (3d, 4d, and 4s orbitals of Cr).³⁵ At the complete basis set limit (CBS), the ${Cr}_{2}$ PEC from sc-DMRG-NEVPT2 agrees well with the dissociation curve from the experiment, yielding spectroscopic parameters for the equilibrium bond length $r_{e}$ = 1.656 Å and a binding energy $D_{e}$ = 1.432 eV that are close to the experimental $r_{e}^{e x p}$ = 1.679 Å and $D_{e}^{e x p}$ = 1.47 $\pm$ 0.1 eV.⁶¹ However, the ${Cr}_{2}$ binding energy in sc-NEVPT2 can be significantly affected by the strong contraction approximation, as we demonstrated in our t-NEVPT2 study of the ${Cr}_{2}$ PEC using the (12e, 12o) active space.³⁷ In this work, we recompute the ${Cr}_{2}$ equilibrium binding energy in the large (12e, 22o) active space using our t-MPS-NEVPT2 algorithm.

We obtained the MPS reference wavefunction as described in the sc-DMRG-NEVPT2 study by Guo et al.³⁵ In short, we used the (12e, 22o) active space and optimized the molecular orbitals using the DMRG-SCF algorithm with $M^{'}$ = 1000 (no frozen core). We then performed a DMRG computation with $M^{'}$ = 6000 using the optimized orbitals to obtain an accurate reference wavefunction. To compute imaginary-time Green’s functions, the reference MPS with the large bond dimension $M^{'}$ was compressed down to an MPS with a smaller bond dimension M. We computed each t-MPS-NEVPT2 energy contribution $E_{a c t}^{[i]}$ by starting with M = 500 and increasing M until convergence. To account for the scalar relativistic effects, we used the cc-pwCV5Z-DK basis set⁶² combined with the spin-free one-electron variant of the X2C Hamiltonian.^63–65 Since strong contraction reduces the binding energy³⁷ in ${Cr}_{2}$ ⁠, the CBS-extrapolated equilibrium bond length r_e computed using the uncontracted t-MPS-NEVPT2 method is expected to be somewhat shorter than r_e = 1.656 Å obtained from sc-DMRG-NEVPT2. Based on the results of our preliminary study,³⁷ we estimate r_e obtained from the uncontracted theory to be $\sim$ 0.005 Å shorter and use the bond distance r(Cr–Cr) = 1.650 Å in our computations.

Table II presents the t-MPS-NEVPT2 correlation energy contributions computed for different values of the bond dimension M. Increasing M from 500 to 800 does not significantly affect the energy components $E_{a c t}^{[i]}$ that depend on the 1-GF (⁠ $[i] \in {[+ 1], [- 1]}$ ⁠) and 2-GF (⁠ $[i] \in {[+ 2], [- 2], [0^{'}]}$ ⁠). Out of five of these contributions, the largest change of +0.3 $m E_{h}$ is observed for $E_{a c t}^{[0^{'}]}$ ⁠, while each of the other four energy terms changes by less than −0.1 $m E_{h}$ ⁠. These changes largely cancel each other, leading to a total energy difference of $\sim$ 0.1 $m E_{h}$ ⁠. Thus, we estimate the total error in the sum of the 1- and 2-GF energy contributions to be less than 0.1 $m E_{h}$ ⁠. However, the remaining energy components $E_{a c t}^{[+ 1^{'}]}$ and $E_{a c t}^{[- 1^{'}]}$ that depend on the contracted 3-GF exhibit slower convergence with M, changing by −1.9 $m E_{h}$ for M = 500 $\to$ 800 and −1.2 $m E_{h}$ for M = 800 $\to$ 1200 (Table II). We increased the bond dimension up to M = 2000, where the remaining error in $E_{a c t}^{[+ 1^{'}]}$ and $E_{a c t}^{[- 1^{'}]}$ is estimated to be less than 0.4 $m E_{h}$ (⁠ $<$ 0.01 eV). The resulting active-space energy contributions $E_{a c t}^{[i]}$ can be used to compute the ${Cr}_{2}$ binding energy for the incomplete cc-pwCV5Z-DK basis set as $D_{e}$ = $E_{e x t}^{(2)} + \sum_{[i]} E_{a c t}^{[i]}$ − $2 \times E_{t - N E V P T 2}^{(2)} (C r)$ ⁠, where $E_{e x t}^{(2)}$ is the ${Cr}_{2}$ correlation energy from the external space [Eq. (15)] and $E_{t - N E V P T 2}^{(2)} (C r)$ is the correlation energy for an isolated Cr atom. To estimate $D_{e}$ at the CBS limit, we first compute the error of strong contraction in the binding energy $Δ E_{c o n t r} (D_{e})$ as shown in Table II, where $E_{s c - N E V P T 2}^{(2)}$ is the sc-DMRG-NEVPT2 correlation energy for ${Cr}_{2}$ reported by Guo et al.,³⁵ while $E_{s c - N E V P T 2}^{(2)} (C r)$ and $E_{t - N E V P T 2}^{(2)} (C r)$ are the correlation energies for a Cr atom computed using the (6e, 11o) CASSCF reference. Assuming that the contraction error does not change significantly from the incomplete (cc-pwCV5Z-DK) to the complete basis set, the CBS-extrapolated ${Cr}_{2}$ binding energy at the t-MPS-NEVPT2 level of theory can be estimated as ${(D_{e})}_{t - N E V P T 2}^{C B S} \approx {(D_{e})}_{s c - N E V P T 2}^{C B S} + Δ E_{c o n t r} (D_{e})$ ⁠, where we use ${(D_{e})}_{s c - N E V P T 2}^{C B S}$ = 1.43 eV from Ref. 35. Although the resulting binding energy ${(D_{e})}_{t - N E V P T 2}^{C B S}$ = 1.52 eV is 0.09 eV larger than ${(D_{e})}_{s c - N E V P T 2}^{C B S}$ ⁠, it is still in a good agreement with the experimental binding energy of 1.47 $\pm$ 0.1 eV reported by Casey and Leopold.⁶¹

TABLE II.

t-MPS-NEVPT2 correlation energy contributions $E_{a c t}^{[i]}$ for Cr₂ at 1.650 Å computed for different values of the bond dimension M. Computations employed the (12e, 22o) active space and the cc-pwCV5Z-DK basis set. Also shown are strong contraction errors $Δ E_{c o n t r}$ in the correlation energy for Cr₂, the Cr atom, and the Cr₂ binding energy (⁠ $D_{e}$ ⁠), evaluated as the difference between the sc-DMRG-NEVPT2 and t-MPS-NEVPT2 energies. The error of strong contraction in the binding energy is shown in bold.

M	$E_{a c t}^{[+ 1]}$	$E_{a c t}^{[- 1]}$	$E_{a c t}^{[+ 2]}$	$E_{a c t}^{[- 2]}$	$E_{a c t}^{[0^{'}]}$	$E_{a c t}^{[+ 1^{'}]}$	$E_{a c t}^{[- 1^{'}]}$
500	−0.1632	−0.2357	−0.0939	−0.1140	−1.6425	−0.0282	−0.0416
800	−0.1633	−0.2357	−0.0940	−0.1140	−1.6422	−0.0286	−0.0431
1200						−0.0288	−0.0441
1600						−0.0290	−0.0446
2000						−0.0291	−0.0449
$Δ E_{c o n t r} ({C r}_{2}) = E_{s c - N E V P T 2}^{(2)} - E_{e x t}^{(2)} - \sum_{[i]} E_{a c t}^{[i]}$ = −1.6351 − 0.6801 + 2.3232 = 0.0080 $E_{h}$
$Δ E_{c o n t r} (C r) = E_{s c - N E V P T 2}^{(2)} (C r) - E_{t - N E V P T 2}^{(2)} (C r)$ = −0.7758 + 0.7781 = 0.0023 $E_{h}$
$Δ E_{c o n t r} (D_{e})$ = $Δ E_{c o n t r} ({C r}_{2})$ − $2 \times Δ E_{c o n t r} (C r)$ = 0.0034 $E_{h}$ = 0.09 eV

M	$E_{a c t}^{[+ 1]}$	$E_{a c t}^{[- 1]}$	$E_{a c t}^{[+ 2]}$	$E_{a c t}^{[- 2]}$	$E_{a c t}^{[0^{'}]}$	$E_{a c t}^{[+ 1^{'}]}$	$E_{a c t}^{[- 1^{'}]}$
500	−0.1632	−0.2357	−0.0939	−0.1140	−1.6425	−0.0282	−0.0416
800	−0.1633	−0.2357	−0.0940	−0.1140	−1.6422	−0.0286	−0.0431
1200						−0.0288	−0.0441
1600						−0.0290	−0.0446
2000						−0.0291	−0.0449
$Δ E_{c o n t r} ({C r}_{2}) = E_{s c - N E V P T 2}^{(2)} - E_{e x t}^{(2)} - \sum_{[i]} E_{a c t}^{[i]}$ = −1.6351 − 0.6801 + 2.3232 = 0.0080 $E_{h}$
$Δ E_{c o n t r} (C r) = E_{s c - N E V P T 2}^{(2)} (C r) - E_{t - N E V P T 2}^{(2)} (C r)$ = −0.7758 + 0.7781 = 0.0023 $E_{h}$
$Δ E_{c o n t r} (D_{e})$ = $Δ E_{c o n t r} ({C r}_{2})$ − $2 \times Δ E_{c o n t r} (C r)$ = 0.0034 $E_{h}$ = 0.09 eV

C. Singlet excited states in all-trans polyenes

1. Background

Conjugated polyenes are important model systems for understanding properties of organic materials (such as polyacetylene), as well as pigments involved in photosynthesis and vision (e.g., carotenoids or retinals). The key to the functionality of these molecules is the ability to absorb and efficiently transfer energy via the low-lying $π – π^{*}$ excited states of the conjugated backbone. One of the simplest examples of polyene systems is the unsubstituted all-trans polyenes that consist of alternating single and double carbon bonds arranged in a linear chain (⁠ $C_{n} H_{n + 2}$ ⁠). The excited states of these molecules are labeled by the symmetry of the $C_{2 h}$ point group (⁠ $A_{g}$ ⁠, $A_{u}$ ⁠, $B_{g}$ ⁠, and $B_{u}$ ⁠), as well as an additional +/− label due to their approximate particle-hole symmetry, which is typically used to distinguish excited states with mainly ionic (+ states) or mainly covalent (− states) excitation character.^66,67

The electronic spectra of all-trans polyenes have been the subject of many experimental and computational studies.^39,67–85 Of particular interest are the excitation energies and ordering of the low-lying singlet electronic states. While the symmetry of the singlet ground state is known to always be $A_{g}^{-}$ ⁠, the experimental and theoretical characterization of the low-lying singlet excited states has been very challenging. In shorter polyenes (⁠ $C_{n} H_{n + 2}$ ⁠, n $\leq$ 8), it has been established that the two lowest-energy singlet transitions are the dipole-allowed excitation to the ionic $1 B_{u}^{+}$ state and the dipole-forbidden excitation to the $2 A_{g}^{-}$ state. Despite the fact that the allowed transition $1 A_{g}^{-}$ $\to$ $1 B_{u}^{+}$ has primarily HOMO $\to$ LUMO excitation character, fluorescence studies suggest that, starting from octatetraene (⁠ $C_{8} H_{10}$ ⁠),⁷⁵ the lowest-energy singlet excitation corresponds to a forbidden transition, such as the $1 A_{g}^{-}$ $\to$ $2 A_{g}^{-}$ transition. Recent high-level theoretical studies of butadiene (⁠ $C_{4} H_{6}$ ⁠) predict the bright $1 B_{u}^{+}$ state vertical excitation to be $\sim$ 0.2 eV lower in energy than the excitation to the dark $2 A_{g}^{-}$ state,⁷³ whereas the two vertical transitions are predicted to be almost degenerate in $C_{8} H_{10}$ ⁠.⁷² In longer polyenes, the low-energy electronic spectrum is more complicated, involving additional dark states $1 B_{u}^{-}$ and $3 A_{g}^{-}$ with energies close to $1 B_{u}^{+}$ and $2 A_{g}^{-}$ ⁠. These dark states have been observed experimentally in all-trans-carotenoids with n = 18–26 $π$ -electrons.^76–78

Early computational studies of long polyenes have been primarily based on model Hamiltonians, such as the Pariser-Parr-Pople (PPP) and Hubbard models.^67,68 In 1986, Tavan and Schulten⁶⁷ performed multi-reference configuration interaction (MRCI) computations using the PPP model, which suggested that at n = 10, the $1 B_{u}^{-}$ state becomes the second excited state, leading to the $1 A_{g}^{-}$ < $2 A_{g}^{-}$ < $1 B_{u}^{-}$ < $1 B_{u}^{+}$ < $3 A_{g}^{-}$ ordering of states for longer polyenes. These computational results were recently reassessed using a combination of MRCI and the semi-empirical OM2 method,⁸⁶ which predicted the inversion of the $1 B_{u}^{-}$ and $1 B_{u}^{+}$ transitions at n = 14.

In 1998, Hirao and co-workers performed the first ab initio computations of the vertical excitation energies in small and medium-size polyenes (n $\leq$ 10) using multi-reference perturbation theory (MRMP) based on the CASSCF reference wavefunctions.⁶⁹ In a later study, Kurashige et al. used a combination of MRMP and CASCI with the (10e, 10o) active space (incomplete $π$ -valence space) for polyene series up to n = 28.⁷⁰ Their computations also suggested that the $1 B_{u}^{-}$ and $1 B_{u}^{+}$ vertical transitions cross at n = 14 and predicted a crossing of $3 A_{g}^{-}$ and $1 B_{u}^{+}$ vertical excitations at n = 22. In 2008, Ghosh et al. performed a DMRG-SCF study of polyenes with 8 $\leq$ n $\leq$ 24 using the complete $π$ -valence active space (ne, no).⁷¹ Although in this study the ionic $1 B_{u}^{+}$ state was not investigated, the authors demonstrated that the self-consistent optimization of orbitals in the DMRG computations of these systems is very important in achieving reasonable agreement with the experiment for the covalent $2 A_{g}^{-}$ ⁠, $1 B_{u}^{-}$ ⁠, and $3 A_{g}^{-}$ states. The work by Ghosh et al., however, did not incorporate dynamic correlation. Although it is generally believed that the dynamic correlation mainly affects the ionic $1 B_{u}^{+}$ state, excitation energies of the covalent states computed using DMRG-SCF overestimate those obtained from the experiment.⁷¹ Using our sc-MPS-NEVPT2 and t-MPS-NEVPT2 methods, we are now in the position to investigate the importance of dynamic correlation effects on the electronic excitations in all-trans polyenes in combination with the complete $π$ -valence space description of static correlation.

2. Details of computations

The ground-state molecular geometries of the all-trans polyenes $C_{n} H_{n + 2}$ (n = 4, 8, 16, and 24) were optimized using the B3LYP method^87,88 implemented in the Psi4 program.⁸⁹ Geometry optimizations were constrained to have $C_{2 h}$ symmetry. For all polyenes, we used the cc-pVDZ basis set.⁴⁸ In addition, to study the effect of the basis set on the excitation energies, we also employed the larger aug-cc-pVTZ basis for the smaller polyenes (n = 4 and 8).

At each optimized geometry, the reference MPS wavefunctions were computed using the DMRG-SCF algorithm for the complete $π$ -valence active space (ne, no). To construct this active space, we first performed a self-consistent field (SCF) computation to obtain a set of canonical Hartree-Fock orbitals. The occupied subset of the SCF orbitals was used to generate an orthonormal set of intrinsic atomic orbitals (IAOs),⁹⁰ from which the $2 p_{z}$ orbitals of carbon atoms (z being the C₂ axis) were easily identified and used as the initial active space orbitals. The initial core and external orbitals were obtained by projecting the active-space orbitals out of the occupied and virtual subspaces of the SCF orbitals, respectively, followed by their symmetric orthogonalization. These initial orbitals were subsequently optimized using the state-averaged DMRG-SCF algorithm with $M^{'} = 250$ ⁠, averaging over the five lowest-energy eigenstates with equal weights. At the end of the DMRG-SCF computations, the active orbitals were relocalized using the Pipek-Mezey procedure,⁹¹ and the core and external orbitals were canonicalized as the eigenvalues of the state-specific generalized Fock operator [Eq. (2)]. To compute the t-MPS-NEVPT2 and sc-MPS-NEVPT2 energies, reference MPS for each eigenstate was first computed with $M^{'} = 500$ and then compressed down to M = 250. For $C_{n} H_{n + 2}$ with n = 4, 8, and 16, the remaining error in the reference and correlation energies was estimated to be less than 0.1 $m E_{h}$ ⁠. To achieve a similar accuracy for $C_{24} H_{26}$ ⁠, we used M = 500 in sc-MPS-NEVPT2 and to compute the $E_{a c t}^{[+ 1^{'}]}$ and $E_{a c t}^{[- 1^{'}]}$ contributions in t-MPS-NEVPT2.

Finally, to assign the spatial symmetries of the excited states, we computed transition dipole moments between the states. Forbidden transitions were assigned to have A_g symmetry, while the allowed transitions were labeled as B_u. The B_u states corresponding to a large transition dipole moment were assigned as $B_{u}^{+}$ ⁠, indicating that the electronic excitation involves a change of particle-hole character, whereas the remaining B_u states were labeled as $B_{u}^{-}$ ⁠. For polyenes with n = 4, 8, 16, and 24, the $1 B_{u}^{+}$ state was found to be the 3rd, 5th, 7th, and 9th singlet eigenstate of the zeroth-order Hamiltonian, respectively. Thus, to compute the $1 B_{u}^{+}$ excitation energies, we increased the number of states evaluated in the state-averaged DMRG computation with $M^{'} = 500$ to 9 and 11 for n = 16 and 24, respectively.

3. Discussion

Table III presents energies, symmetries, and oscillator strengths for the low-lying singlet states of polyenes $C_{n} H_{n + 2}$ (n = 4, 8, 16, and 24) computed using the DMRG-SCF, sc-MPS-NEVPT2, and t-MPS-NEVPT2 methods with the cc-pVDZ basis set. In addition, Fig. 3 plots the relative energies of the excited states obtained from DMRG-SCF and t-MPS-NEVPT2 as a function of the chain length n. For polyenes with n = 8, 16, and 24, our DMRG-SCF computations predict the energies and ordering of the covalent excited states $2 A_{g}^{-}$ < $1 B_{u}^{-}$ < $3 A_{g}^{-}$ in close agreement with the results obtained by Ghosh et al.⁷¹ In particular, for all polyenes (including $C_{4} H_{6}$ ⁠), DMRG-SCF predicts the $2 A_{g}^{-}$ state to be the lowest-energy singlet excited state, while the energy of the ionic $1 B_{u}^{+}$ state dominated by the HOMO $\to$ LUMO transition is computed to be 1.6 eV to 2.3 eV higher.

TABLE III.

Energies, symmetries, and oscillator strengths for the low-lying singlet states in the all-trans conjugated polyenes $C_{n} H_{n + 2}$ ⁠. Entries for the $1 A_{g}^{-}$ ground states give the total energies in $E_{h}$ computed using DMRG-SCF, sc-MPS-NEVPT2, and t-MPS-NEVPT2 with the cc-pVDZ basis set. Entries for the excited states give vertical excitation energies from the ground state in eV. The notation (ne, mo) denotes the active spaces used in the DMRG and NEVPT2 computations. Oscillator strengths in a.u. were computed using DMRG-SCF. The reference numbers in brackets are experimental measurements for substituted polyenes.

Polyene	State	DMRG-SCF	sc-MPS-NEVPT2	t-MPS-NEVPT2	Oscillator strength	Reference
$C_{4} H_{6}$	$1 A_{g}^{-}$	−154.9772	−155.4862	−155.4866
(4e, 4o)	$1 B_{u}^{+}$	8.29	6.27	6.17	2.251	5.92^a, 6.21^b
	$2 A_{g}^{-}$	6.67	6.96	6.94	Forbidden	6.41^b
$C_{8} H_{10}$	$1 A_{g}^{-}$	−308.8259	−309.8154	−309.8168
(8e, 8o)	$1 B_{u}^{+}$	6.74	4.24	4.12	3.345	4.41^c, 4.76^d
	$2 A_{g}^{-}$	4.71	4.77	4.75	Forbidden	3.59^e, 4.81^d
	$1 B_{u}^{-}$	5.91	6.03	6.00	0.056	5.96^d
	$3 A_{g}^{-}$	6.64	6.80	6.76	Forbidden
$C_{16} H_{18}$	$1 A_{g}^{-}$	−616.5142	−618.4719	−618.4754
(16e, 16o)	$1 B_{u}^{+}$	5.51	2.82	2.64	4.964	(2.82)^f
	$2 A_{g}^{-}$	3.24	3.14	3.09	Forbidden	(2.21)^f
	$1 B_{u}^{-}$	4.02	3.95	3.89	0.031
	$3 A_{g}^{-}$	4.77	4.73	4.65	Forbidden
$C_{24} H_{26}$	$1 A_{g}^{-}$	−924.2014	−927.1284	−927.1354
(24e, 24o)	$1 B_{u}^{+}$	5.04	2.25	2.00	6.167	(2.25)^g
	$2 A_{g}^{-}$	2.72	2.54	2.45	Forbidden	(1.51)^g
	$1 B_{u}^{-}$	3.24	3.08	2.94	0.025	(1.80)^g
	$3 A_{g}^{-}$	3.77	3.63	3.49	Forbidden	(2.04)^g

Polyene	State	DMRG-SCF	sc-MPS-NEVPT2	t-MPS-NEVPT2	Oscillator strength	Reference
$C_{4} H_{6}$	$1 A_{g}^{-}$	−154.9772	−155.4862	−155.4866
(4e, 4o)	$1 B_{u}^{+}$	8.29	6.27	6.17	2.251	5.92^a, 6.21^b
	$2 A_{g}^{-}$	6.67	6.96	6.94	Forbidden	6.41^b
$C_{8} H_{10}$	$1 A_{g}^{-}$	−308.8259	−309.8154	−309.8168
(8e, 8o)	$1 B_{u}^{+}$	6.74	4.24	4.12	3.345	4.41^c, 4.76^d
	$2 A_{g}^{-}$	4.71	4.77	4.75	Forbidden	3.59^e, 4.81^d
	$1 B_{u}^{-}$	5.91	6.03	6.00	0.056	5.96^d
	$3 A_{g}^{-}$	6.64	6.80	6.76	Forbidden
$C_{16} H_{18}$	$1 A_{g}^{-}$	−616.5142	−618.4719	−618.4754
(16e, 16o)	$1 B_{u}^{+}$	5.51	2.82	2.64	4.964	(2.82)^f
	$2 A_{g}^{-}$	3.24	3.14	3.09	Forbidden	(2.21)^f
	$1 B_{u}^{-}$	4.02	3.95	3.89	0.031
	$3 A_{g}^{-}$	4.77	4.73	4.65	Forbidden
$C_{24} H_{26}$	$1 A_{g}^{-}$	−924.2014	−927.1284	−927.1354
(24e, 24o)	$1 B_{u}^{+}$	5.04	2.25	2.00	6.167	(2.25)^g
	$2 A_{g}^{-}$	2.72	2.54	2.45	Forbidden	(1.51)^g
	$1 B_{u}^{-}$	3.24	3.08	2.94	0.025	(1.80)^g
	$3 A_{g}^{-}$	3.77	3.63	3.49	Forbidden	(2.04)^g

^a

Experimental adiabatic excitation energy.^79,80

^b

Best theoretical vertical excitation energy.⁷³

^c

Experimental adiabatic excitation energy.^81–83

^d

Best theoretical vertical excitation energy.⁷²

^e

Experimental adiabatic excitation energy.⁸⁴

^f

Experimental adiabatic excitation energy for a substituted polyene.⁸⁵

^g

Experimental adiabatic excitation energy for a substituted polyene.⁷⁸

FIG. 3.

View large Download slide

Excitation energies in eV for the $1 B_{u}^{+}$ ⁠, $2 A_{g}^{-}$ ⁠, $1 B_{u}^{-}$ ⁠, and $3 A_{g}^{-}$ singlet states in the all-trans conjugated polyenes $C_{n} H_{n + 2}$ computed using DMRG-SCF and t-MPS-NEVPT2 with the cc-pVDZ basis set. Computations employed (ne, no) active spaces, where the number of active electrons and orbitals n is equal to the number of C atoms.

Including the dynamic correlation lowers the relative energy of the ionic 1 $B_{u}^{+}$ state drastically. This is seen in the difference of the 1 $B_{u}^{+}$ excitation energies computed using sc-MPS-NEVPT2 and DMRG-SCF, which ranges from 2 eV in $C_{4} H_{6}$ to 2.7 eV in $C_{24} H_{26}$ ⁠, reducing the excitation energy for the longer polyene by more than a factor of two. Similar results are obtained using the uncontracted t-MPS-NEVPT2 method, which lowers the $1 B_{u}^{+}$ energy by an additional 0.10 to 0.25 eV relative to sc-MPS-NEVPT2, due to the removal of the strong contraction approximation. Both methods predict $1 B_{u}^{+}$ to be the lowest-energy singlet excited state. The effect of dynamic correlation on the relative energies of the covalent states $2 A_{g}^{-}$ ⁠, $1 B_{u}^{-}$ ⁠, and $3 A_{g}^{-}$ is much smaller (Fig. 3). Nevertheless, the correlation contributions to the relative energies of these states increase with chain length and become significant for longer polyenes, reaching $\sim$ 0.3 eV in $C_{24} H_{26}$ at the t-MPS-NEVPT2 level of theory. Interestingly, in shorter polyenes (⁠ $C_{4} H_{6}$ and $C_{8} H_{10}$ ⁠), including the dynamic correlation actually leads to an increase of the excitation energies for the covalent states computed using sc-MPS-NEVPT2 and t-MPS-NEVPT2 relative to DMRG-SCF, indicating that the correlation is more important in the ground rather than the excited electronic states for these molecules. As we describe below, a similar trend is observed even when the larger aug-cc-pVTZ basis set is used in the computations. For longer polyenes, dynamic correlation lowers the vertical excitation energies for the covalent transitions. As a result, the sc-MPS-NEVPT2 and t-MPS-NEVPT2 excitation energies decrease with increasing chain length n faster than compared to DMRG-SCF, as shown in Fig. 3.

Figure 4 plots the error of the strong contraction approximation in the excitation energies of polyenes with the number of carbon atoms computed as the difference between the sc-MPS-NEVPT2 and t-MPS-NEVPT2 results. Although in the shorter polyenes (⁠ $C_{4} H_{6}$ and $C_{8} H_{10}$ ⁠), strong contraction mainly affects the energy of the ionic $1 B_{u}^{+}$ state, in the longer polyenes, the errors of strong contraction become significant even for the covalent $2 A_{g}^{-}$ ⁠, $1 B_{u}^{-}$ ⁠, and $3 A_{g}^{-}$ states. In particular, for $C_{24} H_{26}$ ⁠, the strong contraction errors amount to 0.09–0.14 eV in the excitation energy for the covalent transitions, which is comparable to the $\sim$ 0.15 eV lowering of the energies due to the dynamic correlation computed by sc-MPS-NEVPT2.

FIG. 4.

View large Download slide

Strong contraction errors in excitation energies (eV) for the lowest-lying singlet states of the all-trans conjugated polyenes $C_{n} H_{n + 2}$ computed as the difference between the sc-MPS-NEVPT2 and t-MPS-NEVPT2 results.

To assess the relative performance of sc-MPS-NEVPT2 and t-MPS-NEVPT2 for the description of excited states with ionic and covalent character, we computed the $C_{4} H_{6}$ and $C_{8} H_{10}$ vertical excitation energies using sc-NEVPT2 and t-NEVPT2 with the large aug-cc-pVTZ basis set and CASSCF reference wavefunctions (Table IV). Increasing the basis set from cc-pVDZ to aug-cc-pVTZ lowers the t-NEVPT2 relative energies of the $1 B_{u}^{+}$ and $2 A_{g}^{-}$ states in $C_{4} H_{6}$ by 0.34 and 0.32 eV, respectively. For the covalent $2 A_{g}^{-}$ state, both sc-NEVPT2 and t-NEVPT2 predict excitation energies of $\sim$ 6.6 eV, in a reasonable agreement with a reference value of 6.41 eV from highly accurate coupled cluster computations.⁷³ Considering the strong basis set dependence of the $2 A_{g}^{-}$ excitation energy, we expect that this agreement will improve at the CBS limit. On the other hand, the sc-NEVPT2 and t-NEVPT2 vertical excitation energies for the $1 B_{u}^{+}$ state (6.06 and 5.83 eV) are too low relative to the best available theoretical value of 6.21 eV (Table IV), indicating that the two methods will significantly underestimate the excitation energy of the ionic state at the CBS limit. A similar performance of sc-NEVPT2 and t-NEVPT2 is observed for $C_{8} H_{10}$ ⁠. In this case, both methods yield excitation energies for the covalent states $2 A_{g}^{-}$ and $1 B_{u}^{-}$ that are in the excellent agreement with the best available theoretical reference values,⁷² while the relative energy of the ionic $1 B_{u}^{+}$ state is underestimated by more than 0.7 eV. Similar errors for the $1 B_{u}^{+}$ state have been observed by Angeli and Pastore in the study of $C_{8} H_{10}$ using the second-order quasi-degenerate perturbation theory (QD-PT2) with the (8e, 8o) active space.⁷² They demonstrated that increasing the active space from (8e, 8o) to (8e, 16o) results in $\sim$ 0.7 eV increase in the $1 B_{u}^{+}$ excitation energy computed using QD-PT2, leading to a closer agreement with the experiment.

TABLE IV.

Energies for the lowest-lying singlet states in the all-trans conjugated polyenes $C_{4} H_{6}$ and $C_{8} H_{10}$ ⁠. Entries for the $1 A_{g}^{-}$ ground states give the total energies in $E_{h}$ computed using CASSCF, sc-NEVPT2, and t-NEVPT2 with the aug-cc-pVTZ basis set. Entries for the excited states give vertical excitation energies from the ground state in eV. The notation (ne, mo) denotes the active spaces used in the CASSCF computations.

Polyene	State	CASSCF	sc-NEVPT2	t-NEVPT2	Reference
$C_{4} H_{6}$	$1 A_{g}^{-}$	−154.9845	−155.7199	−155.7211
(4e, 4o)	$1 B_{u}^{+}$	7.33	6.06	5.83	5.92^a, 6.21^b
	$2 A_{g}^{-}$	5.40	6.61	6.62	6.41^b
$C_{8} H_{10}$	$1 A_{g}^{-}$	−308.9070	−310.2677	−310.2694
(8e, 8o)	$1 B_{u}^{+}$	6.65	4.02	3.79	4.41^c, 4.76^d
	$2 A_{g}^{-}$	4.79	4.81	4.78	3.59^e, 4.81^d
	$1 B_{u}^{-}$	6.00	6.07	6.02	5.96^d
	$3 A_{g}^{-}$	6.74	6.82	6.75

Polyene	State	CASSCF	sc-NEVPT2	t-NEVPT2	Reference
$C_{4} H_{6}$	$1 A_{g}^{-}$	−154.9845	−155.7199	−155.7211
(4e, 4o)	$1 B_{u}^{+}$	7.33	6.06	5.83	5.92^a, 6.21^b
	$2 A_{g}^{-}$	5.40	6.61	6.62	6.41^b
$C_{8} H_{10}$	$1 A_{g}^{-}$	−308.9070	−310.2677	−310.2694
(8e, 8o)	$1 B_{u}^{+}$	6.65	4.02	3.79	4.41^c, 4.76^d
	$2 A_{g}^{-}$	4.79	4.81	4.78	3.59^e, 4.81^d
	$1 B_{u}^{-}$	6.00	6.07	6.02	5.96^d
	$3 A_{g}^{-}$	6.74	6.82	6.75

^a

Experimental adiabatic excitation energy.^79,80

^b

Best theoretical vertical excitation energy.⁷³

^c

Experimental adiabatic excitation energy.^81–83

^d

Best theoretical vertical excitation energy.⁷²

^e

Experimental adiabatic excitation energy.⁸⁴

The underestimation of the $1 B_{u}^{+}$ relative energies due to insufficiently large active spaces may explain the dependence of the excitation energies with the chain length n computed using t-MPS-NEVPT2 (Fig. 3), where we observe no crossing of $1 B_{u}^{+}$ with the covalent $2 A_{g}^{-}$ and $1 B_{u}^{-}$ transitions, contrary to the experimental results on substituted polyenes^78,85 and some of the calculations.^70,86 Re-investigating the trends in the excitation energies of the long polyenes with a double $π$ -active space, as employed by Angeli and Pastore, will be the subject of future work.

V. CONCLUSIONS

In this work, we developed an efficient algorithm to describe dynamic correlation in strongly correlated systems with large active spaces and basis sets based on a combination of the time-dependent second-order N-electron valence perturbation theory (t-NEVPT2) with matrix product state (MPS) reference wavefunctions. The resulting t-MPS-NEVPT2 approach is equivalent to the fully uncontracted N-electron valence perturbation theory but can efficiently compute correlation energies using the time-dependent density matrix renormalization group algorithm. In addition, we presented a new MPS-based implementation of strongly contracted NEVPT2 (sc-MPS-NEVPT2) that has a lower computational scaling than the commonly used internally contracted NEVPT2 variants. Using these new methods, we computed the dissociation energy of the chromium dimer and investigated the importance of dynamic correlation in the low-lying excited states of all-trans polyenes (⁠ $C_{4} H_{6}$ to $C_{24} H_{26}$ ⁠). For the chromium dimer, the active space used included the “double d” shell of the chromium atoms. Our ${Cr}_{2}$ dissociation energy computed using the uncontracted t-MPS-NEVPT2 method is in a good agreement with the experimental data and the previously reported results from strongly contracted NEVPT2.³⁵ In a study of polyenes, using a complete $π$ -valence active space, our results demonstrate that including dynamic correlation is very important for the ionic $π – π^{*}$ transitions, while it is less important for excitations with covalent character. In particular, for the $C_{24} H_{26}$ polyene, incorporating dynamic correlation at the t-MPS-NEVPT2 level of theory lowers the energy of the ionic transition by $\sim$ 3 eV, whereas only $\sim$ 0.3 eV change is observed for the covalent electronic states. Our results suggest that both sc-MPS-NEVPT2 and t-MPS-NEVPT2 significantly underestimate the relative energy of the ionic state when combined with the complete $π$ -valence active space. In this case, using the “double” $π$ -active space is expected to improve the description of the excitation energies and will be the subject of our future work. Overall, our results demonstrate that t-MPS-NEVPT2 and sc-MPS-NEVPT2 are promising methods for the study of strongly correlated systems that can be applied to a variety of challenging problems.

ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy (DOE), Office of Science through Award No. DE-SC0008624 (primary support for A.Y.S.) and Award No. DE-SC0010530 (additional support for G.K.-L.C.), and used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The authors would like to thank Dr. Zhendong Li for insightful discussions.

APPENDIX A: SPIN-ORBITAL TWO-BODY GREEN’S FUNCTIONS IN SPIN-ADAPTED DMRG

Here , we describe how to compute the spin-orbital 2-GF using the spin-adapted DMRG algorithm. First, we evaluate the elements of $G^{[- 2]} (τ)$ ⁠, which in the spin-orbital basis have the following form: $G_{z_{κ} w_{λ}}^{x_{ρ} y_{σ}} (τ) = ⟨) Ψ_{0} | a_{x ρ}^{†} (τ) a_{y σ}^{†} (τ) a_{w λ} a_{z κ} | Ψ_{0} ⟩) \equiv ⟨) a_{x ρ}^{†} (τ) a_{y σ}^{†} (τ) a_{w λ} a_{z κ} ⟩)$ ⁠, where x, y, w, z denote the spatial orbitals and $ρ, σ, κ, λ$ denote the spin of these orbitals (e.g., $ρ = α, β$ ⁠). As discussed in Sec. II C 3, in spin-adapted DMRG, $a_{x}^{†}$ and $a_{x}$ are the spin tensor creation and annihilation operators. We use a shorthand notation for the product of two tensor operators ${\hat{A}}_{x y}^{S, M} \equiv {[a_{x}^{†} a_{y}^{†}]}_{M}^{S}$ ⁠, where square brackets denote spin coupling (S = 0, 1). The adjoint of a spin tensor operator is defined with an additional sign factor to preserve the Condon-Shortley phase convention: ${\hat{A}}_{x y}^{S, M ‡} = {(- 1)}^{S + M} {\hat{A}}_{x y}^{S, - M †}$ ⁠. In addition, we define the spin-adapted tensor product $\otimes_{S}$ as $Â_{1}^{S_{1}} \otimes_{S} Â_{2}^{S_{2}} = {[{\hat{A}}_{1}^{S_{1}} {\hat{A}}_{2}^{S_{2}}]}^{S}$ ⁠, where ${[{\hat{A}}_{1}^{S_{1}} {\hat{A}}_{2}^{S_{2}}]}_{M}^{S} = \sum_{M_{1} M_{2}} c_{M_{1}, M_{2}, M}^{S_{1}, S_{2}, S} {({\hat{A}}_{1})}_{M_{1}}^{S_{1}} {({\hat{A}}_{2})}_{M_{2}}^{S_{2}}$ and $c_{M_{1}, M_{2}, M}^{S_{1}, S_{2}, S}$ are the Clebsch-Gordan coefficients. Thus, the elements $G_{z_{κ} w_{λ}}^{x_{ρ} y_{σ}} (τ)$ can be computed by solving the linear system of equations

(\begin{matrix} 0 & \frac{1}{2} & \frac{1}{2} & - \frac{1}{2} & - \frac{1}{2} & 0 \\ \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{3}} \\ 0 & - \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & - \frac{1}{2} & 0 \\ 0 & \frac{1}{2} & - \frac{1}{2} & \frac{1}{2} & - \frac{1}{2} & 0 \\ \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 & - \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{6}} & - \frac{1}{\sqrt{6}} & - \frac{1}{\sqrt{6}} & - \frac{1}{\sqrt{6}} & - \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} \end{matrix}) (\begin{matrix} ⟨ a_{x α}^{†} (τ) a_{y α}^{†} (τ) a_{w α} a_{z α} ⟩ \\ ⟨ a_{x α}^{†} (τ) a_{y β}^{†} (τ) a_{w β} a_{z α} ⟩ \\ ⟨ a_{x β}^{†} (τ) a_{y α}^{†} (τ) a_{w α} a_{z β} ⟩ \\ ⟨ a_{x β}^{†} (τ) a_{y α}^{†} (τ) a_{w β} a_{z α} ⟩ \\ ⟨ a_{x α}^{†} (τ) a_{y β}^{†} (τ) a_{w α} a_{z β} ⟩ \\ ⟨ a_{x β}^{†} (τ) a_{y β}^{†} (τ) a_{w β} a_{z β} ⟩ \end{matrix}) = (\begin{matrix} ⟨ {\hat{A}}_{x y}^{0} (τ) \otimes_{0} {\hat{A}}_{z w}^{0 ‡} ⟩ \\ ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{0} {\hat{A}}_{z w}^{1 ‡} ⟩ \\ ⟨ {\hat{A}}_{x y}^{0} (τ) \otimes_{1} {\hat{A}}_{z w}^{1 ‡} ⟩ \\ ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{1} {\hat{A}}_{z w}^{0 ‡} ⟩ \\ ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{1} {\hat{A}}_{z w}^{1 ‡} ⟩ \\ ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{2} {\hat{A}}_{z w}^{1 ‡} ⟩ \end{matrix}) .

(A1)

Equation (A1) can also be used to compute elements of $G^{[+ 2]} (τ)$ by replacing creation operators with annihilation operators and vice versa, e.g., $⟨ a_{x α}^{†} (τ) a_{y β}^{†} (τ) a_{w β} a_{z α} ⟩ \to ⟨ a_{x α} (τ) a_{y β} (τ) a_{w β}^{†} a_{z α}^{†} ⟩$ and $⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{0} {\hat{A}}_{z w}^{1 ‡} ⟩ \to ⟨ {\hat{A}}_{y x}^{1 ‡} (τ) \otimes_{0} {\hat{A}}_{w z}^{1} ⟩$ ⁠. To compute elements of $G^{[0^{'}]} (τ)$ ⁠, we define a product of tensor operators ${\hat{B}}_{x y}^{S, M} \equiv {[a_{x}^{†} a_{y}]}_{M}^{S}$ ⁠. In this case, the system of linear equations has the form

(\begin{matrix} \frac{1}{2} & \frac{1}{2} & 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{\sqrt{12}} & - \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & - \frac{1}{\sqrt{12}} & \frac{1}{\sqrt{12}} \\ - \frac{1}{2} & \frac{1}{2} & 0 & 0 & - \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & - \frac{1}{2} & - \frac{1}{2} \\ 0 & 0 & \frac{1}{\sqrt{2}} & - \frac{1}{\sqrt{2}} & 0 & 0 \\ - \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{6}} & - \frac{1}{\sqrt{6}} \end{matrix}) (\begin{matrix} ⟨ a_{x α}^{†} (τ) a_{y α} (τ) a_{w α}^{†} a_{z α} ⟩ \\ ⟨ a_{x α}^{†} (τ) a_{y α} (τ) a_{w β}^{†} a_{z β} ⟩ \\ ⟨ a_{x α}^{†} (τ) a_{y β} (τ) a_{w β}^{†} a_{z α} ⟩ \\ ⟨ a_{x β}^{†} (τ) a_{y α} (τ) a_{w α}^{†} a_{z β} ⟩ \\ ⟨ a_{x β}^{†} (τ) a_{y β} (τ) a_{w α}^{†} a_{z α} ⟩ \\ ⟨ a_{x β}^{†} (τ) a_{y β} (τ) a_{w β}^{†} a_{z β} ⟩ \end{matrix}) = (\begin{matrix} ⟨ {\hat{B}}_{x y}^{0} (τ) \otimes_{0} {\hat{B}}_{z w}^{0 ‡} ⟩ \\ ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{0} {\hat{B}}_{z w}^{1 ‡} ⟩ \\ ⟨ {\hat{B}}_{x y}^{0} (τ) \otimes_{1} {\hat{B}}_{z w}^{1 ‡} ⟩ \\ ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{1} {\hat{B}}_{z w}^{0 ‡} ⟩ \\ ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{1} {\hat{B}}_{z w}^{1 ‡} ⟩ \\ ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{2} {\hat{B}}_{z w}^{1 ‡} ⟩ \end{matrix}) .

(A2)

Equations (A1) and (A2) can be simplified for a closed-shell reference wavefunction $| Ψ_{0} (⟩$ to obtain compact expressions

⟨ a_{x α}^{†} (τ) a_{y α}^{†} (τ) a_{w α} a_{z α} ⟩ = \frac{1}{\sqrt{3}} ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{0} {\hat{A}}_{z w}^{1 ‡} ⟩,

(A3)

⟨ a_{x α}^{†} (τ) a_{y β}^{†} (τ) a_{w β} a_{z α} ⟩ = \frac{1}{\sqrt{12}} ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{0} {\hat{A}}_{z w}^{1 ‡} ⟩ + \frac{1}{2} ⟨ {\hat{A}}_{x y}^{0} (τ) \otimes_{0} {\hat{A}}_{z w}^{0 ‡} ⟩,

(A4)

⟨ a_{x β}^{†} (τ) a_{y α}^{†} (τ) a_{w β} a_{z α} ⟩ = \frac{1}{\sqrt{12}} ⟨ {\hat{A}}_{x y}^{1} (τ) \otimes_{0} {\hat{A}}_{z w}^{1 ‡} ⟩ - \frac{1}{2} ⟨ {\hat{A}}_{x y}^{0} (τ) \otimes_{0} {\hat{A}}_{z w}^{0 ‡} ⟩,

(A5)

⟨ a_{x α}^{†} (τ) a_{y β} (τ) a_{w β}^{†} a_{z α} ⟩ = \frac{1}{\sqrt{3}} ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{0} {\hat{B}}_{z w}^{1 ‡} ⟩,

(A6)

⟨ a_{x α}^{†} (τ) a_{y α} (τ) a_{w α}^{†} a_{z α} ⟩ = \frac{1}{2} ⟨ {\hat{B}}_{x y}^{0} (τ) \otimes_{0} {\hat{B}}_{z w}^{0 ‡} ⟩ + \frac{1}{\sqrt{12}} ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{0} {\hat{B}}_{z w}^{1 ‡} ⟩,

(A7)

⟨ a_{x α}^{†} (τ) a_{y α} (τ) a_{w β}^{†} a_{z β} ⟩ = \frac{1}{2} ⟨ {\hat{B}}_{x y}^{0} (τ) \otimes_{0} {\hat{B}}_{z w}^{0 ‡} ⟩ - \frac{1}{\sqrt{12}} ⟨ {\hat{B}}_{x y}^{1} (τ) \otimes_{0} {\hat{B}}_{z w}^{1 ‡} ⟩ .

(A8)

The remaining spin-orbital elements can be obtained using the following symmetry relations: $⟨ a_{x α}^{†} (τ) a_{y α}^{†} (τ) a_{w α} a_{z α} ⟩ = ⟨ a_{x β}^{†} (τ) a_{y β}^{†} (τ) a_{w β} a_{z β} ⟩$ ⁠, $⟨ a_{x α}^{†} (τ) a_{y β}^{†} (τ) a_{w β} a_{z α} ⟩ = ⟨ a_{x β}^{†} (τ) a_{y α}^{†} (τ) a_{w α} a_{z β} ⟩$ ⁠, $⟨ a_{x β}^{†} (τ) a_{y α}^{†} (τ) a_{w β} a_{z α} ⟩ = ⟨ a_{x α}^{†} (τ) a_{y β}^{†} (τ) a_{w α} a_{z β} ⟩$ ⁠, $⟨ a_{x α}^{†} (τ) a_{y α} (τ) a_{w α}^{†} a_{z α} ⟩ = ⟨ a_{x β}^{†} (τ) a_{y β} (τ) a_{w β}^{†} a_{z β} ⟩$ ⁠, $⟨ a_{x α}^{†} (τ) a_{y β} (τ) a_{w β}^{†} a_{z α} ⟩ = ⟨ a_{x β}^{†} (τ) a_{y α} (τ) a_{w α}^{†} a_{z β} ⟩$ ⁠, $⟨ a_{x α}^{†} (τ) a_{y α} (τ) a_{w β}^{†} a_{z β} ⟩ = ⟨ a_{x β}^{†} (τ) a_{y β} (τ) a_{w α}^{†} a_{z α} ⟩$ ⁠.

APPENDIX B: STRONGLY CONTRACTED NEVPT2 WITH MPS COMPRESSION (sc-MPS-NEVPT2)

In this section, we describe an efficient implementation of strongly contracted NEVPT2 with MPS reference wavefunctions (sc-MPS-NEVPT2). Although a sc-NEVPT2 implementation on top of MPS wavefunctions was reported in Ref. 35, here we additionally introduce MPS compression, which greatly reduces the computational cost and allows for larger active spaces to be treated. We will not describe the general theory of sc-NEVPT2 in detail but instead refer the reader to Refs. 22 and 23.

In sc-NEVPT2, for each of the 7 subspaces (⁠ $i \in {+ 1, - 1, + 2, - 2, + 1^{'}, - 1^{'}, 0^{'}}$ ⁠) described in Sec. II A, the first-order wavefunction is expanded in a basis consisting of a single perturber function for each core or external index. The perturber function is constructed by fixing (a set of) of core or external indices in $\hat{V}$ and acting on $| Ψ_{0} (⟩$ ⁠. For example, in the $[+ 1^{'}]$ subclass, the perturber functions are $| Ψ_{i}^{[+ 1^{'}]} (⟩ = {\hat{h}}_{i}^{†} a_{i} | Ψ_{0} (⟩$ ⁠. Because the perturber functions are all orthogonal to each other, they lead to very simple expressions for the second-order energy. For example, the perturbation contribution in the $[+ 1^{'}]$ subclass is

\begin{array}{rcl} E^{[+ 1^{'}]} & = & \sum_{i} \frac{| ⟨ Ψ_{0} | a_{i}^{†} {\hat{h}}_{i} | {\tilde{Ψ}}_{i}^{[+ 1^{'}]} ⟩ |^{2}}{E_{{\tilde{Ψ}}_{i}^{[+ 1^{'}]}} - E_{0}} \\ = & \sum_{i} \frac{| | | Ψ_{i}^{[+ 1^{'}]} ⟩) | |^{2}}{E_{{\tilde{Ψ}}_{i}^{[+ 1^{'}]}} - E_{0}}, \end{array}

(B1)

where $| {\tilde{Ψ}}_{i}^{[+ 1^{'}]} (⟩ = \frac{| Ψ_{i}^{[+ 1^{'}]} (⟩}{| | | Ψ_{i}^{[+ 1^{'}]} (⟩ | |}$ and $E_{{\tilde{Ψ}}_{i}^{[+ 1^{'}]}} = \frac{⟨ Ψ_{i}^{[+ 1^{'}]} | {\hat{H}}_{D} | Ψ_{i}^{[+ 1^{'}]} ⟩}{⟨ Ψ_{i}^{[+ 1^{'}]} | Ψ_{i}^{[+ 1^{'}]} ⟩}$ ⁠.

The main bottleneck in obtaining the sc-NEVPT2 energy is computing the term $⟨ Ψ_{i}^{[+ 1^{'}]} | {\hat{H}}_{D} | Ψ_{i}^{[+ 1^{'}]} ⟩$ and the analogous contribution for $| Ψ_{a}^{[- 1^{'}]} (⟩$ ⁠. This corresponds to

\begin{array}{rcl} ⟨ Ψ_{i}^{[+ 1^{'}]} | {\hat{H}}_{D} | Ψ_{i}^{[+ 1^{'}]} ⟩ & = & ⟨ Ψ_{0} | a_{i}^{†} {\hat{h}}_{i} {\hat{H}}_{D} {\hat{h}}_{i}^{†} a_{i} | Ψ_{0} ⟩ \\ = & ⟨ Ψ_{0} | {\hat{h}}_{i} {\hat{H}}_{D} {\hat{h}}_{i}^{†} | Ψ_{0} ⟩ - 𝜖_{i} ⟨ Ψ_{0} | {\hat{h}}_{i} {\hat{h}}_{i}^{†} | Ψ_{0} ⟩ \\ = & ⟨ Ψ_{0} | {\hat{h}}_{i} [{\hat{H}}_{D}, {\hat{h}}_{i}^{†}] | Ψ_{0} ⟩ \\ + (E_{0} - 𝜖_{i}) ⟨ Ψ_{0} | {\hat{h}}_{i} {\hat{h}}_{i}^{†} | Ψ_{0} ⟩ . \end{array}

(B2)

The above object involves 8 active indices, and the expectation value can be computed from the contraction of the 4-RDM and two-electron integrals.^22,23,92

However, similarly to how we avoid calculating the 3-GF in Eqs. (8) and (9) in t-NEVPT2, we can avoid calculating the 4-RDM in sc-NEVPT2. We can represent $| Ψ_{i}^{[+ 1^{'}]} (⟩$ and $| Ψ_{a}^{[- 1^{'}]} (⟩$ as MPS of bond dimension $M_{1}$ ⁠, by carrying out a variational compression as described in Sec. II C 2 to minimize $| | | Ψ_{i}^{[+ 1^{'}]} (⟩ - {\hat{h}}_{i}^{†} a_{i} | Ψ_{0} (⟩ | |$ and $| | | Ψ_{a}^{[- 1^{'}]} (⟩ - a_{a}^{†} {\hat{h}}_{a} | Ψ_{0} (⟩ | |$ ⁠. Assuming $M_{1} \sim M_{0}$ ⁠, the cost of carrying out compressions for all the perturbers in class $[+ 1^{'}]$ and $[- 1^{'}]$ is $O (N_{e x t} N_{a c t}^{2} M_{0}^{3})$ ⁠. The computational cost for the subsequent ${\hat{H}}_{D}$ expectation value is $O (N_{e x t} N_{a c t}^{3} M_{0}^{3})$ ⁠; however, as no sweeps are required, this has a very low prefactor, and in our calculations, the amount of time in this step is negligible. Thus, the computational cost of sc-NEVPT2 with MPS compression (sc-MPS-NEVPT2) scales in practice as $O (N_{a c t}^{4} M_{0}^{3} + N_{a c t}^{6} M_{0}^{2} + N_{e x t} N_{a c t}^{2} M_{0}^{3})$ ⁠. This is significantly lower than the $O (N_{a c t}^{5} M_{0}^{3} + N_{a c t}^{8} M_{0}^{2})$ cost of computing the 4-RDM, which dominates the standard sc-NEVPT2 algorithm presented in Ref. 35.

REFERENCES

1.

M.

Schütz

and

H.-J.

Werner

,

J. Chem. Phys.

114

,

661

(

2001

).

https://doi.org/10.1063/1.1330207

Google Scholar

Crossref

2.

H.-J.

Werner

,

F. R.

Manby

, and

P. J.

Knowles

,

J. Chem. Phys.

118

,

8149

(

2003

).

https://doi.org/10.1063/1.1564816

Google Scholar

Crossref

3.

C.

Riplinger

and

F.

Neese

,

J. Chem. Phys.

138

,

034106

(

2013

).

https://doi.org/10.1063/1.4773581

Google Scholar

Crossref

PubMed

4.

C.

Riplinger

,

B.

Sandhoefer

,

A.

Hansen

, and

F.

Neese

,

J. Chem. Phys.

139

,

134101

(

2013

).

https://doi.org/10.1063/1.4821834

Google Scholar

Crossref

PubMed

5.

S. R.

White

and

R. L.

Martin

,

J. Chem. Phys.

110

,

4127

(

1999

).

https://doi.org/10.1063/1.478295

Google Scholar

Crossref

6.

Ö.

Legeza

,

R.

Noack

,

J.

Sólyom

, and

L.

Tincani

,

Computational Many-Particle Physics

(

Springer

,

2008

), pp.

653

–

664

.

Google Scholar

Crossref

7.

G. H.

Booth

,

A. J. W.

Thom

, and

A.

Alavi

,

J. Chem. Phys.

131

,

054106

(

2009

).

https://doi.org/10.1063/1.3193710

Google Scholar

Crossref

PubMed

8.

Y.

Kurashige

and

T.

Yanai

,

J. Chem. Phys.

130

,

234114

(

2009

).

https://doi.org/10.1063/1.3152576

Google Scholar

Crossref

PubMed

9.

K. H.

Marti

and

M.

Reiher

,

Phys. Chem. Chem. Phys.

13

,

6750

(

2011

).

https://doi.org/10.1039/c0cp01883j

Google Scholar

Crossref

PubMed

10.

G. K.-L.

Chan

and

S.

Sharma

,

Annu. Rev. Phys. Chem.

62

,

465

(

2011

).

https://doi.org/10.1146/annurev-physchem-032210-103338

Google Scholar

Crossref

PubMed

11.

S.

Wouters

and

D.

Van Neck

,

Eur. Phys. J. D

68

,

272

(

2014

).

https://doi.org/10.1140/epjd/e2014-50500-1

Google Scholar

Crossref

12.

Y.

Kurashige

,

G. K.-L.

Chan

, and

T.

Yanai

,

Nat. Chem.

5

,

660

(

2013

).

https://doi.org/10.1038/nchem.1677

Google Scholar

Crossref

PubMed

13.

J.

Chalupský

,

T. A.

Rokob

,

Y.

Kurashige

,

T.

Yanai

,

E. I.

Solomon

,

L.

Rulíšek

, and

M.

Srnec

,

J. Am. Chem. Soc.

136

,

15977

(

2014

).

https://doi.org/10.1021/ja506934k

Google Scholar

Crossref

PubMed

14.

S.

Wouters

,

T.

Bogaerts

,

P.

Van Der Voort

,

V.

Van Speybroeck

, and

D.

Van Neck

,

J. Chem. Phys.

140

,

241103

(

2014

).

https://doi.org/10.1063/1.4885815

Google Scholar

Crossref

PubMed

15.

S.

Sharma

,

K.

Sivalingam

,

F.

Neese

, and

G. K.-L.

Chan

,

Nat. Chem.

6

,

927

(

2014

).

https://doi.org/10.1038/nchem.2041

Google Scholar

Crossref

PubMed

16.

M. K.

MacLeod

and

T.

Shiozaki

,

J. Chem. Phys.

142

,

051103

(

2015

).

https://doi.org/10.1063/1.4907717

Google Scholar

Crossref

PubMed

17.

S.

Evangelisti

,

J. P.

Daudey

, and

J.-P. P.

Malrieu

,

Phys. Rev. A

35

,

4930

(

1987

).

https://doi.org/10.1103/physreva.35.4930

Google Scholar

Crossref

18.

K.

Kowalski

and

P.

Piecuch

,

Int. J. Quantum Chem.

80

,

757

(

2000

).

https://doi.org/10.1002/1097-461x(2000)80:4/5<757::aid-qua25>3.0.co;2-a

Google Scholar

Crossref

19.

K.

Kowalski

and

P.

Piecuch

,

Phys. Rev. A

61

,

052506

(

2000

).

https://doi.org/10.1103/physreva.61.052506

Google Scholar

Crossref

20.

F. A.

Evangelista

,

J. Chem. Phys.

141

,

054109

(

2014

).

https://doi.org/10.1063/1.4890660

Google Scholar

Crossref

PubMed

21.

H.-J.

Werner

and

P. J.

Knowles

,

J. Chem. Phys.

89

,

5803

(

1988

).

https://doi.org/10.1063/1.455556

Google Scholar

Crossref

22.

C.

Angeli

,

R.

Cimiraglia

,

S.

Evangelisti

,

T.

Leininger

, and

J.-P. P.

Malrieu

,

J. Chem. Phys.

114

,

10252

(

2001

).

https://doi.org/10.1063/1.1361246

Google Scholar

Crossref

23.

C.

Angeli

,

R.

Cimiraglia

, and

J.-P. P.

Malrieu

,

Chem. Phys. Lett.

350

,

297

(

2001

).

https://doi.org/10.1016/s0009-2614(01)01303-3

Google Scholar

Crossref

24.

K.

Hirao

,

Chem. Phys. Lett.

190

,

374

(

1992

).

https://doi.org/10.1016/0009-2614(92)85354-d

Google Scholar

Crossref

25.

K.

Andersson

,

P. Å.

Malmqvist

, and

B. O.

Roos

,

J. Chem. Phys.

96

,

1218

(

1992

).

https://doi.org/10.1063/1.462209

Google Scholar

Crossref

26.

H.-J.

Werner

,

Mol. Phys.

89

,

645

(

1996

).

https://doi.org/10.1080/00268979609482499

Google Scholar

Crossref

27.

T.

Yanai

and

G. K.-L.

Chan

,

J. Chem. Phys.

124

,

194106

(

2006

).

https://doi.org/10.1063/1.2196410

Google Scholar

Crossref

PubMed

28.

T.

Yanai

and

G. K.-L.

Chan

,

J. Chem. Phys.

127

,

104107

(

2007

).

https://doi.org/10.1063/1.2761870

Google Scholar

Crossref

PubMed

29.

Y.

Kurashige

and

T.

Yanai

,

J. Chem. Phys.

135

,

094104

(

2011

).

https://doi.org/10.1063/1.3629454

Google Scholar

Crossref

PubMed

30.

M.

Saitow

,

Y.

Kurashige

, and

T.

Yanai

,

J. Chem. Phys.

139

,

044118

(

2013

).

https://doi.org/10.1063/1.4816627

Google Scholar

Crossref

PubMed

31.

M.

Saitow

,

Y.

Kurashige

, and

T.

Yanai

,

J. Chem. Theory Comput.

11

,

5120

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00270

Google Scholar

Crossref

PubMed

32.

D.

Datta

,

L.

Kong

, and

M.

Nooijen

,

J. Chem. Phys.

134

,

214116

(

2011

).

https://doi.org/10.1063/1.3592494

Google Scholar

Crossref

PubMed

33.

F. A.

Evangelista

and

J.

Gauss

,

J. Chem. Phys.

134

,

114102

(

2011

).

https://doi.org/10.1063/1.3559149

Google Scholar

Crossref

PubMed

34.

A.

Köhn

,

M.

Hanauer

,

L. A.

Mück

,

T.-C.

Jagau

, and

J.

Gauss

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

3

,

176

(

2013

).

https://doi.org/10.1002/wcms.1120

Google Scholar

Crossref

35.

S.

Guo

,

M. A.

Watson

,

W.

Hu

,

Q.

Sun

, and

G. K.-L.

Chan

,

J. Chem. Theory Comput.

12

,

1583

(

2016

).

https://doi.org/10.1021/acs.jctc.5b01225

Google Scholar

Crossref

PubMed

36.

L.

Freitag

,

S.

Knecht

,

C.

Angeli

, and

M.

Reiher

,

J. Chem. Theory Comput.

13

,

451

(

2017

).

https://doi.org/10.1021/acs.jctc.6b00778

Google Scholar

Crossref

PubMed

37.

A. Y.

Sokolov

and

G. K.-L.

Chan

,

J. Chem. Phys.

144

,

064102

(

2016

).

https://doi.org/10.1063/1.4941606

Google Scholar

Crossref

PubMed

38.

K. G.

Dyall

,

J. Chem. Phys.

102

,

4909

(

1995

).

https://doi.org/10.1063/1.469539

Google Scholar

Crossref

39.

D.

Zgid

,

D.

Ghosh

,

E.

Neuscamman

, and

G. K.-L.

Chan

,

J. Chem. Phys.

130

,

194107

(

2009

).

https://doi.org/10.1063/1.3132922

Google Scholar

Crossref

PubMed

40.

U.

Schollwöck

,

Rev. Mod. Phys.

77

,

259

(

2005

).

https://doi.org/10.1103/revmodphys.77.259

Google Scholar

Crossref

41.

G. K.-L.

Chan

,

A.

Keselman

,

N.

Nakatani

,

Z.

Li

, and

S. R.

White

,

J. Chem. Phys.

145

,

014102

(

2016

).

https://doi.org/10.1063/1.4955108

Google Scholar

Crossref

PubMed

42.

A. E.

Feiguin

and

S. R.

White

,

Phys. Rev. B

72

,

020404

(

2005

).

https://doi.org/10.1103/physrevb.72.220401

Google Scholar

Crossref

43.

W. H.

Press

,

S. A.

Teukolsky

,

W. T.

Vetterling

, and

B. P.

Flannery

,

Numerical Recipes: The Art of Scientific Computing

(

Cambridge University Press

,

New York

,

1986

).

Google Scholar

44.

M.

Roemelt

,

S.

Guo

, and

G. K.-L.

Chan

,

J. Chem. Phys.

144

,

204113

(

2016

).

https://doi.org/10.1063/1.4950757

Google Scholar

Crossref

PubMed

45.

S.

Sharma

and

G. K.-L.

Chan

,

J. Chem. Phys.

136

,

124121

(

2012

).

https://doi.org/10.1063/1.3695642

Google Scholar

Crossref

PubMed

46.

W.

Tatsuaki

,

Phys. Rev. E

61

,

3199

(

2000

).

https://doi.org/10.1103/physreve.61.3199

Google Scholar

Crossref

47.

Q.

Sun

,

T. C.

Berkelbach

,

N. S.

Blunt

,

G. H.

Booth

,

S.

Guo

,

Z.

Li

,

J.

Liu

,

J.

McClain

,

E. R.

Sayfutyarova

,

S.

Sharma

,

S.

Wouters

, and

G.K.-L.

Chan

, e-print arXiv:1701.08223 (

2017

).

48.

T. H.

Dunning

,

J. Chem. Phys.

90

,

1007

(

1989

).

https://doi.org/10.1063/1.456153

Google Scholar

Crossref

49.

S.

Sharma

and

G. K.-L.

Chan

,

J. Chem. Phys.

141

,

111101

(

2014

).

https://doi.org/10.1063/1.4895977

Google Scholar

Crossref

PubMed

50.

K.

Andersson

,

B. O.

Roos

,

P. Å.

Malmqvist

, and

P. O.

Widmark

,

Chem. Phys. Lett.

230

,

391

(

1994

).

https://doi.org/10.1016/0009-2614(94)01183-4

Google Scholar

Crossref

51.

B. O.

Roos

and

K.

Andersson

,

Chem. Phys. Lett.

245

,

215

(

1995

).

https://doi.org/10.1016/0009-2614(95)01010-7

Google Scholar

Crossref

52.

K.

Andersson

,

Chem. Phys. Lett.

237

,

212

(

1995

).

https://doi.org/10.1016/0009-2614(95)00328-2

Google Scholar

Crossref

53.

H.

Stoll

and

H.-J.

Werner

,

Mol. Phys.

88

,

793

(

1996

).

https://doi.org/10.1080/00268979609482455

Google Scholar

Crossref

54.

H.

Dachsel

,

R. J.

Harrison

, and

D. A.

Dixon

,

J. Phys. Chem. A

103

,

152

(

1999

).

https://doi.org/10.1021/jp982648s

Google Scholar

Crossref

55.

B. O.

Roos

,

Collect. Czech. Chem. Commun.

68

,

265

(

2003

).

https://doi.org/10.1135/cccc20030265

Google Scholar

Crossref

56.

P.

Celani

,

H.

Stoll

,

H.-J.

Werner

, and

P. J.

Knowles

,

Mol. Phys.

102

,

2369

(

2004

).

https://doi.org/10.1080/00268970412331317788

Google Scholar

Crossref

57.

C.

Angeli

,

B.

Bories

,

A.

Cavallini

, and

R.

Cimiraglia

,

J. Chem. Phys.

124

,

054108

(

2006

).

https://doi.org/10.1063/1.2148946

Google Scholar

Crossref

PubMed

58.

F.

Ruipérez

,

F.

Aquilante

,

J. M.

Ugalde

, and

I.

Infante

,

J. Chem. Theory Comput.

7

,

1640

(

2011

).

https://doi.org/10.1021/ct200048z

Google Scholar

Crossref

PubMed

59.

T.

Müller

,

J. Phys. Chem. A

113

,

12729

(

2009

).

https://doi.org/10.1021/jp905254u

Google Scholar

Crossref

PubMed

60.

W.

Purwanto

,

S.

Zhang

, and

H.

Krakauer

,

J. Chem. Phys.

142

,

064302

(

2015

).

https://doi.org/10.1063/1.4906829

Google Scholar

Crossref

PubMed

61.

S. M.

Casey

and

D. G.

Leopold

,

J. Phys. Chem.

97

,

816

(

1993

).

https://doi.org/10.1021/j100106a005

Google Scholar

Crossref

62.

N. B.

Balabanov

and

K. A.

Peterson

,

J. Chem. Phys.

123

,

064107

(

2005

).

https://doi.org/10.1063/1.1998907

Google Scholar

Crossref

63.

W.

Liu

,

Mol. Phys.

108

,

1679

(

2010

).

https://doi.org/10.1080/00268971003781571

Google Scholar

Crossref

64.

T.

Saue

,

ChemPhysChem

12

,

3077

(

2011

).

https://doi.org/10.1002/cphc.201100682

Google Scholar

Crossref

PubMed

65.

D.

Peng

and

M.

Reiher

,

Theor. Chem. Acc.

131

,

1081

(

2012

).

https://doi.org/10.1007/s00214-011-1081-y

Google Scholar

Crossref

66.

R.

Pariser

,

J. Chem. Phys.

24

,

250

(

1956

).

https://doi.org/10.1063/1.1742461

Google Scholar

Crossref

67.

P.

Tavan

and

K.

Schulten

,

J. Chem. Phys.

85

,

6602

(

1986

).

https://doi.org/10.1063/1.451442

Google Scholar

Crossref

68.

P.

Tavan

and

K.

Schulten

,

Phys. Rev. B

36

,

4337

(

1987

).

https://doi.org/10.1103/physrevb.36.4337

Google Scholar

Crossref

69.

K.

Nakayama

,

H.

Nakano

, and

K.

Hirao

,

Int. J. Quantum Chem.

66

,

157

(

1998

).

https://doi.org/10.1002/(sici)1097-461x(1998)66:2<157::aid-qua7>3.0.co;2-u

Google Scholar

Crossref

70.

Y.

Kurashige

,

H.

Nakano

,

Y.

Nakao

, and

K.

Hirao

,

Chem. Phys. Lett.

400

,

425

(

2004

).

https://doi.org/10.1016/j.cplett.2004.10.141

Google Scholar

Crossref

71.

D.

Ghosh

,

J.

Hachmann

,

T.

Yanai

, and

G. K.-L.

Chan

,

J. Chem. Phys.

128

,

144117

(

2008

).

https://doi.org/10.1063/1.2883976

Google Scholar

Crossref

PubMed

72.

C.

Angeli

and

M.

Pastore

,

J. Chem. Phys.

134

,

184302

(

2011

).

https://doi.org/10.1063/1.3585607

Google Scholar

Crossref

PubMed

73.

M. A.

Watson

and

G. K.-L.

Chan

,

J. Chem. Theory Comput.

8

,

4013

(

2012

).

https://doi.org/10.1021/ct300591z

Google Scholar

Crossref

PubMed

74.

E.

Ronca

,

C.

Angeli

,

L.

Belpassi

,

F.

De Angelis

,

F.

Tarantelli

, and

M.

Pastore

,

J. Chem. Theory Comput.

10

,

4014

(

2014

).

https://doi.org/10.1021/ct5004675

Google Scholar

Crossref

PubMed

75.

B. S.

Hudson

and

B. E.

Kohler

,

J. Chem. Phys.

59

,

4984

(

1973

).

https://doi.org/10.1063/1.1680717

Google Scholar

Crossref

76.

T.

Sashima

,

H.

Nagae

,

M.

Kuki

, and

Y.

Koyama

,

Chem. Phys. Lett.

299

,

187

(

1999

).

https://doi.org/10.1016/s0009-2614(98)01278-0

Google Scholar

Crossref

77.

G.

Cerullo

,

Science

298

,

2395

(

2002

).

https://doi.org/10.1126/science.1074685

Google Scholar

Crossref

PubMed

78.

K.

Furuichi

,

T.

Sashima

, and

Y.

Koyama

,

Chem. Phys. Lett.

356

,

547

(

2002

).

https://doi.org/10.1016/s0009-2614(02)00412-8

Google Scholar

Crossref

79.

J. P.

Doering

and

R.

McDiarmid

,

J. Chem. Phys.

75

,

2477

(

1981

).

https://doi.org/10.1063/1.442278

Google Scholar

Crossref

80.

J. P.

Doering

and

R.

McDiarmid

,

J. Chem. Phys.

73

,

3617

(

1980

).

https://doi.org/10.1063/1.440587

Google Scholar

Crossref

81.

L. A.

Heimbrook

,

J. E.

Kenny

,

B. E.

Kohler

, and

G. W.

Scott

,

J. Chem. Phys.

75

,

4338

(

1981

).

https://doi.org/10.1063/1.442640

Google Scholar

Crossref

82.

L. A.

Heimbrook

,

B. E.

Kohler

, and

I. J.

Levy

,

J. Chem. Phys.

81

,

1592

(

1984

).

https://doi.org/10.1063/1.447888

Google Scholar

Crossref

83.

D. G.

Leopold

,

V.

Vaida

, and

M. F.

Granville

,

J. Chem. Phys.

81

,

4210

(

1984

).

https://doi.org/10.1063/1.447452

Google Scholar

Crossref

84.

H.

Petek

,

A. J.

Bell

,

Y. S.

Choi

,

K.

Yoshihara

,

B. A.

Tounge

, and

R. L.

Christensen

,

J. Chem. Phys.

98

,

3777

(

1993

).

https://doi.org/10.1063/1.464056

Google Scholar

Crossref

85.

B. E.

Kohler

,

C.

Spangler

, and

C.

Westerfield

,

J. Chem. Phys.

89

,

5422

(

1988

).

https://doi.org/10.1063/1.455594

Google Scholar

Crossref

86.

M.

Schmidt

and

P.

Tavan

,

J. Chem. Phys.

136

,

124309

(

2012

).

https://doi.org/10.1063/1.3696880

Google Scholar

Crossref

PubMed

87.

A. D.

Becke

,

J. Chem. Phys.

98

,

5648

(

1993

).

https://doi.org/10.1063/1.464913

Google Scholar

Crossref

88.

C.

Lee

,

W.

Yang

, and

R. G.

Parr

,

Phys. Rev. B

37

,

785

(

1988

).

https://doi.org/10.1103/physrevb.37.785

Google Scholar

Crossref

89.

J. M.

Turney

,

A. C.

Simmonett

,

R. M.

Parrish

,

E. G.

Hohenstein

,

F. A.

Evangelista

,

J. T.

Fermann

,

B. J.

Mintz

,

L. A.

Burns

,

J. J.

Wilke

,

M. L.

Abrams

,

N. J.

Russ

,

M. L.

Leininger

,

C. L.

Janssen

,

E. T.

Seidl

,

W. D.

Allen

,

H. F.

Schaefer

,

R. A.

King

,

E. F.

Valeev

,

C. D.

Sherrill

, and

T. D.

Crawford

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

2

,

556

(

2012

).

https://doi.org/10.1002/wcms.93

Google Scholar

Crossref

90.

G.

Knizia

,

J. Chem. Theory Comput.

9

,

4834

(

2013

).

https://doi.org/10.1021/ct400687b

Google Scholar

Crossref

PubMed

91.

J.

Pipek

and

P. G.

Mezey

,

J. Chem. Phys.

90

,

4916

(

1989

).

https://doi.org/10.1063/1.456588

Google Scholar

Crossref

92.

C.

Angeli

,

R.

Cimiraglia

, and

J.-P. P.

Malrieu

,

J. Chem. Phys.

117

,

9138

(

2002

).

https://doi.org/10.1063/1.1515317

Google Scholar

Crossref

2017

Author(s)

Time-dependent N-electron valence perturbation theory with matrix product state reference wavefunctions for large active spaces and basis sets: Applications to the chromium dimer and all-trans polyenes

I. INTRODUCTION

II. THEORY

A. Time-dependent formulation of N-electron valence perturbation theory (t-NEVPT2)

B. Matrix product state (MPS) wavefunctions and the density matrix renormalization group (DMRG) algorithm

C. t-NEVPT2 with MPS reference wavefunctions (t-MPS-NEVPT2)

1. Overview of the algorithm

2. Green’s function evaluation

3. Spin-adaptation

III. COMPUTATIONAL DETAILS

IV. RESULTS

A. Analysis of numerical accuracy: $N_{2}$ dissociation

B. Chromium dimer

C. Singlet excited states in all-trans polyenes

1. Background

2. Details of computations

3. Discussion

V. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: SPIN-ORBITAL TWO-BODY GREEN’S FUNCTIONS IN SPIN-ADAPTED DMRG

APPENDIX B: STRONGLY CONTRACTED NEVPT2 WITH MPS COMPRESSION (sc-MPS-NEVPT2)

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Time-dependent N-electron valence perturbation theory with matrix product state reference wavefunctions for large active spaces and basis sets: Applications to the chromium dimer and all-trans polyenes

I. INTRODUCTION

II. THEORY

A. Time-dependent formulation of N-electron valence perturbation theory (t-NEVPT2)

B. Matrix product state (MPS) wavefunctions and the density matrix renormalization group (DMRG) algorithm

C. t-NEVPT2 with MPS reference wavefunctions (t-MPS-NEVPT2)

1. Overview of the algorithm

2. Green’s function evaluation

3. Spin-adaptation

III. COMPUTATIONAL DETAILS

IV. RESULTS

A. Analysis of numerical accuracy: N2 dissociation

B. Chromium dimer

C. Singlet excited states in all-trans polyenes

1. Background

2. Details of computations

3. Discussion

V. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: SPIN-ORBITAL TWO-BODY GREEN’S FUNCTIONS IN SPIN-ADAPTED DMRG

APPENDIX B: STRONGLY CONTRACTED NEVPT2 WITH MPS COMPRESSION (sc-MPS-NEVPT2)

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

A. Analysis of numerical accuracy: $N_{2}$ dissociation