# iCISCF: An Iterative Configuration Interaction-based Multiconfigurational
Self-consistent Field Theory for Large Active Spaces

Yang Guo  Qingdao Institute for Theoretical and Computational Sciences,
Institute of Frontier and Interdisciplinary Science, Shandong University,
Qingdao, Shandong 266237, China  Ning Zhang  Beijing National Laboratory for
Molecular Sciences, Institute of Theoretical and Computational Chemistry,
College of Chemistry and Molecular Engineering, Peking University, Beijing
100871, China  Yibo Lei  Key Laboratory of Synthetic and Natural Functional
Molecule Chemistry of Ministry of Education, College of Chemistry and
Materials Science, Shaanxi Key Laboratory of Physico-Inorganic Chemistry,
Northwest University, Xi’an, Shaanxi 710127, China  Wenjian Liu
<EMAIL_ADDRESS>Qingdao Institute for Theoretical and Computational Sciences,
Institute of Frontier and Interdisciplinary Science, Shandong University,
Qingdao, Shandong 266237, China

###### Abstract

An iterative configuration interaction (iCI)-based multiconfigurational self-
consistent field (SCF) theory, iCISCF, is proposed to handle systems that
require large complete active spaces (CAS). The success of iCISCF stems from
three ingredients: (1) efficient selection of individual configuration state
functions spanning the CAS, meanwhile maintaining full spin symmetry; (2) the
use of Jacobi rotation for the optimization of active orbitals, in conjunction
with a quasi-Newton algorithm for the core/active-virtual and core-active
orbital rotations; (3) a second-order perturbative treatment of the residual
space left over by the selection procedure (i.e., iCISCF(2)). Just like
selected iCI being a very accurate approximation to CASCI, iCISCF(2) is a very
accurate approximation to CASSCF. Several examples that go beyond the
capability of CASSCF are taken as showcases to reveal the performances of
iCISCF and iCISCF(2).

##  1 Introduction

It is well known that a multiconfigurational wave function is demanded even
for a qualitative description of a chemical system with a dense set of
energetically adjacent frontier molecular orbitals (MO). The complete active
space self-consistent field (CASSCF) theory1, 2, 3, 4, 5, 6, 7 is usually
taken as the first step towards an accurate description of such strongly
correlated systems. The popularity of CASSCF stems from its operational
simplicity: starting with the partition of the in total $M$ MOs into $M_{c}$
core (always doubly occupied), $M_{a}$ active (variably occupied) and
$M_{v}=M-M_{c}-M_{a}$ virtual (always zero occupied) orbtials, all
configurations are then generated by distributing the $N_{a}=N_{e}-2M_{c}$
active electrons in the $M_{a}$ active orbitals in all possible ways, thereby
leading to a subspace full configuration interaction (FCI) problem, usually
denoted as CASCI. While the total energy cannot be altered by arbitrary
rotations (unitary transformations) within each of the three subsets of
orbitals, the rotations in between do lower the energy further. Therefore, a
CASSCF calculation amounts to optimizing the CI and MO coefficients
simultaneously to make the energy stationary. Yet, given so many years of
algorithmic developments (see Refs. 8, 9, 10 for the most recent ones), the
largest CASSCF calculation performed so far involves only 22 active electrons
in 22 active orbitals, which was made possible only by massive
parallelization11. The huge computational cost stems of course from the
exponential growing of the size of CASCI. To break the record, some
approximation must be introduced to simplify the CASCI calculation. As a
matter of fact, as long as the approximation can be controlled so as to
guarantee the quality of the orbitals, the approximation itself is not an
issue at all, simply because CASSCF is not the end of the calculation but
rather the preparation for subsequent treatment of dynamic correlation, which
can fully recover the marginal loss due to the approximate treatment of the CI
coefficients. A large number of approaches have been developed in the last
decades to approximate FCI/CASCI (see Ref. 12 for a complete collection and
classification of such approaches). In particular, some of these13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 have been adapted to
CASSCF, enabling much larger CASSCF calculations. The present work aims to
combine the iterative configuration interaction (iCI) approach29 with CASSCF,
leading to iCISCF as a new member of the static-dynamic-static (SDS)30 family
of methods (SDSPT230, 31, SDSCI29, 30, iCI29, iCIPT232, 12, iVI (iterative
vector interaction)33, 34, $\mathbbm{i}$CAS (imposed automatic selection and
localization of active orbitals)35, and iOI (iterative orbital
interaction)36). Among the various unique features12 of iCI29, the following
ones are particularly relevant to CASSCF:

  1. (1)

Configuration state functions (CSF) instead of Slater determinants (SD) are
taken as the many-electron basis so as to guarantee full spin symmetry, which
is of vital importance for describing low-spin states of general open-shell
systems.

  2. (2)

Individual CSFs can be selected with a very efficient algorithm.

  3. (3)

A particle-hole representation is employed to establish connections between
randomly selected CSFs so as to construct the Hamiltonian matrix very
efficiently in a compressed form. In particular, the connections between hole-
strings and between particle-strings can be shared by CSF spaces of different
spatial and/or spin symmetries, thereby facilitating the simultaneous
calculation of several states of different spatial and/or spin symmetries with
a common set of orthonormal orbitals.

  4. (4)

The diagonalization of the Hamiltonian matrix is done by iVI33, 34, which is
able to access directly the roots of a given energy window or specified
characters, thereby guiding the SCF iterations to converge to excited states
free of variational collapse (a point that is not pursued here though).

The iCISCF algorithm is described in Sec. 2, which is followed by pilot
applications in Sec. 3. The following conventions are to be used throughout:
(1) core, active, virtual, and arbitrary MOs are designated by
$\\{i,j,k,l,\cdots\\}$, $\\{t,u,v,w,\cdots\\}$, $\\{a,b,c,d,\cdots\\}$, and
$\\{o,p,q,r,s,\cdots\\}$, respectively; (2) repeated indices are always summed
up.

##  2 iCISCF

Given a set of MOs that are partitioned into core, active and virtual ones,
the selection of individual CSFs in the CASCI space $P$ is performed
iteratively until convergence, with a single parameter $C_{\mathrm{min}}$ for
controlling the size of the final, selected space $P_{m}$. Specifically,
starting with a guess space $P_{0}$, only those CSFs $\\{|I\rangle\in
P-P_{0}\\}$ satisfying the ranking criterion $\max_{J\in
P_{0}}|\frac{H_{IJ}C_{J}}{E^{(0)}-H_{II}}|\geq C_{\mathrm{min}}$ are put into
$P_{0}$, so as to extend $P_{0}$ to $P_{1}$; the Hamiltonian matrix in $P_{1}$
is then diagonalized33, 34 to prune away those CSFs of coefficients smaller in
magnitude than $C_{\mathrm{min}}$, so as to reduce $P_{1}$ to $P_{m}$ upon
convergence (for more details, see Ref. 12). Finally, the Hamiltonian matrix
$\mathbf{H}^{(0)}$ in $P_{m}$, viz.

$\displaystyle H_{IJ}^{(0)}$ $\displaystyle=$ $\displaystyle\langle
I|\hat{H}|J\rangle,\quad\forall I,J\in P_{m},$ (1) $\displaystyle\hat{H}$
$\displaystyle=$ $\displaystyle
h_{pq}\hat{E}_{pq}+\frac{1}{2}g_{pq,rs}\hat{e}_{pq,rs},$ (2)
$\displaystyle\hat{E}_{pq}$ $\displaystyle=$ $\displaystyle
a_{p\alpha}^{\dagger}a_{q\alpha}+a_{p\beta}^{\dagger}a_{q\beta},$ (3)
$\displaystyle\hat{e}_{pq,rs}$ $\displaystyle=$
$\displaystyle\\{\hat{E}_{pq}\hat{E}_{rs}\\}=\\{\hat{E}_{rs}\hat{E}_{pq}\\}=\hat{E}_{pq}\hat{E}_{rs}-\delta_{qr}\hat{E}_{ps},$
(4)

is diagonalized to yield the normalized iCI wave function

$\displaystyle|0\rangle=\sum_{|I\rangle\in P_{m}}|I\rangle
C_{I}^{(0)},\quad\langle I|J\rangle=\delta_{IJ},\quad\sum_{I}C_{I}^{2}=1,$ (5)

with energy

$\displaystyle E^{(0)}$ $\displaystyle=$ $\displaystyle
tr(\mathbf{h}\mathbf{D})+\frac{1}{2}\sum_{rs}tr(\mathbf{J}^{rs}\mathbf{P}^{sr}),$
(6) $\displaystyle D_{pq}$ $\displaystyle=$ $\displaystyle\langle
0|\hat{E}_{pq}|0\rangle=D_{qp},\quad\Gamma_{pq,rs}=\langle
0|\hat{e}_{pq,rs}|0\rangle=\Gamma_{qp,sr},$ (7) $\displaystyle J^{rs}_{pq}$
$\displaystyle=$ $\displaystyle g_{pq,rs},\quad
P^{rs}_{pq}=\frac{1}{2}(\Gamma_{pq,rs}+\Gamma_{qp,rs})=P^{sr}_{pq}=P^{pq}_{rs}.$
(8)

Note in passing that the basic coupling coefficients $\langle
I|\hat{E}_{pq}|J\rangle$ and $\langle I|\hat{e}_{pq,rs}|J\rangle$ in the
Hamiltonian matrix (1) depend only on the relative occupations of the MOs in
the CSF pairs $\\{|I\rangle,|J\rangle\\}$ but not on the individual MOs. As
such, they can be evaluated and reused very efficiently with the tabulated
unitary group approach32. For later use, we here define $P_{s}=P-P_{m}$ as the
residual of $P$ left over by the selection, whereas $Q=1-P$ as the
complementary space for dynamic correlation.

To further optimize the orbitals, we parameterize the iCISCF wave function
as37

$\displaystyle|\tilde{0}\rangle$ $\displaystyle=$ $\displaystyle\sum_{I\in
P_{m}}exp(-\hat{\kappa})|I\rangle C_{I},$ (9) $\displaystyle\hat{\kappa}$
$\displaystyle=$
$\displaystyle\sum_{pq}\kappa_{pq}\hat{E}_{pq}=\sum_{p>q}\kappa_{pq}\hat{E}_{pq}^{-}=\frac{1}{2}\sum_{p,q}\kappa_{pq}\hat{E}_{pq}^{-},$
(10) $\displaystyle\bm{\kappa}$ $\displaystyle=$
$\displaystyle-\bm{\kappa}^{\dagger},\quad\hat{E}_{pq}^{-}=\hat{E}_{pq}-\hat{E}_{qp},\quad\frac{\partial\hat{\kappa}}{\partial\kappa_{rs}}=\hat{E}_{rs}^{-}.$
(11)

The skew symmetry of $\bm{\kappa}$ implies that only $\kappa_{pq}(p>q)$ are
independent parameters. Moreover, rotations within the core and within the
virtual orbitals are also redundant. Therefore, the $\hat{\kappa}$ operator
actually reads

$\displaystyle\hat{\kappa}=\sum_{ai}\kappa_{ai}\hat{E}_{ai}^{-}+\sum_{ti}\kappa_{ti}\hat{E}_{ti}^{-}+\sum_{at}\kappa_{at}\hat{E}_{at}^{-}+\sum_{t>u}\kappa_{tu}\hat{{E}}_{tu}^{-}.$
(12)

The last term in Eq. (12) vanishes in CASSCF but has to be included in iCISCF
due to the truncation of CASCI. It is precisely this term that renders the
orbital optimization of iCISCF more difficult than that of CASSCF. Minimizing
the Lagrangian

$\displaystyle L[\bm{\kappa},\mathbf{C}]$ $\displaystyle=$
$\displaystyle\langle\tilde{0}|\hat{H}|\tilde{0}\rangle-E(\mathbf{C}^{\dagger}\mathbf{C}-\mathbf{I})$
(13)

leads to the stationary conditions

$\displaystyle G^{\mathrm{c}}_{I}$ $\displaystyle=$ $\displaystyle 2\langle
I|\hat{H}-E^{(0)}|0\rangle=0,$ (14) $\displaystyle G^{\mathrm{o}}_{pq}$
$\displaystyle=$ $\displaystyle\langle
0|[\hat{E}_{pq}^{-},\hat{H}]|0\rangle=2\langle
0|[\hat{E}_{pq},\hat{H}]|0\rangle$ (15) $\displaystyle=$ $\displaystyle
2(F-F^{\dagger})_{pq}=0,$ (16) $\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle\mathbf{D}\mathbf{h}+\sum_{rs}\mathbf{P}^{rs}\mathbf{J}^{sr},$
(17)

where the partial derivatives are evaluated at the expansion point
$\bm{\xi}^{(0)}=(\mathbf{C}^{(0)},\mathbf{0})^{T}$. Eq. (16) can be recast38
into a Hartree-Fock-like equation, such that it can be solved by
diagonalization, just like Eq. (14). However, such first-order orbital
optimization method is hardly useful due to very slow convergence. A more
robust method is the second-order Newton-Raphson (QR) scheme, which amounts to
expanding the Lagrangian $L$ (13) to second order at the expansion point
$\bm{\xi}^{(0)}$,

$\displaystyle
L^{(2)}[\bm{\kappa},\mathbf{C}]=E^{(0)}+\bm{\xi}^{(1)T}\mathbf{G}+\frac{1}{2}\bm{\xi}^{(1)T}\mathbf{E}\bm{\xi}^{(1)},$
(18)

the stationary condition of which is

$\displaystyle\mathbf{E}\bm{\xi}^{(1)}=-\mathbf{G}$ (19)

or in block form

$\displaystyle\begin{pmatrix}\mathbf{E}^{\mathrm{cc}}&\mathbf{E}^{\mathrm{co}}\\\
\mathbf{E}^{\mathrm{oc}}&\mathbf{E}^{\mathrm{oo}}\end{pmatrix}\begin{pmatrix}\mathbf{C}^{(1)}\\\
\bm{\kappa}^{(1)}\end{pmatrix}=-\begin{pmatrix}\mathbf{G}^{\mathrm{c}}\\\
\mathbf{G}^{\mathrm{o}}\end{pmatrix},$ (20)

where the second-order partial derivatives (at the expansion point
$\bm{\xi}^{(0)}$) read

$\displaystyle E^{\mathrm{cc}}_{IJ}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}L}{\partial C_{I}\partial C_{J}}=2\langle
I|\hat{H}-E^{(0)}|J\rangle,$ (21) $\displaystyle E^{\mathrm{oc}}_{pq,I}$
$\displaystyle=$
$\displaystyle\frac{\partial^{2}L}{\partial\kappa_{pq}\partial C_{I}}=2\langle
0|[\hat{E}_{pq}^{-},\hat{H}]|I\rangle=E^{\mathrm{co}}_{I,pq},$ (22)
$\displaystyle E^{\mathrm{oo}}_{pq,rs}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}L}{\partial\kappa_{pq}\partial\kappa_{rs}}=\frac{1}{2}\langle
0|[\hat{E}_{pq}^{-},[\hat{E}_{rs}^{-},\hat{H}]]|0\rangle+\frac{1}{2}\langle
0|[\hat{E}_{rs}^{-},[\hat{E}_{pq}^{-},\hat{H}]]|0\rangle$ (23)
$\displaystyle=$ $\displaystyle\langle
0|[\hat{E}_{pq}^{-},[\hat{E}_{rs}^{-},\hat{H}]]|0\rangle+\frac{1}{2}\langle
0|[[\hat{E}_{rs}^{-},\hat{E}_{pq}^{-}],\hat{H}]|0\rangle.$ (24)

The first row of Eq. (20) can be rearranged to

$\displaystyle
2(\mathbf{H}^{(0)}-E^{(0)}\mathbf{I})\mathbf{C}^{(1)}=-\mathbf{G}^{\mathrm{c}}-2\mathbf{H}^{(1)}\mathbf{C}^{(0)},$
(25)

where

$\displaystyle H^{(1)}_{IJ}$ $\displaystyle=$ $\displaystyle\langle
I|\hat{H}_{\kappa}|J\rangle,$ (26) $\displaystyle\hat{H}_{\kappa}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\kappa_{pq}^{(1)}[\hat{E}_{pq}^{-},\hat{H}]=\hat{H}_{\kappa}^{\dagger}$
(27) $\displaystyle=$
$\displaystyle[\bm{\kappa}^{(1)},\mathbf{h}]_{pq}\hat{E}_{pq}+\frac{1}{2}[\bm{\kappa}^{(1)},\mathbf{J}^{rs}]_{pq}\hat{e}_{pq,rs}+\frac{1}{2}[\bm{\kappa}^{(1)},\mathbf{J}^{pq}]_{rs}\hat{e}_{rs,pq}$
(28) $\displaystyle=$
$\displaystyle[\bm{\kappa}^{(1)},\mathbf{h}]_{pq}\hat{E}_{pq}+[\bm{\kappa}^{(1)},\mathbf{J}^{rs}]_{pq}\hat{e}_{pq,rs},$
(29)

in which $\hat{H}_{\kappa}$ is a bona fide Hamiltonian operator37. Likewise,
the second row of Eq. (20) can be rearranged to

$\displaystyle(\mathbf{E}^{\mathrm{oo}}\bm{\kappa}^{(1)})_{pq}$
$\displaystyle=$ $\displaystyle-\mathbf{G}^{\mathrm{o}}_{pq}-2\langle
0|[\hat{E}_{pq}^{-},\hat{H}|J\rangle C_{J}^{(1)},$ (30)
$\displaystyle(\mathbf{E}^{\mathrm{oo}}\bm{\kappa}^{(1)})_{pq}$
$\displaystyle=$ $\displaystyle\langle
0|[\hat{E}_{pq}^{-},\hat{H}_{\kappa}]|0\rangle+[\mathbf{G^{\mathrm{o}}},\bm{\kappa}^{(1)}]_{pq},$
(31) $\displaystyle\langle 0|[\hat{E}_{pq}^{-},\hat{H}|J\rangle C_{J}^{(1)}$
$\displaystyle=$ $\displaystyle
2[\mathbf{D}^{(1)},\mathbf{h}]_{pq}+2\sum_{rs}[\mathbf{P}^{(1)rs},\mathbf{J}^{sr}]_{pq},$
(32) $\displaystyle D^{(1)}_{pq}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle 0|\hat{E}_{pq}+\hat{E}_{qp}|J\rangle
C_{J}^{(1)},$ (33) $\displaystyle P^{(1)rs}_{pq}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle 0|\hat{e}_{pq,rs}+\hat{e}_{qp,rs}|J\rangle
C_{J}^{(1)}.$ (34)

Eqs. (25) and (30) can be viewed as response equations for the CI
displacements $\mathbf{C}^{(1)}$ and orbital Newton steps $\bm{\kappa}^{(1)}$,
respectively. The former involves the first-order active space Hamiltonian
(26) due to $\bm{\kappa}^{(1)}$, whereas the latter involves the first-order
reduced density matrices (RDM) (33)/(34) due to $\mathbf{C}^{(1)}$. This
particular reformulation21 is advantageous in that it decouples the orbital
optimization from the CI solver implementation in each Newton step, such that
any CI solver can readily be used. An even more sophisticated formulation is
the Werner-Meyer-Knowles (WMK) approach39, 9, 10, where the Lagrangian (13) is
expanded to second order in $\mathbf{T}=e^{-\bm{\kappa}}-\mathbf{I}$, thereby
containing terms infinite order in $\bm{\kappa}$. However, considering that
the simultaneous optimization of the CI coefficients and the orbitals suffers
from severe linear dependence due to the presence of active-active orbital
rotations and that the CI step is rate determining for a large CAS, the second
term on the right-hand side of both Eqs. (25) and (30) can be ignored, thereby
leading to

$\displaystyle(\mathbf{H}^{(0)}-E^{(0)}\mathbf{I})\mathbf{C}^{(1)}$
$\displaystyle=$ $\displaystyle-\mathbf{G}^{\mathrm{c}},$ (35)
$\displaystyle\mathbf{E}^{\mathrm{oo}}\bm{\kappa}^{(1)}$ $\displaystyle=$
$\displaystyle-\mathbf{G}^{\mathrm{o}}.$ (36)

Eq. (35) is the usual iterative partial diagonalization of the CI matrix
$\mathbf{H}^{(0)}$. For the given RDMs, Eq. (36) can, in the spirit of quasi-
Newton (QR) methods, be solved iteratively with the Broyden-Fletcher-Goldfarb-
Shanno (BFGS) algorithm40, 41, which amounts to updating iteratively the
inverse of the orbital Hessian $\mathbf{E}^{\mathrm{oo}}$ with the diagonal
elements (see Appendix A) as the initial guess. While such QN-based CASSCF has
been in use for a long time42, preliminary experimentations showed that the
QN-based iCISCF often does not converge. To handle the troublesome active-
active orbital rotations, we adopt the Jacobi rotation (JR) algorithm, which
is the most robust first-order optimization method. In particular, it has been
shown43 that the convergence pattern of JR-based MCSCF is not much affected
when a not fully converged CI vector is used. In the actual running of iCISCF,
the JR algorithm is first employed to optimize the active orbitals until a
convergence threshold is fulfilled. Then, the QN algorithm is invoked to
determine the remaining core/active-virtual and core-active rotation
parameters. This hybrid algorithm is to be denoted as JR+QN, to distinguish
from the pure QN algorithm. The flowchart of iCISCF is depicted in Fig. 1. If
a state-specific Epstein-Nesbet type of second-order perturbation theory
(ENPT2) is carried out in $P_{s}$, the method will be denoted as iCISCF(2).
Likewise, if the ENPT2 is carried out in the joint space $P_{s}+Q=1-P_{m}$,
the method will be the original iCIPT232, 12. Just like that iCIPT2 is very
close to FCI32, 12, 44, iCISCF(2) is very close to CASSCF. In principle, the
iCISCF(2) energy functional can be minimized for orbital optimization.
However, it has been shown24 that this does not change the orbitals
discernibly.

---
Figure 1: Flowchart of iCISCF using two different active-active rotation
schemes.

##  3 Results and discussion

All calculations were performed with the BDF program package45, 46, 47, 48, 49
on a computer node equipped with Intel(R) Xeon(R) Gold 6240 CPUs (in total 36
physical cores) and 128 GB memory. The convergence threshold for the energy is
10-10 $E_{\mathrm{h}}$. All extrapolations were based on a four-point linear
fit of $E_{\mathrm{iCISCF(2)}}$ versus $E^{(2)}$ calculated with different
Cmin, which is possible because the orbitals with different Cmin are virtually
identical.

###  3.1 Optimization of active orbitals

Compared with CASSCF, the additional couplings between the active-active
orbital rotations and the CI displacements render iCISCF (and other
approximate CASSCF approaches as well) more difficult to converge, especially
when they are optimized simultaneously. To overcome such difficulties, the two
sets of parameters have to be decoupled. Even in this case, special care has
to be taken of. For instance, in the adaptive sampling configuration
interaction-based CASSCF approach (ASCISCF)27, 28, a sequence of supermacro-,
macro-, and micro-iterations is invoked: the selection of a fixed number of
SDs is performed in each supermacroiteration, which consists of a fixed number
of macroiterations, where the list of SDs is kept fixed, such that only the CI
coefficients are allowed to respond to changes in the orbitals resulting from
the microiterations. As such, the macro- and micro-iterations together mimic
CASSCF (where the list of SDs is fixed by construction). In contrast, in
iCISCF, the selection of individual CSFs is performed in each macroiteration.
Since the selection is controlled by a single parameter $C_{\mathrm{min}}$,
both the size and members of the resulting $P_{m}$ space may change along
macroiterations. Therefore, it is of first interest to see how the
optimization of the active orbitals behaves. To this end, the lowest adiabatic
singlet and triplet states of hexacene (six linearly fused benzene rings) are
studied with CAS(26,26)/cc-pVDZ50 and $D_{2h}$ symmetry. As a matter of fact,
identifying the restricted (open-shell) Hartree-Fock (R(O)HF) orbitals as
initial guesses for iCISCF is nontrivial for calculations with large active
space. However, the pre-CMOs obtained by diagonalizing the molecular Fock
matrix projected on the prechosen valence atomic orbitals or the localized
pre-CMOs (pre-LMOs)35, can be employed to construct initial guess CMOs or LMOs
for iCISCF calculations. Thus, in the calculations of hexacene, two different
set of initial guess MOs, CMOs and LMOs, constructed by the $\mathbbm{i}$CAS
method are used as initial guesses.

In both QN and JR+QN algorithms, the optimization of the active orbitals is
turned off when the energy difference between two adjacent macroiterations is
below the threshold $A_{\mathrm{min}}$. The results with different
$A_{\mathrm{min}}$ setting are plotted in Fig. 2, where the energies without
optimizing the active orbitals are taken as the zero points. It can be seen
that both algorithms converge to essentially the same results for the same
$A_{\mathrm{min}}$, and a value of $1.0\times 10^{-6}$ $E_{\mathrm{h}}$ is
required for $A_{\mathbf{min}}$ to get fully converged results. It is also
seen from the Fig. 2 that the initial guess has a grear effect on the
convergence of iCISCF results. By using LMOs as initial guess, internal active
space optimization is not as important as in the calculations with CMOs.
Moreover, the calculations, with different initial MOs, do not always converge
to the same results. For example, by setting C${}_{\mathrm{min}}=5.0\times
10^{-5}$ and $A_{\mathrm{min}}=1.0\times 10^{-8}$, the triplet calculations
with LMOs guesses deliver lower absolute iCISCF energies than those with CMOs
do, which is as large as 6.8 $mE_{\mathrm{h}}$, although the final iCISCF(2)
results are close to each other.

The numbers of macro- and microiterations in these calculations are further
summarized in Table 1. The results show that JR+QN does outperform QN in
general: QN requires somewhat more microiterations than JR+QN and even fails
to converge within 200 macroiterations in two cases. When using a tight
$A_{\mathrm{min}}$ threshold, the JR+QN algorithm needs fewer macroiterations
to converge for most cases. It can also be seen that using LMOs as initial
guesses, the number of macroiteraitons reduces substantially, especially for
calculations with the C${}_{\mathrm{min}}=5.0\times 10^{-4}$ setting. The
triplet state calculations converge faster than that of singlet states.
However, the calculations for singlet states with
C${}_{\mathrm{min}}=5.0\times 10^{-5}$ and $A_{\mathrm{min}}=1.0\times
10^{-8}$ settings still need more than one hundred iterations. The convergence
trends of the calculations with these settings optimized by QN and JR+QN
algorithms are given in Fig. 3. Within the region from $1.0\times 10^{-3}$ to
$1.0\times 10^{-6}E_{\mathrm{h}}$, all calculations converge slowly. This
could be attributed to the coupling between active-active rotation and the CI
coefficients.

In the present work, the JR+QN algorithm is used as the default algorithms for
the iCISCF wave function optimization. The LMOs constructed by
$\mathbbm{i}$CAS from the pre-LMOs are used as initial guesses for all
calculations, if not otherwise specified. By default, $A_{\mathrm{min}}$ is
chosen to be 100 times the energy convergence threshold, such that iCISCF has
only one parameter, Cmin, to play.

|
---|---
(a) C${}_{\mathrm{min}}=5.0\times 10^{-4}$ | (b) C${}_{\mathrm{min}}=5.0\times 10^{-5}$
Figure 2: Convergence patterns of the iCISCF calculations of hexacene. The
energies without active-active rotations are taken as zero points. For
calculations not converged within 200 macroiterations, the iCISCF energies of
the 201st cycle are used. Table 1: Number of macroiterations (MA) of the
iCISCF calculations of hexacene using the quasi-Newton (QN) and hybrid Jacobi
rotation and QN algorithms (JR+QN) with CMOs and LMOs as initial guess,
respectively. Averaged number of microiterations (MI) per macroiteration are
given as well

|  |  C${}_{\mathrm{min}}=5.0\times 10^{-4}$ a |  C${}_{\mathrm{min}}=5.0\times 10^{-5}$ a
---|---|---|---
|  | Singlet | Triplet | Singlet | Triplet
|  | MA | MI | MA | MI | MA | MI | MA | MI
|  |  |  |  | CMOs |  |  |
w/o active-active rotations | 6 | 3.7 | 6 | 3.8 | 7 | 3.3 | 8 | 3.0
|  $A_{\mathrm{min}}$ a |  |  |  |  |  |  |  |
QN | $1.0\times 10^{-2}$ |  n.c.c | 2.2 | 7 | 5.0 | 6 | 6.2 | 8 | 5.1
$1.0\times 10^{-4}$ | 8 | 6.8 | 7 | 6.0 | 7 | 7.5 | 51 | 14.8
$1.0\times 10^{-6}$ | 128 | 11.1 | 133 | 12.0 | 174 | 12.4 | 78 | 14.5
$1.0\times 10^{-8}$ | 168 | 11.1 | 147 | 12.0 |  n.c.c | 12.3 | 80 | 14.6
JR+QN | $1.0\times 10^{-2}$ | 8 | 3.6 | 6 | 4.0 | 7 | 3.6 | 7 | 3.7
$1.0\times 10^{-4}$ | 9 | 3.6 | 8 | 3.5 | 8 | 3.4 | 51 | 3.0
$1.0\times 10^{-6}$ | 125 | 3.4 | 123 | 3.1 | 173 | 2.3 | 72 | 2.7
$1.0\times 10^{-8}$ | 147 | 3.4 | 128 | 3.1 |  n.c.c | 2.3 | 78 | 2.7
|  |  |  |  | LMOs |  |  |
w/o active-active rotations | 6 | 5.3 | 6 | 5.5 | 6 | 6.0 | 12 | 3.3
|  $A_{\mathrm{min}}$ a |  |  |  |  |  |  |  |
QN | $1.0\times 10^{-2}$ | 5 | 8.2 | 6 | 7.5 | 11 | 4.3 | 7 | 6.1
$1.0\times 10^{-4}$ | 9 | 12.0 | 11 | 10.1 | 11 | 7.5 | 14 | 9.0
$1.0\times 10^{-6}$ | 22 | 10.2 | 25 | 11.4 | 87 | 8.2 | 50 | 8.2
$1.0\times 10^{-8}$ | 28 | 8.8 | 34 | 10.4 | 154 | 5.7 | 82 | 5.5
JR+QN | $1.0\times 10^{-2}$ | 6 | 6.2 | 6 | 6.0 | 8 | 5.0 | 7 | 5.3
$1.0\times 10^{-4}$ | 11 | 4.9 | 11 | 5.3 | 13 | 4.1 | 14 | 4.1
$1.0\times 10^{-6}$ | 19 | 3.9 | 13 | 5.2 | 89 | 1.9 | 46 | 2.5
$1.0\times 10^{-8}$ | 21 | 3.6 | 36 | 3.8 | 112 | 1.7 | 46 | 2.5

  * a

$C_{\mathrm{min}}$: threshold for terminating the selection of individual
configuration state functions.

  * b

$A_{\mathrm{min}}$: threshold for terminating the active-active orbital
rotations.

  * c

n.c.: not converged within 200 macroiterations.

---
Figure 3: The convergence trends of hexacene calculations using QN and JR+QN
algorithms, with C${}_{\mathrm{min}}=5.0\times 10^{-5}$ and
$A_{\mathrm{min}}=1.0\times 10^{-8}$ settings, using LMOs as initial guesses.

###  3.2 Adiabatic singlet-triplet gaps of Polyacenes

The singlet-triplet (S-T) gaps of polyacenes have been studied by several
groups using various methods tailored for large active spaces51, 14, 52, 28.
To compare with such calculations, the adiabatic S-T gaps are recalculated
with iCISCF and iCISCF(2) using the cc-pVDZ basis. The geometries of
naphthalene and C4n+2H2n+4 ($n$=3-10) are taken from Refs. 52 and 28,
respectively. The distributions of the orbitals in the irreducible
representations (irrep) of $D_{2h}$ are given in Table 2. Interestingly, the
increments of inactive orbitals in $B_{1g}$ and $B_{2u}$ are different for
even and old n-acenes. As a result, the HOMO and LUMO of even and old n-acenes
have different symmetries. The calculated energies of hexacene are given in
Table 3. To make a close comparison with ASCISCF running with 106 SDs28, we
also report the interpolated iCISCF and iCISCF(2) energies. It can be seen
that iCISCF produces slightly lower absolute energies than ASCISCF does for
both singlet and triplet state employing the same $P_{m}$ size. The deviation
between the two methods at the variational level is about 3.4 kcal/mol for
${}^{1}A_{g}$ state and 1.74 kcal/mol for ${}^{3}B_{u}$ state. Presumably,
this small discrepancy stems from spin contamination in ASCISCF, which is
partly removed by the PT2 correction, thereby leading to a closer agreement
between ASCISCF(2) and iCISCF(2) for the S-T gap, within 1.0 kcal/mol. The
extrapolated S-T gaps by iCISCF(2) and ASCISCF(2) agree within 1.0 kcal/mol as
well.

Table 2:  Distributions of core and active orbitals of n-acenes in irreducible representations of $D_{2h}$ $n$ | active space |  | $A_{g}$ | $B_{3g}$ | $B_{2g}$ | $B_{1g}$ | $A_{u}$ | $B_{3u}$ | $B_{2u}$ | $B_{1u}$
---|---|---|---|---|---|---|---|---|---|---
5 | (22, 22) | core | 18 | 0 | 0 | 13 | 0 | 15 | 16 | 0
active | 0 | 5 | 6 | 0 | 5 | 0 | 0 | 6
6 | (26, 26) | core | 21 | 0 | 0 | 16 | 0 | 17 | 19 | 0
active | 0 | 6 | 7 | 0 | 6 | 0 | 0 | 7
7 | (30, 30) | core | 24 | 0 | 0 | 18 | 0 | 20 | 22 | 0
active | 0 | 7 | 8 | 0 | 7 | 0 | 0 | 8
Table 3: iCISCF and iCISCF(2) energies (+994.0 $E_{\mathrm{h}}$) of hexacene
with CAS(26, 26)/cc-pVDZ

${}^{1}A_{g}$ ${}^{3}B_{3u}$ Gap (kcal/mol) Cmin $N_{\mathrm{CSF}}$ a
$N_{\mathrm{SD}}$ b $E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$
$N_{\mathrm{CSF}}$ a $N_{\mathrm{SD}}$ b $E_{\mathrm{iCISCF}}$
$E_{\mathrm{iCISCF(2)}}$ $E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$
$1.0\times 10^{-4}$ 80618 305351 -0.29320 -0.30699 46057 200173 -0.25549
-0.27317 23.66 21.22 $7.5\times 10^{-5}$ 128424 495560 -0.29675 -0.30901
197610 343609 -0.26043 -0.27586 22.79 20.80 $5.0\times 10^{-5}$ 227406 900944
-0.30084 -0.31119 346020 612068 -0.26517 -0.27837 22.38 20.59 $2.5\times
10^{-5}$ 577599 2405533 -0.30641 -0.31407 875294 1593875 -0.27174 -0.28185
21.76 20.22 Estimated c - 1000000 -0.30108 -0.31126 - 1000000 -0.26847
-0.28012 20.46 19.54 Extrapolated d $3.7\times 10^{12}$ $2.7\times 10^{13}$ -
-0.32300 $8.4\times 10^{12}$ $2.3\times 10^{13}$ - -0.29345 - 18.53 ASCISCF(2)
e - 1000000 -0.29562 -0.30753 - 1000000 -0.26570 -0.27771 18.78 18.72
Extrapolated - - - - - - - - - 19.4

  * a

Number of selected configuration state functions.

  * b

Number of Slater determinants corresponding to $N_{\mathrm{CSF}}$.

  * c

Estimated by linear fit of $E_{\mathrm{iCISCF}}$ or $E_{\mathrm{iCISCF(2)}}$
as function of $log_{10}N_{\mathrm{SD}}$.

  * d

$N_{\mathrm{CSF}}$/$N_{\mathrm{SD}}$: total number of CSFs/SDs in the symmetry
adapted CAS.

  * e

Ref. 28.

The absolute energies of n-acenes computed at iCISCF and iCISCF(2) level have
been provided in the Supporting information (SI). The extrapolated iCISCF(2)
S-T splittings of n-acenes are plotted in Fig. 4. As can be seen, the results
with the cc-pVDZ and cc-pVTZ basis sets are very close, reflecting marginal
basis set incompleteness errors. Except for noacene, the iCISCF(2)/cc-pVDZ
results are very similar to those by ASCISCF(2)/cc-pVDZ28. The obvious
deviations of iCISCF(2) (and ASCISCF(2)) from ACI(2)-DSRG-MRPT2/cc-pVTZ
(adaptive CI with second-order perturbative multireference-driven similarity
renormalization group)52 should then be ascribed to dynamic correlation from
the $Q$ space that is accounted for by DSRG-MRPT2. For more details, the
differences $\Delta E_{\mathrm{ST}}$ between the actually calculated and the
extrapolated S-T gaps are further plotted in Fig. 5, which shows that $\Delta
E_{\mathrm{ST}}$ becomes larger as the system size increases but can be
reduced by decreasing $C_{\mathrm{min}}$. However, even for the smallest
$C_{\mathrm{min}}$ (i.e., $2.5\times 10^{-5}$) used here, $\Delta
E_{\mathrm{ST}}$ for decacene is still as large as 3.0 kcal/mol, which is to
be compared with the extrapolated value of S-T gap, 5.2 kcal/mol. Therefore,
predicting accurately the S-T gaps of decacene and beyond remains a challenge
if the extrapolation step is not taken.

---
Figure 4: Extrapolated iCISCF(2) adiabatic S-T gaps of polyacenes with cc-pVDZ
and cc-pVTZ, to compare with those by ASCISCF(2)/cc-pVDZ28 and ACI(2)-DSRG-
MRPT2/cc-pVTZ52.
---
Figure 5: Differences ($\Delta E_{\mathrm{ST}}$) between the calculated and
extrapolated adiabatic S-T gaps of polyacenes by iCISCF(2)/cc-pVDZ.

###  3.3 Fe-porphyrin complex

Fe(II) porphyrin (FeP) complexes are another set of classic test systems for
strongly correlated methods53, 24, 54, 55, 27, due to both large active spaces
and small energy gaps between different spin states. The free-base porphyrin
already has 24 valence $\pi$ orbitals. Therefore, considering the eight
electrons and five $3d$ orbitals of Fe, a minimal active space of FeP would be
as large as CAS(32, 29). Even so, the correct spin ordering still cannot be
reproduced24, 27. That ${}^{3}B_{1g}$ is lower than ${}^{5}A_{g}$ was only
found by going to CAS(44, 44)24, which consists of additional Fe
$4d4p_{x}4p_{y}$ and $4\times\mathrm{N}$ $2p_{x}2p_{y}$. Realizing that
CAS(44, 44) is not really rational, Levine et al.27 proposed to use CAS(40,
42), which includes Fe $4s4p4d$ and four $\sigma$ lone pair orbitals of N
coordinated to Fe on top of CAS(32, 29). To compare directly with these
results, the iCISCF/cc-pVDZ calculations were performed with the
${}^{3}B_{1g}$ geometry56 for both ${}^{3}B_{1g}$ and ${}^{5}A_{g}$. The
distributions of the orbitals of FeP in the irreps of $D_{2h}$ are given in
Table 4. Upon convergence after 140 macroiterations of iCISCF-CAS(32, 29) with
$C_{\mathrm{min}}=1.0\times 10^{-4}$, close inspections reveal that 24 of the
29 active orbitals are indeed the valence $\pi$ orbitals of porphyrin, but the
remaining 5 are not all the expected Fe 3$d$. It can be seen clearly from Fig.
6 that in the case of ${}^{5}A_{g}$, the nearly doubly occupied 3$d_{z^{2}}$
orbital is pushed out of the active space, whereas a doubly occupied ligand
orbital is rotated into the active space to correlate the singly occupied
3$d_{x^{2}-y^{2}}$ orbital. For the ${}^{3}B_{1g}$ state, the converged active
space does not contain the 3$d_{x^{2}-y^{2}}$ orbital, which is pushed to the
virtual space. Instead, the 4$d_{z^{2}}$ orbital enters the active space to
incorporate the 3$d_{z^{2}}$ orbital. The five 3$d$ orbitals can be kept
inside the active space only by freezing the 3$d_{z^{2}}$ (3$d_{x^{2}-y^{2}}$)
orbital in the calculations of ${}^{5}A_{g}$ (${}^{3}B_{1g}$).

Table 4:  Distributions of orbitals of FeP in irreducible representations of $D_{2h}$ active space |  | $A_{g}$ | $B_{3g}$ | $B_{2g}$ | $B_{1g}$ | $A_{u}$ | $B_{3u}$ | $B_{2u}$ | $B_{1u}$
---|---|---|---|---|---|---|---|---|---
(32, 29) | core | 23 | 0 | 0 | 14 | 0 | 19 | 19 | 2
active | 2 | 7 | 7 | 1 | 5 | 0 | 0 | 7
(40, 42) | core | 21 | 0 | 0 | 14 | 0 | 18 | 18 | 2
active | 7 | 8 | 8 | 2 | 5 | 2 | 2 | 8
(44, 44) | core | 21 | 0 | 0 | 14 | 0 | 17 | 17 | 2
active | 6 | 8 | 8 | 4 | 5 | 3 | 3 | 7
---
Figure 6: The CAS(32, 29) active orbitals of (a) initial guess, (b) converged
${}^{5}A_{g}$, and (c) converged ${}^{3}B_{1g}$. The 24 valence $\pi$ orbitals
of porphyrin are not shown.

The iCISCF-CAS(32, 29) energies are documented in Table 5. As the previous
studies24, 27, the ${}^{5}A_{g}$ state is predicted to be the ground state.
The extrapolated iCISCF(2) results are in good agreement with those from
HCISCF(2)24 for both the total and relative energies. In contrast, the
extrapolated energy gap by ASCISCF(2)27 is just half that by iCISCF(2) or
HCISCF(2). This is mainly because the ASCISCF(2) energy for ${}^{3}B_{1g}$ is
too low as compared with the iCISCF(2) one. Freezing the 3$d_{x^{2}-y^{2}}$
(3$d_{z^{2}}$) orbital inside the active space increases the energy of
${}^{3}B_{1g}$ (${}^{5}A_{g}$) by 3.6 (0.3) kcal/mol, thereby increasing the
gap by 3.3 kcal/mol.

Table 5: iCISCF-CAS(32, 29) energies (+2245.0 $E_{\mathrm{h}}$) of FeP a

${}^{3}B_{1g}$ ${}^{5}A_{g}$ Gap (kcal/mol) $C_{\mathrm{min}}$
$N_{\mathrm{CSF}}$ $N_{\mathrm{SD}}$ $E_{\mathrm{iCISCF}}$
$E_{\mathrm{iCISCF(2)}}$ $N_{\mathrm{CSF}}$ $N_{\mathrm{SD}}$
$E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$ $E_{\mathrm{iCISCF}}$
$E_{\mathrm{iCISCF(2)}}$ frozen b $1.0\times 10^{-4}$ 103238 186387 0.02457
0.00983 108153 161556 -0.00508 -0.02113 -18.61 -19.43 $7.5\times 10^{-5}$
154932 284279 0.02144 0.00850 165285 249441 -0.00860 -0.02283 -18.85 -19.66
$5.0\times 10^{-5}$ 271010 507312 0.01768 0.00691 296625 454698 -0.01231
-0.02413 -18.82 -19.48 $2.5\times 10^{-5}$ 691268 1331168 0.01258 0.00474
774101 1211597 -0.01777 -0.02639 -19.04 -19.53 Extrapolated $1.7\times
10^{14}$ $5.0\times 10^{14}$ - -0.0010 $1.7\times 10^{14}$ $3.4\times 10^{14}$
- -0.0323 - -19.64 fully optimized $1.0\times 10^{-4}$ 106323 193856 0.02069
0.00467 109165 162983 -0.00533 -0.02149 -16.33 -16.41 $7.5\times 10^{-5}$
165932 307657 0.01723 0.00323 165971 250522 -0.00858 -0.02284 -16.19 -16.36
$5.0\times 10^{-5}$ 296081 558728 0.01314 0.00152 302190 463045 -0.01273
-0.02455 -16.23 -16.36 $2.5\times 10^{-5}$ 758072 1469095 0.00763 -0.00077
785566 1228691 -0.01818 -0.02680 -16.19 -16.33 Extrapolated $1.7\times
10^{14}$ $5.0\times 10^{14}$ - -0.0068 $1.7\times 10^{14}$ $3.4\times 10^{14}$
- -0.0328 - -16.36 HCISCF(2) c - 533623 0.0224 - - 379536 0.0020 - -12.80 -
Extrapolated - - - -0.0049 - - - -0.0314 - -16.63 ASCISCF(2) d - 500000
-0.0022 - - 500000 -0.0191 - -10.60 - Extrapolated - - - -0.0186 - - - -0.0324
- -8.63

  * a

See Table 3 for additional explanations.

  * b

The 3$d_{z^{2}}$ and 3$d_{x^{2}-y^{2}}$ orbitals are kept frozen in the active
space when calculating ${}^{5}A_{g}$ and ${}^{3}B_{1g}$, respectively.

  * c

Ref. 24

  * d

Ref. 27

To compare with the HCISCF-CAS(44, 44) results24, the iCISCF-CAS(44, 44)
calculations were also performed. Except for the 24 $\pi$ orbitals of
porphyrin, the remaining active orbitals are plotted in Fig. S1 in the SI. It
turns out that the Fe 4$d_{x^{2}-y^{2}}$ orbital is pushed out of the active
space (in line with the previous observation27), while the 27${}^{th}A_{g}$
orbital, consisting mainly of Fe $4s$ and N $3p$, is pulled into the active
space. The former can be viewed as the “anti-bonding” orbital of
3$d_{x^{2}-y^{2}}$ whereas the latter can be viewed as the “anti-bonding”
orbital of the 22${}^{th}A_{g}$, a lone pair orbital of four N atoms
coordinate to Fe. It is better to put both 4$d_{x^{2}-y^{2}}$ and $4s$ into
the active space as in the case of CAS(40, 42). If only one of them is to be
included, the Fe $4s$ instead of 4$d_{x^{2}-y^{2}}$ orbital is favored for
lower energies. Again, we performed calculations by freezing the
4$d_{x^{2}-y^{2}}$ orbital in the active space, see Fig. S2 for the converged
4$d_{x^{2}-y^{2}}$ and 4$d_{z^{2}}$ orbitals. The latter is slightly mixed
with the 4$s$ orbitals after optimization.

Table 6: iCISCF-CAS(44, 44) energies (+2245.0 $E_{\mathrm{h}}$) of FeP a

${}^{3}B_{1g}$ ${}^{5}A_{g}$ Gap (kcal/mol) Cmin $N_{\mathrm{CSF}}$
$N_{\mathrm{SD}}$ $E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$
$N_{\mathrm{CSF}}$ $N_{\mathrm{SD}}$ $E_{\mathrm{iCISCF}}$
$E_{\mathrm{iCISCF(2)}}$ $E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$ frozen
b $1.0\times 10^{-4}$ 197817 366295 -0.133538 -0.169007 179562 264701
-0.135383 -0.167162 -1.16 1.16 $7.5\times 10^{-5}$ 309455 582411 -0.139264
-0.171517 280883 419065 -0.140490 -0.169338 -0.77 1.37 $5.0\times 10^{-5}$
574704 1103405 -0.146524 -0.174667 527665 800161 -0.147107 -0.172201 -0.37
1.55 $2.5\times 10^{-5}$ 1629683 3222519 -0.156921 -0.179173 1510972 2345179
-0.156565 -0.176247 0.22 1.84 Extrapolated $1.2\times 10^{23}$ $5.1\times
10^{23}$ - -0.1963 $1.4\times 10^{23}$ $3.9\times 10^{23}$ - -0.1911 - 3.28
fully optimized $1.0\times 10^{-4}$ 199609 372869 -0.15415 -0.19010 180638
267708 -0.14558 -0.17787 5.38 7.67 $7.5\times 10^{-5}$ 317682 603672 -0.16021
-0.19264 281562 422181 -0.15082 -0.18003 5.90 7.91 $5.0\times 10^{-5}$ 587999
1138209 -0.16752 -0.19583 543338 828687 -0.15784 -0.18304 6.07 8.02 $2.5\times
10^{-5}$ 1648308 3278017 -0.17797 -0.20037 1541226 2402169 -0.16738 -0.18716
6.64 8.29 Extrapolated $1.2\times 10^{23}$ $5.1\times 10^{23}$ - -0.2174
$1.4\times 10^{23}$ $3.9\times 10^{23}$ - -0.2019 - 9.73 HCISCF(2) c - 2133424
-0.1567 - - 1450271 -0.1457 - 6.90 - Extrapolated - - - -0.1996 - - - -0.1965
- 1.93 ASCISCF(2) d - 500000 -0.2044 - - 500000 -0.2315 - -17.01 -
Extrapolated - - - -0.25617 - - - -0.28688 - -19.27

  * a

See Table 3 for additional explanations.

  * b

The 4$d_{x^{2}-y^{2}}$ orbitals are kept frozen in the active space.

  * c

Ref. 24

  * d

Ref. 27

The iCISCF-CAS(44, 44) energies are documented in Table 6. In all calculations
with different $C_{\mathrm{min}}$, the ${}^{3}B_{1g}$ state is predicted to be
the ground state. However, the presently extrapolated iCISCF(2) energy gap
(9.7 kcal/mol) between ${}^{5}A_{g}$ and ${}^{3}B_{1g}$ is significantly
larger than that (1.9 kcal/mol) by HCISCF(2)24. The latter is even smaller
than the gap (3.3 kcal/mol) obtained by iCISCF(2) with 4$d_{x^{2}-y^{2}}$
frozen. Close inspections reveal that, although the HCISCF gap (6.9 kcal/mol)
is very much the same as the iCISCF one (6.6 kcal/mol) with
$C_{\mathrm{min}}=2.5\times 10^{-5}$, the HCISCF energies for ${}^{5}A_{g}$
and ${}^{3}B_{1g}$ are actually close to the iCISCF ones with
$C_{\mathrm{min}}=1.0\times 10^{-4}$, given dramatic differences in the
numbers of SDs. It can therefore be concluded that the HCISCF calculations24
are far from convergence. Note also that ASCISCF-CAS(44, 44)27 predicts much
lower variational energies but still incorrect spin ordering, which was
ascribed to spatial symmetry breaking. To investigate such symmetry-breaking
effects, we repeated the calculations using C2v instead of D2h. Indeed,
localized symmetry-breaking active orbitals are observed (cf. Fig. S3),
leading to lower energies than those with D2h.

Some representative active orbitals from iCISCF-CAS(40, 42) calculations are
plotted in Fig. S4. As can be seen from the energetics in Table 7, both iCISCF
and iCISCF(2) predict the correct ground state with different
$C_{\mathrm{min}}$. And in this case, there exist perfect agreements between
iCISCF and ASCISCF27 for both the extrapolated total and relative energies.
Yet, rather unexpectedly, the iCISCF-CAS(40, 42) energies are even lower than
the iCISCF-CAS(44, 44) ones (cf. Table 6), reflecting the importance of
rational construction of active spaces (for which the $\mathbbm{i}$CAS method
is a good choice).

Table 7: iCISCF-CAS(40, 42) energies (+2245.0 $E_{\mathrm{h}}$) of FeP a

${}^{3}B_{1g}$ ${}^{5}A_{g}$ Gap (kcal/mol) Cmin $N_{\mathrm{CSF}}$
$N_{\mathrm{SD}}$a $E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$
$N_{\mathrm{CSF}}$ $N_{\mathrm{SD}}$a $E_{\mathrm{iCISCF}}$
$E_{\mathrm{iCISCF(2)}}$ $E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$
$1.0\times 10^{-4}$ 199810 370332 -0.15477 -0.19185 176251 259571 -0.15428
-0.18806 0.31 2.38 $7.5\times 10^{-5}$ 314006 591402 -0.16086 -0.19468 279206
416350 -0.16002 -0.19072 0.52 2.48 $5.0\times 10^{-5}$ 595172 1143717 -0.16882
-0.19838 535154 811146 -0.16738 -0.19412 0.91 2.67 $2.5\times 10^{-5}$ 1686557
3331667 -0.18018 -0.20369 1566899 2429062 -0.17800 -0.19907 1.37 2.90
Estimated - 500000 -0.1587 - - 500000 -0.1617 - -1.87 - Extrapolated
$7.3\times 10^{21}$ $3.0\times 10^{22}$ - -0.2242 $8.5\times 10^{21}$
$2.2\times 10^{22}$ - -0.2173 - 4.34 ASCISCF(2) b - 500000 -0.1671 - - 500000
-0.1699 - -1.76 - Extrapolated - - - -0.2208 - - - -0.2137 - 4.46

  * a

See Table 3 for additional explanations.

  * b

Ref. 27

###  3.4 A rational design of active space for Fe-porphyrin

Although the FeP complexes have been studied by many groups utilizing
different active spaces, there are not many works explaining how their active
spaces for FeP are designed. It has been proved that the $\mathbbm{i}$CAS
method developed by some of us is a useful tool to design and generate initial
guess MOs for MCSCF calculations.35 In this section, we would like to design a
few active spaces from small to large utilizing the $\mathbbm{i}$CAS
algorithm. All initial guess MOs in this section are generated based on the
ROHF wave function of ${}^{5}A_{g}$ state.

One of the smallest active space to describe various spin states of FeP is
CAS(6, 5), which contains only the 3$d$ orbitals and electrons of Fe2+. It has
been proved that this small active space is not able to describe the strong
interaction between Fe2+ and porphyrin anion.57 More orbitals and electrons
have to be included in the active space. By comparing the results of CAS(32,
29) and CAS(40, 42) in Sec. 3.3, it is evident that the $\pi$ orbitals of
porphyrin do not influence the spin splitting of Fe2+ as significantly as the
ligand orbitals and 4$s$4$p$4$d$ orbitals of Fe do. Usually, the ligand
orbitals interact with the metal centre strongly. Thus, in the first step, the
five 3$d$ orbitals of Fe, the 2$s$, 2$p_{x}$, 2$p_{y}$ orbitals of the four N
atoms are used to generate pre-LMOs and the subsequent LMOs guesses for iCISCF
calculations. Besides the five 3$d$ orbitals, 12 LMOs, eight doubly occupied
MOs (DOMOs) and four virtual MOs, are generated by the $\mathbbm{i}$CAS as
active MOs. In fact, not all the 12 LMOs are spatially close to the Fe2+. By
checking their coefficients, we found that only four doubly occupied LMOs are
mixed with the atomic orbitals of Fe. The four ligand LMOs are exhibited in
Fig. 7, together with their significant MOs coefficients. If only the four
LMOs and the Fe 3$d$ orbitals are considered as active MOs, the active space
will be CAS(14,9). However, such an active space only contains DOMOs and
SOMOs, and does not have any virtual orbitals to relax the electrons in the
active space. By analyzing the coefficients of the four LMOs, another 3 atomic
orbitals are added, the 4$s$, 4$p_{x}$ and 4$p_{y}$ orbtial of Fe2+. Thus, a
CAS(14,12) active space is generated by the $\mathbbm{i}$CAS method. Moreover,
to account for the doubles-shell effects, the 4$d$ orbitals of Fe atom could
be considered, resulting in iCISCF-CAS(14,17) calculations.

Although the $\pi$ orbitals of porphyrin are not spatially adjacent to Fe2+,
some studies showed that they do differentially stabilize the triplet states
over quintet states of FeP.58 For the low-lying excited states of FeP, the 24
$\pi$ orbitals are distributed to four different irreps of $D_{2h}$ symmetry,
six $B_{3g}$ MOs (three DOMOs), six $B_{2g}$ MOs (three DOMOs), five $A_{u}$
MOs (two DOMOs), and seven $B_{1u}$ MOs (five DOMOs). To consider the
influence of the $\pi$ orbitals, it may not be necessary to include all 24
$\pi$ MOs in the active space. In principle, the $\pi$ orbitals with close
energies as the HOMO and LUMO, should interact with the Fe2+ more strongly.
Thus, the 24 $\pi$ LMOs, constructed from the pre-LMOs, are recanonicalized
based on their Fock matrix (the Fock matrix of DOMOs or virtual MOs are
diagonalized separately), which is called regional LMOs. After the
canonicalization, the eigenvalues of the Fock matrix within each subspace can
be considered as approximate orbital energies. The 24 $\pi$ regional LMOs with
their irreps and energies are given in Fig. 8. Thus, for the first step, ten
orbitals (from the 9th to 18th MOs in Fig. 8) are added to the CAS(14,17)
active space, including two $B_{3g}$ MOs (one DOMO), two $B_{2g}$ MOs (one
DOMO), two $A_{u}$ MOs (one DOMO), and four $B_{1u}$ MOs (two DOMOs).
Moreover, we found that the orbital dominated by the 4$p_{z}$ orbital of Fe
with $B_{1u}$ symmetry has lower orbital energy than the frontier virtual
$\pi$ orbital in Fig. 8 (the 14th MO), which could be included in the active
space as well. Thus, two sets of guess MOs with and without the 4$p_{z}$
orbital of Fe centre, CAS(14,28) and CAS(14,27), are generated as well. To
consider the influence of more $\pi$ orbitals, another seven orbitals (from
the 5th to 8th and from 19th to 21th in Fig. 8) could be taken into account.
Thus, on top of CAS(24,28), eight more electrons and seven more active MOs are
added to the active space, resulting in iCISCF-CAS(32,35) calculations.
Finally, by adding the remaining valence $\pi$ orbital of porphyrin to the
active space, we arrive at the active space proposed by Levine and coworkers,
CAS(40,42). The distribution of the core and active orbitals of calculations
with the above mentioned seven different active spaces are given in Table 8.

Using all active spaces in Table 8, the ${}^{3}A_{2g}$ (${}^{3}B_{1g}$ in
$D_{2h}$), ${}^{3}E_{g}$ (${}^{3}B_{2g}$+${}^{3}B_{3g}$ in $D_{2h}$) and
${}^{5}A_{1g}$ (${}^{5}A_{g}$ in $D_{2h}$) state of FeP are computed. Note
that for the Cartesian coordinates offered by Groenhof and coworkers, the
$D_{4h}$ symmetry is not applied during the geometry optimization56. Thus, the
degenerate ${}^{3}E_{g}$ state is split to two states in $D_{2h}$ symmetry,
${}^{3}B_{2g}$ and ${}^{3}B_{3g}$. The three triplet and one quintet states
computed using various active space are given in Table 9. The results show
that with the two smallest active space, CAS(6,5) and CAS(14,12), the
${}^{5}A_{g}$ state is predicted to be the ground state. Thus, to include only
the ligand orbitals in the active space, the triplet states are not stabilized
over the quintet state very much. However, by taking into account the double-
shell effects (iCISCF-CAS(14,17) results in Table 8), the
${}^{5}A_{g}-^{3}B_{1g}$ splitting amounts to -0.55 kcal/mol. Compared to that
with CAS(14,12), the ${}^{5}A_{g}-^{3}B_{1g}$ gap increases 13.90 kcal/mol.
Moreover, the excitation energies of ${}^{3}B_{2g}$ and ${}^{3}B_{3g}$ state
with respect to the ${}^{3}B_{1g}$ state are close to that delivered by the
iCISCF-CAS(40,42) calculations as well. Thus, it is the five 4$d$ orbitals
destabilizing the quintet states with respect to the ${}^{3}B_{1g}$ state more
significantly.

To include some $\pi$ orbitals of porphyrin or the 4$p_{z}$ orbital of Fe in
the active space, iCISCF results using four different active spaces are
reported in Table 8 as well. Compared to CAS(24,27), the relative energies
computed by iCISCF-CAS(24,28) are in good agreement with the most accurate
ones, iCISCF-CAS(40,42). The deviations of excitation energies computed by the
extrapolated iCISCF results are less than 1.0 kcal/mol. Further increase the
active space to CAS(32,35), more accurate ${}^{5}A_{g}-^{3}B_{1g}$ gaps are
predicted. However, the relative energies of ${}^{3}B_{2g}$ and ${}^{3}B_{3g}$
states with respect to the ${}^{3}B_{1g}$ state are getting worse.
Nevertheless, the absolute deviations of extrapolated results are still less
than 1.0 kcal/mol with respect to the iCISCF-CAS(40,42) results.

The results in Table 8 show that the 4$d$ orbitals are the most important to
differentially stabilize the triplet states over the ${}^{5}A_{g}$ state.
Although the $\pi$ orbitals of porphyrin are essential to deliver correct spin
ordering, it may not be necessary to take all the 24 $\pi$ orbitals into
consideration. The active space CAS(24,28) can already predict qualitatively
correct results for the four states, which could be used as starting points
for the subsequent dynamic calculations.

Table 8:  Distributions of orbitals of FeP in irreducible representations of $D_{2h}$ for some iCISCF calculations active space |  | $A_{g}$ | $B_{3g}$ | $B_{2g}$ | $B_{1g}$ | $A_{u}$ | $B_{3u}$ | $B_{2u}$ | $B_{1u}$
---|---|---|---|---|---|---|---|---|---
( 6, 5) | core | 23 | 3 | 3 | 14 | 2 | 19 | 19 | 7
active | 2 | 1 | 1 | 1 | 0 | 0 | 0 | 0
(14, 12) | core | 21 | 3 | 3 | 14 | 2 | 18 | 18 | 7
active | 5 | 1 | 1 | 1 | 0 | 2 | 2 | 0
(14, 17) | core | 21 | 3 | 3 | 14 | 2 | 18 | 18 | 7
active | 7 | 2 | 2 | 2 | 0 | 2 | 2 | 0
(24, 27) | core | 21 | 2 | 2 | 14 | 1 | 18 | 18 | 5
active | 7 | 4 | 4 | 2 | 2 | 2 | 2 | 4
(24, 28) | core | 21 | 2 | 2 | 14 | 1 | 18 | 18 | 5
active | 7 | 4 | 4 | 2 | 2 | 2 | 2 | 5
(32, 35) | core | 21 | 1 | 1 | 14 | 0 | 18 | 18 | 4
active | 7 | 6 | 6 | 2 | 4 | 2 | 2 | 6
(40, 42) | core | 21 | 0 | 0 | 14 | 0 | 18 | 18 | 2
active | 7 | 8 | 8 | 2 | 5 | 2 | 2 | 8
---
Figure 7: The four ligand orbitals close to the Fe centre. Their significant
MO coefficients are listed as well.
---
Figure 8: The 24 $\pi$ regional LMOs of porphyrin constructed by the
$\mathbbm{i}$CAS method, ordered by their orbital energies (in
$E_{\mathrm{h}}$). Table 9:  The energies of lowest ${}^{3}B_{1g}$ state of
FeP (in $E_{\mathrm{h}}$) using seven different active spaces, the relative
energies of ${}^{3}B_{2g}$, ${}^{3}B_{3g}$, and ${}^{5}A_{g}$ state w.r.t the
${}^{3}B_{1g}$ state are given in kcal/mol

$E_{\mathrm{iCISCF}}$ $E_{\mathrm{iCISCF(2)}}$ Cmin ${}^{3}B_{1g}$
${}^{3}B_{2g}$ ${}^{3}B_{3g}$ ${}^{5}A_{g}$ ${}^{3}B_{1g}$ ${}^{3}B_{2g}$
${}^{3}B_{3g}$ ${}^{5}A_{g}$ CAS(6,5) $1.0\times 10^{-4}$ -2244.72472 4.33
4.19 -17.55 -2244.72472 4.33 4.19 -17.55 $0.0$ -2244.72472 4.33 4.19 -17.55
-2244.72472 4.33 4.19 -17.55 CAS(14,12) $1.0\times 10^{-4}$ -2244.76054 -0.06
-0.24 -14.45 -2244.76057 -0.08 -0.26 -14.45 $0.0$ -2244.76057 -0.08 -0.26
-14.45 -2244.76057 -0.08 -0.26 -14.45 CAS(14,17) $1.0\times 10^{-4}$
-2244.86467 1.94 1.68 -0.74 -2244.86554 1.80 1.54 -0.58 $7.5\times 10^{-5}$
-2244.86487 1.90 1.63 -0.67 -2244.86556 1.80 1.54 -0.57 $5.0\times 10^{-5}$
-2244.86507 1.86 1.60 -0.65 -2244.86558 1.80 1.53 -0.57 $2.5\times 10^{-5}$
-2244.86532 1.84 1.58 -0.60 -2244.86559 1.80 1.53 -0.57 Extrapolated - - - -
-2244.86562 1.79 1.54 -0.55 CAS(24,27) $1.0\times 10^{-4}$ -2244.98960 2.14
1.76 1.79 -2244.99602 1.95 1.56 2.79 $7.5\times 10^{-5}$ -2244.99090 2.32 1.88
1.98 -2244.99634 2.02 1.61 2.86 $5.0\times 10^{-5}$ -2244.99233 2.16 1.76 2.17
-2244.99668 1.95 1.56 2.91 $2.5\times 10^{-5}$ -2244.99420 2.06 1.65 2.45
-2244.99706 1.89 1.51 2.97 Extrapolated - - - - -2244.99792 1.81 1.44 3.17
CAS(24,28) $1.0\times 10^{-4}$ -2244.99504 2.57 2.08 1.85 -2245.00255 2.36
1.86 2.99 $7.5\times 10^{-5}$ -2244.99656 2.78 2.27 2.09 -2245.00297 2.45 1.94
3.07 $5.0\times 10^{-5}$ -2244.99827 2.62 2.15 2.32 -2245.00343 2.37 1.89 3.15
$2.5\times 10^{-5}$ -2245.00052 2.50 2.02 2.63 -2245.00395 2.31 1.83 3.22
Extrapolated - - - - -2245.00516 2.21 1.75 3.47 CAS(32,35) $1.0\times 10^{-4}$
-2245.07673 3.17 2.83 1.07 -2245.09709 2.65 2.21 2.73 $7.5\times 10^{-5}$
-2245.08070 3.14 2.78 1.57 -2245.09877 2.66 2.20 3.00 $5.0\times 10^{-5}$
-2245.08519 2.86 2.57 1.86 -2245.10065 2.49 2.03 3.15 $2.5\times 10^{-5}$
-2245.09145 2.63 2.29 2.24 -2245.10316 2.28 1.80 3.30 Extrapolated - - - -
-2245.11142 1.71 1.02 4.26 CAS(40,42) $1.0\times 10^{-4}$ -2245.15477 4.21
3.94 0.31 -2245.19185 3.11 2.78 2.38 $7.5\times 10^{-5}$ -2245.16086 3.60 3.29
0.52 -2245.19468 3.05 2.67 2.48 $5.0\times 10^{-5}$ -2245.16882 3.48 3.14 0.91
-2245.19838 2.88 2.50 2.67 $2.5\times 10^{-5}$ -2245.18018 3.29 2.95 1.37
-2245.20369 2.79 2.38 2.90 Extrapolated - - - - -2245.22419 2.57 1.99 4.34

  * a

See Table 3 for additional explanations.

##  4 Conclusions and Outlook

A nearly exact CASSCF approach, iCISCF, has been proposed to treat strongly
correlated systems that are beyond the capability of CASSCF. It combines the
selected iCI for the active space solver and the hybrid JR+QN algorithm for
orbital optimization. Upon convergence, an inner-space second-order
perturbation step (i.e., iCISCF(2)) can further be taken to improve the energy
estimate. The current implementation of iCISCF/iCISCF(2) can handle ca. 50
active orbitals in conjunction with more than 1000 basis functions, as shown
by the showcases of polyacenes and Fe(II)-porphyrin. Further developments of
iCISCF will include (1) the use of $\mathbbm{i}$CAS35 for automatic
construction and localization of CAS, (2) extrapolation/interpolation schemes
for speedup of orbital optimization, (3) analytic energy gradients (which are
operationally precisely the same as those of CASSCF), and (4) stochastic
treatment of dynamic correlation from the complementary space $Q$.

## Acknowledgement

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 21833001 and 21973054), Mountain Tai Climb Program of Shandong
Province, and Key-Area Research and Development Program of Guangdong Province
(Grant No. 2020B0101350001). YG was further supported by the Qilu Young
Scholar Program of Shandong University.

## Data Availability Statement

The data that supports the findings of this study are available within the
article.

##  Appendix A Orbital gradients and diagonal Hessian

By virtue of the basic commutator

$\displaystyle[\hat{E}_{pq},\hat{E}_{rs}]=\delta_{qr}\hat{E}_{ps}-\delta_{ps}\hat{E}_{rq},$
(A.1)

the following commutators can be derived

$\displaystyle[\hat{E}_{pq},\hat{H}]$ $\displaystyle=$
$\displaystyle(1-P_{pq}c.c.)\sum_{\sigma}a_{p\sigma}^{\dagger}[a_{q\sigma},H]$
(A.2) $\displaystyle=$
$\displaystyle\\{\hat{E}_{pr}h_{rq}+\hat{e}_{pr,st}g_{rq,st}\\}-\\{h_{pr}\hat{E}_{rq}+g_{pr,st}\hat{e}_{rq,st}\\},$
(A.3) $\displaystyle[\hat{E}_{pq}^{-},\hat{H}]$ $\displaystyle=$
$\displaystyle[\hat{E}_{pq},\hat{H}]+c.c.$ $\displaystyle=$
$\displaystyle\\{(\hat{E}_{pr}+\hat{E}_{rp})h_{rq}+(\hat{e}_{pr,st}+\hat{e}_{rp,st})g_{rq,st}\\}$
$\displaystyle-$
$\displaystyle\\{h_{pr}(\hat{E}_{rq}+\hat{E}_{qr})+g_{pr,st}(\hat{e}_{rq,st}+\hat{e}_{qr,st})\\},$
(A.5)

in terms of which the orbital gradients (15), the $\hat{H}_{\kappa}$
Hamiltonian (27), and the sigma vectors (31) and (32) can readily be obtained.
Further in view of the following identities

$\displaystyle D_{ip}$ $\displaystyle=$ $\displaystyle 2\delta_{ip},$ (A.6)
$\displaystyle\Gamma_{pq,ij}$ $\displaystyle=$ $\displaystyle
2\delta_{ij}D_{pq}-\delta_{iq}D_{pj},$ (A.7) $\displaystyle\Gamma_{pi,jq}$
$\displaystyle=$ $\displaystyle 4\delta_{ip}\delta_{jq}-\delta_{ij}D_{pq},$
(A.8)

the explicit expressions for the orbital gradients (15) and the diagonal
elements of the orbital Hessian (24) (which are required by the BFGS
algorithm) can readily be derived

$\displaystyle G^{\mathrm{o}}_{ia}$ $\displaystyle=$ $\displaystyle
4f_{ia}^{\mathrm{c}}+4f_{ia}^{\mathrm{a}},$ (A.9) $\displaystyle
G^{\mathrm{o}}_{it}$ $\displaystyle=$ $\displaystyle
4f_{it}^{\mathrm{c}}+4f_{it}^{\mathrm{a}}-2F_{ti},$ (A.10) $\displaystyle
G^{\mathrm{o}}_{ta}$ $\displaystyle=$ $\displaystyle 2F_{ta},$ (A.11)
$\displaystyle G^{\mathrm{o}}_{ut}$ $\displaystyle=$ $\displaystyle
2(F_{ut}-F_{tu}),$ (A.12) $\displaystyle E^{\mathrm{oo}}_{ia,ia}$
$\displaystyle=$ $\displaystyle
4\\{f_{aa}^{\mathrm{c}}+f_{aa}^{\mathrm{a}}-f_{ii}^{\mathrm{c}}-f_{ii}^{\mathrm{a}}+3K^{ii}_{aa}-J^{ii}_{aa}\\},$
(A.13) $\displaystyle E^{\mathrm{oo}}_{it,it}$ $\displaystyle=$ $\displaystyle
4\\{f_{tt}^{\mathrm{c}}+f_{tt}^{\mathrm{a}}-f_{ii}^{\mathrm{c}}-f_{ii}^{\mathrm{a}}+\frac{1}{2}f_{ii}^{\mathrm{c}}D_{tt}-\frac{1}{2}F_{tt}$
(A.14) $\displaystyle+$
$\displaystyle\sum_{u}(\delta_{tu}-D_{tu})(3K^{ii}_{tu}-J^{ii}_{tu})+\frac{1}{2}\sum_{u,v}J^{uv}_{ii}P^{vu}_{tt}+\sum_{u,v}K^{uv}_{ii}Q^{vu}_{tt}\\},$
$\displaystyle E^{\mathrm{oo}}_{ta,ta}$ $\displaystyle=$ $\displaystyle
4\\{\frac{1}{2}D_{tt}f_{aa}^{\mathrm{c}}+\frac{1}{2}\sum_{u,v}P^{uv}_{tt}J^{uv}_{aa}+\sum_{u,v}Q^{uv}_{tt}K^{vu}_{aa}-\frac{1}{2}F_{tt}\\},$
(A.15) $\displaystyle E^{\mathrm{oo}}_{tu,tu}$ $\displaystyle=$ $\displaystyle
2(1+P_{tu})W_{tu,tu}-4W_{tu,ut}-2(F_{tt}+F_{uu}),\quad\forall t>u,$ (A.16)

where

$\displaystyle f_{pq}^{\mathrm{c}}$ $\displaystyle=$ $\displaystyle
h_{pq}+\sum_{i}(2J^{ii}_{pq}-K^{ii}_{pq}),$ (A.17) $\displaystyle
f_{pq}^{\mathrm{a}}$ $\displaystyle=$
$\displaystyle\sum_{u,t}(J^{tu}_{pq}-\frac{1}{2}K^{tu}_{pq})D_{ut},$ (A.18)
$\displaystyle F_{up}$ $\displaystyle=$
$\displaystyle\sum_{v}D_{uv}f^{\mathrm{c}}_{vp}+\sum_{vwt}P^{wt}_{uv}J^{tw}_{vp},$
(A.19) $\displaystyle W_{tu,tu}$ $\displaystyle=$ $\displaystyle
D_{tt}f_{uu}^{\mathrm{c}}+\sum_{rs}P^{rs}_{tt}J^{sr}_{uu}+2\sum_{rs}Q^{rs}_{tt}K^{sr}_{uu},$
(A.20) $\displaystyle W_{tu,ut}$ $\displaystyle=$ $\displaystyle
D_{tu}f_{ut}^{\mathrm{c}}+\sum_{rs}P^{rs}_{tu}J^{sr}_{ut}+2\sum_{rs}Q^{rs}_{tu}K^{st}_{ut},$
(A.21) $\displaystyle Q^{vw}_{tu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\Gamma_{tv,uw}+\Gamma_{vt,uw})=Q^{wv}_{ut},$ (A.22)
$\displaystyle K^{rs}_{pq}$ $\displaystyle=$ $\displaystyle
g_{pr,sq}=K^{sr}_{qp}.$ (A.23)

## References

  * Cheung et al. 1978 Cheung, L.; Sundberg, K.; Ruedenberg, K. Dimerization of carbene to ethylene. _J. Am. Chem. Soc._ 1978, _100_ , 8024–8025
  * Roos et al. 1980 Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. M. A complete active space scf method (casscf) using a density-matrix formulated super-ci approach. _Chem. Phys._ 1980, _48_ , 157–173
  * Werner 1987 Werner, H.-J. Matrix-formulated direct multiconfiguration self-consistent field and multiconfiguration reference configuration-interaction methods. _Adv. Chem. Phys._ 1987, _69_ , 1–62
  * Roos 1987 Roos, B. O. The complete active space self-consistent field method and its applications in electronic structure calculations. _Adv. Chem. Phys._ 1987, _69_ , 399–445
  * Shepard 1987 Shepard, R. The multiconfiguration self-consistent field method. _Adv. Chem. Phys_ 1987, _69_ , 63–200
  * Schmidt and Gordon 1998 Schmidt, M. W.; Gordon, M. S. The construction and interpretation of MCSCF wavefunctions. _Annu. Rev. Phys. Chem._ 1998, _49_ , 233–266
  * Szalay et al. 2012 Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration self-consistent field and multireference configuration interaction methods and applications. _Chem. Rev._ 2012, _112_ , 108–181
  * Sun 2016 Sun, Q. Co-iterative augmented Hessian method for orbital optimization. _arXiv preprint arXiv:1610.08423_ 2016,
  * Kreplin et al. 2019 Kreplin, D. A.; Knowles, P. J.; Werner, H.-J. Second-order MCSCF optimization revisited. I. Improved algorithms for fast and robust second-order CASSCF convergence. _J. Chem. Phys._ 2019, _150_ , 194106
  * Kreplin et al. 2020 Kreplin, D. A.; Knowles, P. J.; Werner, H.-J. MCSCF optimization revisited. II. Combined first- and second-order orbital optimization for large molecules. _J. Chem. Phys._ 2020, _152_ , 074102
  * Vogiatzis et al. 2017 Vogiatzis, K. D.; Ma, D.; Olsen, J.; Gagliardi, L.; de Jong, W. A. Pushing configuration-interaction to the limit: Towards massively parallel MCSCF calculations. _J. Chem. Phys._ 2017, _147_ , 184111
  * Zhang et al. 2021 Zhang, N.; Liu, W.; Hoffmann, M. R. Further Development of iCIPT2 for Strongly Correlated Electrons. _J. Chem. Theory Comput._ 2021, _17_ , 949–964
  * Gidofalvi and Mazziotti 2008 Gidofalvi, G.; Mazziotti, D. A. Active-space two-electron reduced-density-matrix method: Complete active-space calculations without diagonalization of the N-electron Hamiltonian. _J. Chem. Phys._ 2008, _129_ , 134108
  * Fosso-Tande et al. 2016 Fosso-Tande, J.; Nguyen, T.-S.; Gidofalvi, G.; DePrince, A. E. Large-Scale Variational Two-Electron Reduced-Density-Matrix-Driven Complete Active Space Self-Consistent Field Methods. _J. Chem. Theory Comput._ 2016, _12_ , 2260–2271
  * Zgid and Nooijen 2008 Zgid, D.; Nooijen, M. The density matrix renormalization group self-consistent field method: Orbital optimization with the density matrix renormalization group method in the active space. _J. Chem. Phys._ 2008, _128_ , 144116
  * Ghosh et al. 2008 Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K.-L. Orbital optimization in the density matrix renormalization group, with applications to polyenes and beta-carotene. _J. Chem. Phys._ 2008, _128_ , 144117
  * Yanai et al. 2009 Yanai, T.; Kurashige, Y.; Ghosh, D.; Chan, G. K.-L. Accelerating convergence in iterative solution for large-scale complete active space self-consistent-field calculations. _Int. J. Quantum Chem._ 2009, _109_ , 2178–2190
  * Ma and Ma 2013 Ma, Y.; Ma, H. Assessment of various natural orbitals as the basis of large active space density-matrix renormalization group calculations. _J. Chem. Phys._ 2013, _138_ , 224105
  * Wouters et al. 2014 Wouters, S.; Bogaerts, T.; Van Der Voort, P.; Van Speybroeck, V.; Van Neck, D. Communication: DMRG-SCF study of the singlet, triplet, and quintet states of oxo-Mn(Salen). _J. Chem. Phys._ 2014, _140_ , 241103
  * Ma et al. 2017 Ma, Y.; Knecht, S.; Keller, S.; Reiher, M. Second-order self-consistent-field density-matrix renormalization group. _J. Chem. Theory Comput._ 2017, _13_ , 2533–2549
  * Sun et al. 2017 Sun, Q.; Yang, J.; Chan, G. K.-L. A general second order complete active space self-consistent-field solver for large-scale systems. _Chem. Phys. Lett._ 2017, _683_ , 291–299
  * Thomas et al. 2015 Thomas, R. E.; Sun, Q.; Alavi, A.; Booth, G. H. Stochastic Multiconfigurational Self-Consistent Field Theory. _J. Chem. Theory Comput._ 2015, _11_ , 5316–5325
  * Li Manni et al. 2016 Li Manni, G.; Smart, S. D.; Alavi, A. Combining the Complete Active Space Self-Consistent Field Method and the Full Configuration Interaction Quantum Monte Carlo within a Super-CI Framework, with Application to Challenging Metal-Porphyrins. _J. Chem. Theory Comput._ 2016, _12_ , 1245–1258
  * Smith et al. 2017 Smith, J. E. T.; Mussard, B.; Holmes, A. A.; Sharma, S. Cheap and near exact CASSCF with large active spaces. _J. Chem. Theory Comput._ 2017, _13_ , 5468–5478
  * Yao and Umrigar 2021 Yao, Y.; Umrigar, C. Orbital Optimization in Selected Configuration Interaction Methods. _arXiv preprint arXiv:2104.02587_ 2021,
  * Zimmerman and Rask 2019 Zimmerman, P. M.; Rask, A. E. Evaluation of full valence correlation energies and gradients. _J. Chem. Phys._ 2019, _150_ , 244117
  * Levine et al. 2020 Levine, D. S.; Hait, D.; Tubman, N. M.; Lehtola, S.; Whaley, K. B.; Head-Gordon, M. CASSCF with Extremely Large Active Spaces Using the Adaptive Sampling Configuration Interaction Method. _J. Chem. Theory Comput._ 2020, _16_ , 2340–2354
  * Park 2021 Park, J. W. Second-Order Orbital Optimization with Large Active Spaces Using Adaptive Sampling Configuration Interaction (ASCI) and Its Application to Molecular Geometry Optimization. _J. Chem. Theory Comput._ 2021, _17_ , 1522–1534
  * Liu and Hoffmann 2016 Liu, W.; Hoffmann, M. R. iCI: Iterative CI toward full CI. _J. Chem. Theory Comput._ 2016, _12_ , 1169–1178, (E) 2016, _12_ , 3000
  * Liu and Hoffmann 2014 Liu, W.; Hoffmann, M. R. SDS: the ‘static-dynamic-static’ framework for strongly correlated electrons. _Theor. Chem. Acc._ 2014, _133_ , 1481
  * Lei et al. 2017 Lei, Y.; Liu, W.; Hoffmann, M. R. Further development of SDSPT2 for strongly correlated electrons. _Mol. Phys._ 2017, _115_ , 2696–2707
  * Zhang et al. 2020 Zhang, N.; Liu, W.; Hoffmann, M. R. Iterative Configuration Interaction with Selection. _J. Chem. Theory Comput._ 2020, _16_ , 2296–2316
  * Huang et al. 2017 Huang, C.; Liu, W.; Xiao, Y.; Hoffmann, M. R. iVI: An iterative vector interaction method for large eigenvalue problems. _J. Comput. Chem._ 2017, _38_ , 2481–2499, (E) 2018, _39_ , 338
  * Huang and Liu 2019 Huang, C.; Liu, W. iVI-TD-DFT: An iterative vector interaction method for exterior/interior roots of TD-DFT. _J. Comput. Chem._ 2019, _40_ , 1023–1037, (E) 2018, _39_ , 338
  * 35 Lei, Y.; Suo, B.; Liu, W. $\mathbbm{i}$CAS: Imposed Automatic Selection and Localization of Complete Active Spaces. (unpublished).
  * Wang and Liu 2021 Wang, Z.; Liu, W. iOI: an Iterative Orbital Interaction Approach for Solving the Self-Consistent Field Problem. _arXiv e-prints_ 2021, arXiv–2105
  * Helgaker et al. 2014 Helgaker, T.; Jorgensen, P.; Olsen, J. _Molecular electronic-structure theory_ ; John Wiley & Sons, 2014
  * Werner and Meyer 1980 Werner, H.-J.; Meyer, W. A quadratically convergent multiconfiguration–self-consistent field method with simultaneous optimization of orbitals and CI coefficients. _J. Chem. Phys._ 1980, _73_ , 2342–2356
  * Werner and Knowles 1985 Werner, H.; Knowles, P. J. A second order multiconfiguration SCF procedure with optimum convergence. _J. Chem. Phys._ 1985, _82_ , 5053–5063
  * Fletcher 2013 Fletcher, R. _Practical methods of optimization_ ; John Wiley & Sons, 2013
  * Fischer and Almlof 1992 Fischer, T. H.; Almlof, J. General methods for geometry and wave function optimization. _J. Phys. Chem._ 1992, _96_ , 9768–9774
  * Chaban et al. 1997 Chaban, G.; Schmidt, M.; Gordon, M. Approximate second order method for orbital optimization of SCF and MCSCF wavefunctions. _Theor. Chem. Acc._ 1997, _97_ , 88–95
  * Ivanic and Ruedenberg 2003 Ivanic, J.; Ruedenberg, K. A MCSCF method for ground and excited states based on full optimizations of successive Jacobi rotations. _J. Comput. Chem._ 2003, _24_ , 1250–1262
  * Eriksen et al. 2020 Eriksen, J. J.; Anderson, T. A.; Deustua, J. E.; Ghanem, K.; Hait, D.; Hoffmann, M. R.; Lee, S.; Levine, D. S.; Magoulas, I.; Shen, J.; Tubman, N. M.; Whaley, K. B.; Xu, E.; Yao, Y.; Zhang, N.; Alavi, A.; Chan, G. K.-L.; Head-Gordon, M.; Liu, W.; Piecuch, P.; Sharma, S.; Ten-no, S. L.; Umrigar, C. J.; Gauss, J. The ground state electronic energy of benzene. _J. Phys. Chem. Lett._ 2020, _11_ , 8922–8929
  * Liu et al. 1997 Liu, W.; Hong, G.; Dai, D.; Li, L.; Dolg, M. The Beijing 4-component density functional theory program package (BDF) and its application to EuO, EuS, YbO and YbS. _Theor. Chem. Acc._ 1997, _96_ , 75–83
  * Liu et al. 2003 Liu, W.; Wang, F.; Li, L. _J. Theor. Comput. Chem._ 2003, _2_ , 257–272
  * Liu et al. 2004 Liu, W.; Wang, F.; Li, L. In _Recent Advances in Relativistic Molecular Theory_ ; Hirao, K., Ishikawa, Y., Eds.; World Scientific: Singapore, 2004; pp 257–282
  * Liu et al. 2004 Liu, W.; Wang, F.; Li, L. In _Encyclopedia of Computational Chemistry_ ; von Ragué Schleyer, P., Allinger, N. L., Clark, T., Gasteiger, J., Kollman, P. A., Schaefer III, H. F., Eds.; Wiley: Chichester, UK, 2004
  * Zhang et al. 2020 Zhang, Y.; Suo, B.; Wang, Z.; Zhang, N.; Li, Z.; Lei, Y.; Zou, W.; Gao, J.; Peng, D.; Pu, Z.; Xiao, Y.; Sun, Q.; Wang, F.; Ma, Y.; Wang, X.; Guo, Y.; Liu, W. BDF: A relativistic electronic structure program package. _J. Chem. Phys._ 2020, _152_ , 064113
  * Dunning 1989 Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. _J. Chem. Phys._ 1989, _90_ , 1007–1023
  * Hachmann et al. 2007 Hachmann, J.; Dorando, J. J.; Avilés, M.; Chan, G. K.-L. The radical character of the acenes: a density matrix renormalization group study. _J. Chem. Phys._ 2007, _127_ , 134309
  * Schriber et al. 2018 Schriber, J. B.; Hannon, K. P.; Li, C.; Evangelista, F. A. A Combined Selected Configuration Interaction and Many-Body Treatment of Static and Dynamical Correlation in Oligoacenes. _J. Chem. Theory Comput._ 2018, _14_ , 6295–6305
  * Olivares-Amaya et al. 2015 Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. The ab-initio density matrix renormalization group in practice. _J. Chem. Phys._ 2015, _142_ , 034102
  * Li Manni et al. 2016 Li Manni, G.; Smart, S. D.; Alavi, A. Combining the Complete Active Space Self-Consistent Field Method and the Full Configuration Interaction Quantum Monte Carlo within a Super-CI Framework, with Application to Challenging Metal-Porphyrins. _J. Chem. Theory Comput._ 2016, _12_ , 1245–1258
  * Li Manni and Alavi 2018 Li Manni, G.; Alavi, A. Understanding the Mechanism Stabilizing Intermediate Spin States in Fe(II)-Porphyrin. _J. Phys. Chem. A_ 2018, _122_ , 4935–4947
  * Groenhof et al. 2005 Groenhof, A. R.; Swart, M.; Ehlers, A. W.; Lammertsma, K. Electronic Ground States of Iron Porphyrin and of the First Species in the Catalytic Reaction Cycle of Cytochrome P450s. _J. Phys. Chem. A_ 2005, _109_ , 3411–3417
  * Pierloot 2003 Pierloot, K. The CASPT2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes. _Mol. Phys._ 2003, _101_ , 2083–2094
  * Weser et al. 2021 Weser, O.; Freitag, L.; Guther, K.; Alavi, A.; Li Manni, G. Chemical insights into the electronic structure of Fe(II) porphyrin using FCIQMC, DMRG, and generalized active spaces. _Int. J. Quantum Chem._ 2021, _121_ , e26454

For TOC only