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

DMRG-tailored coupled cluster method in the 4c-relativistic domain Jan Brandejs,1, 2, a) Jakub Višňák,1, 2, 3, b) Libor Veis,1, c) Örs Legeza,4, d) and Jiří Pittner1, e) 1)J. Heyrovský Institute of Physical Chemistry, Academy of Sciences of the Czech Republic, v.v.i., Dolejškova 3, 18223 Prague 8, Czech Republic 2)Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic 3)Czech Academic City in Erbil, Yassin Najar street, Kurani Ankawa, Erbil, Kurdistan Region of Iraq 4)Strongly Correlated Systems “Lendület” Research Group, Institute for Solid State Physics and Optics, MTA Wigner Research Centre for Physics, H-1121 Budapest, Konkoly-Thege Miklós út 29-33, Hungary


I. INTRODUCTION
At the turn of the millennium, the density matrix renormalization group method (DMRG) 1 was introduced to quantumchemical community [2][3][4] and since then, it has seen a large surge in the use for multireference systems. The biggest advantage of DMRG method is its capability to treat large active spaces, 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 in 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 pioneering stage. Past efforts were either based on second order perturbation theory 7 , internally contracted MRCI (multireference configuration interaction) 8 , random phase approximation 9 , canonical transformation method 10 , or the perturbation theory with matrix product states 11 . Our group has followed a different pathway to deal with the dynamical correlation, the coupled cluster method externally corrected by DMRG 12 . 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 a) Electronic mail: jan.brandejs@jh-inst.cas.cz b) Electronic mail: jakub.visnak@jh-inst.cas.cz c) Electronic mail: libor.veis@jh-inst.cas.cz d) Electronic mail: legeza.ors@wigner.mta.hu e) Electronic mail: jiri.pittner@jh-inst.cas.cz for the static correlation. Second step is CC analysis of 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 , yields very promising results 12 . 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.
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 results are available 14 . It should be noted that DMRG is best suited for staticcorrelation problems while TlH is dominated mostly by dynamic correlation, for which CC approaches are excellent.

II. THEORY
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 nonrelativistic case where the indices P, Q, R, S run over the positive-energy 4component 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 1/2 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 15 employs a quaternion symmetry approach which combines the Kramers and binary double group symmetry (D * 2h and subgroups) 16 . The double groups can be sorted into three classes based on the application of the Frobenius-Schur indicator to their irreducible representations: "real groups" (D * 2h , D * 2 , and C * 2v ); "complex groups" (C * 2h , C * 2 , and C * s ); and "quaternion groups" (C * i and C * 1 ) 17 . 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 even number of barred indices are nonzero. Finally, in the remaining case of "quaternion groups" all the integrals have to be included and are complex-valued 17,18 .
The idea of externally corrected coupled cluster methods is to take information on static correlation from some non-CC external source, and to include it into the subsequent CC treatment 19 . The conceptually simplest approach is the tailored CC method (TCC) proposed by Bartlett et al. 13,[20][21][22] , which uses the split-amplitude ansatz for the wave function introduced by Piecuch et al. 23,24 |Ψ = e T ext e T cas |Φ where T cas containing amplitudes with all active indices is "frozen" at values obtained from 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. The simplest version of the method truncates both T cas and T ext to single and double excitations. Since there is a singledeterminantal 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 the 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 FCI energy. In general, a quadratic error bound valid for TNS-TCC methods is derived 25 .
In 12,26 we have described how to obtain T cas from the DMRG wave function using concepts of quantum information theory 27 in the non-relativistic case, yielding the DMRG-TCCSD method. The DMRG method 28 is a procedure which variationally optimizes the wave function in the form of the matrix product state (MPS) ansatz 5 . The quantum chemical version of DMRG (QC-DMRG) [29][30][31][32][33][34] 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 5 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 i th orbital), right site (four dimensional tensor space of (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 c a i and c ab i j for a, b, i, j ∈ CAS can be efficiently calculated by contractions of MPS matrices 35,36 . 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 4 assuring the same accuracy along the sweep. Alternatively, one can use the onesite approach in the last sweep 37 .
Once the CI coefficients c a i and c ab i j 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 14,38 . Secondly, 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.

A. Comparison with non-relativistic TCC
In order to compare with non-relativistic version of TCCSD method, the system of hydrogen fluoride was chosen, as it is a biatomic with light nuclei. CC-pVTZ basis was used at the internuclear distance of 89.96 Å. 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.

B. TlH
We have used the computational protocol of Ref. 14 for direct comparison with their energies. Orbitals and MO integrals were generated with the Dirac program package 15 We used the Dirac-Coulomb Hamiltonian and triple-zeta basis sets for Tl (cv3z) and H (cc-pVTZ), which include core-correlating functions for Tl. DMRG calculations were performed with the relativistic branch of the QC-DMRG-Budapest program 39 . C 2v double group symmetry with real irreps was assumed. The 4c-RHF energy was −20275.41661 E h . We used MP2 natural spinors (NS) from the Dirac program 15 as the spinor basis for electroncorrelation 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 1 shows a scheme of embedded active spaces used in the procedure. Initialization of DMRG, i.e., optimal ordering of spinors, was set up as discussed in Ref. 33. The numerical accuracy was controlled by the dynamic block state selection approach (DBSS) 4 keeping up to thousands of block states for the a priory set quantum information loss threshold χ = 10 −6 .  Once we reproduced the MP2 and CCSD energy of TlH in equilibrium geometry from Ref. 14, we applied the 4c-TCCSD method. Obtained energies and their respective deviations from the reference CCSDTQ calculation 14 are listed in Table II. In 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 Figure 2a, TCCSD significantly improves DMRG energy towards FCI, even for the smallest CAS space. In fact, further enlarging of CAS over the size of 14 spinors is counter productive. Although the TCC must reproduce FCI energy when CAS is extended to all spinors, the TCC energy does not approach this limit monotonously 40 . The reason is that the "frozen" T CAS amplitudes cannot reflect the influence of the dynamical correlation in the external space back on the active CAS space, therefore extending CAS space first exacerbates the results. In practice, the optimal CAS size related to the energy minimum is usually independent of M FIG. 2: Equilibrium energy of TlH calculated using the 4c-TCCSD and 4c-DMRG methods with different sizes of DMRG active space, as given in Table II. and therefore can be determined with low bond dimension (M) DMRG calculations 40 . As demonstrated by the chart 2b, 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). Table III shows the obtained spectroscopic properties of TlH. 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 14-spinor space in Table II. Moreover, for 14-spinors, TCCSD with M=512 states is close to exact DBSS calculation, which indicates that the results are 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.

V. CONCLUSION
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 is 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 very promising. Even with a small active space, the new method showed comparable performance to DMRG with large CAS (14,47). The optimal CAS size related to the energy minimum was determined with low cost DMRG calculations. Calculated spectroscopic properties of TlH agree with experimental values within the error bounds.