There are three essential problems in computational relativistic chemistry: Electrons moving at relativistic speeds, close lying states, and dynamical correlation. Currently available quantum-chemical methods are capable of solving systems with one or two of these issues. However, there is a significant class of molecules in which all the three effects are present. These are the heavier transition metal compounds, lanthanides, and actinides with open d or f shells. For such systems, sufficiently accurate numerical methods are not available, which hinders the application of theoretical chemistry in this field. In this paper, we combine two numerical methods in order to address this challenging class of molecules. These are the relativistic versions of coupled cluster methods and the density matrix renormalization group (DMRG) method. To the best of our knowledge, this is the first relativistic implementation of the coupled cluster method externally corrected by DMRG. The method brings a significant reduction of computational costs as we demonstrate on the system of TlH, AsH, and SbH.

At the turn of the millennium, the density matrix renormalization group (DMRG) method1 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,7 internally contracted MRCI (multireference configuration interaction),8 random phase approximation,9 canonical transformation method,10 or the perturbation theory with matrix product states (MPS).11 

Our group has followed a different pathway to deal with the dynamical correlation: the coupled cluster (CC) 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 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.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.

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.15 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 Ms = 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 Ms = 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.

FIG. 1.

Spin–orbit splitting in AsH and SbH.

FIG. 1.

Spin–orbit splitting in AsH and SbH.

Close modal

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,

(1)

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,

(2)

The Kramers symmetry replaces the spin symmetry in the non-relativistic theory; in particular, MS is not a good quantum number and MK 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 MK. Since each creation or annihilation operator in (1) changes MK by ±1/2, the Hamiltonian couples states with |ΔMK| ≤ 2. Moreover, the index permutation symmetry of the 2e-integrals in (1) is lower than in the non-relativistic case.

The Dirac program16 employs a quaternion symmetry approach that combines the Kramers and binary double group symmetry (D2h* and subgroups).17 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*).18 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.20 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

(3)

where Tcas 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 Text 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.26 The simplest version of the method truncates both Tcas and Text to single and double excitations. Since there is a single-determinantal Fermi vacuum, the excitation operators Text and Tcas commute, which keeps the method very simple. TCC can thus use the standard CCSD solver, modified to keep the amplitudes from Tcas fixed. Thanks to the two-body Hamiltonian, tailored CCSD energy with Text = 0 and Tcas 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.27 

In Refs. 12 and 28, we have described how to obtain Tcas from the DMRG wave function using concepts of quantum information theory29 in the non-relativistic case, yielding the DMRG–TCCSD method. The DMRG method30 is a procedure that variationally optimizes the wave function in the form of the matrix product state (MPS) ansatz.5 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,5 

(4)

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 (Ml 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 (Mr 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) procedure4 assuring the same accuracy along the sweep. Alternatively, one can use the one-site approach in the last sweep.39 

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

(5)
(6)

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 ΔMK 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 ΔMK. 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.41 The Dirac program package16 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.16 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,4 keeping up to thousands of block states for the a priori set quantum information loss threshold χ = 10−6.

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 = ErelEnonrel, 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.

TABLE I.

Comparison of energies of the HF molecule obtained from relativistic and non-relativistic methods. The rightmost column shows the difference between the 4c-relativistic energy Erel and non-relativistic energy Enonrel.

MethodErel (Eh)Enonrel (Eh)ΔE (Eh)
RHF −100.149 72 −100.058 46 −0.091 26 
DMRG(6,6) −100.158 68 −100.067 37 −0.091 30 
CCSD −100.423 22 −100.331 72 −0.091 50 
TCCSD(6,6) −100.424 18 −100.332 46 −0.091 72 
MethodErel (Eh)Enonrel (Eh)ΔE (Eh)
RHF −100.149 72 −100.058 46 −0.091 26 
DMRG(6,6) −100.158 68 −100.067 37 −0.091 30 
CCSD −100.423 22 −100.331 72 −0.091 50 
TCCSD(6,6) −100.424 18 −100.332 46 −0.091 72 

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 Eh. We used MP2 natural spinors (NS) from the Dirac program16 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 Eh. 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.

FIG. 2.

Schematic depiction of active spaces used in the 4-TCCSD procedure for TlH.

FIG. 2.

Schematic depiction of active spaces used in the 4-TCCSD procedure for TlH.

Close modal

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 Eh for AsH and −6481.107 75 Eh 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 Eh to +0.25 Eh and +1.01 Eh 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 calculation15 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 mEH. 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).

TABLE II.

Total electronic energy and energy differences ΔEel (in mEh) for various methods with respect to the 4c-CCSDTQ(14, 47) reference energy of −20 275.840 242 33 Eh15 for TlH at the experimental equilibrium internuclear distance 1.872 Å.

MethodEel (Eh)ΔEel (mEh)
4c-MP2(14, 47) −20 275.853 72 −13.49 
4c-CCSD(14, 47) −20 275.829 66 10.58 
4c-CCSD(T) (14, 47) −20 275.840 56 −0.32 
4c-DMRG(14, 47)a,15  −20 275.837 67 2.57 
4c-TCCSD(14, 10) −20 275.830 42 9.83 
4c-TCCSD(14, 11) −20 275.831 70 8.54 
4c-TCCSD(14, 12) −20 275.832 57 7.67 
4c-TCCSD(14, 13) −20 275.833 29 6.95 
4c-TCCSD(14, 14) −20 275.834 30 5.94 
4c-TCCSD(14, 15) −20 275.832 24 8.00 
4c-TCCSD(14, 16) −20 275.829 02 11.23 
4c-TCCSD(14, 17) −20 275.824 05 16.19 
MethodEel (Eh)ΔEel (mEh)
4c-MP2(14, 47) −20 275.853 72 −13.49 
4c-CCSD(14, 47) −20 275.829 66 10.58 
4c-CCSD(T) (14, 47) −20 275.840 56 −0.32 
4c-DMRG(14, 47)a,15  −20 275.837 67 2.57 
4c-TCCSD(14, 10) −20 275.830 42 9.83 
4c-TCCSD(14, 11) −20 275.831 70 8.54 
4c-TCCSD(14, 12) −20 275.832 57 7.67 
4c-TCCSD(14, 13) −20 275.833 29 6.95 
4c-TCCSD(14, 14) −20 275.834 30 5.94 
4c-TCCSD(14, 15) −20 275.832 24 8.00 
4c-TCCSD(14, 16) −20 275.829 02 11.23 
4c-TCCSD(14, 17) −20 275.824 05 16.19 
a

4c-DMRG(14, 47)[4500, 1024, 2048, 10−5] (see Ref. 15).

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.42 The obvious reason is that the “frozen” TCAS 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.42 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.42 

FIG. 3.

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. (a) Comparison of TCCSD and DMRG methods. The horizontal solid line represents the “FCI-limit” from the large 4c-DMRG(14, 47)15 calculation. (b) Details of 4c-TCCSD energies.

FIG. 3.

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. (a) Comparison of TCCSD and DMRG methods. The horizontal solid line represents the “FCI-limit” from the large 4c-DMRG(14, 47)15 calculation. (b) Details of 4c-TCCSD energies.

Close modal

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).

TABLE III.

Spectroscopic constants of 205TlH obtained from 4c-TCCSD, compared with calculations and experimental work from the literature. The spectroscopic constants have been evaluated from the potential energy curve fit with two different methodologies. In the case of the TWOFIT methodology, the number of points has been selected according to mean displacement in the harmonic ground state criterion. In the case of VIBANAL methodology, a wider symmetric interval around equilibrium geometry has been selected. In all cases, internuclear separation axis sampling was chosen to be 0.02 Å.

Methodre (Å)ωe (cm−1)ωexe (cm−1)
Expt.a,43  1.872 1391 22.7 
4c-DMRG(14, 47)15  1.873 1411 26.6 
4c-CCSD(14, 47)15  1.871 1405 19.4 
4c-TCCSD(14, 10)b DBSS 1.874 1404 24.6 
4c-TCCSD(14, 10)b M = 512 1.874 1403 23.4 
4c-TCCSD(14, 14)b DBSS 1.869 1412 22.6 
4c-TCCSD(14, 14)c DBSS 1.869 1411 20.1 
4c-TCCSD(14, 14)b M = 512 1.869 1411 22.6 
4c-TCCSD(14, 14)c M = 512 1.869 1411 20.3 
4c-TCCSD(14, 17)b DBSS 1.859 1426 20.5 
4c-TCCSD(14, 17)c DBSS 1.859 1426 17.5 
4c-TCCSD(14, 17)b M = 512 1.859 1428 22.4 
4c-TCCSD(14, 17)c M = 512 1.859 1428 29.8 
Methodre (Å)ωe (cm−1)ωexe (cm−1)
Expt.a,43  1.872 1391 22.7 
4c-DMRG(14, 47)15  1.873 1411 26.6 
4c-CCSD(14, 47)15  1.871 1405 19.4 
4c-TCCSD(14, 10)b DBSS 1.874 1404 24.6 
4c-TCCSD(14, 10)b M = 512 1.874 1403 23.4 
4c-TCCSD(14, 14)b DBSS 1.869 1412 22.6 
4c-TCCSD(14, 14)c DBSS 1.869 1411 20.1 
4c-TCCSD(14, 14)b M = 512 1.869 1411 22.6 
4c-TCCSD(14, 14)c M = 512 1.869 1411 20.3 
4c-TCCSD(14, 17)b DBSS 1.859 1426 20.5 
4c-TCCSD(14, 17)c DBSS 1.859 1426 17.5 
4c-TCCSD(14, 17)b M = 512 1.859 1428 22.4 
4c-TCCSD(14, 17)c M = 512 1.859 1428 29.8 
a

GRECP spin–orbit MRD-CI (see Ref. 43).

b

TWOFIT fourth order polynomial.

c

VIBANAL 10th order polynomial. Rmin–Rmax (Å): 1.64–2.20 for (14, 14), 1.70–2.04 for (14, 17) DBSS, and 1.72–2.00 for (14, 17) M = 512.

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.

FIG. 4.

Dissociation curve of TlH.

FIG. 4.

Dissociation curve of TlH.

Close modal

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.

TABLE IV.

Three configurations with the highest coefficients (in absolute values) for TlH, SbH, and AsH in the equilibrium internuclear distance, as generated by DMRG with the active space of 14 spinors. In each case, the CC amplitudes were generated with respect to the closed shell reference listed in the first row for given system. Here, 2 is for a doubly occupied Kramers pair and 0 is for an empty Kramers pair.

DeterminantCoefficient
TlH 2 2 2 2 2 2 2 0 0 0 0 0 0 0 Reference 
 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0.978 62 
 2 2 2 2 2 2 0 2 0 0 0 0 0 0 0.059 28 
 2 2 2 2 2 2 0 0 2 0 0 0 0 0 0.055 40 
SbH 2 2 2 2 2 2 2 2 0 0 0 0 0 0 Reference 
 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0.825 18 
 2 2 2 2 2 2 2 0 2 0 0 0 0 0 0.547 48 
 2 2 2 2 2 2 0 2 0 2 0 0 0 0 0.048 72 
AsH 2 2 2 2 2 2 2 2 0 0 0 0 0 0 Reference 
 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0.756 73 
 2 2 2 2 2 2 2 0 2 0 0 0 0 0 0.640 48 
 2 2 2 2 2 2 0 2 0 2 0 0 0 0 0.041 23 
DeterminantCoefficient
TlH 2 2 2 2 2 2 2 0 0 0 0 0 0 0 Reference 
 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0.978 62 
 2 2 2 2 2 2 0 2 0 0 0 0 0 0 0.059 28 
 2 2 2 2 2 2 0 0 2 0 0 0 0 0 0.055 40 
SbH 2 2 2 2 2 2 2 2 0 0 0 0 0 0 Reference 
 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0.825 18 
 2 2 2 2 2 2 2 0 2 0 0 0 0 0 0.547 48 
 2 2 2 2 2 2 0 2 0 2 0 0 0 0 0.048 72 
AsH 2 2 2 2 2 2 2 2 0 0 0 0 0 0 Reference 
 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0.756 73 
 2 2 2 2 2 2 2 0 2 0 0 0 0 0 0.640 48 
 2 2 2 2 2 2 0 2 0 2 0 0 0 0 0.041 23 

The equilibrium energies of AsH in Table V show that 4c-TCCSD improved the 4c-CCSD by 17 mEh, with just 3 mEh 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 re 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.

TABLE V.

Total energy for various methods for AsH at the experimental equilibrium internuclear distance45 of 1.5343 Å.

MethodEel (Eh)
4c-MP2(16, 81) −2260.532 27 
4c-CCSD(16, 81) −2260.532 41 
4c-TCCSD(16, 14) −2260.549 45 
4c-CCSD(T)(16, 81) −2260.551 33 
MethodEel (Eh)
4c-MP2(16, 81) −2260.532 27 
4c-CCSD(16, 81) −2260.532 41 
4c-TCCSD(16, 14) −2260.549 45 
4c-CCSD(T)(16, 81) −2260.551 33 
TABLE VI.

Spectroscopic constants of 75AsH obtained from 4c-TCCSD, compared with calculations and experimental work from the literature.

Methodre (Å)ωe (cm−1)ωexe (cm−1)
Expt.44  1.523 2156 39.2 
4c-SCFa 1.513 2382 33.6 
4c-MP2(16, 81)a 1.503 2256 35.0 
4c-DMRG(16, 14)a DBSS 1.545 2051 42.2 
4c-CCSD(16, 81)a 1.505 2281 40.0 
4c-TCCSD(16, 14)b DBSS 1.517 2154 43.9 
4c-TCCSD(16, 14)a DBSS 1.517 2172 39.6 
4c-CCSD(T) (16, 81)a 1.521 2146 40.6 
Methodre (Å)ωe (cm−1)ωexe (cm−1)
Expt.44  1.523 2156 39.2 
4c-SCFa 1.513 2382 33.6 
4c-MP2(16, 81)a 1.503 2256 35.0 
4c-DMRG(16, 14)a DBSS 1.545 2051 42.2 
4c-CCSD(16, 81)a 1.505 2281 40.0 
4c-TCCSD(16, 14)b DBSS 1.517 2154 43.9 
4c-TCCSD(16, 14)a DBSS 1.517 2172 39.6 
4c-CCSD(T) (16, 81)a 1.521 2146 40.6 
a

TWOFIT fourth order polynomial.

b

VIBANAL 10th order polynomial.

FIG. 5.

Detail of the potential curve of AsH near the equilibrium internuclear distance.

FIG. 5.

Detail of the potential curve of AsH near the equilibrium internuclear distance.

Close modal

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 ωexe. Methods from CC theory succeeded in describing re, ωe, and ωexe, and TCCSD again outputs values between CCSD and CCSD(T). Compared with CCSD, TCCSD improved the vibrational constant ωe.

TABLE VII.

Total energy for various methods for SbH at the equilibrium internuclear distance46 of 1.701 87 Å.

MethodEel (Eh)
4c-MP2(16, 81) −6481.696 51 
4c-CCSD(16, 81) −6481.671 08 
4c-TCCSD(16, 14) −6481.683 80 
4c-CCSD(T)(16, 81) −6481.693 15 
MethodEel (Eh)
4c-MP2(16, 81) −6481.696 51 
4c-CCSD(16, 81) −6481.671 08 
4c-TCCSD(16, 14) −6481.683 80 
4c-CCSD(T)(16, 81) −6481.693 15 
TABLE VIII.

Spectroscopic constants of 121SbH obtained from 4c-TCCSD compared with calculations and experimental work from the literature.

Methodre (Å)ωe (cm−1)ωexe (cm−1)
Expt.46  1.702   
Expt.44  1.711 1897  
Expt.47   1923 34.2 
4c-SCFa 1.705 2024 28.0 
4c-SCFb 1.704 2043 38.1 
4c-MP2(16, 81)a 1.693 2009 30.2 
4c-DMRG(16, 14)a M = 2200 1.737 1839 35.3 
4c-CCSD(16, 81)a 1.706 1945 32.6 
4c-TCCSD(16, 14)a M = 2200 1.706 1937 36.4 
4c-CCSD(T)(16, 81)a 1.710 1916 35.3 
Methodre (Å)ωe (cm−1)ωexe (cm−1)
Expt.46  1.702   
Expt.44  1.711 1897  
Expt.47   1923 34.2 
4c-SCFa 1.705 2024 28.0 
4c-SCFb 1.704 2043 38.1 
4c-MP2(16, 81)a 1.693 2009 30.2 
4c-DMRG(16, 14)a M = 2200 1.737 1839 35.3 
4c-CCSD(16, 81)a 1.706 1945 32.6 
4c-TCCSD(16, 14)a M = 2200 1.706 1937 36.4 
4c-CCSD(T)(16, 81)a 1.710 1916 35.3 
a

TWOFIT fourth order polynomial.

b

VIBANAL 10th order polynomial.

FIG. 6.

Detail of the potential curve of SbH near the equilibrium internuclear distance.

FIG. 6.

Detail of the potential curve of SbH near the equilibrium internuclear distance.

Close modal

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.

1.
2.
S. R.
White
and
R. L.
Martin
,
J. Chem. Phys.
110
,
4127
(
1999
).
3.
G. K.-L.
Chan
and
M.
Head-Gordon
,
J. Chem. Phys.
116
,
4462
(
2002
).
4.
Ö.
Legeza
,
J.
Röder
, and
B. A.
Hess
,
Phys. Rev. B
67
,
125114
(
2003
).
5.
6.
R.
Olivares-Amaya
,
W.
Hu
,
N.
Nakatani
,
S.
Sharma
,
J.
Yang
, and
G. K.-L.
Chan
,
J. Chem. Phys.
142
,
034102
(
2015
).
7.
Y.
Kurashige
and
T.
Yanai
,
J. Chem. Phys.
135
,
094104
(
2011
).
8.
M.
Saitow
,
Y.
Kurashige
, and
T.
Yanai
,
J. Chem. Phys.
139
,
044118
(
2013
).
9.
S.
Wouters
,
N.
Nakatani
,
D.
Van Neck
, and
G. K.-L.
Chan
,
Phys. Rev. B
88
,
075122
(
2013
).
10.
T.
Yanai
and
G. K.-L.
Chan
,
J. Chem. Phys.
124
,
194106
(
2006
).
11.
J.
Ren
,
Y.
Yi
, and
Z.
Shuai
,
J. Chem. Theory Comput.
12
,
4871
(
2016
).
12.
L.
Veis
,
A.
Antalík
,
J.
Brabec
,
F.
Neese
,
Ö.
Legeza
, and
J.
Pittner
,
J. Phys. Chem. Lett.
7
,
4072
(
2016
).
13.
T.
Kinoshita
,
O.
Hino
, and
R. J.
Bartlett
,
J. Chem. Phys.
123
,
074106
(
2005
).
14.
O.
Hino
,
T.
Kinoshita
,
G. K.-L.
Chan
, and
R. J.
Bartlett
,
J. Chem. Phys.
124
,
114311
(
2006
).
15.
S.
Knecht
,
Ö.
Legeza
, and
M.
Reiher
,
J. Chem. Phys.
140
,
041101
(
2014
).
16.
T.
Saue
 et al., DIRAC, A relativistic ab initio electronic structure program (
2018
), http://www.diracprogram.org.
17.
T.
Saue
and
H. J. A.
Jensen
,
J. Chem. Phys.
111
,
6211
(
1999
).
18.
K. G.
Dyall
and
K.
Fægri
, Jr.
,
Introduction to Relativistic Quantum Chemistry
(
Oxford University Press
,
2007
).
19.
J.
Thyssen
, “
Development and applications of methods for correlated relativistic calculations of molecular properties
,” Ph.D. thesis,
Univeristy of Southern Denmark
,
2001
.
20.
X.
Li
and
J.
Paldus
,
J. Comput. Phys.
107
,
6257
(
1997
).
21.
D. I.
Lyakh
,
V. F.
Lotrich
, and
R. J.
Bartlett
,
Chem. Phys. Lett.
501
,
166
(
2011
).
22.
A.
Melnichuk
and
R. J.
Bartlett
,
J. Comput. Phys.
137
,
214103
(
2012
).
23.
A.
Melnichuk
and
R. J.
Bartlett
,
J. Comput. Phys.
140
,
064113
(
2014
).
24.
P.
Piecuch
,
N.
Oliphant
, and
L.
Adamowicz
,
J. Chem. Phys.
99
,
1875
(
1993
).
25.
P.
Piecuch
and
L.
Adamowicz
,
J. Chem. Phys.
100
,
5792
(
1994
).
26.
K.
Kowalski
,
J. Chem. Phys.
148
,
094104
(
2018
).
27.
F. M.
Faulstich
,
A.
Laestadius
,
Ö.
Legeza
,
R.
Schneider
, and
S.
Kvaal
,
SIAM J. Numer. Anal.
57
,
2579
(
2019
).
28.
L.
Veis
,
A.
Antalík
,
J.
Brabec
,
F.
Neese
,
Ö.
Legeza
, and
J.
Pittner
,
J. Phys. Chem. Lett.
8
,
291
(
2017
).
29.
Ö.
Legeza
and
J.
Sólyom
,
Phys. Rev. B
68
,
195116
(
2003
).
30.
U.
Schollwöck
,
Rev. Mod. Phys.
77
,
259
(
2005
).
31.
Ö.
Legeza
,
R.
Noack
,
J.
Sólyom
, and
L.
Tincani
, in
Computational Many-Particle Physics
, Lecture Notes in Physics Vol. 739, edited by
H.
Fehske
,
R.
Schneider
, and
A.
Weisse
(
Springer
Berlin Heidelberg
,
2008
), pp.
653
664
.
32.
K. H.
Marti
and
M.
Reiher
,
Z. Phys. Chem.
224
,
583
(
2010
).
33.
G. K.-L.
Chan
and
S.
Sharma
,
Ann. Rev. Phys. Chem.
62
,
465
(
2011
).
34.
S.
Wouters
and
D.
Van Neck
,
Eur. Phys. J. D
68
,
272
(
2014
).
35.
S.
Szalay
,
M.
Pfeffer
,
V.
Murg
,
G.
Barcza
,
F.
Verstraete
,
R.
Schneider
, and
Ö.
Legeza
,
Int. J. Quantum Chem.
115
,
1342
(
2015
).
36.
T.
Yanai
,
Y.
Kurashige
,
W.
Mizukami
,
J.
Chalupský
,
T. N.
Lan
, and
M.
Saitow
,
Int. J. Quantum Chem.
115
,
283
(
2015
).
37.
G.
Moritz
and
M.
Reiher
,
J. Chem. Phys.
126
,
244109
(
2007
).
38.
K.
Boguslawski
,
K. H.
Marti
, and
M.
Reiher
,
J. Chem. Phys.
134
,
224101
(
2011
).
39.
D.
Zgid
and
M.
Nooijen
,
J. Chem. Phys.
128
,
144115
(
2008
).
40.
S.
Battaglia
,
S.
Keller
, and
S.
Knecht
,
J. Chem. Theory Comput.
14
,
2353
(
2018
).
41.
Ö.
Legeza
,
L.
Veis
, and
T.
Mosoni
, QC-DMRG-Budapest, a Program for Quantum Chemical DMRG Calculations.
42.
F. M.
Faulstich
,
M.
Máté
,
A.
Laestadius
,
M. A.
Csirik
,
L.
Veis
,
A.
Antalik
,
J.
Brabec
,
R.
Schneider
,
J.
Pittner
,
S.
Kvaal
, and
Ö.
Legeza
,
J. Chem. Theory Comput.
15
,
2206
(
2019
).
43.
A. V.
Titov
,
N. S.
Mosyagin
,
A. B.
Alekseyev
, and
R. J.
Buenker
,
Int. J. Quantum Chem.
81
,
409
(
2001
).
44.
K. D.
Hensel
,
R. A.
Hughes
, and
J. M.
Brown
,
J. Chem. Soc., Faraday Trans.
91
,
2999
(
1995
).
45.
K.
Balasubramanian
,
Chem. Rev.
89
,
1801
(
1989
).
46.
M.
Beutel
,
K. D.
Setzer
,
O.
Shestakov
, and
E. H.
Fink
,
J. Mol. Spectrosc.
179
,
79
(
1996
).
47.
R. D.
Urban
,
K.
Essig
, and
H.
Jones
,
J. Chem. Phys.
99
,
1591
(
1993
).