Large-Scale Benchmarking of Multireference Vertical-Excitation Calculations via Automated Active-Space Selection

We have calculated state-averaged complete-active-space self-consistent-field (SA-CASSCF), multiconfiguration pair-density functional theory (MC-PDFT), hybrid MC-PDFT (HMC-PDFT), and n-electron valence state second-order perturbation theory (NEVPT2) excitation energies with the approximate pair coefficient (APC) automated active-space selection scheme for the QUESTDB benchmark database of 542 vertical excitation energies. We eliminated poor active spaces (20–40% of calculations) by applying a threshold to the SA-CASSCF absolute error. With the remaining calculations, we find that NEVPT2 performance is significantly impacted by the size of the basis set the wave functions are converged in, regardless of the quality of their description, which is a problem absent in MC-PDFT. Additionally, we find that HMC-PDFT is a significant improvement over MC-PDFT with the translated PBE (tPBE) density functional and that it performs about as well as NEVPT2 and second-order coupled cluster on a set of 373 excitations in the QUESTDB database. We optimized the percentage of SA-CASSCF energy to include in HMC-PDFT when using the tPBE on-top functional, and we find the 25% value used in tPBE0 to be optimal. This work is by far the largest benchmarking of MC-PDFT and HMC-PDFT to date, and the data produced in this work are useful as a validation of HMC-PDFT and of the APC active-space selection scheme. We have made all the wave functions produced in this work (orbitals and CI vectors) available to the public and encourage the community to utilize this data as a tool in the development of further multireference model chemistries.


Introduction
The accurate treatment of excited states is critical for understanding photochemical phenomena, [1][2][3][4][5][6][7] and it has been a long-standing goal of the electronic structure community. [8][9][10][11][12][13][14][15][16][17][18][19] Although treating excited states is difficult in general, it is particularly challenging when single-determinant methods such as Hartree-Fock or Kohn-Sham density functional theory provide a poor reference state for predicting excited states. This can occur either because the excited states vary greatly from the ground state (e.g., double excitations 20 ) or because the ground state itself is not well-described (e.g., strongly correlated systems [21][22][23][24] ). One can overcome these deficiencies by using multiple-determinant reference states, and the methods that take this approach are called multireference methods.
The most popular multireference method is the complete active space self-consistent field (CASSCF) method, 25 which expresses approximate wave functions in the space of all possible configurations of electrons in an "active space" of orbitals and electrons. These wave functions can then serve as references for perturbation theories such as MC-QDPT, 26,27 CASPT2, 28,29 and NEVPT2. 30,31 Alternatively, quantitative accuracy can be achieved by using a nonclassical-energy functional applied to the the converged wave function in multiconfiguration nonclassical functional theory (MC-NCFT). 32- 36 The total energy is then a sum of the classical portion of the CASSCF energy and nonclassical energy from the functional.
The most common form of MC-NCFT utilizes nonclassical-energy functionals obtained by translating Kohn-Sham exchange-correlation functionals for use with multiconfigurational wave functions via the on-top pair density and is called multiconfiguration pair-density functional theory (MC-PDFT). The translated PBE functional (tPBE) has been used as the functional in the majority of MC-PDFT calculations to date. The nonclassical energy from a density functional can be mixed with the nonclassical part of the CASSCF energy to form a "hybrid" nonclassical functional, for example, using a 0.75:0.25 mixture of tPBE and CASSCF nonclassical energies yields the tPBE0 functional. 37 Difficulties encountered in all such post-CASSCF methods are making the active space large enough and well-balanced enough to converge to a qualitatively accurate description of the underlying wave function(s). The results can depend significantly on the size and nature of the active space and the initial orbital guess. 38 Moreover, in many occasions the orbitals will change character during their optimization. For these reasons, such calculations often require expert human guidance to carefully choose the active space size and composition and monitor them during the calculations.
Recently, Jacquemin and coworkers published the QUESTDB benchmark dataset of 542 vertical excitation energies on a diverse set of small and midsize main-group molecules, calculated via a variety of high-level wave function methods in the aug-cc-pVTZ 60,61 basis. 20,[62][63][64][65][66] In the present paper we have undertaken the automated calculation of these excitation energies with SA-CASSCF, NEVPT2, and MC-PDFT using the APC-ranked-orbital active space selection scheme. To benchmark and analyze the performance of various multireference methods on this diverse set of excitations, we eliminate poor active spaces (20-30% of calculations) by setting an error threshold on the SA-CASSCF excitation energy because that has previously been shown to be good way to judge the quality of the active space. 49 By analyzing results across different active space and basis set size choices, we find different trends in the performance of MC-PDFT and NEVPT2 where the performance of NEVPT2 is overly dependent on the basis set in which the underlying wave function is converged.
Additionally, we are able to produce the first large-scale and robust comparison of MC-PDFT to other single-reference methods such as CC2, and find the CASSCF mixing parameter of 0.25 used in tPBE0 to be optimal. We have made all the wave functions converged in this work available to the public via publication of the converged orbitals and CI vectors and encourage others to use these data in the development of further multireference model chemistries.

Methods
Data Overview. The data we have examined can be found in the QUESTDB dataset, 65 which consists of 542 vertical excitations of small and midsize main-group molecules (molecules with 1-10 non-hydrogenic atoms). Of these excitations, 491 are from singlet ground states and by its spatial and spin symmetries, and benchmark values are reported as "theoretical best estimates" (TBEs) calculated with a variety of high-level methods with the aug-cc-pVTZ 60,61 basis. These TBEs have been used in this work to judge the errors of all computed excitation energies, even those obtained with a different basis set.
Active Space Selection. To obtain orbitals for the active space selection scheme, we started with a restricted Hartree-Fock singlet wave function for closed-shell molecules and a restricted open-shell Hartree-Fock doublet wave function for doublet molecules, as calculated using PySCF. 67 The molecular point group was reduced to the highest available symmetry implemented for the PySCF SA-CASSCF solver: C 2h , C 2v , C s , or D 2h . The APC-rankedorbital active-space-selection scheme 36,59 starts with a set of candidate localized orbitals, ranks them by their approximated orbital entropies, and then eliminates orbitals starting from the lowest-entropy orbitals (those with the highest entropies are considered to be the most important) until the active space size reaches a predetermined maximum number of configuration state functions. We next describe the generation of candidate orbitals, then the ranking scheme, and finally the maximum-size criteria.
Following previous work, 36 up to 23 lowest-energy virtual orbitals of the Hartree-Fock calculation were selected, and orbitals within this subset were grouped by symmetry and Boys-localized 68 within each symmetry. Likewise, up to 23 highest-energy doubly occupied orbitals were also grouped by symmetry and Boys-localized within each symmetry. These two sets of localized orbitals (and the one singly occupied orbital, when present) were then considered as candidates for the active space. Next we describe how we ranked the localized orbitals.
In the originally published APC ranking scheme, given a set of doubly occupied candidate orbitals i and virtual orbitals a, one calculated the approximate-pair-coefficient (APC) matrix C ia as where K aa is the diagonal virtual element of the exchange matrix and F aa and F ii are diagonal virtual and doubly occupied elements of the Fock matrix in the MO basis. The entropies of doubly occupied orbitals i are calculated as and those of virtual orbitals a as However, in the present work we have found that in larger molecules (with >350 augcc-pVTZ basis functions) the APC entropies tend to overestimate the interaction of some virtual orbitals with the doubly occupied orbitals, artificially inflating the entropies of all doubly occupied orbitals and causing the selection of highly imbalanced active spaces (SI Figure 2). To overcome this issue, we propose an algorithmic extension of APC in which high-entropy virtual orbitals are removed from consideration when calculating entropies and then assigned the maximum entropy value (i.e. treated in the same way as singly occupied orbitals). The algorithm takes the following steps: • Provide the sets of candidate doubly occupied/ singly occupied/virtual orbitals and then remove the highest-entropy virtual orbital from L a and put it in L s . Repeat N times.
The above algorithm has a single parameter N (APC-N ) which is the number of times the highest-entropy virtual orbital is removed. In this work we have found a good value of N to be 2 (a scheme we refer to as APC-2), and we find that using APC-2 entropies results in more balanced active spaces and lower SA-CASSCF error than APC (SI Figure 1). The APC-2 entropies are then used to rank and select the orbitals for the active space by dropping the lowest-entropy orbitals until the active space size is lower than some maximum active space size. However, we note that the current scheme should be improved for the treatment of orbitals with degenerate entropies (SI Section 1.4).
In more detail, the active space size for a given set of N orb active orbitals containing N elec active electrons is calculated via the equation: 69 where N α + N β = N elec and N α = N β for even N elec and N α = N β + 1 for odd N elec . The maximum active space size N Max CSF is set via a specification of a maximum number of active electrons and orbitals (N Max elec ,N Max orb ) whose size is calculated via equation 4; this maximum active space choice is notated as max(N Max elec ,N Max orb ). In this work we calculate results at three choices of max(N Max elec ,N Max orb ): max (8,8) (N Max CSF = 1764), max(10,10) (N Max CSF = 19404), and max (12,12) (N Max CSF = 226512). Following this specification, all orbitals are selected and then the lowest-entropy orbital is successively dropped until the size of the active space calculated via equation 4 is less than or equal to N Max CSF .
As an example, we guide the reader through choosing a max(4,4) active space (N Max CSF = 20) from a set of orbitals with occupancies and entropies {(n j , S j )}: Calculation of the Excitation Energies. Calculations of excited-state wave functions were carried out by state-averaged CASSCF, averaging over the ground state and the minimum necessary number of excited states of the symmetry specified by QUESTDB. For example, in a C 2v molecule (e.g., water) if the symmetry of the excited state under consideration is specified to be 1 A 2 with no lower 1 A 2 excitations present, then the state averaging was done evenly over the 1 A 1 ground state and the 1 A 2 excited state. For higher 1 A 2 excitations, however, the state averaging included an additional 1 A 2 state for each 1 A 2 excitation lower in energy (again with weights for state averaging being the same for all states averaged).
Standard convergence parameters were employed, and for a few poor active space choices (0.4% of cases) the calculations failed to converge.
Because the highest available point groups supported by the SA-CASSCF solver in PySCF have lower symmetry than those specified by QUESTDB for single atoms and diatomics, the point groups sometimes had to be reduced to the highest-symmetry subgroup. Additionally, the labeling of different irreps is sometimes a choice of axis convention, such as between B 1 and B 2 in C 2v or between B 1g /B 2g /B 3g and B 1u /B 2u /B 3u in D 2h ; we have done our best to match the irrep we think was used in QUESTDB. Calculations were done with the highest M S allowed by the spin symmetry (e.g., if an excited state has S = 1, then for 8 electrons in the active space the active space would have 5 α and 3 β electrons).
The tPBE and NEVPT2 energies of the converged SA-CASSCF states were then calculated using the implementations of these methods in PySCF. Our implementation of MC-PDFT within PySCF is currently available in the mrh repository. 70 Additionally, tPBE0 energies were calculated by averaging the SA-CASSCF and tPBE energies: 37 The only implementation of NEVPT2 currently in PySCF is strongly contracted NEVPT2 (SC-NEVPT2), 31 and our NEVPT2 calculations use this.
In order to maintain a consistent labeling, the excited state to be compared to the QUESTDB excitation energy was chosen to be the state highest in energy as judged by tPBE. Although this is not a fail-proof scheme in terms of isolating the "same" QUESTDB state of the specified symmetry due to root flipping, we have found it to be satisfactory for our work as the converged QUESTDB wave functions are unavailable and labels such as "n → π * " are ambiguous non-observables. However, because the present work shows that tPBE0 is more accurate than tPBE, we suggest ordering the states by tPBE0 in future work. 3 Results

Eliminating Poor Active Spaces
We calculated excitation energies for all 542 vertical excitations listed in the QUESTDB database with six combinations of active space and basis set: four involving max(12,12) APC-2 active spaces with decreasing basis size (aug-cc-pVTZ, 60,61 jun-cc-pVTZ, 77 cc-pVTZ, 78,79 cc-pVDZ 78,79 ) and two involving jun-cc-pVTZ with decreasing active space size (max (10,10) and max (8,8)). We will refer to these combinations throughout the paper as Aug (12,12), Jun (12,12), TZ (12,12), DZ (12,12), Jun(10,10), and Jun (8,8).  that we obtain for each of the combinations of active space and basis set for all 3237 wave functions that we converged. As expected, we find that the SA-CASSCF error increases when we move from a larger to a smaller basis set with a given active space scheme or when we move from a larger active space to a smaller one with a given basis set. However, in order to reasonably evaluate the accuracies of these methods, we need to eliminate results whose error is driven mainly by poorly chosen active spaces. To analyze only cases with reasonable active spaces we set a threshold T on the SA-CASSCF error of 1.1 eV (T SA−CASSCF = 1.1 eV). That is, we consider that the APC scheme has produced a good active space if the error in the SA-CASSCF excitation energy is less than 1.1 eV.  Further discussion of the error threshold is given in the Supporting Information, which shows the 1.1 eV SA-CASSCF error threshold to be optimal (albeit imperfect) for isolating subsets of automated wave function calculations that reproduce results curated by hand. 84 Additionally, we analyze alternative error thresholds on the NEVPT2 and tPBE0 error. However, an error criterion cannot be used when a benchmark excitation energy or experimental excitation energy is not available. Nevertheless, when an accurate value is not available, one can still use this criterion (although with somewhat less reliably) by comparing to one's best estimate rather than to an accurate value. Clearly, if one's best estimate is good, this will work as well as comparing to an accurate value.

Comparison to Other Methods
Data Overview. In the QUESTDB database, 65    has nearly zero mean signed error. We note that the same good balance seems to occur in ADC(2.5), 85 which averages ADC(3) and ADC(2).  which is only slightly lower than the overall un-dropped-out fraction of 436/542 (80%) in the Aug (12,12) calculations. In keeping with the usual recommendation to use multireference methods for this class of excitation, we find that multireference methods are the only methods that perform consistently well on double excitations. Figure 6: Mean absolute errors (in eV) of tPBE, tPBE0, and NEVPT2 Aug (12,12) calculations on various types of S 0 excitations included by the threshold T SA-CASSCF = 1.1 eV. Figure 6 shows the errors classified by excitation type. In line with the Sarkar study, 84 we find that NEVPT2 is more accurate for triplet excitations than singlet excitations, and tPBE and tPBE0 follow this same trend. The figure shows that, with the exception of Rydberg states, tPBE0 has better performance than tPBE for every excitation category, and therefore we recommend the use of tPBE0 rather then tPBE for calculating excitation energies of valence excitations. We also recommend tPBE0 for calculating a spectrum containing both valence and Rydberg excitations since the performance of the two methods is very similar (on average) for Rydberg states.  correction at different active space/basis set sizes and compares these timings to those for the calculation of the tPBE on-top energy by the methodology in section 2. We find that at the normal grid size (grids_level = 3 in PySCF), tPBE is on average 114× less expensive than NEVPT2 for the large max (12,12) active spaces. This is because -as is well known -the cost of NEVPT2 scales very poorly with the size of the active space, while the cost of tPBE0 remains independent of that. Furthermore, the memory required for NEVPT2 also increases with active space size. It is around the max (12,12) active space size that the compute time for the perturbative correction begins to exceed the compute time of the underlying SA-CASSCF step, while the compute time of tPBE remains low. 86 For smaller active spaces such as max (8,8), the cost of NEVPT2 is comparable to that of tPBE and tPBE0.

Method Timing
Keeping the cost down is important in many applications. Figure 7 shows that the dependence of MC-PDFT compute times on grid size is a significant consideration; we observe a roughly 10× increase in cost from grids_level = 3 (3.4 s) to grids_level = 9 (30.8 s). Our studies find that the standard grids_level = 3 in PySCF is sufficient for excitations such as those we have calculated because we only see a significant change between the maximum and default grid size for a single excitation (SI Figure 22). Therefore we recommend standard grid sizes for most applications involving state-averaged MC-PDFT.

Optimizing the Mixing Parameter in Hybrid tPBE
A major motivation of this work was to generate data for benchmarking and improving MC-PDFT. As a first use of our data to optimize MC-PDFT functionals, we have investigated the optimal mixing parameter λ for hybrid tPBE (htPBE, for which the energy is given by    Figure 20); using the more robust tPBE0 error threshold (discussed in the SI) removes this discrepancy (SI Figure 21). This suggests that a higher value of λ may be optimal for cases in which wave function error dominates. 20 The work presented here is the largest application to date of automated multireference calculations on a broad range of molecules. The generation of 3237 multireference excitation energies has allowed us to gain insight into how to eliminate poorly chosen active spaces and has identified trends in the performance of MC-PDFT and NEVPT2. This work has been possible only through the careful work of Loos, Jacquemin, and coworkers in compiling the QUESTDB dataset 20,62,[64][65][66]87 and the recent work of of Sarkar et. al. 84 which has enabled us to compare our automatically generated results to hand-selected active space calculations.
We see this initial publication as laying the groundwork for several future applications related to MC-PDFT and high-throughput multireference calculations including: • Using the generated data to train and test novel functionals for MC-NCFT, representing a continuation of our initial work that used carbene singlet-triplet excitation energies to train machine-learned functionals. 36 • Improving the active space selection scheme. Our finding that error thresholds can be used to determine the fraction of poor wave functions in the calculated excitation energies can be used as a measure to benchmark the effectiveness of different active space selection schemes.
• Determining if a selected active space is well-chosen without reference to the underlying benchmark values. For specific active spaces and basis sets there appears to be promise in looking at differences between different methods (SI Figure 13), but a method that is generalizable across active spaces and basis sets has yet to be found.
Additionally, we expect that the wave functions converged in this work will be of interest for the development of different post-CASSCF methods such as multireference adiabatic connection (AC) 88 and algebraic diagrammatic construction (ADC). 89 For this reason, we are making all 3237 converged wave functions freely available for public use. We hope that this data will be useful to the electronic structure community both for comparing to the results published here and for developing and testing their own methods.
In summary, we have carried out the largest benchmarking of SA-CASSCF and MC-PDFT to date. This was accomplished by means of an automatic active-space selection scheme and use of a SA-CASSCF error threshold to eliminate poor active-space choices. On a set of 373 aug-cc-pVTZ excitation energies, we find that tPBE0 and NEVPT2 perform with similar accuracy to CC2, while tPBE lags behind. However, the accuracy of NEVPT2 degrades with basis set size even as the quality of the underlying density and on-top pair density appear to remain the same. As expected, we find that tPBE0 is orders of magnitude less expensive than NEVPT2 for larger active spaces, and we recommend its use for the calculation of a broad range of excitation energies, including double excitations.

Acknowledgment
The present work is supported by the National