Numerical and Theoretical Aspects of the DMRG-TCC Method Exemplified by the Nitrogen Dimer

Numerical and Theoretical Aspects of the DMRG-TCC Method Exemplified by the Nitrogen Dimer

Fabian M. Faulstich f.m.faulstich@kjemi.uio.no    Andre Laestadius    Simen Kvaal Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway    Mihály Máté, Mihály András Csirik    Örs Legeza Wigner Research Center for Physics, H-1525, P.O. Box 49, Budapest, Hungary    Libor Veis, Andrej Antalik, Jiří Brabec    Jiří Pittner J. Heyrovský Institute of Physical Chemistry, Academy of Sciences of the Czech Republic, v.v.i., Dolejškova 3, 18223 Prague 8, Czech Republic    Reinhold Schneider Modeling, Simulation and Optimization in Science, Department of Mathematics, Technische Universität Berlin, Sekretariat MA 5-3, Straße des 17. Juni 136, 10623 Berlin,Germany
Abstract

In this article, we investigate the numerical and theoretical aspects of the coupled-cluster method tailored by matrix-product states. We investigate chemical properties of the used method, such as energy size extensivity and the equivalence of linked and unlinked formulation. The existing mathematical analysis is here elaborated in a quantum chemical framework. In particular, we highlight the use of a so-called CAS-ext gap describing the basis splitting between the complete active space and the external part. Moreover, the behavior of the energy error as a function of the optimal basis splitting is discussed. We show numerical investigations on the robustness with respect to the bond dimensions of the single orbital entropy and the mutual information, which are quantities that are used to choose the complete active space. Furthermore, we extend the mathematical analysis with a numerical study on the complete active space dependence of the error.

I. Introduction

The coupled-cluster (CC) theory has played a revolutionary role in establishing a new level of high accuracy in electronic structure calculations and quantum-chemical simulations. Despite the immense progress made in the field, computational schemes aiming at describing quasi-degenerate electronic structures of chemical systems are still unreliable. These multi-configuration systems, also called strongly correlated systems, form one of the most challenging computational problems in quantum chemistry. Since these systems appear in various research areas, a reliable computational scheme is of major interest for natural sciences. Recently, the computational advantages of the novel density matrix renormalization group tailored coupled-cluster (DMRG-TCC) method restricted to single (S) and double (D) excitations were demonstrated on large and statically correlated systems by L.V. et al. in Refs. veis2016coupled (); veis2018intricate (). Furthermore, computations showed that the use of the DMRG-TCCSD method is indispensable in order to determine the proper structure of the low lying energy spectrum in strongly correlated systems veis2018intricate (). In addition to these computational features, the DMRG-TCC approach is a promising candidate for a black-box quasi multi-reference scheme as big parts of the program already provide a black-box implementation up to a numerical threshold. This makes it accessible to a broad class of researchers from various fields of study. For an alternative multi-reference CC method that makes use of matrix product states and a modified DMRG algorithm, we refer to the linearized CCSD theory of Sharma and Alavi sharma2015multireference ().
The mathematical analysis of CC schemes is far from being complete, especially with regard to multi-reference methods, however, many important steps have already been taken. The list of fundamental and mathematical chemistry articles aiming to describe the existence and nature of solutions of CC equations is too long to be summarized here. We will limit our discussion to a short selection of publications addressing such fundamental issues of chemistry.
As a system of polynomial equations the CC equations can have real or, if the cluster operator is truncated, complex solutions. A standard tools to compute a solution of these non-linear equations is the Newton–Raphson and the quasi Newton method. However, these methods may diverge if the Jacobian or the approximated Jacobian become singular. This is in particular the case when strongly correlated systems are considered. These and other related aspects of the CC theory have been addressed by Živković and Monkhorst vzivkovic1977existence (); vzivkovic1978analytic () and Piecuch et al. piecuch1990coupled (). Significant advances in the understanding of the nature of multiple solutions of single-reference CC have been made by Kowalski and Jankowski kowalski1998towards () and by Piecuch and Kowalski piecuch2000computational (). An interesting attempt to address the existence of a cluster operator and cluster expansion in the open-shell case was done by Jeziorski and Paldus jeziorski1989valence (). The first advances within the rigorous realm of local functional analysis were performed by R.S. and Rohwedder, providing analyses of the closed-shell CC method for non multi-configuration systems schneider2009analysis (); rohwedder2013continuous (); rohwedder2013error (). Since then, the local mathematical analysis of CC schemes was extended by A.L. and S.K. analyzing the extended CC method laestadius2018analysis () and revisiting the bivariational approach lowdin1983stability (); arponen1983variational (), and F.M.F. et al. providing the first local mathematical analysis of a multi-reference scheme, namely the CC method tailored by tensor network states (TNS-TCC) faulstich2018analysis (). As this mathematical branch of CC theory is very young, channeling abstract mathematical results into the quantum-chemistry community is a highly interdisciplinary task merging both fields. A first attempt in this direction was done by A.L. and F.M.F. linking the physical assumption of a HOMO-LUMO gap to the somewhat abstract Gårding inequality and in that context presenting key aspects of Refs. schneider2009analysis (); rohwedder2013continuous (); rohwedder2013error (); laestadius2018analysis () from a more quantum chemical point of view laestadius2018coupled ().
With this article, we aim to bridge the mathematical results in Ref. faulstich2018analysis () of the TNS-TCC method into the quantum-chemistry community. Furthermore, we derive chemical properties of the TCC method and extend these results with a numerical study on the complete active space (CAS) dependence of the DMRG-TCCSD error.

II. The DMRG-TCC Method

As a post-Hartree–Fock method, the TCC approach was introduced by Kinoshita et al. in Ref. kinoshita2005coupled () as an alternative to the expensive and knotty conventional multi-reference methods. It divides the cluster operator into a complete active space (CAS) part, denoted , and an external (ext) part , i.e., the wave function is parametrized as

Separating the cluster operator into several parts goes back to Piecuch et al. piecuch1993state (); piecuch1994state (). In this formulation the linked CC equations are given by

(1)

Computing first and keeping it fixed for the dynamical correction via the CCSD method restricts the above equations to not in the CAS, i.e., for all in the CAS (we say that is in the -orthogonal complement of the CAS). We emphasize that this includes mixed states, e.g., where is an element of the CAS but is not. We consider a CAS created by spin-orbitals , forming a subspace of the full configuration interaction (FCI) space created by the entire set of spin-orbitals . We here assume the spin-orbitals to be eigenfunctions of the system’s Fock operator. Note that the following analysis can be applied for any single-particle operator fulfilling the properties used in Ref. faulstich2018analysis () – not only the Fock operator.

Based on the single reference approach, the TCC method needs a large CAS to cover most of the static correlations. As the size of the CAS scales exponentially, an efficient approximation scheme for strongly correlated systems is indispensable for the TCC method to have practical significance. One of the most efficient schemes for static correlation is the DMRG method chan2011density (). Going back to the physicists S. R. White and R. L. Martin white1999ab (), it was introduced to quantum chemistry as an alternative to the CI or CC approach. However, the major disadvantage of the DMRG is that in order to compute dynamical correlation high bond dimensions (tensor ranks) may be necessary, making the DMRG a potentially costly method veis2018intricate (); chan2011density (). Nevertheless, as a TNS-TCC method, the DMRG-TCC approach is an efficient method since the CAS is supposed to cover the statically correlated orbital of the system. This avoids the DMRG method’s weak point and allows to employ a much larger CAS compared to the traditional CI-TCC method.
A notable advantage of the TCC approach over the conventional MRCC methods is that excitation operators commute. This is due to the fact that the Hartee–Fock method yields a single reference solution , which implies that separating the cluster operator corresponds to a partition of excitation operators. Hence, and commute. This makes the DMRG-TCC method’s analysis much more accessible than conventional MRCC methods and therewith facilitates establishing sound mathematical results faulstich2018analysis (). Note also that the dynamical correlation via the CCSD method distinguishes the DMRG-TCC method from standard CAS-SCF methods, which suffer from an imbalance in the considered correlation making these methods inapplicable for chemical phenomena like surface hopping. Consequently, due to its more balanced correlation the DMRG-TCC approach is much better suited for these phenomena, which widens the horizon of possible applications. We remark, however, that the computationally most demanding step of the DMRG-TCC calculation is the DMRG part, and its cost increases rapidly with . Thus the above application requires further investigations. Alternative to the dynamical correction via the CC approach, the DMRG-MRCI method in Ref. saitow2013multireference () utilizes an internally contracted CI algorithm different from a conventional CI calculation.

III. Chemical Properties of The DMRG-TCC

It is desired that quantum-chemical computations possess certain features representing the system’s chemical and physical behavior. Despite their similarity, the CC and TCC method have essentially different properties, which are here elaborated. A basic property of the CC method is the equivalence of linked and unlinked CC equations. We point out that this equivalence is in general not true for the DMRG-TCCSD scheme. This is a consequence of the CAS ansatz since it yields mixed states, i.e., two particle excitations with one excitation into the CAS. The respective overlap integrals in the unlinked CC equations will then not vanish unless the single excitation amplitudes are equal to zero. Generalizing this result for rank complete truncations of order we find that all excitation amplitudes need to be zero but for the -th one. This is somewhat surprising as the equivalence of linked and unlinked CC equations holds for rank complete truncations of the single-reference CC method.
For the sake of simplicity we show this results for the DMRG-TCCSD method. The general case can be proven in similar a fashion. We define the matrix representation with elements for . Note that, as increases the excitation rank, is an atomic lower triangular matrix and therefore not singular. Assuming that the linked CC equations hold, the non-singularity of yields

As the full projection manifold is complete under de-excitation, we obtain that

(2)

Note that the first two terms on the r.h.s. in Eq. (2) describe the unlinked CC equations. We analyze the last term on the r.h.s. in Eq. (2) by expanding the inner products, i.e.,

The first term in this expansion vanishes due to orthogonality. The same holds true for all terms where enters to the power of two or higher. However, as the external space contains mixed states, we find that is not necessarily zero, namely, for and with and . This proves the claim.
Subsequently, we elaborate the size extensivity of the DMRG-TCCSD. Let two DMRG-TCCSD wave functions for the individual subsystems and be

The corresponding energies are given by

and the amplitudes fulfill

in terms of the effective, similarity-transformed Hamiltonians

The Hamiltonian of the compound system of the non-interacting subsystems can be written as . Since the TCC approach corresponds to a partitioning of the cluster amplitude we note that for

With , the energy of the compound systems can be written as

It remains to show that

solves the Schrödinger equation, i.e., for all holds . Splitting the argument into three case, we note that

where . This proves the energy size extensivity for the untruncated TCC method. From this we conclude the energy size extensivity for the DMRG-TCCSD scheme, because the truncation only affects the product states and these are zero in the above projection.

Looking at TCC energy expression we observe that due to the Slater–Condon rules, these equations are independent of CAS excitations higher than order three, i.e., amplitudes of for . More precisely, due to the fact that in the TCCSD case external space amplitudes can at most contain one virtual orbital in the CAS, the TCCSD amplitude expressions become independent of , i.e.,

where the primed variables describe orbitals in the CAS, the non-primed variable describes an orbital in the external part and are occupied orbitals. Note, this does not imply that we can restrict the CAS computation to a manifold characterizing excitations with rank less or equal to three as for strongly correlated systems these can still be relevant. However, it reduces the number of terms entering the DMRG-TCCSD energy computations significantly.
In contrast to the originally introduced CI-TCCSD method taking only for into account kinoshita2005coupled (), the additional consideration of corresponds to an exact treatment of the CAS contributions to the energy. We emphasize that the additional terms do not change the methods complexity. This is due to the fact that including the CAS triple excitation amplitudes will not exceed the dominating complexities of the CCSD myhre2016multilevel () nor of the DMRG. However, the extraction of the CI-triples from the DMRG wave function is costly and a corresponding efficiency investigation is left for future work.

IV. Analysis of the DMRG-TCC

In the sequel we discuss and elaborate mathematical properties of the TCC approach and their influence on the DMRG-TCC method. The presentation here is held brief and the interested reader is referred to Ref. faulstich2018analysis () and the references therein for further mathematical details.

A. The Complete Active Space Choice

As pointed out in the previous section, the TCC method relies on a well-chosen CAS, i.e., a large enough CAS that covers the system’s static correlation. Consequently, we require a quantitative measurement for the quality of the CAS, which presents the first obstacle for creating a non-empirical model since the chemical concept of correlation is not well-defined lyakh2011multireference (). In the DMRG-TCC method, we use a quantum information theory approach to classify the spin-orbital correlation. This classification is based on the mutual-information

This two particle entropy is defined via the von Neumann entropy of the reduced density operators szalay2017correlation (). Note that the mutual-information describes two-particle correlations, for a more general connection between multi-particle correlations and -correlations we refer the reader to Szalay et. al szalay2017correlation (). We emphasize that in practice this is a basis dependent quantity, which is in agreement with the chemical definition of correlation concepts lyakh2011multireference (). We identify pairs of spin-orbitals contributing to a high mutual information value as strongly correlated, the pairs contributing to the plateau region as non-dynamically correlated and the pairs contributing to the mutual information tail as dynamically correlated (see Fig. 3). The mutual-information profile can be well approximated from a prior DMRG computation on the full system. Due to the size of the full system we only compute a DMRG solution of low bond dimension (also called tensor rank). These low-accuracy calculations, however, already provide a good qualitative entropy profile, i.e., the shapes of profiles obtained for low bond dimension, , agree well with the the ones obtained in the FCI limit. Here, we refer to Fig. 2 and Fig. 3 showing the single orbital entropy and mutual information profiles, respectively, for various values and for three different geometries of the N molecule. The orbitals with large entropies can be identified from the low- calculations providing a black-box tool to form the CAS space including the strongly correlated orbitals.

A central observation is that for (i.e. ), the DMRG-TCCSD becomes the CCSD method and for (i.e. ), it is the DMRG method. We recall that the CCSD method can not resolve static correlation and the DMRG method needs high tensor ranks for dynamically correlated systems. This suggests that the error obtains a minimum for some with , i.e., there exists an optimal choice of determining the basis splitting and therewith the choice of the CAS. Note that this feature becomes important for large systems since high bond dimensions become simply impossible to compute with available methods.

B. Local Analysis of the DMRG-TCC Method

The CC method can be formulated as non-linear Galerkin scheme schneider2009analysis (), which is a well-established framework in numerical analysis to convert the continuous Schrödinger equation to a discrete problem. For the DMRG-TCC method a first local analysis was performed in Ref. faulstich2018analysis (). There, a quantitative error estimate with respect to the basis truncation was established. F.M.F. et al. showed under certain assumptions (Assumption A and B in the sequel) that the DMRG-TCC method possesses a locally unique and quasi-optimal solution (cf. Sec. 4.1 in Ref. faulstich2018analysis ()). In case of the DMRG-TCC method the latter means: For a fixed basis set the CC solution tailored by a DMRG solution on a fixed CAS is up to a multiplicative constant the best possible solution in the approximation space defined by the basis set. In other words, the CC method provides the best possible dynamical correction for a given CAS solution such as a DMRG solution.
Note that local uniqueness ensures that for a fixed basis set, the computed DMRG-TCC solution is unique in a neighborhood around the exact solution. We emphasize that this result is derived under the assumption that the CAS solution is fixed. Consequently, for different CAS solutions we obtain in general different TCC solutions, i.e., different cluster amplitudes.
Subsequently, parts of the results in Ref. faulstich2018analysis () are explained in a setting adapted to the theoretical chemistry perspective. The TCC function is given by for not in the CAS. Note that we use the convention where small letters correspond to cluster amplitudes, whereas capital letters describe cluster operators. The corresponding TCC energy expression is given by

Consequently, the linked TCC equations (1) then become

Within this framework the locally unique and quasi-optimal solutions of the TCC method were obtained under two assumptions (see Assumption A and B in Ref. faulstich2018analysis ()).
First, Assumption A requires that the Fock operator is bounded and satisfies a so-called Gårding inequality. For a more detailed description of these properties in this context we refer the reader to Ref. laestadius2018coupled (). Second, it is assumed that there exists a CAS-ext gap in the spectrum of the Fock operator, i.e., there is a gap between the -th and the -st orbital energies. Note, that spectral gap assumptions (cf. HOMO-LUMO gap) are standard in the analysis of dynamically correlated systems. To be applicable to multi-configuration systems, the analysis in Ref. faulstich2018analysis () rests on a CAS chosen such that the -th and the -st spin orbital are non-degenerate. Intuitively, this gap assumption means that the CAS captures the static correlation of the system.
Assumption B is concerned with the fluctuation operator . This operator describes the degeneracy of the system since it is the difference of the Hamiltonian and a single particle operator. Using the similarity transformed and fixing the CAS amplitudes , the map

is assumed to have a small enough Lipschitz-continuity constant (see Eq. (20) in Ref. faulstich2018analysis ()). The physical interpretation of this Lipschitz condition is at the moment unclear.

C. Error Estimates for the DMRG-TCC Method

A major difference between the CI and CC method is that the CC formalism is not variational. Hence, it is not evident that the CC energy error decays quadratically with respect to the error of the wave function or cluster amplitudes. Note that the TCC approach represents merely a partition of the cluster operator, however, its error analysis is more delicate than the traditional CC method’s analysis. The TCC-energy error is measured as difference to the FCI energy. Let describe the FCI solution on the whole space, i.e., . Using the exponential parametrization and the above introduced separation of the cluster operator, we have

(3)

An important observation is that the TCC approach ignores the coupling from the external space into the CAS. It follows that the FCI solution on the CAS is an approximation to the projection of onto the CAS

where is the -orthogonal projection onto the CAS. For a reasonably sized CAS the FCI solution is rarely computationally accessible and we introduce the DMRG solution on the CAS as an approximation of ,

Tailoring the CC method with these different CAS solutions leads in general to different TCC solutions. In the case of , the TCC method yields the best possible solution with respect to the chosen CAS, i.e., . This solution is in general different from fulfilling and its truncated version satisfying , where denotes the -orthogonal projection onto the corresponding Galerkin space. In the context of the DMRG-TCC theory, the Galerkin space represents a truncation in the excitation rank of the cluster operator, e.g., DMRG-TCCD, DMRG-TCCSD, etc.
For the following argument, suppose that an appropriate CAS has been fixed. The total DMRG-TCC energy error can be estimated as faulstich2018analysis ()

(4)

where each term of the r.h.s. in Eq. (4) is now discussed. As a technical remark, the norms on either the Hilbert space of cluster amplitudes or wave functions are here simply denoted . These norms are not just the - or -norm, respectively, but also measure the kinetic energy. It should be clear from context which Hilbert space is in question and we refer to Ref. rohwedder2013continuous () for formal definitions. The first term is defined as

which describes the truncation error of the CCSD method tailored by . We emphasize that the dynamical corrections via the CCSD and the untruncated CC method are here tailored by the same CAS solution. Hence, the energy error corresponds to a single reference CC energy error, which suggests an analysis similar to Refs. schneider2009analysis (); rohwedder2013error (). Indeed, the Aubin–Nitsche duality method aubin1967behavior (); nitsche1968kriterium (); ruchovec1969study () yields a quadratic a priori error estimate in (and in terms of the Lagrange mulitpliers, see Theorem 29 in Ref. faulstich2018analysis ()).
Second, we discuss the term

Here, different CAS solutions with fixed external solutions are used to compute the energies. This suggests that is connected with the error

(5)

describing the approximation error of the DMRG solution on the CAS (see Lemma 27 in Ref. faulstich2018analysis ()). Indeed,

(6)

with , for , where are the orbital energies. The are the (translated) Fock energies, more precisely, , with . Note that the wave function is in general not an eigenfunction of , however, it is an eigenfunction of the projected Hamiltonian . Eq. (5) involves the exponential parametrization. This can be estimated by the energy error of the DMRG wave function, denoted , namely,

(7)

In Sec. V the energy error of the DMRG wave function is controlled by the threshold value , i.e., . Hence, for well chosen CAS the difference is sufficiently small such that holds. This again shows the importance of a well-chosen CAS. Furthermore, the last term in Eq. (6) can be eliminated via orbital rotations, as it is a sum of single excitation amplitudes.
Finally, we consider

(8)

Since is a stationary point of we have . A calculation involving Taylor expanding around (see Lemma 26 in Ref. faulstich2018analysis ()) yields

(9)

Note that the above error is caused by the assumed basis splitting, namely, the correlation from the external part into the CAS is ignored. Therefore, the best possible solution for a given basis splitting differs in general from the FCI solution .
Combining now the three quadratic bounds gives an overall quadratic a priori energy error estimate for the DMRG-TCC method. The interested reader is referred to Ref. faulstich2018analysis () for a more detailed treatment of the above analysis.

D. On the -dependency of the Error Estimates

The error estimate outlined above is for a fixed CAS, i.e., a particular basis splitting, and bounds the energy error in terms of truncated amplitudes. Because the TCC solution depends strongly on the choice of the CAS, it is motivated to further investigate the -dependence of the error . However, the above derived error bound has a highly complicated dependence since not only the amplitudes but also the implicit constants (in ) and norms depend on . Therefore, the analysis in Ref. faulstich2018analysis () is not directly applicable to take the full dependence into account.
In the limit where we obtain that since the TCC method is numerically stable, i.e., a small perturbation in corresponds to a small perturbation in the solution . Furthermore, if we assume that , which is reasonable for the equilibrium bond length of N, the error can be bound as

(10)

Here the subscript on and highlights the -dependence. We remark that we here used the less accurate structure on the amplitude space compared to the structure in Eq. (9). This yields -independent vectors

as well as an -independent norm. The -depenence of will be investigated numerically in more detail in Sec. V B 5.

V. The Splitting Error for N

Including the dependence in the above performed error analysis explicitly is a highly non-trivial task involving many mathematical obstacles and is part of our current research. Therefore, we here extend the mathematical results from Sec. IV with a numerical investigation on this dependence. Our study is presented for the N molecule using the cc-pVDZ basis, which is a common basis for benchmark computations developed by Dunning and coworkers dunning1989gaussian (). We investigate three different geometries of the Nitrogen dimer by stretching the molecule, thus the performance of DMRG-TCCSD method is assessed against DMRG and single reference CC methods for bond lengths , and . In the equilibrium geometry the system is weakly correlated implying that single reference CC methods yield reliable results. For increasing bond length the system shows multi-reference character, i.e., static correlations become more dominant. For this results in the variational breakdown of single reference CC methods kowalski2000renormalized (). This breakdown can be overcome with the DMRG-TCCSD method once a large and well chosen CAS is formed, we therefore refer to the DMRG-TCCSD method as numerically stable with respect to the bond length along the potential energy surface (PES).

As mentioned before, the DMRG method is in general less efficient to recover dynamic correlations since it requires large computational resources. However, due to the specific CAS choice the computational resource for the DMRG part of the TCC scheme is expected to be significantly lower than a pure DMRG calculation for the same level of accuracy.

A. Computational Details

In practice, a routine application of the TCC method to strongly correlated molecular systems, i.e., to multi-reference problems, became possible only recently since it requires a very accurate solution in a large CAS including all static correlations. Tensor network state methods fulfill such a high accuracy criterion, but the efficiency of the TNS-TCCSD method strongly depends on various parameters of the involved algorithms. Some of these are defined rigorously while others are more heuristic from the mathematical point of view. In this section we present the optimization steps for the most important parameters of the DMRG-TCCSD method and outline how the numerical error study in Sec. V B is performed.

As elaborated in Sec. II and IV A, the CAS choice is essential for the computational success of TNS-TCC methods. In addition, the error of the TNS method used to approximate the CAS part depends on various approximations. These include the proper choice of a finite dimensional basis to describe the chemical compound, the tensor network structure, and the mapping of the molecular orbitals onto the given network Szalay-2015 (). Fortunately, all these can be optimized by utilizing concepts of quantum information theory, introduced in Sec. IV A (see also the included references). In the following, we restrict the numerical study to the DMRG-TCCSD method but the results presented here should also hold for other TNS approaches Murg-2010a (); Nakatani-2013 (); Murg-2015 (); Szalay-2015 (); Gunst-2018 ().
In the DMRG-TCCSD case the tensor network topology in the CAS corresponds to a single branched tensor tree, i.e., a one dimensional topology. Thus permutations of orbitals along such an artificial chain effect the convergence for a given CAS choice Chan-2002a (); Legeza-2003b (). This orbital-ordering optimization can be carried out based on spectral graph theory Barcza-2011 (); Fertitta-2014 () by minimizing the entanglement distance Rissler-2006 (), defined as . In order to speed up the convergence of the DMRG procedure the configuration interaction based dynamically extended active space (CI-DEAS) method is applied Legeza-2003b (); Szalay-2015 (). In the course of these optimization steps, the single orbital entropy () and the two-orbital mutual information are calculated iteratively until convergence is reached. The size of the active space is systematically increased by including orbitals with the largest single site entropy values, which at the same time correspond to orbitals contributing to the largest matrix elements of the mutual information. Thus, the decreasingly ordered values of define the so-called CAS vector, CAS, which provides a guide in what order to extend the CAS by including additional orbitals. The bond dimensions (tensor rank) in the DMRG method can be kept fixed or adapted dynamically (Dynamic Block State Selection (DBSS) approach) in order to fulfill an a priori defined error margin Legeza-2003a (); Legeza-2004b (). Accurate extrapolation to the truncation free limit is possible as a function of the truncation error Legeza-1996 (); Legeza-2003a ().
In our DMRG implementation QCDMRG-Budapest () we use a spatial orbital basis, i.e., the local tensor space of a single orbital is dimensional. In this representation an orbital can be empty, singly occupied with either a spin up or spin down electron, or doubly occupied with opposite spins. Note, in contrast to Sec. IV we need spatial orbitals to describe an -electron wave function and similar changes apply to the size of the basis set so that we use from here on. The single orbital entropy therefore varies between 0 and , while the two-orbital mutual information varies between 0 and .
Next we provide a short description how to perform DMRG-TCCSD calculations in practice. Note that we leave the discussion on the optimal choice of for the following sections.
First the CAS is formed from the full orbital space by setting . DMRG calculations are performed iteratively with fixed low bond dimension (or with a large error margin) in order to determine the optimal ordering and the CAS-vector as described above. Thus, the corresponding single-orbital entropy and mutual information are also calculated. These calculations already provide a good qualitative description of the entropy profiles with respect to the exact solution, i.e., strongly correlated orbitals can be identified.
Using a given we form the CAS from the Hartree–Fock orbitals and the first virtual orbitals from the CAS, i.e., orbitals with the largest single orbital entropy values. We emphasize that these orbitals contribute to the largest matrix elements in . We carry out the orbital ordering optimization on the given CAS and perform a large-scale DMRG calculation with a low error threshold margin in order to get an accurate approximation of the . Note, that the DMRG method yields a normalized wave function, i.e., the overlap with the reference determinant is not necessarily equal to one.
Using the matrix product state representation of obtained by the DMRG method we determine the zero reference overlap, single and double CI coefficients of the full tensor representation of the wave function. Next, these are used to calculate the and amplitudes, which form the input of the forthcoming CCSD calculation.
Next, the cluster amplitudes for the external part, i.e., and , are calculated in the course of the DMRG-TCCSD scheme.
As we discus in the next section, finding the optimal CAS, i.e., -splitting, is a highly non-trivial problem, and at the present stage we can only present a solution that is considered as a heuristic approach in terms of rigorous mathematics. In practice, we repeat steps 2-4 for a large DMRG-truncation error as a function of , thus we find local energy minima (see Fig. 4) using a relatively cheap DMRG-TCCSD scheme. Around such a local minimum we perform more accurate DMRG-TCCSD calculations by lowering the DMRG-truncation error in order to refine the optimal . We also monitor the maximum number of DMRG block states required to reach the a priori defined DMRG-error margin as a function of . Since it can happen that several values lead to low error DMRG-TCCSD energies, while the computational effort increases significantly with increasing we select the optimal that leads to low DMRG-TCCSD energy but also minimizes the required DMRG block states. Using the optimal value we perform large-scale DMRG-TCCSD calculation using a relatively tight error bound for the DMRG-truncation error.
We close this section with a brief summary of the numerically accessible error terms and relate them to equations presented in Sec. IV. Note that the error analysis in Sec. IV is presented for a given , thus here the dependence is also omitted.
For a given split, the accuracy of depends on the DMRG truncation error, . As has been shown in Refs. Legeza-1996 (); Legeza-2003a () the relative error, is a liner function of on a logarithmic scale. Therefore, extrapolation to the FCI limit can be carried out as a function of . In addition, the error term appearing in Eq. (7) can be controlled.
Note that terms appearing in Eq. (6) and Eq. (7) include FCI solutions of the considered system. However, in special cases these can be well approximated as follows: Considering a small enough system that is dynamically correlated, like the Nitrogen dimer near the equilibrium geometry with the here chosen basis set. The CI-coefficients are then extractable from the matrix product state representation of a wave function, e.g., or . Note that calculating all CI-coefficients scales exponentially with the size of the CAS. However, since the system is dynamically correlated zeroth order, single and double excitation coefficients are sufficient. Hence the error terms and in Eq. (6) and Eq. (7), respectively, can be well approximated. We remark that this exponential scaling with the CAS size also effects the computational costs of the CAS CI-triples, which are needed for an exact treatment of the TCCSD energy equation. However, investigations of the influence of the CAS CI-triples on the computed energies are left for future work.

B. Results and Discussion

In this section, we investigate the overall error dependence of DMRG-TCCSD as a function of and as a function of the DMRG-truncation error . For our numerical error study we perform steps 1-4 discussed in Sec. V A for each . For each geometry and we also carry out very high accuracy DMRG calculations on the full orbital space, i.e., by setting the truncation error to and . This data is used as a reference for the FCI solution.

1. An Entropy Study on the Full Orbital Space

We start our investigation by showing DMRG results for the full orbital space, i.e., the CAS is formed from orbitals, and for various fixed values and for . In the latter case the maximum bond dimension was set to . In Fig. 1 (a) we show the relative error of the ground-state energy as a function of the DMRG-truncation error on a logarithmic scale. For the FCI energy, , the CCSDTQPH reference energy is used given in Ref. Chan-2004b (). It is visible that the relative error is a linear function of the truncation error on a logarithmic scale, thus extrapolation to the truncation free solution can be carried out according to Refs. Legeza-1996 (); Legeza-2003a ().

Figure 1: (a) Relative error of the ground-state energy as a function of the DMRG-truncation error on a logarithmic scale obtained for the full orbital space () with . (b) Maximum number of block states as a function of for the a priori defined truncation error with (blue), (green), and (red).

In Figs. 2 and 3 we present the sorted values of the single orbital entropy and of the mutual information obtained for fixed and with for the three geometries. As can be seen in the figures the entropy profiles obtained with low-rank DMRG calculations already resemble the main characteristics of the exact profile (). Therefore, orbitals with large single orbital entropies, also contributing to large matrix elements of , can easily be identified from a low-rank computation. The ordered orbital indices define the CAS-vector, and the CAS for the DMRG-TCCSD can be formed accordingly as discussed in Sec. V A.

Figure 2: Single orbital entropy for (blue), (green), (red) obtained for the full orbital space () with DMRG for fixed and for .

Figure 3: Mutual information for (blue), (green), (red) obtained for the full orbital space () with DMRG for fixed and for .

Taking a look at Fig. 2 it becomes apparent that shifts upwards for increasing indicating the higher contribution of static correlations for the stretched geometries. Similarly the first 50-100 matrix elements of also take larger values for larger while the exponential tail, corresponding to dynamic correlations, is less effected. The gap between large and small values of the orbital entropies gets larger and its position shifts rightward for larger . Thus, for the stretched geometries more orbitals must be included in the CAS during the TCC scheme in order to determine the static correlations accurately. We remark here, that the orbitals contributing to the high values of the single orbital entropy and mutual information matrix elements change for the different geometries according to chemical bond forming and breaking processes Boguslawski-2013 ().

2. Numerical Investigation of the Error’s -dependence

In order to obtain in the FCI limit, we perform high-accuracy DMRG calculations with . The CAS was formed by including all Hartee–Fock orbitals and its size was increased systematically by including orbitals with the largest entropies according to the CAS vector. Orbitals with degenerate single orbital entropies, due to symmetry considerations, are added to the CAS at the same time. Thus there are some missing points in the following figures. For each restricted CAS we carry our the usual optimization steps of a DMRG scheme as discussed in Sec. V A, with low bond dimension followed by a high-accuracy calculation with using eight sweeps Szalay-2015 (). Our DMRG ground-state energies for together with the CCSD (corresponding to a DMRG-TCCSD calculation where ), and CCSDTQ reference energies, are shown in Fig. 4 near the equilibrium bond length, . The single-reference coupled cluster calculations were performed in NWChem VALIEV20101477 (), we employed the cc-pVDZ basis set in the spherical representation. For the CCSDTQPH energy was taken as a reference for the FCI energy Chan-2004b ().

Figure 4: Ground-state energy of the N molecule near the equilibrium geometry, , obtained with DMRG-TCCSD for and for various DMRG truncation errors . The CCSD, CCSDT and CCSDTQ reference energies are shown by dotted, dashed and dashed-dotted lines, respectively. The DMRG energy with on the full space, i.e., , is taken as a reference for the FCI energy. For the CAS was additionally formed by taking orbitals according to increasing values of the single-orbital entropy values, i.e., inverse to the other CAS extensions. This is labeled by CAS.

The DMRG energy starts from the Hartree–Fock energy for and decreases monotonically with increasing until the full orbital solution with is reached. It is remarkable, however, that the DMRG-TCCSD energy is significantly below the CCSD energy for all CAS choices, even for a very small . The error, however, shows an irregular behavior taking small values for several different -s. This is due to the fact that the DMRG-TCCSD approach suffers from a methodological error, i.e., certain fraction of the correlations are lost, since the CAS is frozen in the CCSD correction. This supports the hypothesis of a -dependent constant as discussed in Sec. IV D. Therefore, whether orbital is part of the CAS or external part provides a different methodological error. This is clearly seen as the error increases between and although the CAS covers more of the system’s static correlation with increasing . This is investigated in more detail in Sec. V B 4.
Since several -splits lead to small DMRG-TCCSD errors, the optimal value from the computational point of view, is determined not only by the error minimum but also by the minimal computational time, i.e., we need to take the computational requirements of the DMRG into account. Note that the size of the DMRG block states contributes significantly to the computational cost of the DMRG calculation. The connection of the block size to the CAS choice is shown in Fig. 1 (b), where the maximal number of DMRG block states is depicted as a function of for the a priori defined truncation error margin . Note that increases rapidly for . The optimal CAS is therefore chosen such that the DMRG block states are not too large and the DMRG-TCCSD provides a low error, i.e., is a local minimum in the residual with respect to .
It is important to note that based on Fig. 4 the DMRG-TCCSD energy got very close to, or even dropped below, the CCSDT energy for several values. Since close to the equilibrium geometry the wave function is dominated by a single reference character, it is expected that DMRG-TCCSD leads to even more robust improvements for the stretched geometries, i.e., when the multi-reference character of the wave function is more pronounced. Our results for the stretched geometries, and , are shown in Figs. 23 5 and 6. As mentioned in Sec. V B 1, for larger values static correlations gain importance signaled by the increase in the single orbital entropy in Fig. 2. Thus the multi-reference character of the wave function becomes apparent through the entropy profiles. According to Fig. 5 the DMRG-TCCSD energy for all values is again below the CCSD computation and for it is even below the CCSDT reference energy. For the CC computation fluctuates with increasing excitation ranks and CCSDT is even far below the FCI reference energy, revealing the variational breakdown of the single-reference CC method for multi-reference problems. In contrast to this, the DMRG-TCCSD energy is again below the CCSD energy for all , but above the the CCSDT energy. The error furthermore shows a local minimum around .

Figure 5: Ground-state energy of the N molecule with bond length , obtained with DMRG-TCCSD for and for various DMRG truncation errors . The CCSD, CCSDT and CCSDTQ reference energies are shown by dotted, dashed and dashed-dotted lines, respectively. The DMRG energy with on the full space, i.e., , is taken as a reference for the FCI energy. For the CAS was additionally formed by taking orbitals according to increasing values of the single-orbital entropy, i.e., inverse to the other CAS extensions. This is labeled by CAS.

Figure 6: Ground-state energy of the N molecule with bond length , obtained with DMRG-TCCSD for and for various DMRG truncation errors . The CCSD, CCSDT and CCSDTQ reference energies are shown by dotted, dashed and dashed-dotted lines, respectively. The DMRG energy with on the full space, i.e., , is taken as a reference for the FCI energy. For the CAS was additionally formed by taking orbitals according to increasing values of the single-orbital entropy, i.e., inverse to the other CAS extensions. This is labeled by CAS.

For the stretched geometries static correlations are more pronounced, there are more orbitals with large entropies, thus the maximum number of DMRG block states increases more rapidly with compared to the situation near the equilibrium geometry (see Fig. 1 (b)). Thus obtaining an error margin within one milli-Hartree for leads to a significant save in computational time and resources.

3. Effect of on the DMRG-TCCSD

In practice, we do not intend to carry out DMRG calculations in the FCI limit, thus usually a larger truncation error is used. Therefore, we have repeated our calculations for larger truncation errors in the range of and . Our results are shown in Figs. 4, 5, and 6. For small the DMRG solution basically provides the Full-CI limit since the a priori set minimum number of block states already leads to a very low truncation error. Therefore, the error of the DMRG-TCCSD is dominated by the methodological error. For the effect of the DMRG truncation error becomes visible and for large the overall error is basically determined by the DMRG solution. For larger between and the DMRG-TCCSD error shows a minimum with respect to . This is exactly the expected trend, since the CCSD method fails to capture static correlation while DMRG requires large bond dimension to recover dynamic correlations, i.e., a low truncation error threshold. In addition, the error minima for different truncation error thresholds happen to be around the same values. This has an important practical consequence: the optimal -split can be determined by performing cheap DMRG-TCCSD calculations using large DMRG truncation error threshold as a function of .
The figures furthermore indicate that has a high peak for . This can be explained by the splitting of the FCI space since this yields that the correlation from external orbitals with CAS orbitals is ignored. Thus we also performed calculations for using a CAS formed by taking orbitals according to increasing values of the single orbital entropy values in order to demonstrate the importance of the CAS extension. The corresponding error profile as a function of near the equilibrium geometry is shown in Fig. 4 labeled by CAS. As expected, the improvement of DMRG-TCCSD is marginal compared to CCSD up to a very large split since differs only marginally from .

4. Numerical Investigation on CAS-ext correlations

Taking another look at Fig. 2, we can confirm that already for small values the most important orbitals, i.e., those with the largest entropies, are included in the CAS.

Figure 7: (a) Sorted values of the mutual information obtained by DMRG() for on a semi-logarithmic scale for N at . (b) Sorted 40 largest matrix elements of the mutual information obtained by DMRG() for on a lin-lin scale for N at

In Fig. 7 the sorted values of the mutual information obtained by DMRG() for is shown on a semi-logarithmic scale. It is apparent from the figure that the largest values of change only slightly with increasing , thus static correlations are basically included for all restricted CAS. The exponential tail of corresponding to dynamic correlations, however, becomes more visible only for larger values. We conclude, for a given split the DMRG method computes the static correlations efficiently and the missing tail of the mutual information with respect to the full orbital space () calculation is captured by the TCC scheme.
Correlations between the CAS and external parts can also be simulated by a DMRG calculation on the full orbital space using an orbital ordering according to the CAS-vector. In this case, the DMRG left block can be considered as the CAS and the right block as the external part. For a pure target state, for example, the ground state, the correlations between the CAS and external part is measured by the block entropy, as a function of . Here is formed by a partial trace on the external part of . The block entropy is shown in Fig. 8 (a). The block entropy decays monotonically for , i.e, the correlations between the CAS and the external part vanish with increasing . In contrast to this, when an ordering according to CAS is used the correlation between CAS and external part remains always strong, i.e., some of the highly correlated orbitals are distributed among the CAS and the external part. Nevertheless, both curves are smooth and they cannot explain the error profile shown in Fig. 4.

5. Numerical Values for the Amplitude Error Analysis

Since correlation analysis based on the entropy functions cannot reveal the error profile shown in Fig. 4, here we reinvestagte the error behavior as a function of but in terms of the CC amlitudes. Therefore, we also present a more detailed description of Eq. (10) in Sec. IV which includes the following terms:

(11)

Here the valid index-pairs are , with , and . The excitation rank is given by where stands for singles, for doubles, and so on. The -s are the labels of excitation operators , and . The corresponding amplitudes are given as . For invalid index-pairs, i.e., index-pairs that are out of range, the amplitudes are always zero. The various amplitudes appering in Eq. (11) are calculated according to the following rules:
1) The : amplitudes in the CAS() obtained by DMRG() solution (represented by CI coefficients ) for CAS(),

(12)

where and .
2) The : amplitudes not in the CAS() obtained from the DMRG() solution (represented by CI coefficients ) for CAS() projected onto CAS(), i.e., the complement (with respect to valid index-pairs) of .
3) The amplitudes in the CAS() are obtained by the DMRG() solution (represented by CI coefficients ) for CAS(). The amplitudes are the same as Eq. (12), but with , where and .
4) The : amplitudes not in the CAS() obtained by TCCSD, i.e., the complement (with respect to valid index-pairs) of