A general non-Abelian density matrix renormalization group algorithm with application to the C2 dimer Sandeep Sharma Citation: The Journal of Chemical Physics 142, 024107 (2015); doi: 10.1063/1.4905237 View online: http://dx.doi.org/10.1063/1.4905237 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/2?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Linear response theory for the density matrix renormalization group: Efficient algorithms for strongly correlated excited states J. Chem. Phys. 140, 024108 (2014); 10.1063/1.4860375 Second-order perturbation theory with a density matrix renormalization group self-consistent field reference function: Theory and application to the study of chromium dimer J. Chem. Phys. 135, 094104 (2011); 10.1063/1.3629454 The density matrix renormalization group self-consistent field method: Orbital optimization with the density matrix renormalization group method in the active space J. Chem. Phys. 128, 144116 (2008); 10.1063/1.2883981 Quantum chemistry using the density matrix renormalization group II J. Chem. Phys. 119, 4148 (2003); 10.1063/1.1593627 The density matrix renormalization group method: Application to the PPP model of a cyclic polyene chain J. Chem. Phys. 108, 9246 (1998); 10.1063/1.476379
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
THE JOURNAL OF CHEMICAL PHYSICS 142, 024107 (2015)
A general non-Abelian density matrix renormalization group algorithm with application to the C2 dimer Sandeep Sharmaa) Department of Chemistry, Frick Laboratory, Princeton University, Princeton, New Jersey 08544, USA
(Received 25 August 2014; accepted 18 December 2014; published online 9 January 2015) 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 C2 dimer in the cc-pVQZ basis set (with a frozen-core), corresponding to a Hilbert space dimension of 1012 many-body states. While our calculated energy lies within the 0.3 mEh 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 mEh , much more accurate than these previous estimates. Due to the additional efficiency afforded by the algorithm, the excitation energies (Te ) of eight lowest lying excited states: a3Πu , b3Σg−, A1Πu , c3Σu+, B 1∆g , B ′1Σ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Σg+, 1Σu+, 1Σg−, and 1Σu−, to an estimated accuracy of 0.1 mEh of the exact result in this basis. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905237] I. INTRODUCTION
The density matrix renormalization group (DMRG) was first introduced over 20 years ago1,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 Martin3 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 theory16 and has also been incorporated into several major quantum chemistry packages such as M,17 Q-,18 and O.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 work21,22 that introduced the WignerEckart theorem to gain large efficiencies and the “quasi-density matrix” to ensure that only spin multiplet states are retained during renormalization. Zgid and Nooijen23 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
a)Electronic address: [email protected]
0021-9606/2015/142(2)/024107/6/$30.00
Weichselbaum26 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 C2 molecule at various bond lengths in a cc-pVQZ basis set,27 to 10 µEh and 0.1 mEh accuracies, respectively. Further, excitation energies (Te ) of eight lowest excited state of the C2 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 WignerEckart theorem, shown in Eq. (1), where ψ γ i and φ β k are states belonging to the i and k rows of irreducible representations γ
142, 024107-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-2
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
TABLE I. Normal and complementary operators used in the usual DMRG algorithm, where the indices are spin indices. In the definition of complementary operators, the repeated indices are summed over. Normal operator
Complementary operator Hˆ =
Iˆ
t i j a †i a j + ⟨i j |k l⟩a †i a †j a l a k Rˆ i = t i j a j + ⟨i j |k l⟩a †j a l a k Pˆ i j = ⟨i j |k l⟩a l a k Qˆ i j = (⟨ik | jl⟩ − ⟨ik |l j⟩)a †k a l
a †i Aˆ i j = a †i a †j Bˆ i j = a †i a j
and β, respectively, and Oˆ α j is the tensor operator of irreducible representation α that transforms as its jth row. ⟨ψ γ i |Oˆ α j |φ β k ⟩ = ⟨ βk α j |γi ⟩⟨ψ γ ||Oˆ α ||φ β ⟩.
(1)
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.
γ 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 ClebschGordan coefficients as shown in Eq. (2). γ [α β ] γ α † β † CCI Jk = CCI Jk i j = ⟨α i β j |γk ⟩a I i a J j , (2) ij δ CCD I Jl K
A. Blocking
=
ij δ l [α i θ m [β j γ k ]] CCD I J K
i jkm
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. 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
=
α † β † γ
⟨ β j γk |θ m ⟩⟨α i θ m |δl ⟩a I i a J j a Kk .
(3)
i jkm γ [α β ]
α † β †
Here, the intermediate CCI Jk 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 (CCD I Jl K ) 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 (CCI Jk i j ) from Eq. (2) is used on the RHS of
TABLE II. Definitions of the normal and complementary operators used in the non-Abelian symmetry adapted DMRG algorithm. C and D denote the creation and destruction operators, for example, “C DD I J K ” represents the product of three tensor operators, one creation on spatial orbital I and two destructions on spatial orbitals J and K. The superscript A g is used to denote the fully symmetric irreducible representation of the symmetry group of interest. The symmetry label with a negative sign (e.g., −α i ) before it represents the row of the irreducible representation (α) that can couple to the original irreducible representation (α i ) to give a fully symmetric irreducible representation ( A g ), i.e., ⟨α i − α i | Ag ⟩ , 0. Normal operator
Complementary operator αβδ ξ α† β† ξ δ α β α † β Hˆ = α i β j t I Ji j a I i a J j + ν I Ji Kj Lk l a I i a J j a Ll a Kk
Iˆ
aI i
α †
α α −α −α Rˆ I i = I J α I t I Ji i a J i + J K L β j δ k ξ l θ m
γ m [α i β j ] γ Aˆ I m J = CC I J
γ Pˆ I Jm = K Lδ k ξ l α i β j
γ m [α β ] γ Bˆ I Jm = C D I J i j
γ Qˆ I m K Lδ k ξ l α i β j J =
αi β j δ ξ ν I J K L k l ⟨θ m |ξ l δ k ⟩⟨−α i |θ m β j ⟩ −α [β θ m [ξ l δ k ]] C DD J LiK j ⟨α i −α i | Ag ⟩
αi β j δ ξ ν I J K L k l ⟨−γ m |ξ l δ k ⟩⟨γ m |α i β j ⟩ −γ [ξ δ ] DD L Km l k ⟨γ m −γ m | Ag ⟩ αi δ β j ξl
(ν I K JkL
αiδ ξ β j
−ν I K LkJ l )⟨−γ m |ξ l δ k ⟩⟨γ m |α i β j ⟩ −γ [δ ξ ] C D K Lm k l ⟨γ m −γ m | Ag ⟩
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-3
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
the definition. Notation Definition γ[α β] α γ ˆβ ˆ ˆ N = L ⊗ M { Nˆ γ k = Nˆ γ k [α i β j ]}.
(4)
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.
ij
H[S] = H[L] ⊗ A g 1[•] + H[•] ⊗ A g 1[L] + (1/2)
(
Ag γ aγ† RJ [•] + aγJ [L] ⊗ A g Rγ† J [L] ⊗ J [•]
)
I ∈L
+(1/2)
(
Ag aγ† J [•] ⊗
RγJ [L] + aγJ [•] ⊗ A g
( ) ) A g γ† Rγ† [L] + (1/2) AγI J [•] ⊗ A g PγI J [L] + Aγ† P I J [L] J I J [•] ⊗
I ∈•
+(1/2)
I J ∈•
BγI J [•] ⊗ A g QγI J [L] ,
I J ∈•
RαI [S]
=
RαI [L] ⊗ α 1[•] + RαI [•] ⊗ α 1[L] + +
W (α β Ag γ; γα)
(
( ) W (α β Ag γ; γα) 2PγI J [L] ⊗ α aJβ†[•] + QγI J [L] ⊗ α aJβ [•]
J ∈• ) γ α β 2P I J [•] ⊗ aJ [L] + QγI J [•] ⊗ α aJβ [L] ,
J ∈L
PγI J [S]
=
PγI J [L] ⊗γ 1[•] + PγI J [•] ⊗γ 1[L] +
ν α i β j δ k ξ l ⟨−γm |ξl δ k ⟩⟨γm |α i β j ⟩ *. +/ aξ [L] ⊗γ aδ [•], IJKL K L ⟨γ − γ | Ag⟩ m m δ ξ α β K ∈•, k l i j L∈L
QγI J [S] = QγI J [L] ⊗γ 1[•] + QγI J [•] ⊗γ 1[L] +
(ν α i δ k β j ξ l − ν α i δ k ξ l β j )⟨−γm |ξl δ k ⟩⟨γm |α i β j ⟩ *. +/ aδ†[•] ⊗γ aξ [L]. IKJL I K LJ K L ⟨γ − γ | Ag⟩ m m δ ξ α β K ∈•, k l i j L∈L
(5)
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 (α β Ag γ;γα) 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 nonAbelian group of interest, very substantial computational savings are possible.
α
α
β
β
⟨φ1 1 ⊗ α φ2 2||O1γ ⊗ ξ O2δ ∥ψ1 1 ⊗ β ψ2 2⟩ β1 β2 β δ ξ ⟨φ1α1 ∥O1γ ∥ψ1β1⟩⟨φ2α2 ∥O2δ ∥ψ2β2⟩. γ = ψ 1ψ 2φ 1φ 2 α1 α2 α (6) 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 2n α 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 operators29 or by just acting the elements of a tensor operator on the vacuum state. The resulting number of multiplets obtained is less than 2n α . 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
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-4
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
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 γ ⟩ = Cγ α β ∥φ1α ⟩⊗γ ∥φ2β ⟩,
(7)
φ1 φ2
β β α 2 1 δ Ag γ Vα α β = φ1 φ2 ψ 1ψ 2φ 1φ 2 α1 α2 α β α β α ×⟨φ1 1 ∥O1γ ∥ψ1 1⟩⟨φ2 2 ∥O2δ ∥ψ2 2⟩Cαα β .
ψ1 ψ2
(8)
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 “quasidensity 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. A Γφ αgψ α = Cγ α β Cγ α β . (9) 1
1
β
φ2
φ1 φ2
ψ1 φ2
III. CARBON DIMER
C2 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 C2 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 the-
ory,46 and valence bond theory.47 Here, we use our non-Abelian symmetry adapted DMRG program to calculate the excitation energies (Te ) with respect to the ground state X 1Σg+ of eight lowest excited states: a3Πu , b3Σg−, A1Πu , c3Σu+, B1∆g , B ′1Σg+, d 3Π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Σ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 C2 molecule with the cc-pVQZ27 basis set using the M 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 µEh as shown in Table IV. The table also shows the discarded weight and calculated energies of a DMRG calculation when only the SU(2) × D2h group is used. The data demonstrate that roughly twice the number of retained renormalized states (M) are needed when only the D2h 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-FCIQMC42 and the CEEIS45 calculations and well within their reported error bars of 0.3 mEh . 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Σg+. The energies are tabulated in Table V and show
TABLE III. DMRG energy (E + 75.0 in Eh ) of various states of the C2 molecule calculated using the cc-pVQZ basis set with a frozen-core. The second column shows the ground state energies calculated to high accuracy (error ≈ 10 µEh ) with the number of renormalized (multiplet) states M of up to 6000. The last 12 columns show the result of four separate DMRG state-averaged calculations performed for three states each of the irreducible representations 1Σg+, 1Σu+, 1Σg−, and 1Σu− with a maximum M of 4000. The difference between the second and the third column gives us an estimate of the error in our state-averaged calculations, which is less than 0.1 mEh along the entirety of the curve. r 1.1 1.2 1.242 53 1.3 1.4 1.6 2
11Σg+ −0.761 25 −0.799 24 −0.802 69 −0.799 37 −0.779 70 −0.724 05 −0.645 60
1Σ + g
−0.761 24 −0.799 20 −0.802 64 −0.799 33 −0.779 65 −0.724 01 −0.645 52
−0.621 83 −0.694 59 −0.712 08 −0.726 33 −0.732 67 −0.704 87 −0.614 69
1Σ + u
−0.502 28 −0.544 90 −0.549 53 −0.548 71 −0.537 76 −0.510 54 −0.492 90
−0.565 09 −0.600 50 −0.603 38 −0.599 77 −0.580 49 −0.523 67 −0.473 73
−0.294 84 −0.355 81 −0.368 59 −0.378 19 −0.410 85 −0.455 14 −0.423 77
1Σ − g
−0.222 78 −0.294 06 −0.327 73 −0.363 39 −0.374 09 −0.346 27 −0.322 65
−0.286 27 −0.359 14 −0.376 14 −0.390 13 −0.400 43 −0.404 08 −0.382 51
−0.243 95 −0.303 43 −0.326 26 −0.347 67 −0.363 28 −0.338 10 −0.315 99
1Σ − u
−0.214 48 −0.293 63 −0.302 24 −0.305 15 −0.303 54 −0.292 00 −0.276 75
−0.413 25 −0.500 18 −0.522 84 −0.545 39 −0.581 36 −0.618 77 −0.613 74
−0.309 27 −0.445 33 −0.485 20 −0.524 76 −0.552 43 −0.547 43 −0.506 34
−0.192 75 −0.278 46 −0.301 05 −0.320 09 −0.335 51 −0.351 38 −0.431 99
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-5
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
TABLE IV. The table gives the energies (E+75) in Eh and discarded weights of DMRG calculations on the X 1Σg+ state of the C2 molecule at the equilibrium bond length of 1.242 53 Å. Columns 2 and 3 give data for calculation where the full SU(2) × D∞h group is used and columns 4 and 5 give data for calculation where only the SU(2) × D2h group is used. D∞h M 1000 2000 3000 4000 ∞
D 2h
Discarded weight
Energy
Discarded weight
Energy
2.83 × 10−6 4.89 × 10−7 1.41 × 10−7 5.50 × 10−8
−0.802 46 −0.802 64 −0.802 68 −0.802 69 −0.802 69
3.17 × 10−6 9.02 × 10−7 3.33 × 10−7 1.47 × 10−7
−0.801 97 −0.802 49 −0.802 60 −0.802 62 −0.802 67
good agreement with the experimentally measured excitation energies, with all errors being less than 0.03 eV except c3Σ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 ccpVQZ 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 mEh 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. 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 mEh for all the geometries.
FIG. 1. The energies of various states of the C2 molecule relative to complete dissociation, consisting of two carbon atoms in triplet states. The shaded circles show the regions of avoided crossing between states belonging to the same irreducible representation. As expected, we see that the energy of three singlet states, two 1Σg+ and one 1Σu− , approach the value of two carbon triplet atoms at stretched bonds.
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 nonAbelian point group symmetries. We used the resulting implementation to calculate the ground state potential energy curve of the C2 molecule with a cc-pVQZ basis set (and frozen core) to an unprecedented (near-exact) accuracy of 10 µEh . 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
TABLE V. The table shows the calculated and experimental excitation energies, in electron volts (eV), of eight lowest lying excited states of the C2 dimer. The electronic energies of the excited states are calculated at the bond length, in Å, given in the second column of the table. All the bond lengths except for the X 1Σg+ state, which is taken as the equilibrium bond length calculated by Booth et al.,41 are derived from experimental observations.31. Excitation energy Te (eV) State
re Å
Calculated
Experimental
X 1Σg+ a 3Π u b 3Σg− A1Π u c 3Σu+ B 1∆ g B′1Σg+ d 3Π g C 1Π g
1.242 53 1.312 1.369 1.318 1.208 1.385 1.377 1.266 1.255
0 0.07 0.78 1.03 1.17 1.49 1.90 2.51 4.31
0 0.09 0.80 1.04 1.13 1.50 1.91 2.48 4.25
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. 1S.
R. White, Phys. Rev. Lett. 69, 2863 (1992). R. White, Phys. Rev. B 48, 10345 (1993). 3S. R. White and R. L. Martin, J. Chem. Phys. 110, 4127 (1999). 4G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465 (2011). 5Ö. Legeza, R. Noack, J. Sólyom, and L. Tincani, Computational ManyParticle Physics, Lecture Notes in Physics, edited by H. Fehske, R. Schneider, and A. Weiße (Springer, Berlin, Heidelberg, 2008), Vol. 739, pp. 653–664. 6K. H. Marti and M. Reiher, Phys. Chem. Chem. Phys. 13, 6750 (2011). 7Y. Kurashige and T. Yanai, J. Phys. Chem. 130, 234114 (2009). 8S. Wouters and D. Van Neck, Eur. Phys. J. D 68, 272 (2014). 9Y. Kurashige, G. K.-L. Chan, and T. Yanai, Nat. Chem. 5, 660 (2013). 10S. Sharma, K. Sivalingam, F. Neese, and G. K.-L. Chan, Nat. Chem. 6, 927–933 (2014). 11D. Zgid and M. Nooijen, J. Chem. Phys. 128, 144116 (2008). 2S.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-6 12D.
Sandeep Sharma
Ghosh, J. Hachmann, T. Yanai, and G. K. L. Chan, J. Chem. Phys. 128, 144117 (2008). 13Y. Kurashige and T. Yanai, J. Chem. Phys. 135, 094104 (2011). 14F. Liu, Y. Kurashige, T. Yanai, and K. Morokuma, J. Chem. Theory Comput. 9, 4462 (2013). 15M. Saitow, Y. Kurashige, and T. Yanai, J. Chem. Phys. 139, 044118 (2013). 16T. Yanai, Y. Kurashige, E. Neuscamman, and G. K.-L. Chan, J. Chem. Phys. 132, 024105 (2010). 17H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 242 (2012). 18Y. Shao et al., Mol. Phys. (published online 2014); available at http://www. tandfonline.com/doi/abs/10.1080/00268976.2014.952696#.VKVQvWTF_ vQ. 19F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 73 (2012). 20G. Sierra and T. Nishino, Nucl. Phys. B 495, 505 (1997). 21I. P. McCulloch and M. Gulacsi, Aust. J. Phys. 53, 597 (2000). 22I. P. McCulloch and M. Gulacsi, Europhys. Lett. 57, 852 (2002). 23D. Zgid and M. Nooijen, J. Chem. Phys. 128, 14107 (2008). 24S. Wouters, W. Poelmans, P. W. Ayers, and D. Van Neck, Comput. Phys. Commun. 185, 1501 (2014). 25S. Sharma and G. K.-L. Chan, J. Chem. Phys. 136, 124121 (2012). 26A. Weichselbaum, Ann. Phys. 327, 2972 (2012). 27T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). 28D. M. Brink and G. R. Satchler, Angular Momentum, Dover Books on Physics (Dover Publications, 1968). 29H. Morton, Group Theory and Its Application to Physical Problems, Oxford Library of the Physical Sciences (Clarendon Press, 1989). 30W. Weltner and R. J. Van Zee, Chem. Rev. 89, 1713 (1989). 31M. Martin, J. Photochem. Photobiol., A 66, 263 (1992).
J. Chem. Phys. 142, 024107 (2015) 32A.
Van Orden and R. J. Saykally, Chem. Rev. 98, 2313 (1998). D. Watts and R. J. Bartlett, J. Chem. Phys. 96, 6073 (1992). 34C. D. Sherrill and P. Piecuch, J. Chem. Phys. 122, 124104 (2005). 35D. I. Lyakh, V. V. Ivanov, and L. Adamowicz, J. Chem. Phys. 128, 074101 (2008). 36D. Datta, L. Kong, and M. Nooijen, J. Chem. Phys. 134, 214116 (2011). 37M. L. Abrams and C. D. Sherrill, J. Chem. Phys. 121, 9211 (2004). 38U. S. Mahapatra, S. Chattopadhyay, and R. K. Chaudhuri, J. Chem. Phys. 129, 024108 (2008). 39W. Jiang and A. K. Wilson, J. Chem. Phys. 134, 034101 (2011). 40A. D. Pradhan, H. Partridge, and C. W. Bauschlicher, J. Chem. Phys. 101, 3857 (1994). 41G. H. Booth, D. Cleland, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 135, 84104 (2011). 42D. Cleland, G. H. Booth, C. Overy, and A. Alavi, J. Chem. Theory Comput. 8, 4138 (2012). 43S. Ten-no, J. Chem. Phys. 138, 164126 (2013). 44Y. Ohtsuka and S. Ten-no, AIP Conference Proceedings 1456, 97 (2012). 45L. Bytautas and K. Ruedenberg, J. Chem. Phys. 122, 154110 (2005). 46G. Gidofalvi and D. A. Mazziotti, J. Chem. Phys. 122, 194104 (2005). 47S. Shaik et al., Nat. Chem. 4, 195 (2012). 48S. L. Altmann and P. Herzig, Point-Group Theory Tables (Clarendon Press, Oxford, 1994). 49G. H. Booth, D. Cleland, A. Alavi, and D. P. Tew, J. Chem. Phys. 137, 164112 (2012). 50S. Wouters, W. Poelmans, P. W. Ayers, and D. V. Neck, Comput. Phys. Commun. 185, 1501 (2014). 51S. Sharma, T. Yanai, G. H. Booth, C. J. Umrigar, and G. K.-L. Chan, J. Chem. Phys. 140, 104112 (2014). 52T. Yanai and T. Shiozaki, J. Chem. Phys. 136, 84107 (2012). 33J.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
THE JOURNAL OF CHEMICAL PHYSICS 142, 024107 (2015)
A general non-Abelian density matrix renormalization group algorithm with application to the C2 dimer Sandeep Sharmaa) Department of Chemistry, Frick Laboratory, Princeton University, Princeton, New Jersey 08544, USA
(Received 25 August 2014; accepted 18 December 2014; published online 9 January 2015) 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 C2 dimer in the cc-pVQZ basis set (with a frozen-core), corresponding to a Hilbert space dimension of 1012 many-body states. While our calculated energy lies within the 0.3 mEh 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 mEh , much more accurate than these previous estimates. Due to the additional efficiency afforded by the algorithm, the excitation energies (Te ) of eight lowest lying excited states: a3Πu , b3Σg−, A1Πu , c3Σu+, B 1∆g , B ′1Σ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Σg+, 1Σu+, 1Σg−, and 1Σu−, to an estimated accuracy of 0.1 mEh of the exact result in this basis. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905237] I. INTRODUCTION
The density matrix renormalization group (DMRG) was first introduced over 20 years ago1,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 Martin3 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 theory16 and has also been incorporated into several major quantum chemistry packages such as M,17 Q-,18 and O.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 work21,22 that introduced the WignerEckart theorem to gain large efficiencies and the “quasi-density matrix” to ensure that only spin multiplet states are retained during renormalization. Zgid and Nooijen23 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
a)Electronic address: [email protected]
0021-9606/2015/142(2)/024107/6/$30.00
Weichselbaum26 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 C2 molecule at various bond lengths in a cc-pVQZ basis set,27 to 10 µEh and 0.1 mEh accuracies, respectively. Further, excitation energies (Te ) of eight lowest excited state of the C2 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 WignerEckart theorem, shown in Eq. (1), where ψ γ i and φ β k are states belonging to the i and k rows of irreducible representations γ
142, 024107-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-2
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
TABLE I. Normal and complementary operators used in the usual DMRG algorithm, where the indices are spin indices. In the definition of complementary operators, the repeated indices are summed over. Normal operator
Complementary operator Hˆ =
Iˆ
t i j a †i a j + ⟨i j |k l⟩a †i a †j a l a k Rˆ i = t i j a j + ⟨i j |k l⟩a †j a l a k Pˆ i j = ⟨i j |k l⟩a l a k Qˆ i j = (⟨ik | jl⟩ − ⟨ik |l j⟩)a †k a l
a †i Aˆ i j = a †i a †j Bˆ i j = a †i a j
and β, respectively, and Oˆ α j is the tensor operator of irreducible representation α that transforms as its jth row. ⟨ψ γ i |Oˆ α j |φ β k ⟩ = ⟨ βk α j |γi ⟩⟨ψ γ ||Oˆ α ||φ β ⟩.
(1)
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.
γ 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 ClebschGordan coefficients as shown in Eq. (2). γ [α β ] γ α † β † CCI Jk = CCI Jk i j = ⟨α i β j |γk ⟩a I i a J j , (2) ij δ CCD I Jl K
A. Blocking
=
ij δ l [α i θ m [β j γ k ]] CCD I J K
i jkm
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. 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
=
α † β † γ
⟨ β j γk |θ m ⟩⟨α i θ m |δl ⟩a I i a J j a Kk .
(3)
i jkm γ [α β ]
α † β †
Here, the intermediate CCI Jk 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 (CCD I Jl K ) 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 (CCI Jk i j ) from Eq. (2) is used on the RHS of
TABLE II. Definitions of the normal and complementary operators used in the non-Abelian symmetry adapted DMRG algorithm. C and D denote the creation and destruction operators, for example, “C DD I J K ” represents the product of three tensor operators, one creation on spatial orbital I and two destructions on spatial orbitals J and K. The superscript A g is used to denote the fully symmetric irreducible representation of the symmetry group of interest. The symmetry label with a negative sign (e.g., −α i ) before it represents the row of the irreducible representation (α) that can couple to the original irreducible representation (α i ) to give a fully symmetric irreducible representation ( A g ), i.e., ⟨α i − α i | Ag ⟩ , 0. Normal operator
Complementary operator αβδ ξ α† β† ξ δ α β α † β Hˆ = α i β j t I Ji j a I i a J j + ν I Ji Kj Lk l a I i a J j a Ll a Kk
Iˆ
aI i
α †
α α −α −α Rˆ I i = I J α I t I Ji i a J i + J K L β j δ k ξ l θ m
γ m [α i β j ] γ Aˆ I m J = CC I J
γ Pˆ I Jm = K Lδ k ξ l α i β j
γ m [α β ] γ Bˆ I Jm = C D I J i j
γ Qˆ I m K Lδ k ξ l α i β j J =
αi β j δ ξ ν I J K L k l ⟨θ m |ξ l δ k ⟩⟨−α i |θ m β j ⟩ −α [β θ m [ξ l δ k ]] C DD J LiK j ⟨α i −α i | Ag ⟩
αi β j δ ξ ν I J K L k l ⟨−γ m |ξ l δ k ⟩⟨γ m |α i β j ⟩ −γ [ξ δ ] DD L Km l k ⟨γ m −γ m | Ag ⟩ αi δ β j ξl
(ν I K JkL
αiδ ξ β j
−ν I K LkJ l )⟨−γ m |ξ l δ k ⟩⟨γ m |α i β j ⟩ −γ [δ ξ ] C D K Lm k l ⟨γ m −γ m | Ag ⟩
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-3
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
the definition. Notation Definition γ[α β] α γ ˆβ ˆ ˆ N = L ⊗ M { Nˆ γ k = Nˆ γ k [α i β j ]}.
(4)
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.
ij
H[S] = H[L] ⊗ A g 1[•] + H[•] ⊗ A g 1[L] + (1/2)
(
Ag γ aγ† RJ [•] + aγJ [L] ⊗ A g Rγ† J [L] ⊗ J [•]
)
I ∈L
+(1/2)
(
Ag aγ† J [•] ⊗
RγJ [L] + aγJ [•] ⊗ A g
( ) ) A g γ† Rγ† [L] + (1/2) AγI J [•] ⊗ A g PγI J [L] + Aγ† P I J [L] J I J [•] ⊗
I ∈•
+(1/2)
I J ∈•
BγI J [•] ⊗ A g QγI J [L] ,
I J ∈•
RαI [S]
=
RαI [L] ⊗ α 1[•] + RαI [•] ⊗ α 1[L] + +
W (α β Ag γ; γα)
(
( ) W (α β Ag γ; γα) 2PγI J [L] ⊗ α aJβ†[•] + QγI J [L] ⊗ α aJβ [•]
J ∈• ) γ α β 2P I J [•] ⊗ aJ [L] + QγI J [•] ⊗ α aJβ [L] ,
J ∈L
PγI J [S]
=
PγI J [L] ⊗γ 1[•] + PγI J [•] ⊗γ 1[L] +
ν α i β j δ k ξ l ⟨−γm |ξl δ k ⟩⟨γm |α i β j ⟩ *. +/ aξ [L] ⊗γ aδ [•], IJKL K L ⟨γ − γ | Ag⟩ m m δ ξ α β K ∈•, k l i j L∈L
QγI J [S] = QγI J [L] ⊗γ 1[•] + QγI J [•] ⊗γ 1[L] +
(ν α i δ k β j ξ l − ν α i δ k ξ l β j )⟨−γm |ξl δ k ⟩⟨γm |α i β j ⟩ *. +/ aδ†[•] ⊗γ aξ [L]. IKJL I K LJ K L ⟨γ − γ | Ag⟩ m m δ ξ α β K ∈•, k l i j L∈L
(5)
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 (α β Ag γ;γα) 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 nonAbelian group of interest, very substantial computational savings are possible.
α
α
β
β
⟨φ1 1 ⊗ α φ2 2||O1γ ⊗ ξ O2δ ∥ψ1 1 ⊗ β ψ2 2⟩ β1 β2 β δ ξ ⟨φ1α1 ∥O1γ ∥ψ1β1⟩⟨φ2α2 ∥O2δ ∥ψ2β2⟩. γ = ψ 1ψ 2φ 1φ 2 α1 α2 α (6) 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 2n α 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 operators29 or by just acting the elements of a tensor operator on the vacuum state. The resulting number of multiplets obtained is less than 2n α . 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
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-4
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
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 γ ⟩ = Cγ α β ∥φ1α ⟩⊗γ ∥φ2β ⟩,
(7)
φ1 φ2
β β α 2 1 δ Ag γ Vα α β = φ1 φ2 ψ 1ψ 2φ 1φ 2 α1 α2 α β α β α ×⟨φ1 1 ∥O1γ ∥ψ1 1⟩⟨φ2 2 ∥O2δ ∥ψ2 2⟩Cαα β .
ψ1 ψ2
(8)
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 “quasidensity 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. A Γφ αgψ α = Cγ α β Cγ α β . (9) 1
1
β
φ2
φ1 φ2
ψ1 φ2
III. CARBON DIMER
C2 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 C2 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 the-
ory,46 and valence bond theory.47 Here, we use our non-Abelian symmetry adapted DMRG program to calculate the excitation energies (Te ) with respect to the ground state X 1Σg+ of eight lowest excited states: a3Πu , b3Σg−, A1Πu , c3Σu+, B1∆g , B ′1Σg+, d 3Π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Σ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 C2 molecule with the cc-pVQZ27 basis set using the M 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 µEh as shown in Table IV. The table also shows the discarded weight and calculated energies of a DMRG calculation when only the SU(2) × D2h group is used. The data demonstrate that roughly twice the number of retained renormalized states (M) are needed when only the D2h 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-FCIQMC42 and the CEEIS45 calculations and well within their reported error bars of 0.3 mEh . 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Σg+. The energies are tabulated in Table V and show
TABLE III. DMRG energy (E + 75.0 in Eh ) of various states of the C2 molecule calculated using the cc-pVQZ basis set with a frozen-core. The second column shows the ground state energies calculated to high accuracy (error ≈ 10 µEh ) with the number of renormalized (multiplet) states M of up to 6000. The last 12 columns show the result of four separate DMRG state-averaged calculations performed for three states each of the irreducible representations 1Σg+, 1Σu+, 1Σg−, and 1Σu− with a maximum M of 4000. The difference between the second and the third column gives us an estimate of the error in our state-averaged calculations, which is less than 0.1 mEh along the entirety of the curve. r 1.1 1.2 1.242 53 1.3 1.4 1.6 2
11Σg+ −0.761 25 −0.799 24 −0.802 69 −0.799 37 −0.779 70 −0.724 05 −0.645 60
1Σ + g
−0.761 24 −0.799 20 −0.802 64 −0.799 33 −0.779 65 −0.724 01 −0.645 52
−0.621 83 −0.694 59 −0.712 08 −0.726 33 −0.732 67 −0.704 87 −0.614 69
1Σ + u
−0.502 28 −0.544 90 −0.549 53 −0.548 71 −0.537 76 −0.510 54 −0.492 90
−0.565 09 −0.600 50 −0.603 38 −0.599 77 −0.580 49 −0.523 67 −0.473 73
−0.294 84 −0.355 81 −0.368 59 −0.378 19 −0.410 85 −0.455 14 −0.423 77
1Σ − g
−0.222 78 −0.294 06 −0.327 73 −0.363 39 −0.374 09 −0.346 27 −0.322 65
−0.286 27 −0.359 14 −0.376 14 −0.390 13 −0.400 43 −0.404 08 −0.382 51
−0.243 95 −0.303 43 −0.326 26 −0.347 67 −0.363 28 −0.338 10 −0.315 99
1Σ − u
−0.214 48 −0.293 63 −0.302 24 −0.305 15 −0.303 54 −0.292 00 −0.276 75
−0.413 25 −0.500 18 −0.522 84 −0.545 39 −0.581 36 −0.618 77 −0.613 74
−0.309 27 −0.445 33 −0.485 20 −0.524 76 −0.552 43 −0.547 43 −0.506 34
−0.192 75 −0.278 46 −0.301 05 −0.320 09 −0.335 51 −0.351 38 −0.431 99
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-5
Sandeep Sharma
J. Chem. Phys. 142, 024107 (2015)
TABLE IV. The table gives the energies (E+75) in Eh and discarded weights of DMRG calculations on the X 1Σg+ state of the C2 molecule at the equilibrium bond length of 1.242 53 Å. Columns 2 and 3 give data for calculation where the full SU(2) × D∞h group is used and columns 4 and 5 give data for calculation where only the SU(2) × D2h group is used. D∞h M 1000 2000 3000 4000 ∞
D 2h
Discarded weight
Energy
Discarded weight
Energy
2.83 × 10−6 4.89 × 10−7 1.41 × 10−7 5.50 × 10−8
−0.802 46 −0.802 64 −0.802 68 −0.802 69 −0.802 69
3.17 × 10−6 9.02 × 10−7 3.33 × 10−7 1.47 × 10−7
−0.801 97 −0.802 49 −0.802 60 −0.802 62 −0.802 67
good agreement with the experimentally measured excitation energies, with all errors being less than 0.03 eV except c3Σ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 ccpVQZ 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 mEh 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. 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 mEh for all the geometries.
FIG. 1. The energies of various states of the C2 molecule relative to complete dissociation, consisting of two carbon atoms in triplet states. The shaded circles show the regions of avoided crossing between states belonging to the same irreducible representation. As expected, we see that the energy of three singlet states, two 1Σg+ and one 1Σu− , approach the value of two carbon triplet atoms at stretched bonds.
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 nonAbelian point group symmetries. We used the resulting implementation to calculate the ground state potential energy curve of the C2 molecule with a cc-pVQZ basis set (and frozen core) to an unprecedented (near-exact) accuracy of 10 µEh . 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
TABLE V. The table shows the calculated and experimental excitation energies, in electron volts (eV), of eight lowest lying excited states of the C2 dimer. The electronic energies of the excited states are calculated at the bond length, in Å, given in the second column of the table. All the bond lengths except for the X 1Σg+ state, which is taken as the equilibrium bond length calculated by Booth et al.,41 are derived from experimental observations.31. Excitation energy Te (eV) State
re Å
Calculated
Experimental
X 1Σg+ a 3Π u b 3Σg− A1Π u c 3Σu+ B 1∆ g B′1Σg+ d 3Π g C 1Π g
1.242 53 1.312 1.369 1.318 1.208 1.385 1.377 1.266 1.255
0 0.07 0.78 1.03 1.17 1.49 1.90 2.51 4.31
0 0.09 0.80 1.04 1.13 1.50 1.91 2.48 4.25
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. 1S.
R. White, Phys. Rev. Lett. 69, 2863 (1992). R. White, Phys. Rev. B 48, 10345 (1993). 3S. R. White and R. L. Martin, J. Chem. Phys. 110, 4127 (1999). 4G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465 (2011). 5Ö. Legeza, R. Noack, J. Sólyom, and L. Tincani, Computational ManyParticle Physics, Lecture Notes in Physics, edited by H. Fehske, R. Schneider, and A. Weiße (Springer, Berlin, Heidelberg, 2008), Vol. 739, pp. 653–664. 6K. H. Marti and M. Reiher, Phys. Chem. Chem. Phys. 13, 6750 (2011). 7Y. Kurashige and T. Yanai, J. Phys. Chem. 130, 234114 (2009). 8S. Wouters and D. Van Neck, Eur. Phys. J. D 68, 272 (2014). 9Y. Kurashige, G. K.-L. Chan, and T. Yanai, Nat. Chem. 5, 660 (2013). 10S. Sharma, K. Sivalingam, F. Neese, and G. K.-L. Chan, Nat. Chem. 6, 927–933 (2014). 11D. Zgid and M. Nooijen, J. Chem. Phys. 128, 144116 (2008). 2S.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
024107-6 12D.
Sandeep Sharma
Ghosh, J. Hachmann, T. Yanai, and G. K. L. Chan, J. Chem. Phys. 128, 144117 (2008). 13Y. Kurashige and T. Yanai, J. Chem. Phys. 135, 094104 (2011). 14F. Liu, Y. Kurashige, T. Yanai, and K. Morokuma, J. Chem. Theory Comput. 9, 4462 (2013). 15M. Saitow, Y. Kurashige, and T. Yanai, J. Chem. Phys. 139, 044118 (2013). 16T. Yanai, Y. Kurashige, E. Neuscamman, and G. K.-L. Chan, J. Chem. Phys. 132, 024105 (2010). 17H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 242 (2012). 18Y. Shao et al., Mol. Phys. (published online 2014); available at http://www. tandfonline.com/doi/abs/10.1080/00268976.2014.952696#.VKVQvWTF_ vQ. 19F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 73 (2012). 20G. Sierra and T. Nishino, Nucl. Phys. B 495, 505 (1997). 21I. P. McCulloch and M. Gulacsi, Aust. J. Phys. 53, 597 (2000). 22I. P. McCulloch and M. Gulacsi, Europhys. Lett. 57, 852 (2002). 23D. Zgid and M. Nooijen, J. Chem. Phys. 128, 14107 (2008). 24S. Wouters, W. Poelmans, P. W. Ayers, and D. Van Neck, Comput. Phys. Commun. 185, 1501 (2014). 25S. Sharma and G. K.-L. Chan, J. Chem. Phys. 136, 124121 (2012). 26A. Weichselbaum, Ann. Phys. 327, 2972 (2012). 27T. H. Dunning, J. Chem. Phys. 90, 1007 (1989). 28D. M. Brink and G. R. Satchler, Angular Momentum, Dover Books on Physics (Dover Publications, 1968). 29H. Morton, Group Theory and Its Application to Physical Problems, Oxford Library of the Physical Sciences (Clarendon Press, 1989). 30W. Weltner and R. J. Van Zee, Chem. Rev. 89, 1713 (1989). 31M. Martin, J. Photochem. Photobiol., A 66, 263 (1992).
J. Chem. Phys. 142, 024107 (2015) 32A.
Van Orden and R. J. Saykally, Chem. Rev. 98, 2313 (1998). D. Watts and R. J. Bartlett, J. Chem. Phys. 96, 6073 (1992). 34C. D. Sherrill and P. Piecuch, J. Chem. Phys. 122, 124104 (2005). 35D. I. Lyakh, V. V. Ivanov, and L. Adamowicz, J. Chem. Phys. 128, 074101 (2008). 36D. Datta, L. Kong, and M. Nooijen, J. Chem. Phys. 134, 214116 (2011). 37M. L. Abrams and C. D. Sherrill, J. Chem. Phys. 121, 9211 (2004). 38U. S. Mahapatra, S. Chattopadhyay, and R. K. Chaudhuri, J. Chem. Phys. 129, 024108 (2008). 39W. Jiang and A. K. Wilson, J. Chem. Phys. 134, 034101 (2011). 40A. D. Pradhan, H. Partridge, and C. W. Bauschlicher, J. Chem. Phys. 101, 3857 (1994). 41G. H. Booth, D. Cleland, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 135, 84104 (2011). 42D. Cleland, G. H. Booth, C. Overy, and A. Alavi, J. Chem. Theory Comput. 8, 4138 (2012). 43S. Ten-no, J. Chem. Phys. 138, 164126 (2013). 44Y. Ohtsuka and S. Ten-no, AIP Conference Proceedings 1456, 97 (2012). 45L. Bytautas and K. Ruedenberg, J. Chem. Phys. 122, 154110 (2005). 46G. Gidofalvi and D. A. Mazziotti, J. Chem. Phys. 122, 194104 (2005). 47S. Shaik et al., Nat. Chem. 4, 195 (2012). 48S. L. Altmann and P. Herzig, Point-Group Theory Tables (Clarendon Press, Oxford, 1994). 49G. H. Booth, D. Cleland, A. Alavi, and D. P. Tew, J. Chem. Phys. 137, 164112 (2012). 50S. Wouters, W. Poelmans, P. W. Ayers, and D. V. Neck, Comput. Phys. Commun. 185, 1501 (2014). 51S. Sharma, T. Yanai, G. H. Booth, C. J. Umrigar, and G. K.-L. Chan, J. Chem. Phys. 140, 104112 (2014). 52T. Yanai and T. Shiozaki, J. Chem. Phys. 136, 84107 (2012). 33J.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Wed, 11 Mar 2015 07:24:40
