A general non-Abelian density matrix renormalization group algorithm with application to the C₂ dimer

The density matrix renormalization group (DMRG) was first introduced over 20 years ago^1,2 and has since become an indispensable part of the toolbox of chemists and physicists. Its application in quantum chemistry began with the seminal paper by White and Martin³ and was soon followed with work by many other groups.^4–10 Recently, the DMRG has seen many important algorithmic advancements such as orbital optimization,^11,12 combination with perturbation theory,^13,14 configuration interaction,¹⁵ and canonical transformation theory¹⁶ and has also been incorporated into several major quantum chemistry packages such as Molpro,¹⁷ Q-chem,¹⁸ and Orca.¹⁹ 

Significant improvements in DMRG efficiency can be gained by utilizing non-Abelian symmetries such as SU(2) symmetry. This was first attempted by Nishino using interaction-round-a-face DMRG (IRF-DMRG).²⁰ This was later followed by McCulloch’s work^21,22 that introduced the Wigner-Eckart theorem to gain large efficiencies and the “quasi-density matrix” to ensure that only spin multiplet states are retained during renormalization. Zgid and Nooijen²³ later implemented an algorithm for targeting specific spin states with general quantum chemistry Hamiltonians using the quasi-density matrix, but did not take advantage of the Wigner-Eckart theorem, and thus did not provide a boost in efficiency. The first efficient spin adapted code for quantum chemistry was presented independently by Wouters et al.²⁴ and Sharma and Chan.²⁵ Weichselbaum²⁶ then presented a general framework for using non-Abelian symmetries in tensor networks and demonstrated how the so-called “inner” and “outer” multiplicities can be appropriately treated. Here, we present an algorithm for utilizing the most general non-Abelian symmetry found in a molecule: SU(2) × non-Abelian point group symmetry, in conjunction with a quantum chemical Hamiltonian, where the added complication of treating outer multiplicities is not present. This extends our previous work,²⁵ which we refer to for a detailed discussion of the algorithm and notation. We use our new algorithm to calculate frozen core near-exact potential curves for the ground and several excited states of the C₂ molecule at various bond lengths in a cc-pVQZ basis set,²⁷ to 10 μE_h and 0.1 mE_h accuracies, respectively. Further, excitation energies (T_e) of eight lowest excited state of the C₂ dimer are calculated and these show good agreement with the observed experimentally measured values.

The sweep algorithm, introduced by White, is commonly used to optimize the underlying matrix product state (MPS) wavefunction of DMRG. The algorithm can be sub-divided into three steps: blocking, wavefunction solution, and decimation. In the following, we take each step of the algorithm and show how it is modified to use non-Abelian symmetry. However, before this, it is appropriate to reiterate that the central concept on which the algorithm hinges is the Wigner-Eckart theorem, shown in Eq. (1), where ψ^γ_i and ϕ^β_k are states belonging to the i and k rows of irreducible representations γ and β, respectively, and $O ˆ α j$ is the tensor operator of irreducible representation α that transforms as its jth row.

This equation states, in essence, that n_γn_β matrix elements of n_α different operators can be stored using just one “reduced matrix element” (⁠$〈 ψ γ ‖ O ˆ α ‖ ϕ β 〉$⁠), as long as the Clebsch-Gordan coefficients 〈β_kα_j|γ_i〉 of the symmetry group are known. To use the Wigner-Eckart theorem, the design of the algorithm should ensure that only those reduced sets of operators and states that transform as irreducible representation of a given non-Abelian group explicitly appear.

In this step, a single site (•) is added to the left block (⁠$L$⁠) (containing k sites) to generate a resulting system block (⁠$S$⁠) containing k + 1 sites. The system block, in general, will contain O(k²) normal (or complementary) operators, each of which can be multiplied with one of the appropriate O(k²) complementary (or normal operators) on the environment block to form the Hamiltonian (although we note that an explicit matrix representation of this Hamiltonian is never constructed). The definition of these normal operators and corresponding complementary operators with which they are multiplied are given in the two columns of Table I.

Unfortunately, the usual normal and complementary operators cannot be used directly with non-Abelian symmetries because they, in general, do not transform as irreducible representations of the group. To work with tensor operators, we start with some notation. First, we denote the creation/destruction tensor operator $a I γ † / a I γ$ of a “spatial orbital” I as a set of n_γ creation/destruction operators ${ a I γ i † } / { p γ i a I γ i }$⁠, one for each “spin orbital” i of the spatial orbital I. Here, p_γi is a phase factor which, in general, is not 1 and depends on the phase convention used in the construction of Clebsch-Gordan coefficients. It ensures that the destruction operators also transform as a set of tensor operators of irreducible representation γ.^25,28 (A greek letter is used to denote the irreducible representation and the capitalized roman letter denotes the spatial orbital.) Note, here the label “spin” and “spatial” is used more generally than in the context of spin symmetry alone. A spatial orbital is the set of all orbitals that transform as different rows of an irreducible representation. For example, when one is using the SU(2) group then each spatial orbital has only 2 spin orbitals, however, when one is using the SU(2) ×D_∞h group, then a spatial orbital belonging to irrep Π_g will have 4 spin orbitals.

When two creation tensor operators $a I α †$ and $a J β †$ are multiplied, they form a new set of n_αn_β operators that do not themselves form a set of tensor operators, but which can be converted into multiplets of tensor operators using the Clebsch-Gordan coefficients as shown in Eq. (2).

Here, the intermediate $C C I J γ k [ α i β j ] =〈 α i β j | γ k 〉 a I α i † a J β j †$ is defined for convenience. When more than two tensor operators are multiplied (see Eq. (3)), then to uniquely define the final tensor operator (⁠$CC D I J K δ l$⁠) the irreducible representation of the intermediate tensor operator (θ) and the order in which the operators are multiplied (δ_l[α_iθ_m[β_jγ_k]]) need to be specified as well.

As a short-hand for the product of two tensor operators of the form specified in Eq. (2) we define the symbol ⊗^γ, which generates n_γ tensor operators $N ˆ γ$ by taking an outer product of two tensor operators $L ˆ α$ and $M ˆ β$ by using the formula given in Eq. (4), where the same notation as the one used for the intermediate (⁠$C C I J γ k [ α i β j ]$⁠) from Eq. (2) is used on the RHS of the definition.

With these notations in place, we can conveniently define a new set of normal and complementary operators as shown in Table II to replace the operators in Table I.

These operators are constructed during the blocking step using the formulae in Eq. (5) (only formulae for the construction of complementary operators are shown), where the quantity W(αβA_gγ; γα) is the Racah W-coefficient.²⁸ In these formulae, the operators on the left (⁠$L$⁠) and dot (•) blocks are multiplied using the outer product notation given in Eq. (4).

A casual look at the redefinition of the operators in Table II suggests that they provide no computational benefit because, although there are fewer tensor operators, each tensor operator hides within it a multiplet of individual operators. Further, even though all the formulae in Eq. (5) are given in terms of these tensor operators, the definition of the outer product ⊗^γ in Eq. (4) ensures that all the individual operators are involved in its evaluation. However, this is where the “magic” of spin algebra and the Wigner-Eckart theorem enters. We have already seen that due to the Wigner-Eckart theorem, we only ever need store the reduced matrix elements of tensor operators. More importantly, during the blocking step the reduced matrix elements of the tensor operators on the system block, constructed using the formulae in Eq. (5), can always be formed with knowledge only of the reduced matrix elements of the tensor operators on the left and dot block, by making use of Eq. (6). Thus, by working with the tensor operators in Table II and utilizing various predefined coefficients such as the Clebsch-Gordan coefficients, Racah’s coefficient and 9-j coefficient of the non-Abelian group of interest, very substantial computational savings are possible.

We end this section on blocking with a few words about building multiplet states in a block. The full Fock space of an individual spatial orbital I of irreducible representation α contains 2^n_α many body states. These states do not form a basis for an irreducible representation of the symmetry group and have to be transformed using a unitary operator to form a set of multiplet states. This can be done in many ways including by using projection operators²⁹ or by just acting the elements of a tensor operator on the vacuum state. The resulting number of multiplets obtained is less than 2^n_α. For example, in SU(2) symmetry, each spatial orbital gives three set of multiplets, and in SU(2) × D_∞h symmetry, one can get up to 7 sets of multiplets from 16 individual states.

The Hilbert space used to express the DMRG wavefunction of symmetry γ is formed by the tensor product of the renormalized Fock space of the system and the environment block as shown in Eq. (7). In this space, the Hamiltonian is diagonalized using the Davidson algorithm, which requires the ability to perform Hamiltonian-wavefunction multiplication (|v〉 = H|c〉). This operation can be performed just in terms of multiplets as shown in Eq. (8). Note that to reduce the scaling of the algorithm from M⁴ to M³ the explicit matrix representation of the Hamiltonian is never generated.

During the renormalization step, the reduced density matrix of the system is formed and its most significant eigenvalues are retained. When using non-Abelian symmetries in DMRG, instead of the usual reduced density matrix the “quasi-density matrix” which belongs to the fully symmetric irreducible representation is formed, as shown in Eq. (9). Diagonalizing the quasi-density matrix ensures that the multiplets belonging to different irreducible representations do not mix together during renormalization.

C₂ is an important intermediate in combustion processes and its electronic spectrum has been the subject of many experiments.^30–32 Also, due the multi-reference nature of the ground and excited states of C₂ at stretched bond lengths, it has attracted theoretical interest as a benchmark system for methods such as coupled-cluster,^33–36 full configuration interaction,³⁷ multireference perturbation theory,^38,39 multireference configuration interaction,⁴⁰ initiator full configuration interaction quantum Monte Carlo (i-FCIQMC),^41,42 model space quantum Monte Carlo (MSQMC),^43,44 correlation energy extrapolation by intrinsic scaling (CEEIS),⁴⁵ reduced-density-matrix theory,⁴⁶ and valence bond theory.⁴⁷ Here, we use our non-Abelian symmetry adapted DMRG program to calculate the excitation energies (T_e) with respect to the ground state $X 1 Σ g +$ of eight lowest excited states: $a 3 Π u , b 3 Σ g − , A 1 Π u , c 3 Σ u + , B 1 Δ g , B ′ 1 Σ g + , d 3 Π g$⁠, and C¹Π_g, which are then compared against experimental values. We have also calculated the potential energy curves of 12 low energy states—the 3 lowest energy levels each for irreducible representations $1 Σ g + , 1 Σ u + , 1 Σ g −$⁠, and $1 Σ u −$⁠. Note, the states of irreducible representation $1 Σ g −$ and $1 Σ u −$ cannot be easily targeted when only Abelian subgroups of D_∞h are used. Further, the use of the full symmetry of the problem allows us to very efficiently obtain high accuracy, which we exploit to affordably compute not just the ground but also the excited states of each of the irreducible representations.

A ground state Hartree-Fock calculation was performed on the C₂ molecule with the cc-pVQZ²⁷ basis set using the Molpro ¹⁷ quantum chemistry package. The orbitals obtained from the Hartree-Fock calculation were first transformed to make them compatible with the Clebsch-Gordan coefficients for the D_∞h point group, as given in Altmann and Herzig.⁴⁸ The core orbitals and electrons were frozen, but the remaining 8 electrons in 108 orbitals were fully correlated. First, only the ground state calculation was performed using DMRG; the resulting energy is tabulated in the second column of Table III. By fitting the DMRG energy versus the discarded weight to a linear curve, the remaining energy error was estimated to be less than 10 μE_h as shown in Table IV. The table also shows the discarded weight and calculated energies of a DMRG calculation when only the SU(2) × D_2h group is used. The data demonstrate that roughly twice the number of retained renormalized states (M) are needed when only the D_2h subgroup of the full D_∞h is utilized. We find that the cpu time per sweep remains practically unchanged. The theoretical scaling of DMRG algorithm is O(M³), although in practice for medium sized problems often O(M²) scaling is observed, which implies that we obtain between a 4-fold and a 8-fold efficiency gain by using the full non-Abelian symmetry of the problem. The calculated DMRG energy is more accurate than the i-FCIQMC⁴² and the CEEIS⁴⁵ calculations and well within their reported error bars of 0.3 mE_h.

Next, electronic energies of eight low lying states are calculating at their respective equilibrium bond lengths,³¹ to obtain their excitation energies with respect to the ground state $X 1 Σ g +$⁠. The energies are tabulated in Table V and show good agreement with the experimentally measured excitation energies, with all errors being less than 0.03 eV except $c 3 Σ u +$ and C¹Π_g states where the errors are 0.04 eV and 0.06 eV, respectively. These errors are not entirely unexpected given the fact that there are still basis set incompleteness errors even cc-pVQZ basis set as demonstrated by Booth et al.⁴⁹ using the explicitly correlated basis sets. Also, we have not take into account the core correlation, which can result in errors in the ground state energy of up to 2 mE_h as shown by Wouters et al.⁵⁰ and most likely higher errors in the calculated excited state energies. As we have recently demonstrated elsewhere,⁵¹ the remaining basis set error can be eliminated by combining efficient DMRG algorithms for large bases as described in this work, with the transcorrelated Hamiltonian approach of Shiozaki and Yanai,⁵² thus enabling solutions of the molecular Schödinger equation to spectroscopic accuracy.

Next, four separate DMRG state-averaged calculations were performed on three states each of irreducible representations $1 Σ g + , 1 Σ u + , 1 Σ g −$⁠, and $1 Σ u −$ with a maximum M (number of retained renormalized (multiplet) states) of 4000. The energies of all the states are shown in columns 3–14 of Table III. The difference between the energies of the first and the second column gives an estimate of error of our state-averaged calculations and is less than 0.1 mE_h for all the geometries. These energies are plotted in Figure 1 and show several curve crossings and avoided curve crossings.

In summary, we have presented an efficient extension of our spin-adapted DMRG algorithm to utilize general non-Abelian point group symmetries. We used the resulting implementation to calculate the ground state potential energy curve of the C₂ molecule with a cc-pVQZ basis set (and frozen core) to an unprecedented (near-exact) accuracy of 10 μE_h. Several excited state potential energy curves were also computed to very high accuracy to demonstrate the large efficiency gains in this new formulation. The excitation energies of some of these excited states show good agreement with the experimental results. These calculations may, for practical chemical purposes, be considered exact within the given basis set.

This work was supported by the U.S. National Science Foundation (NSF) through Grant No. NSF-CHE-1265277. Additional support for software development was provided through Grant No. NSF-OCI-1265278. S.S. would also like to thank Garnet Chan for many helpful discussions.