View PDF Version
 
DOI: 10.1039/D4FD00054D (Paper) Faraday Discuss., 2024, Advance Article

Spinless formulation of linearized adiabatic connection approximation and its comparison with the second order N-electron valence state perturbation theory

Yang Guo*a and Katarzyna Pernal*b
aQingdao Institute for Theoretical and Computational Sciences, Institute of Frontier Chemistry, School of Chemistry and Chemical Engineering, Shandong University, Qingdao, Shandong 266237, China. E-mail: yang.guo@sdu.edu.cn
bInstitute of Physics, Lodz University of Technology, ul. Wolczanska 217/221, 93-005 Lodz, Poland. E-mail: pernalk@gmail.com

Received 9th March 2024 , Accepted 2nd April 2024

First published on 3rd April 2024

Abstract

The adiabatic connection (AC) approximation, along with its linearized variant AC0, was introduced as a method of obtaining dynamic correlation energy. When using a complete active space self-consistent field (CASSCF) wave function as a reference, the AC0 approximation is considered one of the most efficient multi-reference perturbation theories. It only involves the use of 1st- and 2nd-order reduced density matrices. However, some numerical results have indicated that the excitation energies predicted by AC0 are not as reliable as those from the second-order N-electron valence state perturbation theory (NEVPT2). In this study, we develop a spinless formulation of AC0 based on the Dyall Hamiltonian and provide a detailed comparison between AC0 and NEVPT2 approaches. We demonstrate the components within the correlation energy expressions that are common to both methods and those unique to either AC0 or NEVPT2. We investigate the role of the terms exclusive to NEVPT2 and explore the possibility of enhancing AC0’s performance in this regard.

1 Introduction

The multi-configurational self-consistent field (MCSCF) method was proposed in the 1960s by several research groups, to capture electron correlations beyond the Hartree–Fock approximation.1–5 A special variant of MCSCF, the complete active space self-consistent field (CASSCF) method, developed by Roos and colleagues, is one of the most popular quantum chemistry methods for strongly correlated systems.6,7 In CASSCF calculations, active spaces, comprising n active electrons and m active spatial molecular orbitals (MOs), must be defined, typically referred to as CAS(n, m). Within the active space, full configuration interaction (FCI) calculations are performed. Thus, the size of the active space is the bottleneck in CASSCF calculations. As far as we know, without any approximation, the largest active space employed by CAS configuration interaction (CI) and FCI is CAS(22,22)[thin space (1/6-em)]8 and FCI(23,26),9 respectively. As a matter of fact, the conventional CASSCF/CASCI algorithms can only deal with systems with about 15 strongly correlated electrons and MOs. Note that the CAS method is also known as the full optimized reaction spaces (FORS) method.10,11

To extend CASSCF/CASCI methods to systems with more strongly correlated electrons or orbitals, various approximate FCI methods have been developed. In the 1970s, Malrieu and Peyerimhoff introduced their own selected CI (sCI) methods,12–14 while Goddard and colleagues developed the generalized valence bond (GVB) method.15 In the 1980s, Roos and colleagues extended CASSCF to constrained CAS and restricted active space SCF (RASSCF) methods.16,17 After 2000, Olsen and Gagliardi proposed the generalized active space (GAS) method,18–20 and Ivanic developed the occupation-restricted-multiple-active-space (ORMAS) method.21 In recent decades, sCI has been revived by various research groups worldwide, including Head-Gordon,22,23 Evangelista,24–26 Neese,27,28 Liu,29–34 and Sharma,35 among others.36–39 The fundamental idea behind sCI, RAS, GAS, and other similar methods is to reduce the dimension of the FCI space by introducing various approximations. Thus, these methods could be collectively referred to as reduced CI (RCI) methods. Alternatively, FCI quantum Monte Carlo (FCIQMC)40–44 and density matrix renormalization group (DMRG) theories45–51 have been developed as approximate FCI methods as well. All these methods have been used to address static correlations in systems with about one hundred strongly correlated electrons.

To capture the missing dynamic correlations of electrons, multi-reference (MR) perturbation theory (PT), MRCI, and MR coupled cluster (CC) methods based on RCI, FCIQMC, or DMRG wave functions have been developed. Using the DMRG reference, Yanai and colleagues developed CASPT2,52–54 MRCISD,55 as well as MRCC[thin space (1/6-em)]56 methods for strongly correlated systems requiring large active spaces. The CASPT2 methods based on RCI[thin space (1/6-em)]57,58 and DMRG[thin space (1/6-em)]59–61 are also reported by various groups. As one of the most popular MRPT2 methods, n-electron valence state perturbation theories (NEVPT2)62,63 based on DMRG,64–68 RCI,69 as well as FCIQMC,70 have been developed by many groups. Various MRCI and MRCC methods based on these approximate FCI wave functions have been developed by many groups as well.71,72 In the fully internally contracted (FIC) MRCISD or NEVPT2 method,62,63 the 5th-order reduced density matrices (RDMs) are required in principle. To reduce computational and storage costs, the rank reduction trick could be employed to reduce the requirement of 5th-order RDMs to 4th-order, provided that CASSCF or CASCI references are used.73 However, when utilizing approximate CAS references, the rank reduction trick could introduce unexpected intruder states.74 To achieve a stable FIC NEVPT2 or MRCISD algorithm with approximate CAS references, the 5th-order RDMs must be evaluated without additional approximations.75 Although the storage bottleneck of the 5th-order RDM could be circumvented by defining intermediates,68,76 their computational costs remain demanding.

To extend FIC MR dynamic correlation methods to systems with large active spaces, it is imperative to develop methods that require only low-order RDMs. Employing approximate cumulant expansion to lower the maximum rank of RDMs entering NEVPT2 equations leads to the appearance of false intruder states, which severely deteriorates the performance of the method.77 Evangelista and colleagues have developed various MR driven similarity renormalization group (DSRG) theories, including DSRG-PT and DSRG-CC, which only necessitate up to 3rd-order RDMs.78 Their extensions to systems with large active spaces have been explored as well.79,80 In 2018, Pernal proposed the adiabatic connection (AC) theory for MR.81,82 Based on transition density matrices from the extended random phase approximation (ERPA), AC theory can be considered a FIC MR method that only requires 1st- and 2nd-order RDMs. The linearized AC approximation (AC0) method could be considered as an MRPT2 method, which reverts to MP2 with a Hartree–Fock reference.82

Many applications of AC and AC0 indicate that they can deliver potential energy curves of ground states as accurately as CASPT2 and NEVPT2.83–85 However, for excited states, the AC0 approximation tends to overestimate the excitation energies.83 This is because, in the present AC0, the transition density matrices are computed from the quasiparticle ERPA, which lacks negative excitation transition density matrices. To address this issue without solving more complicated ERPA equations, the AC0D model has been developed to recover the missing correlation energies due to negative excitation transition density matrices for excited states.86 Most recently, Pernal and Veis extended AC0 to DMRG[thin space (1/6-em)]87 and explored further improvements in its accuracy by removing fixed-density matrix approximation.88

Although the theory of AC0 has been reported elsewhere, its relationship with conventional MRPT2, like NEVPT2, has not been discussed. In the present work, the equations of AC0 employing the Dyall Hamiltonian62 as the 0th-order Hamiltonian are presented and compared with NEVPT2 in detail. Like NEVPT2, the correlation energies computed by AC0 are decomposed into eight different subspaces. Detailed comparisons of correlation energies between AC0 and NEVPT2 are provided, and the importance of D corrections to AC0 for excited states is also discussed.

2 Theory

Throughout this work, the following conventions are used: {i, j, k, l, ⋯ } represent doubly occupied (occ), {t, u, v, w, ⋯ } denote active (act), {a, b, c, d, ⋯ } refer to virtual (vir), and {p, q, r, s} stand for unspecified orbitals. The spin-free excitation operator image file: d4fd00054d-t1.tif is employed for deriving all equations and matrix elements, avoiding the transformation from spin orbitals to spatial orbitals.

2.1 NEVPT2

In MRPT2, it is essential to define the full Hamiltonian (Ĥ), the 0th-order Hamiltonian (Ĥ0) and wave function (|Ψ0〉). The relationship between Ĥ and Ĥ0 is defined as:
 
Ĥ = [P with combining circumflex]Ĥ0[P with combining circumflex] + [Q with combining circumflex]Ĥ0[Q with combining circumflex] + [V with combining circumflex]. (1)
Here, [P with combining circumflex] and [Q with combining circumflex] are projection operators to the configuration state functions (CSFs) of the reference space and the 1st-order interaction space (FOIS), respectively. In the present work, the CASSCF or CASCI wave function serves as the 0th-order wave function. The [Q with combining circumflex] operator can be decomposed into eight different operators depending on the excitation type (Sk),
 
[Q with combining circumflex]Sijab = ÊbjÊai (2)
 
[Q with combining circumflex]Sija = ÊtjÊai (3)
 
[Q with combining circumflex]Siab = ÊbtÊai (4)
 
[Q with combining circumflex]Sij = ÊujÊti (5)
 
[Q with combining circumflex]Sab = ÊbuÊat (6)
 
[Q with combining circumflex]Si = ÊwuÊti (7)
 
[Q with combining circumflex]Sa = ÊwuÊat (8)
 
[Q with combining circumflex]Sia = ÊutÊai, ÊauÊti. (9)

In NEVPT2, the Dyall Hamiltonian is employed as the 0th-order Hamiltonian,62

 
image file: d4fd00054d-t2.tif(10)
 
image file: d4fd00054d-t3.tif(11)
where
 
êprqs = ÊpqÊrs − δqrÊps, (12)
 
image file: d4fd00054d-t4.tif(13)

Fii, defined as

 
image file: d4fd00054d-t5.tif(14)
and Faa defined analogously to Fii, represent the orbital energies of doubly occupied and virtual MOs, respectively. Notably, the utilization of the Dyall Hamiltonian ensures the absence of off-diagonal contributions in eqn (1),
 
image file: d4fd00054d-t6.tif(15)

Once the 1st-order wave function of NEVPT2 is obtained, the correlation energies could be computed as (from now on the acronym NEV will be used interchangeably with NEVPT2 for brevity)

 
image file: d4fd00054d-t7.tif(16)
where
 
|ΨSk〉 = [Q with combining circumflex]Sk|Ψ0〉. (17)

The amplitudes TSk are usually computed by solving the projected 1st-order equations within each subspace Sk,

 
〈[capital Psi, Greek, tilde]Sk|(Ĥ0 − E0)|ΨSk〉TSk + 〈[capital Psi, Greek, tilde]Sk|[V with combining circumflex]|Ψ0〉 = 0, (18)
where [capital Psi, Greek, tilde]Sk is the contravariant projector of ΨSk, 〈[capital Psi, Greek, tilde]Sk|ΨSk〉 = 1.89,90 The usage of [capital Psi, Greek, tilde]Sk could lead to simpler projected equations. As 〈ΨSk|Ĥ0|Ψ0〉 = 0, the 1st-order equations can be re-written as,
 
〈[capital Psi, Greek, tilde]Sk|(Ĥ0 − E0)|ΨSk〉TSk + 〈[capital Psi, Greek, tilde]Sk|Ĥ|Ψ0〉 = 0. (19)

Furthermore, with CASCI or CASSCF reference, the rank reduction trick can be used to reduce one order of RDMs,73

 
〈[capital Psi, Greek, tilde]Sk|[Ĥ0,[Q with combining circumflex]Sk]|Ψ0〉TSk + 〈[capital Psi, Greek, tilde]Sk|Ĥ|Ψ0〉 = 0. (20)

Consequently, RDMs up to 4th-order are needed for CAS reference functions. In particular, for the subspaces Sij, Sab, Sia, two active indices are present in the 1st-order wave function and eqn (20) involves up to 3rd-order RDMs in the active space. Equations for the subspaces Si and Sa are most demanding and they require up to 4th-order RDMs. The other subspaces Sk do not involve higher than 2nd-order RDMs in the active space. The NEVPT2 equations derived using eqn (20) have been reported elsewhere.73,90 Note that when using an approximate CASCI reference, the rank-reduction trick does not hold75 and the original equation, eqn (19), should be employed. Thus, in the Si and Sa subspaces, the 5th-order RDMs are required.

Since, as it has been shown, in NEVPT2 there is no coupling between different subspaces, the amplitude equations of different subspaces could be solved independently. Therefore, the total correlation energy of NEVPT2 can be computed individually within each subspace and subsequently summed up

 
image file: d4fd00054d-t8.tif(21)
where
 
image file: d4fd00054d-t9.tif(22)
Various eigenvalues (εμ) are obtained by solving the extended Koopmans’ equations,73 which have clear physical meanings.91 Specifically, εEA1μ(εEA2μ) and εIP1μ(εIP2μ) are single (double) electron affinities and single (double) ionization potentials in the active space, respectively. image file: d4fd00054d-t10.tif represent ionization potentials as well as shake-up state energies (electron affinities as well as shake-off state energies). Additionally, εEEμ denote electronic excitation energies within the active space. The explicit forms of the matrix elements OSk can be found in ref. 73 and 90.

2.2 AC0

The entire AC0 algorithm consists of two major steps: solving ERPA problems and computing correlation energies. While the explicit equations of AC0 with CAS reference have been published,82 it should be noted that they were formulated based on spin orbitals and were not decomposed into subspace contributions as in NEVPT2. In this section, the spinless AC0 equations based on the Dyall Hamiltonian are derived. The correlation energies of AC0 are separated into eight different subspaces as NEVPT2, which allows for a direct comparison of these two methods. It is also shown that the AC0 correlation energy expression is derivable from the second-order energy expression in the perturbation theory avoiding the adiabatic connection technique used in ref. 82.
2.2.1 Extended random phase approximation. In AC, the generalized eigenvalue equation of ERPA reads,
 
image file: d4fd00054d-t11.tif(23)
where Aα, Bα and M are defined as,
 
image file: d4fd00054d-t12.tif(24)
 
image file: d4fd00054d-t13.tif(25)

The Ĥα in eqn (24) is the linearly interpolated Hamiltonian between Ĥ0 and the exact Hamiltonian Ĥ,

 
Ĥα = Ĥ0 + αĤ′ = Ĥ0 + α(Ĥ − Ĥ0), (26)
corresponding, respectively, to α = 0 and α = 1 limits. ERPA has been derived from Rowe’s equation of motion.92 By assuming that an excitation operator is truncated to single excitations
 
image file: d4fd00054d-t14.tif(27)
the μth eigenstate of Ĥα could be generated from a ground state |Ψ0α〉,
 
|Ψμα〉 = Ôμα|Ψ0α〉. (28)

By assuming the consistency (killer) condition,

 
〈Ψ0α|Ôμα = 0, (29)
the matrices in the resulting equation of motion are expressible in terms of α-dependent RDMs up to the 2nd-order. In ERPA, see eqn (23), the RDMs from the reference wave function are employed. In equations of motion derived for the excitation operator as in eqn (27), the eigenvectors are directly related to one-electron transition density matrices
 
image file: d4fd00054d-t15.tif(30)
and in the ERPA model one finds
 
γμα = MXμα + MYμα. (31)

In the present work, the Dyall Hamiltonian, see eqn (10), is chosen as Ĥ0, which is different from the group Hamiltonian93 adopted in previous AC studies. Despite this difference in Ĥ0 selection, the explicit matrix elements engaged in AC0 remain exactly the same. The 1st- and 2nd-order RDMs of CASSCF involved in eqn (23) are symmetrized,

 
Γtu = 〈Ψ0|Êtu|Ψ0〉, ([capital Gamma, Greek, tilde]tu = 2δtu − Γtu), (32)
 
image file: d4fd00054d-t16.tif(33)

In eqn (32), the [capital Gamma, Greek, tilde]tu is the 1st-order hole density matrix.

Only two types of elements, A0 (B0) and A1 (B1), corresponding to α = 0 and α = 1 limits are evaluated in the AC0 method, respectively. At α = 0, the EPRA equation displays a block-diagonal structure depending on the single excitation types. This allows for its decoupling into four blocks: occupied–virtual (occ–vir), occupied–active (occ–act), active-virtual (act–vir), and active–active (act–act):

 
image file: d4fd00054d-t17.tif(34)
 
image file: d4fd00054d-t18.tif(35)
 
image file: d4fd00054d-t19.tif(36)
 
image file: d4fd00054d-t20.tif(37)
respectively. Various matrices in the above equations are defined as,
 
A0,occ–viriajb = δijδab(Faa − Fii), (38)
 
Mocc–viriajb = δijδab, (39)
 
A0,occ–actitju = δij(Kocc–acttu − [capital Gamma, Greek, tilde]tuFii), (40)
 
Mocc–actitju = δij(2δtu − Γtu) = δij([capital Gamma, Greek, tilde]tu), (41)
 
A0,act–virtaub = δab(Kact–virtu + ΓtuFaa), (42)
 
Mact–virtaub = δab(Γtu), (43)
 
image file: d4fd00054d-t21.tif(44)
 
B0,act–acttuwv = A0,act–acttuvw, (45)
 
image file: d4fd00054d-t22.tif(46)

Since there is no off-diagonal block in eqn (34)–(36), they could be further simplified to,

 
A0,occ–virXocc–virμ = wocc–virμMocc–virXocc–virμ, Yocc–virμ = 0, (47)
 
A0,occ–actXocc–actμ = wocc–actμMocc–actXocc–actμ, Yocc–actμ = 0, (48)
 
A0,act–virXact–virμ = wact–virμMact–virXact–virμ, Yact–virμ = 0. (49)
In eqn (47), Aocc–vir is diagonal, and Mocc–vir represents a unit matrix. Consequently, there is no need to solve the generalized eigenvalue (GE) equation in this block. For eqn (48) and (49), it suffices to solve the GE problems involving solely the active indexes,
 
Kocc–actCocc–actμ = εocc–actμ[capital Gamma, Greek, tilde]Cocc–actμ, (50)
 
Kact–virCact–virμ = εact–virμΓCact–virμ. (51)

Since the matrices

 
image file: d4fd00054d-t23.tif(52)
 
image file: d4fd00054d-t24.tif(53)
are the same as the Koopmans’ matrices of Sija and Siab subspace in NEVPT2,73 the eigenvalues and corresponding eigenvectors of the above two GE equations are exactly the same as the GE equations in NEVPT2 for the Sija and Siab subspaces,
 
εocc–actμ = εEA1μ, (54)
 
εact–virμ = εIP1μ. (55)

The final eigenvalues and eigenvectors of eqn (48) and (49) can be restored for each doubly occupied MO i or virtual MO a,

 
image file: d4fd00054d-t25.tif(56)
 
image file: d4fd00054d-t26.tif(57)

Since A0,act–act (B0,act–act) is exactly the same as A1,act–act (B1,act–act), i.e.

 
A1,act–acttuvw = B1,act–acttuwv = A0,act–acttuvw = B0,act–acttuwv (58)
the superscripts 0 and 1 for act–act matrices are omitted in the rest of this work. To address eqn (37), a transformation is required,
 
image file: d4fd00054d-t27.tif(59)
where
 
image file: d4fd00054d-t28.tif(60)
 
image file: d4fd00054d-t29.tif(61)
 
image file: d4fd00054d-t30.tif(62)
 
image file: d4fd00054d-t31.tif(63)

Then, the transformed eqn 59 can be solved as in conventional time-dependent HF,

 
(Āact–act + [B with combining macron]act–act)(Āact–act − [B with combining macron]act–act)([X with combining macron]act–actμ − Ȳact–actμ) = (wact–actμ)2([X with combining macron]act–actμ − Ȳact–actμ) (64)
and the vectors ([X with combining macron]act–actμ + Ȳact–actμ) are computed as
 
image file: d4fd00054d-t32.tif(65)

Assuming a properly-selected active space, all eigenvalues [wact–actμ]2 correspond to squared excitation energies in the active space and they should be positive for ground state calculations. However, when dealing with excited states, a few eigenvalues may turn negative or approach zero. In this work, if any [wact–actμ]2 is smaller than 10−6, the eigenvalue and its corresponding eigenvector are disregarded in the computations.

2.2.2 Energy expressions of AC0. The derivation of AC0 energy expression has been reported elsewhere.81,82 To compare it with NEVPT2, a new derivation based on the transition density matrices is given below. Without loss of generality, Ĥ′ = Ĥ − Ĥ0 is assumed to consist of one- and two-electron operators
 
image file: d4fd00054d-t33.tif(66)

The 2nd-order in α correction to the ground state energy,

 
E(2)0 = 〈Ψ(0)0|Ĥ′|Ψ(1)0〉, (67)
can be expressed in terms of 0th- and 1st-order one-electron transition density matrices,
 
[γ(0)μ]pq = 〈Ψ(0)0|Êpq|Ψ(0)μ〉, (68)
and
 
[γ(1)μ]pq = 〈Ψ(0)0|Êpq|Ψ(1)μ〉 + 〈Ψ(1)0|Êpq|Ψ(0)μ〉, (69)
respectively. In eqn (67), Ψ(1)0 is the 1st-order correction to the ground state wave function.

The formal derivation begins with noticing that from the completeness of eigenstates of Ĥα at any value of α

 
image file: d4fd00054d-t34.tif(70)
one obtains a resolution of identity
 
image file: d4fd00054d-t35.tif(71)
satisfied by 0th-order states and a null projector
 
image file: d4fd00054d-t36.tif(72)

By inserting eqn (71) into the two-body term of eqn (66), the expression for the 2nd-order correlation energy reads,

 
image file: d4fd00054d-t37.tif(73)

By adding image file: d4fd00054d-t38.tif, the final correlation energy could be computed from transition density matrices defined in eqn (68) and (69),

 
image file: d4fd00054d-t39.tif(74)

The above expression could be further re-written as,

 
image file: d4fd00054d-t40.tif(75)
where the 0th- and 1st-order RDMs of the ground state are defined as,
 
[γ(0)]pq = 〈Ψ(0)0|Êpq|Ψ(0)0〉, (76)
 
[γ(1)]pq = 〈Ψ(0)0|Êpq|Ψ(1)0〉 + 〈Ψ(1)0|Êpq|Ψ(0)0〉. (77)

Eqn (75) turns into the AC0 energy expression after assuming that the 1-RDM is fixed (γ(1) = 0) and the transition density matrices in the 0th- and 1st-order are obtained from the ERPA equation of motion, eqn (23),

 
image file: d4fd00054d-t41.tif(78)
The spin-free ERPA transition density matrices, eqn (31), are given by the eigenvectors [X(0)μ, Y(0)μ] corresponding to α = 0 and the 1st -order ones [X(1)μ, Y(1)μ] of ERPA, respectively. The latter follows from a standard perturbation theory treatment94 applied to GE in eqn (23), which, ultimately, yields the working formula for AC0 reading
 
image file: d4fd00054d-t42.tif(79)
where A(1) and B(1) are the 1st-order perturbing matrices simply obtained as A(1) = A1 − A0, similarly for B(1). The vectors Z+ and Z− in eqn (79) are 0th-order (α = 0) ERPA eigenvectors. Taking into account the decoupling of ERPA equations at α = 0, cf. eqn (34)–(37), the Z± vectors could be factorized as,
 
image file: d4fd00054d-t43.tif(80)
 
image file: d4fd00054d-t44.tif(81)

Depending on the μ and λ indices, the AC0 correlation energies can be classified into ten distinct components, see Table 1, which allows one to write the total AC0 energy as

 
image file: d4fd00054d-t45.tif(82)
where matrices OAC0μλ defined in eqn (79) are given in the Appendix and the eigenvalues ωμ take forms as shown in eqn (83), depending on the type of μ block,
 
image file: d4fd00054d-t46.tif(83)

Table 1 The decomposition of AC0 correlation energies into different subspaces depending on the μ and λ in eqn (79) and their correspondence with the subspaces of NEVPT2
image file: d4fd00054d-u1.tif

When both μ and λ are in the act–act block, the final correlation energy is always zero. This is due to the fact that,

 
image file: d4fd00054d-t47.tif(84)

Among the remaining nine components (μ, λ), two of them, (occ–act, act–vir) and (occ–vir, act–act), can be regarded as belonging to Sia subspace. Consequently, akin to NEVPT2, the correlation energy of AC0 can be computed across eight different subspaces, presented in Table 1. The important result of this section is that by showing that the correlation energy in the 2nd-order as given in eqn (67), is equivalent to the formula in eqn (75) expressed by the 1st-order TRDMs. This justifies a comparison of ESk terms contributing to AC0 with their counterpart contributing to NEVPT2 energy. The explicit energy expressions of ESkAC0 are given subspace by subspace in the Appendix. Remarkably, the energy expressions of the Sijab, Sija and Siab subspaces are exactly the same as those of FIC-NEVPT2 (ref. 73)

 
image file: d4fd00054d-t48.tif(85)
 
image file: d4fd00054d-t49.tif(86)
 
image file: d4fd00054d-t50.tif(87)

This means that for subspaces Sk involving only one active index, contributions to NEVPT2 and AC0 energies are identical.

2.2.3 D correction. It has been found that the AC0 method can produce ground state energies as accurate as NEVPT2. However, it often overestimates the excitation energies, even for low-lying excited states. Recently, one of us proposed the D correction to recover the missing correlation energies in AC0 computed for an excited state by considering the negative-transition density matrices.86 In this subsection, the explicit expression of the D correction is recapitulated, and the final D correction is separated into three Sk subspaces.

For a given excited state J, computation of the AC0 correlation energy proceeds according to the same protocol as described in Sections 2.2.1 and 2.2.2, the difference being that the ERPA GEs, eqn (34)–(36), are solved for ΨJ CAS reference wave function, instead of that of ground state Ψ0. It is known that all of the eigenvalues wact–act, corresponding to excitations within the active orbitals space, are positive in the quasiparticle ERPA, see eqn (37). Thus, only the positive-energy transitions are considered in the calculation of excited states. However, by enlarging the manifold of the excitation operators in ERPA, there should be, as well, negative eigenvalues wact–act corresponding to transitions from state J to lower states, at least to the ground state. To account for the missing correlation energies due to the de-excitation in the AC0 calculation for state J, we could evaluate the correction from all states I(I < J) with lower excitation energies instead. In the AC0 calculations of state I, there should be a one-to-one correspondence between wact–act and ECASK − ECASI, where K is a singly excited state higher than I. In principle, each state I with lower energy than state J should contribute D corrections to the correlation energies of state J. For state J, the D correction from state I is computed as follows:

 
image file: d4fd00054d-t51.tif(88)
(see the Appendix for explicit expression of the matrices O′AC0 and OAC0). Only the Sia, Si, and Sa subspaces, involving excitations between states I and J in the active space, wact–actλ of state I, contribute to the D correction of state J. The final correlation energy of AC0D for state J is computed as follows:
 
image file: d4fd00054d-t52.tif(89)

Although there is a one-to-one correspondence between wact–act of state I and CASSCF energy differences of states J and I, ECASJ − ECASI, the assignment is not straightforward, especially when there are many states with energies close to state J. To assist in the assignment, for molecules with spatial symmetries, the irreducible representation (irrep) of each wact–actλ is determined from the irreps of orbital pairs with non-zero elements in Yact–actλ.

3 Computational details

All calculations were performed using a development version of ORCA.95–97 For the potential energy surface (PES) calculations of N2, to ensure that the Sij subspace has more significant correlation energies, all orbitals were correlated. Therefore, the cc-pwcVQZ basis set was employed.98 For the calculations of Cr2 and selected organic molecules from Thiel’s test set,99 the cc-pVQZ basis set was used with the frozen core approximation.100 The resolution of identity (RI) approximation was utilized to reduce the computational costs involving two-electron integrals in all correlation calculations, with auxiliary basis sets appropriate for the cc-pwcVQZ or cc-pVQZ basis sets, respectively. For Cr2 calculations, only the 1s2s2p orbitals are frozen.

For the calculations of excited states of organic molecules, we adopted the active spaces reported by Thiel and coworkers, which are detailed in Table 2. Specifically, for furan, pyrrole, benzene, and octatetraene, only the π electrons and orbitals are included in the active space. In the case of imidazole, in addition to the π electrons and orbitals, the lone-pair orbital of nitrogen is also included. We also adopted the computational procedure suggested by Thiel and coworkers.99 The ground state energy using the state-specific (SS) CASSCF reference was computed to calculate the excitation energies of states with different symmetries compared to the ground state. For the D correction, it is necessary to determine the energy order of excited states at the CASSCF level beforehand. In some cases, the order of excited states obtained from state-averaged (SA) CASSCF calculations may differ from that obtained from SS-CASSCF calculations. The different orders of excited states could lead to slightly different D corrections. In this work, the order from SA-CASSCF calculations for states with the same spin symmetry was used, which is consistent with a previous work published by one of the authors.86 The D corrections were computed using the same SA-CASSCF reference as well.

Table 2 The correlation energies (in eV) of subspaces Sk and the total correlation energies of AC0, AC0D and NEVPT2 (NEV) for selected organic molecules. For the AC0(D) results, only the deviations with respect to NEVPT2 are given. The correlation energies of Sijab, Siab and Sija subspaces are exactly the same for AC0(D) and NEVPT2, therefore they are not shown but they are included in the total correlation energies
Molecule (active space) State Sab Sij Sia Sa Si Total correlation energy
AC0 NEV AC0 NEV AC0 AC0D NEV AC0 AC0D NEV AC0 AC0D NEV AC0 AC0D NEV
Furan CAS(6,5) 11A1(SS) −0.14 −1.05 −0.01 −0.05 0.21 0.21 −2.99 0.04 0.04 −0.34 0.00 0.00 0.00 0.09 0.09 −23.87
11A1(SA) −0.15 −1.12 −0.01 −0.05 0.23 0.23 −2.77 0.07 0.07 −0.51 0.00 0.00 0.00 0.13 0.13 −24.27
21A1 −0.11 −0.99 −0.01 −0.04 0.57 0.39 −2.52 0.14 0.11 −0.50 0.00 0.00 0.00 0.59 0.39 −24.13
31A1 −0.12 −1.16 −0.01 −0.07 2.24 1.01 −3.96 0.57 0.28 −0.85 0.00 0.00 0.00 2.68 1.16 −26.16
11B2 −0.17 −1.30 −0.01 −0.07 1.23 0.54 −3.64 0.34 0.11 −0.92 0.00 0.00 0.00 1.39 0.47 −25.95
13A1 −0.10 −0.95 0.00 −0.04 0.26 0.07 −2.30 0.03 −0.02 −0.38 0.00 0.00 0.00 0.18 −0.06 −23.66
13B2 −0.12 −0.97 −0.01 −0.04 0.18 0.18 −2.26 0.01 0.01 −0.27 0.00 0.00 0.00 0.06 0.06 −23.43
Pyrrole CAS(6,5) 11A1(SS) −0.15 −1.06 −0.01 −0.05 0.21 0.21 −3.00 0.07 0.07 −0.42 0.00 0.00 0.00 0.12 0.12 −23.30
11A1(SA) −0.16 −1.12 −0.01 −0.05 0.21 0.21 −2.73 0.06 0.06 −0.50 0.00 0.00 0.00 0.10 0.10 −23.53
21A1 −0.11 −1.02 0.00 −0.04 0.54 0.35 −2.43 0.17 0.13 −0.52 0.00 0.00 0.00 0.59 0.37 −23.44
31A1 −0.12 −1.14 −0.01 −0.06 1.80 0.82 −3.46 0.54 0.30 −0.80 0.00 0.00 0.00 2.21 0.99 −24.96
11B2 −0.16 −1.25 −0.01 −0.06 0.86 0.14 −3.04 0.25 0.02 −0.80 0.00 0.00 0.00 0.93 −0.01 −24.64
13A1 −0.11 −0.99 0.00 −0.03 0.29 0.05 −2.32 0.05 −0.02 −0.47 0.00 0.00 0.00 0.23 −0.09 −23.19
13B2 −0.13 −0.99 −0.01 −0.04 0.17 0.17 −2.21 0.01 0.01 −0.29 0.00 0.00 0.00 0.05 0.05 −22.78
Imidazole CAS(8,7) 11A′(SS) −0.22 −1.74 −0.01 −0.04 0.20 0.20 −3.71 0.08 0.08 −0.76 0.00 0.00 −0.01 0.04 0.04 −23.69
11A′(SA) −0.25 −1.89 −0.01 −0.05 0.23 0.23 −3.42 0.10 0.10 −0.89 0.00 0.00 −0.01 0.07 0.07 −24.12
21A′ −0.21 −1.81 −0.01 −0.04 0.60 0.21 −3.25 0.18 0.07 −0.91 0.00 0.00 −0.01 0.57 0.07 −24.12
31A′ −0.26 −2.04 −0.01 −0.06 1.12 0.65 −3.90 0.41 0.23 −1.16 0.00 −0.01 −0.01 1.25 0.59 −25.31
41A′ −0.22 −1.90 −0.01 −0.05 1.49 0.51 −3.95 0.48 0.07 −1.13 0.00 −0.01 −0.01 1.74 0.35 −25.22
11A′′ −0.23 −1.84 0.00 −0.03 0.30 0.27 −3.29 0.06 −0.06 −0.75 0.03 0.02 −0.09 0.15 0.00 −23.82
21A′′ −0.23 −1.86 −0.01 −0.03 0.35 0.35 −3.41 0.06 0.04 −0.83 0.04 0.03 −0.11 0.21 0.18 −24.16
13A′ −0.20 −1.71 −0.01 −0.03 0.32 0.32 −3.19 0.08 0.08 −0.77 0.00 0.00 −0.01 0.20 0.20 −23.65
23A′ −0.20 −1.75 −0.01 −0.04 0.36 0.33 −3.13 0.12 0.11 −0.82 0.00 0.00 −0.01 0.28 0.23 −23.76
33A′ −0.23 −1.87 −0.01 −0.04 0.50 0.29 −3.41 0.21 0.15 −1.06 0.00 0.00 −0.01 0.47 0.20 −24.46
43A′ −0.21 −1.77 0.00 −0.03 0.77 0.15 −3.62 0.35 0.01 −1.11 0.00 −0.01 −0.01 0.91 −0.07 −24.67
13A′′ −0.22 −1.81 −0.01 −0.03 0.31 0.29 −3.30 0.09 −0.01 −0.80 0.02 0.01 −0.08 0.18 0.06 −23.81
23A′′ −0.22 −1.78 0.00 −0.03 0.36 0.34 −3.45 0.10 0.02 −0.89 0.03 0.03 −0.11 0.26 0.17 −24.15
Benzene CAS(6,6) 11Ag −0.10 −0.73 −0.02 −0.09 0.25 0.25 −3.64 0.04 0.04 −0.28 0.00 0.00 0.00 0.17 0.17 −25.49
11E2g −0.06 −0.57 −0.01 −0.07 0.41 −0.59 −2.79 0.06 −0.09 −0.29 0.00 0.00 0.00 0.41 −0.75 −25.15
11B2u −0.07 −0.64 −0.01 −0.09 0.34 0.29 −2.99 0.01 0.00 −0.27 0.00 0.00 0.00 0.26 0.21 −25.25
11B1u −0.14 −1.07 −0.02 −0.12 1.22 0.65 −3.92 0.24 0.16 −0.71 0.00 0.00 0.00 1.30 0.65 −27.29
11E1u −0.11 −1.00 −0.02 −0.12 1.69 0.46 −4.42 0.35 0.14 −0.83 0.00 0.00 0.00 1.92 0.47 −27.90
13E2g −0.05 −0.55 −0.01 −0.07 0.34 −0.39 −2.75 0.05 −0.09 −0.28 0.00 0.00 0.00 0.33 −0.54 −25.06
13E1u −0.07 −0.70 −0.01 −0.08 0.38 0.03 −3.11 0.11 0.03 −0.47 0.00 0.00 0.00 0.40 −0.02 −25.67
13B2u −0.10 −0.95 −0.02 −0.11 0.89 0.11 −4.23 0.21 0.06 −0.88 0.00 0.00 0.00 0.98 0.05 −27.50
13B1u −0.06 −0.59 −0.01 −0.08 0.26 0.26 −2.93 0.02 0.02 −0.23 0.00 0.00 0.00 0.21 0.21 −25.02
Octatetraene CAS(8,8) 11Ag(SS) −0.13 −0.87 −0.03 −0.11 0.31 0.31 −4.41 0.01 0.01 −0.22 0.00 0.00 0.00 0.16 0.16 −34.34
11Ag(SA) −0.14 −0.91 −0.02 −0.12 0.30 0.30 −4.18 0.00 0.00 −0.38 0.00 0.00 0.00 0.15 0.15 −34.67
21Ag −0.09 −0.81 −0.02 −0.11 0.41 0.07 −4.07 0.03 −0.01 −0.40 0.00 0.00 0.00 0.34 −0.05 −34.64
31Ag −0.13 −1.20 −0.02 −0.16 1.98 −0.07 −5.92 0.66 0.37 −1.21 0.00 0.00 0.00 2.49 0.15 −37.87
41Ag −0.08 −0.77 −0.01 −0.10 0.47 −0.04 −3.85 0.17 0.08 −0.49 0.00 0.00 0.00 0.55 −0.05 −34.62
11Bu −0.14 −1.21 −0.02 −0.15 1.32 0.61 −5.66 0.41 0.31 −1.09 0.00 0.00 0.00 1.57 0.75 −37.44
21Bu −0.08 −0.79 −0.01 −0.10 0.44 0.06 −3.89 0.08 0.02 −0.46 0.00 0.00 0.00 0.43 −0.01 −34.64
31Bu −0.06 −0.73 −0.01 −0.10 0.54 −0.98 −3.62 0.16 −0.08 −0.49 0.00 0.00 0.00 0.62 −1.14 −34.49
13Ag −0.09 −0.79 −0.02 −0.10 0.35 0.29 −4.05 0.04 0.02 −0.34 0.00 0.00 0.00 0.28 0.21 −34.40
13Bu −0.10 −0.82 −0.02 −0.11 0.33 0.33 −4.28 0.03 0.03 −0.32 0.00 0.00 0.00 0.24 0.24 −34.48
ME −0.14 — −0.01 — 0.58 0.22 — 0.15 0.06 — 0.00 0.00 — 0.58 0.12 —
RMSE 0.15 — 0.01 — 0.75 0.41 — 0.21 0.10 — 0.01 0.01 — 0.84 0.42 —

4 Results

The accuracy of AC0 in describing bond breaking processes and singlet–triplet gaps has been thoroughly studied by Pernal and coworkers.82,85 AC0 has been shown to produce results as accurate as NEVPT2 in most systems at the ground state. To provide a more detailed comparison between AC0 and NEVPT2, the energy differences within subspaces Sk for AC0 are presented. For the Sijab, Siab, and Sija subspaces, we have demonstrated that their energy expressions in NEVPT2 and AC0 are exactly the same, which has been further confirmed numerically. Therefore, for the comparison between AC0 and NEVPT2, only the remaining five subspaces are considered.

The energy deviation of AC0 with respect to NEVPT2 for the N2 molecule is depicted in Fig. 1. Energies from both methods have been obtained with CAS(10,8) wave function. Along the entire PES, the correlation energies of the Sab subspace are significantly overestimated by AC0, while those of the Sa subspace are underestimated. The errors from the Sab and Sa subspaces slightly cancel each other out. Since there are only two doubly occupied MOs in the calculations, the correlation energies of the Sij, Si, and Sia subspaces are not very significant. The overestimated energies of Sab in AC0 can be attributed to the violation of Pauli’s exclusion principle, as discussed in the literature.101 In the AC0 wave function of the Sab subspace, one active electron is allowed to be excited to the same virtual molecular orbital (MO) twice. Throughout the PES of N2, AC0 overestimates the correlation energy by an average of 0.25 eV. Although the binding energies predicted by AC0 are close to those by NEVPT2, there is a hump at the bond dissociation region (1.6 Å). This could be attributed to the missing high-order particle-hole contribution in the Sa subspace of AC0, which is important for the strongly correlated region.


image file: d4fd00054d-f1.tif
Fig. 1 The deviations (in mEh) of AC0 correlation energies with respect to NEVPT2 for N2.

The AC0 correlation energies of Cr2 with respect to those of NEVPT2, utilizing CAS(12,12) as reference, are shown in Fig. 2. Similar to N2, the AC0 correlation energies of the Sab subspace are overestimated by about 25 mEh throughout the curve. Moreover, since there are more doubly occupied MOs, the correlation energies from the Sij subspace are also overestimated. In contrast, the AC0 energies of the Sia, Si, and Sa subspaces, which involve wact–act, are underestimated with respect to those of NEVPT2. There is a hump located at the equilibrium region for Sa. It was reported by Roos and coworkers that CASPT2 suffers from the intruder state problem around the equilibrium bond length of Cr2.102 For Cr2, the equilibrium bond distance region might be more strongly correlated. The difference of Sa between AC0 and NEVPT2 could be an indicator of strong correlations and the need to go beyond single excitations in the active space (or even double excitations, as in this region also NEVPT2 struggles) in the excitation operator in ERPA. The PES computed by AC0 and NEVPT2 is shown in Fig. 3. Due to the hump introduced by the Sa subspace, the binding energy of Cr2 is underestimated by AC0 compared to NEVPT2.


image file: d4fd00054d-f2.tif
Fig. 2 The deviations (in mEh) of AC0 correlation energies with respect to NEVPT2 for Cr2.

image file: d4fd00054d-f3.tif
Fig. 3 The potential energy curves of Cr2 computed by AC0 and NEVPT2. The newly fitted experimental results are given for comparison.103

The excited states of five organic molecules selected from Thiel’s test set99 are computed using AC0(D) and NEVPT2. The contributions to correlation energies from Sk subspaces and the total correlation energies of NEVPT2 together with the corresponding deviations between AC0(D) and NEVPT2 are presented in Table 2. For almost all states, similar to the calculations of diatomic molecules, the AC0 correlation energies of Sab and Sij subspaces are overestimated, whereas the energies of the remaining three subspaces are underestimated, using NEVPT2 results as the reference. For Sab and Sij subspaces, the errors of AC0 between ground states and excited states are insignificant. AC0 and NEVPT2 correlation energies of Si subspaces are practically the same. However, for the other two subspaces involving wact–act, Sia and Sa, the AC0 approximation significantly underestimates the correlation energies of excited states. The higher the excitation energies, the greater the deviations. These observations confirm the important role of the negative transitions in the active space, which ERPA lacks. The errors between NEVPT2 and AC0 are mainly due to the Sia subspace. For the Sia subspace, the mean error (ME) and root mean squared error (RMSE) of AC0 with respect to NEVPT2 are 0.58 eV and 0.75 eV, respectively. These errors are close to the overall errors between NEVPT2 and AC0. By adding the D corrections to the subspace, the ME of Sia with respect to the NEVPT2 results is reduced to 0.22 eV. For the Sa space the ME error is reduced to 0.06 eV. Significant reductions of the AC0 correlation energy errors in Sia and Sa spaces after adding the D correction, proves that the lowering of the error is not coincidental. It is a consequence of effectively adding contributions to the correlation energies of excited states from negative transitions in the active space by the correction.

The results of excitation energies are summarized in Table 3. The AC0 approximation tends to overestimate the excitation energies compared to NEVPT2, which is a consequence of underestimating the correlation energies of Sia and Sa subspaces, as already discussed. By including the D correction the errors in these subspaces are reduced and the results are significantly improved for most low-lying excited states. The RMSE of excitation energies computed by the AC0D is reduced to 0.42 eV, which is close to that of NEVPT2, employing CC3 results as reference. For excited states with excitation energies higher than 6.0 eV, the AC0D corrections may not always improve the results systematically. For example, AC0 delivers more accurate excitation energies for the 13B2u state of benzene than that of NEVPT2. By including the D correction, the excitation energy is reduced by about 0.93 eV. To systematically improve the accuracy of AC0 for excited states, one should construct more accurate transition densities from high-order ERPA, using projection spaces beyond the one-particle-one-hole space. Note that for the 11Bu state of octatetraene, the NEVPT2 method underestimates the excitation energies significantly, by 1.19 eV, due to the mixing of Rydberg and valence states. Previous studies using triple–ζ basis sets also showed underestimations of about 0.9 eV.104,105 To achieve more accurate results for the 11Bu state, it may be necessary to increase the size of the active space.106

Table 3 The excitation energies (in eV) of selected organic molecules computed by AC0, AC0D and NEVPT2 with respect to CC3 results.99
Molecule State AC0 AC0D NEVPT2 CC3
Furan 21A1 0.55 0.35 0.10 6.62
31A1 2.20 0.46 0.68 8.53
11B2 0.90 −0.02 −0.40 6.60
13A1 0.23 −0.02 0.13 5.48
13B2 0.12 0.12 0.15 4.17
Pyrrole 21A1 0.62 0.40 0.14 6.40
31A1 1.99 0.77 −0.12 8.17
11B2 0.71 −0.24 −0.10 6.71
13A1 0.25 −0.07 0.14 5.51
13B2 0.18 0.18 0.25 4.48
Imidazole 21A′ 0.70 0.20 0.21 6.58
31A′ 0.81 0.15 −0.37 7.10
41A′ 1.68 0.29 0.02 8.45
11A′′ 0.12 −0.03 0.02 6.82
21A′′ 0.08 0.05 −0.08 7.93
13A′ 0.23 0.23 0.07 4.69
23A′ 0.29 0.25 0.06 5.79
33A′ 0.35 0.07 −0.08 6.55
43A′ 0.47 −0.51 −0.39 7.42
13A′′ 0.12 0.01 −0.01 6.37
23A′′ 0.14 0.05 −0.08 7.51
Benzene 11E2g 0.20 −0.96 −0.03 8.43
11B2u 0.24 0.18 0.15 5.07
11B1u 0.57 −0.08 −0.56 6.68
11E1u 1.14 −0.31 −0.61 7.45
13E2g 0.27 −0.60 0.11 7.49
13E1u 0.22 −0.20 −0.01 4.90
13B2u −0.05 −0.98 −0.86 6.04
13B1u 0.25 0.25 0.22 4.12
Octatetraene 21Ag −0.05 −0.44 −0.25 4.97
31Ag 1.58 −0.76 −0.76 6.50
41Ag 0.28 −0.32 −0.12 6.81
11Bu 0.21 −0.60 −1.19 4.94
21Bu 0.04 −0.41 −0.23 6.06
31Bu 0.89 −0.87 0.43 7.91
13Ag 0.17 0.10 0.05 3.67
13Bu 0.09 0.09 0.01 2.30
ME 0.51 −0.08 −0.12 —
MUE 0.51 0.32 0.24 —
RMSE 0.75 0.42 0.35 —

5 Conclusions

In this study, we have revisited the working equations of AC0 using spin-free operators and introduced a factorization of the correlation energies into eight distinct subspaces, akin to FIC-NEVPT2. This has facilitated a systematic comparison of the AC0 and NEVPT2 methods and identified subspaces in which AC0 and NEVPT2 deviate most significantly. We have demonstrated that, despite the fact that AC0 is based on equations of motion with the excitation operator constrained to single excitations, while FOIS in NEVPT2 is generated by doubly-exciting operators, the equations governing the Sijab, Siab, and Sija subspaces in AC0 are identical to those in FIC-NEVPT2. However, we observed that the correlation energies of the Sij and Sab subspaces in AC0 tend to be overestimated due to violations of Pauli’s exclusion principle. For strongly correlated systems, the Sia, Si, and Sa subspaces play a crucial role in producing reliable results in any MRPT2 method. In AC0, as the ERPA equation is solved within the one-particle-one-hole space, it may not yield results as accurate as NEVPT2. This problem is particularly acute in exited state calculations. We have shown that the overestimation of excited state energies by AC0 is due to positive errors in Sia and Sa subspaces, calculated with respect to their NEVPT2 counterparts. To address this limitation, the AC0D method approximates negative excitation corrections and has been shown to produce improved results compared to AC0. To systematically enhance the accuracy of the AC0 method, further incorporating larger excitation operator manifolds in the ERPA, such as the particle–particle ERPA, is essential. We are currently engaged in projects in this direction.

Conflicts of interest

There are no conflicts to declare.

Appendix

Appendix: AC0 energy expressions of eight subspaces

Notice that in the following expression, the upper index AC0 has been removed from OSk matrices for brevity.

Sijab:

 
Oijab = (ia|jb)[2(ia|jb) − (ib|ja)] (A.1)
 
image file: d4fd00054d-t53.tif(A.2)
Sija:
 
image file: d4fd00054d-t54.tif(A.3)
 
Oija = (ia|jμ)[2(ia|jμ) − (iμ|ja)] (A.4)
 
image file: d4fd00054d-t55.tif(A.5)
Siab:
 
image file: d4fd00054d-t56.tif(A.6)
 
Oiab = (ia|μb)[2(ia|μb) − (ib|μa)] (A.7)
 
image file: d4fd00054d-t78.tif(A.8)
Sij:
 
image file: d4fd00054d-t57.tif(A.9)
 
image file: d4fd00054d-t58.tif(A.10)
 
Oij = (iμ|jλ)Niμ,jλ (A.11)
 
image file: d4fd00054d-t59.tif(A.12)
Sab:
 
image file: d4fd00054d-t60.tif(A.13)
 
image file: d4fd00054d-t61.tif(A.14)
 
Oab = (μa|λb)Nμa,λb (A.15)
 
image file: d4fd00054d-t62.tif(A.16)
Si:
 
image file: d4fd00054d-t63.tif(A.17)
 
image file: d4fd00054d-t64.tif(A.18)
 
Oi = (iμ|λ)Niμ,λ (A.19)
 
image file: d4fd00054d-t65.tif(A.20)
Sa:
 
image file: d4fd00054d-t66.tif(A.21)
 
image file: d4fd00054d-t67.tif(A.22)
 
Oa = (μa|λ)Nμa,λ (A.23)
 
image file: d4fd00054d-t68.tif(A.24)
Sia:
 
image file: d4fd00054d-t69.tif(A.25)
 
image file: d4fd00054d-t70.tif(A.26)
 
image file: d4fd00054d-t71.tif(A.27)
 
image file: d4fd00054d-t72.tif(A.28)
 
Oia = (iμ|λa)Niμ,λa (A.29)
 
image file: d4fd00054d-t73.tif(A.30)
 
image file: d4fd00054d-t74.tif(A.31)
 
image file: d4fd00054d-t75.tif(A.32)

In the above expressions, the A1 and B1 matrices are closely related to the orbital Hessian of the reference wave function,

 
image file: d4fd00054d-t76.tif(A.33)
where
 
Gpq = ΓprFcqr + Γpr,sp′(qr|sp′), (A.34)
 
image file: d4fd00054d-t77.tif(A.35)

Acknowledgements

This work received support from the National Natural Science Foundation of China (grant no. 22273052), the Qilu Young Scholar Program of Shandong University, and the National Science Center of Poland under grant no. 2019/35/B/ST4/01310. YG expresses gratitude to KP for the encouragement. Additionally, YG acknowledges stimulating discussions with Kantharuban Sivalingam, Frank Neese, and Wenjian Liu.

References

  1. G. Das and A. C. Wahl, J. Chem. Phys., 1966, 44, 87–96 CrossRef CAS .
  2. E. Clementi and A. Veillard, J. Chem. Phys., 1966, 44, 3050–3053 CrossRef CAS .
  3. J. Hinze and C. C. Roothaan, Prog. Theor. Phys., Suppl., 1967, 40, 37–51 CrossRef .
  4. B. Levy, Int. J. Quantum Chem., 1970, 4, 297–313 CrossRef CAS .
  5. K. Ruedenberg, L. Cheung and S. Elbert, Int. J. Quantum Chem., 1979, 16, 1069–1101 CrossRef CAS .
  6. B. O. Roos, P. R. Taylor and P. E. M. Siegbahn, Chem. Phys., 1980, 48, 157–173 CrossRef CAS .
  7. B. O. Roos, Adv. Chem. Phys., 1987, 69, 399–445 CrossRef CAS .
  8. K. D. Vogiatzis, D. Ma, J. Olsen, L. Gagliardi and W. A. de Jong, J. Chem. Phys., 2017, 147, 184111 CrossRef PubMed .
  9. H. Gao, S. Imamura, A. Kasagi and E. Yoshida, J. Chem. Theory Comput., 2024, 20, 1185–1192 CrossRef CAS PubMed .
  10. L. Cheung, K. Sundberg and K. Ruedenberg, J. Am. Chem. Soc., 1978, 100, 8024–8025 CrossRef CAS .
  11. K. Ruedenberg, M. W. Schmidt, M. M. Gilbert and S. Elbert, Chem. Phys., 1982, 71, 41–49 CrossRef CAS .
  12. B. Huron, J. P. Malrieu and P. Rancurel, J. Chem. Phys., 1973, 58, 5745–5759 CrossRef CAS .
  13. R. J. Buenker and S. D. Peyerimhoff, Theor. Chim. Acta, 1968, 12, 183–199 CrossRef CAS .
  14. R. J. Buenker and S. D. Peyerimhoff, Theor. Chim. Acta, 1974, 35, 33–58 CrossRef CAS .
  15. W. A. Goddard III, T. H. Dunning Jr, W. J. Hunt and P. J. Hay, Acc. Chem. Res., 1973, 6, 368–376 CrossRef .
  16. S. P. Walch, C. W. Bauschlicher, B. O. Roos and C. J. Nelin, Chem. Phys. Lett., 1983, 103, 175–179 CrossRef CAS .
  17. J. Olsen, B. O. Roos, P. Jørgensen and H. J. A. Jensen, J. Chem. Phys., 1988, 89, 2185–2192 CrossRef CAS .
  18. T. Fleig, J. Olsen and C. M. Marian, J. Chem. Phys., 2001, 114, 4775–4790 CrossRef CAS .
  19. T. Fleig, J. Olsen and L. Visscher, J. Chem. Phys., 2003, 119, 2963–2971 CrossRef CAS .
  20. D. Ma, G. Li Manni and L. Gagliardi, J. Chem. Phys., 2011, 135, 044128 CrossRef PubMed .
  21. J. Ivanic, J. Chem. Phys., 2003, 119, 9364–9376 CrossRef CAS .
  22. N. M. Tubman, C. D. Freeman, D. S. Levine, D. Hait, M. Head-Gordon and K. B. Whaley, J. Chem. Theory Comput., 2020, 16, 2139–2159 CrossRef CAS PubMed .
  23. D. S. Levine, D. Hait, N. M. Tubman, S. Lehtola, K. B. Whaley and M. Head-Gordon, J. Chem. Theory Comput., 2020, 16, 2340–2354 CrossRef CAS PubMed .
  24. F. A. Evangelista, J. Chem. Phys., 2014, 140, 124114 CrossRef PubMed .
  25. J. B. Schriber and F. A. Evangelista, J. Chem. Phys., 2016, 144, 161106 CrossRef PubMed .
  26. J. B. Schriber and F. A. Evangelista, J. Chem. Theory Comput., 2017, 13, 5354–5366 CrossRef CAS PubMed .
  27. V. G. Chilkuri and F. Neese, J. Comput. Chem., 2021, 42, 982–1005 CrossRef CAS PubMed .
  28. V. G. Chilkuri and F. Neese, J. Chem. Theory Comput., 2021, 17, 2868–2885 CrossRef CAS PubMed .
  29. W. Liu and M. R. Hoffmann, J. Chem. Theory Comput., 2016, 12, 1169–1178 CrossRef CAS PubMed .
  30. N. Zhang, W. Liu and M. R. Hoffmann, J. Chem. Theory Comput., 2020, 16, 2296–2316 CrossRef CAS PubMed .
  31. N. Zhang, W. Liu and M. R. Hoffmann, J. Chem. Theory Comput., 2021, 17, 949–964 CrossRef CAS PubMed .
  32. Y. Guo, N. Zhang, Y. Lei and W. Liu, J. Chem. Theory Comput., 2021, 17, 7545–7561 CrossRef CAS PubMed .
  33. N. Zhang, Y. Xiao and W. Liu, J. Phys.: Condens. Matter, 2022, 34, 224007 CrossRef CAS PubMed .
  34. Y. Guo, N. Zhang and W. Liu, J. Chem. Theory Comput., 2023, 19, 6668–6685 CrossRef CAS PubMed .
  35. J. E. T. Smith, B. Mussard, A. A. Holmes and S. Sharma, J. Chem. Theory Comput., 2017, 13, 5468–5478 CrossRef CAS PubMed .
  36. Y. Yao and C. J. Umrigar, J. Chem. Theory Comput., 2021, 17, 4183–4194 CrossRef CAS PubMed .
  37. P. M. Zimmerman and A. E. Rask, J. Chem. Phys., 2019, 150, 244117 CrossRef PubMed .
  38. J. W. Park, J. Chem. Theory Comput., 2021, 17, 1522–1534 CrossRef CAS PubMed .
  39. Y. Garniron, A. Scemama, E. Giner, M. Caffarel and P.-F. Loos, J. Chem. Phys., 2018, 149, 064103 CrossRef PubMed .
  40. G. H. Booth, A. J. W. Thom and A. Alavi, J. Chem. Phys., 2009, 131, 054106 CrossRef PubMed .
  41. R. E. Thomas, Q. Sun, A. Alavi and G. H. Booth, J. Chem. Theory Comput., 2015, 11, 5316–5325 CrossRef CAS PubMed .
  42. W. Dobrautz, S. D. Smart and A. Alavi, J. Chem. Phys., 2019, 151, 094104 CrossRef PubMed .
  43. E. Giner, A. Scemama and M. Caffarel, Can. J. Chem., 2013, 91, 879–885 CrossRef CAS .
  44. E. Giner, R. Assaraf and J. Toulouse, Mol. Phys., 2016, 114, 910–920 CrossRef CAS .
  45. D. Zgid and M. Nooijen, J. Chem. Phys., 2008, 128, 144116 CrossRef PubMed .
  46. D. Ghosh, J. Hachmann, T. Yanai and G. K.-L. Chan, J. Chem. Phys., 2008, 128, 144117 CrossRef PubMed .
  47. T. Yanai, Y. Kurashige, D. Ghosh and G. K.-L. Chan, Int. J. Quantum Chem., 2009, 109, 2178–2190 CrossRef CAS .
  48. Y. Ma and H. Ma, J. Chem. Phys., 2013, 138, 224105 CrossRef PubMed .
  49. S. Wouters, T. Bogaerts, P. Van Der Voort, V. Van Speybroeck and D. Van Neck, J. Chem. Phys., 2014, 140, 241103 CrossRef PubMed .
  50. Y. Ma, S. Knecht, S. Keller and M. Reiher, J. Chem. Theory Comput., 2017, 13, 2533–2549 CrossRef CAS PubMed .
  51. Q. Sun, J. Yang and G. K.-L. Chan, Chem. Phys. Lett., 2017, 683, 291–299 CrossRef CAS .
  52. Y. Kurashige and T. Yanai, J. Chem. Phys., 2011, 135, 094104 CrossRef PubMed .
  53. Y. Kurashige, J. Chalupský, T. N. Lan and T. Yanai, J. Chem. Phys., 2014, 141, 174111 CrossRef PubMed .
  54. T. Yanai, M. Saitow, X.-G. Xiong, J. Chalupský, Y. Kurashige, S. Guo and S. Sharma, J. Chem. Theory Comput., 2017, 13, 4829–4840 CrossRef CAS PubMed .
  55. M. Saitow, Y. Kurashige and T. Yanai, J. Chem. Phys., 2013, 139, 044118 CrossRef PubMed .
  56. T. Yanai, Y. Kurashige, E. Neuscamman and G. K. Chan, J. Chem. Phys., 2010, 132, 024105 CrossRef PubMed .
  57. P. Celani and H.-J. Werner, J. Chem. Phys., 2000, 112, 5546–5557 CrossRef CAS .
  58. P.-Å. Malmqvist, K. Pierloot, A. R. M. Shahi, C. J. Cramer and L. Gagliardi, J. Chem. Phys., 2008, 128, 204109 CrossRef PubMed .
  59. Q. M. Phung, S. Wouters and K. Pierloot, J. Chem. Theory Comput., 2016, 12, 4352–4361 CrossRef CAS PubMed .
  60. S. Wouters, V. Van Speybroeck and D. Van Neck, J. Chem. Phys., 2016, 145, 054120 CrossRef PubMed .
  61. N. Nakatani and S. Guo, J. Chem. Phys., 2017, 146, 094102 CrossRef .
  62. K. G. Dyall, J. Chem. Phys., 1995, 102, 4909–4918 CrossRef CAS .
  63. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J.-P. Malrieu, J. Chem. Phys., 2001, 114, 10252–10264 CrossRef CAS .
  64. S. Guo, M. A. Watson, W. Hu, Q. Sun and G. K.-L. Chan, J. Chem. Theory Comput., 2016, 12, 1583–1591 CrossRef CAS PubMed .
  65. M. Roemelt, S. Guo and G. K.-L. Chan, J. Chem. Phys., 2016, 144, 204113 CrossRef PubMed .
  66. S. Sharma, G. Knizia, S. Guo and A. Alavi, J. Chem. Theory Comput., 2017, 13, 488–498 CrossRef CAS PubMed .
  67. L. Freitag, S. Knecht, C. Angeli and M. Reiher, J. Chem. Theory Comput., 2017, 13, 451–459 CrossRef CAS PubMed .
  68. A. Y. Sokolov, S. Guo, E. Ronca and G. K.-L. Chan, J. Chem. Phys., 2017, 146, 244102 CrossRef PubMed .
  69. E. Xu and S. Li, J. Chem. Phys., 2013, 139, 174111 CrossRef PubMed .
  70. R. J. Anderson, T. Shiozaki and G. H. Booth, J. Chem. Phys., 2020, 152, 054101 CrossRef CAS PubMed .
  71. S. Sharma and A. Alavi, J. Chem. Phys., 2015, 143, 102815 CrossRef PubMed .
  72. Z. Luo, Y. Ma, X. Wang and H. Ma, J. Chem. Theory Comput., 2018, 14, 4747–4755 CrossRef CAS PubMed .
  73. C. Angeli, R. Cimiraglia and J.-P. Malrieu, J. Chem. Phys., 2002, 117, 9138–9153 CrossRef CAS .
  74. Y. Guo, K. Sivalingam and F. Neese, J. Chem. Phys., 2021, 154, 214111 CrossRef CAS PubMed .
  75. Y. Guo, K. Sivalingam, C. Kollmar and F. Neese, J. Chem. Phys., 2021, 154, 214113 CrossRef CAS PubMed .
  76. C. Kollmar, K. Sivalingam, Y. Guo and F. Neese, J. Chem. Phys., 2021, 155, 234104 CrossRef CAS PubMed .
  77. D. Zgid, D. Ghosh, E. Neuscamman and G. K. Chan, J. Chem. Phys., 2009, 130, 194107 CrossRef PubMed .
  78. F. A. Evangelista, J. Chem. Phys., 2014, 141, 054109 CrossRef PubMed .
  79. J. B. Schriber, K. P. Hannon, C. Li and F. A. Evangelista, J. Chem. Theory Comput., 2018, 14, 6295–6305 CrossRef CAS PubMed .
  80. J. W. Park, J. Chem. Theory Comput., 2023, 19, 6263–6272 CrossRef CAS PubMed .
  81. K. Pernal, Phys. Rev. Lett., 2018, 120, 013001 CrossRef CAS PubMed .
  82. E. Pastorczak and K. Pernal, J. Chem. Theory Comput., 2018, 14, 3493–3503 CrossRef CAS PubMed .
  83. K. Pernal, J. Chem. Phys., 2018, 149, 204101 CrossRef PubMed .
  84. E. Pastorczak, M. Hapka, L. Veis and K. Pernal, J. Phys. Chem. Lett., 2019, 10, 4668–4674 CrossRef CAS PubMed .
  85. D. Drwal, P. Beran, M. Hapka, M. Modrzejewski, A. Sokół, L. Veis and K. Pernal, J. Phys. Chem. Lett., 2022, 13, 4570–4578 CrossRef CAS PubMed .
  86. D. Drwal, E. Pastorczak and K. Pernal, J. Chem. Phys., 2021, 154, 164102 CrossRef CAS PubMed .
  87. P. Beran, M. Matoušek, M. Hapka, K. Pernal and L. Veis, J. Chem. Theory Comput., 2021, 17, 7575–7585 CrossRef CAS PubMed .
  88. M. Matoušek, M. Hapka, L. Veis and K. Pernal, J. Chem. Phys., 2023, 158, 054105 CrossRef PubMed .
  89. P. Pulay, S. Saebø and W. Meyer, J. Chem. Phys., 1984, 81, 1901–1905 CrossRef CAS .
  90. Y. Guo, K. Sivalingam, E. F. Valeev and F. Neese, J. Chem. Phys., 2016, 144, 094111 CrossRef PubMed .
  91. C. Angeli and R. Cimiraglia, J. Phys. Chem. A, 2014, 118, 6435–6439 CrossRef CAS PubMed .
  92. D. J. Rowe, Rev. Mod. Phys., 1968, 40, 153–166 CrossRef .
  93. E. Rosta and P. R. Surján, J. Chem. Phys., 2002, 116, 878–890 CrossRef CAS .
  94. K. Pernal, J. Chem. Theory Comput., 2014, 10, 4332–4341 CrossRef CAS PubMed .
  95. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73 CAS .
  96. F. Neese, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2022, 12, e1606 CrossRef .
  97. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108 CrossRef CAS PubMed .
  98. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023 CrossRef .
  99. M. Schreiber, M. R. Silva-Junior, S. P. A. Sauer and W. Thiel, J. Chem. Phys., 2008, 128, 134110 CrossRef PubMed .
  100. N. B. Balabanov and K. A. Peterson, J. Chem. Phys., 2005, 123, 064107 CrossRef PubMed .
  101. G. E. Scuseria, T. M. Henderson and I. W. Bulik, J. Chem. Phys., 2013, 139, 104113 CrossRef PubMed .
  102. B. O. Roos and K. Andersson, Chem. Phys. Lett., 1995, 245, 215–223 CrossRef CAS .
  103. H. R. Larsson, H. Zhai, C. J. Umrigar and G. K.-L. Chan, J. Am. Chem. Soc., 2022, 144, 15932–15937 CrossRef CAS PubMed .
  104. I. Schapiro, K. Sivalingam and F. Neese, J. Chem. Theory Comput., 2013, 9, 3567–3580 CrossRef CAS PubMed .
  105. Y. Song, Y. Guo, Y. Lei, N. Zhang and W. Liu, Top. Curr. Chem., 2021, 379, 1–56 CrossRef PubMed .
  106. S. Manna, R. K. Chaudhuri and S. Chattopadhyay, J. Chem. Phys., 2020, 152, 244105 CrossRef CAS PubMed .
This journal is © The Royal Society of Chemistry 2024