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

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

arXiv:1703.10830·physics.chem-ph·April 9, 2018

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

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

PDF

TL;DR

This paper introduces a novel time-dependent NEVPT2 method combined with matrix product state wavefunctions, enabling efficient dynamic correlation calculations for large active spaces and basis sets, demonstrated on chromium dimer and polyenes.

Contribution

It develops a new t-MPS-NEVPT2 approach and a low-scaling MPS-based strongly-contracted NEVPT2 method for large-scale electronic structure calculations.

Findings

01

Accurate dissociation energy for chromium dimer

02

Low-lying excited states in polyenes

03

Efficient handling of large active spaces

Abstract

In earlier work [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 to compute uncontracted dynamic correlation energies for large active spaces and basis sets, using the time-dependent density matrix renormalization group (td-DMRG) 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…

Tables4

Table 1. Table 1: Equations for the elements of tensors ϵ [ 𝒊 ] ( τ ) superscript bold-italic-ϵ delimited-[] 𝒊 𝜏 \boldsymbol{\epsilon^{[i]}}(\tau) and 𝐆 ( τ ) [ 𝒊 ] \mathbf{G}\boldsymbol{{}^{[i]}}(\tau) used in the evaluation of the t-NEVPT2 active-space energy contributions E a c t [ i ] superscript subscript 𝐸 𝑎 𝑐 𝑡 delimited-[] 𝑖 E_{act}^{[i]} in Eq. 16 .

$[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} ⟩$

Table 2. Table 2: t-MPS-NEVPT2 correlation energy contributions E a c t [ i ] superscript subscript 𝐸 𝑎 𝑐 𝑡 delimited-[] 𝑖 E_{act}^{[i]} for Cr 2 at 1.650 Å computed for different values of the bond dimension M 𝑀 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 Δ subscript 𝐸 𝑐 𝑜 𝑛 𝑡 𝑟 \Delta E_{contr} in the correlation energy for Cr 2 , the Cr atom, and the Cr 2 binding energy ( D e subscript 𝐷 𝑒 D_{e} ), evaluated as the difference between the sc-DMRG-NEVPT2 and t-MPS-NEVPT2 energies.

$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} ({Cr}_{2}) = E_{sc - NEVPT2}^{(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} (Cr) = E_{sc - NEVPT2}^{(2)} (Cr) - E_{t - NEVPT2}^{(2)} (Cr)$ = $-$ 0.7758 $+$ 0.7781 = 0.0023 $E_{h}$
$Δ E_{c o n t r} (D_{e})$ = $Δ E_{c o n t r} ({Cr}_{2})$ $-$ $2 \times Δ E_{c o n t r} (Cr)$ = 0.0034 $E_{h}$ = 0.09 eV

Table 3. Table 3: 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 − 1 superscript subscript 𝐴 𝑔 {1}A_{g}^{-} ground states give the total energies in E h subscript 𝐸 ℎ 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 ( n 𝑛 n e, m 𝑚 m o) 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
\ce C4H6	$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
\ce C8H10	$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
\ce C16H18	$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
\ce C24H26	$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

Table 4. Table 4: Energies for the lowest-lying singlet states in the all-trans conjugated polyenes \ce C4H6 and \ce C8H10. Entries for the 1 A g − 1 superscript subscript 𝐴 𝑔 {1}A_{g}^{-} ground states give the total energies in E h subscript 𝐸 ℎ 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 ( n 𝑛 n e, m 𝑚 m o) denotes the active spaces used in the CASSCF computations.

Polyene	State	CASSCF	sc-NEVPT2	t-NEVPT2	Reference
\ce C4H6	$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
\ce C8H10	$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

Equations57

\hat{H}_{D} =

\hat{H}_{D} =

f_{p}^{q}

f_{p}^{q}

\hat{H}_{a c t}

\hat{H}_{a c t}

E^{(2)}

E^{(2)}

= - \int_{0}^{\infty} d τ ⟨ Ψ_{0} ⟩ \hat{V}^{' †} e^{- (\hat{H}_{D} - E_{D}) τ} \hat{V}^{'} Ψ_{0},

E^{(2)}

E^{(2)}

+ E^{[+ 1^{'}]} + E^{[- 1^{'}]} + E^{[0^{'}]} .

\hat{h}_{p}^{†}

\hat{h}_{p}^{†}

\tilde{h}_{p}^{q}

E^{[+ 1^{'}]}

E^{[+ 1^{'}]}

E^{[- 1^{'}]}

∣ Ψ ⟩ = n_{1} \dots n_{k} \sum A^{n_{1}} A^{n_{2}} \dots A^{n_{p}} \dots A^{n_{k}} ∣ n_{1} n_{2} \dots n_{k} ⟩,

∣ Ψ ⟩ = n_{1} \dots n_{k} \sum A^{n_{1}} A^{n_{2}} \dots A^{n_{p}} \dots A^{n_{k}} ∣ n_{1} n_{2} \dots n_{k} ⟩,

∣ l_{α_{p}} ⟩

∣ l_{α_{p}} ⟩

∣ r_{α_{p + 1}} ⟩

∣ Ψ ⟩ = l_{α_{p}} n_{p} r_{α_{p + 1}} \sum [A^{n_{p}}]_{α_{p} α_{p + 1}} ∣ l_{α_{p}} n_{p} r_{α_{p + 1}} ⟩ .

∣ Ψ ⟩ = l_{α_{p}} n_{p} r_{α_{p + 1}} \sum [A^{n_{p}}]_{α_{p} α_{p + 1}} ∣ l_{α_{p}} n_{p} r_{α_{p + 1}} ⟩ .

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

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

E_{e x t}^{(2)}

E_{e x t}^{(2)}

\displaystyle+2\sum_{ixya}\frac{\tilde{h}_{{a}}^{{i}}{v}_{{ix}}^{{ay}}\gamma_{{y}}^{{x}}}{\varepsilon_{i}-\varepsilon_{a}}\

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

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

∣ Φ_{μ}^{[i]} (τ^{'}) ⟩ \approx ∣ Φ_{μ}^{[i]} (τ) ⟩ + n = 1 \sum 6 c_{n} ∣ k_{n} (δ τ) ⟩,

∣ Φ_{μ}^{[i]} (τ^{'}) ⟩ \approx ∣ Φ_{μ}^{[i]} (τ) ⟩ + n = 1 \sum 6 c_{n} ∣ k_{n} (δ τ) ⟩,

∣ k_{n} (δ τ) ⟩ = - δ τ (\hat{H}_{D} - E_{D}) [∣ Φ_{μ}^{[i]} (τ) ⟩ + m = 1 \sum n - 1 b_{nm} ∣ k_{m} (δ τ) ⟩] .

∣ k_{n} (δ τ) ⟩ = - δ τ (\hat{H}_{D} - E_{D}) [∣ Φ_{μ}^{[i]} (τ) ⟩ + m = 1 \sum n - 1 b_{nm} ∣ k_{m} (δ τ) ⟩] .

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

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

0 \frac{1}{3} 00 \frac{1}{2} \frac{1}{6} \frac{1}{2} \frac{1}{12} - \frac{1}{2} \frac{1}{2} 0 - \frac{1}{6} \frac{1}{2} \frac{1}{12} \frac{1}{2} - \frac{1}{2} 0 - \frac{1}{6} - \frac{1}{2} \frac{1}{12} \frac{1}{2} \frac{1}{2} 0 - \frac{1}{6} - \frac{1}{2} \frac{1}{12} - \frac{1}{2} - \frac{1}{2} 0 - \frac{1}{6} 0 \frac{1}{3} 00 - \frac{1}{2} \frac{1}{6} ⟨ 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 β} ⟩ = ⟨ \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 ‡} ⟩ .

0 \frac{1}{3} 00 \frac{1}{2} \frac{1}{6} \frac{1}{2} \frac{1}{12} - \frac{1}{2} \frac{1}{2} 0 - \frac{1}{6} \frac{1}{2} \frac{1}{12} \frac{1}{2} - \frac{1}{2} 0 - \frac{1}{6} - \frac{1}{2} \frac{1}{12} \frac{1}{2} \frac{1}{2} 0 - \frac{1}{6} - \frac{1}{2} \frac{1}{12} - \frac{1}{2} - \frac{1}{2} 0 - \frac{1}{6} 0 \frac{1}{3} 00 - \frac{1}{2} \frac{1}{6} ⟨ 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 β} ⟩ = ⟨ \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 ‡} ⟩ .

\frac{1}{2} \frac{1}{12} - \frac{1}{2} \frac{1}{2} 0 - \frac{1}{6} \frac{1}{2} - \frac{1}{12} \frac{1}{2} \frac{1}{2} 0 \frac{1}{6} 0 \frac{1}{3} 00 \frac{1}{2} \frac{1}{6} 0 \frac{1}{3} 00 - \frac{1}{2} \frac{1}{6} \frac{1}{2} - \frac{1}{12} - \frac{1}{2} - \frac{1}{2} 0 \frac{1}{6} \frac{1}{2} \frac{1}{12} \frac{1}{2} - \frac{1}{2} 0 - \frac{1}{6} ⟨ 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 β} ⟩ = ⟨ \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 ‡} ⟩ .

\frac{1}{2} \frac{1}{12} - \frac{1}{2} \frac{1}{2} 0 - \frac{1}{6} \frac{1}{2} - \frac{1}{12} \frac{1}{2} \frac{1}{2} 0 \frac{1}{6} 0 \frac{1}{3} 00 \frac{1}{2} \frac{1}{6} 0 \frac{1}{3} 00 - \frac{1}{2} \frac{1}{6} \frac{1}{2} - \frac{1}{12} - \frac{1}{2} - \frac{1}{2} 0 \frac{1}{6} \frac{1}{2} \frac{1}{12} \frac{1}{2} - \frac{1}{2} 0 - \frac{1}{6} ⟨ 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 β} ⟩ = ⟨ \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 ‡} ⟩ .

⟨ 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 β} ⟩

E^{[+ 1^{'}]}

E^{[+ 1^{'}]}

= i \sum \frac{∣∣ ∣ Ψ _{i}^{[+ 1^{'}]} ⟩ ∣ ∣ ^{2}}{E _{\tilde{Ψ}_{i}^{[+ 1^{'}]}} - E _{0}},

⟨ Ψ_{i}^{[+ 1^{'}]} ⟩ \hat{H}_{D} Ψ_{i}^{[+ 1^{'}]}

⟨ Ψ_{i}^{[+ 1^{'}]} ⟩ \hat{H}_{D} Ψ_{i}^{[+ 1^{'}]}

= ⟨ Ψ_{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} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

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

Alexander Yu. Sokolov

[email protected]

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA

Sheng Guo

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA

Enrico Ronca

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA

Garnet Kin-Lic Chan

[email protected]

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA

Abstract

In earlier work [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 to compute uncontracted dynamic correlation energies for large active spaces and basis sets, using the time-dependent density matrix renormalization group (td-DMRG) 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 (\ceC4H6 to \ceC24H26), 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.Schutz:2001p661 ; Werner:2003p8149 ; Riplinger:2013p034106 ; Riplinger:2013p134101 Similarly, efficient methods for static correlationWhite:1999p4127 ; Legeza2008 ; Booth:2009p054106 ; Kurashige:2009p234114 ; Marti:2011p6750 ; Chan:2011p465 ; Wouters:2014p272 now exist that allow us to simulate systems with a large number of strongly correlated open-shell orbitals.Kurashige:2013p660 ; Chalupsky:2014p15977 ; Wouters:2014p241103 ; Sharma:2014p927 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.Kurashige:2013p660 ; Sharma:2014p927 Technically, the challenges are the high computational cost of the existing canonical algorithms, algebraic complexity of the underlying equations,MacLeod:2015p051103 as well as numerical instabilities in the simulations.Evangelisti:1987p4930 ; Kowalski:2000p757 ; Kowalski:2000p052506 ; Evangelista:2014p054109

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,Werner:1988p5803 ; Angeli:2001p10252 ; Angeli:2001p297 which represents the complicated active-space wavefunction in terms of simpler quantities, such as the reduced density matrices (RDM). While this approximation has been enormously useful in quantum chemistry, Hirao:1992p374 ; Werner:1988p5803 ; Andersson:1992p1218 ; Werner:1996p645 ; Angeli:2001p10252 ; Angeli:2001p297 ; Yanai:2006p194106 ; Yanai:2007p104107 ; Kurashige:2011p094104 ; Saitow:2013p044118 ; Saitow:2015p5120 ; Datta:2011p214116 ; Evangelista:2011p114102 ; Kohn:2012p176 ; Guo:2016p1583 ; Freitag:2017p451 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 theorySokolov:2016p064102 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). 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 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) algorithmWhite:1999p4127 ; Kurashige:2009p234114 ; Marti:2011p6750 ; Chan:2011p465 ; Wouters:2014p272 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),Angeli:2001p10252 ; Angeli:2001p297 ; Guo:2016p1583 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 (Section II.1) and matrix product state wavefunctions (Section II.2). We then describe the details of our t-MPS-NEVPT2 implementation (Section II.3) and outline the computational details (Section III). In Section IV, we use t-MPS-NEVPT2 to compute energies along the \ceN2 dissociation curve (Section IV.1), the dissociation energy of the chromium dimer (Section IV.2), and the vertical excitation energies in all-trans polyenes (Section IV.3). In Section V, we present the conclusions of our work. The details of the efficient sc-MPS-NEVPT2 implementation are described in the Appendix B.

II Theory

II.1 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 $\ket{\Psi_{0}}$ computed in a complete active space (CAS) of molecular orbitals. We require $\ket{\Psi_{0}}$ to be an eigenfunction of a zeroth-order Hamiltonian $\hat{H}^{(0)}$ . Following our previous work,Sokolov:2016p064102 we assume that $\hat{H}^{(0)}$ is the Dyall Hamiltonian,Dyall:1995p4909 $\hat{H}^{(0)}=\hat{H}_{D}$ , defined as:

[TABLE]

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 $\varepsilon_{i}$ and $\varepsilon_{a}$ are defined as the core and external eigenvalues of the generalized Fock operator,

[TABLE]

where $h_{{p}}^{{q}}$ and ${v}_{{pq}}^{{rs}}$ are the usual one- and antisymmetrized two-electron integrals, and $\gamma_{{q}}^{{p}}=\braket{\Psi_{0}}{a^{\dagger}_{p}a^{\ }_{q}}{\Psi_{0}}$ is the one-particle density matrix of $\ket{\Psi_{0}}$ . The $\hat{H}_{act}$ operator describes the two-electron interaction in the active space:

[TABLE]

Having specified $\ket{\Psi_{0}}$ and $\hat{H}^{(0)}$ , we now consider an expansion of the energy with respect to the perturbation $\lambda\hat{V}=\lambda(\hat{H}-\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 expansionSokolov:2016p064102 with respect to the perturbation operator $\hat{V}(\tau)=e^{(\hat{H}_{D}-E_{D})\tau}\,\hat{V}\,e^{-(\hat{H}_{D}-E_{D})\tau}$ . The second-order energy contribution can be written as:

[TABLE]

where $\tau$ is imaginary time and $\hat{V}^{\prime}$ denotes the part of $\hat{V}$ that contributes to the first-order wavefunction $\ket{\Psi^{(1)}}$ ( $\hat{V}^{\prime}=\hat{Q}\hat{V}$ , $\hat{Q}=1-\ket{\Psi_{0}}\bra{\Psi_{0}}$ ). Section II.1 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).Angeli:2001p10252 ; Angeli:2001p297 However, unlike the standard uncontracted theory, the time-dependent energy expression (II.1) 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):

[TABLE]

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

[TABLE]

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

[TABLE]

The bottleneck of the time-dependent NEVPT2 (t-NEVPT2) algorithm, when $\ket{\Psi_{0}}$ is expanded in determinants, is the evaluation of the 2-GF that has $\mathcal{O}(N_{\tau}\times N_{det}\times N_{act}^{6})$ computational cost, where $N_{act}$ is the number of active orbitals, $N_{det}$ is the dimension of the CAS Hilbert space, and $N_{\tau}$ is the number of time steps ( $N_{\tau}$ $\sim$ 10-20). As a result, t-NEVPT2 has a significantly lower scaling with $N_{act}$ than the internally-contracted NEVPT2 approaches,Angeli:2001p10252 ; Angeli:2001p297 which require computation of up to the four-particle reduced density matrix (4-RDM) with $\mathcal{O}(N_{det}\times N_{act}^{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_{det}\sim\mathcal{O}(e^{N_{act}})$ , and solving for $\ket{\Psi_{0}}$ becomes a bottleneck for large active spaces. In the following sections we will describe how this limitation can be ameliorated by expressing $\ket{\Psi_{0}}$ in the matrix product state (MPS) representation.

II.2 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:

[TABLE]

where $\ket{n_{1}n_{2}\ldots n_{k}}$ is a Slater determinant, $\mathbf{A}^{n_{p}}$ is a tensor that corresponds to only one orbital $p$ with occupancy $n_{p}\in\{\ket{},\ket{\uparrow},\ket{\downarrow},\ket{\uparrow\downarrow}\}$ (also referred as the site $p$ ), and the total number of tensors $\mathbf{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,\ldots,k-1\}$ , $\mathbf{A}^{n_{p}}$ is a matrix with fixed dimensions $M\times M$ . The dimensions of $\mathbf{A}^{n_{1}}$ and $\mathbf{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}\ldots n_{k}$ , the product of tensors in Eq. 10 yields a scalar corresponding to the coefficient of the determinant $\ket{n_{1}n_{2}\ldots n_{k}}$ in the expansion of the wavefunction $\ket{\Psi}$ . 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\ldots 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

[TABLE]

and the total wavefunction becomes

[TABLE]

In this context, we can think of $\mathbf{A}^{n_{p}}$ as representing a “site wavefunction”, and operators projected into the renormalized space $\{\ket{l_{\alpha_{p}}n_{p}r_{\alpha_{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.

II.3 t-NEVPT2 with MPS reference wavefunctions (t-MPS-NEVPT2)

II.3.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

[TABLE]

where we combine terms involving the core and external orbitals

[TABLE]

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

[TABLE]

where $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ is the Green’s function tensor and $\boldsymbol{\epsilon^{[i]}}(\tau)$ is the integral prefactor tensor that can be computed using the one- and two-electron integrals and orbital energies ( $\varepsilon_{i}$ and $\varepsilon_{a}$ ). Explicit equations for the elements of tensors $\boldsymbol{\epsilon^{[i]}}(\tau)$ and $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ are given in Table 1.

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 ( $\tau$ ), and (iii) integration in imaginary time. In step (i), the elements of $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ are computed for a set of $\tau$ values (time steps, $\tau_{k}$ ). In step (ii), the prefactors $\boldsymbol{\epsilon^{[i]}}(\tau)$ are evaluated using the equations in Table 1, and the energy contributions at each time step are computed as $E_{act}^{[i]}(\tau_{k})=\boldsymbol{\epsilon^{[i]}}(\tau_{k})\mathbf{G}\boldsymbol{{}^{[i]}}(\tau_{k})$ . For example, $E_{act}^{[-2]}(\tau_{k})$ is computed as the dot product of two tensors: $\epsilon^{[-2]}(\tau_{k})=\epsilon_{zw}^{xy}(\tau_{k})=\frac{1}{8}\sum_{ab}{v}_{{zw}}^{{ab}}{v}_{{ab}}^{{xy}}e^{-(\varepsilon_{a}+\varepsilon_{b})\tau_{k}}$ and $G^{[-2]}(\tau_{k})=G_{zw}^{xy}(\tau_{k})=\braket{\Psi_{0}}{a^{\dagger}_{x}(\tau_{k})a^{\dagger}_{y}(\tau_{k})a^{\ }_{w}a^{\ }_{z}}{\Psi_{0}}$ . Finally, in step (iii), each energy contribution $E_{act}^{[i]}$ is evaluated by fitting the computed values $E_{act}^{[i]}(\tau_{k})$ to an exponential function $\sum_{i}\mathcal{E}(\tau)=a_{i}e^{-b_{i}\tau}$ and integrating the obtained result. We note that the above algorithm outline can be used to implement t-NEVPT2 for reference wavefunctions in any representation (e.g., determinant-basedSokolov:2016p064102 or MPS-based). Only the details of step (i) depend on the explicit form of $\ket{\Psi_{0}}$ . In Section II.3.2, we will describe how imaginary-time Green’s functions $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ are evaluated in the MPS representation.

II.3.2 Green’s function evaluation

To organize our MPS implementation of the imaginary-time Green’s functions, we express the elements of $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ (Table 1) in the general form $G_{\mu\nu}^{[i]}(\tau)=\braket{\Phi_{\mu}^{[i]}(\tau)}{\Phi_{\nu}^{[i]}}$ , where $\ket{\Phi_{\nu}^{[i]}}$ is obtained by applying all creation and annihilation operators at $\tau=0$ on $\ket{\Psi_{0}}$ (e.g., $a^{\ }_{x}\ket{\Psi_{0}}$ , $a^{\dagger}_{x}a^{\ }_{y}\ket{\Psi_{0}}$ ), while the application of operators at $\tau=\tau^{\prime}$ yields $\ket{\Phi_{\mu}^{[i]}(\tau^{\prime})}$ (e.g., $a^{\ }_{x}(\tau^{\prime})\ket{\Psi_{0}}$ , $a^{\dagger}_{x}(\tau^{\prime})a^{\ }_{y}(\tau^{\prime})\ket{\Psi_{0}}$ ). In our notation, indices $\mu$ and $\nu$ run over the total number of states $\ket{\Phi_{\mu}^{[i]}}$ (i.e., $N_{act}$ for $a^{\ }_{x}\ket{\Psi_{0}}$ or $N^{2}_{act}$ for $a^{\dagger}_{x}a^{\ }_{y}\ket{\Psi_{0}}$ ).

We begin the computation of $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ by optimizing the zeroth-order MPS wavefunction $\ket{\Psi_{0}}$ using the DMRG algorithm with bond dimension $M_{0}$ . To evaluate the elements of $\mathbf{G}\boldsymbol{{}^{[i]}}(\tau)$ we employ the algorithm below.

Compute initial wavefunctions $\ket{\Phi_{\mu}^{[i]}}$ for every $\mu$ . Each wavefunction $\ket{\Phi_{\mu}^{[i]}}$ is computed as an individual MPS and stored on disk. Evaluation of $\ket{\Phi_{\mu}^{[i]}}$ ( $[i]$ $\in$ $\{[+1]$ , $[-1]$ , $[+2]$ , $[-2]$ , $[0^{\prime}]\}$ ) requires applying the operator $a^{\dagger}_{x}$ or $a^{\ }_{x}$ on $\ket{\Psi_{0}}$ one or more times. The resulting MPS is of exactly the same bond dimension as the original MPS. For $[i]$ $\in$ $\{[+1^{\prime}]$ , $[-1^{\prime}]\}$ , the wavefunctions $\ket{\Phi_{\mu}^{[+1^{\prime}]}}$ and $\ket{\Phi_{\mu}^{[-1^{\prime}]}}$ are computed by compressing the MPS form of $\hat{h}_{i}^{\dagger}\ket{\Psi_{0}}$ and $\hat{h}_{a}\ket{\Psi_{0}}$ , where $\hat{h}_{i}^{\dagger}$ and $\hat{h}_{a}$ are defined as in Eq. 6. Since applying $\hat{h}_{i}^{\dagger}$ or $\hat{h}_{a}$ involves a summation over the active-space labels $x$ , the resulting MPS is of larger bond dimension than $\ket{\Psi_{0}}$ . We use variational compression to obtain $\ket{\Phi_{\mu}^{[+1^{\prime}]}}$ and $\ket{\Phi_{\mu}^{[-1^{\prime}]}}$ MPS by minimizing $||\ket{\Phi_{i}^{[+1^{\prime}]}}-\hat{h}_{i}^{\dagger}\ket{\Psi_{0}}||$ and $||\ket{\Phi_{a}^{[-1^{\prime}]}}-\hat{h}_{a}\ket{\Psi_{0}}||$ , where $\ket{\Phi_{i}^{[+1^{\prime}]}}$ and $\ket{\Phi_{a}^{[-1^{\prime}]}}$ have bond dimension $M_{1}>M_{0}$ . To maximize the efficiency of computing the overlaps appearing in the compression (e.g. $\bra{\Phi_{i}^{[+1^{\prime}]}}\hat{h}_{i}^{\dagger}\ket{\Psi_{0}}$ ) we use the standard DMRG formalism of “normal” and “complementary” operators, which requires building $\mathcal{O}(k)$ such operators from the left and right blocks. In practice, we find that $M_{1}\approx 2M_{0}$ usually provides a sufficiently accurate compression. Overall, this step of the t-MPS-NEVPT2 algorithm has $\mathcal{O}(N_{act}^{2}M_{0}^{3})+\mathcal{O}(N_{ext}N_{act}^{2}M_{1}^{2}M_{0})$ scaling for computing all $\ket{\Phi_{\mu}^{[i]}}$ wavefunctions with $[i]$ $\in$ $\{[+2]$ , $[-2]$ , $[0^{\prime}]\}$ and $[i]$ $\in$ $\{[+1^{\prime}]$ , $[-1^{\prime}]\}$ , respectively ( $N_{ext}$ is the number of external orbitals). 2. 2.

Propagate wavefunctions $\ket{\Phi_{\mu}^{[i]}}$ in imaginary time $\tau$ according to the time-dependent Schrödinger equation $\ket{\Phi_{\mu}^{[i]}(\tau)}=e^{-(\hat{H}_{D}-E_{D})\tau}\ket{\Phi_{\mu}^{[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,Press:2007 which allows to automatically control the time step $\delta\tau$ based on the error estimate of the fifth-order approximation of the propagator $e^{-(\hat{H}_{D}-E_{D})\tau}$ . In the ERK (4,5) method, the wavefunction at time $\tau^{\prime}=\tau+\delta\tau$ is approximated to $\mathcal{O}((\delta\tau)^{5})$ as

[TABLE]

where states $\ket{k_{n}(\delta\tau)}$ are obtained by successively applying the zeroth-order Hamiltonian six times:

[TABLE]

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 $\ket{\Phi_{\mu}^{[i]}(\tau^{\prime})}$ , which we denote as $\ket{\Phi_{\mu}^{[i]}(\tau^{\prime})[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, $\Delta=\sqrt{\left|\left|\ket{\Phi_{\mu}^{[i]}(\tau^{\prime})}-\ket{\Phi_{\mu}^{[i]}(\tau^{\prime})[4]}\right|\right|}$ .

As mentioned in Section II.2, 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 $\ket{\Phi_{\mu}^{[i]}}$ is now the site wavefunction, the states $\ket{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 $\Delta_{p}$ . As in Ref. 42, the site wavefunctions $\ket{\Phi_{\mu}^{[i]}(\tau)}$ , $\ket{\Phi_{\mu}^{[i]}(\tau+\delta\tau/3)}$ , $\ket{\Phi_{\mu}^{[i]}(\tau+2\delta\tau/3)}$ , $\ket{\Phi_{\mu}^{[i]}(\tau+\delta\tau)}$ 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 $\delta\tau\sim 10^{-3}$ and determine the new time-step after each propagation as $\delta\tau^{\prime}=\min(2\times\delta\tau,\delta\tau_{emb})$ , where $\delta\tau_{emb}=\delta\tau\left|\frac{\Delta E_{conv}}{\max_{p}(\{\Delta_{p}\})}\right|^{1/5}$ and $\Delta E_{conv}$ is the specified accuracy threshold. Note that if $\max_{p}(\{\Delta_{p}\})>\Delta E_{conv}$ , 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 $\ket{\Phi_{\mu}^{[i]}}$ has $\mathcal{O}(N_{act}^{5}M_{0}^{3})+\mathcal{O}(N_{act}^{6}M_{0}^{2})+\mathcal{O}(N_{ext}N_{act}^{3}M_{1}^{3})+\mathcal{O}(N_{ext}N_{act}^{4}M_{1}^{2})$ scaling. 3. 3.

Compute Green’s function elements at every time step $\tau_{k}$ as the overlap between two MPS $G_{\mu\nu}^{[i]}(\tau_{k})=\braket{\Phi_{\mu}^{[i]}(\tau_{k})}{\Phi_{\nu}^{[i]}}$ . The computational cost of computing all $G_{\mu\nu}^{[i]}(\tau_{k})$ elements has $\mathcal{O}(N_{act}^{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 $\mathcal{O}(N_{act}^{2})$ MPS wavefunctions $\ket{\Phi_{\mu}^{[i]}}$ and evaluation of $\mathcal{O}(N_{act}^{4})$ matrix elements $G_{\mu\nu}^{[i]}(\tau_{k})$ , leading to a total $\mathcal{O}(N_{\tau}N_{act}^{5}M_{0}^{3})+\mathcal{O}(N_{\tau}N_{act}^{6}M_{0}^{2})$ computational scaling, where $N_{\tau}$ $\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 $\mathcal{O}(N_{act}^{4}M_{0}^{3})+\mathcal{O}(N_{act}^{6}M_{0}^{2})$ scaling, is of lower cost.Guo:2016p1583 The higher scaling of our 2-GF implementation is because each MPS $\ket{\Phi_{\mu}^{[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 have 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 $\ket{\Phi_{\mu}^{[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 $\mathcal{O}(N_{act}^{2})$ components:

[TABLE]

where each component $E_{\mu}^{[i]}$ can be evaluated independently with $\mathcal{O}(N_{\tau}N_{act}^{3}M_{0}^{3})+\mathcal{O}(N_{\tau}N_{act}^{4}M_{0}^{2})$ cost. In addition, as we will demonstrate in Section 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.

II.3.3 Spin-adaptation

In Section II.3.2 we discussed the evaluation of the $G_{\mu\nu}^{[i]}(\tau)$ matrix elements in terms of creation and annihilation operators $a^{\dagger}_{p}(\tau)$ and $a^{\ }_{p}(\tau)$ . In a non-spin-adapted DMRG algorithm, $a^{\dagger}$ and $a$ are 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_{\lambda}w_{\kappa}}^{x_{\rho}y_{\sigma}}(\tau)=\braket{\Psi_{0}}{a^{\dagger}_{x_{\rho}}(\tau)a^{\dagger}_{y_{\sigma}}(\tau)a^{\ }_{w_{\kappa}}a^{\ }_{z_{\lambda}}}{\Psi_{0}}$ , where we use the Roman characters $x,y,w,z$ to denote spatial orbitals and Greek characters $\rho,\sigma,\kappa,\lambda$ to denote the spin of these orbitals. Here, $a^{\ }_{w_{\kappa}}$ and $a^{\ }_{z_{\lambda}}$ are simple operators, and application of a pair of operators $a^{\ }_{w_{\kappa}}a^{\ }_{z_{\lambda}}\ket{\Psi_{0}}$ results in a single MPS wavefunction $\ket{\Phi_{w_{\kappa}z_{\lambda}}^{[-2]}}$ .

In a spin-adapted DMRG implementation,Sharma:2012p124121 $a^{\dagger}$ 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. $\braket{\Psi_{0}}{\{[a^{\dagger}_{x}(\tau)a^{\dagger}_{y}(\tau)]^{S_{1}}_{M_{1}}[a^{\ }_{w}a^{\ }_{z}]^{S_{2}}_{M_{2}}\}^{S_{3}}_{M_{3}}}{\Psi_{0}}$ , where we use brackets to denote the coupling of spins $S_{1}$ , $S_{2}$ , and $S_{3}$ for different tensor products. As an example, we consider the evaluation of $G_{\mu\nu}^{[-2]}(\tau)$ for a reference state $\ket{\Psi_{0}}$ with $S=0$ . (Generalization to reference states with $S\neq 0$ is straightforward and is simplified using singlet embedding).Tatsuaki:2000p3199 ; Sharma:2012p124121 Applying a pair of tensor operators on the reference state $[a^{\ }_{w}a^{\ }_{z}]^{S_{1}}_{M_{1}}\ket{\Psi_{0}}$ results in two sets of eigenstates $\ket{{\Phi_{wz}^{[-2]}}_{M_{1}}^{S_{1}}}$ with $S_{1}$ $=$ $0,1$ . The wavefunctions $\ket{{\Phi_{wz}^{[-2]}}_{M_{1}}^{S_{1}}}$ are propagated in imaginary time for each value of $S_{1}$ , 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_{zw}^{xy})^{S_{1}}(\tau)=\sum_{M_{1}}(-1)^{S_{1}+M_{1}}c_{M_{1},-M_{1},0}^{S_{1},S_{1},0}\braket{\Psi_{0}}{[a^{\dagger}_{x}(\tau)a^{\dagger}_{y}(\tau)]^{S_{1}}_{M_{1}}[a^{\ }_{w}a^{\ }_{z}]^{S_{1}}_{M_{1}}}{\Psi_{0}}$ , where we used the fact that $S_{1}=S_{2}$ and $S_{3}=0$ for a closed-shell $\ket{\Psi_{0}}$ . Expanding the spin tensor operators $[a^{\ }_{w}a^{\ }_{z}]^{S_{1}}$ as a combination of spin-orbital operators $a^{\ }_{w_{\rho}}$ and $a^{\ }_{z_{\sigma}}$ leads to a system of linear equations that can be solved to obtain the spin-orbital 2-GF elements $G_{z_{\lambda}w_{\kappa}}^{x_{\rho}y_{\sigma}}(\tau)$ from the matrix elements $(G_{zw}^{xy})^{S_{1}}(\tau)$ (see the 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.Sharma:2012p124121 The t-MPS-NEVPT2 and sc-MPS-NEVPT2 correlation energies were computed using Pyscf,Sun:ArXiv08223 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 ( $n$ e, $m$ o), 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*-9* $E_{h}$ ) and orbital optimization iterations (10*-6* $E_{h}$ ). In t-MPS-NEVPT2, the accuracy of time propagation was controlled by setting $\Delta E_{conv}=10^{-4}$ $E_{h}$ for all systems but the N2 molecule, where the tighter threshold $\Delta E_{conv}=10^{-6}$ $E_{h}$ was used (see Section II.3.2 for details). Specific details pertaining to our computations of the all-trans polyenes are described in Section IV.3.

As described in Section II.3, 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_{0}$ and $M_{1}$ to the same value, and whenever we refer to bond dimension $M$ we imply that $M_{0}=M_{1}=M$ .

IV Results

IV.1 Analysis of numerical accuracy: N2 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.Sokolov:2016p064102 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 \ceN2 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^{\prime}$ = 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.Dunning:1989p1007

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 \ceN2 potential energy curve (PEC). As the bond dimension $M$ increases, the t-MPS-NEVPT2 errors decrease, becoming smaller than 0.1 $mE_{h}$ at $M$ = 200 for all energy points. As can be seen from Figure 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 (Figure 1). For all geometries with $r$ $\leq$ 2.1 Å, our t-MPS-NEVPT2 algorithm achieves an accuracy of better than $10^{-6}$ $E_{h}$ , while errors in the range of $\sim$ $10^{-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 \ceN2 dissociation using the cc-pVDZ basis set and the (10e,8o) active space.Sharma:2014p111101

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_{act}^{[i]}$ ( $[i]$ $\in$ $\{[+1]$ , $[-1]$ , $[+2]$ , $[-2]$ , $[0^{\prime}]\}$ ), and (ii) the sum of energy terms that depend on the contracted 3-GF ( $E_{act}^{[+1^{\prime}]}$ + $E_{act}^{[-1^{\prime}]}$ ). These errors are plotted in Figures 2a and 2b 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 contribution (ii) from the contracted 3-GF. At $M$ = 50, contribution (i) shows errors of $\sim$ $10^{-4}$ $E_{h}$ for a large range of geometries, while the errors in contribution (ii) are an order of magnitude larger. Comparing Figures 2a and 2b, 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. Figure 1 and Figure 2b). As discussed in Section II.3.2, the observed slower convergence of contribution (ii) is the result of the MPS compression used to evaluate $E_{act}^{[+1^{\prime}]}$ and $E_{act}^{[-1^{\prime}]}$ , 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$ $2M$ to obtain contribution (ii).

IV.2 Chromium dimer

The chromium dimer (\ceCr2) is a notoriously challenging system for multi-reference methods. Computational studiesAndersson:1994p391 ; Roos:1995p215 ; Andersson:1995p212 ; Stoll:1996p793 ; Dachsel:1999p152 ; Roos:2003p265 ; Celani:2004p2369 ; Angeli:2006p054108 ; Ruiperez:2011p1640 ; Muller:2009p12729 ; Kurashige:2011p094104 ; Purwanto:2015p064302 ; Sokolov:2016p064102 ; Guo:2016p1583 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 PEC of this molecule. In 2011, Kurashige and Yanai showed that the \ceCr2 dissociation curve can be accurately described by combining the CASPT2 variant of internally-contracted perturbation theory with DMRG in a large (12e,28o) active space,Kurashige:2011p094104 which included the $3d$ , $4d$ , $4s$ , and $4p$ orbitals of the Cr atoms. Very recently (2016), Guo et al. computed \ceCr2 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).Guo:2016p1583 At the complete basis set limit (CBS), the \ceCr2 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}^{exp}$ = 1.679 Å and $D_{e}^{exp}$ = 1.47 $\pm$ 0.1 eV.Casey:1993p816 However, the \ceCr2 binding energy in sc-NEVPT2 can be significantly affected by the strong contraction approximation, as we demonstrated in our t-NEVPT2 study of the \ceCr2 PEC using the (12e,12o) active space.Sokolov:2016p064102 In this work, we recompute the \ceCr2 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.Guo:2016p1583 In short, we used the (12e,22o) active space and optimized the molecular orbitals using the DMRG-SCF algorithm with $M^{\prime}$ = 1000 (no frozen core). We then performed a DMRG computation with $M^{\prime}$ = 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^{\prime}$ was compressed down to an MPS with a smaller bond dimension $M$ . We computed each t-MPS-NEVPT2 energy contribution $E_{act}^{[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 setBalabanov:2005p064107 combined with the spin-free one-electron variant of the X2C Hamiltonian.Liu:2010p1679 ; Saue:2011p3077 ; Peng:2012p1081 Since strong contraction reduces the binding energySokolov:2016p064102 in \ceCr2, 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,Sokolov:2016p064102 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 2 presents 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_{act}^{[i]}$ that depend on the 1-GF ( $[i]$ $\in$ $\{[+1]$ , $[-1]\}$ ) and 2-GF ( $[i]$ $\in$ $\{[+2]$ , $[-2]$ , $[0^{\prime}]\}$ ). Out of five of these contributions, the largest change of +0.3 $mE_{h}$ is observed for $E_{act}^{[0^{\prime}]}$ , while each of the other four energy terms changes by less than $-$ 0.1 $mE_{h}$ . These changes largely cancel each other, leading to a total energy difference of $\sim$ 0.1 $mE_{h}$ . Thus, we estimate the total error in the sum of the 1- and 2-GF energy contributions to be less than 0.1 $mE_{h}$ . However, the remaining energy components $E_{act}^{[+1^{\prime}]}$ and $E_{act}^{[-1^{\prime}]}$ that depend on the contracted 3-GF exhibit slower convergence with $M$ , changing by $-$ 1.9 $mE_{h}$ for $M$ = 500 $\rightarrow$ 800 and $-$ 1.2 $mE_{h}$ for $M$ = 800 $\rightarrow$ 1200 (Table 2). We increased the bond dimension up to $M$ = 2000, where the remaining error in $E_{act}^{[+1^{\prime}]}$ and $E_{act}^{[-1^{\prime}]}$ is estimated to be less than 0.4 $mE_{h}$ ( $<$ 0.01 eV). The resulting active-space energy contributions $E_{act}^{[i]}$ can be used to compute the \ceCr2 binding energy for the incomplete cc-pwCV5Z-DK basis set as $D_{e}$ = $E^{(2)}_{ext}$ $+$ $\sum_{[i]}E_{act}^{[i]}$ $-$ $E^{(2)}_{\mathrm{t-NEVPT2}}(\mathrm{Cr})$ , where $E^{(2)}_{ext}$ is the \ceCr2 correlation energy from the external space (Section II.3.1) and $E^{(2)}_{\mathrm{t-NEVPT2}}(\mathrm{Cr})$ 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 $\Delta E_{contr}(D_{e})$ as shown in Table 2, where $E^{(2)}_{\mathrm{sc-NEVPT2}}$ is the sc-DMRG-NEVPT2 correlation energy for \ceCr2 reported by Guo et al. ( $M$ = 1200),Guo:2016p1583 while $E^{(2)}_{\mathrm{sc-NEVPT2}}(\mathrm{Cr})$ and $E^{(2)}_{\mathrm{t-NEVPT2}}(\mathrm{Cr})$ 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 \ceCr2 binding energy at the t-MPS-NEVPT2 level of theory can be estimated as $(D_{e})_{\mathrm{t-NEVPT2}}^{\mathrm{CBS}}$ $\approx$ $(D_{e})_{\mathrm{sc-NEVPT2}}^{\mathrm{CBS}}$ + $\Delta E_{contr}(D_{e})$ , where we use $(D_{e})_{\mathrm{sc-NEVPT2}}^{\mathrm{CBS}}$ = 1.43 eV from Ref. 35. Although the resulting binding energy $(D_{e})_{\mathrm{t-NEVPT2}}^{\mathrm{CBS}}$ = 1.52 eV is 0.09 eV larger than $(D_{e})_{\mathrm{sc-NEVPT2}}^{\mathrm{CBS}}$ , it is still in a good agreement with the experimental binding energy of 1.47 $\pm$ 0.1 eV reported by Casey and Leopold.Casey:1993p816

IV.3 Singlet excited states in all-trans polyenes

IV.3.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 $\pi$ - $\pi^{*}$ excited states of the conjugated backbone. One of the simplest examples of polyene systems are the unsubstituted all-trans polyenes that consist of alternating single and double carbon bonds arranged in a linear chain (CnHn+2). The excited states of these molecules are labeled by the symmetry of the $C_{2h}$ 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.Pariser:1956p250 ; Tavan:1986p6602

The electronic spectra of all-trans polyenes have been the subject of many experimental and computational studies.Tavan:1986p6602 ; Tavan:1987p4337 ; Nakayama:1998p157 ; Kurashige:2004p425 ; Ghosh:2008p144117 ; Zgid:2009p194107 ; Angeli:2011p184302 ; Watson:2012p4013 ; Ronca:2014p4014 ; Hudson:1973p4984 ; Sashima:1999p187 ; Cerullo:2002p2395 ; Furuichi:2002p547 ; Doering:1981p2477 ; Doering:1980p3617 ; Heimbrook:1981p4338 ; Heimbrook:1984p1592 ; Leopold:1984p4210 ; Petek:1993p3777 ; Kohler:1988p5422 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 (CnHn+2, $n$ $\leq$ 8), it has been established that the two lowest-energy singlet transitions are the dipole-allowed excitation to the ionic $1B_{u}^{+}$ state and the dipole-forbidden excitation to the ${2}A_{g}^{-}$ state. Despite the fact that the allowed transition ${1}A_{g}^{-}$ $\rightarrow$ $1B_{u}^{+}$ has primarily HOMO $\rightarrow$ LUMO excitation character, fluorescence studies suggest that, starting from octatetraene (\ceC8H10),Hudson:1973p4984 the lowest-energy singlet excitation corresponds to a forbidden transition, such as the ${1}A_{g}^{-}$ $\rightarrow$ ${2}A_{g}^{-}$ transition. Recent high-level theoretical studies of butadiene (\ceC4H6) predict the bright $1B_{u}^{+}$ state vertical excitation to be $\sim$ 0.2 eV lower in energy than the excitation to the dark ${2}A_{g}^{-}$ state,Watson:2012p4013 whereas the two vertical transitions are predicted to be almost degenerate in \ceC8H10.Angeli:2011p184302 In longer polyenes, the low-energy electronic spectrum is more complicated, involving additional dark states $1B_{u}^{-}$ and ${3}A_{g}^{-}$ with energies close to $1B_{u}^{+}$ and ${2}A_{g}^{-}$ . These dark states have been observed experimentally in *all-trans-*carotenoids with $n$ = 18–26 $\pi$ -electrons.Sashima:1999p187 ; Cerullo:2002p2395 ; Furuichi:2002p547

Early computational studies of long polyenes have been primarily based on model Hamiltonians, such as the Pariser-Parr-Pople (PPP) and Hubbard models.Tavan:1986p6602 ; Tavan:1987p4337 In 1986, Tavan and SchultenTavan:1986p6602 performed multi-reference configuration interaction (MRCI) computations using the PPP model, which suggested that at $n$ $=$ 10 the $1B_{u}^{-}$ state becomes the second excited state, leading to the ${1}A_{g}^{-}$ $<$ ${2}A_{g}^{-}$ $<$ $1B_{u}^{-}$ $<$ $1B_{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,Schmidt:2012p124309 which predicted the inversion of the $1B_{u}^{-}$ and $1B_{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.Nakayama:1998p157 In a later study, Kurashige et al. used a combination of MRMP and CASCI with the (10e,10o) active space (incomplete $\pi$ -valence space) for polyene series up to $n$ = 28.Kurashige:2004p425 Their computations also suggested that the $1B_{u}^{-}$ and $1B_{u}^{+}$ vertical transitions cross at $n$ = 14 and predicted a crossing of ${3}A_{g}^{-}$ and $1B_{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 $\pi$ -valence active space ( $n$ e, $n$ o).Ghosh:2008p144117 Although in this study the ionic $1B_{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}^{-}$ , $1B_{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 $1B_{u}^{+}$ state, excitation energies of the covalent states computed using DMRG-SCF overestimate those obtained from the experiment.Ghosh:2008p144117 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 $\pi$ -valence space description of static correlation.

IV.3.2 Details of computations

The ground-state molecular geometries of the all-trans polyenes CnHn+2 ( $n$ = 4, 8, 16, and 24) were optimized using the B3LYP methodBecke:1993p5648 ; Lee:1988p785 implemented in the Psi4 program.Turney:2012p556 Geometry optimizations were constrained to have $C_{2h}$ symmetry. For all polyenes, we used the cc-pVDZ basis set.Dunning:1989p1007 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 $\pi$ -valence active space ( $n$ e, $n$ o). 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 (IAO),Knizia:2013p4834 from which the $2p_{z}$ orbitals of carbon atoms ( $z$ being the $C_{2}$ 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^{\prime}$ = 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,Pipek:1989p4916 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 were first computed with $M^{\prime}$ = 500 and then compressed down to $M$ = 250. For CnHn+2 with $n$ = 4, 8, and 16, the remaining error in the reference and correlation energies was estimated to be less than 0.1 $mE_{h}$ . To achieve a similar accuracy for \ceC24H26, we used $M$ = 500 in sc-MPS-NEVPT2 and to compute the $E_{act}^{[+1^{\prime}]}$ and $E_{act}^{[-1^{\prime}]}$ 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 $1B_{u}^{+}$ state was found to be the 3rd, 5th, 7th, and 9th singlet eigenstate of the zeroth-order Hamiltonian, respectively. Thus, to compute the $1B_{u}^{+}$ excitation energies, we increased the number of states evaluated in the state-averaged DMRG computation with $M^{\prime}$ = 500 to 9 and 11 for $n$ = 16 and 24, respectively.

IV.3.3 Discussion

Table 3 presents energies, symmetries, and oscillator strengths for the low-lying singlet states of polyenes CnHn+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, Figure 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}^{-}$ $<$ $1B_{u}^{-}$ $<$ ${3}A_{g}^{-}$ in close agreement with the results obtained by Ghosh et al.Ghosh:2008p144117 In particular, for all polyenes (including \ceC4H6), DMRG-SCF predicts the ${2}A_{g}^{-}$ state to be the lowest-energy singlet excited state, while the energy of the ionic $1B_{u}^{+}$ state dominated by the HOMO $\rightarrow$ LUMO transition is computed to be 1.6 eV to 2.3 eV higher.

Including the dynamic correlation lowers the relative energy of the ionic $1B_{u}^{+}$ state drastically. This is seen in the difference of the $1B_{u}^{+}$ excitation energies computed using sc-MPS-NEVPT2 and DMRG-SCF, which ranges from 2 eV in \ceC4H6 to 2.7 eV in \ceC24H26, 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 $1B_{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 $1B_{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}^{-}$ , $1B_{u}^{-}$ , and ${3}A_{g}^{-}$ is much smaller (Figure 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 \ceC24H26 at the t-MPS-NEVPT2 level of theory. Interestingly, in shorter polyenes (\ceC4H6 and \ceC8H10), 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 Figure 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 (\ceC4H6 and \ceC8H10) strong contraction mainly affects the energy of the ionic $1B_{u}^{+}$ state, in the longer polyenes the errors of strong contraction become significant even for the covalent ${2}A_{g}^{-}$ , $1B_{u}^{-}$ , and ${3}A_{g}^{-}$ states. In particular, for \ceC24H26, the strong contraction errors amount to 0.09 to 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 dynamic correlation computed by sc-MPS-NEVPT2.

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 \ceC4H6 and \ceC8H10 vertical excitation energies using sc-NEVPT2 and t-NEVPT2 with the large aug-cc-pVTZ basis set and CASSCF reference wavefunctions (Table 4). Increasing the basis set from cc-pVDZ to aug-cc-pVTZ lowers the t-NEVPT2 relative energies of the $1B_{u}^{+}$ and ${2}A_{g}^{-}$ states in \ceC4H6 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.Watson:2012p4013 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 $1B_{u}^{+}$ state (6.06 and 5.83 eV) are too low relative to the best available theoretical value of 6.21 eV (Table 4), 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 \ceC8H10. In this case, both methods yield excitation energies for the covalent states ${2}A_{g}^{-}$ and $1B_{u}^{-}$ that are in the excellent agreement with the best available theoretical reference values,Angeli:2011p184302 while the relative energy of the ionic $1B_{u}^{+}$ state is underestimated by more than 0.7 eV. Similar errors for the $1B_{u}^{+}$ state have been observed by Angeli and Pastore in the study of \ceC8H10 using the second-order quasi-degenerate perturbation theory (QD-PT2) with the (8e,8o) active space.Angeli:2011p184302 They demonstrated that increasing the active space from (8e,8o) to (8e,16o) results in $\sim$ 0.7 eV increase in the $1B_{u}^{+}$ excitation energy computed using QD-PT2, leading to a closer agreement with the experiment.

The underestimation of the $1B_{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 (Figure 3), where we observe no crossing of $1B_{u}^{+}$ with the covalent ${2}A_{g}^{-}$ and $1B_{u}^{-}$ transitions, contrary to experimental results on substituted polyenesFuruichi:2002p547 ; Kohler:1988p5422 and some of calculations.Kurashige:2004p425 ; Schmidt:2012p124309 Re-investigating the trends in the excitation energies of the long polyenes with a double $\pi$ -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 (\ceC4H6 to \ceC24H26). For the chromium dimer, the active space used included the “double d” shell of the chromium atoms. Our \ceCr2 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.Guo:2016p1583 In a study of polyenes, using a complete $\pi$ -valence active space, our results demonstrate that including dynamic correlation is very important for the ionic $\pi$ - $\pi^{*}$ transitions, while it is less important for excitations with covalent character. In particular, for the \ceC24H26 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 $\pi$ -valence active space. In this case, using the “double” $\pi$ -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.

VI Acknowledgements

This work was supported by the U.S. Department of Energy (DOE), Office of Science through Award DE-SC0008624 (primary support for A.Y.S.), and Award 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.

VII 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 $\mathbf{G}\boldsymbol{{}^{[-2]}}(\tau)$ , which in the spin-orbital basis have the following form: $G_{z_{\kappa}w_{\lambda}}^{x_{\rho}y_{\sigma}}(\tau)=\braket{\Psi_{0}}{a^{\dagger}_{x_{\rho}}(\tau)a^{\dagger}_{y_{\sigma}}(\tau)a^{\ }_{w_{\lambda}}a^{\ }_{z_{\kappa}}}{\Psi_{0}}\equiv\braket{a^{\dagger}_{x\rho}(\tau)a^{\dagger}_{y\sigma}(\tau)a^{\ }_{w\lambda}a^{\ }_{z\kappa}}$ , where $x,y,w,z$ denote spatial orbitals and $\rho,\sigma,\kappa,\lambda$ denote the spin of these orbitals (e.g., $\rho=\alpha,\beta$ ). As discussed in Section II.3.3, in spin-adapted DMRG $a^{\dagger}_{x}$ and $a^{\ }_{x}$ are spin tensor creation and annihilation operators. We use a shorthand notation for the product of two tensor operators $\hat{A}_{xy}^{S,M}\equiv[a^{\dagger}_{x}a^{\dagger}_{y}]^{S}_{M}$ , 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}_{xy}^{S,M\ddagger}=(-1)^{S+M}\hat{A}_{xy}^{S,-M\dagger}$ . In addition, we define the spin-adapted tensor product $\otimes_{S}$ as $\hat{A}_{1}^{S_{1}}\otimes_{S}\hat{A}_{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}}]^{S}_{M}=\sum_{M_{1}M_{2}}c_{M_{1},M_{2},M}^{S_{1},S_{2},S}(\hat{A}_{1})^{S_{1}}_{M_{1}}(\hat{A}_{2})^{S_{2}}_{M_{2}}$ and $c_{M_{1},M_{2},M}^{S_{1},S_{2},S}$ are Clebsch-Gordan coefficients. Thus, the elements $G_{z_{\kappa}w_{\lambda}}^{x_{\rho}y_{\sigma}}(\tau)$ can be computed by solving the linear system of equations:

[TABLE]

Eq. 20 can also be used to compute elements of $\mathbf{G}\boldsymbol{{}^{[+2]}}(\tau)$ by replacing creation operators with annihilation operators and vice versa, e.g. $\braket{a^{\dagger}_{x\alpha}(\tau)a^{\dagger}_{y\beta}(\tau)a^{\ }_{w\beta}a^{\ }_{z\alpha}}\rightarrow\braket{a^{\ }_{x\alpha}(\tau)a^{\ }_{y\beta}(\tau)a^{\dagger}_{w\beta}a^{\dagger}_{z\alpha}}$ and $\braket{\hat{A}_{xy}^{1}(\tau)\otimes_{0}\hat{A}_{zw}^{1\ddagger}}\rightarrow\braket{\hat{A}_{yx}^{1\ddagger}(\tau)\otimes_{0}\hat{A}_{wz}^{1}}$ . To compute elements of $\mathbf{G}\boldsymbol{{}^{[0^{\prime}]}}(\tau)$ , we define a product of tensor operators $\hat{B}_{xy}^{S,M}\equiv[a^{\dagger}_{x}a^{\ }_{y}]^{S}_{M}$ . In this case, the system of linear equations has the form:

[TABLE]

Eqs. 20 and 21 can be simplified for a closed-shell reference wavefunction $\ket{\Psi_{0}}$ to obtain compact expressions:

[TABLE]

The remaining spin-orbital elements can be obtained using the following symmetry relations: $\braket{a^{\dagger}_{x\alpha}(\tau)a^{\dagger}_{y\alpha}(\tau)a^{\ }_{w\alpha}a^{\ }_{z\alpha}}=\braket{a^{\dagger}_{x\beta}(\tau)a^{\dagger}_{y\beta}(\tau)a^{\ }_{w\beta}a^{\ }_{z\beta}}$ , $\braket{a^{\dagger}_{x\alpha}(\tau)a^{\dagger}_{y\beta}(\tau)a^{\ }_{w\beta}a^{\ }_{z\alpha}}=\braket{a^{\dagger}_{x\beta}(\tau)a^{\dagger}_{y\alpha}(\tau)a^{\ }_{w\alpha}a^{\ }_{z\beta}}$ , $\braket{a^{\dagger}_{x\beta}(\tau)a^{\dagger}_{y\alpha}(\tau)a^{\ }_{w\beta}a^{\ }_{z\alpha}}=\braket{a^{\dagger}_{x\alpha}(\tau)a^{\dagger}_{y\beta}(\tau)a^{\ }_{w\alpha}a^{\ }_{z\beta}}$ , $\braket{a^{\dagger}_{x\alpha}(\tau)a^{\ }_{y\alpha}(\tau)a^{\dagger}_{w\alpha}a^{\ }_{z\alpha}}=\braket{a^{\dagger}_{x\beta}(\tau)a^{\ }_{y\beta}(\tau)a^{\dagger}_{w\beta}a^{\ }_{z\beta}}$ , $\braket{a^{\dagger}_{x\alpha}(\tau)a^{\ }_{y\beta}(\tau)a^{\dagger}_{w\beta}a^{\ }_{z\alpha}}=\braket{a^{\dagger}_{x\beta}(\tau)a^{\ }_{y\alpha}(\tau)a^{\dagger}_{w\alpha}a^{\ }_{z\beta}}$ , $\braket{a^{\dagger}_{x\alpha}(\tau)a^{\ }_{y\alpha}(\tau)a^{\dagger}_{w\beta}a^{\ }_{z\beta}}=\braket{a^{\dagger}_{x\beta}(\tau)a^{\ }_{y\beta}(\tau)a^{\dagger}_{w\alpha}a^{\ }_{z\alpha}}$ .

VIII 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^{\prime}$ , $-1^{\prime}$ , $0^{\prime}\}$ ) described in Section II.1, 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 $\ket{\Psi_{0}}$ . For example, in the $[+1^{\prime}]$ subclass, the perturber functions are $\ket{\Psi_{i}^{[+1^{\prime}]}}=\hat{h}^{{\dagger}}_{i}a_{i}\ket{\Psi_{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^{\prime}]$ subclass is

[TABLE]

where $\ket{\tilde{\Psi}_{i}^{[+1^{\prime}]}}=\frac{\ket{\Psi_{i}^{[+1^{\prime}]}}}{||\ket{\Psi_{i}^{[+1^{\prime}]}}||}$ , $E_{\tilde{\Psi}_{i}^{[+1^{\prime}]}}=\frac{\braket{\Psi_{i}^{[+1^{\prime}]}}{\hat{H}_{D}}{\Psi_{i}^{[+1^{\prime}]}}}{\braket{\Psi_{i}^{[+1^{\prime}]}}{\Psi_{i}^{[+1^{\prime}]}}}$ .

The main bottleneck in obtaining the sc-NEVPT2 energy is computing the term $\braket{\Psi_{i}^{[+1^{\prime}]}}{\hat{H}_{D}}{\Psi_{i}^{[+1^{\prime}]}}$ and the analogous contribution for $\ket{\Psi_{a}^{[-1^{\prime}]}}$ . This corresponds to

[TABLE]

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. Angeli:2001p10252 ; Angeli:2001p297 ; Angeli:2002ik

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 $\ket{\Psi_{i}^{[+1^{\prime}]}}$ and $\ket{\Psi_{a}^{[-1^{\prime}]}}$ as MPS of bond dimension $M_{1}$ , by carrying out a variational compression as described in Section II.3.2 to minimize $||\ket{\Psi_{i}^{[+1^{\prime}]}}-\hat{h}_{i}^{\dagger}a_{i}\ket{\Psi_{0}}||$ and $||\ket{\Psi_{a}^{[-1^{\prime}]}}-a_{a}^{{\dagger}}\hat{h}_{a}\ket{\Psi_{0}}||$ . Assuming $M_{1}\sim M_{0}$ , the cost of carrying out compressions for all the perturbers in class $[+1^{\prime}]$ and $[-1^{\prime}]$ is $\mathcal{O}(N_{ext}N^{2}_{act}M_{0}^{3})$ . The computational cost for the subsequent $\hat{H}_{D}$ expectation value is $\mathcal{O}(N_{ext}N^{3}_{act}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 $\mathcal{O}(N_{act}^{4}M_{0}^{3}+N_{act}^{6}M_{0}^{2}+N_{ext}N^{2}_{act}M_{0}^{3})$ . This is significantly lower than the $\mathcal{O}(N_{act}^{5}M_{0}^{3}+N_{act}^{8}M_{0}^{2})$ cost of computing the 4-RDM, which dominates the standard sc-NEVPT2 algorithm presented in Ref. 35.

Bibliography92

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) M. Schütz and H.-J. Werner, J. Chem. Phys. 114 , 661 (2001).
2(2) H.-J. Werner, F. R. Manby, and P. J. Knowles, J. Chem. Phys. 118 , 8149 (2003).
3(3) C. Riplinger and F. Neese, J. Chem. Phys. 138 , 034106 (2013).
4(4) C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, J. Chem. Phys. 139 , 134101 (2013).
5(5) S. R. White and R. L. Martin, J. Chem. Phys. 110 , 4127 (1999).
6(6) Ö. Legeza, R. Noack, J. Sólyom, and L. Tincani, in Computational Many-Particle Physics (Springer, 2008), pp. 653–664.
7(7) G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 131 , 054106 (2009).
8(8) Y. Kurashige and T. Yanai, J. Chem. Phys. 130 , 234114 (2009).

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Time-dependent NNN-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

Abstract

I Introduction

II Theory

II.1 Time-dependent formulation of NNN-electron valence perturbation theory (t-NEVPT2)

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

II.3 t-NEVPT2 with MPS reference wavefunctions (t-MPS-NEVPT2)

II.3.1 Overview of the algorithm

II.3.2 Green’s function evaluation

II.3.3 Spin-adaptation

III Computational Details

IV Results

IV.1 Analysis of numerical accuracy: N2 dissociation

IV.2 Chromium dimer

IV.3 Singlet excited states in all-trans polyenes

IV.3.1 Background

IV.3.2 Details of computations

IV.3.3 Discussion

V Conclusions

VI Acknowledgements

VII Appendix A: Spin-orbital two-body Green’s functions in spin-adapted DMRG

VIII Appendix B: Strongly-contracted NEVPT2 with MPS compression (sc-MPS-NEVPT2)

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

II.1 Time-dependent formulation of $N$ -electron valence perturbation theory (t-NEVPT2)