Communication: Four-component density matrix renormalization group

We present the first implementation of the relativistic quantum chemical two- and four-component density matrix renormalization group algorithm that includes a variational description of scalar-relativistic effects and spin–orbit coupling. Numerical results based on the four-component Dirac–Coulomb Hamiltonian are presented for the standard reference molecule for correlated relativistic benchmarks: thallium hydride.

Owing to remarkable advances in the past decades, relativistic quantum chemical methods have become a routinely applicable and indispensable tool for the accurate description of the chemistry and spectroscopy of heavy-element compounds^1–3 and even of first- and second-row molecules.⁴ Major challenges for relativistic quantum chemistry originate from (i) the reduction of non-relativistic (spin and spatial) symmetries caused by magnetic couplings that lead to in general complex wave functions and require the use of double-group symmetry as well as (ii) the large number of (unpaired) valence electrons to be correlated (in particular for heavy elements) and (iii) the occurrence of near-degeneracies of electronic states. Popular quantum chemical methods such as CASSCF/CASPT2/SO-CASPT2⁵ assume an additivity of electron correlation and spin–orbit effects or a weak polarization of orbitals due to spin–orbit interaction, or both. Hence, for heavy-element compounds accuracy is inevitably limited as relativistic effects and static or dynamic electron correlation are often not only large but also counteracting.^1,6

To address the latter issue adequately, a number of genuine relativistic multiconfigurational and multireference approaches have been proposed.^7–9 In this Communication we merge the strengths of the density matrix renormalization group (DMRG) algorithm,¹⁰ which has been successfully introduced to the field of non-relativistic quantum chemistry,^11–13 with a variational description of all relativistic effects in the orbital basis. This new four-component (4c) DMRG ansatz goes beyond preceding scalar-relativistic DMRG approaches^14,15 and allows us to efficiently describe first and foremost non-dynamic correlation (or strong correlations) in heavy-element complexes by means of extensive active orbital spaces which would surmount capabilities of any to-date available relativistic multiconfigurational approach.

Our relativistic DMRG implementation rests on the time-independent, first-quantized 4c-Dirac–Coulomb(–Breit) Hamiltonian¹⁶ – any two-component (2c) Hamiltonian is also directly usable – (with positive-energy projectors omitted for brevity),^2,3

where

$ĥD(i)$ is the one-electron Dirac Hamiltonian for electron i,

$\hat{g}(i,j)$

$ĝ(i,j)$ is a two-electron operator describing the interaction between electrons i and j, and V_NN is the classical nuclear repulsion energy operator. In the absence of any external magnetic field it can be shown that Eq. (1) is symmetric under time-reversal³ from which follows that a fermion four-component spinor function ϕ_i occurs in Kramers pairs

$\lbrace \phi _i,\bar{\phi }_i\rbrace$

${ϕi,ϕ¯i}$⁠. A spinor

$\bar{\phi }_i$

$ϕ¯i$ can thus be obtained from the action of the time-reversal operator

$\hat{\mathcal {K}} = -i\Sigma _y {\hat{\mathcal {K}}}_0$

$K̂=−iΣyK̂0$ on ϕ_i, that is

${\hat{\mathcal {K}}}\phi _i = \bar{\phi }_i$

$K̂ϕi=ϕ¯i$⁠. Hence, our 4c- (or 2c-)spinor basis is comprised of Kramers pairs.

In the no-pair approximation, we can formulate the resulting Hamiltonian in second-quantized and normal ordered form,

where the summation indices P, Q, R, S strictly refer to positive-energy spinors, and

$FPQ$ and

$V_{PR}^{QS} = (G_{PR}^{QS} - G_{PR}^{SQ} )$

$VPRQS=(GPRQS−GPRSQ)$ are Fock-matrix elements and antisymmetrized two-electron integrals

$G_{PR}^{QS}$

$GPRQS$⁠, respectively.

We benefit from a quaternion symmetry scheme¹⁷ that has been implemented for the binary double groups

$D2h*$ and subgroups thereof in the Dirac program package which our DMRG program is interfaced. In this scheme, point group symmetry and quaternion operator algebra are combined advantageously such that the eigenvalue equation, Eq. (1), can be solved either using real (double groups

${\rm D}_{2h}^{\ast }$

$D2h*$⁠,

${\rm D}_{2}^{\ast }$

$D2*$⁠, and

${\rm C}_{2v}^{\ast }$

$C2v*$⁠; resulting number of non-zero real matrices of a quaternion operator matrix representation: NZ = 1), complex (⁠

${\rm C}_{2h}^{\ast }$

$C2h*$⁠,

${\rm C}_{2}^{\ast }$

$C2*$⁠, and

${\rm C}_{s}^{\ast }$

$Cs*$⁠; NZ = 2), or quaternion algebra (⁠

${\rm C}_{i}^{\ast }$

$Ci*$ and

${\rm C}_{1}^{\ast }$

$C1*$⁠; NZ = 4). Working in a Kramers-paired spinor basis, one can then show that all operator matrix elements

$t_{p\bar{q}}$

$tpq¯$ of a time-symmetric one-electron operator

$\hat{t}$

$t̂$ are zero by symmetry. Furthermore, the complete set of two-electron integrals

$G_{PR}^{QS}$

$GPRQS$ of the two-electron (Coulomb) operator

$\hat{g}$

$ĝ$ in molecular spinor (MS) basis can be cast into a 4 × 3 ((NZ,3)) matrix representation (see also Appendix B.3 page 161ff of Ref. 18),

where

$R$ and

$\mathcal {I}$

$I$ denote the real and complex parts of a two-electron integral in MS representation, respectively. The number of nonzero rows for a given binary double group thus corresponds to the NZ rank as given above. Important symmetry reductions for both the one- and two-electron integrals are therefore being taken into account in our relativistic Kramers-unrestricted DMRG implementation.

In a Kramers-restricted spinor basis all one-electron matrix elements

$FPQ$ (see Eq. (2)) among barred and unbarred components will be identical while matrix elements between barred and unbarred are non-zero only in the NZ = 4 case. In contrast, a two-electron integral

$G_{PR}^{QS}$

$GPRQS$ may generally be comprised of barred and unbarred spinors. As illustrated by Eq. (3) for NZ = 1 and NZ = 2, respectively, only an even number (n_barred = 0, 2, 4) of barred spinors yields a non-vanishing two-electron integral whereas for NZ = 4 all combinations are contributing. Even though integrals can be made real-valued (NZ = 1), permutational symmetry is reduced by a factor two compared to the 8-fold permutational symmetry in the non-relativistic case since orbitals are complex in a relativistic framework.

In DMRG, electron–electron correlation is taken into account by an iterative procedure that minimizes the Rayleigh quotient corresponding to the electronic Hamiltonian

$H̃̂$ and eventually converges a full-CI-type wave function within the selected active orbital space. The full configuration Hilbert space of a finite system comprising N MSs, Λ^(N), is built from tensor product spaces of local orbital (tensor) spaces Λ_i,¹⁹ which can be written as

$\Lambda ^{(N)}=\otimes _{i=1}^N \Lambda _i$

$Λ(N)=⊗i=1NΛi$⁠. Our implementation exploits a two-dimensional local Hilbert-space representation, q = 2, where each spinor can be either empty or singly occupied. The tensor space dimension is then 2^N with N being the number of spinors.

In the two-site DMRG variant,¹⁰ that is the basis for our relativistic DMRG implementation, Λ^(N) is approximated by a tensor product space of four tensor spaces, i.e.,

$Λ DMRG (N)=Λ(l)⊗Λl+1⊗Λl+2⊗Λ(r)$⁠. The dimensions of the corresponding local left (l) and right (r) spaces are denoted as

$M_l={\rm \dim }\ \Lambda ^{(l)}$

$Ml=dimΛ(l)$ and

$M_r={\rm \dim }\ \Lambda ^{(r)}$

$Mr=dimΛ(r)$⁠, respectively. With

$q\break={\rm \dim }\ \Lambda _{l+1}={\rm \dim }\ \Lambda _{l+2}$

$q=dimΛl+1=dimΛl+2$ the resulting dimensionality of the DMRG wave function is

${\rm \dim }\ \Lambda ^{(N)}_{\rm DMRG}=q^2M_lM_r\ll q^N$

$dimΛ DMRG (N)=q2MlMr≪qN$⁠. The number of block states, M_l and M_r, required to achieve sufficient convergence can be regarded as a function of the level of entanglement among the molecular orbitals. Hence the maximum number of block states M_max = max (M_l, M_r) determines the accuracy of a DMRG calculation.²⁰ 

The success and numerical efficiency of the DMRG algorithm rely on a subsequent application of the singular value decomposition (SVD) theorem^19,21 while the performance depends on the level of entanglement encoded in the wave function.²² During an SVD step, the finite system is divided into two parts

$Λ DMRG (N)=Λ(s)⊗Λ(e)$⁠, namely into a system and an environment block. In each DMRG step, the basis states of the system block are then transformed to a new truncated basis set by a unitary transformation based on the preceding SVD.²³ This transformation depends therefore on how accurately the environment is represented²⁴ as well as on the level of truncation.²⁰ As a consequence the accuracy of the DMRG method is governed by the truncation error, δɛ_TR, as well as by the environmental error, δɛ_sweep.²⁵ The latter is minimized in each DMRG macro-iteration by a successive application of the SVD going through the system back and forth (“sweeping”).

In order to minimize δɛ_sweep, which is usually largest during the initial sweep of the DMRG approach because of a poor representation of the environment, we take advantage of the Configuration Interaction based Extended Active Space procedure (CI-DEAS)^26,27 to efficiently construct the environmental basis states by means of an orbital entropy profile.²⁸ The latter is dependent on the orbital ordering along a (fictitious) one-dimensional chain^22,29 and determines the maximum number of block states M_max = max (M) that is needed to satisfy an a priori defined accuracy threshold given by a value χ.

The truncation error δɛ_TR is a function of the total number of block states M. Assuming M_l = M_r = M we can exploit a second-order polynomial fit as a function of 1/M by taking the limit of zero energy change between two sweeps E_el(M, δɛ_sweep = 0) for a given M to provide a good estimate for the truncation-free solution.^25,30

We demonstrate the capabilities of our 4c-DMRG implementation at the example of the thallium hydride molecule since this system has become a standard benchmark molecule for a plethora of relativistic methods.^31–40,58 Spinors and MS integrals were computed with a development version of the Dirac12 program package⁴¹ using the Dirac–Coulomb (DC) Hamiltonian and triple-ζ basis sets for Tl (cv3z)^42,43 and H (cc-pVTZ),⁴⁴ which include core-correlating functions for Tl. All DMRG calculations were performed with the relativistic development branch of the QC-DMRG-BUDAPEST program.⁴⁵ 

$C2v*$ double group symmetry (NZ = 1) was assumed throughout all calculations for TlH. MP2 natural spinors (NSs),⁴⁶ correlating the Tl 5s5p4f5d6s6p and H 1s electrons while keeping the remaining core electrons of Tl frozen, served as the orbital basis for all electron-correlation calculations. Since Dirac12 requires to use uncontracted basis sets in a four-component framework, a virtual orbital threshold was set at 135 hartree, such that the initial virtual correlation space in the MP2 calculation comprised all recommended core-valence and valence-correlation functions. The final active space was then chosen to include all occupied spinors that have MP2-NS occupancies less than 1.98 as well as all virtuals up to a cutoff of ≈0.001 in the MP2-NS occupation numbers. Given this criterion, an active space of 14 electrons – the occupied Tl 5d6s6p plus H 1s shells – in 47 Kramers pairs (94 spinors) was used in the CI^47–49 MP2, CC,^50,51 and DMRG calculations. The latter are further characterized by the choice of

$\smash{M_{\rm max}, M_{\rm min}, M_{\rm min}^{\rm DEAS}}$

$M max ,M min ,M min DEAS$⁠, and χ, denoted in the following as DMRG(14,94)[

$\smash{M_{\rm max},M_{\rm min},M_{\rm min}^{\rm DEAS},\chi}$

$M max ,M min ,M min DEAS ,χ$].

Figure 1 depicts the one-orbital s_i and two-orbital I_ij entropy profiles^22,52,53 at the experimental internuclear distance

$re exp$ = 1.872 Å computed from an initial DMRG[256,256,256,10⁻⁵] calculation. We first note that the one-orbital entropy profile (left-hand side of Figure 1) is nearly perfectly symmetric with respect to the unbarred (#1–#47) and barred (#48–#94) spinors where any slight deviation is an artefact of the preset low M_min, M_max values. The total quantum information I_tot encoded in the wave function, defined as the sum of one-orbital entropies, I_tot = ∑_is_i, can be taken as a measure of the importance of dynamic (weak) electron correlation. The lower I_tot (compared to

$I_{\rm tot}^{\rm max}=\sum _i s_i^{\rm max}= N\ln (2)=65.15$

$I tot max =∑isi max =Nln(2)=65.15$⁠), the more important will be an appropriate account of dynamic electron correlation in order to grasp all important correlation effects. In the present case of TlH we have

$I_{\rm tot}\simeq 2.23\ll I_{\rm tot}^{\rm max}$

$I tot ≃2.23≪I tot max$ which points to the fact that TlH is a predominantly single-reference close to its equilibrium structure.

The two-orbital mutual information, I_ij, confirms this qualitative picture. I_ij values are visualized in the right panel of Figure 1, where the degree of entanglement between spinors is marked by a color-coded connecting line. While few spinors are weakly entangled (red) the majority is entangled with even smaller strengths (green). Since several spinors are mutually entangled with the same order of magnitude, we expect that large M_min, M_max values combined with a low quantum information loss threshold χ are required to reach a fully converged DMRG wave function.

To corroborate this hypothesis we compiled in Table I total energy differences for various standard wave-function-expansion methods as well as for our 4c-DMRG(14,94)[4500,1024,2048,10⁻⁵] model with respect to a chosen 4c-CCSDTQ reference at

$re exp$ = 1.872 Å. The 4c-DMRG wave function was built from an optimized ordering of orbitals based on the entropy profiles given in Figure 1 and by applying high accuracy settings in the initial CI-DEAS sweep (with CI_level = 4 and χ_CI = 10⁻⁸). These initial conditions ensured both a rapid elimination of the environmental error and a fast total convergence towards the global minimum. The 4c-CISDTQ energy is reached after no more than six sweeps of the 4c-DMRG wave function optimization procedure. Inspection of Table I furthermore reveals that the 4c-DMRG energy is, although being below our best variational 4c-CISDTQ energy, still 2.57 mH higher than the reference 4c-CCSDTQ as well as 2.89 mH higher than the single-reference 4c-CCSD(T) energies. We recall that DMRG is best suited for static-correlation problems, whereas TlH is dominated by dynamic correlation. The missing dynamic-correlation contributions may be captured by multireference perturbation theory on top of the relativistic DMRG wave function. However, extrapolating the DMRG energy for a given M to the limit E_el(M, δɛ_sweep = 0) provides an effective means to eliminate the truncation error. The resulting best estimate of E_el(M → ∞) = −20275.8395 H at

$r_e^{\rm exp}$

$re exp$=1.872 Å is then as close as +0.7 mH to the CCSDTQ reference energy. Moreover, we found that even a reduced-M 4c-DMRG[512, δɛ_sweep = 0] potential energy curve does not only effectively reproduce the shape of the 4c-CCSDTQ potential energy curve but also yields accurate spectroscopic constants—extracted from a fourth-order polynomial fit and compiled in Table II. The 4c-CCSDTQ data are in excellent agreement with experiment for the equilibrium internuclear distance r_e, harmonic frequency ω_e, and for the anharmonicity constant ω_ex_e. Our DMRG results also show an excellent agreement with experiment for r_e, while predicting slightly too high values for ω_e (+20 cm⁻¹) and ω_ex_e (+4 cm⁻¹), respectively.

We conclude with a note on the computational demands of our 4c-DMRG implementation in comparison to CCSDTQ. The benchmark DMRG[4500,1024,2048,10⁻⁵] calculation required ≈50 GB of core memory to represent all operators of the left and right blocks of the most expensive local optimization step while the relativistic MRCC code⁵⁴ had similar memory requirements for the optimization of the various t amplitudes.

Our new 4c- and 2c-DMRG approach (no 2c-results shown here) bears the potential to become a new powerful tool for the theoretical chemistry and photophysics of heavy-element molecules dominated by strong static electron correlation as observed, for instance, in lanthanide and actinide complexes.

M.R. has been financially supported by SNF (No. 200020_144458/1) and O.L. by OTKA (Nos. K100908 and NN110360).