Communication: A flexible multi-reference perturbation theory by minimizing the Hylleraas functional with matrix product states

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 H₀. Changing H₀ (as in CASPT2 versus NEVPT2) usually requires re-deriving and re-implementing non-trivial intermediates. In our formulation we can easily change H₀ 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 H[Ψ₁], Eq. (1), which is minimized with respect to the first order wavefunction |Ψ₁⟩,

Here H₀, E₀, and |Ψ₀⟩ are, respectively, the zeroth order Hamiltonian, energy, and wavefunction, V is the perturbation (H − H₀), and Q the projector onto the space orthogonal to |Ψ₀⟩. It can be easily verified that at the minimum, the wavefunction |Ψ₁⟩ satisfies the familiar equation

In MPS-PT2 we minimize the Hylleraas functional, while expressing |Ψ₁⟩ 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 (⟨Ψ₁|Ψ₂⟩) and transition elements of operators between two MPSs (⟨Ψ₁|O|Ψ₂⟩) can be evaluated in O(M³) 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.⁴⁰ 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

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

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,

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.

In the sweep algorithm, each tensor [Aⁿ] of the MPS wavefunction |Ψ₁⟩ 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

$Ail−1ilnl$ of |Ψ₁⟩, 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}$

$Ail−1ilnl$ (the quantity being solved for, shown with a red circle) with a 6-index tensor, formed by contracting all tensors in the operator H₀ − E₀ and |Ψ₁⟩ with

$A^{n_l}_{i_{l-1} i_l}$

$Ail−1ilnl$ removed. Similarly, the second and third terms involve the tensors of V, |Ψ₀⟩ and all tensors of |Ψ₁⟩ with

$A^{n_l}_{i_{l-1} i_l}$

$Ail−1ilnl$ removed. The resulting linear equation for

$A^{n_l}_{i_{l-1} i_l}$

$Ail−1ilnl$ is in O(dM²) unknowns. Solving it by conjugate gradient requires only O(M³) 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}$

$Ail−1ilnl$⁠. It is clear that the tensor so obtained depends on the values of all the other tensors [Aⁿ] in |Ψ₁⟩. Thus we sweep through all the orbitals of |Ψ₁⟩ 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 |Ψ₁⟩ and |Ψ₀⟩ (M₁ and M₀, respectively) and the number of active and external orbitals (k_a and k_v, respectively). For concreteness, we consider H₀ to be Dyall's Hamiltonian,⁴¹ as used in NEVPT2. All the contractions to be performed at each sweep iteration, together with their respective (leading order in M, k_a, k_v) 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 H₀, 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(ka2)$ auxiliary bond dimension, the same as the Dyall H₀. 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 k_v. 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 k_v (the number of virtual orbitals) to reflect the total number of steps in the sweep, assuming k_v + k_a ∼ O(k_v). 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 k_v is very large. For example, forming the MPO for V (through DMRG complementary operator techniques⁴⁵) costs

$O(k_v^2k_a^2M_0M_1)$

$O(kv2ka2M0M1)$/sweep, which becomes important when k_v > M.)

The above costs assume the Dyall H₀. A similar analysis applies to the CASPT2 H₀, 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 ⟨Ψ₀|Ψ₁⟩, but these only cost

$O(kaM02M1+kaM0M12)$ for the entire sweep, which is independent of k_v. The one-body Fock operator leads to a lower overall computational cost for the first term in the Hylleraas functional involving H₀. 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 H₀ in NEVPT and for the CASPT2 H₀ 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.³⁸ 

We have implemented the above MPS-PT sweep algorithm to evaluate the NEVPT2 and NEVPT3 energy in the BLOCK code.⁴⁶ 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 M₁, 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 H₀ and V is

$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 V_i. Each V_i produces a mutually orthogonal state when acting on |Ψ₀⟩ and thus the Hylleraas functional can be computed and minimized for each V_i 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 N₂ 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.⁴⁸ Table I shows the energies of the N₂ 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 μE_h with M₁ = 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 Molpro⁴⁹ and MPS-PT2 calculations with smaller M₁ = 100–300, are plotted against increasing N₂ bond-length in Figure 2. We see that the value of M₁ required to accurately describe the first order wavefunction |Ψ₁⟩ increases somewhat with increasing bond length. This can be understood by realising that the first-order wavefunction contains contributions from

$O(Naka2kv2)$ determinants (where N_a is the number of significant determinants in the active space). As N_a 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$

$M≪Naka2kv2$ and except for M = 100, the MPS-PT error remains small at all bond-lengths.

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

$|Ψ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.

Figure 4 shows the relationship between the size of the basis set and the value of M₁ required to obtain a second-order energy within 1 mE_h of the exact fully uncontracted MPS-PT2 energy. We see that, encouragingly, the M₁ 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, M₁ must asymptotically be independent of k_v.

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.