We describe a formulation of multi-reference perturbation theory that obtains a rigorous upper bound to the second order energy by minimizing the Hylleraas functional in the space of matrix product states (MPS). The first order wavefunctions so obtained can also be used to compute the third order energy with little overhead. Our formulation has several advantages including (i) flexibility with respect to the choice of zeroth order Hamiltonian, (ii) recovery of the exact uncontracted multi-reference perturbation theory energies in the limit of large MPS bond dimension, (iii) no requirement to compute high body density matrices, (iv) an embarrassingly parallel algorithm (scaling up to the number of virtual orbitals, squared, processors). Preliminary numerical examples show that the MPS bond dimension required for accurate first order wavefunctions scales sub-linearly with the size of the basis.

Commonly, electron correlation is classified into static (or strong) correlation and dynamic (or weak) correlation. The former is essential to capture the qualitative electronic structure, and the latter to generate quantitatively accurate results. In most chemical systems, static correlation involves a small subset of orbitals with degenerate orbital energies on the energy scale of the Coulomb repulsion. In this subspace, also known as the active space, multi-configuration self-consistent field (MCSCF) calculations have traditionally been performed. However, the cost of (numerically exact) MCSCF scales exponentially with the number of active orbitals and is thus restricted to active spaces with about 16 electrons in 16 orbitals. In the last decade, this active space restriction has effectively been removed, without significant numerical errors, with the advent of new near-exact methods such as the density matrix renormalization group (DMRG)1–8 and full configuration interaction quantum Monte Carlo (FCIQMC).9–11 (More approximate techniques for large active spaces, such as restricted active space and general active space methods,12,13 high-order active-space coupled cluster,14–17 and variational reduced density matrix methods,18–20 have also been advanced.) However, for quantitative chemical accuracy, these new active space methods must still be augmented with techniques to include dynamic correlation, and this remains an important frontier.21–25 

The primary technical difficulty in treating dynamic correlation is the large number of orbitals require to converge the short-range electron-electron cusp. Explicit correlation (through R12 and F12 factors)26–31 significantly ameliorates, but does not eliminate, the need for large basis sets. One affordable approach to include dynamic correlation from a large number of external orbitals is through perturbation theory. Examples of such multi-reference perturbation methods include complete active space perturbation theory (CASPT2)32,33 and n-electron valence state perturbation theory (NEVPT2).34,35 However, even when recently combined with DMRG reference functions,21,36 these methods have only been used with moderate active space sizes. The key bottleneck is the need to construct intermediates involving three- and four-body reduced density matrices (RDMs) in the active space. As the number of active orbitals increases, calculating and storing these RDMs become prohibitively expensive.

In this Communication, we describe an alternate formulation of a multi-reference perturbation theory, which we call matrix product state perturbation theory (MPS-PT). The new formulation possesses substantial advantages as follows:

  1. It bypasses the need for high-body RDMs in the active space, potentially allowing very large active spaces to be used.

  2. Perturbation theory depends strongly on the zeroth-order Hamiltonian H0. Changing H0 (as in CASPT2 versus NEVPT2) usually requires re-deriving and re-implementing non-trivial intermediates. In our formulation we can easily change H0 without significantly changing the implementation.

  3. To reduce the cost of multi-reference perturbation theory, the first order wavefunction is often determined in a restricted space, as in internally contracted CASPT2 or partially and strongly contracted NEVPT2 (PC-NEVPT2, SC-NEVPT2). Aggressive contraction can lead to additional errors which are difficult to estimate a priori. We also contract the first order space via the MPS bond dimension, but do so in an automatic and optimal way. This contraction systematically and rapidly converges to the uncontracted result with increasing bond dimension, and the errors of contraction can be robustly estimated.

  4. The algorithm is highly parallel and has a modest scaling with the MPS bond dimension and size of basis.

Our current implementation is significantly slower than existing NEVPT2 and CASPT2 implementations for small active spaces but in contrast allows very large active spaces to be used.

It is well known that perturbation theory can be formulated as a variational problem.37,39 For the second order energy, the variational functional is the Hylleraas functional H1], Eq. (1), which is minimized with respect to the first order wavefunction |Ψ1⟩,

\begin{align}H[\Psi _1] = \langle \Psi _1| H_0 - E_0|\Psi _1\rangle + 2 \langle \Psi _1|QV|\Psi _0\rangle .\end{align}
H[Ψ1]=Ψ1|H0E0|Ψ1+2Ψ1|QV|Ψ0.
(1)

Here H0, E0, and |Ψ0⟩ are, respectively, the zeroth order Hamiltonian, energy, and wavefunction, V is the perturbation (HH0), and Q the projector onto the space orthogonal to |Ψ0⟩. It can be easily verified that at the minimum, the wavefunction |Ψ1⟩ satisfies the familiar equation

\begin{align}(H_0 - E_0)|\Psi _1\rangle = -QV|\Psi _0\rangle .\end{align}
(H0E0)|Ψ1=QV|Ψ0.
(2)

In MPS-PT2 we minimize the Hylleraas functional, while expressing |Ψ1⟩ as a MPS. MPS forms a complete variational space, and the quality of the MPS solution is controlled by a single parameter “M,” the dimension of the auxiliary bond in the MPS (see Figure 1). Overlaps (⟨Ψ12⟩) and transition elements of operators between two MPSs (⟨Ψ1|O2⟩) can be evaluated in O(M3) central processing unit (CPU) time. The Hylleraas functional can be minimized with respect to a MPS of arbitrary bond dimension using a sweep algorithm analogous to that in the DMRG.40 In the limit of large M, the solution converges to the exact (uncontracted) first order wavefunction (and second order energy). From the first order wavefunction, the third order correction can also be calculated due to Wigner's 2n + 1 rule, as

\begin{align}E_3 = \langle \Psi _1|V|\Psi _1\rangle - E_1 \langle \Psi _1|\Psi _1\rangle .\end{align}
E3=Ψ1|V|Ψ1E1Ψ1|Ψ1.
(3)

A general MPS representing |Ψ1⟩ is shown in Eq. (4), where ni is the occupation of orbital i,

\begin{align}|\Psi _1\rangle = \sum _{\lbrace n\rbrace , i_i\cdots i_{k-1}} A^{n_1}_{i_1} A^{n_2}_{i_1 i_2} \ldots A^{n_k}_{i_{k-1}} |n_1 n_2 \ldots n_k\rangle .\end{align}
|Ψ1={n},iiik1Ai1n1Ai1i2n2...Aik1nk|n1n2...nk.
(4)

Eq. (4) is conveniently expressed graphically as a tensor network, as shown in panel (a) of Figure 1 (see Refs. 42–44 for a detailed explanation of the graphical notation). Similarly, any operator Ω can be written in matrix product operator (MPO) form,

\begin{align}\Omega \!=\!\!\! \sum _{\lbrace n\rbrace \lbrace n^{\prime }\rbrace , i_i\cdots i_{k-1}} \!\!\! W^{n_1 n_1^{\prime }}_{i_1} W^{n_2 n_2^{\prime }}_{i_1 i_2} \ldots W^{n_k n_k^{\prime }}_{i_{k-1}} |n_1 n_2 \ldots n_k\rangle \langle n_1^{\prime } n_2^{\prime } \ldots n_k^{\prime }|,\end{align}
Ω={n}{n},iiik1Wi1n1n1Wi1i2n2n2...Wik1nknk|n1n2...nkn1n2...nk|,
(5)

and this is also expressed graphically in panel (a) of Figure 1. With this notation, the Hylleraas functional is shown in panel (b) of Figure 1, where the operators and wavefunctions have been represented as MPOs and MPSs, respectively, and the corresponding tensor networks are contracted to obtain the final expression.

FIG. 1.

(a) A matrix product state (MPS) can be represented graphically using a series of three-dimensional tensors, shown with circles, each of which is associated with an orbital. The free index, also known as the physical index (pointing upwards) of the tensors, denotes the occupation of the orbital and the other two indices, known as auxiliary indices, are sequentially contracted. The dimension d of the physical index is 4 for a spatial orbital, and the dimension M of the auxiliary index can be increased to make the MPS arbitrarily flexible. Similarly, a matrix product operator (MPO) can be represented as a series of four-dimensional tensors, with two physical indices and auxiliary indices contracted sequentially. Due to the 2-body nature of the Hamiltonians in quantum chemistry, the auxiliary bond dimension of a MPO is always less than k2, where k is the number of orbitals. (b) and (c) The tensors in wavefunctions Ψ1 and Ψ0 are represented by circles and triangles, respectively. The zeroth order Hamiltonian H0 and perturbation operator V are represented by squares and diamonds, respectively. Panel (b) shows the three terms in the Hylleraas functional. Panel (c) shows the partial derivative of the Hylleraas functional

$\partial H[\Psi _1]/\partial A^{n_l}_{i_{l-1}i_l}$
H[Ψ1]/Ail1ilnl that is set to zero in the optimization. The red tensor,
$A^{n_l}_{i_{l-1}i_l}$
Ail1ilnl
is the quantity being solved for. (d) Contractions that need to be carried out at each step of the sweep iteration. The corresponding costs (assuming H0 is Dyall's Hamiltonian41) are shown. For each term, individual contractions are numbered according to the order in which they are performed. If two contractions have the same number, it means that the two indices on the same tensor are fused to form a larger composite index and then the contraction is performed. The solid tensors represent the results of contracting all the open tensors in panel (c) that are not shown in this figure.

FIG. 1.

(a) A matrix product state (MPS) can be represented graphically using a series of three-dimensional tensors, shown with circles, each of which is associated with an orbital. The free index, also known as the physical index (pointing upwards) of the tensors, denotes the occupation of the orbital and the other two indices, known as auxiliary indices, are sequentially contracted. The dimension d of the physical index is 4 for a spatial orbital, and the dimension M of the auxiliary index can be increased to make the MPS arbitrarily flexible. Similarly, a matrix product operator (MPO) can be represented as a series of four-dimensional tensors, with two physical indices and auxiliary indices contracted sequentially. Due to the 2-body nature of the Hamiltonians in quantum chemistry, the auxiliary bond dimension of a MPO is always less than k2, where k is the number of orbitals. (b) and (c) The tensors in wavefunctions Ψ1 and Ψ0 are represented by circles and triangles, respectively. The zeroth order Hamiltonian H0 and perturbation operator V are represented by squares and diamonds, respectively. Panel (b) shows the three terms in the Hylleraas functional. Panel (c) shows the partial derivative of the Hylleraas functional

$\partial H[\Psi _1]/\partial A^{n_l}_{i_{l-1}i_l}$
H[Ψ1]/Ail1ilnl that is set to zero in the optimization. The red tensor,
$A^{n_l}_{i_{l-1}i_l}$
Ail1ilnl
is the quantity being solved for. (d) Contractions that need to be carried out at each step of the sweep iteration. The corresponding costs (assuming H0 is Dyall's Hamiltonian41) are shown. For each term, individual contractions are numbered according to the order in which they are performed. If two contractions have the same number, it means that the two indices on the same tensor are fused to form a larger composite index and then the contraction is performed. The solid tensors represent the results of contracting all the open tensors in panel (c) that are not shown in this figure.

Close modal

In the sweep algorithm, each tensor [An] of the MPS wavefunction |Ψ1⟩ is optimized in sequence, keeping the other tensors fixed. This converts the nonlinear optimization of the MPS (which is multi-linear in its parameters) to a series of linear equations that we solve using the conjugate gradient algorithm. At step l of the sweep, we take the partial derivative of the Hylleraas functional with respect to the set of tensor elements

$A^{n_l}_{i_{l-1} i_l}$
Ail1ilnl of |Ψ1⟩, and set it to zero. The resulting linear equation is shown graphically in panel (c) of Figure 1 and has 3 terms, corresponding to the partial derivatives of the three terms in the Hylleraas functional. The first involves a contraction of
$A^{n_l}_{i_{l-1} i_l}$
Ail1ilnl
(the quantity being solved for, shown with a red circle) with a 6-index tensor, formed by contracting all tensors in the operator H0E0 and |Ψ1⟩ with
$A^{n_l}_{i_{l-1} i_l}$
Ail1ilnl
removed. Similarly, the second and third terms involve the tensors of V, |Ψ0⟩ and all tensors of |Ψ1⟩ with
$A^{n_l}_{i_{l-1} i_l}$
Ail1ilnl
removed. The resulting linear equation for
$A^{n_l}_{i_{l-1} i_l}$
Ail1ilnl
is in O(dM2) unknowns. Solving it by conjugate gradient requires only O(M3) time (due to the special structure of the operator-vector product) and yields the current best estimate of the tensor
$A^{n_l}_{i_{l-1} i_l}$
Ail1ilnl
. It is clear that the tensor so obtained depends on the values of all the other tensors [An] in |Ψ1⟩. Thus we sweep through all the orbitals of |Ψ1⟩ solving a linear algebra problem at each step l = 1…k to optimize each tensor, while keeping the others fixed. This algorithm is essentially the same as the standard DMRG sweep algorithm, the only difference being that instead of an eigenvalue problem, a linear equation is solved at each step.

The computational cost of the above sweep can be expressed in detail in terms of the auxiliary bond dimensions of |Ψ1⟩ and |Ψ0⟩ (M1 and M0, respectively) and the number of active and external orbitals (ka and kv, respectively). For concreteness, we consider H0 to be Dyall's Hamiltonian,41 as used in NEVPT2. All the contractions to be performed at each sweep iteration, together with their respective (leading order in M, ka, kv) computational costs, are shown in panel (d) of Figure 1. Contractions 1, 2, and 3 of the first two terms have the same leading scaling with respect to the number of orbitals even though the operator in the first term (denoted by squares) is a Dyall Hamiltonian H0, which does not couple the active and external spaces, and the operator V in the second term (denoted by diamonds) couples the two spaces. The low scaling of the second term is because components of V with more than 2 external indices do not contribute to the second order energy and can be neglected. The resulting MPO for V, under this restriction, can then always be represented with

$O(k_a^2)$
O(ka2) auxiliary bond dimension, the same as the Dyall H0. The costs for contractions 1–5 in the third term are not shown because they are used to construct the bottom diagram, which is a constant that does not change during the sweep, and thus has a subleading dependence on kv. The final computational cost of a single sweep is
$O(k_vk_a^2M_1^3 + k_vk_a^2M_0M_1^2 + k_vk_a^2M_0^2M_1)$
O(kvka2M13+kvka2M0M12+kvka2M02M1)
, where we have multiplied the most expensive costs in each individual step shown in panel (d) of Figure 1 by kv (the number of virtual orbitals) to reflect the total number of steps in the sweep, assuming kv + kaO(kv). We have also dropped the dependence on d, which is a fixed constant (usually 4). (For brevity, we have not considered the cost of terms with subleading (lower than cubic) dependence on M. These only become important when kv is very large. For example, forming the MPO for V (through DMRG complementary operator techniques45) costs
$O(k_v^2k_a^2M_0M_1)$
O(kv2ka2M0M1)
/sweep, which becomes important when kv > M.)

The above costs assume the Dyall H0. A similar analysis applies to the CASPT2 H0, which takes the form of a one-body Fock operator multiplied by many-body projectors. The many-body projectors lead to additional overlap computations of the form ⟨Ψ01⟩, but these only cost

$O(k_aM_0^2M_1 + k_aM_0M_1^2)$
O(kaM02M1+kaM0M12) for the entire sweep, which is independent of kv. The one-body Fock operator leads to a lower overall computational cost for the first term in the Hylleraas functional involving H0. However, the second term involving V results in the same leading cost as for NEVPT2 above. Additionally, note that because we do not proceed via special reduced density matrix intermediates, the computations for both the Dyall H0 in NEVPT and for the CASPT2 H0 involve the same basic contractions of MPOs and MPSs and thus it is simple to change the zeroth order Hamiltonian in the implementation. This is similar to the simplicity of implementing fully uncontracted determinant algorithms, except here contraction is automatically provided by the finite bond dimension of the MPS. Here, we would also like to point out that this approach can be extended to excited state calculations by a small modification to the current procedure.38 

We have implemented the above MPS-PT sweep algorithm to evaluate the NEVPT2 and NEVPT3 energy in the BLOCK code.46 For computational robustness, we have implemented the two-site version of the sweep algorithm (see Refs. 40 and 42 for a discussion of the difference between one-site and two-site sweep algorithms). In the limit of large M1, our MPS-PT algorithm will yield the uncontracted NEVPT2 and NEVPT3 energies. We have parallelized the algorithm exactly in the manner of the DMRG module of the BLOCK code, where each auxiliary bond contributing to the MPO (arising from a term in the operator/complementary operator product) is evaluated on a different processor. As explained previously, the auxiliary bond dimension of both H0 and V is

$O(k_a^2)$
O(ka2)⁠, thus the code is expected to efficiently parallelize over
$O(k_a^2)$
O(ka2)
processors. We refer the reader to the detailed description of DMRG parallelization in Ref. 47. There is an additional “embarrassing” parallelization over
$O(k_v^2)$
O(kv2)
processors which we mention (but have not implemented). This arises by breaking the summation in V into
$k_v^2$
kv2
individual terms Vi. Each Vi produces a mutually orthogonal state when acting on |Ψ0⟩ and thus the Hylleraas functional can be computed and minimized for each Vi separately, giving a final second-order energy that is just the sum of these terms.

Symmetries of the Hamiltonian are easily utilized in an analogous fashion to the DMRG module of the BLOCK code, which currently implements many symmetries including particle number, Abelian and non-Abelian point groups, and SU(2) (spin symmetry).

We next illustrate various aspects of the theory by performing calculations on the N2 bond dissociation benchmark using a full-valence CAS of 10 electrons in 8 orbitals and the cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets.48 Table I shows the energies of the N2 molecule with various NEVPT2 and MPS-PT variants at several bond-lengths in the cc-pVDZ basis set. The MPS-PT2 energy is converged to a few tens of μEh with M1 = 800 and always gives an energy lower than the partially contracted NEVPT2 theory (since it is converging to the fully uncontracted result) although the energy difference is not large. The error of various standard flavors of contracted NEVPT2 theory, PC-NEVPT2 and SC-NEVPT2 (evaluated using Molpro49 and MPS-PT2 calculations with smaller M1 = 100–300, are plotted against increasing N2 bond-length in Figure 2. We see that the value of M1 required to accurately describe the first order wavefunction |Ψ1⟩ increases somewhat with increasing bond length. This can be understood by realising that the first-order wavefunction contains contributions from

$O(N_a k_a^2k_v^2)$
O(Naka2kv2) determinants (where Na is the number of significant determinants in the active space). As Na increases at increased bond-lengths, the first-order wavefunction becomes more complicated and requires slightly larger M. Nonetheless,
$M\ll N_a k_a^2 k_v^2$
MNaka2kv2
and except for M = 100, the MPS-PT error remains small at all bond-lengths.

Table I.

Energy (E+108.0 in Eh) for a MPS-PT calculation performed with a cc-pVDZ basis set for the 1Σg state of N2 molecule at various bond-lengths. A (10e, 8o) active space was used. PC-NEVPT2, and SC-NEVPT2 results are given for comparison. The last four columns show the MPS-PT2 energy calculated using different auxiliary bond dimensions M1.

   MPS-PT3M1 used in MPS-PT2
re (Å)PC-NEVPT2SC-NEVPT2M1 = 800800300200100
1.8 −1.1430 −1.1396 −1.1626 −1.1432 −1.1431 −1.1428 −1.1403 
2.0 −1.2420 −1.2389 −1.2620 −1.2422 −1.2421 −1.2418 −1.2382 
2.2 −1.2487 −1.2458 −1.2689 −1.2490 −1.2487 −1.2483 −1.2432 
2.4 −1.2129 −1.2103 −1.2333 −1.2132 −1.2129 −1.2123 −1.2066 
2.7 −1.1342 −1.1318 −1.1549 −1.1345 −1.1341 −1.1332 −1.1258 
3.0 −1.0591 −1.0568 −1.0807 −1.0594 −1.0589 −1.0578 −1.0484 
3.6 −0.9648 −0.9631 −0.9896 −0.9651 −0.9641 −0.9624 −0.9430 
4.2 −0.9350 −0.9342 −0.9621 −0.9352 −0.9341 −0.9324 −0.9055 
   MPS-PT3M1 used in MPS-PT2
re (Å)PC-NEVPT2SC-NEVPT2M1 = 800800300200100
1.8 −1.1430 −1.1396 −1.1626 −1.1432 −1.1431 −1.1428 −1.1403 
2.0 −1.2420 −1.2389 −1.2620 −1.2422 −1.2421 −1.2418 −1.2382 
2.2 −1.2487 −1.2458 −1.2689 −1.2490 −1.2487 −1.2483 −1.2432 
2.4 −1.2129 −1.2103 −1.2333 −1.2132 −1.2129 −1.2123 −1.2066 
2.7 −1.1342 −1.1318 −1.1549 −1.1345 −1.1341 −1.1332 −1.1258 
3.0 −1.0591 −1.0568 −1.0807 −1.0594 −1.0589 −1.0578 −1.0484 
3.6 −0.9648 −0.9631 −0.9896 −0.9651 −0.9641 −0.9624 −0.9430 
4.2 −0.9350 −0.9342 −0.9621 −0.9352 −0.9341 −0.9324 −0.9055 
FIG. 2.

Error in the second-order energy (in Eh) for various flavors of NEVPT2 and MPS-PT2 calculations performed with cc-pVDZ basis sets for the 1Σg state of N2 molecule at several bond lengths. A (10e, 8o) active space was used for these calculations. The graph shows that the value of M1 required to capture the correct second-order energy using the MPS-PT2 theory increases with the N2 bond length, most likely due to increase in the multi-reference character of the wavefunction. Nonetheless, for M1 > 100, the error remains quite small at all bond-lengths.

FIG. 2.

Error in the second-order energy (in Eh) for various flavors of NEVPT2 and MPS-PT2 calculations performed with cc-pVDZ basis sets for the 1Σg state of N2 molecule at several bond lengths. A (10e, 8o) active space was used for these calculations. The graph shows that the value of M1 required to capture the correct second-order energy using the MPS-PT2 theory increases with the N2 bond length, most likely due to increase in the multi-reference character of the wavefunction. Nonetheless, for M1 > 100, the error remains quite small at all bond-lengths.

Close modal

When one views the DMRG algorithm in the traditional renormalization group language, each tensor in the MPS (in canonical form) can be viewed as mapping from 4M states (for d = 4) to the M most significant renormalized states. In the two-site variant of the algorithm, the weight of these discarded 3M states measures the severity of the truncation performed in the renormalization step. The discarded weight goes to zero in the limit of large M and for smaller M usually exhibits a linear relationship with the energy error.40,50 Figure 3 shows that after an initial rapid relaxation in the energy, the error in the MPS-PT2 indeed scales linearly with the corresponding discarded weight in

$\mathinner {|{\Psi _1}\rangle }$
|Ψ1⁠. Thus the relationship between the energy error and discarded weight can be used to gauge the convergence of the MPS-PT2 energy with M, just as in the DMRG algorithm.

FIG. 3.

Error in the second-order energy (in Eh) versus the discarded weight for a MPS-PT2 calculation performed with cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets for the 1Σg state of N2 molecule at a bond length of 2 Å. A (10e, 8o) active space was used for these calculations. After an initial rapid relaxation of the energy one recovers the linear relationship between the energy and the discarded weight.

FIG. 3.

Error in the second-order energy (in Eh) versus the discarded weight for a MPS-PT2 calculation performed with cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets for the 1Σg state of N2 molecule at a bond length of 2 Å. A (10e, 8o) active space was used for these calculations. After an initial rapid relaxation of the energy one recovers the linear relationship between the energy and the discarded weight.

Close modal

Figure 4 shows the relationship between the size of the basis set and the value of M1 required to obtain a second-order energy within 1 mEh of the exact fully uncontracted MPS-PT2 energy. We see that, encouragingly, the M1 required scales sub-linearly with the size of the basis. Indeed, the entanglement in the basis set limit is bounded by a constant (the dimension of the active orbital Hilbert space). Thus, M1 must asymptotically be independent of kv.

FIG. 4.

The size of the auxiliary bond dimension M of the first order wavefunction required to obtain a second-order energy within 1 mEh of the exact fully uncontracted MPS-PT2 energy.

FIG. 4.

The size of the auxiliary bond dimension M of the first order wavefunction required to obtain a second-order energy within 1 mEh of the exact fully uncontracted MPS-PT2 energy.

Close modal

To summarize, we have shown that minimizing the Hylleraas functional in the space of wavefunctions described by a MPS can be used to obtain accurate upper bounds to the exact multi-reference second order perturbation theory energy. With little additional computational overhead, third order corrections to the energies can also be obtained. No active space density matrices are required, allowing large active spaces to be treated, and our implementation can be easily modified to accommodate any reasonable zeroth order hamiltonian. Parallelism is also easily exploited. In benchmark calculations we have shown that this approach efficiently recovers the uncontracted perturbation theory energy, beyond existing internally contracted multi-reference perturbation results. Finally, we observe that the bond dimension of the first order wavefunction required to obtain chemical accuracy scales sub-linearly with the size of the basis set, which is promising for future applications with very large basis sets.

This work was supported by the US National Science Foundation (NSF) through Grant No. NSF-CHE-1265277. Additional support for software development was provided through Grant No. NSF-OCI-1265278. S.S. would also like to thank Brecht Verstichel for many helpful discussions.

1.
S. R.
White
and
R. L.
Martin
,
J. Chem. Phys.
110
,
4127
(
1999
).
2.
G. K.-L.
Chan
and
S.
Sharma
,
Annu. Rev. Phys. Chem.
62
,
465
(
2011
).
3.
Ö.
Legeza
,
R.
Noack
,
J.
Sólyom
, and
L.
Tincani
,
Computational Many-Particle Physics
,
Lecture Notes in Physics
Vol.
739
, edited by
H.
Fehske
,
R.
Schneider
, and
A.
Weie
(
Springer
,
Berlin, Heidelberg
,
2008
), pp.
653
664
.
4.
K. H.
Marti
and
M.
Reiher
,
Phys. Chem. Chem. Phys.
13
,
6750
(
2011
).
5.
Y.
Kurashige
and
T.
Yanai
,
J. Phys. Chem.
130
,
234114
(
2009
).
6.
S.
Wouters
and
D.
Van Neck
, “
The density matrix renormalization group for ab initio quantum chemistry
,” e-print arXiv:1407.2040v1 (
2014
).
7.
Y.
Kurashige
,
G. K.-L.
Chan
, and
T.
Yanai
,
Nat. Chem.
5
,
660
(
2013
).
8.
S.
Sharma
,
K.
Sivalingam
,
F.
Neese
, and
G. K.-L.
Chan
,
Nat. Chem.
(published online
2014
).
9.
G. H.
Booth
,
D.
Cleland
,
A.
Alavi
, and
D. P.
Tew
,
J. Chem. Phys.
137
,
164112
(
2012
).
10.
G. H.
Booth
,
A. J. W.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
131
,
54106
(
2009
).
11.
G. H.
Booth
,
S. D.
Smart
, and
A.
Alavi
,
Mol. Phys.
112
,
1855
(
2014
).
12.
J.
Olsen
,
B. O.
Roos
,
P.
Jørgensen
, and
H. J. A.
Jensen
,
J. Chem. Phys.
89
,
2185
(
1988
).
13.
P. A.
Malmqvist
,
A.
Rendell
, and
B. O.
Roos
,
J. Phys. Chem.
94
,
5477
(
1990
).
14.
M.
Kállay
and
P. R.
Surján
,
J. Chem. Phys.
115
,
2945
(
2001
).
15.
D. W.
Small
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
13
,
19285
(
2011
).
16.
P.
Piecuch
,
S. A.
Kucharski
, and
R. J.
Bartlett
,
J. Chem. Phys.
110
,
6103
(
1999
).
17.
J.
Olsen
,
J. Chem. Phys.
113
,
7140
(
2000
).
18.
H.
Nakatsuji
and
K.
Yasuda
,
Phys. Rev. Lett.
76
,
1039
(
1996
).
19.
D. A.
Mazziotti
,
Chem. Rev.
112
,
244
(
2012
).
20.
B.
Verstichel
,
H.
van Aggelen
,
D.
Van Neck
,
P. W.
Ayers
, and
P.
Bultinck
,
Phys. Rev. A
80
,
032508
(
2009
).
21.
Y.
Kurashige
and
T.
Yanai
,
J. Chem. Phys.
135
,
094104
(
2011
).
22.
S.
Sharma
,
T.
Yanai
,
G. H.
Booth
,
C. J.
Umrigar
, and
G. K.-L.
Chan
,
J. Chem. Phys.
140
,
104112
(
2014
).
23.
T.
Yanai
,
Y.
Kurashige
,
E.
Neuscamman
, and
G. K.-L.
Chan
,
J. Chem. Phys.
132
,
024105
(
2010
).
24.
M.
Saitow
,
Y.
Kurashige
, and
T.
Yanai
,
J. Chem. Phys.
139
,
044118
(
2013
).
25.
D.
Ma
,
G.
Li Manni
, and
L.
Gagliardi
,
J. Chem. Phys.
135
,
044128
(
2011
).
26.
S.
Ten-no
,
Theor. Chem. Acc.
131
,
1070
(
2012
).
27.
M.
Torheyden
and
E. F.
Valeev
,
Phys. Chem. Chem. Phys.
10
,
3410
(
2008
).
28.
L.
Kong
,
F. A.
Bischoff
, and
E. F.
Valeev
,
Chem. Rev.
112
,
75
(
2011
).
29.
W.
Klopper
,
F. R.
Manby
,
S.
Ten-No
, and
E. F.
Valeev
,
Int. Rev. Phys. Chem.
25
,
427
(
2006
).
30.
T.
Yanai
and
T.
Shiozaki
,
J. Chem. Phys.
136
,
84107
(
2012
).
31.
T. B.
Adler
,
H.-J.
Werner
, and
F. R.
Manby
,
J. Chem. Phys.
130
,
054106
(
2009
).
32.
K.
Andersson
,
P. A.
Malmqvist
,
B. O.
Roos
,
A. J.
Sadlej
, and
K.
Wolinski
,
J. Phys. Chem.
94
,
5483
(
1990
).
33.
B. O.
Roos
,
Acc. Chem. Res.
32
,
137
(
1999
).
34.
C.
Angeli
,
R.
Cimiraglia
,
S.
Evangelisti
,
T.
Leininger
, and
J.-P.
Malrieu
,
J. Chem. Phys.
114
,
10252
(
2001
).
35.
I.
Schapiro
,
K.
Sivalingam
, and
F.
Neese
,
J. Chem. Theory Comput.
9
,
3567
(
2013
).
36.
F.
Liu
,
Y.
Kurashige
,
T.
Yanai
, and
K.
Morokuma
,
J. Chem. Theory Comput.
9
,
4462
(
2013
).
37.
O.
Sinanoǧlu
,
J. Chem. Phys.
34
,
1237
(
1961
).
38.
O.
Sinanoǧlu
,
Phys. Rev.
122
,
491
(
1961
).
39.
E.
Hylleraas
,
Z. Phys.
65
,
209
(
1930
).
40.
G. K.-L.
Chan
and
M.
Head-Gordon
,
J. Chem. Phys.
116
,
4462
(
2002
).
41.
K. G.
Dyall
,
J. Chem. Phys.
102
,
4909
(
1995
).
42.
43.
G. K.-L.
Chan
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
907
(
2012
).
45.
46.
S.
Sharma
and
G. K.-L.
Chan
,
J. Chem. Phys.
136
,
124121
(
2012
).
47.
G. K.-L.
Chan
,
J. Chem. Phys.
120
,
3172
(
2004
).
48.
T. H.
Dunning
,
J. Chem. Phys.
90
,
1007
(
1989
).
49.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
, and
M.
Schütz
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
242
(
2012
).
50.
Ö.
Legeza
and
G.
Fáth
,
Phys. Rev. B
53
,
14349
(
1996
).