Driven similarity renormalization group: Third-order multireference perturbation theory

TL;DR

This paper introduces DSRG-MRPT3, a third-order multireference perturbation theory that is intruder free, size consistent, non-iterative, and includes reference relaxation, showing comparable accuracy to established methods.

Contribution

The paper presents the development and benchmarking of DSRG-MRPT3, a novel third-order multireference perturbation theory with improved accuracy and efficiency over previous methods.

Findings

DSRG-MRPT3 has nonparallelism errors similar to CASPT3 and MRCISD.

It shows significant improvements over DSRG second-order perturbation theory.

Accurately predicts singlet-triplet splitting within 0.1 kcal/mol of coupled cluster results.

Abstract

A third-order multireference perturbation theory based on the driven similarity renormalization group approach (DSRG-MRPT3) is presented. The DSRG-MRPT3 method has several appealing features: a) it is intruder free, b) it is size consistent, c) it leads to a non-iterative algorithm with ${\cal O}(N^6)$ scaling, and d) it includes reference relaxation effects. The DSRG-MRPT3 scheme is benchmarked on the potential energy curves of F$_2$, H$_2$O$_2$, C$_2$H$_6$, and N$_2$ along the F-F, O-O, C-C, and N-N bond dissociation coordinates, respectively. The nonparallelism errors of DSRG-MRPT3 are consistent to those of CASPT3 and MRCISD, and show significant improvements over those obtained from DSRG second-order multireference perturbation theory. Our efficient implementation of the DSRG-MRPT3 based on factorized electron repulsion integrals enables studies of medium-sized open-shell organic compounds. This point is demonstrated with computations of the singlet-triplet splitting ($\Delta_{\rm ST} = E_{\rm T} - E_{\rm S}$) of 9,10-anthracyne. At the DSRG-MRPT3 level of theory, our best estimate of the adiabatic $\Delta_{\rm ST}$ is 3.9 kcal mol$^{-1}$, a value that is within 0.1 kcal mol$^{-1}$ from multireference coupled cluster results.

Tables

Table 1. Table 1: Maximum error (MAX) and nonparallelism error (NPE) for the ground-state potential energy curves of \ce F-F, \ce HO-OH, \ce H3C-CH3, and \ce N ≡ \equiv N computed with various methods (reported in units of m E h subscript 𝐸 h E_{\rm h} ). All DSRG methods employ a value of the flow variable s = 0.5 𝑠 0.5 s=0.5 E h − 2 superscript subscript 𝐸 h 2 E_{\rm h}^{-2} . The last column shows the average NPE.

	\ceF2		\ceH2O2		\ceC2H6		\ceN2		Average
Method	MAX	NPE	MAX	NPE	MAX	NPE	MAX	NPE	NPE
pc-NEVPT2	$23.20$	$8.08$	$32.21$	$7.71$	$51.36$	$2.97$	$34.17$	$1.41$	$5.04$
CASPT2	$16.46$	$1.69$	$24.69$	$1.66$	$46.25$	$3.22$	$22.87$	$9.49$	$4.02$
uDSRG-MRPT2	$25.23$	$10.24$	$31.14$	$7.23$	$52.14$	$5.50$	$32.91$	$18.87$	$10.46$
DSRG-MRPT2	$20.33$	$5.34$	$29.38$	$5.46$	$51.06$	$4.42$	$32.35$	$18.34$	$8.39$
CASPT3	$9.29$	$5.38$	$12.51$	$5.23$	$13.76$	$3.08$	$7.77$	$6.31$	$5.00$
uDSRG-MRPT3	$9.46$	$5.03$	$15.50$	$6.04$	$19.15$	$2.73$	$15.19$	$6.05$	$4.96$
DSRG-MRPT3	$6.81$	$2.39$	$15.04$	$5.59$	$19.13$	$2.71$	$13.88$	$4.78$	$3.87$
MRCISD	$20.53$	$4.47$	$30.73$	$5.24$	$35.40$	$3.89$	$10.67$	$3.27$	$4.22$
MRCISD+Q	$1.60$	$1.35$	$3.02$	$1.94$	$4.29$	$2.27$	$1.96$	$1.50$	$1.77$
uMR-LDSRG(2)	$3.83$	$4.86$	$3.89$	$3.95$	$6.08$	$3.43$	$9.30$	$5.25$	$4.37$
MR-LDSRG(2)r1	$-1.03$	$1.34$	$0.78$	$0.86$	$5.28$	$2.87$	$8.25$	$4.98$	$2.51$
MR-LDSRG(2)	$-5.19$	$4.24$	$-3.87$	$3.75$	$4.92$	$3.12$	$7.61$	$4.81$	$3.98$

Table 2. Table 2: Adiabatic singlet-triplet splitting ( Δ ST = E T − E S subscript Δ ST subscript 𝐸 T subscript 𝐸 S \Delta_{\rm ST}=E_{\rm T}-E_{\rm S} , in kcal mol -1 ) of 9,10-anthracyne computed using various methods with the cc-pVTZ basis set. All values are corrected by the zero-point vibrational energy ( + 0.264 0.264 +0.264 kcal mol -1 ) obtained at the CASSCF(2,2)-DSRG-MRPT3/cc-pVDZ level of theory.

	CAS
Method	(2,2)	(4,4)	(8,8)	(12,12)	(16,16)
CASSCF	$0.61$	$2.99$	$4.52$	$2.64$	$2.42$
CASPT2	$4.08$	$5.47$	$3.24$	$4.78$
sc-NEVPT2	$2.90$	$3.79$	$0.78$	$4.02$
pc-NEVPT2	$3.08$	$3.86$	$0.90$	$4.01$
uDSRG-MRPT2a	$2.15$	$3.82$	$3.20$	$3.72$	$3.66$
DSRG-MRPT2a	$2.58$	$3.89$	$3.21$	$3.86$	$3.78$
DSRG-MRPT2b	$3.19$	$4.00$	$2.52$	$3.96$	$3.85$
uDSRG-MRPT3a	$1.74$	$3.50$	$5.36$	$3.46$	$3.40$
DSRG-MRPT3a	$1.98$	$3.70$	$5.50$	$3.66$	$3.61$
DSRG-MRPT3b	$2.38$	$3.86$	$5.58$	$3.93$	$3.90$
MRCISD	$1.23$	$3.23$	$4.30$	$3.00$
MRCISD+Q	$1.71$	$3.41$	$4.23$	$3.25$
Mk-MRCCSD	$4.44$
Mk-MRCCSD(T)	$3.98$

Equations

$$|{\Psi }_0⟩=\sum _{\mu =1}^d{c}_{\mu }|{\Phi }^{\mu }⟩,$$

$$\hat{H}\to \bar{H}(s)={\hat{U}}^{†}(s)\hat{H}\hat{U}(s),s\ge 0.$$

$${\bar{H}}^{\mathrm{N}}(s)=\hat{R}(s).$$

$$\hat{U}(s)=\exp [\hat{A}(s)]=\exp [\hat{T}(s)-{\hat{T}}^{†}(s)].$$

$$\bar{H}(s)=\hat{H}+\sum _{k=1}^{\infty }\frac{1}{k!}\underset{k\text{ nested commutators}}{\underbrace{[\cdots [[\hat{H},\hat{A}(s)],\hat{A}(s)],\cdots ]}}.$$

$$E_{\mathrm{u}}(s)=⟨{\Psi }_0⟩\bar{H}(s){\Psi }_0.$$

$$\sum _{\mu }^d⟨{\Phi }_{\nu }⟩\bar{H}(s){\Phi }^{\mu }{c}_{\mu }=E(s)c_{\nu }.$$

$${\bar{H}}^{(0)}(s)=$$

$${\bar{H}}^{(1)}(s)=$$

$${\bar{H}}^{(2)}(s)=$$

$${\bar{H}}^{(3)}(s)=$$

$$+\frac{1}{2}[{\tilde{H}}^{(2)}(s),{\hat{A}}^{(1)}(s)],$$

$${\hat{R}}^{(1)}(s)$$

$${\hat{R}}^{(2)}(s)$$

$${\hat{H}}^{(0)}=E_0+{\hat{F}}^{(0)},$$

$${\hat{F}}^{(0)}$$

$$=\sum _p{ϵ}_p\{{\hat{a}}_p^p\}.$$

$$f_p^q=h_p^q+\sum _{ij}{v}_{pi}^{qj}{\gamma }_j^i.$$

$$E_{\mathrm{u}}^{(0)}(s)=$$

$$E_{\mathrm{u}}^{(1)}(s)=$$

$$E_{\mathrm{u}}^{(2)}(s)=$$

$$E_{\mathrm{u}}^{(3)}(s)=$$

$$+⟨{\Psi }_0⟩[{\tilde{H}}^{(2)}(s),{\hat{T}}^{(1)}(s)]{\Psi }_0,$$

$$\begin{array}{rl}{\bar{H}}^{(3)}(s)\approx & [{\hat{H}}^{(0)},{\hat{A}}^{(3)}(s)]+\frac{1}{2}[{\tilde{H}}^{(1)}(s),{\hat{A}}^{(2)}(s)]\\ & +\frac{1}{2}[{\tilde{H}}_{1,2}^{(2)}(s),{\hat{A}}^{(1)}(s)],\end{array}$$

$$t_a^{i,(1)}(s)$$

$$t_{ab}^{ij,(1)}(s)$$

$$t_a^{i,(2)}(s)$$

$$t_{ab}^{ij,(2)}(s)$$

$$E_{\mathrm{u}}^{[3]}=\sum _{n=0}^3{E}_{\mathrm{u}}^{(n)}.$$

$$\sum _{\mu }^d⟨{\Phi }_{\nu }⟩{\bar{H}}_{1,2}^{[3]}(s){\Phi }^{\mu }{c}_{\mu }=E^{[3]}(s)c_{\nu },$$

$${\hat{C}}_k(s)=\frac{1}{k}[{\hat{C}}_{k-1}(s),\hat{A}(s){]}_{1,2},k=1,2,3,\cdots ,$$

$${[{\tilde{H}}^{(1)}(s)]}_i^a=$$

$${[{\tilde{H}}^{(1)}(s)]}_{ij}^{ab}=$$

$$(pq|rs)\approx \sum _Q^M{B}_{pq}^Q{B}_{rs}^Q,$$

$${\bar{H}}_{ef}^{ij,(2)}(s)\leftarrow \sum _{gh}^{\mathbf{V}}{v}_{ef}^{gh,(1)}{t}_{gh}^{ij,(1)}(s).$$

$$E_{\mathrm{u}}^{(2)}\leftarrow \frac{1}{4}\sum _{\genfrac{}{}{0ex}{}{xyz}{uvw}}^{\mathbf{A}}\sum _e^{\mathbf{V}}{v}_{xy}^{we,(1)}(s)t_{ez}^{uv,(1)}(s){\lambda }_{uvw}^{xyz},$$

$$\text{NPE}=\underset{r\in \mathbf{R}}{\max }\Delta E(r)-\underset{r\in \mathbf{R}}{\min }\Delta E(r),$$

Full text

I Introduction

Multireference perturbation theory (MRPT) based on a complete active space (CAS)Roos, Taylor, and Siegbahn (1980); Ruedenberg et al. (1982); Werner and Knowles (1985) wave function is one of the simplest and most popular quantum chemistry approaches for studying near-degenerate electronic states. A number of MRPTs have been proposed, among which the second-order complete active space perturbation theory (CASPT2) by Roos and co-workers is perhaps the most successful one.Andersson et al. (1990); Andersson, Malmqvist, and Roos (1992) The partially- and strongly-contracted variants of second-order $n$-electron perturbation theory (pc- and sc-NEVPT2), introduced more recently, are also growing in popularity.Dyall (1995); Angeli et al. (2001); Angeli, Cimiraglia, and Malrieu (2002); Angeli, Pastore, and Cimiraglia (2007) One of the advantages of CASPT2, pc- and sc-NEVPT2 is that they employ an internally contracted formalism, whereby the first-order correction to the wave function is generated by excitation operators acting on the entire reference wave function. Therefore, the computational cost of these methods scales only polynomially with the number of active space orbitals—a tremendous reduction in cost compared to exponentially-scaling uncontracted multireference formalisms.

Despite their success, CASPT2 and NEVPT2 have some crucial limitations. Perhaps the most unsettling feature of CASPT2—and most other MRPT methods—is the intruder-state problem.Evangelisti, Daudey, and Malrieu (1987); Paldus et al. (1993); Kowalski and Piecuch (2000); Roos and Andersson (1995) Intruder states are encountered when the zeroth-order energy of the reference and excited configurations are near degenerate. Intruder states introduce singularities in the energy denominators and yield excitation amplitudes with unphysically large values. Several approaches have been proposed to address intruders,Dyall (1995); Andersson (1995); Ghigo, Roos, and Malmqvist (2004); Witek et al. (2002); Taube and Bartlett (2009); Camacho, Witek, and Yamamoto (2009) yet the most straightforward is level (denominator) shifting.Roos and Andersson (1995); Forsberg and Malmqvist (1997) When intruder states are weakly coupled to the reference, level shifting works well for both ground and excited states. However, in computations on multiple excited states, finding a ubiquitous level shift that removes intruders for all states may be challenging.Roos et al. (1996) In addition, level shifting introduces some arbitrariness in CASPT2 results.Camacho, Witek, and Yamamoto (2009); Vancoillie, Malmqvist, and Veryazov (2016) Another less severe problem is the small size-consistency error carried by CASPT2 due to the use of projectors in the zeroth-order Hamiltonian [$\hat{H}^{(0)}$].Rintelman et al. (2005); Helgaker, Jørgensen, and Olsen (2000) These projectors are introduced because the reference is not an eigenfunction of the average Fock operator.

In NEVPT2 the above two issues are neatly solved by employing Dyall’s zeroth-order Hamiltonian that includes two-electron interactions within the active orbitals.Dyall (1995) However, both CASPT2 and NEVPT2 have computational bottlenecks that prevent computations with large active spaces. For a complete active space, the energy expressions of CASPT2 and NEVPT2 demand the four-particle reduced density matrix (4RDM) of the reference. The memory cost of storing this quantity grows as the eighth power of the number of active orbitals ($N_{\bf A}$) and quickly becomes the Achilles’ heel of these methods when $N_{\bf A}>16$. The computational cost of CASTP2 and NEVPT2 is then dominated by tensor contractions that scale as ${\cal O}(N_{\bf A}^{8})$ and ${\cal O}(N_{\bf A}^{9})$, respectively. Several attempts have been made to avoid the 4RDM via cumulant decompositions,Kutzelnigg and Mukherjee (1997, 1999); Yanai and Chan (2007) and encouraging results have been obtained for CASPT2.Kurashige and Yanai (2011) Further approximations to the three-particle density matrix are less promising as “false intruders” appear on the potential energy curves.Kurashige and Yanai (2011); Zgid et al. (2009) Unfortunately, even when the 4RDM is neglected, CASPT2 and pc-NEVPT2 still require removing linear dependencies in the excitation manifold by diagonalizing the overlap metric.Andersson, Malmqvist, and Roos (1992); Angeli et al. (2001) This step also has a computational costs proportional to ${\cal O}(N_{\bf A}^{9})$ and restricts state-of-the-art CASPT2 computations to $N_{\bf A}\approx 30$.Kurashige and Yanai (2011)

Multireference perturbation theories based on the driven similarity renormalization groupEvangelista (2014) (DSRG) provide a solution to both the intruder-state problem and the computational scaling limitations of CASPT2 and NEVPT2. The DSRG is a many-body approach closely related to the in-medium similarity renormalization group (IM-SRG),Głazek and Wilson (1993); Wegner (2000); Tsukiyama, Bogner, and Schwenk (2011); Hergert et al. (2013, 2016) coupled cluster,Čížek (1966); Bartlett (1981); Crawford and Schaefer (2000); Bartlett and Musiał (2007) and canonical transformation theories.Yanai and Chan (2006, 2007); Neuscamman, Yanai, and Chan (2010) Like in the IM-SRG, the DSRG separates excitation energy scales via a continuous unitary transformation of the Hamiltonian controlled by a flow parameter $s$. A perturbative analysis of the DSRG shows that this transformation folds in correlation effects from excited configurations that correspond to energy denominators larger than a cutoff $\Lambda=s^{-1/2}$, while it leaves untouched those excitations for which the denominators are smaller than $\Lambda$.Kehrein (2006); Hergert (2017) As such, the DSRG avoids intruders for finite values of $s$ and yields a transformed (renormalized) Hamiltonian with modified many-body interactions. Nonetheless, the DSRG is distinct from the IM-SRG. While the IM-SRG directly determines the renormalized Hamiltonian by solving a collection of ordinary differential equations, the DSRG obtains it from a set of coupled nonlinear equations. Another critical ingredient of the DSRG is the use of Fock-space many-body conditions,Kutzelnigg (1982); *Kutzelnigg:1983dr; *Kutzelnigg:1984eg; *Kutzelnigg:1985fj; Stolarczyk and Monkhorst (1985a); *Stolarczyk:1985ct; *Stolarczyk:1988ci; *Stolarczyk:1988cv; Lindgren (1978); Nooijen and Bartlett (1996) which lead to equations in terms of normal-ordered second-quantized operators. For multireference theories, the many-body approach effectively avoids the need to orthogonalize the excitation manifold.Datta, Kong, and Nooijen (2011); Datta and Nooijen (2012); Li and Evangelista (2015, 2016)

In our previous work, we have explored the multireference DSRGLi and Evangelista (2015, 2016) (MR-DSRG) and its second-order perturbation theory (DSRG-MRPT2).Li and Evangelista (2015); Hannon, Li, and Evangelista (2016) The MR-DSRG formalism is built upon the algebra of Mukherjee and Kutzelnigg’s generalized normal ordering and Wick’s theorem,Kutzelnigg and Mukherjee (1997); Mukherjee (1997); Mahapatra et al. (1998); Shamasundar (2009); Kong, Nooijen, and Mukherjee (2010); Kutzelnigg, Shamasundar, and Mukherjee (2010); Sinha, Maitra, and Mukherjee (2013) where operator contractions lead to density cumulants.Kutzelnigg and Mukherjee (1997); Kutzelnigg, Shamasundar, and Mukherjee (2010); Hanauer and Köhn (2012) The advantage of this scheme is that the DSRG-MRPT2 energy requires at most the three-particle density cumulant. Consequently, the DSRG-MRPT2 approach has a scaling that is proportional to $\mathcal{O}(N_{\bf A}^{6})$, and could potentially be applied to systems with large active spaces.

Our previous work has shown that the accuracy of the DSRG-MRPT2 is similar to that of other second-order MRPTs.Li and Evangelista (2015) Furthermore, the linearized MR-DSRG with one- and two-body operators [MR-LDSRG(2)] greatly improves upon the accuracy of DSRG-MRPT2, but requires a recursive evaluation of the Hamiltonian and an iterative update of the cluster amplitudes.Li and Evangelista (2016) The cost of MR-LDSRG(2) computations currently limits applications of this method to systems with 200–300 orbitals. In this work we propose to overcome this limitation by developing a third-order multireference perturbation theory (MRPT3) based on the DSRG. A noniterative DSRG-MRPT3 would be less expensive than the MR-LDSRG(2) scheme and likely to be more accurate than second-order perturbation theory. Indeed, several third-order MRPTs have been formulated and they were found to be superior with respect to the corresponding second-order MRPTs.Wolinski, Sellers, and Pulay (1987); Werner (1996); Angeli and Cimiraglia (2002); Angeli et al. (2006); Jiang, Khait, and Hoffmann (2006); Khait, Jiang, and Hoffmann (2009) For example, the CASPT3 implementation of Werner provides geometries and harmonic frequencies of small molecules that are as accurate as those of multireference configuration interaction with singles and doubles (MRCISD), but only costs as one iteration of MRCISD.Werner (1996)

In this work, we derive and implement a third-order DSRG-MRPT (DSRG-MRPT3). The zeroth-order Hamiltonian is chosen to contain only the diagonal blocks of the Fock operator, an identical choice made in the CASPT2D approach.Andersson et al. (1990) When applied to the DSRG-MRPT3, this choice of $\hat{H}^{(0)}$ leads to an efficient non-iterative formalism that is free from the intruder-state problem,Li and Evangelista (2015) is rigorously size extensive,Bartlett (1981) and may be interfaced with any reference wave function for which the one-, two-, and three-body density cumulants are computed. In this study we also develop relaxed DSRG-MRPT2 and MRPT3 approaches in which the reference wave function is optimized under the effects of dynamic electron correlation. Most state-specific MRPT2 approaches, including CASPT2 and NEVPT2, do not account for reference relaxation effects, with the notable exception of Mukherjee’s state-specific MRPT2,Mahapatra, Datta, and Mukherjee (1999a, 1998); Sinha Mahapatra, Datta, and Mukherjee (1999); Chaudhuri et al. (2005); Chattopadhyay et al. (2016) generalized Van Vleck PT2,Shavitt and Redmon (1980); Kirtman (1981); Khait, Song, and Hoffmann (2002); Hoffmann et al. (2009) and multiconfigurational PT2.Rolik, Szabados, and Surján (2003); Szabados et al. (2005); Chaudhuri et al. (2005)

We begin our discussion of DSRG-MRPT3 by providing a brief overview of the general MR-DSRG ansatz in Sec. II.1, followed in Sec. II.2 by a detailed perturbative analysis. Then, in Sec. III, we provide details of the DSRG-MRPT3 implementation in our open-source code Forte.FOR (2016) Section V consists of two parts. In the first one we study the potential energy curves of \ceF2, \ceHO-OH, \ceH3C-CH3, and \ceN2 to assess the accuracy of the DSRG-MRPT3 method and compare its performance with that of the DSRG-MRPT2 and the MR-LDSRG(2) approaches. In the second part we compute the singlet-triplet gap of 9,10-anthracyne and compare DSRG-MRPT3 results with other multireference methods including CASPT2,Celani and Werner (2000) MRCISD,Knowles and Werner (1988); Werner and Knowles (1988); Shamasundar, Knizia, and Werner (2011) and Mukherjee multireference coupled cluster theory.Mahapatra et al. (1998); Mahapatra, Datta, and Mukherjee (1999b); Evangelista, Allen, and Schaefer (2006); Mahapatra, Datta, and Mukherjee (2010); Evangelista et al. (2010) Finally, in Sec. VI we conclude and discuss future applications of the DSRG-MRPT3 scheme.

II Theory

II.1 An overview of MR-DSRG

In this section we briefly review the MR-DSRG formalism. Readers who are interested in the details regarding the operator parameterization should refer to the original DSRG (Ref. 29) and MR-DSRG (Refs. 55 and 56) papers. In MR-DSRG theory, we employ the generalized normal ordering of Mukherjee and Kutzelnigg (MK-GNO) to deal with the algebra of second-quantized operators.Kutzelnigg and Mukherjee (1997); Mukherjee (1997); Mahapatra et al. (1998); Shamasundar (2009); Kong, Nooijen, and Mukherjee (2010); Kutzelnigg, Shamasundar, and Mukherjee (2010); Sinha, Maitra, and Mukherjee (2013) Under MK-GNO, operator contractions are associated with density cumulants,Kutzelnigg and Mukherjee (1997); Kutzelnigg, Shamasundar, and Mukherjee (2010); Hanauer and Köhn (2012) which embody all information of the reference wave function. The reference wave function used to define the MK-GNO Fermi vacuum is a multideterminantal wave function:

[TABLE]

where each determinant $\Phi^{\mu}$ is weighted by the coefficient $c_{\mu}$. In this work we further assume that the set of determinants $\{\Phi^{\mu},\mu=1,2,\dots,d\}$ forms a complete active space (CAS), although this is not generally required by the MR-DSRG formalism. The coefficients $c_{\mu}$ and the molecular orbitals are determined by the CAS self-consistent field (CASSCF) procedure.Roos, Taylor, and Siegbahn (1980); Ruedenberg et al. (1982) The set of spin orbitals $\{\phi^{p},p=1,2,\dots,N\}$ then falls into three subsets: core ($\bf C$, with indices $m,n$), active ($\bf A$, with indices $u,v,w,x,y,z$), and virtual ($\bf V$, with indices $e,f,g,h$) of dimension $N_{\bf C}$, $N_{\bf A}$, and $N_{\bf V}$, respectively. Two composite orbital sets are introduced: hole ($\bf H=C\cup A$, with indices $i,j,k,l$) and particle ($\bf P=A\cup V$, with indices $a,b,c,d$) of size $N_{\bf H}=N_{\bf C}+N_{\bf A}$ and $N_{\bf P}=N_{\bf A}+N_{\bf V}$, respectively. General orbitals ($\bf G$) are denoted by indices $p,q,r,s$.

The philosophy of the DSRG ansatz is to define a continuous ($s$-dependent) unitary operator $\hat{U}(s)$ that transforms the bare Hamiltonian $\hat{H}$ to a band-diagonal operator $\bar{H}(s)$, namely:

[TABLE]

When the flow variable $s$ approaches infinity, we require this transformation to exactly zero the couplings between the reference state and its internally-contracted excited configurations. In the many-body formalism,Lindgren (1978); Nooijen and Bartlett (1996); Datta, Kong, and Nooijen (2011); Datta and Nooijen (2012); Demel, Datta, and Nooijen (2013) this coupling is conveniently represented by the non-diagonal terms of the Hamiltonian $[\bar{H}^{\rm N}(s)]$. The DSRG assumes that the unitary transformation is controlled by a Hermitian source operator $\hat{R}(s)$. The corresponding many-body DSRG flow equationEvangelista (2014); Li and Evangelista (2015, 2016) realizes this idea by equating the non-diagonal terms of the Hamiltonian to the source operator:

[TABLE]

The source operator $\hat{R}(s)$ renormalizes the Hamiltonian in such a way that excited configurations that are energetically separated from the reference by at least $s^{-1/2}$ are decoupled to each other.Kehrein (2006); Evangelista (2014) Note, the parametrization of $\hat{R}(s)$ that can achieve this renormalization transformation is not unique. In this work we use a source operator (see Ref. 29) that reproduces the transformation of second-order perturbation theory based on the single-reference similarity renormalization group.Tsukiyama, Bogner, and Schwenk (2011); Hergert et al. (2016)

Once $\hat{R}(s)$ is defined, the unitary operator $\hat{U}(s)$ can be determined via Eqs. (2) and (3). In order to set up Eq. (3), we need to write the DSRG transformed Hamiltonian $\bar{H}(s)$ as a sum of second-quantized operators. The unitary operator $\hat{U}(s)$ is expressed as the exponential of an anti-Hermitian operator $\hat{A}(s)$, and $\hat{A}(s)$ is further related to the coupled cluster excitation operator $\hat{T}(s)$,

[TABLE]

Note that internal amplitudes ${t}_{uv\cdots}^{xy\cdots}(s)$ with $u,v,\dots,x,y,\dots\in\bf A$ are redundant and hence set to zero. Using the Baker–Campbell–Hausdorff (BCH) expansion, we write the renormalized Hamiltonian $\bar{H}(s)$ as a series of commutators of $\hat{H}$ and $\hat{A}(s)$,

[TABLE]

The many-body expression of $\bar{H}(s)$ is obtained by evaluating the commutators using the MK-GNO Wick’s theorem and subsequently collecting the same rank of normal-ordered operators.Evangelista (2014); Li and Evangelista (2015, 2016) Because the BCH formula in Eq. (5) does not terminate, we truncate each commutator $[\hat{H},\hat{A}(s)]$ to contain at most two-body operators,Yanai and Chan (2006, 2007); Li and Evangelista (2016) and consider the series converged when the Frobenius norm of the $k$-nested commutator is less than $10^{-12}$ $E_{\rm h}$.Evangelista and Gauss (2012); Evangelista (2014); Li and Evangelista (2016) We shall use the subscript “1,2” whenever the operator is truncated to contain at most two-body terms.

The MR-DSRG energy can be computed in three ways. The first approach simply computes the expectation value of $\bar{H}(s)$ using the reference wave function $\Psi_{0}$,

[TABLE]

The energy obtained via Eq. (6) is unrelaxed because the weight of each reference determinant ($c_{\mu}$) is fixed and equal to the reference CASCI/CASSCF wave function. Alternatively, the coefficient $c_{\mu}$ may be optimized by diagonalizing $\bar{H}(s)$ within the set of reference determinants $\{\Phi^{\mu}\}$:

[TABLE]

We consider two approaches to treat reference relaxation effects. The first one consists in the iterative solutions of Eqs. (7) and (3) until self consistency is reached. Achieving self consistency of the amplitude and eigenvalue equations is a necessary condition to guarantee that when $\hat{A}(s)$ is not truncated, the BCH expansion of $\bar{H}(s)$ is not approximated, and $s\rightarrow\infty$, then the full MR-DSRG is equivalent to full configuration interaction. In the second approach, we approximate the relaxed energy with a single-step correction which consist in solving Eq. (7) only once and take the eigenvalue of the matrix $\braket{\Phi_{\nu}}{\bar{H}(s)}{\Phi^{\mu}}$ as the relaxed energy. This single-step relaxation scheme is economical since it requires only one $N^{6}$ step. Both relaxation schemes are compared and benchmarked in Sec. V.

II.2 The DSRG-MRPT3 method

We first partition the normal-ordered bare Hamiltonian into a zeroth-order term $\hat{H}^{(0)}$ and a first-order fluctuation potential $\hat{H}^{(1)}$. The DSRG perturbation theory is derived from an order-by-order expansion of the source operator $\hat{R}(s)$, the anti-Hermitian operator $\hat{A}(s)$, and the unrelaxed energy $E_{\rm u}(s)$, while the coefficients $c_{\mu}$ in Eq. (7) are not treated perturbatively.Li and Evangelista (2015) The zeroth- till third-order DSRG transformed Hamiltonian are given by

[TABLE]

where we have introduced the combined first-order Hamiltonian $\tilde{H}^{(1)}(s)=\hat{H}^{(1)}+\bar{H}^{(1)}(s)$ and the second-order counterpart $\tilde{H}^{(2)}(s)=\bar{H}^{(2)}(s)-\frac{1}{6}[[\hat{H}^{(0)},\hat{A}^{(1)}(s)],\hat{A}^{(1)}(s)]$. The first- and second-order cluster amplitudes are determined by:

[TABLE]

To simplify the structure of Eqs. (8)–(11), the zeroth-order Hamiltonian $\hat{H}^{(0)}$ is chosen to include the reference energy and the diagonal blocks of the Fock operator

[TABLE]

where

[TABLE]

Here the string of normal-ordered creation ($\hat{a}^{\dagger}{}$) and annihilation ($\hat{a}{}$) operators are compactly written as $\{{\hat{a}_{rs\dots}^{pq\dots}}\}=\{{\hat{a}^{\dagger}_{p}\hat{a}^{\dagger}_{q}\dots\hat{a}_{s}\hat{a}_{r}}\}$. The generalized Fock matrix element ${f}_{p}^{q}$ is expressed in terms of one-electron (${h}_{p}^{q}$) and antisymmetrized two-electron (${v}_{pq}^{rs}$) integrals as well as the one-particle density matrix $\gamma_{q}^{p}=\braket{\Psi_{0}}{\hat{a}^{\dagger}_{p}\hat{a}_{q}}{\Psi_{0}}$,

[TABLE]

As indicated by the last line of Eq. (II.2), molecular orbitals are semicanonicalized such that the core, active and virtual blocks of the generalized Fock matrix are diagonal ($\epsilon_{p}={f}_{p}^{p}$).

With this zeroth-order Hamiltonian, the commutator $[\hat{H}^{(0)},\hat{A}^{(n)}(s)]$ yields only one- and two-body non-diagonal terms (that is, $\{{\hat{a}_{i}^{a}}\}$, $\{{\hat{a}_{ij}^{ab}}\}$, $\{{\hat{a}_{a}^{i}}\}$, and $\{{\hat{a}_{ab}^{ij}}\}$). As a result, expectation values of the form $\braket{\Psi_{0}}{[\hat{H}^{(0)},\hat{A}^{(n)}(s)]}{\Psi_{0}}=0$ for any $n$. The unrelaxed zeroth- till third-order energies written in terms of the cluster operators are thus:

[TABLE]

where the prefactors that enter in Eqs. (10) and (11) cancel since $\braket{\Psi_{0}}{[\hat{H},\hat{A}(s)]}{\Psi_{0}}=2\braket{\Psi_{0}}{[\hat{H},\hat{T}(s)]}{\Psi_{0}}$.

Note that the contribution $\braket{\Psi_{0}}{[\tilde{H}^{(2)}(s),\hat{T}^{(1)}(s)]}{\Psi_{0}}$ to $E^{(3)}_{\rm u}(s)$ reported in Eq. (20) contains contractions that involve the four-body cumulant of the reference. To avoid the cost of computing and storing the four-body cumulant we approximate the third-order transformed Hamiltonian with:

[TABLE]

where $\tilde{H}^{(2)}_{1,2}(s)$ includes only the one- and two-body components of the operator $\tilde{H}^{(2)}(s)$ and neglects the remaining three-body components. This approximation is consistent with the linear MR-DSRG truncated to two-body operators [MR-LDSRG(2)].Li and Evangelista (2016) By this we mean that if the perturbative series for $\bar{H}^{(n)}$ were to be resummed to infinite order, consistently truncating all commutators to one- and two-body operators, it would yield (if convergent) the unrelaxed MR-LDSRG(2) transformed Hamiltonian.

From Eq. (12), we can derive the explicit expressions of the first-order amplitudes,Li and Evangelista (2015)

[TABLE]

where ${\Delta}_{ab\dots}^{ij\dots}=\epsilon_{i}+\epsilon_{j}+\dots-\epsilon_{a}-\epsilon_{b}-\dots$ is a general Møller–Plesset denominator. Analogous expressions for the second-order amplitudes can be derived from Eq. (13),

[TABLE]

Here ${f}_{a}^{i,(2)}(s)$ and ${v}_{ab}^{ij,(2)}(s)$ are one- and two-body non-diagonal elements of $\frac{1}{2}[\tilde{H}^{(1)}(s),\hat{A}^{(1)}(s)]$, respectively.

Once the first- and second-order amplitudes are known, the unrelaxed DSRG-MRPT3 energy is obtained by summing the zeroth- to third-order scalar terms of Eqs. (17)–(20),

[TABLE]

The relaxed energy is obtained by a one-step diagonalization [Eq. (7)] using the summed third-order DSRG transformed Hamiltonian:

[TABLE]

where we have applied the operator truncation $\bar{H}^{[3]}_{1,2}(s)=\sum_{n=0}^{3}\bar{H}^{(n)}_{1,2}(s)$. As such, both relaxed and unrelaxed DSRG-MRPT3 energies are computed in a non-iterative fashion. In Appendix A, we show that solving Eq. (26) only requires minuscule amount of work comparing to computing Eq. (25). For brevity, we shall take the energy of Eq. (26) as our DSRG-MRPT3 energy, otherwise as uDSRG-MRPT3 for the unrelaxed energy.

We note that other more sophisticated zeroth-order Hamiltonian such as the one proposed by DyallDyall (1995) or the retaining-excitation HamiltonianFink (2006, 2009) can also be straightforwardly applied to the DSRG perturbation theory. However as discussed in Appendix B, using these Hamiltonians mostly lead to unnecessary complications and bring only small advantages.

III Implementation

To obtain the DSRG-MRPT3 working equations, the core task is to evaluate the commutators in Eqs. (9)–(11). Essentially the $k$-th nested commutator [$\hat{C}_{k}(s)$] is computed from the following recursive relation:

[TABLE]

where $\hat{C}_{0}(s)=\hat{H}$. In Ref. 56 we have reported all terms of $\hat{C}_{1}(s)=[\hat{H},\hat{A}(s)]_{1,2}$ and we shall not repeat them here. Instead, we would like to discuss those terms that are necessary to implement the DSRG-MRPT3 scheme.

Computing the scalar term of commutator $[\hat{H},\hat{T}(s)]$ for example requires only the non-diagonal terms of $\hat{H}$ (${f}_{i}^{a}$ and ${v}_{ij}^{ab}$). This statement suggests that we need elements such as ${[\tilde{H}^{(n)}(s)]}_{i}^{a}$ and ${[\tilde{H}^{(n)}(s)]}_{ij}^{ab}$ for $n=1,2$ to compute the second- and third-order energies. For $n=1$ we are able to write out explicit expressions due to the simple structure of $\hat{H}^{(0)}$:

[TABLE]

When $n=2$ we call for the non-diagonal elements of $\frac{1}{2}[\tilde{H}^{(1)}(s),\hat{A}^{(1)}(s)]_{1,2}$, computing which necessitate all elements of $\tilde{H}^{(1)}(s)$. Fortunately except for non-diagonal terms of $\tilde{H}^{(1)}(s)$, the remainders are equivalent to those in the first-order bare Hamiltonian. We therefore factorize the two-electron integrals using density fitting (DF) or Cholesky decomposition (CD) techniquesWhitten (1973); Dunlap, Connolly, and Sabin (1979); Werner, Manby, and Knowles (2003); Beebe and Linderberg (1977); Koch, Sánchez de Merás, and Pedersen (2003); Aquilante et al. (2008, 2011); Weigend, Kattannek, and Ahlrichs (2009); DePrince and Sherrill (2013); Hannon, Li, and Evangelista (2016) to avoid storing the four-index tensor of size $N^{4}$. Specifically two-electron integrals are approximated as a contraction of two three-index tensors:

[TABLE]

where $Q$ is the auxiliary index. The quantity $M$ is maximum value of $Q$ and it is roughly three times as large as the regular basis set ($M\approx 3N$).

For small active spaces ($N_{\bf A}\ll N_{\bf C}$) the asymptotic scaling of the DSRG-MRPT3 method is ${\cal O}(N_{\bf V}^{4}N_{\bf C}^{2})$, resulting from the contraction:

[TABLE]

This tensor contraction is also contained in the coupled cluster equations and it is a scourge because its computational cost cannot be reduced by factorizing the two-electron integrals ${v}_{ef}^{gh,(1)}$ via DF or Cholesky decomposition,Schutz and Manby (2003) instead it requires the use of alternative factorizations.Parrish et al. (2014) As a compromise, we form ${v}_{ef}^{gh,(1)}$ in batches of compound indices $gh$ according to Eq. (30), with the size of each batch automatically determined by available memory. For large active spaces, the computational cost of the DSRG-MRPT3 is constrained by terms such as:

[TABLE]

where $\lambda_{uvw}^{xyz}$ is the three-particle density cumulant of the reference. Although the term given in Eq. (32) scales as ${\cal O}(N_{\bf A}^{6}N_{\bf V})$, the effort made to compute it is still significantly less than the most expensive step in CASPT2 or NEVPT2 [${\cal O}(N_{\bf A}^{9})$].

The DSRG-MRPT3 method is implemented in our open-source code Forte,FOR (2016) a plugin to the Psi4 packageTurney et al. (2012) that specializes on multireference methods. All tensor contractions are written using the syntax provided by the open-source tensor library Ambit.AMB (2016) The three-index DF integrals generated by Psi4 are read using a general interface developed by Hannon et al.Hannon, Li, and Evangelista (2016) The DSRG-MRPT3 equations are spin integrated such that one-body terms are decomposed into $\alpha$ and $\beta$ blocks and two-body terms are divided into $\alpha\alpha$, $\alpha\beta$, and $\beta\beta$ contributions. For convenience, the current implementation constantly stores nine tensors of size $N_{\bf H}^{2}N_{\bf P}^{2}/16$ coming from the spin integration of the $t_{2}$ amplitudes (either first- or second-order) and two two-body intermediates. This storage requirement limits practical applications to about 650 basis functions of 100 correlated electrons on a computer with 120 GB of memory.

IV Computational Details

We tested the DSRG-MRPT3 method on the ground-state potential energy curves (PECs) of \ceF2, \ceH2O2, \ceC2H6, and \ceN2 by breaking the respective F–F, O–O, C–C, and N–N bonds. The DSRG-MRPT3 PECs were compared to those of other multireference methods including DSRG-MRPT2,Li and Evangelista (2015) MR-LDSRG(2),Li and Evangelista (2016) NEVPT2,Angeli et al. (2001); Angeli, Cimiraglia, and Malrieu (2002); Angeli, Pastore, and Cimiraglia (2007) CASPT2,Andersson, Malmqvist, and Roos (1992); Werner (1996) CASPT3,Werner (1996) MRCISD,Knowles and Werner (1988); Werner and Knowles (1988) and MRCISD with Davidson correction (MRCISD+Q).Langhoff and Davidson (1974) All multireference computations are based on minimal CASSCF references, specifically, CAS(2,2) for \ceF2, \ceH2O2, and \ceC2H6, while CAS(6,6) for \ceN2. Full configuration interaction (FCI) data served as the benchmark for \ceF2 and \ceN2.Bytautas et al. (2007); Li and Evangelista (2016) For \ceH2O2 and \ceC2H6, we took the reference data from Ref. 110, which were computed using coupled cluster singles, doubles and triples augmented with second-order perturbative quadruples corrections [CCSDT(2)${}_{\text{Q}}$].Hirata et al. (2004) Dunning’s correlation-consistent double-$\zeta$ (cc-pVDZ) basis setDunning (1989) was used and the molecular orbitals constructed mainly from the 1$s$ orbitals of C, N, O, and F atoms were frozen in all post-CASSCF computations.

To measure the quality of the potential energy curves we use the nonparallelism error (NPE) computed with respect to a reference method [FCI or CCSDT(2)Q] over a range of bond lengths $\bf R$, defined as:

[TABLE]

where $\Delta E(r)$ is the error with respect to the reference energy at bond length $r$.

As an application to medium-sized molecules, we studied the singlet-triplet splitting ($\Delta_{\rm ST}=E_{\rm T}-E_{\rm S}$) of 9,10-anthracyne (9,10-didehydroanthracene, see Fig. 1), which is recently found capable of retro-Bergman cyclization on a NaCl/Cu(111) surface when manipulated with the CO tip of an atomic force microscope.Schuler et al. (2016) We first optimized the geometries of singlet and triplet 9,10-anthracyne at the CASSCF(2,2)-DSRG-MRPT3/cc-pVDZ level of theory using gradients from 3-point finite-difference computations. These geometries were characterized as minima by finite-difference harmonic vibrational analyses. Then we computed the $\Delta_{\rm ST}$ of 9,10-anthracyne using various multireference methods including DSRG-MRPT2/3, NEVPT2, the partially contracted version of CASPT2Celani and Werner (2000) and MRCISDShamasundar, Knizia, and Werner (2011) as implemented in the RS2C and CIC modules of MOLPRO,Werner et al. (2012) Mukherjee multireference coupled cluster theory with singles and doubles (Mk-MRCCSD),Mahapatra et al. (1998); Mahapatra, Datta, and Mukherjee (1999b); Evangelista, Allen, and Schaefer (2006); Mahapatra, Datta, and Mukherjee (2010) and Mk-MRCCSD with perturbative triples [Mk-MRCCSD(T)]Evangelista et al. (2010) based on a CASSCF reference. Both cc-pVDZ and cc-pVTZ basis sets were adopted and the $1s$-like orbitals on carbon atoms were excluded for dynamic correlations. The two-electron integrals in DSRG-MRPT2/3 computations were approximated by Cholesky decompositionBeebe and Linderberg (1977); Koch, Sánchez de Merás, and Pedersen (2003); Aquilante et al. (2008, 2011) with a threshold of $10^{-8}$ a.u. for geometry optimizations and $10^{-6}$ a.u. for single points.

The NEVPT2, CASPT2, CASPT3, and MRCISD energies were computed using the Molpro 2015.1 programWerner et al. (2012, 2015) and the remaining were obtained from Psi4.Turney et al. (2012) The $xyz$ coordinates of the optimized 9,10-anthracyne, as well as all energies on the PECs of \ceF2, \ceH2O2, \ceC2H6, and \ceN2 are available in the supplementary material.SI For convenience, the supplementary materialSI also includes the internal coordinates of \ceH2O2 and \ceC2H6, which are taken from Refs. 110 and 118, respectively.

V Results

V.1 Potential energy curves

In this section we assess the accuracy of DSRG-MRPT3 by investigating four bond breaking processes. We use a recent benchmark set by Yang, Jalan, Green, and Truhlar (YJGT) that was used to compare the performance of numerous single-reference coupled cluster and multireference methods.Yang et al. (2013) The YJGT set contains \ceF2, \ceH2O2 and \ceC2H6 scanned along the \ceF-F, \ceO-O, and \ceC-C bonds, respectively. To test a multiple-bond breaking process, we consider the dissociation curve of \ceN2, a routine benchmark for multireference methods.Yanai and Chan (2007); Das, Mukherjee, and Kállay (2010); Kurashige et al. (2014); Fosso-Tande et al. (2016)

The errors in the computed ground-state potential energy curves of \ceF2, \ceH2O2, \ceC2H6 and \ceN2 are shown in Fig. 2 and Table 1. Fig. 2 shows that third-order MRPTs recover a larger fraction of electron correlation than the corresponding second-order methods. For example, going from DSRG-MRPT2 to DSRG-MRPT3, the maximum error is reduced by an amount between one third (\ceF2) and one half (\ceH2O2). Despite the fact that the DSRG-MRPT3 does not always yield the curve with the smallest absolute error, we find that the average NPE of the DSRG-MRPT3 (3.87 m$E_{\rm h}$) is smaller than that of CASPT3 (5.00 m$E_{\rm h}$), MRCISD (4.22 m$E_{\rm h}$), and the MR-LDSRG(2) (3.98 m$E_{\rm h}$). Interestingly, the DSRG-MRPT3 significantly improves upon the PEC of \ceN2 computed with DSRG-MRPT2. For this molecule, the DSRG-MRPT3 gives a NPE equal to 4.78 m$E_{\rm h}$ vs. 18.34 m$E_{\rm h}$ for the DSRG-MRPT2. Overall, a comparison of the various MR-DSRG methods considered here suggests that the accuracy of these methods follows the trend: DSRG-MRPT2 $\ll$ DSRG-MRPT3 $<$ MR-LDSRG(2).

Figure 3 shows a comparison of different reference-relaxation (RR) approaches applied to perturbative and nonperturbative MR-DSRG methods. For all combinations of molecules and methods the one-step RR curves are more accurate than the corresponding unrelaxed curves. For example, the average NPE of the unrelaxed DSRG-MRPT3 is 1.09 m$E_{\rm h}$ higher than the one-step RR approach. For the MR-LDSRG(2) series, the one-step RR results (labeled with the subscript “r1”) show energy errors that are halfway between those from unrelaxed and fully relaxed computations. These findings suggest that at the perturbative level it is generally advantageous to introduce reference relaxation effects and that a one-step approach is a good compromise between accuracy and cost. Beyond perturbation theory we observe that full relaxation of the reference does in certain cases give results inferior to the one-step correction (see the \ceF2 and \ceH2O2 curves). However, in our opinion the fully relaxed MR-LDSRG(2) is the most rigorous formulation of the nonpertubative MR-DSRG approach and should be the default choice in applications to chemical problems.

V.2 Singlet-triplet splittings of 9,10-anthracyne

The simplest active space for 9,10-anthracyne is a CAS(2,2) that consists of two electrons in two $\sigma$ orbitals that belong to the dehydrogenated carbon atoms. From our previous experience with p-benzyne,Li and Evangelista (2015, 2016) adding six $\pi$-type orbitals to the active space increases the singlet-triplet splitting ($\Delta_{\rm ST}=E_{\rm T}-E_{\rm S}$) by up to 2 kcal mol*-1*. As noted before,Li and Evangelista (2016) this energy shift likely results from an improved treatment of static correlation effects by the larger active space, since both the CASSCF and DSRG-MRPT2 singlet-triplet splittings are shifted by a similar amount. For 9,10-anthracyne, there are fourteen $\pi$ orbitals, so that a full treatment of the $\sigma$ and $\pi$ orbitals leads to a CAS(16,16) reference wave function. The sixteen orbitals of the CAS(16,16) space are shown in Figure 4 along with their corresponding occupation numbers (ONs). The ON for each of the sixteen orbitals lies in the range of 0.02–1.98, justifying the inclusion of all $\pi$ orbitals to the active space.Bofill and Pulay (1989); Keller et al. (2015) Since—with the exception of the DSRG-MRPT2/3 approaches—the CAS(16,16) is too large for the methods considered here, we also report results for three smaller active spaces that include a subset of the $\pi$ space: CAS(4,4), CAS(8,8), and CAS(12,12). Specifically, the CAS(4,4) active space includes the $2b_{1g}$ and $3b_{3u}$ orbitals, while the CAS(8,8) further adds $1a_{u},2b_{2g},2a_{u}$, and $3b_{2g}$ orbitals. The CAS(12,12) active space is obtained by augmenting the CAS(8,8) reference with the $2b_{3u},1b_{1g},4b_{3u}$, and $3b_{1g}$ orbitals.

The CASSCF(2,2)-DSRG-MRPT3/cc-pVDZ optimized geometries of singlet and triplet 9,10-anthracyne are shown in Fig. 1. All \ceC-C bond lengths are consistent with those predicted by density functional theory (DFT) except for the \ceC_4a-C_9a bond of the singlet ground state.Schuler et al. (2016) Interestingly, the singlet \ceC_4a-C_9a bond length reported in Ref. 114 (1.578 Å) is not only 0.11 Å longer than our prediction (1.466 Å), but also larger than a typical $sp^{3}$ hybridized \ceC-C single bond like in ethane (1.536 Å, from Ref. 124).

In Table 2, we list the singlet-triplet splittings ($\Delta_{\rm ST}=E_{\rm T}-E_{\rm S}$) of 9,10-anthracyne computed using numerous multireference methods and the cc-pVTZ basis set. All values are adjusted by $+0.264$ kcal mol*-1* to account for zero-point harmonic vibrational energy corrections computed at the CASSCF(2,2)-DSRG-MRPT3/cc-pVDZ level of theory. Because the singlet-triplet splitting computed with the cc-pVDZ and cc-pVTZ basis are in very good agreement (within 0.5 kcal mol*-1*), we report results only for the latter basis set and provide data for the former in the supplementary material.SI Benchmark computations on p-benzyne show that Mk-MRCCSD(T) yields an accurate singlet-triplet splitting,Evangelista et al. (2010, 2012) therefore, in the absence of an experimental value for 9,10-anthracyne, we take the $\Delta_{\rm ST}$ from Mk-MRCCSD(T) as our reference.

In contrast to the DFT prediction of a triplet ground state,Schuler et al. (2016) all methods reported in Table 2 favor a singlet ground state. For the minimal active space, the CASPT2 $\Delta_{\rm ST}$ (4.1 kcal mol*-1*) is in excellent agreement with our reference value [4.0 kcal mol*-1* from Mk-MRCCSD(T)], while the remaining methods underestimate the $\Delta_{\rm ST}$ by 1–2 kcal mol*-1*. As the active-space size increases, better agreements to the Mk-MRCCSD(T) value are observed for all methods. For example, when using a CASSCF(16,16) reference, singlet-triplet splittings of DSRG-MRPT2 and -MRPT3 (s = 0.5 $E_{\rm h}^{-2}$) are 3.8 and 3.6 kcal mol*-1*, respectively. For both DSRG-MRPT2 and -MRPT3 approaches, reference relaxation shifts the $\Delta_{\rm ST}$ to higher values, improving the agreement with the Mk-MRCCSD(T) result. The effect of reference relaxation is largest for the minimal active space, where in the case of the DSRG-MRPT2 ($s=0.5$ $E_{\rm h}^{-2}$) it amounts to an increase in the singlet-triplet splitting by 0.4 kcal mol*-1*. We also observe that the $\Delta_{\rm ST}$ obtained from DSRG-MRPT2 and -MRPT3 approaches does not depend significantly on the value of the flow variable: as $s$ increases from 0.5 to 1.0 $E_{\rm h}^{-2}$, the largest variation in $\Delta_{\rm ST}$ is less than 0.4 kcal mol*-1*.

We note that results for the CAS(8,8) reference are somewhat erratic and they deviate significantly from those obtained with CAS(4,4), CAS(12,12), and CAS(16,16) references. A plausible reason for this inconsistency is that some of the singlet and triplet semicanonical orbitals of CASSCF(8,8) are considerably different from those obtained with other active spaces (see supplementary materialSI ).

It is instructive to compare the $\Delta_{\rm ST}$ between 9,10-anthracyne and p-benzyne. The $\Delta_{\rm ST}$ of p-benzyne is measured to be $3.8\pm 0.4$ kcal mol*-1* by photoelectron spectroscopy,Wenthold, Squires, and Lineberger (1998) and Mk-MRCCSD(T)/cc-pVTZ with a CASSCF(2,2) predicts $\Delta_{\rm ST}$ = 4.45 kcal mol*-1*.Evangelista et al. (2012) Therefore, we would expect the experimental $\Delta_{\rm ST}$ of 9,10-anthracyne to fall in the range 2.9–3.7 kcal mol*-1*. The predictions of NEVPT2 and DSRG-MRPT2 always fall in this range for all the active spaces considered here, while CASPT2 constantly overestimate the $\Delta_{\rm ST}$ of 9,10-anthracyne. The DSRG-MRPT3 results are always in line with those of MRCISD, and the minimal active space seems to be insufficient in this case.

VI Conclusions

We have introduced the DSRG-MRPT3 scheme, a renormalized third-order multireference perturbation theory derived from the MR-DSRG approach.Li and Evangelista (2015, 2016); Hannon, Li, and Evangelista (2016) By working in a semicanonical basis and using a zeroth-order Hamiltonian [$\hat{H}^{(0)}$] that contains only the diagonal blocks of the Fock matrix, the DSRG-MRPT3 energy may be obtained by a non-iterative procedure that scales as ${\cal O}(N_{\bf C}^{2}N_{\bf V}^{4})$. The cost of evaluating the DSRG-MRPT3 energy is reduced by truncating each commutator in the effective Hamiltonian with only one- and two-body operators.Yanai and Chan (2006, 2007) As a result, the DSRG-MRPT3 energy and amplitude equations require up to three-body density cumulants. Furthermore, we show that the third-order effective Hamiltonian may be easily formed from the first- and second-order amplitudes, and that it may be diagonalized to obtain a relaxed model space.

We have benchmarked the DSRG-MRPT3 method on the ground-state potential energy curves of \ceF2, \ceH_2O_2, \ceC_2H_6 and \ceN2 and found that, on average, the DSRG-MRPT3 nonparallelism error (3.9 m$E_{\rm h}$) is comparable to that of CASPT3 (5.0 m$E_{\rm h}$) and MRCISD (4.2 m$E_{\rm h}$). We have also shown that accounting for reference relaxation effects with a one step diagonalization of the effective Hamiltonian does generally improve the quality of DSRG-MRPT2 and -MRPT3 results. For example, in the case of DSRG-MRPT3, the one-step relaxation approach reduces the average NPE by 1.1 m$E_{\rm h}$. Since the cost of evaluating the one-step relaxed energy is a small fraction of the cost of an unrelaxed DSRG-MRPT3 computation, we recommend the former approach as the default method.

In our opinion, the DSRG-MRPT3 is best viewed as an economical multireference method that offers a good compromise between accuracy and computational cost. Compared to more expensive nonperturbative (i.e. coupled cluster like) multireference methods, the DSRG-MRPT3 has the advantage that it does not rely on an iterative procedure and has reduced I/O costs. In its current implementation, the DSRG-MRPT3 may be routinely applied in computations with a hundred correlated electrons and 500–600 basis functions. This point is illustrated with computations of the singlet-triplet splitting ($\Delta_{\rm ST}$) of 9,10-anthracyne (C14H8). The DSRG-MRPT3 based on a CASSCF(16,16) reference with $s=1.0$ $E_{\rm h}^{-2}$ predicts the $\Delta_{\rm ST}$ to be 3.9 kcal mol*-1*, only 0.1 kcal mol*-1* smaller than our best estimate from Mk-MRCCSD(T). In our experience, the CAS(16,16) active space is too large for practical CASPT2 and NEVPT2 computations on a single computer node. With Cholesky decomposed integrals (using a $10^{-6}$ a.u. threshold), the CAS(16,16) DSRG-MRPT3 energy of 9,10-anthracyne may be computed in about 12.5 hours using 16 threads on two Intel Xeon E5-2650 v2 processors and 128GB of memory.

In conclusion, we have shown that the DSRG framework may be used to formulate a third-order multireference perturbation theory that avoids some of the major limitations of other MRPT approaches and has a favorable accuracy/cost ratio. These results suggest that it might be worthwhile to explore multi-state generalizations of the DSRG-MRPT2 and DSRG-MRPT3 methods to compute excited state energies and their analytic gradients. This work also presents additional benchmark results that validate the accuracy of perturbative and nonperturbative computational methods based on the MR-DSRG formalism.

Acknowledgements.

C.L. would like to thank Dr. Zhi Sun for performing the CASSCF(12,12)-NEVPT2 computations of 9,10-anthracyne. This work was supported by Emory University and the U.S. Department of Energy under Award No. DE-SC0016004.

Appendix A Reference relaxation

In this section we provide details of the reference relaxation procedure used in the DSRG-MRPT2/3 and MR-LDSRG(2) methods. We first introduce the many-body expression of the DSRG transformed Hamiltonian

[TABLE]

where $\bar{H}_{0}(s)$ is the scalar term obtained by summing the reference energy ($E_{0}$) and all the fully contracted contributions from $\bar{H}(s)$. The quantity $\bar{H}_{k}(s)$ contains the $k$-body contributions to $\bar{H}(s)$:

[TABLE]

The approximated Hamiltonian $\bar{H}_{1,2}(s)$ is thus obtained when setting $n=2$ in Eq. (34).

In order to solve the eigenvalue problem of Eq. (26), we express the transformed Hamiltonian using operators normal ordered with respect to the true vacuum. Defining the two-particle density matrix $\gamma_{rs}^{pq}=\braket{\Psi_{0}}{\hat{a}^{\dagger}_{p}\hat{a}^{\dagger}_{q}\hat{a}_{s}\hat{a}_{r}}{\Psi_{0}}$, the $\bar{H}_{1,2}(s)$ operator is written as [“$(s)$” is dropped for brevity],Li and Evangelista (2016)

[TABLE]

The quantities $[{o}_{p}^{q}={\bar{H}}_{p}^{q}-\sum_{ij}{\bar{H}}_{pi}^{qj}\gamma_{j}^{i}]$ and $[{w}_{pq}^{rs}={\bar{H}}_{pq}^{rs}]$ may be considered as MR-DSRG dressed one- and two-electron integrals, respectively.

A common strategy to compute the CASCI energy is to fold the contribution of core orbitals into a scalar term ($c_{0}$) and the one-body term labeled by all active indices (${c}_{u}^{v}$). Using the dressed integrals (${o}_{p}^{q},{w}_{pq}^{rs}$), the additional scalar term coming from the core orbitals is given by

[TABLE]

and the corresponding modified one-body operator is

[TABLE]

Obviously, Eq. (A) cancels the second line of Eq. (A). Thus we prove that only those elements of $\bar{H}(s)$ labeled by active indices are necessary to obtain the relaxed MR-DSRG energy. The cost to evaluate ${\bar{H}}_{uv}^{xy}$ scales as ${\cal O}(N_{\bf V}^{2}N_{\bf A}^{4})$ and is significantly smaller than then the cost to evaluate the DSRG-MRPT3 energy, which scales as ${\cal O}(N_{\bf V}^{4}N_{\bf C}^{2})$.

Appendix B Other choices of zeroth-order Hamiltonian

The DSRG-MRPT3 approach proposed in this work uses a diagonal one-body zeroth-order Hamiltonian. However, partitioning schemes that include two-body operators may also be used with the DSRG-MRPT3, including the Dyall HamiltonianDyall (1995) and the retaining-excitation (RE) Hamiltonian.Fink (2006, 2009) In this section we analyze the implications of including two-body operators in the zeroth-order Hamiltonian, focusing on the computational efficiency and the avoidance of intruders in the DSRG-MRPT3. If these partitionings led to an intruder-free DSRG-MRPT approach, then one could avoid the $s$-dependence of the energy by taking the limit $s\rightarrow\infty$. Unfortunately, our analysis shows that even when two-body interactions are included in the zeroth-order Hamiltonian, some form of renormalization is necessary to avoid intruders.

When written in normal-ordered form with respect to the reference wave function, the zeroth-order Hamiltonian of Dyall may be written as:

[TABLE]

where $\hat{F}^{(0)}$ is defined as in Eq. (II.2). A first important consequence of the presence of the two-body term is the fact that already at first-order the commutator $[\hat{H}^{(0)},\hat{A}^{(1)}(s)]$ produces three-body terms. Consequently, evaluating the second-order energy requires computing of a subset of the three-body operators. This problem may be avoided by invoking the “1,2” operator approximation, which introduces truncation errors already in the second-order energy.

Another complication resulting from the two-body terms in $\hat{H}^{(0)}$ is the need for an iterative solution of the amplitudes. Again, we take the zeroth-order Hamiltonian to be $\hat{H}^{(0)}_{\rm Dyall}$ and consider the limit $s\rightarrow\infty$. The first-order amplitudes corresponding to the promotion two electrons from active to virtual orbitals [${t}_{ef}^{uv,(1)}(s)$] are determined by the equation:

[TABLE]

where for clarity we dropped the symbol “(s)” from the amplitudes. In Eq. (B) the amplitude ${t}_{ef}^{uv,(1)}$ is coupled to all amplitudes of the form ${t}_{ef}^{xy,(1)}$ where $x,y\in\mathbf{A}$. Consequently, the DSRG-MRPT based on the Dyall Hamiltonian requires an iterative solution of a set of linear equations for certain classes of amplitudes. Nevertheless, the most numerous amplitudes, double excitations from core orbitals to virtual orbitals, ${t}_{ef}^{mn,(1)}$, may still be computed via a noniterative procedure.

Last, we consider whether the use of Dyall’s Hamiltonian in the DSRG-MRPT formalism may avoid the intruder state problem, as it is know to be the case in NEVPT.Angeli et al. (2001, 2006); Angeli, Pastore, and Cimiraglia (2007) To this end, we evaluate the diagonal preconditioner or “shifted” denominator corresponding to certain classes of excitations prone to give small denominators and determine if extra terms that arise from Dyall’s Hamiltonian may help avoid divergences. We first go back to Eq. (B), and identify the shifted denominator (${D}_{ef}^{uv}$) for the ${t}_{ef}^{uv,(1)}$ amplitudes with the expression:

[TABLE]

The analysis of this expression is simpler in the basis of natural orbitals. Defining the natural occupation of orbital $\phi^{u}$ as $n_{u}=\gamma^{u}_{u}$ and assuming a natural orbital basis, we may simplify Eq. (41) to:

[TABLE]

The term $\epsilon_{u}+\epsilon_{v}-\epsilon_{e}-\epsilon_{f}$ is the standard Møller–Plesset denominator, and in general it is negative since virtual orbitals are assumed to lie higher in energy than active orbitals. The term $\frac{1}{2}v^{uv,(0)}_{uv}(1-n_{u}-n_{v})$ arises from the two-body active part of the Hamiltonian and may be positive (when $n_{u}+n_{v}<1$), negative (when $n_{u}+n_{v}>1$), or zero. Consequently, the additional terms that arise from the active part of the two-body operator cannot guarantee that ${D}_{ef}^{uv}\neq 0$. Excitations most prone to small denominators are singles coupled with a spectator excitation (from active to active orbitals), for example, ${t}_{ex}^{uv,(1)}(s)$. In this case we also find that the shifted denominators may accidentally be zero since the contribution of the two-body Hamiltonian does not have a well defined sign.

Next, we consider single excitations. Interestingly, in this case we find that the Dyall Hamiltonian does not lead to modified energy denominators. Hence, intruder states may still arise when singles denominators that involve one active orbital approach zero. This is in contrast with NEVPT2, where single denominators are shifted by two-body integrals. This difference arises from the use of many-body conditions in DSRG-MRPT and projective conditions in the case of NEVPT2.

In summary, our analysis shows that the use of Dyall’s partitioning of the Hamiltonian cannot avoid the intruder state problem in perturbation theories like the DSRG-MRPT, which are derived from a set of many-body conditions. In other words, renormalization of small denominators is necessary even when the zeroth-order Hamiltonian contains two-electron terms. From another perspective, the difficulty in converging the nonperturbative MR-LDSRG(2) equationsLi and Evangelista (2016) when $s\rightarrow\infty$, also suggests that using other zeroth-order Hamiltonians with two-body contributions is unlikely solution to the intruder state problem in DSRG multireference perturbation theory.