We extend our previous work [S. Sharma and G. K.-L. Chan, J. Chem. Phys. **136**, 124121 (2012)], which described a spin-adapted (SU(2) symmetry) density matrix renormalization group algorithm, to additionally utilize general non-Abelian point group symmetries. A key strength of the present formulation is that the requisite tensor operators are not hard-coded for each symmetry group, but are instead generated on the fly using the appropriate Clebsch-Gordan coefficients. This allows our single implementation to easily enable (or disable) any non-Abelian point group symmetry (including SU(2) spin symmetry). We use our implementation to compute the ground state potential energy curve of the C_{2} dimer in the cc-pVQZ basis set (with a frozen-core), corresponding to a Hilbert space dimension of 10^{12} many-body states. While our calculated energy lies within the 0.3 mE_{h} error bound of previous initiator full configuration interaction quantum Monte Carlo and correlation energy extrapolation by intrinsic scaling calculations, our estimated residual error is only 0.01 mE_{h}, much more accurate than these previous estimates. Due to the additional efficiency afforded by the algorithm, the excitation energies (*T _{e}*) of eight lowest lying excited states:

*a*

^{3}Π

_{u}, $ b 3 \Sigma g \u2212 $,

*A*

^{1}Π

_{u}, $ c 3 \Sigma u + $,

*B*

^{1}Δ

_{g}, $ B \u2032 1 \Sigma g + $,

*d*

^{3}Π

_{g}, and

*C*

^{1}Π

_{g}are calculated, which agree with experimentally derived values to better than 0.06 eV. In addition, we also compute the potential energy curves of twelve states: the three lowest levels for each of the irreducible representations $ 1 \Sigma g + $, $ 1 \Sigma u + $, $ 1 \Sigma g \u2212 $, and $ 1 \Sigma u \u2212 $, to an estimated accuracy of 0.1 mE

_{h}of the exact result in this basis.

## I. INTRODUCTION

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^{3} 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,^{15} and canonical transformation theory^{16} and has also been incorporated into several major quantum chemistry packages such as Molpro,^{17} Q-chem,^{18} and Orca.^{19}

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).^{20} 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^{23} 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.*^{24} and Sharma and Chan.^{25} Weichselbaum^{26} 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,^{25} 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_{2} molecule at various bond lengths in a cc-pVQZ basis set,^{27} to 10 *μ*E_{h} and 0.1 mE_{h} accuracies, respectively. Further, excitation energies (*T _{e}*) of eight lowest excited state of the C

_{2}dimer are calculated and these show good agreement with the observed experimentally measured values.

## II. NON-ABELIAN SYMMETRIES IN DMRG

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 \u02c6 \alpha j $ is the tensor operator of irreducible representation *α* that transforms as its *j*th row.

This equation states, in essence, that *n*_{γ}*n*_{β} matrix elements of *n*_{α} different operators can be stored using just one “reduced matrix element” ($\u3008 \psi \gamma \Vert O \u02c6 \alpha \Vert \varphi \beta \u3009$), as long as the Clebsch-Gordan coefficients 〈*β _{k}α_{j}*|

*γ*〉 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.

_{i}### A. Blocking

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*^{2}) normal (or complementary) operators, each of which can be multiplied with one of the appropriate *O*(*k*^{2}) 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.

Normal operator . | Complementary operator . |
---|---|

$ I \u02c6 $ | $ H \u02c6 = t i j a i \u2020 a j +\u3008ij|kl\u3009 a i \u2020 a j \u2020 a l a k $ |

$ a i \u2020 $ | $ R \u02c6 i = t i j a j +\u3008ij|kl\u3009 a j \u2020 a l a k $ |

$ A \u02c6 i j = a i \u2020 a j \u2020 $ | $ P \u02c6 i j =\u3008ij|kl\u3009 a l a k $ |

$ B \u02c6 i j = a i \u2020 a j $ | $ Q \u02c6 i j = ( \u3008 i k | j l \u3009 \u2212 \u3008 i k | l j \u3009 ) a k \u2020 a l $ |

Normal operator . | Complementary operator . |
---|---|

$ I \u02c6 $ | $ H \u02c6 = t i j a i \u2020 a j +\u3008ij|kl\u3009 a i \u2020 a j \u2020 a l a k $ |

$ a i \u2020 $ | $ R \u02c6 i = t i j a j +\u3008ij|kl\u3009 a j \u2020 a l a k $ |

$ A \u02c6 i j = a i \u2020 a j \u2020 $ | $ P \u02c6 i j =\u3008ij|kl\u3009 a l a k $ |

$ B \u02c6 i j = a i \u2020 a j $ | $ Q \u02c6 i j = ( \u3008 i k | j l \u3009 \u2212 \u3008 i k | l j \u3009 ) a k \u2020 a l $ |

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 \gamma \u2020 / a I \gamma $ of a “spatial orbital” *I* as a set of *n*_{γ} creation/destruction operators $ { a I \gamma i \u2020 } / { p \gamma i a I \gamma 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 \alpha \u2020 $ and $ a J \beta \u2020 $ 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 \gamma k [ \alpha i \beta j ] =\u3008 \alpha i \beta j | \gamma k \u3009 a I \alpha i \u2020 a J \beta j \u2020 $ 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 \delta l $) the irreducible representation of the intermediate tensor operator (*θ*) and the order in which the operators are multiplied (*δ _{l}*[

*α*[

_{i}θ_{m}*β*]]) need to be specified as well.

_{j}γ_{k}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 \u02c6 \gamma $ by taking an outer product of two tensor operators $ L \u02c6 \alpha $ and $ M \u02c6 \beta $ by using the formula given in Eq. (4), where the same notation as the one used for the intermediate ($C C I J \gamma k [ \alpha i \beta 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.

Normal operator . | Complementary operator . |
---|---|

$ I \u02c6 $ | $ H \u02c6 = \u2211 \alpha i \beta j t I J \alpha i \beta j a I \alpha i \u2020 a J \beta j +\u2211 \nu I J K L \alpha i \beta j \delta k \xi l a I \alpha i \u2020 a J \beta j \u2020 a L \xi l a K \delta k $ |

$ a I \alpha i \u2020 $ | $ R \u02c6 I \alpha i = \u2211 I J \alpha I t I J \alpha i \u2212 \alpha i a J \u2212 \alpha i + \u2211 J K L \beta j \delta k \xi l \theta m \nu I J K L \alpha i \beta j \delta k \xi l \u3008 \theta m | \xi l \delta k \u3009 \u3008 \u2212 \alpha i | \theta m \beta j \u3009 \u3008 \alpha i \u2212 \alpha i | A g \u3009 CD D J L K \u2212 \alpha i [ \beta j \theta m [ \xi l \delta k ] ] $ |

$ A \u02c6 I J \gamma m =C C I J \gamma m [ \alpha i \beta j ] $ | $ P \u02c6 I J \gamma m = \u2211 K L \delta k \xi l \alpha i \beta j \nu I J K L \alpha i \beta j \delta k \xi l \u3008 \u2212 \gamma m | \xi l \delta k \u3009 \u3008 \gamma m | \alpha i \beta j \u3009 \u3008 \gamma m \u2212 \gamma m | A g \u3009 D D L K \u2212 \gamma m [ \xi l \delta k ] $ |

$ B \u02c6 I J \gamma m =C D I J \gamma m [ \alpha i \beta j ] $ | $ Q \u02c6 I J \gamma m = \u2211 K L \delta k \xi l \alpha i \beta j ( \nu I K J L \alpha i \delta k \beta j \xi l \u2212 \nu I K L J \alpha i \delta k \xi l \beta j ) \u3008 \u2212 \gamma m | \xi l \delta k \u3009 \u3008 \gamma m | \alpha i \beta j \u3009 \u3008 \gamma m \u2212 \gamma m | A g \u3009 C D K L \u2212 \gamma m [ \delta k \xi l ] $ |

Normal operator . | Complementary operator . |
---|---|

$ I \u02c6 $ | $ H \u02c6 = \u2211 \alpha i \beta j t I J \alpha i \beta j a I \alpha i \u2020 a J \beta j +\u2211 \nu I J K L \alpha i \beta j \delta k \xi l a I \alpha i \u2020 a J \beta j \u2020 a L \xi l a K \delta k $ |

$ a I \alpha i \u2020 $ | $ R \u02c6 I \alpha i = \u2211 I J \alpha I t I J \alpha i \u2212 \alpha i a J \u2212 \alpha i + \u2211 J K L \beta j \delta k \xi l \theta m \nu I J K L \alpha i \beta j \delta k \xi l \u3008 \theta m | \xi l \delta k \u3009 \u3008 \u2212 \alpha i | \theta m \beta j \u3009 \u3008 \alpha i \u2212 \alpha i | A g \u3009 CD D J L K \u2212 \alpha i [ \beta j \theta m [ \xi l \delta k ] ] $ |

$ A \u02c6 I J \gamma m =C C I J \gamma m [ \alpha i \beta j ] $ | $ P \u02c6 I J \gamma m = \u2211 K L \delta k \xi l \alpha i \beta j \nu I J K L \alpha i \beta j \delta k \xi l \u3008 \u2212 \gamma m | \xi l \delta k \u3009 \u3008 \gamma m | \alpha i \beta j \u3009 \u3008 \gamma m \u2212 \gamma m | A g \u3009 D D L K \u2212 \gamma m [ \xi l \delta k ] $ |

$ B \u02c6 I J \gamma m =C D I J \gamma m [ \alpha i \beta j ] $ | $ Q \u02c6 I J \gamma m = \u2211 K L \delta k \xi l \alpha i \beta j ( \nu I K J L \alpha i \delta k \beta j \xi l \u2212 \nu I K L J \alpha i \delta k \xi l \beta j ) \u3008 \u2212 \gamma m | \xi l \delta k \u3009 \u3008 \gamma m | \alpha i \beta j \u3009 \u3008 \gamma m \u2212 \gamma m | A g \u3009 C D K L \u2212 \gamma m [ \delta k \xi l ] $ |

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.

^{28}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^{29} 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.

### B. Wavefunction solution

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*^{4} to *M*^{3} the explicit matrix representation of the Hamiltonian is never generated.

### C. Renormalization

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.

## III. CARBON DIMER

C_{2} 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_{2} at stretched bond lengths, it has attracted theoretical interest as a benchmark system for methods such as coupled-cluster,^{33–36} full configuration interaction,^{37} multireference perturbation theory,^{38,39} multireference configuration interaction,^{40} 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),^{45} reduced-density-matrix theory,^{46} and valence bond theory.^{47} 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 \Sigma g + $ of eight lowest excited states: $ a 3 \Pi u , b 3 \Sigma g \u2212 , A 1 \Pi u , c 3 \Sigma u + , B 1 \Delta g , B \u2032 1 \Sigma g + , d 3 \Pi g $, and

*C*

^{1}Π

_{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 \Sigma g + , 1 \Sigma u + , 1 \Sigma g \u2212 $, and $ 1 \Sigma u \u2212 $. Note, the states of irreducible representation $ 1 \Sigma g \u2212 $ and $ 1 \Sigma u \u2212 $ 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_{2} molecule with the cc-pVQZ^{27} basis set using the Molpro ^{17} 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.^{48} 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*^{3}), although in practice for medium sized problems often *O*(*M*^{2}) 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^{42} and the CEEIS^{45} calculations and well within their reported error bars of 0.3 mE_{h}.

r . | $ 1 1 \Sigma g + $ . | $ 1 \Sigma g + $ . | $ 1 \Sigma u + $ . | $ 1 \Sigma g \u2212 $ . | $ 1 \Sigma u \u2212 $ . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1.1 | −0.761 25 | −0.761 24 | −0.621 83 | −0.502 28 | −0.565 09 | −0.294 84 | −0.222 78 | −0.286 27 | −0.243 95 | −0.214 48 | −0.413 25 | −0.309 27 | −0.192 75 |

1.2 | −0.799 24 | −0.799 20 | −0.694 59 | −0.544 90 | −0.600 50 | −0.355 81 | −0.294 06 | −0.359 14 | −0.303 43 | −0.293 63 | −0.500 18 | −0.445 33 | −0.278 46 |

1.242 53 | −0.802 69 | −0.802 64 | −0.712 08 | −0.549 53 | −0.603 38 | −0.368 59 | −0.327 73 | −0.376 14 | −0.326 26 | −0.302 24 | −0.522 84 | −0.485 20 | −0.301 05 |

1.3 | −0.799 37 | −0.799 33 | −0.726 33 | −0.548 71 | −0.599 77 | −0.378 19 | −0.363 39 | −0.390 13 | −0.347 67 | −0.305 15 | −0.545 39 | −0.524 76 | −0.320 09 |

1.4 | −0.779 70 | −0.779 65 | −0.732 67 | −0.537 76 | −0.580 49 | −0.410 85 | −0.374 09 | −0.400 43 | −0.363 28 | −0.303 54 | −0.581 36 | −0.552 43 | −0.335 51 |

1.6 | −0.724 05 | −0.724 01 | −0.704 87 | −0.510 54 | −0.523 67 | −0.455 14 | −0.346 27 | −0.404 08 | −0.338 10 | −0.292 00 | −0.618 77 | −0.547 43 | −0.351 38 |

2 | −0.645 60 | −0.645 52 | −0.614 69 | −0.492 90 | −0.473 73 | −0.423 77 | −0.322 65 | −0.382 51 | −0.315 99 | −0.276 75 | −0.613 74 | −0.506 34 | −0.431 99 |

r . | $ 1 1 \Sigma g + $ . | $ 1 \Sigma g + $ . | $ 1 \Sigma u + $ . | $ 1 \Sigma g \u2212 $ . | $ 1 \Sigma u \u2212 $ . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1.1 | −0.761 25 | −0.761 24 | −0.621 83 | −0.502 28 | −0.565 09 | −0.294 84 | −0.222 78 | −0.286 27 | −0.243 95 | −0.214 48 | −0.413 25 | −0.309 27 | −0.192 75 |

1.2 | −0.799 24 | −0.799 20 | −0.694 59 | −0.544 90 | −0.600 50 | −0.355 81 | −0.294 06 | −0.359 14 | −0.303 43 | −0.293 63 | −0.500 18 | −0.445 33 | −0.278 46 |

1.242 53 | −0.802 69 | −0.802 64 | −0.712 08 | −0.549 53 | −0.603 38 | −0.368 59 | −0.327 73 | −0.376 14 | −0.326 26 | −0.302 24 | −0.522 84 | −0.485 20 | −0.301 05 |

1.3 | −0.799 37 | −0.799 33 | −0.726 33 | −0.548 71 | −0.599 77 | −0.378 19 | −0.363 39 | −0.390 13 | −0.347 67 | −0.305 15 | −0.545 39 | −0.524 76 | −0.320 09 |

1.4 | −0.779 70 | −0.779 65 | −0.732 67 | −0.537 76 | −0.580 49 | −0.410 85 | −0.374 09 | −0.400 43 | −0.363 28 | −0.303 54 | −0.581 36 | −0.552 43 | −0.335 51 |

1.6 | −0.724 05 | −0.724 01 | −0.704 87 | −0.510 54 | −0.523 67 | −0.455 14 | −0.346 27 | −0.404 08 | −0.338 10 | −0.292 00 | −0.618 77 | −0.547 43 | −0.351 38 |

2 | −0.645 60 | −0.645 52 | −0.614 69 | −0.492 90 | −0.473 73 | −0.423 77 | −0.322 65 | −0.382 51 | −0.315 99 | −0.276 75 | −0.613 74 | −0.506 34 | −0.431 99 |

. | D_{∞h}
. | D_{2h}
. | ||
---|---|---|---|---|

M . | Discarded weight . | Energy . | Discarded weight . | Energy . |

1000 | 2.83 × 10^{−6} | −0.802 46 | 3.17 × 10^{−6} | −0.801 97 |

2000 | 4.89 × 10^{−7} | −0.802 64 | 9.02 × 10^{−7} | −0.802 49 |

3000 | 1.41 × 10^{−7} | −0.802 68 | 3.33 × 10^{−7} | −0.802 60 |

4000 | 5.50 × 10^{−8} | −0.802 69 | 1.47 × 10^{−7} | −0.802 62 |

∞ | −0.802 69 | −0.802 67 |

. | D_{∞h}
. | D_{2h}
. | ||
---|---|---|---|---|

M . | Discarded weight . | Energy . | Discarded weight . | Energy . |

1000 | 2.83 × 10^{−6} | −0.802 46 | 3.17 × 10^{−6} | −0.801 97 |

2000 | 4.89 × 10^{−7} | −0.802 64 | 9.02 × 10^{−7} | −0.802 49 |

3000 | 1.41 × 10^{−7} | −0.802 68 | 3.33 × 10^{−7} | −0.802 60 |

4000 | 5.50 × 10^{−8} | −0.802 69 | 1.47 × 10^{−7} | −0.802 62 |

∞ | −0.802 69 | −0.802 67 |

Next, electronic energies of eight low lying states are calculating at their respective equilibrium bond lengths,^{31} to obtain their excitation energies with respect to the ground state $ X 1 \Sigma 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 \Sigma u + $ and *C*^{1}Π_{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.*^{49} 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.*^{50} and most likely higher errors in the calculated excited state energies. As we have recently demonstrated elsewhere,^{51} 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,^{52} thus enabling solutions of the molecular Schödinger equation to spectroscopic accuracy.

. | r
. _{e} | Excitation energy T (eV)
. _{e} | |
---|---|---|---|

State . | Å . | Calculated . | Experimental . |

$ X 1 \Sigma g + $ | 1.242 53 | 0 | 0 |

a^{3}Π_{u} | 1.312 | 0.07 | 0.09 |

$ b 3 \Sigma g \u2212 $ | 1.369 | 0.78 | 0.80 |

A^{1}Π_{u} | 1.318 | 1.03 | 1.04 |

$ c 3 \Sigma u + $ | 1.208 | 1.17 | 1.13 |

B^{1}Δ_{g} | 1.385 | 1.49 | 1.50 |

$ B \u2032 1 \Sigma g + $ | 1.377 | 1.90 | 1.91 |

d^{3}Π_{g} | 1.266 | 2.51 | 2.48 |

C^{1}Π_{g} | 1.255 | 4.31 | 4.25 |

. | r
. _{e} | Excitation energy T (eV)
. _{e} | |
---|---|---|---|

State . | Å . | Calculated . | Experimental . |

$ X 1 \Sigma g + $ | 1.242 53 | 0 | 0 |

a^{3}Π_{u} | 1.312 | 0.07 | 0.09 |

$ b 3 \Sigma g \u2212 $ | 1.369 | 0.78 | 0.80 |

A^{1}Π_{u} | 1.318 | 1.03 | 1.04 |

$ c 3 \Sigma u + $ | 1.208 | 1.17 | 1.13 |

B^{1}Δ_{g} | 1.385 | 1.49 | 1.50 |

$ B \u2032 1 \Sigma g + $ | 1.377 | 1.90 | 1.91 |

d^{3}Π_{g} | 1.266 | 2.51 | 2.48 |

C^{1}Π_{g} | 1.255 | 4.31 | 4.25 |

Next, four separate DMRG state-averaged calculations were performed on three states each of irreducible representations $ 1 \Sigma g + , 1 \Sigma u + , 1 \Sigma g \u2212 $, and $ 1 \Sigma u \u2212 $ 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.

## IV. CONCLUSION

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_{2} 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.

## Acknowledgments

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.