Toward DMRG-tailored coupled cluster method in the 4c-relativistic domain

At the turn of the millennium, the density matrix renormalization group (DMRG) method¹ was introduced to the quantum-chemical community,^2–4 and since then, it has seen a large surge in the use for multireference systems. The biggest advantage of the DMRG method is its capability to treat large active spaces, and current implementations can go to about 50 active space spinors.^5,6 However, a major drawback of DMRG is its inability to capture dynamical correlation, since it cannot include all virtual spinors. This correlation has a strong influence on the target systems of this project, which thus aims to address this problem. The DMRG method is already well established and computational chemists started to use it; however, the methods for treating the dynamical correlations on top of DMRG are still in the pioneering stage. Past efforts were either based on second order perturbation theory,⁷ internally contracted MRCI (multireference configuration interaction),⁸ random phase approximation,⁹ canonical transformation method,¹⁰ or the perturbation theory with matrix product states (MPS).¹¹ 

Our group has followed a different pathway to deal with the dynamical correlation: the coupled cluster (CC) method externally corrected by DMRG.¹² As the name suggests, this is a combination of DMRG and the coupled cluster (CC) method. The CC method is known for its ability to describe dynamical correlation. In the externally corrected approach, first a DMRG calculation is performed on the strongly correlated active space, keeping the rest of the system fixed. This accounts for the static correlation. The second step is the CC analysis of the matrix product state (MPS) wave function, obtained from DMRG. Then, a CC calculation is performed on the rest of the system, keeping, in turn, the active space amplitudes fixed, which captures the dynamical correlation. Already the simplest version thereof, the tailored CCSD (CC with single and double excitations) approach,^13,14 yields very promising results.¹² Remarkably, all previous approaches based on the use of DMRG output in another method have so far been non-relativistic, leaving the relativistic domain unexplored. This is the focus of this paper.

First, we demonstrate the capabilities of our relativistic 4c-TCCSD implementation on the example of the thallium hydride (TlH) molecule, which has become a standard benchmark molecule for relativistic methods and most importantly large-scale DMRG and up to CCSDTQ (CC with single, double, triple, and quadruple excitations) results are available.¹⁵ It should be noted that DMRG is best suited for static-correlation problems, while TlH is dominated mostly by dynamic correlation, for which CC approaches are excellent. In order to study the behavior of the TCCSD method for more multireference systems, we performed tests on AsH and SbH molecules. The multireference character in the AsH and SbH ground state arises from the fact that two determinants are needed to describe the M_s = 0 triplet component in the spin-free case and therefor is somewhat “artificial.” Figure 1 depicts the $π1/22$ and $π3/22$ determinants arising from M_s = 0 triplet component’s determinants due to the spin–orbit splitting for AsH and SbH. The heavier the atom is, the greater the splitting, and then, $π1/22$ is more dominant and the X0⁺ ground state becomes less multireference. Hence, we expect AsH to be of stronger multireference nature than SbH.

Present-day relativistic calculations are often carried out within the no-pair approximation, where the Dirac–Coulomb Hamiltonian is embedded by projectors eliminating the troublesome negative-energy solutions, which yields a second quantized Hamiltonian formally analogous to the non-relativistic case,

where the indices P, Q, R, and S run over the positive-energy four-component spinors spanning the one-electron basis. The barred spinors (⁠$φp¯$⁠) and unbarred spinors (φ_p) form Kramers pairs related to each other by action of the time-reversal operator K,

The Kramers symmetry replaces the spin symmetry in the non-relativistic theory; in particular, M_S is not a good quantum number and M_K projection is defined instead, which is $12$ for unbarred spinors (A) and −1/2 for spinors with barred indices (B). The capital indices in (1) run over both spinors of a Kramers pair. In contrast to the non-relativistic case, the Hamiltonian (1) is, in general, not block-diagonal in M_K. Since each creation or annihilation operator in (1) changes M_K by ±1/2, the Hamiltonian couples states with |ΔM_K| ≤ 2. Moreover, the index permutation symmetry of the 2e-integrals in (1) is lower than in the non-relativistic case.

The Dirac program¹⁶ employs a quaternion symmetry approach that combines the Kramers and binary double group symmetry (⁠$D2h*$ and subgroups).¹⁷ The double groups can be sorted into three classes based on the application of the Frobenius–Schur indicator to their irreducible representations: “real groups” (⁠$D2h*$⁠, $D2*$⁠, and $C2v*$⁠), “complex groups” (⁠$C2h*$⁠, $C2*$⁠, and $Cs*$⁠), and “quaternion groups” (⁠$Ci*$ and $C1*$⁠).¹⁸ Generalization of non-relativistic methods is simplest in the “real groups” case, where the integrals are real-valued and the ones with odd number of barred (B) indices vanish. In practice, it means that additional “spin cases” of integrals (AB|AB) and (AB|BA) (in Mulliken notation) have to be included. For the complex groups, the integrals are complex-valued, but still, only integrals with an even number of barred indices are non-zero. Finally, in the remaining case of “quaternion groups,” all the integrals have to be included and are complex-valued.^18,19

The idea of externally corrected coupled cluster methods is to take information on the static correlation from some non-CC external source and to include it into the subsequent CC treatment.²⁰ The conceptually simplest approach is the tailored CC (TCC) method proposed by Bartlett et al.,^13,21–23 which uses the split-amplitude ansatz for the wave function introduced by Piecuch et al.,^24,25

where T_cas containing amplitudes with all active indices is “frozen” at values obtained from the complete active space configuration interaction (CASCI) or, in our case, from DMRG. The external cluster operator T_ext is composed of amplitudes with at least one index outside the CAS space. Another way to justify this ansatz is the formulation of CC equations based on excitation subalgebras recently introduced by Kowalski.²⁶ The simplest version of the method truncates both T_cas and T_ext to single and double excitations. Since there is a single-determinantal Fermi vacuum, the excitation operators T_ext and T_cas commute, which keeps the method very simple. TCC can thus use the standard CCSD solver, modified to keep the amplitudes from T_cas fixed. Thanks to the two-body Hamiltonian, tailored CCSD energy with T_ext = 0 and T_cas from CASCI (complete active space configuration interaction) reproduces the CASCI energy. In the limit of CAS space, including all MOs, TCC thus recovers the full configuration interaction (FCI) energy. In general, error bound valid for tensor network state TCC (TNS-TCC) methods is derived.²⁷ 

In Refs. 12 and 28, we have described how to obtain T_cas from the DMRG wave function using concepts of quantum information theory²⁹ in the non-relativistic case, yielding the DMRG–TCCSD method. The DMRG method³⁰ is a procedure that variationally optimizes the wave function in the form of the matrix product state (MPS) ansatz.⁵ The quantum chemical version of DMRG (QC-DMRG)^31–36 eventually converges to the FCI solution in a given orbital space, i.e., to CASCI. The practical version of DMRG is the two-site algorithm, which provides the wave function in the two-site MPS form,⁵ 

where α_i ∈ {|0⟩, |↓⟩, |↑⟩, |↓↑⟩}, and for a given pair of adjacent indices [i, (i + 1)], W is a four index tensor, which corresponds to the eigenfunction of the electronic Hamiltonian expanded in the tensor product space of four tensor spaces defined on an ordered orbital chain, so called left block (M_l dimensional tensor space), left site (four dimensional tensor space of the ith orbital), right site [four dimensional tensor space of the (i + 1)th orbital], and right block (M_r dimensional tensor space).

When employing the two-site MPS wave function [Eq. (4)] for the purposes of the TCCSD method, the CI expansion coefficients $cia$ and $cijab$ for a, b, i, j ∈ CAS can be efficiently calculated by contractions of MPS matrices.^37,38 We would like to note that using the two-site DMRG approach, in practice, means using the wave-function calculated at different sites, and it can only be employed together with the dynamical block state selection (DBSS) procedure⁴ assuring the same accuracy along the sweep. Alternatively, one can use the one-site approach in the last sweep.³⁹ 

Once the CI coefficients $cia$ and $cijab$ have been obtained, the standard CC analysis is performed to convert them to the CC amplitudes,

The generalization of the DMRG–TCCSD method to the relativistic 4c case has to consider several points. First of all, the additional integral classes with nonzero ΔM_K have to be implemented in the DMRG Hamiltonian.^15,40 Second, there will be more CI coefficients and, subsequently, CC amplitudes to be obtained from the MPS wave function, corresponding to excitations with nonzero ΔM_K. Finally, except for the “real groups,” the DMRG procedure has to work with complex matrices, and the resulting cluster amplitudes will also be complex-valued. In the present work, we have selected numerical examples with “real groups” symmetry, while the complex generalization of the DMRG code is in progress.

In the present work, we have used the two-site DMRG variant together with DBSS through the course of the whole DMRG procedure and obtained the CC amplitudes from the resulting two-site form of the MPS wave function. For all systems, DMRG calculations were performed exclusively with the QC-DMRG-Budapest program.⁴¹ The Dirac program package¹⁶ was used for the remaining relativistic calculations, whereas the Orca program was used for remaining non-relativistic calculations. Orbitals and MO integrals were generated with the Dirac program package.¹⁶ We used the Dirac–Coulomb Hamiltonian and triple-zeta basis sets for the heavier of the two atoms (cv3z for Tl, As, and Sb and cc-pVTZ for F) as well as for hydrogen (cc-pVTZ), which include core-correlating functions for the heavier atom. Initialization of DMRG, i.e., optimal ordering of spinors, was set up, as discussed in Ref. 35. The numerical accuracy was controlled by DBSS,⁴ keeping up to thousands of block states for the a priori set quantum information loss threshold χ = 10⁻⁶.

In order to compare with the non-relativistic version of the TCCSD method, the system of hydrogen fluoride was chosen as it is biatomic with light nuclei. CC-pVTZ basis was used at the internuclear distance of 0.8996 Å.

Consistent methods should exhibit a constant shift of relativistic and non-relativistic energy ΔE = E_rel − E_nonrel, given by a different Hamiltonian. Table I shows that TCCSD is consistent with RHF, CCSD, and DMRG methods in terms of ΔE up to a millihartree.

We have used the computational protocol of Ref. 15 for direct comparison with their energies. $C2v⋆$ double group symmetry with real irreps was assumed. The 4c-RHF energy was −20 275.416 61 E_h. We used MP2 natural spinors (NS) from the Dirac program¹⁶ as the spinor basis for electron-correlation calculations, correlating the Tl 5s, 5p, 4f, 5d, 6s, 6p, and H 1s electrons, while keeping the remaining core electrons of Tl frozen. Using uncontracted basis, a virtual spinor threshold was set at 135 E_h. The resulting space (14, 47) was chosen by ordering MP2 NS by their occupations and taking those with values between 1.98 and 0.001. In this space, the 4c-TCCSD was performed, with DMRG calculations in the procedure limited to subspaces of (14, 10), (14, 14), (14, 17), (14, 25), and (14, 29), with spinors sorted by MP2 occupations. Figure 2 shows a scheme of embedded active spaces used in the procedure.

Since SbH is a heavier homolog of AsH (which is itself a homolog of nitrene, NH), the procedure was very similar for both of them. In contrast, with TlH, instead of MP2 NS, we used average-of-configuration SCF spinors. The 4c-SCF energy was −2260.042 61 E_h for AsH and −6481.107 75 E_h for SbH. The DMRG calculation in the 4c-TCCSD procedure was limited to subspaces of (16, 14) and (16, 23).

Dominant contributions to active spaces for the AsH active space (16, 14) are As (3d, 4s, 4p, 5s, 5p and H: 1s, and for active space (16, 23), we add As: 4d and H: 2s, 2p to the former. Energies of MOs are from −2.1 E_h to +0.25 E_h and +1.01 E_h for (16, 14) and (16, 23), respectively (for internuclear distance 1.52 Å). Dominant contributions to active spaces for SbH are analogous, except for principal quantum numbers, which are higher by 1.

Once we reproduced the MP2 and CCSD energy of TlH in equilibrium geometry from Ref. 15, we applied the 4c-TCCSD method. Obtained energies and their respective deviations from the reference CCSDTQ calculation¹⁵ are listed in Table II. In the case of the optimal selection of active space of 14-spinors, the TCCSD method improved the CCSD energy by 4.94 mE_H. While TCCSD introduces only a minor computational cost increase over CCSD, it cuts the energy error in half. This shows the practical advantage of the method. The energy obtained by TCCSD is comparable even with large-scale DMRG in the full CAS(14, 47).

As we can see from the high accuracy of the 4c-CCSD(T) energy, the system does not exhibit a considerable multireference character. Therefore, even a rather small CAS of 14 spinors is sufficient for a good description of the system.

As shown on the chart in Fig. 3(a), TCCSD significantly improves the DMRG energy toward FCI, even for the smallest CAS space. In fact, further enlarging of CAS over the size of 14 spinors is counterproductive. Although the TCC must reproduce the FCI energy when CAS is extended to all spinors, the TCC energy does not approach this limit monotonically.⁴² The obvious reason is that the “frozen” T_CAS amplitudes cannot reflect the influence of the dynamical correlation on the external space back on the active CAS space; therefore, extending CAS space first exacerbates the results. Unfortunately, more detailed understanding of this highly nonlinear behavior is still an open problem. Despite a considerable effort, so far, no quantities able to predict the optimal CAS-EXT split a priori have been identified.⁴² Nevertheless, we have found that an error minimum can be obtained by sweeping through the entire orbital space with low cost DMRG calculations, and this minimum does not shift or shifts only a little when more accurate calculations are performed. Therefore, in practice, the optimal CAS size related to the energy minimum is usually independent of M and can be determined with low bond dimension (M) DMRG calculations.⁴² 

As demonstrated by the chart in Fig. 3(b), the optimal CAS size is 14 spinors for the equilibrium energy calculation. This CAS size is optimal not only for energies but also for the calculation of spectroscopic properties, including the low bond dimension calculations with M = 512 (see Table III).

The obtained spectroscopic properties of TlH are listed in Table III and the corresponding dissociation curve is depicted in Fig. 4. Even for a small active space of 10 spinors, TCCSD shows an agreement with the experiment comparable with the large DMRG(14, 47) calculation. For the 14-spinor space, spectroscopic constants obtained by TCCSD exhibit the best agreement with the experiment, thus being consistent with the lowest energy single point result of the 14-spinor space in Table II. Moreover, TCCSD based on DMRG with M = 512 states or on DMRG with DBSS yields very similar results, indicating that the underlying DMRG is well converged. However, for 17-spinors, there is a bigger difference and M = 512 might not be accurate enough. This is in accordance with the previous findings of dynamical correlation effects.

In order to further assess the feasibility of the method for multireference systems, we studied AsH and SbH molecules. Table IV compares the determinants with highest coefficients as calculated by 4c-DMRG with the space of (14, 14) for TlH and (16, 14) for AsH and SbH. The coefficients confirm that the need arises for a multireference description of the ground state of AsH and that SbH is between AsH and TlH in terms of its multireference character.

The equilibrium energies of AsH in Table V show that 4c-TCCSD improved the 4c-CCSD by 17 mE_h, with just 3 mE_h difference from 4c-CCSD(T). In this case, a more accurate theoretical reference energy is unavailable, and therefore, we turned to spectroscopic constants, as shown in Table VI, to enable a comparison with accurate IR spectra obtained by CO laser magnetic resonance in Ref. 44. The respective potential curve is plotted in Fig. 5. The comparison of both the internuclear distance r_e and the vibrational constant ω_e shows a clear advantage of 4c-TCCSD over 4c-CCSD in this case, which we attribute to the multireference nature of AsH. Despite the multireference character, the 4c-CCSD(T) with perturbative triples still prevails over 4c-TCCSD. This shows that at this range of internuclear distances, the multireference character is not strong enough to cause CCSD(T) to fail. Nevertheless, in AsH, it is strong enough that TCCSD provides a major improvement over CCSD at the same computational scaling.

Calculated equilibrium energies for SbH are listed in Table VII. Compared with the higher order methods, the 4c-MP2 method outputs lower energy, which might be due to a failure to describe the multireference character of this system. As with the previous systems, the TCCSD energy is between CCSD and CCSD(T). However, in constrast, with TlH, where static correlation was not essential for a good description, here it is not clear if CCSD(T) is a good benchmark, since for multireference systems, CCSD(T) tends to output too low energy. Unfortunately, we miss a more accurate calculation to compare with; hence, we again turn to spectroscopic constants, which are listed in Table VIII, with the respective potential curve in Fig. 6. Considering the calculated internuclear distance, MP2 and DMRG are inaccurate, as they were for AsH. Oddly enough, plain DMRG with a small active space is in both cases very accurate for ω_e and ω_ex_e. Methods from CC theory succeeded in describing r_e, ω_e, and ω_ex_e, and TCCSD again outputs values between CCSD and CCSD(T). Compared with CCSD, TCCSD improved the vibrational constant ω_e.

We have implemented the relativistic tailored coupled clusters method, which is capable of treating relativistic, strongly correlated systems both in terms of static and dynamical correlation. The aim was to show that compared with the previously published calculations, we can obtain results of equal quality with much smaller active space, i.e., at a fraction of computational cost. The results presented are promising. Even with a small active space, the new method showed comparable performance for TlH to DMRG with large CAS(14, 47). The optimal CAS size related to the energy minimum was determined with low cost DMRG calculations. The calculated spectroscopic properties of TlH are in agreement with the experimental values within the error bounds. The comparison with experimental spectroscopic constants for AsH, which has a stronger multireference character, has shown that TCCSD is able to describe such systems more accurately that CCSD, with a computational cost lower than CCSD(T).

We thank the developers of the Dirac program, in particular, Dr. Lucas Visscher and Dr. Stefan Knecht for providing access to the development version of the code and for helpful discussions. The work of the Czech team has been supported by the Czech Science Foundation (Grant No. 18-24563S). Ö. Legeza has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) through Grant No. K120569 and by the Hungarian Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017-00001). The development of the relativistic DMRG libraries was supported by the Center for Scalable and Predictive methods for Excitation and Correlated phenomena (SPEC), which is funded as part of the Computational Chemical Sciences Program by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences at Pacific Northwest National Laboratory. M. Máté was supported by the ÚNKP-19-3 New National Excellence Program of the Ministry for Innovation and Technology. Mutual visits with the Hungarian group have partly been supported by the Hungarian–Czech Joint Research (Project No. MTA/19/04). Part of the CPU time for the numerical computations was supported by the Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project “IT4Innovations National Supercomputing Center – LM2015070.” Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum provided under the program “Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042) is appreciated.