Coupled cluster downfolding methods: The effect of double commutator terms on the accuracy of ground-state energies

Over the last few decades, the coupled cluster (CC) theory^1–9 has evolved into one of the most accurate and dominant theories to describe various many-body systems across spatial scales and therefore addressing fundamental problems in quantum chemistry,^10–17 material sciences,^18–25 and nuclear structure theory.^26–28 Many strengths of the single-reference CC formalism (SR-CC) originate in the exponential representation of the ground-state wave function |Ψ⟩,

where T and |Φ⟩ correspond to the cluster operator and a reference function that provides an approximation to the exact ground-state wave function. For practical applications, one can define a hierarchy of CC approximations by increasing the excitation level in the cluster operator. Another important feature of the CC formalism stems from the linked cluster theorem,^29,30 which allows one to build efficient algorithms for the inclusion of the higher-rank excitations. When these two features are combined, they give rise to efficient and accurate methodologies that account for higher-order correlation effects and that have been widely used in physics and quantum chemistry. The diverse manifold of higher-order approximations include, among others, categories of active-space, multi-reference, externally corrected, and adaptive approaches as well as non-iterative corrections and approaches that combine stochastic techniques (for key reviews and papers, see Refs. 7, 9, and 31–39 and references therein).

Recently, it was shown using the sub-system embedding sub-algebras (SES)⁴⁰ approach that the CC formalism is a natural renormalization procedure that allows one to construct effective Hamiltonians in reduced-dimensionality active spaces for classes of sub-systems of the quantum system of interest. This formalism provides a mathematically rigorous procedure to construct the many-body form of the active-space effective Hamiltonians using only external (out of active space) fermionic degrees of freedom that can be extracted from single reference CC wave function expansions. In contrast, the internal degrees of freedom (inside active space) are determined by the diagonalization of effective Hamiltonians in the active space. The SES-CC formalism is also valid for any sub-system, and the resulting sub-problems (or computational blocks) can be integrated into the computational flows (or quantum flow algorithms described in Ref. 41), which can traverse large sub-spaces of the corresponding Hilbert space. For example, using this approach, one can define flows involving occupied localized orbitals, bypassing certain problems of traditional local CC formulations. Additionally, SES-CC can be extended to the time domain, leading to time-dependent Hamiltonians that describe the dynamics of the entire system in the active space (provided that time-dependent external amplitudes are known or can be efficiently approximated/evaluated).

In light of the above discussion, high-accuracy CC-based techniques for novel representations of quantum many-body problems in reduced-dimensionality spaces may play an important role not only in solving complex problems but also in enabling new forms of simulations associated with the emergence of the quantum computing (QC) technologies. Although reduced-dimensionality methods have been existing for a long time and are usually identified with Complete Active Space Self-Consistent Field (CASSCF) methods (see, for example, Refs. 42 and 43 and references therein), density matrix renormalization group (DMRG) methods,^44–47 local CC formulations,^48–60 and correlated approaches where the energies are obtained through the diagonalization of the effective Hamiltonians in reduced-size sub-spaces (or equivalently-model spaces),^61–91 the SES-CC formulation allows one to construct active-space effective Hamiltonians using single-reference concepts. A problem with SES-CC effective Hamiltonians in certain applications, such as quantum computing, is their non-Hermitian character. In light of this problem, we have recently explored the double unitary coupled cluster (DUCC) ansatz,^92–96 which, in analogy to the SES-CC techniques, results in an active-space representation of quantum problems with the added benefit that the resulting DUCC effective Hamiltonians are Hermitian.

In analogy to the canonical SES-CC approach and various effective Hamiltonian formulations based on the wave operator formalism and/or similarity transformations, the DUCC formalism also decouples the external and internal fermionic degrees of freedom. These Hamiltonians have been recently studied in the context of quantum simulations using various quantum solvers, including Quantum Phase Estimator (QPE)^97–103 and Variational Quantum Eigensolvers (VQEs),^104–118 where representations of the DUCC effective Hamiltonians used in these studies were defined by the lowest-order contributions to the corresponding commutator expansion. The area of similarity transformations of electronic Hamiltonians induced by exponential operators involving cluster operators expressed in terms of non-commutative algebras is under active development (for recent developments, see Refs. 119–121).

In this manuscript, we will extend this analysis to DUCC effective Hamiltonians that include higher-order correlations effects. In particular, we will focus on approximations involving double commutators. The main focus will be on the approximate schemes where the class of external excitations used to define DUCC effective Hamiltonians includes singly and doubly excited amplitudes that are used to construct one- and two-body effective interactions. Within this approximation, we consider algorithms that are correct up to the third order of the standard Møller–Plesset perturbation theory. We illustrate the performance of these methods on the example of several benchmark systems used in recent studies of quantum algorithms based on the downfolding formalisms. The paper focuses on examining the role of higher-order terms in the commutator expansion. Therefore, we chose model systems that reduce concern for missing three-body effects and higher-order excitations. This includes the Be, Li₂, and H₂O systems. Additionally, the Li₂ system epitomizes the situation where the coupled cluster single double (CCSD) level of theory provides an accurate approximation to the exact [full configuration interaction (FCI)] formalism, and the active space can be properly defined. Using this example, we demonstrate that in this situation, the approximations defined by double commutators can provide a systematic improvement in the accuracy of energies obtained by the diagonalization of the effective Hamiltonians defined by single-commutator-based approaches. For the H₂O system, we discuss the behavior of the DUCC energies as functions of active space size. These results provide yet another illustration of the feasibility of the CC downfolding techniques.

The DUCC effective Hamiltonian formalism has been amply discussed in recent papers.^41,92,93,95 Here we only describe basic tenets of this approach. The DUCC formalism utilizes a composite unitary CC ansatz to represent the exact wave function |Ψ⟩, i.e.,

where σ_ext and σ_int are general-type anti-Hermitian operators,

All cluster amplitudes defining the σ_int cluster operator carry active spin-orbital indices only and define all possible excitations within the corresponding active space. The external cluster operator σ_ext is defined by parameters carrying at least one inactive spin-orbital index (for more details, see the discussion of DUCC expansion in Ref. 93).

Using the DUCC representation [Eq. (2)], it can be shown that the energy of the entire system (once the exact form of the σ_ext operator is known) can be calculated through the diagonalization of the effective/downfolded Hamiltonian in the active space, i.e.,

where

and

In Eq. (6), P and Q_int are the projection operators on the reference function |Φ⟩ and all excited Slater determinants (defined with respect to the |Φ⟩ reference) in the active space, respectively. Once the external cluster amplitudes are known (or can be effectively approximated), the energy of the system can be evaluated by diagonalizing Hermitian effective/downfolded Hamiltonian in the active space using various quantum/conventional-computing solvers. Important steps toward developing practical computational schemes based on the utilization of the DUCC downfolded Hamiltonians are (1) efficient approximation of the σ_ext operator and (2) approximate form of the non-terminating commutator expansions defining downfolded Hamiltonians. A legitimate approximation for σ_ext is through the utilization of the external part of the standard cluster operator

which has been discussed in Ref. 92. In particular, T_ext can be approximated by SR-CCSD amplitudes that carry at least one external spin-orbital index. Other possible sources for evaluating external cluster amplitudes are higher-rank single-reference CC methods, truncated expansions,^119,122 and approximate unitary CC formulations such as UCC(n) methods.^123,124

Given the fact that the $H̄ext$ operator is defined as an infinite expansion in terms of the σ_ext operator, we cannot calculate its exact form. For practical reasons, in this paper, we implement a three-step simplification procedure. The first step, denoted as a simplification S1, is the estimation of σ_ext through the canonical CCSD external amplitudes, i.e., $σext≃Text(CCSD)−Text(CCSD)†$⁠. Although approximating cluster operators involving non-commuting components is a challenging task, the analysis of Ref. 123 suggests that the above approximation is justified. As possible alternatives to the CCSD amplitudes, one can consider amplitudes obtained with various UCC(n) models (amply discussed in Ref. 123) based on the perturbative analysis of the unitary CC energy. An interesting feature of the UCC(n) formalisms is that the UCC(3) equations for cluster amplitudes are equivalent to the linear CCD (CC with doubles) equations. The second approximation, denoted as S2, is to limit $H̄ext$ to one- and two-body contributions. The last simplification, which we will explore in more detail in this work, is the truncation of the commutator expansion of $H̄ext$ (denoted as S3).

In this work, in addition to the single commutator expansion used in the previous analysis of the CC downfolding,⁹² we include non-trivial contributions stemming from the double commutators. We collected all the approximate versions of the downfolded Hamiltonians in Table I, where we analyze seven approximate schemes A(1)–A(7). In the expressions for A(1)–A(7),

and

where p, q, r, and s are spin-orbital labels, $hqp$ and $vrspq$ represent one-electron and antisymmetric two-electron integrals, and

where F_N and V_N are one- and two-body components of the H_N operator. The commutators terms in Table I are in the particle–hole normal product form and we refer the reader to Ref. 92 for the procedure on how to obtain the corresponding Hamiltonian in physical-vacuum normal product form. The first approximation, A(1), utilizes bare Hamiltonian in the active space while approximations A(2)–A(4) are driven by single commutators [in approximations A(2) and A(4), F_N-dependent double commutator terms are added to provide perturbative consistency to the second order of many-body perturbation theory (MBPT)]. Approximations A(5)–A(7) introduce non-trivial terms due to the second commutators. In approximations A(5) and A(7), in analogy to A(2) and A(4), F_N-dependent triple commutators are introduced to provide perturbative consistency at the third order of perturbation theory [MBPT(3)].

The perturbative analysis of the downfolded Hamiltonians involves the perturbative expansion for the T_ext operator. It should be, however, remembered that due to the construction of the active space and the definition of the external part of the cluster operator, the riskiest amplitudes characterized by large values are not explicitly present in T_ext. Therefore, the perturbative analysis with T_ext may not involve typical problems encountered by the full T operator for the challenging cases defined by the presence of strong correlation effects.

The efficiency of commutator expansion also hinges upon the judicious choice of the active space. When the active space adequately captures the static correlation effects leading to small values of σ_ext amplitudes; the commutator expansion should provide a hierarchical class of approximations, otherwise, infinite commutator expansion involving large values of external amplitudes may lead to divergent nature of expansions based on this scheme.

The non-Hermitian and Hermitian formulations of the CC downfolding can also be used as tools in analyzing techniques focusing on the partitioning of the CC correlation effects into internal and external parts.^125,126 For simplicity, we will assume that the external part of the cluster operator is defined by doubly excited external amplitudes (T_{ext, 2}) taken, for example, from the CCSD calculations. Let us also focus on the A(2) approximation, where we retain only fully contracted diagrams stemming from [H_N, σ_ext] and $12[[FN,σext],σext]$ terms, where $σext≃Text, 2−Text, 2†$⁠. These diagrams are (for simplicity, we assumed HF orbitals) given by the following expressions:

and

where $t̄ijab$s are external cluster amplitudes, which are equal to zero if all indices i, j, a, b are active, and ɛ_p’s are the HF orbital energies. The expressions in Eqs. (12) and (13) can be identified with the external correlation energies (it is also the primary source of double-counting of correlation effects if we apply approximations relying solely on terminating at a particular rank of the commutator expansion). The expression in Eq. (13) is equivalent up to the second-order to the negative value of the external correlation energy. Therefore, assuming the cancellation of the first term in (12) and the term in Eq. (13) leads to the following form of the downfolded Hamiltonian:

where $ΔEext(CC)$ is the external correlation energy $(ΔEext(CCSD))$⁠. In this case, the eigenvalue of the active-space Hamiltonian (14) is a sum of the CAS energy and the external part of the CCSD energy. The same result can be obtained using non-Hermitian CC downfolded Hamiltonian (see Ref. 40),

where $e−TextHeText$ can be expressed using finite commutator expansion as

Again, focusing on single and double commutators and retaining only their scalar components leads to the following form of the effective Hamiltonian:

which again furnishes the same energy as in the Hermitian case. In the following part of the discussion, we will refer to this model as A(1) + $ΔEext(CC)$⁠. In particular, we will focus on the case where the CCSD formalism is used to evaluate the external part of the correlation energy [A(1) + $ΔEext(CCSD)$].

A simple approximation discussed above does not describe the impact of the coupling between internal and external cluster amplitudes (T_int/T_ext or σ_int/σ_ext) as defined by the full form of many-body expansions for the effective Hamiltonians. Another class of approximations can be associated with using approximate forms of the external part of cluster amplitudes at the MBPT(2) level as discussed in Ref. 126. By (1) using the approximate form of the T_ext operator in Eq. (15), (2) taking into account the equivalence at the solution between the connected form of the CC equations and effective Hamiltonian [Eq. (15)] eigenvalue problem,^40,41 (3) approximating internal cluster amplitudes by singles and doubles, and (4) solving equations corresponding to active-space singles and doubles, one arrives at the energy expansion identical with the one obtained in Ref. 126.

Although there are alternative ways for analyzing the internal and external parts of the CC energy, the CC downfolding offers a simple way of introducing these approximations using the language of effective Hamiltonians.

The implementations of A(2)–A(7) formulations have been integrated with the Tensor Contraction Engine (TCE)¹²⁷ environment of NWChem,^128,129 which allows us to use a variety of formulations to extract the external cluster operator. In the present studies, we use the TCE CCSD module to provide singly and doubly excited parts of the external cluster operator (T_ext). Calculations were performed using restricted Hartree–Fock (RHF) and natural orbitals at the MBPT(2) level (for the Li₂ system only). Additionally, in all calculations, all occupied spin orbitals were considered active.

To evaluate the quality of downfolding CC methods, we compare their energies with the results of the coupled cluster single double triple (CCSDT) and coupled cluster single double triple quadruple (CCSDQT) methods correlating only active orbitals. The CCSDT and CCSDTQ results in the active space correspond to having only internal amplitudes without accounting for the external excitations through an effective downfolded Hamiltonian. There are a couple of reasons for this:

• To measure how much correlation is captured by the downfolding procedure. Since some works, especially in quantum computing, truncate higher-energy orbitals, these numbers illustrate the effect of negating the dynamical correlation and why one should not negate the higher-energy orbitals.

• To measure the effect of higher-order excitations in the active space as the diagonalization of the DUCC effective Hamiltonian will capture those excitations. This helps to rule out any fortuitous agreement with full CC methods due to the approximations used.

The calculations for the beryllium atom were performed using cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets in active spaces composed of 5, 6, and 9 lowest-lying RHF orbitals. The results of our simulations are shown in Fig. 1 and Tables II and III. In Fig. 1, we show the errors of approximations A(1)–A(7) obtained with various active spaces with respect to the FCI results obtained when correlating all RHF orbitals in cc-pVXZ (X = D, T, Q) basis sets. One can observe a consistent pattern emerging: (1) the behavior of downfolded methods is similar for all basis sets employed, (2) the accuracy of the energies obtained with the effective Hamiltonians increases as the size of the active space increases and as higher-order commutators are included in the effective Hamiltonian expansion. These observations are corroborated by the results in Table II, where total energies are compared, and Table III, where the percentage of total correlation energy reproduced by a given downfolded Hamiltonian is shown. The importance of the higher-order terms stemming from the double commutator is well illustrated by the A(5)–A(7) methods for small (five orbitals) active space in the cc-pVQZ basis set (Table III), where these approaches can reproduce the exact correlation energy. For this specific example, the performance of the single-commutator-driven approaches A(2)–A(4) is significantly worse. This is associated with the quality of the active space, which is reflected by the fact that diagonalizing bare Hamiltonian in the active space defined by five orbitals reproduces only 18% of the total ground-state correlation energy. Using larger active space for the cc-pVQZ basis set, defined by nine lowest-lying RHF orbitals, results in single-commutator-driven approaches A(2) and A(4) performing much better than in the case of small five-orbital active space. As discussed in Refs. 92 and 94, approximation A(3) is not fully consistent at the second order of perturbation theory and may result in doubling the external correlation effects, which is reflected by the results shown in Fig. 1. One should also observe systematic improvements of the quality of the A(2) and A(4) formulations with the size of the active space employed, which is a consequence of the inclusion of important F_N dependent double commutator terms. For larger active spaces, double-commutator-driven approximations A(5)–A(7) approach the FCI level of accuracy. The utilization of the larger active spaces is also advantageous from the point of minimizing the effects of additional approximations such as neglecting three-body interactions in the effective Hamiltonians or using an approximate form of the σ_ext operator. These approximations are what may account for why the results obtained with five orbitals are in a slightly better agreement with FCI than some larger active space results, but this discrepancy takes a backseat to the importance of higher-order commutators.

In this subsection, we will compare the performance of the A(i) (i = 2, …, 7) approximations for the Li₂ system in the cc-pVTZ basis set (spherical harmonic representation of d orbitals was employed) containing 60 basis set functions. In analogy to previous studies⁹⁴ that demonstrated advantages of using correlated molecular orbital basis, we will use the natural orbitals calculated at the second order of Møller–Plesset perturbation theory [MBPT(2)]. To define active spaces, we used ten MBPT(2) natural orbitals characterized by the largest occupancies. The results of calculations are shown in Fig. 2 and Table IV. In Fig. 2, we show the potential energy surfaces A(1), A(4), and A(7) and CCSDT formulations (the values of calculated energies are shown in Table IV). It should be stressed that the full CCSDT formalism for the Li₂ system provides nearly FCI level of accuracy in the calculated ground-state energies. We also included an additional model in which the external correlation energy computed using external amplitudes from CCSD, $ΔEext(CCSD)$⁠, is added to the energy from the bare Hamiltonian in the active space [A(1)]. From Fig. 2, one can see that energies obtained with the A(1) approach significantly underestimate dissociation energy (errors for large separations are significantly larger than the errors near the equilibrium geometry). While these errors are improved upon by adding the external correlation energy $ΔEext(CCSD)$⁠, such disjoint extrapolation techniques can lead to an imbalance in correlation energy, which has been observed in previous studies as well.¹³⁰ This is exemplified by the change from overestimation of the correlation energy near equilibrium to underestimation of the correlation energy at stretch bond lengths. One should notice that the A(7) formalism provides a very balanced description of correlation effects and provides further improvement of the A(4) results. It also provides smaller errors that the A(1) + $ΔEext(CCSD)$ extrapolated results for most geometries. This fact is well illustrated in Table IV, where one can observe that the A(7) errors with respect to the full CCSDT energies do not exceed 2.5 mhartree. The same error bound for A(4) formulations is larger and amounts to 6.9 mhartree. It is also worth emphasizing that the high quality of the A(7) results was achieved despite the fact that the MBPT(2) formalism may not be the method of choice for large Li–Li separations and MBPT(2) natural orbitals may significantly differ from those defined by non-perturbative formalisms such as the CCSD approach.

Breaking a single bond in the water molecule provides yet another example for illustrating the performance of the downfolding methods in describing chemical transformations. We also use this benchmark to analyze the convergence of the commutator expansions as a function of the active space size. In the present studies, we use the equilibrium geometry of the water molecule defined by the following parameters: R_OH = R_e = 0.961 83 Å and ∠_HOH = 103.9215°. In Table V, we compare the energies obtained with various downfolding methods [A(2)–A(7)] with CCSDTQ and CCSDTQ-in-active-space formulations for active space defined by 12 active RHF orbitals. The CCSDTQ-in-active-space helps in calibrating the importance of the correlation effects originating in the orthogonal complement of the active space. For all geometries considered here, we used the cc-pVTZ basis set (employing the spherical representation of the d orbitals), which yields 58 RHF orbitals. From Table V, one can see that for all geometries, the A(6) and A(7) approximations provide energies of similar values, which are consistently located above the CCSDTQ ones. For each geometry considered in Table V, A(6) and A(7) provide significant improvements in the CCSDTQ-in-active-space energies. For example, the CCSDTQ-in-active-space errors with respect to the CCSDTQ energies where all orbitals are correlated amount to 234, 230, and 227 mhartree for 1.0R_e, 1.5R_e, and 2.0R_e, respectively. For comparison, the analogous A(7) errors take values 3, 10, and 22 mhartree. It should also be stressed that the perturbative A(5) approximation consistently underestimates the CCSDTQ energies by 5–7 mhartree. The double commutator expansions A(5)–A(7) also improve the quality of single-commutator formulations [including the most complete A(4) formulation].

To explore the impact of the size of the active space on the quality of downfolding formulations, we performed a series of calculations with the A(2)–A(7) approaches for the equilibrium structure of the water molecule for active spaces defined by 6, 7, …, 16 active orbitals corresponding to the lowest-lying RHF molecular orbitals (see Table VI). For the water structure where one OH bond is stretched to 2.0R_e, we perform similar studies with representative A(4) and A(7) approaches for active spaces defined by 12–16 active orbitals (see Table VII). It is justified to expect that the increase in the active space size can compensate for simplifications S1–S3 (discussed in Sec. II), which are indispensable in dealing with non-terminating series. For example, part of the effects due to external triply (T_{ext, 3}) and quadruply (T_{ext, 4}) excited amplitudes for small active spaces, which are not included in the discussed approximations, can be accommodated in the FCI-type solvers (quantum and classical) as internal excitations for larger active spaces. Indeed, by scrutinizing Table VI and Fig. 3, one can observe that for larger active space comprising 16 orbitals (which still provides a nearly four-fold reduction in the total size of the orbital space), the A(7) results provide nearly the CCSDTQ level of accuracy. At the same time, A(7) results are more accurate than the energies obtained with single commutator expansion A(4) (see Fig. 3) and converge quickly to the CCSDTQ level of accuracy. This is happening despite the fact that a very simple form of external T_ext is used and that the CCSDTQ-in-active-space energy for active space defined by 16 orbitals is characterized by 193 mhartree of error with respect to the full CCSDTQ energy.

In Table VII, we can observe a similar behavior for the 2.0R_e case. In this situation, the error of the A(7) energies with respect to the CCSDTQ data for 16 orbital active space is characterized by a larger error (6.746 mhartree) compared to the equilibrium case. However, given the reduction of 184.505 mhartree error obtained with the CCSDTQ-in-active-space approach, one should consider A(7) accuracy as quite satisfying given the number of simplifications (S1–S3) used to define the corresponding A(7) downfolded Hamiltonian. It is also instructive to observe that the quality of the 16-active orbital A(7) energy is significantly better than the one obtained with the full CCSD approach (17.839 mhartree of error), which provides the source for the external amplitudes.

In this paper, we discuss the level of accuracy provided by the active-space downfolded Hamiltonians, including variants involving non-trivial terms stemming from double commutators. In analogy to previous studies, for the sake of perturbative consistency, in some approximations, we also included Fock-operator-dependent terms stemming from triple commutators. In all approximations analyzed in this paper, we used a subset of CCSD cluster amplitudes to evaluate the so-called external part of the anti-Hermitian cluster operator σ_ext. To assess the performance of these approximations, we performed a series of DUCC simulations for typical benchmark systems such as the beryllium atom, the Li₂ dimer, and the water molecule. Using these examples, we demonstrated that the inclusion of double-commutator terms results in consistent improvements in the DUCC approaches based on a single-commutator expansion. This was best seen on the example of a bond breaking in the Li₂ system when the integration of the A(7) approach with natural orbitals for modest size active space resulted in nearly the CCSDT level of accuracy for all geometries considered. Using the H₂O example, we also demonstrated that the effect of S1–S3 approximations needed to deal with the non-terminating series can be, to a large extent, accommodated by enlarging the size of the active space. At the same time, we demonstrated that A(7) energies quickly approach the full CCSDTQ results as the size of the active space increases. The double-commutator-based approximation A(7) not only outperforms the single-commutator-based A(4) approach but can also provide CCSDTQ level quality for active spaces, which offer nearly four-fold reduction in the size of the orbital space employed.

A natural extension of the current studies is to include higher-many-body components in the downfolded Hamiltonians. In particular, in the forthcoming papers, we will include (1) three-body interactions (in the particle–hole representation) stemming from contractions of the two-body part of the electronic Hamiltonian with σ_{ext 2} operators and (2) terms due to the σ_{ext 3} and/or σ_{ext 4} operators.

We believe that energy accuracies obtainable with downfolding techniques are getting us close to quantum simulations of realistic chemical processes. The downfolding Hamiltonians also open an opportunity for utilizing machine learning methods to extract the analytical form of the effective inter-electron interactions in small-dimensionality sub-spaces of the entire Hilbert space.

This work was supported by the “Embedding QC into Many-body Frameworks for Strongly Correlated Molecular and Materials Systems” project, which is funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences (BES), the Division of Chemical Sciences, Geosciences, and Biosciences. Part of this work was supported by the Quantum Science Center (QSC), a National Quantum Information Science Research Center of the U.S. Department of Energy (DOE).

The authors have no conflicts of interest to declare.

The data that support the findings of this study are available from the corresponding author upon reasonable request.