An alternative choice of the zeroth-order Hamiltonian in CASPT2 theory.

A zeroth-order Hamiltonian based on Koopmans matrices for complete active space second-order perturbation theory (CASPT2) is presented. This Hamiltonian involves three types of Fock matrices. The original CASPT2 Fock matrix is retained for excitation classes where the excitation does not change the number of electrons in the complete active space (CAS). For excitation classes involving a change in the number of electrons in the CAS, two alternative Fock matrices corresponding to either positive or negative ionization of the CAS are introduced. These are constructed such that they exactly reproduce the Koopmans matrices for a singly ionized CAS. Test calculations indicate that the method gives better excitation energies than CASPT2 without using empirical parameters, for example, the ionization potential-electron affinity shift, which is also designed to improve excitation energies. The method is also less prone to intruder states than conventional CASPT2. Moreover, the dissociation curve for the chromium dimer looks much more reasonable than the one obtained with conventional CASPT2.


I. INTRODUCTION
The complete active space self-consistent-field (CASSCF) method 1 describing static electron correlation represents the starting point for obtaining electronic wave functions with multiconfigurational character. These become important for nearly degenerate ground states occurring in the formation or breaking of chemical bonds or transition metal complexes. As in the single reference case, various methods exist for adding dynamic electron correlation on top of a reference wave function. Multireference configuration interaction (MRCI) both uncontracted 2,3 and contracted 4-7 is based on the addition of doubly excited configurations to the reference wave function but lacks size consistency. Multireference coupled cluster (MRCC) theory, [8][9][10][11][12][13][14] on the other hand, is both complex and costly and not yet ready for application to larger systems. Thus, approaches based on perturbation theory represent a viable alternative for the treatment of dynamic electron correlation in the multireference case. The most important of these are complete active space second-order perturbation theory (CASPT2) 15,16 and n-electron valence space perturbation theory (NEVPT2).. [17][18][19] In contrast to NEVPT2, CASPT2 is plagued by the problem of intruder states arising from near-zero energy denominators in the perturbation expansion with Cr 2 being a prominent example where these states lead to several discontinuities in the dissociation curve. 20 The most obvious remedy for avoiding intruder states is level shifts increasing the energy denominators. 21,22 Another modification of the CASPT2 zeroth-order Hamiltonian aims at a more balanced treatment of closed-shell and open-shell configurations. 23 A more specific level shift technique adapted to the ionic nature of the CAS for perturber functions belonging to excitation classes where electrons are added to or removed from the CAS has also been suggested. 24 These shifts are denoted as ionization potential-electron affinity (IPEA) shifts and depend on the occupation numbers of the active orbitals and involve an empirical parameter. Based on minimization of errors for dissociation and excitation energies, a default value of 0.25 a.u. has been suggested for this parameter. 24 However, there is no unique choice fitting all purposes. Studies on spin crossover phenomena in organometallic complexes are controversial, suggesting a much larger value in the interval 0.5-0.7 a.u. [25][26][27] or the default value. 28 Evaluation of magnetic coupling constants in heterobinuclear complexes found a value of zero, i.e., the original CASPT2 method to be most appropriate. 29 The results for the spin energetics of Fe II and Fe III heme complexes employing IPEA shifts of 0.25 a.u. and 0.5 a.u. are ambiguous. 30 The effect of the IPEA shift on the dissociation curve of the chromium dimer has also been studied. 31,32 While Ruipérez et al. 31 found a value of 0.45 a.u. to be best The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp suited for the reproduction of the experimental values in the case of a (12,12) CAS, Vancoillie et al. 32 observed that the dependence on the IPEA shift is strongly reduced if the CAS is increased. A critical analysis 33 also demonstrates that the application of IPEA shifts is not always appropriate.
In the single reference case, the diagonal elements of the Fock operator used in MP2 theory correspond to negative ionization potentials and electron affinities for doubly occupied and empty orbitals, respectively. Analogously, one would expect the energies ascribed to active orbitals to correspond to electron affinities if an electron is excited from the doubly occupied to the active space and to ionization potentials if an electron is excited from the active to the virtual space. This, however, is not the case for the Fock matrix applied in the CASPT2 zeroth-order Hamiltonian. Considering an active orbital with an occupation number of one, the corresponding MO energy is approximately related to the arithmetic mean of an ionization potential and an electron affinity, 24 leading to a lowering of the energy denominators and thus possibly to intruder states in both cases. The IPEA shifts are chosen such that they shift the active MO energies toward electron affinities if an electron is excited into the CAS and toward ionization potentials if an excitation removes an electron from the CAS. The approach suggested in this contribution is, in principle, based on the same idea but is intended to avoid the use of a more or less empirical parameter, which has to be varied within a wide range depending on the nature of the experimental values one wants to obtain. It is well-known that electron affinities and ionization potentials can be evaluated on the basis of extended Koopmans' theorem. [34][35][36][37] Starting from a CASSCF wave function, it has been demonstrated that such an approach gives vertical ionization potentials in good agreement with the experimental values in the context of NEVPT2. 38 Thus, we think that it could be useful to construct a single particle operator with the elements of the Fock matrix chosen such that this operator reproduces the Koopmans matrices for corresponding excitations of an electron into or out of the CAS. The two Fock matrices for positive and negative ionization of the CAS are then used in the CASPT2 formalism for the corresponding excitations. The details of this approach are described in Sec. II.
Test calculations have been performed to demonstrate the usefulness of the new approach. We have used a subset of the Thiel benchmark set 39 for the evaluation of excitation energies. As the chromium dimer is frequently used as a test case for any multireference method, the dissociation curve for this molecule has also been evaluated. Finally, we compare the spectroscopic constants for some simple diatomic molecules as obtained from several CASPT2 variants and NEVPT2.

II. THE CONSTRUCTION OF FOCK MATRICES FOR THE CASPT2 ZEROTH-ORDER HAMILTONIAN FROM KOOPMANS MATRICES
Since the Koopmans matrices play a central role in the development of the present formalism, we would like to briefly remind the reader to their definition. Given the wave function |Ψ N ⟩ as a solution of the Schrödinger equation for a system with N electrons, extended Koopmans' theorem [34][35][36][37] is motivated by the idea that the wave functions and corresponding energies |Ψ N +1 ⟩, E N +1 or |Ψ N −1 ⟩, E N −1 for the same system with, respectively, one electron added to or removed from it might be calculated using only the first-order and second-order density matrices of the N-electron system without solving the Schrödinger equation for N + 1-or N − 1-electron systems from scratch. Given an orbital basis set ψq, the wave functions of these systems are represented as (1) σ is the spin orientation of the electron to be added or removed. The variational parameters to be optimized are the coefficients c + μpσ and c − μpσ , which are assumed to be real. They are obtained according to the stationary conditions δ ± μσ = 0 with where the wave function |Ψ N ⟩ is assumed to be normalized. Inserting Eq. (1) in Eq. (2) and assuming that |Ψ N ⟩ is an exact eigenfunction of the Schrödinger equation, i.e.,Ĥ|Ψ N ⟩ = E N |Ψ N ⟩, leads to with the Koopmans matrices and the metric Note that the Hermiticity of the Koopmans matrices is not apparent from their formal definition given in Eq. (4). In fact, it depends on the assumption that |Ψ N ⟩ is an exact eigenfunction of the Schrödinger equation, which only holds for a full CI wave function.
Writing the Born-Oppenheimer Hamiltonian aŝ one finally obtains for the Koopmans matrices Summing the Koopmans matrices over spin and introducing the particle excitation operatorÊ one obtains In the following, we would like to give the Koopmans matrices for a CASSCF wave function where one has to distinguish three orbital subspaces. Indices i, j, . . ., t, u, . . ., and a, b, . . . are used for internal (doubly occupied), active, and virtual orbitals, respectively. Using p, q, . . . as general indices, one obtains with is the CASSCF wave function. Note that the Koopmans matrices, as given by Eq. (11), are not Hermitian, e.g., K − it ≠ K − ti , because the CASSCF wave function is of course not an exact eigenfunction of the Born-Oppenheimer Hamiltonian. One may note, however, that the matrix elements of their active-active block, i.e., K − tu and K + tu , appear in some of the matrix elements of NEVPT2 theory. Using the Dyall Hamiltonian 40 with the constant C in Eq. (13) chosen such that the zeroth-order energy corresponds to the energy of the CASSCF wave function, i.e., with ⟨Ψ t i =Êti Ψ (0) ⟩ and ⟨Ψ a t =Êat Ψ (0) ⟩. Note that the Dyall Hamiltonian contains a bielectronic term, which leads to the appearance of Koopmans matrices in the matrix elements between excited configurations in Eq. (14). CASPT2, on the other hand, employs a zeroth-order Hamiltonian based on the single particle Fock operator where the matrix F is given by Eq. (12). Writing the CASPT2 zerothorder Hamiltonian aŝ withQ andP sd representing projection operators on the orthogonal complement of |Ψ (0) ⟩ in the CAS configuration interaction (CASCI) space and the space spanned by internally contracted configuration state functions obtained from the application of single and double excitation operators to the CAS-CI wave function |Ψ (0) ⟩, respectively, one obtains with Comparing Eqs. (14) and (17), one might wish to replace the Fock matrix elements in the latter by other ones chosen such that the lefthand sides of these equations become identical. Denoting these Fock matrices as F + and F − , one obtains Equation (19) represents two systems of linear equations for the determination of the matrix elements of F + and F − , which are thus chosen such that they reproduce the Koopmans matrices in the matrix elements of excited configurations. Here, we have chosen singly excited configurations for the sake of simplicity, but exactly the same Koopmans matrices appear also in the matrix elements of doubly excited configurations: ⟩. Since the eigenvalues of the Koopmans matrices K + and K − are related to electron affinities and ionization potentials, respectively, it seems obvious that using F + for excitation classes where electrons are added to the CAS and F − for those where electrons are removed from the CAS leads to increased energy denominators and thus possibly to the avoidance of intruder states. Remember that the original Fock matrix represents a kind of compromise between electron affinities and ionization potentials because it results from the arithmetic mean of the Koopmans matrix K + with eigenvalues representing negative ionization potentials and the negative of the Koopmans matrix K − , the eigenvalues of which correspond to electron affinities, The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp Thus, the energy levels of the active space resulting from this choice can be expected to be too low if an electron is excited into the active space and too high if an electron is excited from the active to the virtual space. This clearly favors the appearance of intruder states. Taking into account the Hermiticity of the Koopmans and Fock matrices, the system of linear equations [Eq. (19)] can be written as with In general, the matrix X is not definite, and its eigenvalues xi may become near singular. This is addressed by a process known as Tikhonov regularization 41 where the inverse eigenvalues are modified using a small parameter α = 10 −2 as follows: Thus, numerical instabilities are avoided, and the solution of Eq. (21) for the matrix elements of F + and F − by inversion of the matrix X does not present any difficulties.

III. NUMERICAL RESULTS
The modified CASPT2 approach described in Sec. II has been implemented in the ORCA program package, 42 which has been used for all results presented in the following. It was our main goal to test the new method denoted as CASPT2-K in the following in cases that are sensitive to intruder states.

A. Excitation energies for a subset of the Thiel benchmark set
A benchmark set of medium-sized organic molecules for testing vertical excitation energies has been compiled by Schreiber et al. 39 This so-called Thiel benchmark set has also been used by Schapiro et al. 43 to probe the excitation energies obtained from statespecific NEVPT2 calculations in comparison with CASPT2. In our calculations, we have a chosen subset of the complete Thiel benchmark set based on a preference for molecules where the CASPT2 method tends to be sensitive to intruder states. These molecules have been identified using the data given by Zobel et al. 33 in the supplementary material of their paper. As we wanted to avoid the use of arbitrary parameters as far as possible, level shifts are only applied if absolutely necessary, i.e., in the case of intruder states that are identified by their reference weight. The setup of our calculations (basis set, frozen core, active space, geometry, symmetry, and weighting factors of state averaging) is identical to the one given by Schreiber et al. 39 and also used by Shapiro et al. 43 In particular, the same TZVP basis set has been applied. 44 In analogy to the work of Shapiro et al., 43 results from an iterative approximate coupledcluster triples model (CC3) are used as a reference because it has been demonstrated that excitation energies obtained from CC3 are close to full CI values. 45 Table I shows results for NEVPT2, CASPT2, the present method denoted as CASPT2-K, and CASPT2-IPEA with the standard value of 0.25 for the shift parameter. A comment on our implementation of the IPEA shift might be appropriate at this point. In the original work of Ghigo et al. 24 the IPEA shift is added to the diagonal elements of the Fock matrix. We followed a different method given by Kepenekian et al. 25 These authors apply the shift given in Table I of their paper for various excitation classes to the diagonal matrix elements ⟨Ψ rs pq Ĥ 0 − E (0) Ψ rs pq ⟩ of the perturber functions, withĤ 0 being the CASPT2 zeroth-order Hamiltonian. Thus, our IPEA results deviate slightly from those obtained from the MOL-CAS program package. 46 The literature on this subject is unclear to us.
Signed deviations from the CC3 excitation energies are given in Table I for the methods mentioned above. The reference weight of the ground state amounts roughly to 0.7-0.8 for the molecules considered. If the reference weight of an excited state has been found to be more than 10% lower than that of the ground state, an imaginary level shift 22 of 0.2 a.u. has been added for the affected states. The corresponding excitation energies are given in italics in Table I.    Table I. It can be seen that the mean signed error is much smaller than the mean unsigned error for NEVPT2, which clearly performs better than the CASPT2 variants. The absolute values of the maximum errors are 0.92 eV, 3.14 eV, 1.92 eV, and 0.87 eV for NEVPT2, CASPT2, CASPT2-K, and CASPT2-IPEA, respectively. The corresponding states are plagued by intruders in of Chemical Physics the case of CASPT2 and CASPT2-K. The CASPT2 variants, especially CASPT2 and CASPT2-K, tend to overestimate the excitation energies. CASPT2-K performs better than CASPT2 but worse than CASPT2-IPEA. However, in contrast to CASPT2-IPEA, it does not involve any empirical parameter. Remember that the default value of the IPEA shift has been chosen to reproduce excitation energies and ionization potentials. 24 Other properties require an adjustment of this parameter. It has already been discussed in the Introduction that the variation can be rather large.
We are confident that our conclusions also hold for the remainder of Thiel's benchmark set, which will be the subject of future studies.

B. The chromium dimer
The chromium dimer is considered as a benchmark molecule for the performance of multireference methods. Various studies using multiconfiguration perturbation theory have been reported. 20,31,32,[47][48][49][50][51] It is therefore of interest to investigate the performance of the new method for this molecule. The dissociation curve for the chromium dimer has been calculated for NEVPT2 and several CASPT2 variants including conventional CASPT2, CASPT2-IPEA with IPEA shifts of 0.25 and 0.5, and CASPT2-K. The minimal (12,12) CAS consists of the 3d and 4s orbitals. The DKH-def2-TZVP and DKH-def2-QZVPP basis sets 44 with the DKH2 Hamiltonian 52,53 have been used. The reported results have been counterpoise corrected for the basis set superposition error (BSSE). An imaginary level shift of 0.2 a.u. has been employed in all CASPT2 calculations to avoid intruder states. The frozen core consists of all orbitals up to and including the 2p shells. The results for the DKH-def2-TZVP and DKH-def2-QZVPP basis sets have been extrapolated to the complete TABLE II. Binding energy De and equilibrium bond distance r 0 for the chromium dimer as obtained from various perturbational multireference methods for the DKH-def2-TZVP and DKH-def2-QZVPP basis sets and the extrapolated basis set limit.  The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp basis set (CBS) limit 54,55 using extrapolation parameters from the def2 basis. 56 The results for the calculated energy surface are shown in Figs. 1-3 for the DKH-def2-TZVP basis set, DKH-def2-QZVPP basis set, and CBS, respectively. These plots show the binding energy, i.e., the difference between the energy of the molecule and two single chromium atoms rather than the total electronic energy. This simplifies the comparison with the experimental values taken from the work of Casey and Leopold, 57 which are also shown. The experimental binding energy of De = 1.56 eV has been adopted from the work of Simaud et al. 58 Looking at the dissociation curves, one observes that, in contrast to conventional CASPT2, the qualitative shape of the curve obtained from CASPT2-K looks reasonable. The CASPT2 curve can be improved by varying the magnitude of the IPEA shift with an optimum value near 0.5. However, this value differs significantly from the standard value of 0.25. Note that CASPT2-K does not rely on such a parameter and is hence preferable. Moreover, the curve obtained with CASPT2-K for the DKH-def2-TZVP basis without any level shift (not shown here) is practically indistinguishable from the one shown in Fig. 1 and thus, in contrast to the other CASPT2 TABLE III. Spectroscopic constants for some diatomic molecules as obtained from various methods with the QZVP basis set. Experimental values are also shown. The energies at the equilibrium bond length are given in atomic units, and bond lengths re are given in Ångström, whereas the equilibrium rotational constant Be, the vibration-rotation interaction constant αe, the harmonic frequency ωe, and the anharmonicity ωeχe are given in cm −1 . The dissociation energies De are given in eV.

Method
Energy The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp variants, free of discontinuities resulting from intruder states. This indicates once more that CASPT2-K is less sensitive to intruder states. Table II shows the binding energy and equilibrium bond lengths for various perturbation methods. The quantitative agreement with the experimental values of r 0 = 1.679 Å 57 and De = 1.56 eV 58 is not particularly impressive, but at least the equilibrium bond lengths obtained for NEVPT2, CASPT2-K, and CASPT2-IPEA with IPEA = 0.5 are quite reasonable. Comparing NEVPT2 and CASPT2-IPEA with IPEA = 0.5, the binding energy of the latter is closer to the experimental value. However, the qualitative shape of the curve is correct for NEVPT2, whereas CASPT2-IPEA with IPEA = 0.5 shows a second minimum instead of a shoulder. The curve for CASPT2-K also has the correct shape but strongly overestimates the binding energy. It is known from previous NEVPT2 and CASPT2 studies with the minimal active space that the results are strongly basis set dependent. Better agreement with the experimental curve is achieved by extending the active space with the second d-shell, which leads to a CAS (12,22). 32,59 Such large active spaces are not tractable with the exact CASSCF reference wave function and are thus not pursued here.

C. Other diatomic molecules
Potential energy surfaces have been computed for the diatomic molecules N 2 , CO, F 2 , O 2 , HF, and BH with various methods using the def2-QZVP basis set. 44 No level shift is applied in the CASPT2 variants except for the standard IPEA shift of 0.25 in CASPT2-IPEA. The 1s orbitals are treated as a frozen core. Table III shows calculated spectroscopic constants obtained from CASSCF, NEVPT2, CASPT2, CASPT2-K, and CASPT2-IPEA. The experimental values shown were taken from the work of Herzberg 60 except for the dissociation energies, which are tabulated in the handbook of chemistry and physics. 61 It can be seen that all dynamic correlation methods lead to better agreement with the experimental results as compared to CASSCF. Deviations of the bond lengths of more than 0.01 Å from the experimental values have only been observed for F 2 where the energy curve is rather shallow, as indicated by the low variational frequency. Most interesting are the systematic trends in the dissociation energy. NEVPT2 yields larger values of the total energy than the CASPT2 variants in all cases. The effect of the IPEA shift on the CASPT2 results is, in general, small in the vicinity of the minimum but shows up clearly near the dissociation limit where it leads to an upshift of the energy so that the IPEA shift always increases the binding energy, as has already been observed by Ghigo et al. 24 It brings the calculated values closer to the experimental ones for N 2 and CO. For all other molecules, however, the CASPT2 value is already too large so that the upshift of the energy worsens the results and leads to an overestimation of the binding energies. This contradicts the findings of Ghigo et al. 24 who observed an improvement for almost all molecules of their larger test set. This may be due to the larger basis set used by these authors. The binding energies obtained from CASPT2-K are, in general, even larger than the ones resulting from CASPT2-IPEA. The same holds for NEVPT2. Thus, both methods tend to overestimate the binding energy. The trends just described can be verified by looking at the dissociation curves and energy tables given in the supplementary material of this article.

IV. CONCLUSIONS
A modified CASPT2 formalism has been presented. Based on the idea that energy levels in the active space should resemble electron affinities if an excitation adds electrons to the CAS and ionization potentials if electrons are removed from the CAS, the Fock matrix elements needed for the CASPT2 zeroth-order Hamiltonian are constructed from the corresponding Koopmans matrices. As a consequence, the correlation energies for the excitation classes corresponding to states ÊaiÊtj|Ψ (0) ⟩ and ÊaiÊ bt |Ψ (0) ⟩ are the same for NEVPT2 and CASPT2D, i.e., CASPT2 without interclass interaction. This has been confirmed numerically with small deviations arising from application of the Tikhonov regularization. Our calculations indicate that this method is less sensitive to intruder states than the standard CASPT2 formalism with and without IPEA shift. We also found that our approach leads to an improvement of excitation energies, as indicated by calculations on a subset of the Thiel benchmark set 39 as compared to conventional CASPT2. The dissociation curve for the chromium dimer represents a challenge for any multireference method. The shape of this curve is considerably improved if the present method is used instead of conventional CASPT2.
Since the modifications of the Fock matrix in CASPT2-K and CASPT2-IPEA are motivated basically by the same idea, it is perhaps not surprising that the changes brought about by these manipulations show similar trends in the results. Both methods improve the excitation energies with respect to conventional CASPT2 and lead to an increase in the dissociation energies for diatomic molecules. The excitation energies obtained from CASPT2-IPEA with the default value for the IPEA shift are still better than those resulting from CASPT2-K, but it should be noted that CASPT2-K is free of adjustable parameters. Remember also that the IPEA shift had to be doubled with respect to its standard value to give a dissociation curve for Cr 2 , showing at least some resemblance to the experimental one, whereas the CASPT2-K curve has a qualitatively correct shape without such manipulations. Thus, overall, we consider CASPT2-K to be a modest but clear improvement over the more empirical CASPT2-IPEA.

SUPPLEMENTARY MATERIAL
The following data are accessible in the supplementary material to this article: (a) a table giving the weight of the reference state of the perturbed wave function for all ground and excited states of the molecules appearing in Table I, (b) tables containing the bond distances and corresponding total electronic energies and binding energies for the chromium dimer as obtained from various perturbation methods for the DKH-def2-TZVP and DKH-def2-QZVPP basis sets and the complete basis set (CBS) limit, (c) two figures showing the influence of the imaginary level shift on the dissociation curves of the chromium dimer for the CASPT2 variants, and (d) tables and plots giving the total electronic energy as a function of the bond distance for all diatomic molecules described in Sec. III C.