# Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). III. Large active space calculations with selected configuration interaction reference

Yang Guo, Kantharuban Sivalingam, Vijay Gopal Chilkuri, Frank Neese

► To cite this version:

Yang Guo, Kantharuban Sivalingam, Vijay Gopal Chilkuri, Frank Neese. Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). III. Large active space calculations with selected configuration interaction reference. The Journal of Chemical Physics, 2025, 162 (14), 10.1063/5.0262473. hal-05081709

HAL Id: hal-05081709

https://hal.science/hal-05081709v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). III. Large active space calculations with selected configuration interaction reference

Yang Guo ID; Kantharuban Sivalingam ID; Vijay Gopal Chilkuri ID; Frank Neese ID

<figure><img src="image_2.png" /></figure>

Check for updates

J. Chem. Phys. 162, 144110 (2025)

https://doi.org/10.1063/5.0262473

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

## Articles You May Be Interested In

Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). II. The full rank NEVPT2 (FR-NEVPT2) formulation

J. Chem. Phys. (June 2021)

Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). I. Revisiting the NEVPT2 construction

J. Chem. Phys. (June 2021)

Analytic first-order derivatives of partially contracted *n*-electron valence state second-order perturbation theory (PC-NEVPT2)

J. Chem. Phys. (September 2019)

<figure><img src="image_5.png" /></figure>

AIP
Publishing

The Journal of Chemical Physics

Special Topics Open
for Submissions

Learn More

<figure><img src="image_6.png" /></figure>

AIP
Publishing
# Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). III. Large active space calculations with selected configuration interaction reference

Cite as: J. Chem. Phys. 162, 144110 (2025); doi: 10.1063/5.0262473

Submitted: 1 February 2025 • Accepted: 24 March 2025 •

Published Online: 11 April 2025

<figure><img src="image_7.png" /></figure>

View Online

<figure><img src="image_8.png" /></figure>

Export Citation

<figure><img src="image_9.png" /></figure>

CrossMark

Yang Guo,<fnref n="1" /> ID Kantharuban Sivalingam,<fnref n="2" /> ID Vijay Gopal Chilkuri,<fnref n="3" /> ID and Frank Neese<fnref n="2" /> a) ID

## AFFILIATIONS

## ABSTRACT

In multi-reference (MR) methods, addressing systems with large active spaces remains a challenge in the field. In Papers I and II of this series, we demonstrated that full rank N-electron valence state second-order perturbation theory (FR-NEVPT2) is a robust MR perturbation theory capable of computing strongly correlated systems with approximate density matrices. However, the previous FR-NEVPT2 implementation requires the computation and storage of fifth-order reduced density matrices (RDMs), limiting the usage of FR-NEVPT2 for systems with large active spaces. In the present work, as Paper III of the series, we report a new FR-NEVPT2 algorithm to handle systems with large active spaces. In the new algorithm, an approximate complete active space (CAS) self-consistent field (SCF) method, iterative configuration expansion (ICE) SCF, is employed to compute the reference wave functions for FR-NEVPT2. Then, the necessary Koopmans matrices of FR-NEVPT2 involving various RDMs are constructed using the intermediates designed by Kollmar et al. [J. Chem. Phys. 155, 234104 (2021)] to avoid storage bottlenecks. The performance of the new FR-NEVPT2 algorithm for systems with large active spaces is evaluated. Our results show that even with aggressive truncation parameters to truncate the ICE-SCF reference wave function, FR-NEVPT2 effectively recovers the missing static correlations of ICE-SCF. Several interesting systems with active spaces up to CAS(34,34) are studied using FR-NEVPT2 with ICE-SCF reference.

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0262473

## I. INTRODUCTION

It is well-known that the computational costs of the complete active space (CAS) self-consistent field (SCF) as well as the post-CASSCF methods increase exponentially with the number of active electrons and orbitals. To enhance the applicability of the CASSCF method to systems with a larger number of strongly correlated active electrons and orbitals, various approximate CAS configuration interaction (CASCI) methods have been developed. In 1969, Whitten and Hackmeyer introduced the first selected CI (sCI) approach to investigate the electronic structure of formaldehyde.<fnref n="1" /> During

the 1970s, Malrieu and Peyerimhoff introduced their own sCI methods,<fnref n="2" />⁻⁴ while Goddard and his co-workers developed the generalized valence bond (GVB) method.⁵ The sCI method, developed by Malrieu and co-workers, is known as configuration interaction using a perturbative selection done iteratively (CIPSI). Its performance has been further enhanced by improving convergence,⁶,⁷ accuracy,⁸ and efficiency.⁹,¹⁰ In the 1980s, building upon the CASSCF theory,¹¹ Roos and co-workers developed the constrained CAS and restricted active space SCF (RASSCF) methods.¹²,¹³ Following the year 2000, Olsen and Gagliardi proposed the generalized active space (GAS) method,¹⁴⁻¹⁶ and Ivanic introduced the occupation-restricted

<footnote n="1">¹ Institute of Frontier Chemistry, School of Chemistry and Chemical Engineering, Shandong University, Qingdao, Shandong 266237, China</footnote>

<footnote n="2">²Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mülheim an der Ruhr, Germany</footnote>

<footnote n="3">³Aix Marseille Univ., CNRS, Centrale Méditerranée, iSm2, Marseille, France</footnote>

<footnote>a) Author to whom correspondence should be addressed: frank.neese@kofo.mpg.de</footnote>
multiple-active-space (ORMAS) method.¹⁷ In recent years, sCI has seen a revival with contributions from various research groups worldwide, including those led by Head-Gordon et al.,¹⁸,¹⁹ Schriber and Evangelista,²⁰,²¹ Chilkuri and Neese,²²,²³ Zhang et al.,²⁴⁻²⁸ and Smith et al.,²⁹ among others.³⁰⁻³³ The fundamental concept underlying all these methods is utilizing a subset of the Hilbert space as an approximation to the complete CASCI space. In the recent decades, full configuration interaction quantum Monte Carlo (FCIQMC)³⁴⁻³⁷ and density matrix renormalization group (DMRG)³⁸⁻⁴⁵ theories have been developed as alternative approximate FCI methods as well. These approaches have been successfully applied to tackle static correlations in systems with approximately one hundred active electrons and MOs.

To account for the dynamic correlations of electrons on top of the approximate CASSCF/CASCI reference, a variety of multi-reference (MR) methods have been developed. At MRCI or MR coupled-cluster (CC) levels, Yanai and his team developed fully internally contracted (FIC) MRCISD and MR canonical transformation methods using the DMRG reference.⁴⁶,⁴⁷ Ma and co-workers reported an externally contracted (ec) MRCISD method with a DMRG reference as well.⁴⁸ Most recently, to further reduce the computational costs of ec-MRCISD, they introduced the renormalized residue technique.⁴⁹ Sharma also reported linearized coupled cluster and strongly contracted MRCI methods based on DMRG-SCF.⁵⁰,⁵¹ Among various dynamic correlation methods, MR perturbation theory (PT) stands out as the most computationally cost-effective. Up to now, the CASPT2 method, proposed by Roos and co-workers, remains the most widely used MRPT2 method.⁵² Based on RAS and selected active space reference wave functions, Celani and Werner introduced a CASPT2 method for large active spaces in 2000.⁵³ Gagliardi and colleagues also developed an alternative RASPT2 method.⁵⁴ In recent decades, CASPT2 methods utilizing DMRG wave functions as references have been independently reported by Yanai et al.,⁵⁵⁻⁵⁷ Wouters et al.,⁵⁸,⁵⁹ and Nakatani and Guo.⁶⁰ Besides CASPT2, several NEVPT2 methods to treat systems with large active spaces were developed as well. Research groups, including those led by Chan et al.,⁶¹⁻⁶³ Freitag et al.,⁶⁴ Sharma et al.,⁵¹,⁶⁵ and Anderson et al.,⁶⁶ have documented NEVPT2 algorithms employing DMRG or FCIQMC references for large active spaces. Xu and Li reported a NEVPT2 method based on the GVB wave function.⁶⁷,⁶⁸ In addition, Evangelista advanced the field with driven similarity renormalization group (DSRG) MRPT2 methods with the sCI reference.⁶⁹ Recently, Park has developed analytical gradients for DSRG-MRPT2 based on the sCI method developed by Head-Gordon and co-workers,¹⁹ adaptive sampling CI, as well.⁷⁰ Most recently, Malrieu and co-workers have introduced a new MRPT2 method for handling sCI references in large active spaces.⁷¹ Alternative MRPT2 approaches based on DMRG references have been explored by Ma and Pernal as well.⁴⁹,⁷²

Thanks to the advantageous features of FIC-NEVPT2,⁷³,⁷⁴ such as its free-from intruder states and unitary invariance, we have made significant efforts over the past decades to extend the capabilities of NEVPT2. In addition to the linear-scaling NEVPT2 for large molecules with moderate CAS spaces⁷⁵ as well as the explicitly correlated NEVPT2 (NEVPT2-F12),⁷⁶ we have developed a linear-scaling version of NEVPT2-F12.⁷⁷ The new method allows us to compute the basis set limit for large molecules at the NEVPT2 level.

To treat systems with large active spaces, in the first paper of this series (referred to as Paper I), we accessed various approximations to the reduced density matrices (RDMs) in FIC-NEVPT2.⁷⁸ However, none of these approximations could produce reliable results due to the presence of false intruder states. In the second paper (Paper II), we identified that intruder states in NEVPT2 calculations using approximate CASCI references primarily arise from the rank reduction trick.⁷⁹ While the rank reduction trick can reduce one order of RDMs to decrease computational costs with a CAS reference, it should not be applied to approximate CASSCF/CASCI references. By incorporating the missing residuals in the rank reduced NEVPT2 formalism, the full-rank NEVPT2 (FR-NEVPT2) method can yield stable results with approximate CAS references. Note that the early NEVPT2 implementations by Dyall⁸⁰ and Angeli et al.⁸¹ should be equivalent to FR-NEVPT2, where the rank reduction trick is not utilized. The FR-NEVPT2 approach requires fifth-order RDMs in the correlation energy calculations of $S_i$ and $S_a$ subspaces. To circumvent the storage bottleneck for fourth- and fifth-order RDMs in FR-NEVPT2, we defined new intermediates for NEVPT2 or FR-NEVPT2 algorithms.⁸² The FR-NEVPT2 algorithm utilizing these new intermediates should be able to treat systems with large active spaces based on the sCI reference.

In the present work, we first introduce the ICE-SCF method, which is based on the sCI method, ICE, developed by some of us.²²,²³ The FR-NEVPT2 algorithm built upon the ICE-SCF reference is then reported. These novel ICE-SCF and FR-NEVPT2 methods are applied to study systems that require larger active spaces, including iron-porphyrin complexes and polyacene molecules. The new algorithms successfully handle systems with an active space of CAS(34,34).

## II. THEORY

### A. ICE-SCF

To implement an ICE-SCF algorithm, one can substitute the CASCI/FCI methods with ICE²³ within any CASSCF orbital optimization drivers. However, due to the use of incomplete CI in the CAS space, the orbital optimization in ICE-SCF must consider not only the typical occupied-virtual (occ-vir), occupied-active (occ-act), and active-virtual (act-vir) rotation parameters but should also include active-active (act-act) rotations.¹⁹,²⁷ Since these rotations will become redundant as the approximate solution is approaching the exact FCI solution, inclusion of such rotations will introduce convergence problems. Hence, these rotations are, in keeping with other approximate FCI methods, neglected. Below, we will investigate the impact of this approximation on the final NEVPT2 results by comparing several ways to canonicalize the final ICE-SCF orbitals.

The ICE method consists of two major parameters, TGEN and TVAR, which control the accuracy and size of the CI space. In each ICE iteration, there are two major steps. Each parameter governs one step. Configurations at the nth step are generated from the singles and doubles (S+D) of selected configurations at the (n-1)th step using the TGEN parameter. However, not all configurations determined by TGEN at the nth step are included in the CI space. The Epstein-Nesbet perturbative energy of each configuration is calculated. If the energy contribution of a configuration exceeds TVAR, the configuration is selected into the variational CI
space. Otherwise, the configuration is treated perturbatively. Thus, the TGEN parameter is responsible for the “expansion” in ICE and is slightly more important than the TVAR parameter. Our previous study shows that by setting TVAR = TGEN × 10⁻⁷, the ICE could always produce stable results and rely solely on one parameter.²³ A test set comprising the FCI reference of 21 diatomic atoms, referred to as FCI21 test set, was used to access the performance of ICE. It has been observed that setting TGEN to 10⁻⁴ results in ICE results deviating from FCI results by less than 1.0 mEH.²³ When TGEN is sufficiently small, the optimization of MOs within the active space can be neglected. However, with large TGEN settings, the absolute deviations of ICE-SCF energies with respect to FCI results can be substantial. Nevertheless, as demonstrated in Sec. III, these errors, also referred to as reference relaxation energies, can be recovered by the subsequent FR-NEVPT2 calculations. Thus, in the present work, the active-active rotation parameters κ<sub>tu</sub> are excluded from the orbital optimization. The orbital rotation operators in ICE-SCF consist solely of κ<sub>ai</sub> (occ-vir), κ<sub>ti</sub> (occ-act), and κ<sub>at</sub> (act-vir), similar to those in CASSCF, expressed as

$$ \hat{\kappa} = \sum_{ai} \kappa_{ai} \hat{E}_{ai} + \sum_{ti} \kappa_{ti} \hat{E}_{ti} + \sum_{at} \kappa_{at} \hat{E}_{at}. \quad (1) $$

The SuperCIPT algorithm is utilized for orbital optimization in ICE-SCF.⁸³ In ORCA, three ICE drivers are available, based on configuration (CFG), configuration state function (CSF), and determinant (DET), leading to three distinct ICE-SCF algorithms. In this study, only the CFG-based ICE-SCF is employed. The default value for TGEN, which governs the accuracy of ICE-SCF results, is set to 10⁻⁴.

## B. FR-NEVPT2 with ICE-SCF reference

In Paper II, we demonstrated that FR-NEVPT2 consistently produces stable results with approximated CASCI reference.⁷⁹ To implement FR-NEVPT2, the residual term of each subspace (excluding the S<sub>ijab</sub> subspace) must be computed and added to the Koopmans matrix based on the reference ICE-SCF/ICE-CI wave function |Ψ<sup>ICE</sup>⟩. Taking the S<sub>a</sub> subspace as an example, the residual of the active space can be expressed as follows:

$$ \tilde{R}_{tuv}^{\overline{tuv}} = \langle \Psi^{\text{ICE}} | \tilde{E}_{\tilde{v}}^{\overline{u}} \tilde{E}_{\tilde{t}}^{\overline{t}} \tilde{E}_{\tilde{u}}^{\overline{v}} H_{\text{Dyall}}^{\text{Act}} | \Psi^{\text{ICE}} \rangle - (E_0 - C') \langle \Psi^{\text{ICE}} | \tilde{E}_{\tilde{v}}^{\overline{u}} \tilde{E}_{\tilde{t}}^{\overline{t}} \tilde{E}_{\tilde{u}}^{\overline{v}} | \Psi^{\text{ICE}} \rangle, \quad (2) $$

where the constant C' represents the ICE-SCF/ICE-CI energy contribution from the doubly occupied MOs, given by $\sum_i 2h_{ii} + \sum_{ij} 2(|ii|jj) - (ij|ij)$. As discussed in Paper II, calculating the first term in $\tilde{R}_{tuv}^{\overline{tuv}}$ involves in the fifth-order RDMs. However, by using properly designed intermediates, the need to store fourth- and fifth-order RDMs can be avoided.⁸² The computational costs as well as the storage of all matrices involved in the FR-NEVPT2 are similar to those of third-order RDMs, with a larger prefactor. In FR-NEVPT2, the residuals of seven different subspaces can be derived from three matrices,

$$ R_1 \tilde{t} = \langle \Psi^{\text{ICE}} | \tilde{E}_t^I H_{\text{Dyall}}^{\text{Act}} | \Psi^{\text{ICE}} \rangle, \quad (3) $$

$$ R_2 \tilde{t}_{tu} = \langle \Psi^{\text{ICE}} | \tilde{E}_t^I \tilde{E}_u^I H_{\text{Dyall}}^{\text{Act}} | \Psi^{\text{ICE}} \rangle, \quad (4) $$

$$ R_3 \tilde{t}_{tuv} = \langle \Psi^{\text{ICE}} | \tilde{E}_{\tilde{v}}^{\overline{u}} \tilde{E}_{\tilde{t}}^{\overline{t}} \tilde{E}_{\tilde{u}}^{\overline{v}} H_{\text{Dyall}}^{\text{Act}} | \Psi^{\text{ICE}} \rangle. \quad (5) $$

Assuming that the three matrices based on ICE-SCF/ICE-CI wave functions are computed, the FR-NEVPT2 algorithm can be implemented straightforwardly on top of the conventional FIC-NEVPT2 algorithm. For the calculation of intermediates R1 to R3, the resolution of identity (RI) CSF space |Ψ<sup>RI</sup>⟩ is also employed as in the evaluation of RDMs.⁷⁸ The intermediate R3 can be computed as

$$ \begin{aligned} R_3 \tilde{t}_{tuv} &= \sum_J \langle \Psi^{\text{ICE}} | \tilde{E}_{\tilde{v}}^{\overline{u}} \tilde{E}_{\tilde{t}}^{\overline{t}} \tilde{E}_{\tilde{u}}^{\overline{v}} | \Psi_J^{\text{RI}} \rangle \langle \Psi_J^{\text{RI}} | H_{\text{Dyall}}^{\text{Act}} | \Psi^{\text{ICE}} \rangle \\ &= \sum_J \langle \Psi^{\text{ICE}} | \tilde{E}_{\tilde{v}}^{\overline{u}} \tilde{E}_{\tilde{t}}^{\overline{t}} \tilde{E}_{\tilde{u}}^{\overline{v}} | \Psi_J^{\text{RI}} \rangle Z_J, \quad (Z_J = \langle \Psi_J^{\text{RI}} | H_{\text{Dyall}}^{\text{Act}} | \Psi^{\text{ICE}} \rangle), \end{aligned} \quad (6) $$

where J runs over all CSFs in the |Ψ<sup>RI</sup>⟩ space. It should be noted that the dimension of the |Ψ<sup>RI</sup>⟩ space depends not only on the number of CSFs in |Ψ<sup>ICE</sup>⟩ but also on the number of active electrons and active MOs.

In the present FR-NEVPT2 algorithm, only RDMs and Koopmans matrices with six or fewer active indices (such as the 3rd-order RDM) are stored in random access memory (RAM). Although the fourth- and fifth-order RDMs are not explicitly computed and stored, intermediates defined in our previous study,⁸² such as

$$ X_{tu}^{I} = \sum_{vwX} \langle \Psi^{RI} | \hat{E}_{t}^{v} \hat{E}_{x}^{w} | \Psi^{ICE} \rangle (vw|ux), \quad (7) $$

could be bottlenecks as well. In the calculation of octacene triplet state, with CAS(34,34) and TVAR = TGEN × 10⁻³ = 10⁻⁶ settings, the ICE-SCF reference consists of 34 736 CSFs, whereas the |Ψ<sup>RI</sup>⟩ space consists of 868 627 838 CSFs. Thus, intermediates such as $X_{tu}^I$ could require up to 7.3 TB of memory, which may not be possible to store in RAM. To address the storage bottleneck, the present FR-NEVPT2 algorithm evaluates these intermediates on-the-fly during the computation of the Koopmans and R1 to R3 matrices. Only vectors associated with the |Ψ<sup>RI</sup>⟩ space, such as Z<sub>J</sub> in Eq. (6), are stored in RAM. In the calculation of octacene, the memory needed to store Z<sub>J</sub> (6.5 GB) is less than that required for storing third-order RDMs (11.5 GB).

After getting the ICE-SCF reference wave functions, the subsequent FR-NEVPT2 calculations could be performed without any truncation parameters. However, the computational costs of FR-NEVPT2 are significantly higher than those of ICE-SCF. Typically, the reference wave function generated by ICE-SCF is too large to be directly used in FR-NEVPT2. Therefore, the ICE-SCF reference must be further truncated before proceeding with FR-NEVPT2 calculations to reduce the computational costs. In the present work, this is achieved by truncating the ICE-SCF wave function utilizing the same TGEN parameter but a more aggressive TVAR parameter,

$$ |\Psi^{\text{ICE}}\rangle = |\bar{\Psi}^{\text{ICE}}\rangle + |\tilde{\Psi}^{\text{ICE}}\rangle = \sum_{\bar{I}} T_{\bar{I}} |\bar{I}\rangle + \sum_{\bar{j}} T_{\bar{j}} |\bar{j}\rangle (T_{\bar{I}}^2 > \text{TVAR}, T_{\bar{j}}^2 \le \text{TVAR}). \quad (8) $$

The truncation of |Ψ<sup>ICE</sup>⟩ is performed in the same way as that of |Ψ<sup>CAS</sup>⟩ using the prescreening (PS) approximation reported in
Papers I and II.⁷⁸,⁷⁹ Note that when applying the PS approximation to truncate the ICE-SCF wave function, $E_0$ in Eq. (2) should be replaced by the new energy

$$ \bar{E}_0 = \langle \bar{\Psi}^{\text{ICE}} | H | \bar{\Psi}^{\text{ICE}} \rangle / \langle \bar{\Psi}^{\text{ICE}} | \bar{\Psi}^{\text{ICE}} \rangle. \quad (9) $$

This is because that the final residuals, as well as various RDMs and Koopmans matrices, are evaluated with the $|\bar{\Psi}^{\text{ICE}}\rangle$ wave functions. For example, the residual of the $S_a$ subspace is computed as

$$ \bar{R}_{tuv}^{\bar{\psi} \bar{\psi} \bar{\psi}} = \langle \bar{\Psi}^{\text{ICE}} | \hat{E}_v^{\bar{\psi}} \hat{E}_t^{\bar{\psi}} \hat{E}_u^{\bar{\psi}} H_{\text{Dyall}}^{\text{Act}} | \bar{\Psi}^{\text{ICE}} \rangle - (\bar{E}_0 - C') \langle \bar{\Psi}^{\text{ICE}} | \hat{E}_v^{\bar{\psi}} \hat{E}_t^{\bar{\psi}} \hat{E}_u^{\bar{\psi}} | \bar{\Psi}^{\text{ICE}} \rangle. \quad (10) $$

### III. RESULTS

The new FR-NEVPT2 algorithm is implemented in a development version based on ORCA 6.0. The CFG-based state-specific ICE-SCF, which is available in the released ORCA 6.0, is employed to compute the reference wave function with large active spaces. In this work, MP2 natural MOs (NMOs) are used as the initial guess for the ICE-SCF calculations. The setting TVAR = TGEN × 10⁻⁷ was used to generate the ICE-SCF MOs. The cc-pVXZ (X = D, T, Q) basis sets, along with their corresponding auxiliary basis sets, are used and are referred to as DZ, TZ, and QZ, respectively. The geometries of the n-acenes (n = 3–7) singlet and triplet states are optimized at the UB3LYP/6-31G(d) level by Evangelista and colleagues,⁶⁴ whereas the Cartesian coordinates of octacene are obtained via ASCI-SCF by Park.³³ For iron porphyrin, the geometries reported by Pierloot and colleagues are adopted. The geometries of all three low-lying states are optimized using PBE0/TZVP.⁸⁴

#### A. ICE-SCF

While the accuracy of ICE has been verified and employed extensively in previous research,²³ the performance of the ICE-SCF

method remains to be thoroughly evaluated. In this study, we investigated the singlet and triplet states of anthracene and tetracene using the ICE-SCF method with double-zeta (DZ) basis sets. The ICE-SCF results, calculated with varying TGEN values, are presented in Table I. The CASSCF results serve as the reference. For anthracene, with a TGEN of 10⁻⁴, the absolute deviations of ICE-SCF are less than 1.1 kcal/mol for both singlet and triplet states, yielding a singlet–triplet (S–T) gap that is as accurate as that obtained from CASSCF. By contrast, for tetracene, the absolute errors of ICE-SCF compared to the CASSCF results exceed 6.0 kcal/mol for both states with TGEN = 5 × 10⁻⁴. However, the adiabatic S–T gaps are accurately predicted by ICE-SCF, with a deviation of only 0.75 kcal/mol. By tightening the TGEN threshold to 10⁻⁵, the absolute errors for both states of tetracene decrease to less than 1.3 kcal/mol. Conversely, using a crude TGEN parameter of 10⁻² can result in absolute deviations exceeding 32.3 kcal/mol, leading to overestimations of the adiabatic S–T gaps by 4.9 and 7.3 kcal/mol for anthracene and tetracene, respectively. Notably, the molecular orbitals (MOs) optimized by ICE-SCF with loose TGEN parameters are as accurate as those obtained with tight TGEN settings. The ICE-CI results, calculated with MOs derived from various TGEN parameters, are shown in Table S1 of the supplementary material. The ICE-CI results with the same TGEN indicate minimal discrepancies, less than 1.0 kcal/mol, across three sets of MOs optimized with different TGEN settings. This suggests that the ICE-SCF MOs optimized with coarse thresholds are nearly indistinguishable from those optimized with tight TGEN settings. Nonetheless, to ensure the reliability of relative energies, we recommend using a TGEN setting no coarser than 10⁻³, provided that such calculations are feasible.

#### B. FR-NEVPT2 with ICE-SCF reference

As discussed in Sec. II B, the computational demands of FR-NEVPT2 are significantly greater than those of ICE-SCF, and the parameters used in ICE-SCF may not be directly used in subsequent FR-NEVPT2 calculations. Therefore, to ensure computational

<table><caption>TABLE I. ICE-SCF results for the singlet and triplet states of anthracene and tetracene using various TGEN parameters and DZ basis sets. The ICE-SCF energies relative to the CASSCF energies (in hartree) for singlet and triplet states of polyacenes, along with the errors in adiabatic S–T gaps (in kcal/mol) produced by ICE-SCF compared to those by CASSCF, are given. The active spaces CAS(14,14) and CAS(18,18) are employed for anthracene and tetracene, respectively.</caption><thead><tr><th rowspan="3">TGEN</th><th colspan="5">Anthracene</th><th colspan="3">Tetracene</th></tr><tr><th colspan="2">Singlet</th><th colspan="3">Triplet</th><th colspan="2">Singlet</th><th colspan="2">Triplet</th></tr><tr><th>No. CSFs</th><th>ICE-SCF<br>(hartree)</th><th>No. CSFs</th><th>ICE-SCF<br>(hartree)</th><th>S–T gap<br>(kcal/mol)</th><th>No. CSFs</th><th>ICE-SCF<br>(hartree)</th><th>No. CSFs</th><th>ICE-SCF<br>(hartree)</th><th>S–T gap<br>(kcal/mol)</th></tr></thead><tbody><tr><td>1 × 10<sup>-2</sup></td><td>1882</td><td>0.022 08</td><td>2610</td><td>0.029 93</td><td>4.93</td><td>5560</td><td>0.039 95</td><td>5604</td><td>0.051 56</td><td>7.28</td></tr><tr><td>1 × 10<sup>-3</sup></td><td>19 500</td><td>0.004 48</td><td>58 131</td><td>0.004 93</td><td>0.28</td><td>94 314</td><td>0.013 92</td><td>304 308</td><td>0.015 57</td><td>1.04</td></tr><tr><td>5 × 10<sup>-4</sup></td><td>23 697</td><td>0.003 55</td><td>76 004</td><td>0.003 65</td><td>0.06</td><td>142 259</td><td>0.009 73</td><td>476 529</td><td>0.010 93</td><td>0.75</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>73 273</td><td>0.001 36</td><td>174 986</td><td>0.001 69</td><td>0.21</td><td>437 167</td><td>0.005 05</td><td>1 123 919</td><td>0.006 26</td><td>0.76</td></tr><tr><td>5 × 10<sup>-5</sup></td><td>100 118</td><td>0.000 99</td><td>251 282</td><td>0.001 11</td><td>0.08</td><td>755 784</td><td>0.003 96</td><td>1 875 736</td><td>0.005 00</td><td>0.65</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>208 235</td><td>0.000 28</td><td>470 732</td><td>0.000 29</td><td>0.01</td><td>2 305 880</td><td>0.001 80</td><td>6 401 907</td><td>0.002 04</td><td>0.15</td></tr><tr><td>CASSCF</td><td>691 160</td><td>-536.212 22</td><td>1 252 895</td><td>-536.140 79</td><td>44.82</td><td>112 298 248</td><td>-688.915 71</td><td>225 503 748</td><td>-688.863 16</td><td>32.98</td></tr></tbody></table>
feasibility, it is necessary to employ somewhat more aggressive truncation parameters when computing the reference of FR-NEVPT2. In the ICE-SCF algorithm, the TVAR parameter serves as the truncation threshold for the PS approximation, as reported in Paper II.⁷⁹ In the present work, subsequent to the convergence of the ICE-SCF calculation with TVAR = TGEN × 10⁻⁷, an ICE-CI calculation is performed using the identical TGEN but with a larger TVAR. This additional ICE-CI step is intended to refine the ICE-SCF reference for FR-NEVPT2 applications. To assess the influence of TVAR,

FR-NEVPT2 calculations are carried out using different sets of ICE-SCF MOs and various TVAR parameters.

The adiabatic singlet–triplet (S–T) gaps of anthracene and tetracene are calculated using FR-NEVPT2 with ICE-SCF MOs generated with TGEN values of 1 × 10⁻³ and 1 × 10⁻⁴. Based on these two sets of MOs, the ICE-CI and FR-NEVPT2 results for anthracene with varying TVAR parameters are presented in Table II. The results indicate that by utilizing TVAR ≤1 × 10⁻⁶ settings, the absolute deviations of FR-NEVPT2 with respect to the canonical

<table><caption>**TABLE II.** ICE-CI and FR-NEVPT2 results of anthracene singlet and triplet states with different TGEN and TVAR parameters. The active space CAS(14,14) is employed with DZ basis sets. Besides the deviations with respect to the canonical CASSCF and NEVPT2 energies (in hartree), the errors of adiabatic S–T gaps (in kcal/mol) with respect to those of CASSCF and NEVPT2 and the number of CSFs in the reference (No. CSF) and the RI space (No. RI_CSF) are given as well.</caption><thead><tr><th rowspan="2">TGEN</th><th rowspan="2">TVAR</th><th colspan="3">Singlet</th><th colspan="3">Triplet</th><th colspan="2">S-T gap (Kcal/mol)</th></tr><tr><th>No. CSF</th><th>No. RI_CSF</th><th>ICE-CI (hartree)</th><th>FR-NEVPT2 (hartree)</th><th>No. CSF</th><th>No. RI_CSF</th><th>ICE-CI (hartree)</th><th>FR-NEVPT2 (hartree)</th><th>ICE-CI</th><th>FR-NEVPT2</th></tr></thead><tbody><tr><td rowspan="4">1 × 10<sup>-3</sup></td><td>1 × 10<sup>-5</sup></td><td>705</td><td>356 841</td><td>0.023 24</td><td>-0.002 49</td><td>1510</td><td>653 076</td><td>0.030 93</td><td>-0.003 76</td><td>4.82</td><td>-0.80</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>4360</td><td>832 102</td><td>0.012 90</td><td>-0.001 67</td><td>9212</td><td>1 798 382</td><td>0.016 74</td><td>-0.002 46</td><td>2.41</td><td>-0.50</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>12 552</td><td>981 061</td><td>0.006 37</td><td>-0.001 74</td><td>29 050</td><td>2 374 613</td><td>0.008 44</td><td>-0.002 39</td><td>1.30</td><td>-0.41</td></tr><tr><td>1 × 10<sup>-11</sup></td><td>19 500</td><td>990 305</td><td>0.004 48</td><td>-0.001 70</td><td>58 131</td><td>2 463 189</td><td>0.004 93</td><td>-0.002 13</td><td>0.28</td><td>-0.27</td></tr><tr><td rowspan="4">1 × 10<sup>-4</sup></td><td>1 × 10<sup>-5</sup></td><td>783</td><td>396 811</td><td>0.022 05</td><td>-0.002 37</td><td>1848</td><td>810 952</td><td>0.029 07</td><td>-0.003 56</td><td>4.41</td><td>-0.74</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>6623</td><td>956 591</td><td>0.009 63</td><td>-0.002 08</td><td>13 375</td><td>2 049 559</td><td>0.013 38</td><td>-0.002 66</td><td>2.35</td><td>-0.36</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>20 584</td><td>1 260 487</td><td>0.003 30</td><td>-0.001 81</td><td>51 415</td><td>2 680 360</td><td>0.004 62</td><td>-0.002 29</td><td>0.83</td><td>-0.30</td></tr><tr><td>1 × 10<sup>-11</sup></td><td>73 273</td><td>2 040 190</td><td>0.001 36</td><td>-0.000 84</td><td>174 986</td><td>4 050 638</td><td>0.001 69</td><td>-0.001 08</td><td>0.21</td><td>-0.15</td></tr><tr><td colspan="2">0.0</td><td>691 160</td><td>2 760 615</td><td>-536.212 22</td><td>-537.849 98</td><td>1 252 895</td><td>5 010 005</td><td>-536.140 79</td><td>-537.776 25</td><td>44.82</td><td>46.26</td></tr></tbody></table>

<table><caption>**TABLE III.** ICE-CI and FR-NEVPT2 results of tetracene singlet and triplet states with different TGEN and TVAR parameters. The active space CAS(18,18) is employed with DZ basis sets. Besides the deviations with respect to the most accurate ICE-CI and FR-NEVPT2 energies (in hartree), the errors of adiabatic S–T gaps (in kcal/mol) and the number of CSFs in the reference (No. CSF) and the RI space (No. RI_CSF) are given as well.</caption><thead><tr><th rowspan="2">TGEN</th><th rowspan="2">TVAR</th><th colspan="3">Singlet</th><th colspan="3">Triplet</th><th colspan="2">S-T gap (Kcal/mol)</th></tr><tr><th>No. CSF</th><th>No. RI_CSF</th><th>ICE-CI (hartree)</th><th>FR-NEVPT2 (hartree)</th><th>No. CSF</th><th>No. RI_CSF</th><th>ICE-CI (hartree)</th><th>FR-NEVPT2 (hartree)</th><th>ICE-CI</th><th>FR-NEVPT2</th></tr></thead><tbody><tr><td rowspan="4">1 × 10<sup>-3</sup></td><td>1 × 10<sup>-5</sup></td><td>1084</td><td>1 843 314</td><td>0.041 97</td><td>-0.001 03</td><td>2915</td><td>3 909 399</td><td>0.046 77</td><td>0.000 17</td><td>3.01</td><td>0.75</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>6036</td><td>9 153 689</td><td>0.029 72</td><td>-0.000 05</td><td>12 836</td><td>20 114 588</td><td>0.032 46</td><td>0.001 09</td><td>1.72</td><td>0.71</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>33 906</td><td>19 412 354</td><td>0.017 52</td><td>-0.000 29</td><td>70 925</td><td>48 728 138</td><td>0.017 26</td><td>0.000 95</td><td>-0.17</td><td>0.78</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>80 124</td><td>22 167 302</td><td>0.010 86</td><td>-0.000 36</td><td>223 359</td><td>68 867 148</td><td>0.008 79</td><td>0.000 85</td><td>-1.30</td><td>0.76</td></tr><tr><td rowspan="4">1 × 10<sup>-4</sup></td><td>1 × 10<sup>-5</sup></td><td>1313</td><td>2 569 830</td><td>0.040 08</td><td>-0.001 48</td><td>3 128</td><td>4 793 844</td><td>0.045 36</td><td>-0.000 27</td><td>3.31</td><td>0.76</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>10 953</td><td>13 261 163</td><td>0.025 07</td><td>-0.000 56</td><td>20 191</td><td>27 747 513</td><td>0.028 14</td><td>0.000 49</td><td>1.93</td><td>0.66</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>76 453</td><td>25 913 513</td><td>0.009 04</td><td>-0.001 56</td><td>165 262</td><td>73 008 727</td><td>0.009 33</td><td>0.000 01</td><td>0.18</td><td>0.99</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>170 456</td><td>37 872 264</td><td>0.002 85</td><td>-0.001 50</td><td>505 827</td><td>90 307 883</td><td>0.000 28</td><td>-0.000 05</td><td>-1.61</td><td>0.91</td></tr><tr><td colspan="2">Reference<br>(TGEN = 5 × 10<sup>-5</sup>)</td><td>755 784</td><td>97 954 321</td><td>-688.911 75<sup>a</sup></td><td>-691.017 29<sup>a</sup></td><td>553 099</td><td>99 305 678</td><td>-688.854 77<sup>b</sup></td><td>-690.968 14<sup>b</sup></td><td>...</td><td>30.84</td></tr></tbody></table>

<footnote><sup>a</sup>For the singlet state, the ICE-CI with TVAR = 10⁻¹¹ is computed as the reference for FR-NEVPT2.</footnote>

<footnote><sup>b</sup>For the triplet state, the ICE-CI with TVAR = 10⁻⁸ is computed as the reference for FR-NEVPT2.</footnote>
NEVPT2 results for singlet and triplet states are less than 2.0 mEh and 3.0 mEh, respectively, even though the reference ICE-CI energies may differ from the CASSCF energies by more than 16.0 mEh. For tetracene, the most precise FR-NEVPT2 results, as seen in Table III, are achieved with TVAR settings of $10^{-11}$ for singlets and $10^{-8}$ for triplets. Since canonical NEVPT2 results are not available, these numbers will serve as the reference values in the following comparisons. The results in Table III reveal that using a crude parameter ($TVAR = 10^{-5}$) can result in absolute errors of ICE-CI energies exceeding 25.0 kcal/mol compared to the reference. Even with the smallest TGEN and TVAR settings, the absolute deviations of ICE-CI for both singlet and triplet states are above 3.0 kcal/mol when compared to the CASSCF results in Table I. To achieve converged S–T gaps at the ICE-CI level, TVAR must be set to $10^{-7}$. Despite the significant deviations in the ICE-CI results with crude TVAR settings, FR-NEVPT2 consistently provides accurate absolute energies. This suggests that FR-NEVPT2 can compensate for deficiencies in the reference wave functions.

These findings align with results reported in Paper II, underscoring the robustness of the FR-NEVPT2 method. The data in Tables II and III demonstrate that setting TGEN to $1 \times 10^{-4}$ and TVAR to $1 \times 10^{-7}$ allows the FR-NEVPT2 method to deliver both accurate absolute and relative energies for anthracene and tetracene. However, for calculations involving more than 30 active electrons and MOs, the settings become unfeasible due to the size of the $|\Psi^{RI}\rangle$ space. In such cases, looser TGEN and TVAR parameters must be employed. To avoid adjusting two parameters, the remainder of this work will use the setting $TVAR = TGEN \times 10^{-3}$ in all NEVPT2 calculations. Thus, the size of the ICE-CI reference and the accuracy of FR-NEVPT2 depend solely on the TGEN parameter. Unless otherwise specified, the TGEN is set to $1 \times 10^{-4}$ for all FR-NEVPT2 calculations.

To further understand the reference relaxation effect of FR-NEVPT2, a detailed analysis of energy contributions for both anthracene and tetracene molecules is conducted using different TGEN parameters. The errors in energy contributions from ICE-CI

<figure><img src="image_10.png" /><figcaption>FIG. 1. FR-NEVPT2 energy contributions of singlet and triplet states of (a) anthracene and (b) tetracene employing different TGEN parameters with respect to the reference results in Tables II and III.</figcaption></figure>
and eight first order interaction subspaces relative to the anthracene reference results are presented in Fig. 1(a). For the FR-NEVPT2 results with TGEN = 1 × 10⁻³, the poor ICE-CI energies are significantly relaxed by the Sᵢₐ subspace, contributing 7.3 and 9.2 kcal/mol for the singlet and triplet states, respectively. In addition, the Sᵢₐ and Sᵢₐb subspaces contribute ~1.0 kcal/mol each. Theoretically, the Sᵢₐ, Sᵢ, and Sᵢₐ subspaces, which involve spectator excitations, should exhibit more pronounced reference relaxation effects. However, for planar polyacene molecules, the correlation energies of the Sᵢ subspace are consistently zero due to spatial symmetry, resulting in no contribution from the Sᵢ subspace. With more accurate ICE-CI references, setting TGEN to 1 × 10⁻⁴ and 1 × 10⁻⁵, the Sᵢₐ subspace continues to be the primary contributor in recovering the reference ICE-CI results. This could be attributed to the wave function of the Sᵢₐ subspace, which comprises single excitations within the σ MO manifold coupled with single excitations within the π MO manifold,⁸⁵ thereby contributing significantly to

the absolute and relative correlation energies. The absolute FR-NEVPT2 results predicted with the three different TGEN parameters are nearly identical, with the largest error being less than 1.5 kcal/mol. By using TGEN = 1 × 10⁻⁶ settings, all contributions of FR-NEVPT2 converge to the canonical results. Similar results were observed for tetracene as well. Although the deviations in the ICE-CI reference exceed 18.0 kcal/mol for both states, the FR-NEVPT2 results with crude TVAR parameters are nearly identical to those obtained using the most accurate parameters. In fact, we performed FR-NEVPT2 calculations with even more lenient settings, such as TGEN = 0.01, where the ICE-CI reference for the singlet state of tetracene comprises only 892 CSFs. However, the deviation of the absolute FR-NEVPT2 result is less than 1.3 kcal/mol compared to the reference in Table III.

For comparison, NEVPT2 calculations are also performed using the same settings as those in Fig. 1. In these NEVPT2 calculations, the residuals of the Koopmans matrices are absent.

<figure><img src="image_11.png" /><figcaption>**FIG. 2.** Rank reduced NEVPT2 energy contributions of singlet and triplet states of (a) anthracene and (b) tetracene employing different TGEN parameters with respect to the reference results in Tables II and III.</figcaption></figure>
Consequently, with crude TGEN parameters, intruder states are encountered, leading to poor absolute and relative results (Fig. 2). When TGEN is set to $1 \times 10^{-3}$, the singlet state energy of anthracene is significantly underestimated, which is primarily attributed to the emergence of false intruder states. For both the singlet and triplet states of anthracene, there are deviations exceeding 4.0 kcal/mol from the canonical NEVPT2 values when TGEN is set to $1 \times 10^{-5}$. The NEVPT2 energies only converge to the reference values when $TGEN \le 1 \times 10^{-6}$ settings are used. Similarly, for tetracene, the NEVPT2 performance is unsatisfactory. The absolute deviations for the singlet state exceed 40 kcal/mol even with the most accurate settings. As highlighted in Paper I, conventional NEVPT2 with an approximate CAS reference, such as ICE-CI, fails to yield reliable results. The residuals in FR-NEVPT2 are crucial not only for generating accurate correlation energies but also for correcting the reference energies, thereby ensuring the precision of the absolute energies.

### C. The influence of active MOs

In most DMRG-SCF methods, the active-active rotation para-meters are not considered in orbital optimization, utilizing localized MOs (LMOs) instead. Typically, the LMOs yield lower absolute energies compared to canonical MOs (CMOs) or NMOs. It has been observed that in sCI-based MCSCF methods, active-active rotations are crucial for obtaining stable relative energies.²⁷ However, as mentioned above, considering orbital rotations within the active space makes the convergence of sCI-based MCSCF methods more challenging. In Sec. III B, we demonstrated that errors at the ICE-SCF level can be effectively recovered by FR-NEVPT2. Consequently, the choice of active MOs should not impact the FR-NEVPT2 results. Nonetheless, in this subsection, we numerically assess the influence of active MOs on FR-NEVPT2.

For the optimization of ICE-SCF orbitals, we utilize the TVAR = $TGEN * 10^{-7} = 1 \times 10^{-11}$ settings (i.e., default settings), with MP2 natural MOs as the initial guess. Employing the converged

<figure><img src="image_12.png" /><figcaption>FIG. 3. FR-NEVPT2 energy contributions of tetracene (a) singlet and (b) triplet states utilizing different types of active MOs with $TGEN = 1 \times 10^{-4}$, compared with the reference results in Table III.</figcaption></figure>
ICE-SCF MOs, FR-NEVPT2 calculations, with TVAR = TGEN
* 10⁻³ = 1 × 10⁻⁷, are conducted using three different types of active
MOs: CMOs, NMOs, and LMOs. The active CMOs are derived
by diagonalizing the generalized Fock matrix of the active space,
whereas NMOs are obtained by diagonalizing the first-order RDMs.
The active LMOs are generated by separately localizing nearly dou-
bly occupied MOs and nearly unoccupied MOs to circumvent dif-
ficulties in generating initial guesses for ICE-CI calculations. The
FR-NEVPT2 contributions using different active MOs are presented
in Fig. 3, with the most accurate FR-NEVPT2 results in Table III
serving as the reference. For the singlet state, the ICE-CI results uti-
lizing LMOs are slightly lower. However, the overall FR-NEVPT2
energies calculated with the three types of active MOs are nearly
identical, with the maximum absolute error in total energy being
less than 1.1 kcal/mol. For the triplet state of tetracene, the ICE-
CI results with LMOs are 1.5 kcal/mol higher than those obtained
with CMOs. Nevertheless, the final FR-NEVPT2 results using differ-
ent active MOs are also close to each other. Therefore, the S–T gaps
computed with different types of active MOs are similar. Regardless
of the type of active MOs used, the FR-NEVPT2 method is capa-
ble of recovering the remaining static correlations within the active
space. Therefore, provided that FR-NEVPT2 calculations are per-
formed, the optimization of orbitals within the active space may be
neglected.

D. Iron porphyrin (FeP)

The spin-state energy order of the iron porphyrin (FeP) com-
plex serves as a critical benchmark for MR methods and attracts
the interest of many theoretical groups.<sup>86–88</sup> Over recent decades,
Pierloot and colleagues have extensively studied FeP using various
MR methods.<sup>84,89–91</sup> They utilized CASPT2 and NEVPT2 methods
with CAS(8,11) active spaces, which encompass the 3d4d orbitals of
iron and one σ-ligand MO, to investigate the spin states of FeP.<sup>84</sup>
In a subsequent CASPT2/CC study, they highlighted the signifi-
cant role of semi-core correlation in the spin state energetics of
transition metal complexes.<sup>92</sup> By incorporating 3s3p correlations at
the CCSD(T) level with CASPT2, the CASPT2/CC approach was
able to predict <sup>3</sup>A<sub>2g</sub> as the ground state. Li Manni and co-workers
suggested that, compared to the 3s3p orbitals, the π-orbitals play
a more important role in stabilizing the triplet states of FeP.<sup>86,93</sup>
However, to correctly describe the low-lying triplet and quintet
states, σ-ligand orbitals and the entire valence π orbitals of the por-
phyrin have to be included simultaneously in the active space in MR
calculations.<sup>94</sup> Recently, one of us designed various active spaces
to examine the energy differences among the <sup>3</sup>A<sub>2g</sub>, <sup>3</sup>E<sub>g</sub>, and <sup>5</sup>A<sub>1g</sub>
states using the iCISCF method.<sup>27</sup> It was found that, at the MCSCF
level, the choice of active space significantly influences the spin
state ordering. The active space CAS(40,42), including all π-orbitals
of the porphyrin, the 3d4s4p4d orbitals of iron, and four σ-ligand
MOs, could produce qualitatively correct energy order. This active
space CAS(40,42), also favored by Head-Gordon and co-workers,<sup>14</sup>
should be an ideal active space for both MCSCF and MRPT2 cal-
culations. The electronic structures of catalytically active species
[FeP]<sup>2−</sup> and FeP were investigated experimentally and theoretically
by our group as well.<sup>88,95</sup> The spectroscopic results indicated that
the ground state of the planar FeP complex is <sup>3</sup>A<sub>2g</sub>. However, this
state is nearly degenerate with a <sup>3</sup>E<sub>g</sub> state, and there is significant

scrambling due to spin-orbit coupling (SOC).<sup>88</sup> Thus, the ground
state of FeP is best interpreted as a mixture of three nearly degener-
ate triplet states.<sup>88</sup> The experiments did not find signs of low-lying
quintet states.

In this work, the adiabatic energy gaps between the three states,
³A₂g, ³Eg, and ⁵A₁g, are computed with ICE-SCF and FR-NEVPT2.
For the calculations with less than 14 active MOs, canonical CASSCF
and NEVPT2 calculations are performed. In calculations with active
MOs ranging from 15 to 20, the TGEN = 1 × 10⁻⁴ settings are
used in ICE-SCF and FR-NEVPT2 calculations. The TGEN = 5
× 10⁻⁴ settings are used in the FR-NEVPT2 calculations involv-
ing more than 20 active MOs, whereas TGEN = 10⁻³ settings are
used for the calculations with the active space CAS(30,29). To
check the influence of the basis set, DZ, TZ, and QZ basis sets are
employed. Calculations with the QZ basis set consist of 1784 basis
functions.

To validate the importance of the 3s3p orbitals, three differ-
ent active spaces are employed with the active space CAS(8,11) as
the starting point. The inclusion of 3s3p orbitals of iron results in
CAS(16,15). Further consideration of the “double sp shell” effect
leads to the active space CAS(16,19). The 19 active MOs involved
in these calculations are depicted in Fig. 4. The energy gaps rela-
tive to the ³A₂g state are presented in Table IV. For the NEVPT2
results based on the CAS(8,11) reference with QZ basis sets, the
energy gaps of ³Eg⁻³A₂g and ⁵A₁g⁻³A₂g are 5.2 and −5.3 kcal/mol,
respectively. The ⁵A₁g state is predicted to be the ground state, which
is consistent with the results reported by Pierloot and colleagues.⁸⁴
Although the inclusion of the 3s3p and 4s4p orbitals in the active
space decreases the ⁵A₁g⁻³A₂g gap at the SCF level, the ⁵A₁g⁻³A₂g
gaps predicted by FR-NEVPT2 are worse than those with CAS(8,11).
The inclusion of 3s3p4s4p orbitals in FR-NEVPT2 calculations stabi-
lizes the quintet states more than the triplet states. Thus, at NEVPT2
levels, to include 3s3p4s4p orbitals to active space may not be rec-
ommended for the prediction of the ⁵A₁g⁻³A₂g gap at the NEVPT2
level. In contrast, by including the 3s3p4s4p orbitals of iron, the
two triplet states ³A₂g and ³Eg become nearly degenerate. The
energy gaps between ³Eg and ³A₂g are less than 0.3 kcal/mol in all
cases, which is consistent with experimental and theoretical findings
reported in our work.⁸⁸ Comparing results with different basis sets,
the energy gaps by FR-NEVPT2 are overestimated with DZ basis
sets. Thus, for FeP, the results computed with a double-ξ basis set
are not reliable.

In the CASSCF results with CAS(8,11), the two doubly occu-
pied orbitals with the highest orbital energies and the two virtual
orbitals with the lowest energies are the four “Gouterman orbitals.”<sup>96</sup>
In order to test the influence of “Gouterman orbitals,” we have
undertaken calculations with an active space of CAS(12,15) (Fig. 5).
As might have been anticipated, inclusion of these orbitals has small
effects at the ICE-SCF level, where the quintet state remains the
ground state by 17–20 kcal/mol (Table V). However, the effect of
including “Gouterman orbitals” in the dynamic correlation calcula-
tion is rather profound and preferentially stabilizes the triplet states
to the extent that they become near-degenerate with the quintet
state. This may be a result of achieving a better balance in the lig-
and field since the degenerate E<sub>g</sub> orbitals are antibonding between
the ligand and the iron d<sub>xy</sub> and d<sub>yz</sub> orbitals. Further including
six additional π-orbitals (Fig. 5), both ICE-SCF and FR-NEVPT2
results with the CAS(20,21) active space are similar to those with
<figure><img src="image_13.png" /><figcaption>FIG. 4. Active MOs used in the calculations of FeP.</figcaption></figure>

<table><caption>TABLE IV. Relative energies of $^3$E$_g$ and $^5$A$_{1g}$ states with respect to the $^3$A$_{2g}$ state of FeP, employing different active spaces in Fig. 4 with different basis sets.</caption><thead><tr><th rowspan="2">Active space</th><th rowspan="2">Method</th><th colspan="3"><sup>3</sup>E<sub>g</sub>-<sup>3</sup>A<sub>2g</sub> (kcal/mol)</th><th colspan="3"><sup>5</sup>A<sub>1g</sub>-<sup>3</sup>A<sub>2g</sub> (kcal/mol)</th></tr><tr><th>DZ</th><th>TZ</th><th>QZ</th><th>DZ</th><th>TZ</th><th>QZ</th></tr></thead><tbody><tr><td rowspan="2">CAS(8,11)</td><td>CASSCF</td><td>2.91</td><td>3.03</td><td>3.08</td><td>-18.33</td><td>-17.43</td><td>-17.57</td></tr><tr><td>NEVPT2</td><td>6.54</td><td>5.44</td><td>5.20</td><td>-11.85</td><td>-6.95</td><td>-5.39</td></tr><tr><td rowspan="2">CAS(16,15)</td><td>ICE-SCF</td><td>0.00</td><td>0.01</td><td>0.01</td><td>-15.14</td><td>-13.82</td><td>-13.80</td></tr><tr><td>FR-NEVPT2</td><td>0.04</td><td>0.03</td><td>0.05</td><td>-15.39</td><td>-10.06</td><td>-9.14</td></tr><tr><td rowspan="2">CAS(16,19)</td><td>ICE-SCF</td><td>0.02</td><td>0.00</td><td>0.01</td><td>-13.94</td><td>-13.57</td><td>-14.15</td></tr><tr><td>FR-NEVPT2</td><td>-0.10</td><td>-0.23</td><td>0.20</td><td>-16.80</td><td>-12.58</td><td>-10.89</td></tr></tbody></table>

CAS(12,15). Utilizing the QZ basis set, the FR-NEVPT2 results for all states are nearly degenerate as well. Based on the CAS(20,21) active space, we attempted to include two σ-ligand MOs. However, we were unable to achieve convergence for all calculations. At convergence, the $3p_x3p_y$ orbitals of Fe enter the active space, while the two σ-ligand MOs are rotated out (they possess identical irreducible representations). To incorporate the two σ-ligand MOs shown in

Fig. 5, the $3s3p4s4p$ orbitals of Fe must also be included in the active space, leading to CAS(30,29). In Table V, compared with the results employing CAS(12,15) and CAS(20,21), the FR-NEVPT2 method with the active space CAS(30,29) delivers slightly worse $^5$A$_{1g}$-^3A$_{2g}$ gaps. The $^5$A$_{1g}$ state is predicted to be 7.7 (4.9) kcal/mol more stable than the $^3$A$_{2g}$ state with QZ (TZ) basis set. Again, the results indicate that the inclusion of $3s3p4s4p$ could stabilize the triplet states
<figure><img src="image_14.png" /><figcaption>FIG. 5. Active MOs used in the calculations of FeP.</figcaption></figure>

more at the ICE-SCF level but stabilizes the quintet state more at the FR-NEVPT2 level.

It is noteworthy that the $ ^5$A$_1$g state still has about 7.7 kcal/mol lower energy than the triplet states even with the QZ basis set and the CAS(30,29) active space. We believe that this is mainly due to the Hartree-Fock "curse" of grossly overestimating the stability of high-spin states. The energy difference between two states with different multiplicity and different electronic configuration contains contributions from the ligand-field splitting and the changes in the electron-electron repulsion. The latter is largely dominated by the fact that in the low multiplicity, one has more opposite spin electron pairs. These electron pairs have much larger pair-correlation energies than same-spin pairs. Consequently, a physically correct prediction of such spin-state energies requires the full recovery of all differential correlation energies at the basis set limit, a task that is exceedingly difficult to achieve. From Table IV, it is evident that the ICE-SCF results are rather insensitive to the basis set extension, presumably because the chosen active spaces still recover rather little dynamic correlation. It is nevertheless gratifying to see that the

FR-NEVPT2 results do get better with increasing basis set size in most cases, as one would expect, since one is getting closer to the basis set limit. However, even at with a quadruple-ξ basis, one is still not at the basis set limit, and in addition, even with a large active space, NEVPT2 is still not at the FCI limit. Consequently, it is to be expected that the quintet state is still overestimated in our calculation. Even larger active space and basis sets are required to deliver qualitatively correct spin-state energy order. To further illustrate the challenge, we note that unpublished results from our group indicate that the strong bias for high-spin states (in FeTPP, but we believe in general) is still not cured at the DLPNO-CCSD(T)/def2-QZVPP level.⁹⁷,⁹⁸ More theoretical and experimental research studies about FeP are ongoing in our group and will be reported in due course.

### E. Polyacene

The S-T gaps of polyacene molecules have been investigated using sCI based DSRG-MRPT2 methods,⁶⁹,⁷⁰ DMRG-based
<table><caption>TABLE V. Relative energies of $^3E_g$ and $^5A_{1g}$ states with respect to the $^3A_{2g}$ state of FeP, employing different active spaces in Fig. 5 with TZ and QZ basis sets. The main difference to Table IV is the addition of the ligand orbitals, e.g., four "Gouterman orbitals."</caption><thead><tr><th rowspan="2">Active space</th><th rowspan="2">Method</th><th colspan="2">$^3E_g-^3A_{2g}$ (kcal/mol)</th><th colspan="2">$^5A_{1g}-^3A_{2g}$ (kcal/mol)</th></tr><tr><th>TZ</th><th>QZ</th><th>TZ</th><th>QZ</th></tr></thead><tbody><tr><td rowspan="2">CAS(12,15)</td><td>ICE-SCF</td><td>0.02</td><td>0.04</td><td>-17.74</td><td>-17.87</td></tr><tr><td>NEVPT2</td><td>-0.01</td><td>-0.04</td><td>-1.16</td><td>0.43</td></tr><tr><td rowspan="2">CAS(20,21)</td><td>ICE-SCF</td><td>0.04</td><td>0.13</td><td>-18.80</td><td>-19.41</td></tr><tr><td>FR-NEVPT2</td><td>-0.16</td><td>0.40</td><td>-2.35</td><td>-0.38</td></tr><tr><td rowspan="2">CAS(30,29)</td><td>ICE-SCF</td><td>0.18</td><td>2.09</td><td>-12.30</td><td>-13.26</td></tr><tr><td>FR-NEVPT2</td><td>-0.01</td><td>1.28</td><td>-4.94</td><td>-7.69</td></tr></tbody></table>

CASPT²⁹⁹ and MRCI⁴⁸ methods, as well as DMRG pair density functional theory.¹⁰⁰ In this subsection, the adiabatic S–T gaps of n-acenes (n = 3–8) are studied by FR-NEVPT2. For n-acenes (n = 3–4), with fewer than 20 active MOs, TGEN = 1 × 10⁻⁴ settings are used, whereas for pentacene and hexacene, TGEN is set as 5 × 10⁻⁴. For the calculations of heptacene and octacene, which involve more than 30 active electrons and orbitals, TGEN is set to 1 × 10⁻³. Despite aggressive parameters are used for heptacene and octacene, the ICE-CI reference wave functions of the triplet state calculations already include 23 836 and 34 860 CSFs, respectively. The RI CSF space comprises over 380 (866) million CSFs for heptacene (octacene). For the triplet state of octacene with the TZ basis set, the computation required 33 days to complete utilizing 20 cores, performed on a cluster equipped with an AMD EPYC 7H12 CPU. In the FR-NEVPT2 calculation, about 95.6% of the computational time was dedicated to the calculation of various RDMs and Koopmans matrices. In particular, 70.1% of the cost was spent on the computation of matrix elements in rank reduced NEVPT2, while 25.5% was spent on calculating the residuals in FR-NEVPT2, which implicitly involves the fifth-order RDMs. The FR-NEVPT2 calculations of triplet heptacene and octacene with tighter TGEN thresholds may require even more time. To assess the convergence of FR-NEVPT2 results with respect to TGEN, the calculations of hexacene with TGEN set to 1 × 10⁻⁴, 5 × 10⁻⁴ and 1 × 10⁻³ are further conducted. Detailed information is provided in the supplementary material (Table S3). Even with TGEN = 1 × 10⁻⁴ settings, the reference ICE-CI energies are far from convergence. However, as in the case of tetracene, the absolute FR-NEVPT2 energies with the different TGEN thresholds are almost the same. The absolute deviations are less than 1.5 (0.5) milli-hartree for the singlet (triplet) state. Given that the FR-NEVPT2 calculations with TVAR = TGEN * 10⁻³ = 1 × 10⁻⁶ settings provide accurate absolute and relative energies for hexacene, the heptacene and octacene results using TGEN = 1 × 10⁻³ are also expected to be reliable.

The adiabatic S–T gaps predicted by ICE-SCF and FR-NEVPT2 with both DZ and TZ basis sets are given in Fig. 6. The S–T gaps determined by ICE-SCF, which lacks dynamic correlation, are overestimated compared with those by FR-NEVPT2. In the case of octacene, the largest system studied, the inclusion of dynamic correlations reduces the energy gap by ~7 kcal/mol. The FR-NEVPT2

<figure><img src="image_15.png" /><figcaption>FIG. 6. Adiabatic S–T gaps of n-acenes (n = 3–7) computed by ICE-SCF and FR-NEVPT2 with DZ and TZ basis sets compared with DSRG-MRPT2,⁶⁹,⁷⁰ DMRG-CASPT²⁹⁹ and DMRG-ec-MRCISD+Q⁴⁸ results.</figcaption></figure>

results obtained with both DZ and TZ basis sets are found to be nearly identical. The basis set has little influence for the adiabatic S–T gaps of polyacenes. However, the gaps delivered by the TZ basis sets are slightly lower than that with the DZ basis set, except for octacene. For anthracene and tetracene, the S–T gaps calculated by FR-NEVPT2 are larger than those obtained from DMRG-CASPT2 and DMRG-ec-MRCISD methods. For anthracene, the NEVPT2 calculations can be performed without any approximation. The S–T gap predicted by the canonical NEVPT2 method matches with that from FR-NEVPT2 using TGEN = 1 × 10⁻⁴ settings. Thus, the overestimation of the S–T gap for anthracene by FR-NEVPT2 stems from the NEVPT2 theory. The sCI based DSRG-MRPT2 methods⁶⁹,⁷⁰ and DMRG-CASPT2 method³⁹ yield lower energy gaps than FR-NEVPT2 does for all n-acenes. However, from pentacene to octacene, the FR-NEVPT2 results are in good agreement
with those from DMRG-ec-MRCISD. For octacene, the S–T gap by
FR-NEVPT2 with TZ basis set is only 0.2 kcal/mol lower than that
of DMRG-ec-MRCISD; the present FR-NEVPT2 results should be
reliable.

IV. CONCLUSION

The development of computationally efficient and robust mul-
tireference methods for strongly correlated systems with large active
spaces is one of the most active topics in contemporary quantum
chemical research. Although studying systems with an active space
around CAS(60,60) is feasible at the MCSCF level, the post-MCSCF
methods, such as MRPT2, with such large active spaces remain a
field of active investigation. In Paper I of the series, various approxi-
mations to reduce the computational costs of rank reduced NEVPT2
for systems with large active spaces are evaluated. In Paper II of
this series, the FR-NEVPT2 theory was proposed and proven to
be a stable method. In the present work, we implemented an effi-
cient FR-NEVPT2 algorithm for a selected CI method developed by
some of us, ICE-SCF. The current FR-NEVPT2 method, together
with the ICE-SCF method, is capable of delivering reliable static
and dynamic correlation energies for systems with active spaces
as large as CAS(34,34). Interesting systems, such as Fe porphyrin
and polyacenes, have been studied using the present FR-NEVPT2
method. We found that although the static correlation predicted
by ICE-SCF may be far from convergence, the remaining static
correlation energies within the active space can be recovered by
FR-NEVPT2. Thus, even with an ICE-SCF reference consisting of
about one thousand CSFs, the FR-NEVPT2 can deliver reliable
absolute energies at the MRPT2 level. In the future, the domain-
based local pair natural orbital method, as well as explicitly cor-
related methods, will be integrated into the present FR-NEVPT2
algorithm. The new methods, capable of producing basis set lim-
its for large molecules with large active spaces, will be reported
in due course.

SUPPLEMENTARY MATERIAL

See the supplementary material for detailed absolute ICE-SCF
and FR-NEVPT2 energies of iron-porphyrin and polyacenes. Exam-
ple input files to perform ICE-SCF and FR-NEVPT2 calculations are
given as well.

ACKNOWLEDGMENTS

The research of this work was supported by the National
Natural Science Foundation of China (Grant Nos. 22273052 and
22433001). K.S. and F.N. gratefully acknowledge the financial sup-
port by the Max Planck Society and the Cluster of Excellence
(RESOLV, University of Bochum, EXC 1069).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

**Yang Guo:** Conceptualization (equal); Software (equal); Writing – original draft (lead); Writing – review & editing (equal). **Kantharuban Sivalingam:** Conceptualization (equal); Software (equal); Writing – review & editing (equal). **Vijay Gopal Chilkuri:** Software (equal); Writing – review & editing (equal). **Frank Neese:** Conceptualization (lead); Funding acquisition (lead); Software (lead); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available
within the article and its [supplementary material](#).

REFERENCES

¹J. L. Whitten and M. Hackmeyer, *J. Chem. Phys.* **51**, 5584 (1969).

²B. Huron, J. P. Malrieu, and P. Rancurel, *J. Chem. Phys.* **58**, 5745 (1973).

³R. J. Buenker and S. D. Peyerimhoff, *Theor. Chim. Acta* **35**, 33 (1974).

⁴R. J. Buenker and S. D. Peyerimhoff, *Theor. Chim. Acta* **39**, 217 (1975).

⁵W. A. Goddard, T. H. Dunning, W. J. Hunt, and P. J. Hay, *Acc. Chem. Res.* **6**, 368 (1973).

⁶S. Evangelisti, J.-P. Daudey, and J.-P. Malrieu, *Chem. Phys.* **75**, 91 (1983).

⁷C. Angeli and R. Cimiraglia, *Theor. Chem. Acc.* **105**, 259 (2001).

⁸C. Angeli, R. Cimiraglia, M. Persico, and A. Toniolo, *Theor. Chem. Acc.* **98**, 57 (1997).

⁹C. Angeli and M. Persico, *Theor. Chem. Acc.* **98**, 117 (1997).

¹⁰C. Angeli, R. Cimiraglia, and M. Persico, *Theor. Chem. Acc.* **100**, 324 (1998).

¹¹B. O. Roos, P. R. Taylor, and P. E. M. Sigbahn, *Chem. Phys.* **48**, 157 (1980).

¹²S. P. Walch, C. W. Bauschlicher, B. O. Roos, and C. J. Nelin, *Chem. Phys. Lett.* **103**, 175 (1983).

¹³J. Olsen, B. O. Roos, P. Jørgensen, and H. J. A. Jansen, *J. Chem. Phys.* **89**, 2185 (1988).

¹⁴T. Fleig, J. Olsen, and C. M. Marian, *J. Chem. Phys.* **114**, 4775 (2001).

¹⁵T. Fleig, J. Olsen, and L. Visscher, *J. Chem. Phys.* **119**, 2963 (2003).

¹⁶D. Ma, G. Li Manni, and L. Gagliardi, *J. Chem. Phys.* **135**, 044128 (2011).

¹⁷J. Ivanic, *J. Chem. Phys.* **119**, 9364 (2003).

¹⁸N. M. Tubman, D. S. Levine, D. Hait, M. Head-Gordon, and K. B. Whaley, *arXiv:1808.02049* (2018).

¹⁹D. S. Levine, D. Hait, N. M. Tubman, S. Lehtola, K. B. Whaley, and M. Head-Gordon, *J. Chem. Theory Comput.* **16**, 2340 (2020).

²⁰F. A. Evangelista, *J. Chem. Phys.* **140**, 124114 (2014).

²¹J. B. Schriber and F. A. Evangelista, *J. Chem. Phys.* **144**, 161106 (2016).

²²V. G. Chilkuri and F. Neese, *J. Comput. Chem.* **42**, 982 (2021).

²³V. G. Chilkuri and F. Neese, *J. Chem. Theory Comput.* **17**, 2868 (2021).

²⁴N. Zhang, W. Liu, and M. R. Hoffmann, *J. Chem. Theory Comput.* **16**, 2296 (2020).

²⁵N. Zhang, W. Liu, and M. R. Hoffmann, *J. Chem. Theory Comput.* **17**, 949 (2021).

²⁶N. Zhang, Y. Xiao, and W. Liu, *J. Phys.: Condens. Matter* **34**, 224007 (2022).

²⁷Y. Guo, N. Zhang, Y. Lei, and W. Liu, *J. Chem. Theory Comput.* **17**, 7545 (2021).

²⁸Y. Guo, N. Zhang, and W. Liu, *J. Chem. Theory Comput.* **19**, 6668 (2023).

²⁹J. E. T. Smith, B. Mussard, A. A. Holmes, and S. Sharma, *J. Chem. Theory Comput.* **13**, 5468 (2017).

³⁰Y. Garniron, A. Scemama, E. Giner, M. Caffarel, and P.-F. Loos, *J. Chem. Phys.* **149**, 064103 (2018).

³¹P. M. Zimmerman and A. E. Rask, *J. Chem. Phys.* **150**, 244117 (2019).

³²Y. Yao and C. J. Umrigar, *J. Chem. Theory Comput.* **17**, 4183 (2021).

³³J. W. Park, *J. Chem. Theory Comput.* **17**, 1522 (2021).

³⁴G. H. Booth, A. J. W. Thom, and A. Alavi, *J. Chem. Phys.* **131**, 054106 (2009).
<footnote n="35">³⁵R. E. Thomas, Q. Sun, A. Alavi, and G. H. Booth, *J. Chem. Theory Comput.* **11**, 5316 (2015).</footnote>

<footnote n="36">³⁶W. Dobrautz, S. D. Smart, and A. Alavi, *J. Chem. Phys.* **151**, 094104 (2019).</footnote>

<footnote n="37">³⁷E. Giner, A. Scemama, and M. Caffarel, *Can. J. Chem.* **91**, 879 (2013).</footnote>

<footnote n="38">³⁸G. K.-L. Chan and M. Head-Gordon, *J. Chem. Phys.* **116**, 4462 (2002).</footnote>

<footnote n="39">³⁹D. Ghosh, J. Hachmann, T. Yanai, and G. K.-L. Chan, *J. Chem. Phys.* **128**, 144117 (2008).</footnote>

<footnote n="40">⁴⁰D. Zgid and M. Nooijen, *J. Chem. Phys.* **128**, 144116 (2008).</footnote>

<footnote n="41">⁴¹T. Yanai, Y. Kurashige, D. Ghosh, and G. K.-L. Chan, *Int. J. Quantum Chem.* **109**, 2178 (2009).</footnote>

<footnote n="42">⁴²S. Wouters, T. Bogaerts, P. Van Der Voort, V. Van Speybroeck, and D. Van Neck, *J. Chem. Phys.* **140**, 241103 (2014).</footnote>

<footnote n="43">⁴³H. Ma and Y. Ma, *J. Chem. Phys.* **138**, 224105 (2013).</footnote>

<footnote n="44">⁴⁴Y. Ma, S. Knecht, S. Keller, and M. Reiher, *J. Chem. Theory Comput.* **13**, 2533 (2017).</footnote>

<footnote n="45">⁴⁵Q. Sun, J. Yang, and G. K.-L. Chan, *Chem. Phys. Lett.* **683**, 291 (2017).</footnote>

<footnote n="46">⁴⁶T. Yanai, Y. Kurashige, E. Neuscamman, and G. K.-L. Chan, *J. Chem. Phys.* **132**, 024105 (2010).</footnote>

<footnote n="47">⁴⁷M. Saitow, Y. Kurashige, and T. Yanai, *J. Chem. Phys.* **139**, 044118 (2013).</footnote>

<footnote n="48">⁴⁸Z. Luo, Y. Ma, X. Wang, and H. Ma, *J. Chem. Theory Comput.* **14**, 4747 (2018).</footnote>

<footnote n="49">⁴⁹Y. Song, Y. Chen, Y. Ma, and H. Ma, *Electron. Struct.* **2**, 014002 (2020).</footnote>

<footnote n="50">⁵⁰S. Sharma and A. Alavi, *J. Chem. Phys.* **143**, 102815 (2015).</footnote>

<footnote n="51">⁵¹A. Mahajan, N. S. Blunt, I. Sabzevari, and S. Sharma, *J. Chem. Phys.* **151**, 211102 (2019).</footnote>

<footnote n="52">⁵²K. Andersson, P. A. Malmqvist, B. O. Roos, A. J. Sadlej, and K. Wolinski, *J. Phys. Chem.* **94**, 5483 (1990).</footnote>

<footnote n="53">⁵³P. Celani and H.-J. Werner, *J. Chem. Phys.* **112**, 5546 (2000).</footnote>

<footnote n="54">⁵⁴P. Å. Malmqvist, K. Pierloot, A. R. M. Shahi, C. J. Cramer, and L. Gagliardi, *J. Chem. Phys.* **128**, 204109 (2008).</footnote>

<footnote n="55">⁵⁵Y. Kurashige and T. Yanai, *J. Chem. Phys.* **135**, 094104 (2011).</footnote>

<footnote n="56">⁵⁶Y. Kurashige, J. Chalupský, T. N. Lan, and T. Yanai, *J. Chem. Phys.* **141**, 174111 (2014).</footnote>

<footnote n="57">⁵⁷T. Yanai, M. Saitow, X.-G. Xiong, J. Chalupský, Y. Kurashige, S. Guo, and S. Sharma, *J. Chem. Theory Comput.* **13**, 4829 (2017).</footnote>

<footnote n="58">⁵⁸S. Wouters, V. Van Speybroeck, and D. Van Neck, *J. Chem. Phys.* **145**, 054120 (2016).</footnote>

<footnote n="59">⁵⁹Q. M. Phung, S. Wouters, and K. Pierloot, *J. Chem. Theory Comput.* **12**, 4352 (2016).</footnote>

<footnote n="60">⁶⁰N. Nakatani and S. Guo, *J. Chem. Phys.* **146**, 094102 (2017).</footnote>

<footnote n="61">⁶¹S. Guo, M. A. Watson, W. Hu, Q. Sun, and G. K.-L. Chan, *J. Chem. Theory Comput.* **12**, 1583 (2016).</footnote>

<footnote n="62">⁶²M. Roemelt, S. Guo, and G. K.-L. Chan, *J. Chem. Phys.* **144**, 204113 (2016).</footnote>

<footnote n="63">⁶³A. Y. Sokolov, S. Guo, E. Ronca, and G. K.-L. Chan, *J. Chem. Phys.* **146**, 244102 (2017).</footnote>

<footnote n="64">⁶⁴L. Freitag, S. Knecht, C. Angeli, and M. Reiher, *J. Chem. Theory Comput.* **13**, 451 (2017).</footnote>

<footnote n="65">⁶⁵S. Sharma, G. Knizia, S. Guo, and A. Alavi, *J. Chem. Theory Comput.* **13**, 488 (2017).</footnote>

<footnote n="66">⁶⁶R. J. Anderson, T. Shiozaki, and G. H. Booth, *J. Chem. Phys.* **152**, 054101 (2020).</footnote>

<footnote n="67">⁶⁷E. Xu and S. Li, *J. Chem. Phys.* **139**, 174111 (2013).</footnote>

<footnote n="68">⁶⁸E. Xu, D. Zhao, and S. Li, *J. Chem. Theory Comput.* **11**, 4634 (2015).</footnote>

<footnote n="69">⁶⁹J. B. Schriber, K. P. Hannon, C. Li, and F. A. Evangelista, *J. Chem. Theory Comput.* **14**, 6295 (2018).</footnote>

<footnote n="70">⁷⁰J. W. Park, *J. Chem. Theory Comput.* **19**, 6263 (2023).</footnote>

<footnote n="71">⁷¹J.-P. Malrieu and J.-L. Heully, *J. Chem. Phys.* **159**, 034110 (2023).</footnote>

<footnote n="72">⁷²P. Beran, M. Matoušek, M. Hapka, K. Pernal, and L. Veis, *J. Chem. Theory Comput.* **17**, 7575 (2021).</footnote>

<footnote n="73">⁷³C. Angeli, R. Cimiraglia, and J.-P. Malrieu, *J. Chem. Phys.* **117**, 9138 (2002).</footnote>

<footnote n="74">⁷⁴C. Angeli, M. Pastore, and R. Cimiraglia, *Theor. Chem. Acc.* **117**, 743 (2007).</footnote>

<footnote n="75">⁷⁵Y. Guo, K. Sivalingam, E. F. Valeev, and F. Neese, *J. Chem. Phys.* **144**, 094111 (2016).</footnote>

<footnote n="76">⁷⁶Y. Guo, K. Sivalingam, E. F. Valeev, and F. Neese, *J. Chem. Phys.* **147**, 064110 (2017).</footnote>

<footnote n="77">⁷⁷Y. Guo, F. Pavošević, K. Sivalingam, U. Becker, E. F. Valeev, and F. Neese, *J. Chem. Phys.* **158**, 094111 (2023).</footnote>

<footnote n="78">⁷⁸Y. Guo, K. Sivalingam, and F. Neese, *J. Chem. Phys.* **154**, 214111 (2021).</footnote>

<footnote n="79">⁷⁹Y. Guo, K. Sivalingam, C. Kollmar, and F. Neese, *J. Chem. Phys.* **154**, 214113 (2021).</footnote>

<footnote n="80">⁸⁰K. G. Dyall, *J. Chem. Phys.* **102**, 4909 (1995).</footnote>

<footnote n="81">⁸¹C. Angeli, R. Cimiraglia, and J.-P. Malrieu, *Chem. Phys. Lett.* **350**, 297 (2001).</footnote>

<footnote n="82">⁸²C. Kollmar, K. Sivalingam, Y. Guo, and F. Neese, *J. Chem. Phys.* **155**, 234104 (2021).</footnote>

<footnote n="83">⁸³C. Kollmar, K. Sivalingam, B. Helmich-Paris, C. Angeli, and F. Neese, *J. Comput. Chem.* **40**, 1463 (2019).</footnote>

<footnote n="84">⁸⁴K. Pierloot, Q. M. Phung, and A. Domingo, *J. Chem. Theory Comput.* **13**, 537 (2017).</footnote>

<footnote n="85">⁸⁵C. Angeli, *J. Comput. Chem.* **30**, 1319 (2009).</footnote>

<footnote n="86">⁸⁶G. Li Manni, S. D. Smart, and A. Alavi, *J. Chem. Theory Comput.* **12**, 1245 (2016).</footnote>

<footnote n="87">⁸⁷M. Radoń, *J. Chem. Theory Comput.* **10**, 2306 (2014).</footnote>

<footnote n="88">⁸⁸M. Tarrago, C. Römelt, J. Nehrkorn, A. Schnegg, F. Neese, E. Bill, and S. Ye, *Inorg. Chem.* **60**, 4966 (2021).</footnote>

<footnote n="89">⁸⁹B. J. Persson, B. O. Roos, and K. Pierloot, *J. Chem. Phys.* **101**, 6810 (1994).</footnote>

<footnote n="90">⁹⁰M. Feldt, Q. M. Phung, K. Pierloot, R. A. Mata, and J. N. Harvey, *J. Chem. Theory Comput.* **15**, 922 (2019).</footnote>

<footnote n="91">⁹¹S. Vancoillie, H. Zhao, V. T. Tran, M. F. A. Hendrickx, and K. Pierloot, *J. Chem. Theory Comput.* **7**, 3961 (2011).</footnote>

<footnote n="92">⁹²Q. M. Phung, M. Feldt, J. N. Harvey, and K. Pierloot, *J. Chem. Theory Comput.* **14**, 2446 (2018).</footnote>

<footnote n="93">⁹³G. Li Manni and A. Alavi, *J. Phys. Chem. A* **122**, 4935 (2018).</footnote>

<footnote n="94">⁹⁴G. Li Manni, D. Kats, D. P. Tew, and A. Alavi, *J. Chem. Theory Comput.* **15**, 1492 (2019).</footnote>

<footnote n="95">⁹⁵M. Tarrago, S. Ye, and F. Neese, *Chem. Sci.* **13**, 10029 (2022).</footnote>

<footnote n="96">⁹⁶M. Gouterman, *J. Mol. Spectrosc.* **6**, 138 (1961).</footnote>

<footnote n="97">⁹⁷M. Saitow, U. Becker, C. Riplinger, E. F. Valeev, and F. Neese, *J. Chem. Phys.* **146**, 164105 (2017).</footnote>

<footnote n="98">⁹⁸Y. Guo, C. Riplinger, D. G. Liakos, U. Becker, M. Saitow, and F. Neese, *J. Chem. Phys.* **152**, 024116 (2020).</footnote>

<footnote n="99">⁹⁹Y. Kurashige and T. Yanai, *Bull. Chem. Soc. Jpn.* **87**, 1071 (2014).</footnote>

<footnote n="100">¹⁰⁰P. Sharma, V. Bernales, S. Knecht, D. G. Truhlar, and L. Gagliardi, *Chem. Sci.* **10**, 1716 (2019).</footnote>