Dynamically weighted multireference perturbation theory: Combining the advantages of multi-state and state-averaged methods

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 ¹Σ⁺ symmetry.^1,7 Multi-state (MS) MRPTs that follow the “diagonalize then perturb then diagonalize” route^1–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) formalism¹² 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) method^17,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),¹⁷ square of hyperbolic secant,¹⁷ and splines.¹⁸ 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.¹⁹ 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).¹⁴ 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

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.²⁰ A generic k-body cluster operator $T^α,k(s)$ is defined as

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 ordered²¹ 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)

where $R^(s)$ is the “source operator” and the superscript “N” indicates the nondiagonal part²² 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.¹⁴ In the unrelaxed variant of the DSRG-PT2 (u-DSRG-PT2), the energy of each state is computed as the expectation value

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

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

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,

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

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

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.¹⁵ 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öhn²⁶ 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

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

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,

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

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−Âα$⁠,²⁷ reaching in the end a generalized eigenvalue problem

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

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,

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

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 ¹Σ⁺ states of LiF.^{1–3,10–12} A two-state computation shows that both methods correctly predict an avoided crossing around r_Li–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-CASPT2²⁹ and QD-NEVPT2.¹² It has been shown that the hump can be cured by an extended multi-state treatment.¹¹ 

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 (r_Li–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 (E₁ + E₂) with respect to the SA-DSRG-PT2 values. Since E₁ + E₂ is equal to the trace of $Hαβeff$⁠, the large errors in E₁ + E₂ near r_Li–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 r_Li–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:

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

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

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

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α(ζ)$].²¹ 

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

The corresponding DW-DSRG-PT effective Hamiltonian is then equivalent to the MS-DSRG-PT formalism with state-specific diagonal elements.¹⁴ 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)¹² 

In this case, since ω⁽¹⁾ = ω⁽²⁾, 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 $Eh−2$⁠) 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 mE_h). 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 NH₃ along the set of reaction coordinates introduced by Truhlar and co-workers.³² In particular, we consider a trigonal planar configuration in which one of the N-H bond lengths (r₁) 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 (A₁ and B₂), 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 NH₃ computed with the SA-, MS-, and DW-DSRG-PT2 methods using three different values of ζ. As shown in Fig. 2(a), using ζ = 50 $Eh−2$ 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.

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 $Eh−2$⁠. 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 $Eh−2$⁠) 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.³³ Electron repulsion integrals were approximated using density fitting where the def2-TZVP-JKFIT³⁴ and def2-TZVP-RI³⁵ auxiliary basis sets were utilized for the RHF and DSRG-MRPT2³⁶ 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.³⁸ 

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.³⁹ 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 (¹E_2g/¹E_1u). All excitation energies are calculated using the ground-state energy of the corresponding theory.

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 ¹B_3u states of anthracene (2¹B_3u, 3¹B_3u). These two states are strongly coupled, and their ordering is predicted differently by CASCI and u-DSRG-PT2. As a result, when the third ¹B_3u state is included, the coupling between the second and third ¹B_3u 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 theory²⁸ and nonperturbative truncation schemes.¹⁴ 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, NH₃, 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

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 $(Δα)abij≠0$ and $(Δβ)abij≠0$] and compute the difference in the value of the amplitudes for references $Ψ0α$ and $Ψ0β$

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

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 H^eff, it is convenient to normal order the operators $H¯α(s)$ with respect to the true physical vacuum,²⁸ 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

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α$ (⁠$γ¯α$⁠)

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