Improving the efficiency of the multireference driven similarity renormalization group via sequential transformation, density fitting, and the non-interacting virtual orbital approximation
Tianyuan Zhang Chenyang Li Francesco A. Evangelista



[
\begin@twocolumnfalse






This study examines several techniques to improve the efficiency of the linearized multireference driven similarity renormalization group truncated to one- and two-body operators MR-LDSRG(2).
We propose a sequential MR-LDSRG(2) sq-MR-LDSRG(2) approach, in which one-body rotations are folded exactly into the Hamiltonian.
This new approach is combined with density fitting (DF) to reduce the storage cost of two-electron integrals. 
To further avoid the storage of large four-index intermediates, we propose a non-interacting virtual orbit (NIVO) approximation in which tensor blocks labeled by three and four virtual indices are neglected.
The NIVO approximation reduces the computational cost prefactor of the MR-LDSRG(2) bringing it closer to that of coupled cluster with singles and doubles (CCSD).
We test the effect of the DF and NIVO approximations on the MR-LDSRG(2) and sq-MR-LDSRG(2) methods by computing properties of eight diatomic molecules. The diatomic constants obtained by DF-sq-MR-LDSRG(2)+NIVO are found to be as accurate as those from the original MR-LDSRG(2) and coupled cluster theory with singles, doubles, and perturbative triples.
Finally, we demonstrate that the DF-sq-MR-LDSRG(2)+NIVO scheme can be applied to chemical systems with more than 550 basis functions by computing the automerization energy of cyclobutadiene with a quintuple-ζ basis set.
The predicted automerization energy is found similar to the values computed with Mukherjee's state-specific multireference coupled cluster theory with singles and doubles.



\end@twocolumnfalse
]





























§ INTRODUCTION


The failure of conventional many-body methods to describe near-degenerate electronic states has motivated the development of many efficient and practical multireference approaches, including perturbation theories (MRPTs)<cit.> and multireference configuration interaction (MRCI) schemes.<cit.>
Considerable efforts have been dedicate to the development of multireference coupled cluster (MRCC) methods,<cit.>
with the goal of creating nonperturbative theories that are both size extensive and systematically improvable.
Analogous many-body methods based on unitary transformations have received considerably less attention.<cit.>
Unitary theories have, in principle, two advantages over traditional coupled cluster approaches: 1) the energy satisfies the variational principle, and 2) the transformed Hamiltonian is Hermitian.
The latter property is an important advantage in new applications of unitary methods to quantum computing, both in quantum algorithms<cit.> and downfolding approaches aimed at reducing the number of orbitals in quantum computations.<cit.>

One of the main obstacles in the formulation of both single- and multi-reference unitary coupled cluster theories is that they lead to nonterminating equations.
The central quantity evaluated in these approaches is the similarity transformed Hamiltonian  (H̅) defined as

    Ĥ→H̅ = Û^†ĤÛ= e^-ÂĤ e^Â

where (Ĥ) is the bare Hamiltonian and Û is a unitary operator. In writing this transformation, we have expressed Û as the exponential of the anti-Hermitian operator Â (Â^† = - Â), which is commonly written in terms of the coupled cluster excitation operator T̂ as Â = T̂ - T̂^†.
Using the Baker–Campbell–Hausdorff (BCH) identity,<cit.> the transformed Hamiltonian may be computed as the following commutator series

    H̅ = Ĥ + [ Ĥ, Â] + 1/2![ [ Ĥ, Â], Â] + 1/3![ [ [ Ĥ, Â], Â], Â]….

Since the operator Â contains both excitations and de-excitations, contractions are possible among components of Â, and as a consequence, the BCH series given in Eq. (<ref>) is nonterminating.

Various approximations have been proposed to evaluate the unitarily transformed Hamiltonian.
Perhaps the simplest way to approximate the nonterminating unitary series is to truncate the BCH expansion after a certain number of commutators.<cit.>
Proof-of-principle studies on unitary coupled cluster (CC) theory<cit.> suggest that for a a series containing up to n-nested commutators, the error decays as 10^-n, and about four commutators are necessary to achieve sub-milliHartree accuracy.
Taube and Bartlett<cit.> have suggested tractable approximations to unitary CC theory based on the Zassenhaus expansion that are exact for a given number of electrons.
A common way to truncate the unitary BCH series is to use a recursive approximation of the commutator [ · , Â], as suggested by Yanai and Chan.<cit.>
In their linear truncation scheme, these authors proposed to approximate each single commutator [ · , Â] with its scalar and one- and two-body components, which we indicate as [ · , Â]_0,1,2.
Since in this truncation scheme the commutator [ · , Â]_0,1,2 preserves the many-body rank (number of creation and annihilation operators) of the Hamiltonian, the full BCH series can then be evaluated via a recursive relation.
An advantage of this approach is that closed-form expressions for terms like [Ô, Â]_0,1,2, where Ô is an operator containing up to two-body terms  can be easily derived.
This truncation scheme has been employed in canonical transformation (CT) theory<cit.> and has been used to truncate normal-ordered  equations in the flow renormalization group of Wegner.<cit.>

We have recently developed a multireference driven similarity renormalization group (MR-DSRG)<cit.> approach that avoids the multiple-parentage problem<cit.> and numerical instabilities<cit.> encountered in other nonperturbative multireference methods.
In the MR-DSRG, we perform a unitary transformation of the Hamiltonian controlled by a flow parameter, which determines to which extend the resulting effective Hamiltonian is  band diagonal.<cit.>
The simplest nonperturbative approximation, the linearized MR-DSRG truncated to two-body operators [MR-LDSRG(2)],<cit.> assumes that Â contains up to two-body operators (singles and doubles) and employs the linear commutator approximation of Yanai and Chan.
Preliminary benchmarks indicate that the MR-LDSRG(2) method is more accurate than CCSD around equilibrium geometries, and that this accuracy is preserved along  potential energy curves, especially for single-bond breaking processes.<cit.>
The cost to evaluate a single commutator in the MR-LDSRG(2) scales as O (N_ C^2  N_ V^2 N^2) =  O (N_ C^2 N_ V^4 + N_ C^3 N_ V^3 + …) where N_ C, N_ V, and N are the numbers of core, virtual, and total orbitals, respectively.
This scaling is identical to that of CC with singles and doubles (CCSD) [O (N_ C^2 N_ V^4)]. However, while in CCSD the most expensive term is evaluated only once, in the MR-LDSRG(2) the same term must be evaluated for each nested commutator.
Consequently, if the BCH series is truncated after n + 1 terms, the computational cost of the MR-LDSRG(2) is roughly n times that of CCSD, where approximately ten or more terms are usually required to convergence the energy error in the BCH series to 10^-12 .<cit.>
A second reason is that computing the BCH series requires storing large intermediate tensors with memory costs that scale as 𝒪(N^4).
When these intermediates are stored in memory, practical computations are limited to 200–300 basis functions on a single computer node.

In this work, we combine a series of improvements and approximations to reduce the computational and memory requirements of the MR-LDSRG(2) down to a small multiple of the cost of CCSD.
To begin, we consider an alternative ansatz for the MR-DSRG based on a sequential similarity transformation.<cit.>
This sequential ansatz reduces the complexity of the MR-DSRG equations and, when combined with integral factorization techniques, reduces significantly the cost to evaluate singles contributions.
Second, we apply density fitting (DF)<cit.> to reduce the memory requirements and the I/O cost by avoiding the storage of of two-electron repulsion integrals.
Together with Cholesky decomposition (CD)<cit.> and other tensor decomposition schemes,<cit.> these techniques have been crucial in enabling computations with 1000 or more basis functions and found application in numerous electronic structure methods,<cit.>
including coupled cluster methods.<cit.>
Third, we reduce the cost of MR-LDSRG(2) computations by neglecting operators that involve three or more virtual electrons.
We term this truncation scheme the non-interacting virtual orbital (NIVO) approximation.
A perturbative analysis of the NIVO approximation shows that the errors introduced appear at third order.
To the best of our knowledge, such approximation has not been introduced in multireference theories, but it is analogous to other truncation schemes used in CCSD in which certain diagrams are have modified coefficients<cit.> or are completely removed.<cit.>

This paper is organized in the following way.
In Sec. <ref> we present an overview of the MR-DSRG theory, discuss the sequential MR-DSRG, and introduced the NIVO approximation.
Details of the implementation together with a discussion of timings are given in Sec. <ref>.
In Sec. <ref> we assess the accuracy of several MR-LDSRG(2) schemes on a benchmark set of several diatomic molecules and determine the automerization barrier of cyclobutadiene.
Finally, in Sec. <ref> we conclude this work with a discussion of the main results.



§ THEORY



We first define the orbital labeling convention employed in this work.
The set of molecular spin orbitals G≡{ϕ^p, p = 1, 2, …, N} is partitioned into core (C), active (A), and virtual (V) components of sizes N_ C, N_ A, and N_ V, respectively.
We use indices m,n to label core orbitals, u,v,x,y to label active orbitals, and e,f,g,h to label virtual orbitals.

For convenience, we also define the set of hole (H = C ∪ A) and particle (P = A ∪ V) orbitals with dimensions N_ H = N_ C + N_ A and N_ P = N_ A + N_ V, respectively.
Hole orbitals are denoted by indices i,j,k,l and those of particle by indices a,b,c,d.
General orbitals are labeled by indices p,q,r,s.








The MR-DSRG assumes a multideterminantal reference wave function ():

    |⟩ = ∑_μ = 1^d c_μ|Φ^μ⟩.

In this work, the set of determinants {Φ^μ; μ = 1, 2, …, d} in Eq. (<ref>) is assumed to form a complete active space (CAS) and the corresponding coefficients {c_μ; μ = 1, 2, …, d} are obtained from a CAS self-consistent field (CASSCF) computation.<cit.> 

All operators are then normal ordered with respect to  according to the scheme of Mukherjee and Kutzelnigg.<cit.>
For example, the bare Hamiltonian is expressed as

    Ĥ = E_0 + ∑_pq f_p^q pq + 1/4∑_pqrs v_pq^rspqrs ,

where E_0=⟨|Ĥ||$⟩ is the reference energy andpq…rs… =  p q …s r is a product of creation (p ≡p^†) and annihilation (p) operators in its normal-ordered form, as indicated by the curly braces.
The generalized Fock matrix (f_p^q) introduced in Eq. (<ref>) is defined as


    f_p^q = h_p^q + ∑_rsvprqsrs,

wherehpq = ⟨ϕ_p | ĥ | ϕ^q|$⟩ and vpqrs = ϕ_p ϕ_qϕ^r ϕ^s are the one-electron and anti-symmetrized two-electron integrals, respectively.
Here, we have also used the one-particle reduced density matrix (1-RDM) defined as pq = ⟨| pq ||$⟩.



 §.§ Review of the MR-DSRG method

  
The MR-DSRG performs a parametric unitary transformation of the bare Hamiltonian analogous to Eq. (<ref>), whereby the anti-Hermitian operatorÂ (s)depends on the so-called flow parameter, a real numbersdefined in the range of[0,∞). The resulting transformed Hamiltonian [H̅(s)] is a function ofsdefined as

    H̅ (s) = e^-Â (s)Ĥ e^Â (s).

The operatorÂ(s)is a sum of many-body operators with rank ranging from one up to the total number of electrons (n),

    Â(s) = ∑_k=1^nÂ_k(s),

whereÂ_k(s)is thek-body component ofÂ(s).
In the DSRG the operatorsÂ_k (s)are written asÂ_k(s) = T̂_k (s) - T̂^†_k (s), whereT̂_k (s)is ans-dependent cluster operator defined as

    T̂_k(s) = 1/(k!)^2∑_ij⋯^𝐇∑_ab⋯^𝐏tab⋯ij⋯ (s) ab⋯ij⋯.

Note that the cluster amplitudestab⋯ij⋯ (s)exclude internal excitations, which are labeled only with active orbital indices.



The DSRG transformed Hamiltonian has a many-body expansion similar to Eq. (<ref>),

    H̅(s)=E̅_0(s)+∑_pqH̅_q^p(s){â_p^q}+1/4∑_pqrsH̅_rs^pq(s) {â_pq^rs}+⋯,

where,

    E̅_0(s) = ⟨| H̅(s) ||,⟩

is the DSRG energy and the tensorsH̅_rs…^pq…(s)are analogous to one- and two-electron integrals but dressed with dynamical correlation effects.

The goal of the DSRG transformation is to decouple the interactions between the reference wave function () and its excited configurations.
Such interactions are the couplings between hole and particle orbitals represented by generalized excitation [H̅ab…ij… (s) ab…ij…] and de-excitation [H̅ij…ab… (s) ij…ab…] operators, whereij⋯∈𝐇andab⋯∈𝐏, excluding cases where all the indices are active orbitals.
These terms ofH̅(s)that the DSRG transformation aims to suppress are called the off-diagonal components and will be denoted asH̅_od(s).
Instead of achieving a full decoupling of the off-diagonal components [i.e.,H̅_od(s) = 0], we demand that the DSRG transformation achieves a partial decoupling, avoiding the components ofH̅_od(s)with small or vanishing Møller–Plesset energy denominators. 
This condition is imposed via the DSRG flow equation, a nonlinear implicit equation of the form

    H̅_ od(s) = R̂ (s),

where the source operatorR̂(s)is Hermitian and continuous ins.
The role of the source operator in the DSRG flow equation is to drive the off-diagonal elements ofH̅(s)to zero, a goal achieved by an appropriate parameterization ofR̂(s).<cit.>

In the linearized MR-DSRG approximation [MR-LDSRG(2)], the BCH series [Eq. (<ref>)] is evaluated by keeping up to two-body normal-ordered operators of each commutator. The transformed Hamiltonian can then be evaluated by the following recursive equations

    {Ĉ^(k+1)(s) =    1/k + 1[Ĉ^(k)_1,2 (s), Â(s)]_0,1,2,
    H̅^(k+1)(s) =    H̅^(k)(s) +  Ĉ^(k+1) (s),
    .

starting fromĈ^(0)(s) = H̅^(0) = Ĥand iterating until the norm ofĈ^(k+1)_1,2 (s)is less than a given convergence threshold.












The solution of Eq. (<ref>) yields a set of amplitudestab…ij… (s)that define the operatorÂ(s)and the DSRG transformed HamiltonianH̅(s).
From this latter quantity, the MR-DSRG electronic energy is computed as the expectation value with respect to the reference

    E̅_0(s) = ⟨| H̅(s) ||.⟩

We refer the energy computed using Eq. (<ref>) as the unrelaxed energy since the reference coefficients are not optimized.
To include reference relaxation effects, we require thatis an eigenstate ofH̅(s)within the space of reference determinants, a condition that is equivalent to solving the eigenvalue problem

    ∑_μ=1^d⟨Φ_ν | H̅(s) | Φ^μ|c⟩'_μ = E(s) c'_ν.

Equation (<ref>) defines a new reference'with expansion coefficientsc'_μ.
This new reference may be used as a starting point for a subsequent MR-DSRG transformation and this procedure can be repeated until convergence and such converged energy is referred as fully relaxed energy.
For the nonperturbative MR-DSRG schemes discussed in this work, we use the fully relaxed energy by default, unless otherwise noted.




 §.§ Simplifying the MR-DSRG equations: Sequential transformation


















Our first modification to the MR-DSRG approach is an alternative way to transform the bare Hamiltonian via a sequence of unitary operators with increasing particle rank

    H̅(s) = e^-Â_n(s)⋯ e^-Â_2(s) e^-Â_1(s)Ĥ e^Â_1(s) e^Â_2(s)⋯ e^Â_n(s).

We term the MR-DSRG approach based on Eq. (<ref>) the sequential MR-DSRG (sq-MR-DSRG), while we refer to the original formalism based on Eq. (<ref>) as the traditional MR-DSRG.
Note that in the limit ofs →∞and no truncation ofÂ(s), both the traditional and sequential MR-DSRG can approach the full configuration interaction limit.<cit.>
However, these schemes are not equivalent for truncatedÂ (s)[for example,n = 2in Eqs. (<ref>) and (<ref>)] due to the fact that operators of different particle rank do not commute, that is,[Â_i (s), Â_j (s)] ≠0fori ≠j.

An advantage of the sq-MR-DSRG approach is thatÂ_1 (s)can be exactly folded into the Hamiltonian via a unitary transformation.
The resultingÂ_1 (s)-dressed Hamiltonian[H̃(s)],

    H̃(s) = e^-Â_1(s)Ĥ e^Â_1(s),

preserves the particle rank of the bare Hamiltonian [Eq. (<ref>)]. The corresponding scalar and tensor components ofH̃(s)can be obtained by a simple unitary transformation of the one- and two-electron integrals (f_p^qandv_pq^rs) and update of the scalar energy.
As will be discussed in Sec. <ref>, theÂ_1 (s)-dressed Hamiltonian can be computed very efficiently when the two-electron integrals are approximated with DF or CD.

The transformed Hamiltonian for the sq-MR-DSRG truncated to one- and two-body operators is given by

    H̅(s)  = e^-Â_2(s)H̃(s) e^Â_2(s).

In the linear approximation, the evaluation of Eq. (<ref>) is simpler than in the traditional MR-LDSRG(2) since the total number of tensor contractions is reduced from 39 to 30.<cit.>
Another advantage of the sequential approach is thatÂ_1 (s)is treated exactly, while in the traditional scheme some contractions involving singles are neglected.
To appreciate this point, consider all the contributions to the double-commutator term in the MR-LDSRG(2) that depend onÂ_1(s)[[Ĥ, Â(s)]_1,2, Â(s)]_0,1,2←       [[Ĥ, Â_1(s)]_1,2, Â_1(s)]_0,1,2
    
    +    [[Ĥ, Â_1(s)]_1,2, Â_2(s)]_0,1,2
    
    +    [[Ĥ, Â_2(s)]_1,2, Â_1(s)]_0,1,2.

The first term on the r.h.s. of Eq. (<ref>) is treated exactly in the MR-LDSRG(2).
However, since contractions involvingÂ_2(s)generate three-body terms  (truncated in the linearized approximation), the contribution ofÂ_1(s)in the second and third terms are not included exactly in the MR-LDSRG(2) transformed Hamiltonian.
In the sequential approach, all contributions fromÂ_1(s)are treated by forming the operatorH̃(s), and since the BCH expansion for such transformation does not generate intermediates with rank greater than two, all terms involvingÂ_1(s)are treated exactly in the linearized approximation.



 §.§ Alleviating the memory bottleneck: The non-interacting virtual orbital (NIVO) approximation


In both the traditional and sequential MR-DSRG approaches, the DF approximation reduces the cost to store both the bare andÂ_1 (s)-dressed Hamiltonian fromO(N^4)toO(N^2M).
However, in the evaluation of the recursive commutator approximation ofH̅(s), two-body operators are generated during the evaluation of each intermediate commutator [Ĉ^(k) (s)] andH̅(s).
These quantities haveO(N^4)storage cost and, thus, reintroduce the bottleneck avoided with DF.

In order to reduce the cost to storeH̅_2(s)andĈ^(k)_2 (s), we shall neglect certain tensor blocks of these operators.
By partitioning of orbitals into core (𝐂), active (𝐀) and virtual (𝐕) spaces, each general 4-index tensor may be subdivided into 81 blocks according to the combination of orbital indices, for example𝐂𝐂𝐂𝐂,𝐀𝐀𝐕𝐕,𝐂𝐀𝐕𝐀, etc.
We propose a non-interacting virtual orbital (NIVO) approximation, which neglects the operator components ofĈ^(k+1)_2 (s) = [Ĉ^(k)_1,2 (s),Â(s)]_2,k ≥0, with three or more virtual orbital indices (𝐕𝐕𝐕𝐕,𝐕𝐂𝐕𝐕,𝐕𝐕𝐕𝐀, etc.) in the recursive definition of the linearized BCH series.
Neglecting these blocks, the number of elements in each NIVO-approximated tensor is reduced fromO(N^4)toO(N^2N_𝐇^2), a size comparable to that of theÂ_2(s)tensor.
For instance, in the cyclobutadiene computation using a quintuple-ζbasis set reported in Sec. <ref>, the memory requirements ofH̅_2(s)orĈ^(k)_2 (s)are reduced from 2.7 TB to 6.8 GB by the NIVO approximation.

To justify the NIVO approximation we analyze its effect on the energy.
The first term in the BCH series that is approximated in the sq-MR-LDSRG(2)+NIVO scheme is the commutatorĈ_2^(1)(s) = [H̃ (s), Â_2 (s)]_2.
Indicating the terms neglected in NIVO asδĈ_2^(1)(s), we see that the first energy contribution affected by the NIVO approximation comes from the expectation value of the triple commutator term [δĈ_0^(3) (s)]

    δĈ_0^(3) (s) = 1/6[[δĈ_2^(1)(s),  Â_2(s)]_1,2, Â_2(s)]_0,

whose contributions are shown as diagrams in Fig. <ref>.
From a perturbation theory perspective, these diagrams are of order four or higher [assumingÂ_2(s)to be of order one] and, therefore, are negligible compared to the leading error (third order) of the linearized commutator approximation.


Hereafter, we shall append “+NIVO” at the end of the method name to indicate the use of NIVO approximation.
For example, the density fitted MR-LDSRG(2) method in the sequential transformation ansatz with NIVO approximation is termed “DF-sq-MR-LDSRG(2)+NIVO”.




§ IMPLEMENTATION


The sq-MR-LDSRG(2) method combined with DF and the NIVO approximation was implemented in Forte,<cit.> an open-source suite of multireference theories for molecular computations.
This implementation reuses several components of our previous MR-LDSRG(2) code based on conventional four-index two-electron integrals.<cit.>
The DSRG equations were implemented as tensor contractions using the Ambit tensor library,<cit.> while integrals were generated using the Psi4 package.<cit.>
In the following, we provide the details of our implementation of the sequential ansatz in combination with DF.<cit.>

 §.§ Sequential transformation


TheÂ_1 (s)-dressed Hamiltonian [H̃(s), Eq. (<ref>)] can be obtained by a unitary transformation ofĤvia the operatorÛ (s) = exp[Â_1(s)].
For clarity, we shall drop the label “(s)” for alls-dependent quantities [H̃ (s), Â_1(s), andÛ (s)] in this section.
The one- and two-body components ofH̃(h̃^q'_p'andṽ_p'q'^r's') are given by

    h̃^q'_p'   = ∑_pq U_q'^q  h^q_p U_p^p',
    ṽ_p'q'^r's'   = ∑_pqrs U_r^r' U_s^s' v_pq^rs U_p'^p U_q'^q.

Here, the unitary matrixUpp'=(𝕌)_p'pand its inverseUp'p = (𝕌)_p'p^∗are given by𝕌 = e^𝔸, where the matrix𝔸is composed of elements of theÂ_1tensor,(𝔸)_ia=taiand(𝔸)_ai=-tai.
Note that we use primed indices only as a way to distinguish labels, yet these indices by no means imply a new set of orbitals.

TheÂ_1-dressed Hamiltonian written in normal ordered form with respect tois given by

    H̃ = Ẽ_0 + ∑_pqf̃_p^q{â_q^p} + 1/4∑_pqrsṽ_pq^rs{â_rs^pq},

where the transformed energy (Ẽ_0) is given by 

    Ẽ_0 = ∑_i'j'^𝐇h̃_i'^j'γ_j'^i' + 1/4∑_i'j'k'l'^𝐇ṽ^k'l'_i'j'γ^i'j'_k'l',

and the Fock matrix elements (f̃^q'_p') are defined as

    f̃^q'_p' = h̃^q'_p' + ∑_i'j'^𝐇ṽ^q'j'_p'i'i'j'.

Note that the quantitiesi'j'andi'j'k'l'in Eqs. (<ref>) and (<ref>) are the untransformed 1- and 2-RDMs of the reference|⟩defined aspq = ⟨| pq ||$⟩ and pqrs = ⟨| pqrs ||$⟩, respectively.

The two-electron integral transformation [Eq. (<ref>)] has a noticeable cost [O(N^5)] and must be repeated each time theÂ_1operator is updated.
However, in the implementation based on DF integrals, this transformation may be performed in a significantly more efficient way.
In DF, the four-index electron repulsion integral tensor as a contraction of a three-index auxiliary tensor (B_pq^Q),

    ⟨pq | rs|≈⟩∑_Q^M B^Q_pr B^Q_qs,

whereMis the dimension of the fitting basis in DF.
Using this decomposition, the unitary transformation may be directly applied to each auxiliary tensor,

    B̃^Q_p'q' = ∑_pq  B^Q_pq  U^p_p' U^q_q',

reducing the cost to evaluateH̃down toO (N^3 M).

Equations (<ref>)–(<ref>) specify the procedures to obtainH̃as a unitary transformation ofĤ.
SinceH̃retains the structure ofĤ, we can reuse most of our previous MR-LDSRG(2) code<cit.> to implement sq-MR-LDSRG(2) by employingH̃(instead ofĤ) and removing terms involvingÂ_1.

As described in Ref. , we evaluate the  commutatorĈ^(k+1) = 1/k + 1[Ĉ^(k)_1,2 , Â]_0,1,2in Eq. (<ref>) using the following recursive system of equations


since[Ĉ^(k)_1,2 , T̂^†] = - [Ĉ^(k)_1,2 , T̂]^†, computing the

    {Ô^(k+1)    = 1/k + 1[Ĉ^(k)_1,2 , T̂]_0,1,2,
    Ĉ^(k+1)   = Ô^(k+1) +  [Ô^(k+1)]^†,
    .

whereÔ^(k+1)is an intermediate containing up to two-body components.
The iteration starts from eitherĈ^(0) = Ĥ, in traditional MR-LDSRG(2), orĈ^(0) = H̃in the sequential version, optionally applying the NIVO approximation to the two-body intermediate tensorsÔ^(k)_2,Ĉ^(k)_2andH̅^(k)_2fork ≥1.



 §.§ Batched tensor contraction for the DF algorithm




Despite the storage cost reduction of the DF and NIVO approximations, another potential memory bottleneck is the size of the intermediate tensors generated during the evaluation of commutators.
For example, consider the following contraction, 

    Orsij←∑_ab^𝐏rsabtabij  ∀ i,j ∈𝐇, ∀ r,s ∈𝐆,

which is also found in the CCSD equations.
In the DF case, Eq. (<ref>) is written as two contractions involving auxiliary tensors,

    Orsij←   ∑_ab^𝐏∑_Q^M B^Q_ar B^Q_bstabij
      - ∑_ab^𝐏∑_Q^M B^Q_as B^Q_brtabij.

The most efficient way to evaluate the first term of Eq. (<ref>) is to introduce the intermediate tensorI_arbs = ∑_Q^M B^Q_ar B^Q_bsof sizeO(N^2 N_𝐏^2).
To avoid storage of these large intermediates, it is common to evaluate Eq. (<ref>) using a batched algorithm, whereby a slice of the tensorI_arbsis computed and contracted on the fly with the amplitudestabij.
To automate this optimization of the tensor contraction we have coded a generic batching algorithm in the tensor library Ambit.<cit.>
Whereas the Ambit code for the first term in Eq. (<ref>) is written as
[fontsize=]
O["ijrs"] += B["Qar"] * B["Qbs"] * t["ijab"];

our new implementation allows batching over the indexrby simply surrounding the contraction with the  function decorator
[fontsize=]
O["ijrs"] += batched("r", B["Qar"] * B["Qbs"] * t["ijab"]);


 §.§ Computational cost reduction


Here we discuss timings for all the MR-LDSRG(2) variants introduced in this work.
In MR-LDSRG(2) theory, the computational bottleneck is forming the DSRG transformed HamiltonianH̅.
Timings for computingH̅ in the case of cyclobutadiene (see Sec. <ref> for details) are summarized in Fig. <ref>.
Detailed timings for the evaluation ofH̅using all the combinations of the approximations considered here are reported in the Supplementary Material.


The total timing (t_tot) for computing the transformed Hamiltonian inniterations is partitioned according to 

    t_tot = t_1 + t_2 + t_misc,

wheret_1andt_2are the timings to evaluate the commutators involvingT̂_1andT̂_2, respectively. In the sequential transformation approach,t_1is instead defined as the time for forming theÂ_1-transformed Hamiltonian.
The termt_miscaccounts for the cost to sort and accumulate the results of contractions withT̂_1andT̂_2, as shown in the second line of Eq. (<ref>).
Figure <ref> shows that the timing for the conventional MR-LDSRG(2) is dominated by contractions involvingT̂_1andT̂_2.
The cost of the singles contractions can be reduced significantly (3–5 times) by employing the sequentially transformed approach, even though at each iteration of the sq-MR-LDSRG(2) equations it is necessary to build the operatorH̃.

Applying the NIVO approximation to the original MR-LDSRG(2) leads to a drastic reduction of the total computational time (×18 speedup).
This reduction in timings is due to several contributing factors.
First, the evaluation of theT̂_1contractions in NIVO is sped up by a factor of𝒪(n N / N_𝐇), wherenis the number of commutators included in the BCH series.
Second, the contributions due to doubles,Ô^(k+1)←[Ĉ^(k), T̂_2], have identical scaling for the first commutator, but fork ≥1they can be evaluated with a speedup of a factor of𝒪(N^2 / N_𝐇^2).
Third, the NIVO approximation also reducest_miscsignificantly because the tensors transpose and accumulation operations costs are reduced fromO(N^4)to𝒪(N^2N_𝐇^2).
For largeN / N_𝐇ratios, the cost to evaluateH̅in the NIVO approximation is dominated by the commutator[Ĥ,T̂_2], with scaling identical to that of CCSD.
For comparison, the similarity-transformed Hamiltonian can be evaluated in 24 s with 's CCSD, in 121 s with our NIVO-MR-LDSRG(2) code, and 2208 s with the original MR-LDSRG(2) code (in both cases employing an unrestricted implementation and C_1symmetry).



In general, we observe an increase int_2due to the extra cost to build two-body intermediates from the auxiliary tensors for methods combined with DF.
The traditional and sequential approaches using the DF/NIVO approximations have similar costs, with the latter being slightly faster due to the efficient transformation of the auxiliary tensors [B̃, Eq. (<ref>)] afforded by the DF approximation.
For this example, the DF-sq-MR-LDSRG(2)+NIVO computation ran 12 times faster than the one using the original approach.
As we will demonstrate in the next section, this method is as accurate as the MR-LDSRG(2) and, therefore, the method we recommend for large-scale computations.



§ RESULTS AND DISCUSSION


 §.§ First row diatomic molecules


We first benchmark the effect of DF and the NIVO approximation on the traditional and sequential versions of the MR-LDSRG(2).
Our test set consists of eight diatomic molecules: BH, HF, LiF, BeO, CO, C2, N2, and F2.
Specifically, we computed equilibrium bond lengths (r_e), harmonic vibrational frequencies (ω_e), anharmonicity constants (ω_e x_e), and dissociation energies (D_0) and compare those to experimental data taken from Ref. .
The dissociation energyD_0includes zero-point vibrational energy corrections that account for anharmonicity effects and is computed asD_0 = D_e - ω_e /2 + ω_e x_e /4(in a.u.), whereD_eis the dissociation energy with respect to the bottom of the potential.<cit.>
Since our current implementation of the MR-DSRG cannot handle half-integer spin states, the energies of the atoms Li, B, C, N, O, and F were computed as half of the energy of the stretched homonuclear diatomic molecule at a distance of 15 Å.
All spectroscopic constants were obtained via a polynomial fit of the energy using nine equally spaced points centered around the equilibrium bond length and separated by 0.2 Å.
For all eight molecules, we adopted a full-valence active space where the 1s orbital of hydrogen, and the 2s and 2p orbitals of first-row elements are considered as active orbitals.
No orbitals were frozen in the CASSCF optimization procedure.
The flow parameter for all DSRG computations was set tos=0.5, as suggested by our previous work.<cit.>
All computations utilized the cc-pVQZ basis set<cit.> and 1s-like orbitals of the first-row elements were frozen in the CC and MR-DSRG treatments of electron correlation.
In DF computations, we employed a mixed flavor of the auxiliary basis sets.
For CASSCF, the cc-pVQZ-JKFIT auxiliary basis set<cit.> was used for H, B, C, N, O and F atoms, and the def2-QZVPP-JKFIT basis set<cit.> was used for Li and Be atoms.
The cc-pVQZ-RI basis set<cit.> was applied to all atoms in DSRG computations.

d[1]D..#1


Figure <ref> and Table <ref> report a comparison of second- and third-order DSRG multireference perturbation theory (DSRG-MRPT2/3), the original MR-LDSRG(2), DF-sq-LDSRG(2)+NIVO, CCSD, and CCSD(T). 
The mean absolute error (MAE) and standard deviation (SD) reported in Table <ref> show that MR-LDSRG(2) method is as accurate as CCSD(T) in predictingr_e,ω_e x_eandD_0, while it predictsω_ethat are of accuracy between that of CCSD and CCSD(T).

The fact that the MR-LDSRG(2) results are more accurate than those from CCSD suggests that the full-valence treatment of static correlation leads to a more accurate treatment of correlation. It is also rewarding to see that in many instances the MR-LDSRG(2) has an accuracy similar to that of CCSD(T), despite the fact that the former does not include triples corrections.


To analyze the impact of each approximation of the MR-LDSRG(2) method, in Table <ref> we report the mean absolute difference between properties computed with and without each approximation.
The use of a sequential ansatz has a modest effect on all properties, with the largest mean absolute differences observed forω_e(1.6 ) andD_0(0.3 ). Nevertheless, the MAE with respect experimental results is nearly unchanged, if not slightly improving.
The DF and NIVO approximations have an effect on molecular properties that is comparable in magnitude and smaller than the deviation introduced by the sequential ansatz.
When these three approximations are combined together, the resulting method shows errors with respect to experimental values that are nearly identical to those from the conventional MR-LDSRG(2).
The only noticeable deviations are found forω_e(MAE 13.9 vs. 14.3 ) andD_0(MAE 2.6 vs. 2.5 ).





In this study, we also investigate the effect of combining the DF-sq-LDSRG(2)+NIVO method with truncation of the BCH expansion, i.e., terminatingH̅(s) = H̃ (s) + ∑_n=1^k 1/n! Ĉ^(n) (s)at a given integerk.
The recursive evaluation ofH̅(s)via Eq. (<ref>) usually requires 10–12.
Truncation of the BCH series to a few terms may therefore introduce speedups of up to 3–4 times.
In Table <ref> and Table <ref> we report statistics computed by approximating the BCH expansion up to  2, 3, and 4 commutators.
The use of only two commutator introduces noticeable deviations with respect to experiments forω_eandD_0. Compared to the full BCH series, this truncation level increases the MAE ofω_eandD_0by 12.6  and0.5, respectively.
The inclusion of the triply-nested commutator significantly reduces these deviations to only 0.1  and 0.2  , respectively.
The four-fold commutator term yieldsr_e,ω_e, andω_e x_ethat are nearly identical to those from the untruncated BCH series, while the MAE ofD_0deviates only by 0.3 .

Since the error introduced by neglecting the four commutator term is smaller or comparable to the other approximations considered here, our results suggest that a BCH series truncated to three commutators may offer a good compromise between accuracy and speed.



 §.§ Cyclobutadiene


Next, we consider the automerization reaction of cyclobutadiene (CBD, C_4H_4). We study the energy difference between the rectangular (D_2h) energy minimum and the square transition state (D_4h).<cit.>
This reaction is a challenging chemistry problem for both experiment and theory.<cit.>
Due to its instability, there are no direct measurements of the reaction barrier, and  experiments performed on substituted cyclobutadienes  suggest the barrier height falls in the range 1.6–10 .<cit.>
In this work, we optimized the equilibrium and transition state geometries using finite differences of energies to compute the barrier height.
Specifically, we compare both DF-MR-LDSRG(2)+NIVO and DF-sq-MR-LDSRG(2)+NIVO optimized geometries to those obtained from the state-specific MRCC of Mukherjee and co-workers (Mk-MRCC)<cit.> as implemented in .
<cit.>

To reduce computational cost, all MR-DSRG calculations performed two steps of the reference relaxation procedure discussed in Sec. <ref>.
A comparison of this procedure with full reference relaxation using the cc-pVDZ basis set shows errors of ca. 0.01  for absolute energies, 0.0001 Å for bond lengths, and 0.001^∘for bond angles.
We applied a Tikhonov regularization denominator shift<cit.> of 1 m in all Mk-MRCC calculations to guarantee convergence.
The Mk-MRCC implementation used in this work neglects effective Hamiltonian couplings between reference determinants that differ by three or more spin orbitals, and therefore yield approximate results when applied to the CAS(4e,4o) reference considered here.
All computations utilized the cc-pVXZ (X=D, T, Q, 5) basis set<cit.>, and the corresponding cc-pVXZ-JKFIT<cit.> and cc-pVXZ-RI<cit.> auxiliary basis sets for DF-CASSCF and DF-DSRG computations, respectively.
The 1s core electrons of carbon atoms were frozen in all post-CASSCF methods.
All results were computed using semi-canonical CASSCF orbitals.



Preliminary computations using the cc-pVDZ basis using the CAS(2e,2o) and CAS(4e,4o) active spaces revealed an interesting aspect of this system.
As shown in Fig. <ref>, thes-dependency of the automerization barrier displays significantly different behavior for these two active spaces.
In both cases, the predicted activation energies change significantly for small values ofs(<0.2), a normal trend observed for all DSRG computations and due to the increased recovery of dynamical correlation energy.
Interestingly, while the CAS(4e,4o) curve flattens out for larger values ofs, the CAS(2e,2o) curve shows a significants-dependence in the ranges ∈[0.5,8].



To understand the origin of this difference we analyze the double substitution amplitudes [tabij(s)] as a function ofsfor both the equilibrium and transition state geometries, as shown in Fig. <ref>.
In the CAS(2e,2o) case, we notice some abnormally large amplitudes (indicated in red), some of which are as large as 0.1. These amplitudes correspond to excitations within the fourπorbitals of CBD, and suggests that the CAS(2e,2o) space is insufficient to capture all static correlation effects in CBD.



The offending amplitudes converge at different rates assincreases, and introduce a strongs-dependence in the energy barrier. Note also that in the limit ofs →∞there is a significant difference in the barrier for the CAS(2e,2o) and CAS(4e,4o) spaces.

In contrast, in the CAS(4e,4o) computations all excitations within theπorbitals are included in the active space and the resulting DSRG amplitudes have absolute values less than 0.05.
Diverging amplitudes in computations with CAS(4e,4o) reference wave functions, corresponding to intruder states, can also be seen in Fig. <ref>.
Our results reported in Table <ref> are all based on the flow parameter values=1.0, which is significantly far from the region  (s > 5.0) where amplitudes begin to diverge, and at the same time leads to well converged absolute energies (see Figure S3 in the Supporting Information).
We also performed computations usings= 0.5  (reported in Table S4) to verify that the automerization energies computed with different values of the flow parameter are consistent.
In general, the difference in automerization energies computed withs = 0.5and 1.0  is of the order of 0.6–0.7 .
Note that intruder states are also encountered in Mk-MRCCSD computations based on the CAS(4e,4o) reference, and lead to convergence issues that could be avoided only via Tikhonov regularization.



Geometric parameters for the optimized structures and energy barriers of CBD computed with the CAS(4e,4o) reference are reported in Table <ref>.
A comparison the energy barrier computed at the CASSCF and correlated levels shows that dynamical correlation is important in this system as it increases is by about 1–3 .
Our best estimate for the automerization barrier of CBD is 10.3  at the DF-sq-MR-LDSRG(2)+NIVO/cc-pV5Z level of theory.
This value is likely to be slightly higher then the exact results since in the Mk-MRCC results perturbative triples corrections contribute to lowering the barrier by ca. 1.5 .
In general, the MR-LDSRG(2) results are between those of Mk-MRCCSD and Mk-MRCCSD(T), reinforcing the same observation we made in the benchmark of diatomic molecules.
For instance, the MR-LDSRG(2) predicted C-C bond length at theD_4hgeometry is 1.4446 Å, which is almost midway between the Mk-MRCCSD (1.4406 Å) and Mk-MRCCSD(T) (1.4480 Å) values.
As expected, the differences between the conventional and sequentially transformed MR-DSRG(2) results are negligible.
Our results in cc-pVDZ and cc-pVTZ bases also agree well with other reported theoretical values, especially those computed with multireference methods,<cit.> and the experimental range reported in Ref. .
Using our new implementation, we were able to perform, for the first time, nonperturbative multireference computations on cyclobutadiene using the cc-pV5Z basis (580 correlated orbitals) on a single node with 128 GB of memory.




§ CONCLUSION


In this work, we describe a strategy to reduce the computational and memory costs of the multireference driven similarity renormalization group (MR-DSRG). We demonstrate that the cost of the linear MR-DSRG with singles and doubles [MR-LDSRG(2)] can be lowered substantially without compromising its accuracy by using a combination of: 1) a sequential unitary transformation, 2) density fitting (DF) of the  two-electron integrals, and 3) the non-interacting virtual orbital (NIVO) operator approximation.
The sequential MR-DSRG scheme introduced in this work [sq-MR-DSRG] reduces the cost of evaluating single-excitations and allows to treat them exactly. Like in the case of Brueckner coupled cluster theory,<cit.> this approach reduces the number of algebraic terms in the DSRG equations because there are no terms (diagrams) containing single excitations.
The use of DF integrals reduces the memory requirements of the original MR-DSRG(2) from𝒪(N^4)toO(N^2M), whereNis the number of basis functions.
Density fitting is particularly convenient when combined with the sequential approach because the contributions of singles can be directly included in the DF auxiliary three-index integrals, reducing the integral transformation cost fromO (N^5)toO (N^3 M), whereMis the number of auxiliary basis functions.

The NIVO approximation neglects the operator components of a commutator with three or more virtual indices.
A formal analysis of this approximation showed that the leading error is of fourth order in perturbation theory.
In practice, NIVO is crucial to both avoiding the memory bottleneck of the MR-DSRG(2) and reducing the computational cost to evaluate the transformed Hamiltonian.

To benchmark the MR-LDSRG(2) and sq-MR-LDSRG(2) approaches and assess the impact of the DF and  NIVO approximations, we computed the spectroscopic constants of eight diatomic molecules using the full-valence active space and the cc-pVQZ basis set.
Compared to experimental data, both MR-DSRG methods yield results that are as accurate as those obtained with CCSD(T).
Moreover, the DF-sq-MR-LDSRG(2)+NIVO results are almost identical to those computed without the NIVO approximation: the harmonic vibrational frequencies, anharmonicity constants, and dissociation energies only differ by, on average, 0.1 , 0.1 , and 0.2 , respectively.
These results supports our claim that the speedup brought by the NIVO approximation does not sacrifice the accuracy of both variants of the MR-LDSRG(2).

Combining DF and the NIVO approximation, both the traditional and sequential MR-LDSRG(2) can be routinely applied to chemical systems with more than 500 basis function.
We demonstrate this point by studying the automerization reaction of cyclobutadiene using a quintuple-ζbasis set (584 basis functions).
Our best estimate of the reaction barrier from DF-sq-MR-LDSRG(2)+NIVO/cc-pV5Z is 10.3 . However, we expect that this result is likely overestimated due to the lack of three-body corrections in the MR-LDSRG(2) theory.
Our results agree well with Mk-MRCCSD predictions and multireference coupled cluster reported in the literature.

In conclusion, we have shown that it is possible to significantly reduce the cost of MR-LDSRG(2) computations without reducing the accuracy of this approach. The sequential approach and NIVO approximations are general, and can be applied to improve the efficiency of other unitary nonperturbative methods (e.g., unitary coupled cluster theory) and downfolding schemes for classical-quantum hybrid algorithms.<cit.>































§ SUPPORTING INFORMATION



Detailed timings and formal scaling required to evaluate the transformed Hamiltonian, error statistics for the eight diatomic molecules, cyclobutadiene automerization energies reported in the literature, and convergence of DF-sq-MR-LDSRG(2)+NIVO absolute energy with respect to flow parameter.





§ ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy under Award No. DE-SC0016004, a Research Fellowship of the Alfred P. Sloan Foundation, and a Camille Dreyfus Teacher-Scholar Award.