We introduce two new approaches to compute near-degenerate electronic states based on the driven similarity renormalization group (DSRG) framework. The first approach is a unitary multi-state formalism based on the DSRG (MS-DSRG), whereby an effective Hamiltonian is built from a set of state-specific solutions. The second approach employs a dynamic weighting parameter to smoothly interpolate between the multi-state and the state-averaged DSRG schemes. The resulting dynamically weighted DSRG (DW-DSRG) theory incorporates the most desirable features of both multi-state approaches (ability to accurately treat many states) and state-averaged methods (correct description of avoided crossings and conical intersections). We formulate second-order perturbation theories (PT2) based on the MS- and DW-DSRG and study the potential energy curves of LiF, the conical intersection of the two lowest singlet states of NH3, and several low-lying excited states of benzene, naphthalene, and anthracene. The DW-DSRG-PT2 predicts the correct avoided crossing of LiF and avoids artifacts produced by the corresponding state-specific and multi-state theories. Excitation energies of the acenes computed with the DW-DSRG-PT2 are found to be more accurate than the corresponding state-averaged values, showing a small dependence on the number of states computed.

Multireference perturbation theories (MRPTs) based on a complete active space (CAS) reference are often the methods of choice for studying near-degenerate electronic states.1–6 Among the various flavors of MRPTs, state-specific schemes belong to the “diagonalize then perturb” philosophy and are well suited to study electronic states that are energetically separated. However, since these approaches neglect the coupling between different perturbed references, they yield qualitatively incorrect potential energy surfaces (PESs) close to points of near degeneracy. A classical example is the dissociation curve of the LiF molecule where the non-crossing rule is violated for the lowest two singlet states of 1Σ+ symmetry.1,7 Multi-state (MS) MRPTs that follow the “diagonalize then perturb then diagonalize” route1–3,8,9 circumvent this problem by explicitly accounting for the coupling between perturbed references. However, MS-MRPTs may still yield incorrect PESs near degenerate regions unless special care is taken to enforce that the zeroth-order Hamiltonian is invariant with respect to unitary rotations of the model space, as done in extended MS (XMS) approaches.10,11

Recently, two of us have proposed a different route to compute multiple near-degenerate excited states based on a state-averaged (SA) formalism12 derived from the multireference driven similarity renormalization group (DSRG).13,14 The SA-DSRG approach follows the “diagonalize then perturb then diagonalize” strategy, in which a single unitary transformation is applied to the Hamiltonian.15,16 This formalism leads to a Hermitian effective Hamiltonian, which, upon diagonalization, gives the energies of a manifold of excited states. The original SA-DSRG formalism has several appealing features: (i) it can deal with near-degenerate electronic states, (ii) it can treat several excited states at the cost of one state-specific computation, and (iii) it includes reference relaxation effects. The SA-DSRG is particularly suited to compute near-degenerate electronic states since it treats dynamical correlation effects in a state-averaged way, avoiding the need for the XMS formalism. However, a state-averaged formalism is likely to degrade in accuracy as the number of states computed increases or when treating states with very different electronic character (e.g., neutral and ionic).

The main goal of this study is to formulate a generalized SA-DSRG scheme that in the limit of well-separated electronic states is equivalent to a MS treatment while for near-degenerate states corresponds to a SA approach. Our approach dynamically interpolates between these two cases, and henceforth, we refer to it as “dynamically weighted” DSRG (DW-DSRG). We have also formulated, for the first time, a multi-state DSRG (MS-DSRG) approach based on unitary transformations. Furthermore, we have investigated a self-consistent variant of the MS-DSRG (sc-MS-DSRG) in which the basis of states is optimized in the presence of dynamical correlation effects.

We note that our approach is similar to the dynamically weighted CAS self-consistent field (DW-CASSCF) method17,18 in which the weight of a state in a state-averaged CASSCF calculation is adjusted depending on its energy. The DW-CASSCF approach successfully avoids energy discontinuities that arise from the crossing of states in and outside the model space.17,18 Different weighting functions have also been proposed, including Gaussian damping (as in this work),17 square of hyperbolic secant,17 and splines.18 However, the purpose of dynamically weighting in our treatment of dynamical electron correlation is to improve the description of states that become near degenerate, and therefore, it is substantially different from that of the DW-CASSCF. Another method relevant to our work is the dynamically weighted polarizable quantum mechanics/molecular mechanics scheme by Hagras and Glover, introduced to avoid artifacts and/or poor convergence near state crossings.19 This approach tunes the quantum mechanical part of the energy between state-specific and state-averaged values using the energy separation between the state of interest and all other states.

We will begin by summarizing the main features of the state-specific driven similarity renormalization group perturbation theory (DSRG-PT).14 In the DSRG-PT approach, we consider a set (or ensemble) of n zeroth-order states, E0{Ψ0α,α=1,2,,n}, where the states Ψ0α are CASSCF or CAS configuration interaction (CASCI) wave functions assumed to be orthonormal. For each ensemble state Ψ0α, the Hamiltonian is transformed via a unitary operator Ûα(s) that decouples Ψ0α from its internally contracted excited configurations. Note that since this transformation is assumed to be different for each state, we label the operator Ûα(s) with the subscript α. The unitary operator Ûα(s) is written in an exponential form as eÂα(s), where Âα(s) is an anti-Hermitian operator that depends on a continuous time-like variable s ∈ [0, ). Consequently, the unitarily transformed DSRG Hamiltonian for state Ψ0α [H¯α(s)] is given by

H¯α(s)=eÂα(s)ĤeÂα(s).
(1)

It is convenient to express this operator as Âα(s)=T^α(s)T^α(s), where T^α(s)=T^α,1(s)+T^α,2(s)+ is an s-dependent generalization of the single-reference cluster operator.20 A generic k-body cluster operator T^α,k(s) is defined as

T^α,k(s)=1(k!)2ijHabPtabij{âijab}α,
(2)

where the orbital spaces H (holes) and P (particles) span the core/active and the active/virtual spaces, respectively, and the notation {Ô}α is used to denote the operator Ôα normal ordered21 with respect to the state Ψ0α, i.e., Ψ0α|{Ô}α|Ψ0α=0.

The anti-Hermitian operator Âα(s) is implicitly determined by a nonlinear many-body equation (DSRG flow equation)

[H¯α(s)]N=R^α(s),
(3)

where R^(s) is the “source operator” and the superscript “N” indicates the nondiagonal part22 of H¯α(s). The flow equation realizes a renormalization transformation of the Hamiltonian. That is, as s goes from zero to infinity, the transformation gradually zeroes the nondiagonal elements of H¯(s), starting with elements that correspond to large energy denominators and then proceeding to those with zero denominators. This feature allows DSRG-based methods to avoid the intruder state problem.13,14,23–25 In the second-order DSRG-PT method (DSRG-PT2), the DSRG transformation and flow equations [Eqs. (1) and (3)] are realized approximately up to second order in perturbation theory.14 In the unrelaxed variant of the DSRG-PT2 (u-DSRG-PT2), the energy of each state is computed as the expectation value

Euα=Ψ0α|H¯α(s)|Ψ0α.
(4)

Reference relaxation effects can be introduced in the state-specific DSRG-PT2,14 but these are not relevant to this work.

In the SA-DSRG-PT2,12 a single unitary transformation of the Hamiltonian is performed for all the states in the ensemble

H¯(s)=eÂ(s)ĤeÂ(s),
(5)

where Â(s) is a state-averaged operator defined in a way analogous to the state-specific operators Âα(s). To produce a transformation that accounts for dynamical correlation in an average way for all the states in the ensemble, an equal-weight density matrix (ρ^) is introduced,

ρ^=α=1nωα|Ψ0αΨ0α|,
(6)

and all operators are normal ordered with respect to ρ^, i.e., {Ô}ρ is normal ordered if Tr(ρ^{Ô}ρ)=0. In the SA-DSRG-PT2, only one set of DSRG equations [Eq. (3)] is solved, using operators normal ordered with respect to ρ^. The contracted SA-DSRG-PT2 (SA-DSRG-PT2c) energy is then obtained by forming the matrix element of H¯(s) in the basis of ensemble states, H¯αβ(s)=Ψ0α|H¯(s)|Ψ0β, and diagonalizing it.

In this section, we introduce a multi-state formulation of the DSRG-PT2 (MS-DSRG-PT2). To the best of our knowledge, multi-state formalisms based on unitary transformations have been largely unexplored. Therefore, we briefly discuss some of the challenges we encountered in the development of the present approach. The DSRG transformation [Eq. (1)] suggests that to each zeroth-order state Ψ0αE0, we may associate a correlated state (Ψ̃α) defined as

|Ψ̃α=eÂα|Ψ0α.
(7)

In the MS-DSRG-PT2, we postulate that the wave function for excited state γγ) is written as a linear combination of correlated states as

|Θγ=β=1neÂβ|Ψ0βCβγ,
(8)

where the operators Âβ are obtained by solving the state-specific DSRG equations and Cβγ are variationally optimized coefficients. This wave function form is similar to the ansatz proposed by Jeziorski and Monkhorst to formulate the state-universal multireference coupled cluster method.15 However, the MS-DSRG differs from the Jeziorski–Monkhorst approach in three ways: (i) it uses a unitary transformation, (ii) it employs a basis of CASSCF/CASCI states instead of model space determinants, and (iii) it uses state-specific cluster operators. The multi-state formalism of Aoto and Köhn26 also considers linear combinations of correlated states; however, these authors consider an exponential transformation that contains only excitation operators and the zeroth-order states in E0 are not assumed to be mutually orthogonal.

The states |Ψ̃α are normalized

Ψ̃α|Ψ̃α=Ψ0α|eÂαeÂα|Ψ0α=Ψ0α|Ψ0α=1;
(9)

however, pairs of states |Ψ̃α, |Ψ̃β with αβ are not orthogonal. It is convenient to introduce the overlap matrix (Sαβ) between correlated states

Sαβ=Ψ0α|eÂαeÂβ|Ψ0β=Ψ0α|eX^αβ|Ψ0β.
(10)

In Eq. (10), we have introduced the operator X^αβ(s)=log[eÂα(s)eÂβ(s)]=Âβ(s)Âα(s)+12[Âβ(s),Âα(s)]+, which can be evaluated using the Baker–Campbell–Hausdorff formula. Note that contrary to the case of coupled cluster theory, the correlated states do not satisfy the intermediate normalization condition, that is,

Ψ0α|Ψ̃α=Ψ0α|eÂα|Ψ0α1, ifÂα0.
(11)

Imposing orthonormality of the correlated wave functions leads to the following condition:

Θγ|Θλ=αβ=1n(Cαγ)*SαβCβλ=δγλ,
(12)

which implies that the coefficient vectors Cβγ are not orthonormal.

To obtain equations for the MS-DSRG approach, we insert |Θγ⟩ into the Schrödinger equation and left-project with Ψ0α|eÂα,27 reaching in the end a generalized eigenvalue problem

β=1nHαβeffCβγ=Eγβ=1nSαβCβγ,
(13)

where the MS-DSRG effective Hamiltonian Hαβeff is given by

Hαβeff=Ψ0α|eÂα(s)ĤeÂβ(s)|Ψ0β=12Ψ0α|H¯α(s)eX^αβ(s)+eX^αβ(s)H¯β(s)|Ψ0β.
(14)

As shown in  Appendix A, when X^αβ(s) is evaluated using first-order operators [Âα(1)(s) and Âβ(1)(s)], then the quantity Ŷαβ(s)=exp[X^αβ(1)(s)]=1+O[Ŷαβ(2)(s)], where Ŷαβ(2)(s) indicates a second order quantity. Consequently, in the evaluation of the effective Hamiltonian and overlap matrices, we retain contributions up to order two and one, respectively,

Hαβeff,[2](s)=12Ψ0α|H¯α[2](s)+H¯β[2](s)|Ψ0β,
(15)
Sαβ[1]=δαβ,
(16)

where superscript “[2]” indicates the sum of zeroth-, first-, and second-order terms. Note that the DSRG transformed Hamiltonian H¯α(s) in Eq. (15) generally contains three- and higher-body contributions

H¯α(s)=H¯α,0(s)+H¯α,1(s)+H¯α,2(s)+,
(17)

where H¯α,k(s) is a k-body operator normal ordered with respect to the reference Ψ0α. In our previous studies,12,14,28 we kept at most the two-body terms of H¯α(s) and this approximation is also assumed in this work. However, we also report, for the first time, results where the elements Hαβeff contain contributions from three-body components of H¯α(s) and H¯β(s). The evaluation of the matrix elements Hαβeff with H¯α(s) truncated to two- and three-body operators is discussed in  Appendix B.

In Fig. 1(a), we show the performance of the SA- and MS-DSRG-PT2 methods on a common testbed for excited-state theories: the avoided crossing of the lowest two 1Σ+ states of LiF.1–3,10–12 A two-state computation shows that both methods correctly predict an avoided crossing around rLi–F = 6.5 Å. However, the MS-DSRG-PT2 excited state curve displays a small “hump” around 5.5 Å, in the region where the SA2-CASSCF(6e, 7o) reference predicts an avoided crossing. This hump is absent from the contracted SA-DSRG-PT2 (SA-DSRG-PT2c) approach but was previously observed for MS-CASPT229 and QD-NEVPT2.12 It has been shown that the hump can be cured by an extended multi-state treatment.11 

FIG. 1.

(a) The lowest two 1Σ+ states of lithium fluoride near the avoided crossing region. The circles show the average energy of these two states at the bond length with the minimum energy gap. (b) Energy differences of the ground (solid lines) and excited (dashed lines) states with respect to the corresponding SA-DSRG-PT2c energies. (c) Difference of the sum of the ground and excited state energies (E1 + E2) with respect to the SA-DSRG-PT2c value. The DW-DSRG-PT2 curves are obtained by setting ζ = 50 Eh2, while the MS-DSRG-PT2 results are obtained by setting ζ to 1014Eh2 in the DW-DSRG-PT2. The SA-DSRG-PT2c method is equivalent to DW-DSRG-PT2 with ζ = 0, and SA-DSRG-PT2c(H3) includes the three-body term H¯α,3(s) in the DSRG transformed Hamiltonian. This example employed a CASSCF reference with six-electron-in-seven-orbital averaged over the lowest two states [SA2-CASSCF(6e, 7o)] using the aug(F)-cc-pVTZ basis set. This basis set consists of the cc-pVTZ basis set for Li30 and the aug-cc-pVTZ basis set for F.31 All 1s-like orbitals of Li, F, and C were not correlated in the DSRG-PT2 part of the computation, and the flow parameter was set to 0.5 Eh2.

FIG. 1.

(a) The lowest two 1Σ+ states of lithium fluoride near the avoided crossing region. The circles show the average energy of these two states at the bond length with the minimum energy gap. (b) Energy differences of the ground (solid lines) and excited (dashed lines) states with respect to the corresponding SA-DSRG-PT2c energies. (c) Difference of the sum of the ground and excited state energies (E1 + E2) with respect to the SA-DSRG-PT2c value. The DW-DSRG-PT2 curves are obtained by setting ζ = 50 Eh2, while the MS-DSRG-PT2 results are obtained by setting ζ to 1014Eh2 in the DW-DSRG-PT2. The SA-DSRG-PT2c method is equivalent to DW-DSRG-PT2 with ζ = 0, and SA-DSRG-PT2c(H3) includes the three-body term H¯α,3(s) in the DSRG transformed Hamiltonian. This example employed a CASSCF reference with six-electron-in-seven-orbital averaged over the lowest two states [SA2-CASSCF(6e, 7o)] using the aug(F)-cc-pVTZ basis set. This basis set consists of the cc-pVTZ basis set for Li30 and the aug-cc-pVTZ basis set for F.31 All 1s-like orbitals of Li, F, and C were not correlated in the DSRG-PT2 part of the computation, and the flow parameter was set to 0.5 Eh2.

Close modal

To further analyze the origin of this hump in the MS-DSRG-PT2, in Fig. 1(b), we plot the energy deviation with respect to SA-DSRG-PT2. Note that in correspondence to the hump (rLi–F = 5.5 Å), the MS-DSRG-PT2 energy deviation is positive for both states. In Fig. 1(c), we plot the deviation of the sum of the ground and excited states energies (E1 + E2) with respect to the SA-DSRG-PT2 values. Since E1 + E2 is equal to the trace of Hαβeff, the large errors in E1 + E2 near rLi–F = 5.5 Å suggest that it is necessary to improve the diagonal elements of the effective Hamiltonian to avoid the hump. One possibility is using zeroth-order states that display the avoided crossing at the correct bond length. This condition can be achieved by formulating a self-consistent MS-DSRG-PT2 (sc-MS-DSRG-PT2) scheme, whereby at each iteration, Hαβeff is built and its eigenvectors are used to define a new set of zeroth-order states Ψ0α, repeating this process until self consistency is reached. Unfortunately, as shown in Fig. 1(c), even the sc-MS-DSRG-PT2 scheme shows a hump, although in this case it is found shifted to around rLi–F = 6.5 Å. In summary, the persistence of the hump, even after improving the zeroth-order states, suggests that in the multi-state formalism, the diagonal elements of Hαβeff are very sensitive to the sudden switch in wave function character that occurs near an avoided crossing.

In DW-DSRG, we propose to evaluate the energy of a manifold of excited states using the MS-DSRG effective Hamiltonian formalism and to determine the operators Âα(s) using a scheme that interpolates between a state-specific and a state-averaged one. One practical consequence of this choice is that near an avoided crossing, the diagonal elements of the effective Hamiltonian become less dependent on the zeroth-order reference states. We will first discuss the procedure to dynamically weight electronic states, and then, we will discuss the construction of the effective Hamiltonian.

We begin by generalizing the MS-DSRG formalism and associating to each target state α an ensemble density matrix (ρ^α) defined by the weights (ωβα) as follows:

ρ^α=β=1nωβα|Ψ0βΨ0β|.
(18)

This expression is a state-specific generalization of the SA-DSRG density matrix [Eq. (6)].

For a zeroth-order state Ψ0α well separated energetically from other states in the ensemble, we require that the density matrix is pure, that is, ρ^α|Ψ0αΨ0α|. On the contrary, when the state Ψ0α is degenerate with one or more other states, we require the density matrix ρ^α to be an equal-weight average of this set of degenerate states. One way to satisfy these two conditions is to define the weight for target state α using a Gaussian function

ωβα(ζ)=eζ(Eα(0)Eβ(0))2γ=1neζ(Eα(0)Eγ(0))2,
(19)

where Eα(0) is the zeroth-order energy of state α and ζ is a parameter. When interpreted as a function of Eβ(0), ωβα(ζ) is peaked around Eα(0). Furthermore, it satisfies the normalization condition βωβα(ζ)=1. Similarly, we define the average k-body density matrix for target state α as

γ¯kα(ζ)=α=1nωβα(ζ)γkββ.
(20)

The quantities γkββ are diagonal elements of the k-body transition density matrices (γkαβ), defined as

[γkαβ]rspq=Ψ0α|âpâqâsâr|Ψ0β.
(21)

In the DW-DSRG-PT, we compute the diagonal matrix elements of the effective Hamiltonian by solving the DSRG flow equation using many-body operators normal ordered with respect to the ensemble density matrix ρ^α. Practically, this amounts to replacing the state-specific reduced-density matrices used in the generalized Wick theorem with the dynamically weighted ones [γ¯kα(ζ)].21 

To illustrate the properties of the DW-DSRG-PT scheme, we consider an ensemble of two reference states, Ψ01 and Ψ02. When the energy separation of these two states is large (|E1(0)E2(0)|), the weights used to define the target state densities are unit vectors

ω(1)10,ω(2)01.
(22)

The corresponding DW-DSRG-PT effective Hamiltonian is then equivalent to the MS-DSRG-PT formalism with state-specific diagonal elements.14 In the limit of degenerate zeroth-order states (|E1(0)E2(0)|0), the weights are equivalent to those used in the equal-weighted contracted SA-DSRG-PT (SA-DSRG-PT2c)12 

ω(1)1212,ω(2)1212.
(23)

In this case, since ω(1) = ω(2), the amplitude equations for both states are identical, which implies that Âα(s) = Âβ(s). Consequently, the correlated states Ψ̃α and Ψ̃β are orthogonal and the overlap matrix Sαβ reduces to the identity. For intermediate cases, the density matrices will smoothly interpolate between these two limits. Note that when the energy separation satisfies |Eα(0)Eβ(0)|ζ1/2, then weights for target state α are approximately diagonal, i.e., ωβα(ζ)δαβ (unless two or more states are exactly degenerate), and the DW-DSRG-PT approach approximates the MS-DSRG-PT. Therefore, the energy ζ−1/2 may be interpreted as a cutoff that determines the transition between the MS and the SA treatment. Note that the DW-DSRG-PT2 is equivalent to the MS-DSRG-PT2 in the limit of ζ.

As shown in Fig. 1(a), the DW-DSRG-PT2 (with ζ = 50 Eh2) correctly reproduces the avoided crossing of LiF and the curves are free from the hump observed for the MS-DSRG-PT2 [see Fig. 1(b)]. Figure 1(b) also demonstrates that the three-body term H¯α,3(s) in the DSRG transformed Hamiltonian [see Eq. (17)] can be safely ignored due to its negligible energy contribution (<0.5 mEh). Thus, the results presented in the following text were computed without the three-body contributions in the effective Hamiltonian.

To test the DW-DSRG-PT2 ability to describe near-degenerate states, we consider the conical intersection between the lowest two singlet states of NH3 along the set of reaction coordinates introduced by Truhlar and co-workers.32 In particular, we consider a trigonal planar configuration in which one of the N-H bond lengths (r1) is allowed to vary, while all angles and the other two N–H bond lengths (1.039 Å) are fixed. In this configuration, the lowest two singlet states belong to different irreducible representations (A1 and B2), and as such, the effective Hamiltonians for all of the theories considered here are diagonal. However, due to the different way the diagonal elements are built, the MS-DSRG-PT2, SA-DSRG-PT2c, and DW-DSRG-PT2 methods will yield different energies.

In Fig. 2, we show the potential energy curves of the first two states of NH3 computed with the SA-, MS-, and DW-DSRG-PT2 methods using three different values of ζ. As shown in Fig. 2(a), using ζ = 50 Eh2 yields DW-DSRG-PT2 results that closely follow the SA-DSRG-PT2c curves near the conical intersection. Note also that the MS-DSRG-PT2 curves parallel the SA-DSRG-PT2c results.

FIG. 2.

(a) Potential energy curves for the lowest two singlet states of planar NH3 along the r1 coordinate computed using various DSRG-PT2 schemes. (b) Energy deviations of the A1 (solid lines) and B2 (dashed lines) states obtained by SA- and DW-DSRG-PT2 methods compared to the MS-DSRG-PT2 values. The vertical dotted line indicates the location of the conical intersection obtained by the equal-weighted SA2-CASSCF(8e, 7o)/aug-cc-pVTZ theory.

FIG. 2.

(a) Potential energy curves for the lowest two singlet states of planar NH3 along the r1 coordinate computed using various DSRG-PT2 schemes. (b) Energy deviations of the A1 (solid lines) and B2 (dashed lines) states obtained by SA- and DW-DSRG-PT2 methods compared to the MS-DSRG-PT2 values. The vertical dotted line indicates the location of the conical intersection obtained by the equal-weighted SA2-CASSCF(8e, 7o)/aug-cc-pVTZ theory.

Close modal

Since the DW-DSRG-PT2 approach is essentially interpolating between the MS- and SA-DSRG schemes, a sudden switch between the non-degenerate and degenerate case can introduce artifacts in the potential energy curves. This problem can be seen in Fig. 2(b), where we plot the DW-DSRG-PT2 curves computed with ζ = 500 and 5000 Eh2. Using larger values of ζ leads to “wiggles” in the potential energy curves that originate from the sudden change in the weights of the zeroth-order states that define the state-averaged density matrix [Eq. (20)]. In particular, it can be seen that away from the conical intersection, the DW-DSRG-PT2 (ζ = 5000 Eh2) curves approach the corresponding MS-DSRG-PT2 results, while near the SA2-CASSCF conical intersection (located around 1.95 Å), they approach the SA-DSRG-PT2c curves. This example shows that one potential issue with the DW scheme is selecting an appropriate value of the parameter ζ. More extensive benchmarks of conical intersections will be required to identify an optimal value of ζ that minimizes errors and at the same time does not introduce artifacts in the description of conical intersections.

Next, we examine the dependence of DW-DSRG-PT2 excitation energies with respect to the number of computed roots for low-lying singlet states of benzene, naphthalene, and anthracene. For each molecule, we employed a CASCI reference that treats all π and π* orbitals as active orbitals, leading to CAS(6, 6), CAS(10, 10), and CAS(14, 14) active spaces for benzene, naphthalene, and anthracene, respectively. Molecular orbitals were computed using restricted Hartree–Fock (RHF) theory and the def2-TZVP basis set.33 Electron repulsion integrals were approximated using density fitting where the def2-TZVP-JKFIT34 and def2-TZVP-RI35 auxiliary basis sets were utilized for the RHF and DSRG-MRPT236 computations, respectively. Following Ref. 37, the ground-state geometries of benzene, naphthalene, and anthracene were optimized by second-order Møller–Plesset perturbation theory with the 6-31G* basis set.38 

In Figs. 3–5, we compare the vertical excitation energies of SA-, MS-, and DW-DSRG-PT2 with various number of averaged states against those of state-specific unrelaxed DSRG-PT2 (u-DSRG-PT2). In our computations on benzene and naphthalene, we add states according to the order of CASPT2 results as reported in Ref. 37, while for anthracene, we follow the ordering from multireference Møller–Plesset perturbation theory.39 Consequently, this order may not follow that obtained from u-DSRG-PT2. In the case of benzene, we avoid symmetry breaking by including all states of a degenerate irreducible representation (1E2g/1E1u). All excitation energies are calculated using the ground-state energy of the corresponding theory.

FIG. 3.

Deviations of benzene vertical excitation energies (in eV) computed using the SA-, MS-, and DW-DSRG-PT2 methods relative to the corresponding state-specific u-DSRG-PT2 values (shown on the right). The MS- and DW-DSRG-PT2 values may be overlapped and difficult to resolve. The SA scheme is the SA-DSRG-PT2c variant in Ref. 12.

FIG. 3.

Deviations of benzene vertical excitation energies (in eV) computed using the SA-, MS-, and DW-DSRG-PT2 methods relative to the corresponding state-specific u-DSRG-PT2 values (shown on the right). The MS- and DW-DSRG-PT2 values may be overlapped and difficult to resolve. The SA scheme is the SA-DSRG-PT2c variant in Ref. 12.

Close modal
FIG. 4.

Deviations of naphthalene vertical excitation energies (in eV) computed using the SA-, MS-, and DW-DSRG-PT2 methods relative to the corresponding state-specific u-DSRG-PT2 values (shown on the right). The MS- and DW-DSRG-PT2 values may be overlapped and difficult to resolve. The SA scheme is the SA-DSRG-PT2c variant in Ref. 12.

FIG. 4.

Deviations of naphthalene vertical excitation energies (in eV) computed using the SA-, MS-, and DW-DSRG-PT2 methods relative to the corresponding state-specific u-DSRG-PT2 values (shown on the right). The MS- and DW-DSRG-PT2 values may be overlapped and difficult to resolve. The SA scheme is the SA-DSRG-PT2c variant in Ref. 12.

Close modal
FIG. 5.

Deviations of anthracene vertical excitation energies (in eV) computed using the SA-, MS-, and DW-DSRG-PT2 methods relative to the corresponding state-specific u-DSRG-PT2 values (shown on the right). The MS- and DW-DSRG-PT2 values may be overlapped and difficult to resolve. The SA scheme is the SA-DSRG-PT2c variant in Ref. 12.

FIG. 5.

Deviations of anthracene vertical excitation energies (in eV) computed using the SA-, MS-, and DW-DSRG-PT2 methods relative to the corresponding state-specific u-DSRG-PT2 values (shown on the right). The MS- and DW-DSRG-PT2 values may be overlapped and difficult to resolve. The SA scheme is the SA-DSRG-PT2c variant in Ref. 12.

Close modal

In general, the MS-DSRG-PT2 excitation energies are in very good agreement with those computed using the u-DSRG-PT2 method. Larger deviations are observed for the second and third 1B3u states of anthracene (21B3u, 31B3u). These two states are strongly coupled, and their ordering is predicted differently by CASCI and u-DSRG-PT2. As a result, when the third 1B3u state is included, the coupling between the second and third 1B3u states pushes the lower state down and the higher state up by roughly 0.15 eV.

The excitation energies from the SA-DSRG-PT2c approach are systematically off from the state-specific u-DSRG-PT2 values by 0.18 eV on average, and deviations as high as 0.24 eV are seen for high-lying states of naphthalene and anthracene. In all cases, the DW-DSRG-PT2 excitation energies lie in between the SA and MS values and do not show systematic deviations as large as those observed for the state-averaged approach. Contrary to our initial expectations, although the SA- and DW-DSRG-PT2 theories show a dependence of the excitation energies with respect to the number of averaged states (increase in error), this effect is not very pronounced for the examples considered here. For a given excited state, variations in the predicted excitation energies with the number of states computed are very stable and vary at most by 0.04 eV.

In this work, we have derived and implemented two multi-state extensions of the second-order perturbation theory based on the driven similarity renormalization group. The first approach is a multi-state DSRG-PT2 (MS-DSRG-PT2) scheme analogous to conventional MS formalisms,1–3,8,9 in which a second-order Hamiltonian is built from a basis of state-specific states. The second approach employs a new strategy in which the diagonal elements of the effective Hamiltonian are obtained from a state-averaged DSRG-PT2 scheme using densities with dynamically adjusted weights (DW-DSRG-PT2). In DW-DSRG-PT2, the weights that determine the average density are automatically determined using the energy separations of the zeroth-order states. This approach smoothly interpolates between a state-specific density (for states that are energetically well separated) and a state-averaged density (for near-degenerate states). This property allows the DW-DSRG-PT2 to avoid artifacts displayed by MS methods near avoided crossing and conical intersections. At the same time, for well-separated states, the DW approach maintains the state-specific character of the MS formalism and avoids the danger of treating states with different chemical nature (e.g., diradical, charge separated, Rydberg) in an averaged way. These benefits are illustrated for two examples: the avoided crossing of LiF and the low-lying excited states of benzene, naphthalene, and anthracene. In the avoided crossing region of LiF, the MS-DSRG-PT2 curves show wiggles, an artifact that is avoided in the DW-DSRG-PT2 treatment. Excitation energies for the three acenes computed at the DW-DSRG-PT2 level are rather insensitive to the number of states considered. By contrast, a state-averaged treatment (SA-DSRG-PT2c) shows more pronounced systematic deviations from state-specific excitation energies and a slight degradation of the low-lying states when more states are computed.

In our opinion, the DW approach is an interesting alternative to extended multi-state (XMS) MRPTs,10,11 which also avoid some of the artifacts of the MS formalism. Both approaches try to improve the standard MS description by a suitable modification of the effective Hamiltonian. In the XMS scheme, the description of conical intersections is improved by introducing a zeroth-order Hamiltonian that is invariant to rotations among the CASCI solutions. The DW-DSRG-PT2 approach follows a different route, whereby dynamical correlation effects are treated appropriately depending if a state is energetically isolated or it belongs to a group of near-degenerate states.

The dynamically weighted scheme for excited states introduced here is not specific to the DSRG-PT2 method. The formalism may be straightforwardly extended to more accurate treatments of electron correlation such as DSRG third-order perturbation theory28 and nonperturbative truncation schemes.14 With small modifications, the DW scheme should be also applicable to other MRPTs, such as the CASPT2 approach. In this case, it would be interesting to examine the relative performance of the XMS vs the DW extensions of CASPT2. Finally, we point out that it may be advantageous to combine the DW-DSRG-PT and DW-CASSCF schemes to obtain smooth potential energy surfaces where conventional treatments based on equally weighted SA-CASSCF fail.

See supplementary material for the energies of LiF, NH3, benzene, naphthalene, and anthracene.

C.L. and F.A.E. were supported by the U.S. Department of Energy under Award No. DE-SC0016004 and by a Research Fellowship of the Alfred P. Sloan Foundation. R.L. acknowledges support by the Swedish Research Council (Grant No. 2016-03398), the National Science Foundation (Grant No. CHE-1464862), and the Wenner-Gren Foundation.

In this section, we show that when the operator X^αβ(s)=log[eÂα(s)eÂβ(s)] is evaluated using first-order operators [Âα(1)(s) and Âβ(1)(s)], exp[X^αβ(1)(s)]1 up to first order in perturbation theory. Expanding the operator Ŷαβ(s)=exp[X^αβ(1)(s)], we can write

Ŷαβ(s)=exp[X^αβ(1)(s)]=1+Âβ(1)(s)Âα(1)(s)+O[Ŷαβ(2)(s)],
(A1)

where O[Ŷαβ(2)(s)] collects all terms of order two and higher. Next, we show that the difference Âβ(1)(s)Âα(1)(s) is a second order quantity. For convenience, we take the limit s [assuming the Møller–Plesset denominators (Δα)abij0 and (Δβ)abij0] and compute the difference in the value of the amplitudes for references Ψ0α and Ψ0β

(tβ)abij,(1)(tα)abij,(1)=vabij,(1)1(Δβ)abij1(Δα)abij=vabij,(1)(Δα)abij(Δβ)abij(Δα)abij(Δβ)abij,
(A2)

where vabij,(1)=ϕaϕb|  |ϕiϕj are antisymmetrized two-electron integrals. To estimate the order of the quantity (Δα)abij(Δβ)abij, we note that the difference between the orbital energies of the two references (ϵα)p(ϵβ)p can be written as

(ϵα)p(ϵβ)p=rsvprps,(1)(γαα)sr(γββ)sr.
(A3)

Consequently, (Δα)abij(Δβ)abij is proportional to vprps,(1) times the difference γααγββ, giving this term an overall order of one (assuming the density matrices are zeroth-order quantities). When (Δα)abij(Δβ)abij is introduced in Eq. (A2), we find out that the total order of the amplitude difference (tβ)abij,(1)(tα)abij,(1) is equal to two, which implies that Âβ(1)(s)Âα(1)(s) is a second-order quantity.

In this section, we discuss the procedure used to evaluate the matrix elements of the effective Hamiltonian Hαβeff [see Eq. (15)]. To compute Hαβeff, we first solve the DSRG equation [Eq. (3)] for each state and evaluate H¯α(s) truncating the equations to the desired level (two- or three-body). By default, the operators H¯α(s) obtained by solving the DSRG equation are normal ordered with respect to the their corresponding density matrix ρ^α. To evaluate the matrix elements of Heff, it is convenient to normal order the operators H¯α(s) with respect to the true physical vacuum,28 i.e., in terms of bare creation and annihilation operators. After this step, a generic contribution of the form Ψ0α|H¯α(s)|Ψ0β can be obtained by contracting the vacuum-ordered H¯α(s) operators with the transition density matrices (γαβ). To illustrate this point, we will evaluate the contribution of the one-body operator H¯α,1(s) to Hαβeff. This matrix element is given by

Ψ0α|H¯α,1(s)|Ψ0β=ijHΨ0α|[H¯α(s)]ij{âji}ρα|Ψ0β
(B1)

and can be simplified by expressing the ρα-ordered operator {âqp}ρα in terms of bare operators using the definition of normal ordered operator, {âqp}ρα=âqp(γ¯α)qp. After introducing this expression in Eq. (B1), we obtain a formula that depends only on the transition density matrix (γαβ) and the average density matrix for state Ψ0α (γ¯α)

Ψ0α|H¯α,1(s)|Ψ0β=ijH[H¯α(s)]ijΨ0α|âji|Ψ0β(γ¯α)ji=ijH[H¯α(s)]ij(γαβγ¯α)ji.
(B2)

Expressions for higher-body contributions can be easily derived following the same procedure.

1.
H.
Nakano
,
J. Chem. Phys.
99
,
7983
(
1993
).
2.
J.
Finley
,
P.-Å.
Malmqvist
,
B. O.
Roos
, and
L.
Serrano-Andrés
,
Chem. Phys. Lett.
288
,
299
(
1998
).
3.
C.
Angeli
,
S.
Borini
,
M.
Cestari
, and
R.
Cimiraglia
,
J. Chem. Phys.
121
,
4043
(
2004
).
4.
Z.
Rolik
,
Á.
Szabados
, and
P. R.
Surján
,
J. Chem. Phys.
119
,
1922
(
2003
).
5.
R. K.
Chaudhuri
,
K. F.
Freed
,
G.
Hose
,
P.
Piecuch
,
K.
Kowalski
,
M.
Włoch
,
S.
Chattopadhyay
,
D.
Mukherjee
,
Z.
Rolik
,
Á.
Szabados
,
G.
Tóth
, and
P. R.
Surján
,
J. Chem. Phys.
122
,
134105
(
2005
).
6.
M. R.
Hoffmann
,
D.
Datta
,
S.
Das
,
D.
Mukherjee
,
Á.
Szabados
,
Z.
Rolik
, and
P. R.
Surján
,
J. Chem. Phys.
131
,
204104
(
2009
).
7.
J.-P.
Malrieu
,
J.-L.
Heully
, and
A.
Zaitsevskii
,
Theor. Chim. Acta
90
,
167
(
1995
).
8.
W.
Liu
and
M. R.
Hoffmann
,
Theor. Chem. Acc.
133
,
1481
(
2014
);
Y.
Lei
,
W.
Liu
, and
M. R.
Hoffmann
,
Mol. Phys.
115
,
2696
(
2017
).
9.
S.
Sharma
,
G.
Jeanmairet
, and
A.
Alavi
,
J. Chem. Phys.
144
,
034103
(
2016
).
10.
A. A.
Granovsky
,
J. Chem. Phys.
134
,
214113
(
2011
).
11.
T.
Shiozaki
,
W.
Győrffy
,
P.
Celani
, and
H.-J.
Werner
,
J. Chem. Phys.
135
,
081106
(
2011
).
12.
C.
Li
and
F. A.
Evangelista
,
J. Chem. Phys.
148
,
124106
(
2018
).
13.
F. A.
Evangelista
,
J. Chem. Phys.
141
,
054109
(
2014
).
14.
C.
Li
and
F. A.
Evangelista
,
J. Chem. Theory Comput.
11
,
2097
(
2015
);
[PubMed]
C.
Li
and
F. A.
Evangelista
,
J. Chem. Phys.
144
,
164114
(
2016
);
[PubMed]
C.
Li
and
F. A.
Evangelista
,
J. Chem. Phys.
148
,
079903
(
2018
).
[PubMed]
15.
B.
Jeziorski
and
H. J.
Monkhorst
,
Phys. Rev. A
24
,
1668
(
1981
).
16.
U. S.
Mahapatra
,
B.
Datta
, and
D.
Mukherjee
,
Chem. Phys. Lett.
299
,
42
(
1999
).
17.
M. P.
Deskevich
,
D. J.
Nesbitt
, and
H.-J.
Werner
,
J. Chem. Phys.
120
,
7281
(
2004
).
18.
W. J.
Glover
,
J. Chem. Phys.
141
,
171102
(
2014
).
19.
M. A.
Hagras
and
W. J.
Glover
,
J. Chem. Theory Comput.
14
,
2137
(
2018
).
20.
R. J.
Bartlett
and
M.
Musiał
,
Rev. Mod. Phys.
79
,
291
(
2007
).
21.
D.
Mukherjee
,
Chem. Phys. Lett.
274
,
561
(
1997
);
W.
Kutzelnigg
and
D.
Mukherjee
,
J. Chem. Phys.
107
,
432
(
1997
);
D.
Sinha
,
R.
Maitra
, and
D.
Mukherjee
,
Comput. Theor. Chem.
1003
,
62
(
2013
).
22.

The nondiagonal part of an operator  is defined by the components of  that couple a reference to its excited configurations.

23.
S.
Evangelisti
,
J. P.
Daudey
, and
J. P.
Malrieu
,
Phys. Rev. A
35
,
4930
(
1987
).
24.
J.
Paldus
,
P.
Piecuch
,
L.
Pylypow
, and
B.
Jeziorski
,
Phys. Rev. A
47
,
2738
(
1993
).
25.
K.
Kowalski
and
P.
Piecuch
,
Phys. Rev. A
61
,
052506
(
2000
).
26.
Y. A.
Aoto
and
A.
Köhn
,
J. Chem. Phys.
144
,
074103
(
2016
).
27.
Note that projecting the Schrödinger equation with Ψ0α| instead of a correlated state leads to the following less symmetric (and less appealing) eigenvalue equation
β=1nΨ0α|ĤeÂβ|Ψ0βCβγ=Eγβ=1nΨ0α|eÂβ|Ψ0βCβγ.
Contrary to the case of a pure excitation operator for which the matrix elements are manifestly connected, e.g.,
Ψ0α|ĤeT^β|Ψ0β=Ψ0α|eT^βĤeT^β|Ψ0β,
in the unitary case, it is not possible to eliminate the exponential term since
Ψ0α|ĤeÂβ|Ψ0β=Ψ0α|eÂβH¯β|Ψ0βΨ0α|H¯β|Ψ0β.
28.
C.
Li
and
F. A.
Evangelista
,
J. Chem. Phys.
146
,
124132
(
2017
);
[PubMed]
C.
Li
and
F. A.
Evangelista
,
J. Chem. Phys.
148
,
079902
(
2018
).
[PubMed]
29.
T.
Yanai
,
M.
Saitow
,
X.-G.
Xiong
,
J.
Chalupsky
,
Y.
Kurashige
,
S.
Guo
, and
S.
Sharma
,
J. Chem. Theory Comput.
13
,
4829
(
2017
).
30.
B. P.
Prascher
,
D. E.
Woon
,
K. A.
Peterson
,
T. H.
Dunning
, and
A. K.
Wilson
,
Theor. Chem. Acc.
128
,
69
(
2010
).
31.
R. A.
Kendall
,
T. H.
Dunning
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
32.
S.
Nangia
and
D. G.
Truhlar
,
J. Chem. Phys.
124
,
124309
(
2006
);
[PubMed]
Z. H.
Li
,
R.
Valero
, and
D. G.
Truhlar
,
Theor. Chem. Acc.
118
,
9
(
2007
).
33.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
34.
F.
Weigend
,
J. Comput. Chem.
29
,
167
(
2007
).
35.
C.
Hättig
,
Phys. Chem. Chem. Phys.
7
,
59
(
2005
).
36.
K. P.
Hannon
,
C.
Li
, and
F. A.
Evangelista
,
J. Chem. Phys.
144
,
204111
(
2016
).
37.
M.
Schreiber
,
M. R.
Silva-Junior
,
S. P. A.
Sauer
, and
W.
Thiel
,
J. Chem. Phys.
128
,
134110
(
2008
).
38.
M. M.
Francl
,
W. J.
Pietro
,
W. J.
Hehre
,
J. S.
Binkley
,
M. S.
Gordon
,
D. J.
DeFrees
, and
J. A.
Pople
,
J. Chem. Phys.
77
,
3654
(
1982
).
39.
Y.
Kawashima
,
T.
Hashimoto
,
H.
Nakano
, and
K.
Hirao
,
Theor. Chem. Acc.
102
,
49
(
1999
).

Supplementary Material