Automated construction of molecular active spaces from atomic valence orbitals
Elvira R. Sayfutyarova Qiming Sun Garnet K. -L. Chan Gerald Knizia




We introduce the atomic valence active space (AVAS), a simple and well-defined automated technique for constructing active orbital spaces for use in multi-configuration and 
multi-reference (MR) electronic structure calculations.
Concretely, the technique constructs active molecular orbitals capable of describing all relevant electronic configurations emerging from a targeted set of atomic valence orbitals (e.g., the metal d orbitals in a coordination complex).
This is achieved via a linear transformation of the occupied and unoccupied orbital spaces from an easily obtainable single-reference wavefunction (such as from a Hartree-Fock or Kohn-Sham calculations) based on projectors to targeted atomic valence orbitals.
We discuss the premises, theory, and implementation of the idea, and several of its variations are tested.
To investigate the performance and accuracy, we calculate the excitation energies for various transition metal complexes in typical application scenarios.
Additionally, we follow the homolytic bond breaking process of a Fenton reaction along its reaction coordinate.
While the described AVAS technique is not an universal solution to the active space problem, its premises are fulfilled in many application scenarios of transition metal chemistry and bond dissociation processes.
In these cases the technique makes MR calculations easier to execute, easier to reproduce by any user, and simplifies
the determination of the appropriate size of the active space required for accurate results.




Manifestations of minimum-bias dijets in high-energy nuclear collisions

    Thomas A. Trainor
    December 30, 2023
========================================================================





§ INTRODUCTION

Multiconfigurational and multireference (MR) methods remain indispensable in the treatment of 
challenging electronic structure problems. 
Transition metal complexes provide a rich source of examples, as they often feature strongly correlated electronic degrees of freedom, which render density functional theory (DFT) calculations unreliable.
Unfortunately, MR methods require an a priori choice of a suitable set of active molecular orbitals (an active space), which critically determines the quality of the results.
The choice of active space is non-trivial and often represents a major challenge in practical computations.

The standard way to choose an active space is as follows: First, a full set of molecular orbitals of the molecule is computed with a simple electronic structure method, such as Hartree-Fock (HF) or Kohn-Sham DFT (KS-DFT). Second, these molecular orbitals, both occupied and unoccupied, are  visually inspected, and based on their shape, energy, occupation numbers, etc., one selects into the active space the set of orbitals which are expected to be chemically most relevant
(e.g., with significant transition metal d-electron character or ligand character in the case of transition metal complexes—regions which are empirically known to be important).

But despite the existence of general advice on how to select a good set of starting orbitals <cit.>, there are some problems with this approach.
First, the selection of molecular orbitals for the active space is performed by the user, normally based on personal experience.

This makes MR methods hard to apply, and gives results that are complex to reproduce and hard to judge in terms of quality.
This stands in contrast to single-reference calculations which do not require active spaces and therefore do not have this level of arbitrariness. 
Second, the molecular orbitals mix together valence orbitals of different character; 
for example, typically a metal's d atomic valence orbitals contribute to a very large number of molecular orbitals, and it is often not easy to truncate these to a small active subset in such a way that they retain the capability of describing all the right physics in complicated systems.
In the case of large complexes, especially with multiple metal ions, this procedure commonly becomes a matter of trial and error.

Techniques to aid in the construction of high quality active spaces are therefore highly desirable, particularly as more  powerful electronic structure methods are becoming available which are capable of treating increasingly large active spaces efficiently.<cit.>
There have been a number of contributions in this area, and the general procedures can loosely be summarized as based on
 estimating the (correlated) occupation numbers of the orbitals (including the full single orbital density matrix in Ref. Keller:ActiveSpaceSelection,Stein:ActiveSpaceSelection), followed by selecting the active orbitals based on partial occupancy according to an input threshold. The various techniques differ primarily 
in how the occupation number information is obtained. For example, correlated occupation numbers
can be estimated from unrestricted Hartree-Fock <cit.> or Kohn-Sham calculations or from correlated calculations,
such as MP2 <cit.> or approximate DMRG calculations <cit.>.

Nonetheless, while these existing automated approaches are advancements from typical ad hoc active space constructions, there is still room for further improvement.
For example, an obvious drawback of the above procedures is that they all require a non-trivial preliminary calculation:
either a correlated calculation must be performed, or a suitable broken symmetry solution must be found.
This is not always possible: for example, there may not always be a broken symmetry solution involving the region of interest, or
the preliminary correlated calculation may simply be too costly. 
A perhaps less obvious drawback, but one which one encounters in practice, is that there is no guarantee that the active orbitals
found in these automated procedures are actually spatially located in the region of chemical interest. For example, in 
models of enzymatic binding sites, particularly those with charged ligands, the unpaired electrons may lie in functional groups which are spatially 
far from and irrelevant to the chemistry of the metal center. In this case, an additional inspection is once again required to choose the subset of
active orbitals of chemical interest.

Here we propose an alternative approach to construct molecular active spaces for multireference problems automatically and systematically.
This approach does not suffer from the drawbacks mentioned above, and is particularly well-suited for practical calculations involving 
transition metal complexes. In its simplest formulation, the procedure  requires only an easily obtained single-determinant reference function, together with a choice of target atomic valence orbitals.
The technique is based on the following idea: It is empirically known that one can typically identify a small set of atomic valence orbitals which give rise to the strong correlation effects (for example, d atomic valence orbitals in transition metal complexes). 
We therefore aim to construct a set of active molecular orbitals by defining them in terms of these atomic valence orbitals.
Concretely, using simple linear algebra, we can define mathematical rotations of the occupied and virtual molecular orbitals of the single-determinant reference function which maximizes their given atomic valence character (e.g. 3d character). 
After this rotation, the relevant molecular orbitals to include into the active space can be selected automatically. 
The ideas employed here are closely related to earlier work of Iwata<cit.> and Schmidt and coworkers<cit.>  regarding the construction of approximate valence spaces;
however these earlier works focused on constructing molecular orbitals capable of spanning the entire valence space (not only the valence space of selected atomic orbitals, with the goal of identifying active spaces for processes involving said atomic orbitals), and the latter work also markedly differs in intent.
Similar mathematical techniques were also used in the construction of molecule-intrinsic minimal basis sets capable of spanning a set of given occupied orbital spaces<cit.>.


Sec. <ref> explains the motivation, background, theory, and implementation of the construction proposed here.
Sec. <ref> then describes criteria for judging the quality of the constructed active spaces, and discusses its application to a large number of prototypical MR calculations of metal complexes.
Sec. <ref> describes conclusions and possible implication for future work.





§ THEORY





 §.§ What is an active space?


Many naturally occuring stable molecules have an electronic structure which is qualitatively well described by a single-determinant self-consistent field (SCF) wave function, 
such as used in Kohn-Sham Density Functional Theory (KS-DFT) or Hartree-Fock (HF).
In this case, the entire space of molecular orbitals (MO) is strictly divided into fully occupied and 
fully unoccupied molecular (spin-) orbitals.
However, there are important classes of chemical systems in which this picture breaks down and where a superposition of multiple N-electron determinants is required to describe the electronic structure even qualitatively.
This phenomenon is sometimes called strong correlation; prototypical cases in which it occurs are (a)  the process of homolytic bond breaking, and (b) various kinds of transition metal complexes—particularly when the complex is in an overall low-spin state generated
by coupling of the metal to other metals or redox non-innocent ligands, as frequently encountered in catalysis and in bio-inorganic systems.

Both of these prototypical cases share the same root cause: The occurrence of valence atomic orbitals with energy levels similar to other valence orbitals, but with poor orbital overlap, 
giving rise to small energy splittings between bonding and anti-bonding linear combinations.
In these cases quantum resonances between both bonding and anti-bonding orbitals must be explicitly considered to qualitatively describe the electronic structure—and single-determinant wave functions are incapable of doing so, because in these each MO either is occupied or is not, but not both.
The core idea developed in this manuscript is that in these two most important cases, the emergence of strong correlation is tightly linked to specific valence atomic orbitals, which are easy to identify. In the case of transition metals, these are the compact d orbitals (and possibly some specific ligand orbitals), and in the case of bond dissociation, these are the valence atomic orbitals of the dissociating atom(s).
Thus, in the following we will assume that a small number of specific valence atomic orbitals are explicitly selected by the user (e.g., the d orbitals of metal centers), and that the goal is to construct an active space suitable for describing all highly relevant determinants that they give rise to.
We stress that this initial selection is easy in practice—the core problem is how to use the information to build a suitable active space.

We thus briefly consider what an active space is. 
An active space wavefunction is one where the superpositions
of determinants are restricted so that varied occupations are found only within the active orbitals {ϕ_1 …ϕ_m}. Technically, this means that the wavefunction can be written as 
a second quantized product

    |Ψ⟩ = {ϕ_1 …ϕ_m } | core⟩

where {ϕ_1 …ϕ_m} denotes a general occupancy wavefunction within the active orbitals, and |core⟩
denotes a single determinant. 
It is clear that the active orbitals must span the space of our chosen specific valence atomic orbitals. However, 
Eq. (<ref>) additionally implies that the rest of the molecule must be well described by the single core determinant.
To achieve this, the active space must contain orbitals additional to our set of specific valence atomic orbitals, and
to remain compact, we need to define the minimal additional set. 


To this end, we here employ the basic observation in density matrix embedding theory (DMET),<cit.> which 
describes how to construct such an active space explicitly.
In particular, DMET tells us that the active space with the above properties is at most twice the size of the initial set of chosen valence atomic orbitals. 
However, in this work we  modify the presentation and practice of the DMET procedure
to make it more natural in the active space setting.
In particular, in contrast to the original presentation in referencesGK_GKC1,GK_GKC2 (but as described in the appendix of Ref. zheng2016ground)
we do not introduce separate “fragment” and “bath” orbitals, but rather retain the occupied and virtual character of the constructed active orbitals.
This has the important benefit that it leads to a natural truncation procedure, which allows us to further reduce the size of the active space in
the most chemically meaningful way.
In Sec. <ref> we describe the isolation of entangled orbitals for the active space construction, and Sec. <ref> will provide technical details and discuss various practical aspects relevant in the active space case.
Finally, Sec. <ref> outlines how the presented approach may be extended to more complex cases, such as double-shell effects or complex metal/ligand interactions.




 §.§ Isolating target-overlapping orbitals for the active space


Let A={|p⟩} denote the (small) set of chosen target valence
atomic orbitals (not necessarily orthonormal), which we expect to be
responsible for strong correlations (e.g. the five 3d orbitals of a
third row transition metal atom of a complex, details on their
selection and representation will follow).  Let |Φ⟩ denote
a Slater determinant, which represents the electronic structure of our
system of interest at the SCF level (HF or KS-DFT).

Our method for active space construction is built around the following
physical assumptions, which are those used in 
DMET<cit.>:(a) a
SCF wave function |Φ⟩ may be unable to describe
precisely how our target atomic orbitals are bonded with the rest of
the molecule; however, it will generally describe the rest of
the molecule reasonably well (experience in DMET suggest that this works
even in cases where the SCF wave function as a whole is qualitatively
very wrong<cit.>), and (b) we can isolate the part of
|Φ⟩ which involves our target AOs from the part which does
not, and employ the former part as the active space (therefore allowing it
to be replaced by a more powerful wave function description), and retain the
simple determinantal description of rest.


To this end, we employ a rotation within the set of occupied molecular
orbitals of |Φ⟩ which splits them into two groups: one group
which has overlap with our target AOs, and one group which does not.
A simple dimensional counting argument will show that for a set of
|A| selected target AOs, there is a rotation of the occupied
molecular orbitals such that at most |A| of them have non-zero
weight on the target AOs.  We similarly split  the set of virtual
molecular orbitals into one group of at most |A| virtual orbitals
which have weight on the target atomic valence orbitals, and the other
group which does not.  The idea is now to explicitly construct these
rotated orbital groups, and then employ the at most 2|A| combined
occupied and virtual orbitals with target overlap as active orbitals,
while leaving the other occupied and virtual orbitals without target
overlap as inactive (closed-shell) or virtual orbitals, respectively,
in the following multi-configuration treatment.  As all of the
selected AOs in A then lie within the span of this active space, the
resulting multi-configurational wave function is then capable of
representing arbitrary quantum resonances involving the target AOs.

We first discuss the occupied case.  Let i = 1… N_occ
denote N_occ occupied molecular orbitals (MOs) of
|Φ⟩.  The projector P̂ onto the space of atomic
orbitals in A, is given by

    P̂   = ∑ _p,q∈ A |p ⟩ [σ^-1]_pq⟨ q |.

Here σ denotes the |A|×|A| target AO overlap matrix with
elements [σ]_pq = ⟨ p|q⟩, and σ^-1 its
matrix inverse.  
Employing these projectors, we construct the first set of active
orbitals by rotating |Φ⟩'s occupied MOs {|i⟩} as
follows.
First, we calculate the N_occ× N_occ overlap
matrix of occupied orbitals projected onto (A), the space of
selected target atomic orbitals:

    [ S^A]_ij = ⟨ i|P̂|j ⟩,

where i,j are occupied orbital indices.  Next we compute the
N_occ× N_occ (unitary) matrix of eigenvectors
[ U]_ij of S^A, such that

    S^A U =  U
       diag(σ_1,…,σ_N_occ)

(where diag(…) denotes a diagonal matrix of the
given elements), or, written in component form,

    ∀ i,j: ∑_k [ S^A]_ik [ U]_kj = [ U]_ijσ_j.

There are at most |A| non-zero eigenvalues {σ_i}, because
(A) is a |A|-dimensional space, and [ S]^A_ij
involves a projection onto it.  Furthermore, the eigenvectors [
  U]_ij of S^A_ij define a rotation on the occupied orbitals:

    |i⟩  ↦  |ĩ⟩ = ∑_k|k⟩ [
         U]_ki,

which clearly separates them into two groups: The at most |A|
rotated occupied orbitals |ĩ⟩ with σ_i≠ 0,
which have non-vanishing overlap with our target atomic orbitals (and
which therefore should go into the active space), and
the remaining |ĩ⟩ with σ_i= 0 which have no
overlap with our target atomic orbitals, and therefore can stay as
inactive (inner closed shell) orbitals in the subsequent
multiconfigurational methods.


Note that the rotated occupied orbitals {|ĩ⟩} in
(<ref>) are obtained as a unitary
transformation of |Φ⟩'s original occupied orbitals
{|i⟩}.  Consequently, a determinantal wave function
|Φ̃⟩ built from the {|ĩ⟩} is physically
equivalent (differs by at most a phase factor) from the original
determinant |Φ⟩.  That is, so far we have done nothing to
|Φ⟩ except for splitting its occupied orbitals into a convenient
set of at most |A| orbitals related to our |A| target AOs and the
remaining set we can treat as inactive.

We then proceed similarly for the virtual orbitals {|a⟩, 1 =
1…,N_vir} of |Φ⟩: We form the
N_vir× N_vir projected overlap matrix

    [S̅]^A_ab = ⟨ a|P̂|b ⟩,

where a,b are virtual orbital indices, then find its unitary matrix
of eigenvectors U̅ such that

    S̅^AU̅ = U̅diag(σ_1,…,σ_N_vir)
           ⇔    ∀ a,b: ∑_c [S̅^A]_ac [U̅]_cb = [U̅]_abσ_b,

and use U̅ to rotate the virtual orbitals via

    |a⟩  ↦  |ã⟩ = ∑_c|c⟩
           [U̅]_ca.

Again, the at most |A| of the new virtual orbitals {|ã⟩} with eigenvalues σ_a≠ 0 are selected for the active
space, while the remaining orbitals will stay unoccupied in the
subsequent multi-configuration treatment.

Finally, having active orbitals with overlap with A from both sides,
occupied and unoccupied orbitals in |Φ⟩, we can form the
total active space by combining the sets of {|ĩ⟩} and
{|ã⟩} with non-zero projected overlap eigenvalues.
Since the combined set includes all orbitals which have
non-vanishing overlap with our target space (A), all of the
selected AOs in A then lie within the span of this active space.
Therefore, a multi-configurational wave function with this active
space will be capable of representing arbitrary quantum resonances
involving the target AOs—which was the goal of our construction. 




 §.§ Technical details of the construction

Sec. <ref> discusses the formal framework of the active space construction.
However, several practical aspects still need to be discussed:
(a) How are the target AOs A={|p⟩} chosen and represented?
(b) How are the actual rotation matrices U (eq. (<ref>)) and U̅ (eq. (<ref>)) computed in practice?
(c) Can the active space (formally twice the number of the target AOs) be further reduced in size?
(d) How should open-shell systems be handled? (In particular, what to do for restricted open-shell functions |Φ⟩?)
We will discuss these questions in the current and next subsections.

Let us first assume that |Φ⟩ is a closed-shell Slater determinant obtained from an SCF calculation. 
Its occupied and virtual molecular orbitals are expressed as:

    |i ⟩   = ∑_μ∈ B_1  |μ⟩ C^μ_i
     |a ⟩   = ∑_μ∈ B_1  |μ⟩C̅^μ_a,

where μ are basis functions from the (large) computational basis set B_1 (e.g., cc-pVTZ or def2-TZVPP), and C^μ_i=[ C_occ]_μ i and C̅^μ_a=[ C_vir]_μ a are the coefficients of the basis function μ in the expansion of the occupied orbital i and virtual orbital a, respectively.
C_occ and C_vir denote the |B_1|× N_occ occupied and |B_1|× N_vir virtual sub-matrices of the |B_1|× |B_1| SCF orbital matrix C (note that N_occ+N_vir=|B_1|—each orbital is either occupied or virtual).


In general, computational basis sets such as B_1 do not contain basis functions directly corresponding to AOs of any sort.
For this reason, we here select our target AOs A={|p⟩} based on a second auxiliary basis set B_2. 
This is a minimal basis set of tabulated free-atom AOs (here MINAO is used<cit.>; but other choices, such as subsets of ANO-RCC,<cit.> ano-pVnZ,<cit.> or ANO-VT-XZ<cit.>, could be considered; see also Sec. <ref>).
This choice leads to  simple expressions for the projected overlap matrices

    S^A_ij = ⟨i|P̂|j|=⟩∑_μμ'  C^μ_i P_μμ' C^μ'_j 
    S̅^A_ab = ⟨a|P̂|b|=⟩∑_μμ'C̅^μ_a P_μμ'C̅^μ'_b,

where the matrix elements of the projector are

    P_μμ' =∑_pp'⟨μ|p ⟩ [σ^-1]_pp'⟨ p'| μ' ⟩.

Combining all formulas into a numerical algorithm, the rotated orbitals are constructed as follows:

 
  * Let A ⊂ B_2 denote the subset of AOs we choose as target AOs for the active space construction (for example, the five 3d AOs in a transition metal complex with one metal center). 
 
  * Form the overlap matrix σ with elements σ_pp' = ⟨ p | p' ⟩, where p, p' ∈ A, as well as
 σ's inverse matrix, with elements σ^pp'=[σ^-1]_pp'. Both matrices have dimension of |A| × |A|.
 
  * Form the overlap matrix 𝐒_21 between the functions p of A ⊂ B_2 and the functions μ of the large basis set B_1, with elements [ S_21]_pμ= ⟨ p | μ⟩.
 
  * Form the projector P_μμ' =∑_pp'⟨μ|p ⟩σ^pp'⟨ p'| μ' ⟩ , or 𝐏=𝐒_21^†σ^-1𝐒_21.
 
  * Form the projected overlap matrices 𝐒^A= C_occ^† P  C_occ for the occupied orbitals (eq. (<ref>)), and S̅^A= C_vir^† P  C_vir for the virtual orbitals  (eq. (<ref>)).

  * Finally, diagonalize both projected overlap matrices to obtain the transformation matrices separating the MO sets by overlap with (A). 
Concretely, diagonalize 𝐒^A to obtain the eigenvector matrix U, and use it to find the transformed occupied orbital matrix C̃_occ =  C_occ U. 
Then diagonalize 𝐒̅^A to obtain the eigenvector matrix U̅, and use it to find the transformed virtual orbital matrix C̃_vir =  C_virU̅. These are the expansion coefficients of {|ĩ⟩} (eq. (<ref>)) and {|ã⟩} (eq. (<ref>)), respectively.



Rather than using the minimal basis B_2 directly, one could consider choosing the target AOs from a set of polarized AOs which take the molecular environment into account, such as the Intrinsic Atomic Orbitals (IAOs).<cit.>
If the IAOs are given as

    |p̃⟩ = ∑_μ |μ⟩ T_μ p,

where T_μ p denotes the elements of the |B_1|× |B_2| IAO transformation matrix,<cit.> this can be incorporated by updating the projection matrix in S^A_ij and S^A_ab in eqs. (<ref>) and (<ref>) as

    𝐏 = 𝐒𝐓(𝐓^†𝐒𝐓)^-1𝐓^†𝐒.

For simplicity, we do not follow this approach in this work, and choose the target AOs directly from the minimal basis B_2 as described above. 

Additionally, in some cases one may wish to extend the target space by some AOs which lie outside the valence space. This scenario is discussed in Sec. <ref>.









 §.§ Truncating the active space


The eigenvalues σ_i and σ_a of the projected overlap matrices 
in eqs. (<ref>) and (<ref>) 
reflect the degree to which the transformed orbitals |ĩ⟩  and |ã⟩ overlap with the space of our target AOs.
If we include every such transformed orbital with σ_a ≠ 0 and
σ_i ≠ 0 into our active space, 

the resulting CAS space will exactly include all electronic configurations which can be formed over the given AOs and the maximum size of the CAS space is twice
that of the target AO space. However, this CAS space may be too large and we may need to truncate it.
Here we can use the fact that often many of the σ_i and σ_a are small.
As a practical measure, we can set a threshold, such as 0.05–0.1, to exclude MOs with negligible overlap with (A).
In addition to reducing the size of the active space, this can further 
improve the reproducibility of calculations in the case of very small eigenvalues. 
This threshold becomes the only numerical parameter to be chosen by the user, and together with the selection of target AOs and the type of the SCF wave function |Φ⟩, it fully determines the active space.

Of course, if truncation is used, the active space no longer captures all possible configurations involving the target AOs which can be formed. However, the truncation does not affect the quality
of the description of the rest of the molecule by the core determinant in Eq. (<ref>). This guarantees
that the CASCI energy lies below the variational HF energy. We can also imagine the opposite tradeoff,
where one obtains a truncated active space which retains the ability to capture all possible configurations involving the target AOs, at the cost of worsening
the quality of the core determinant which describes the rest of the molecule.
(In DMET language, this would correspond to truncating the “bath” orbitals,
which is considered in Refs. Wouters2016 and qmmmdmet).
However, this may be a worse truncation procedure in the current setting, as the energy gained by treating the fluctuations in occupation number in the target AOs (such as a TM
3d shell) may not make up for the energy lost in incompletely describing the mean-field hybridization between the target AOs and the rest of the molecule.  
In particular, this second bath truncation procedure can, in principle, lead to a CASCI energy above the variational HF energy.




 §.§ Treatment of open-shell systems



If |Φ⟩ is a closed-shell determinant, then the active space construction algorithm can be directly used as described in Sec. <ref>.
However, in the case of open shell determinants, several choices can be considered:

  
  * One may perform the algorithm separately for alpha and beta orbitals, thus creating active orbitals with different spatial parts for alpha and beta-spin electrons.
While this choice is the most straight-forward and, arguably, creates the best initial active orbitals in the open-shell case, this option is not directly feasible if a spin-adapted multiconfigurational calculation will follow—most existing MCSCF programs cannot use such unrestricted orbitals (although this is implemented in the code we use here <cit.>).
A possible remedy for this problem would be to construct a single set of “corresponding orbitals”<cit.> from the separate alpha- and beta-orbital sets, but this has not been tested here.

  *  One may use exclusively the alpha
orbitals to construct the active space orbitals (and inactive orbitals determining the core determinant).
This treatment can be applied to both restricted and unrestricted SCF functions |Φ⟩ in a simple manner.
The rationale for this is that in the restricted open-shell case, the occupied beta orbitals lie entirely within the linear span of the occupied alpha orbitals, so one can argue that this choice takes care of both spin cases.
However, this argument is somewhat misleading because it may lead to some unoccupied beta orbitals being transformed into the core space, therefore enforcing their occupation with two electrons.
This error in the core means that the CASCI energy may be higher than the variational HF energy.
If the singly-occupied MOs in |Φ⟩ have only small components on the target valence AOs, this can lead to very bad CASCI wavefunctions.


  *  If a ROHF determinant |Φ⟩ is used, one may apply the construction of Sec. <ref> exclusively to the doubly-occupied and fully unoccupied orbitals of |Φ⟩, to form the core determinant and initial part of the active space, and then include additionally all the singly-occupied orbitals of |Φ⟩ into the active space.
The CASCI energy is then guaranteed to be below the variational HF energy, and further spin-adaptation can be used.
The main drawback is that the active space is usually larger in this procedure.


By default, we use method <ref>, however, we compare the different schemes in one of the systems below.








§ COMPUTATIONAL DETAILS



We implemented the atomic valence active space (AVAS) construction within the PySCF<cit.> package. 
All Complete Active Space Self-Consistent Field (CASSCF), Complete Active Space Configuration Interaction (CASCI), and strongly contracted N-electron valence state perturbation theory (NEVPT2)<cit.>
calculations were carried out using PySCF. 
For active spaces with more than 16–17 orbitals we used the Block code<cit.> through the  PySCF-Block interface to perform DMRG calculations in the active space.
The AVAS construction is also being implemented in a development version of Molpro.

For simplicity, we used all-electron cc-pVTZ-DK<cit.> basis sets for all systems, apart from the Fenton reaction (vide infra).
For the auxiliary minimal basis B_2 used to choose the target AOs, we employed the MINAO basis<cit.>, which is a truncated subset of the cc-pVTZ basis; for most atoms, this set 
consists of spherically averaged ground-state Hartree-Fock orbitals for the free atoms.
We did not use point group symmetry in the present calculations, since in a straight-forward implementation, the MOs do not necessarily retain symmetry-adaption after rotation; however, if symmetry-adaptable sets of target AOs are chosen, symmetry respecting orbital rotations can in principle be constructed.

For completeness, scalar relativistic effects were included using the exact-two-component (X2C) approach,<cit.> but this did not lead to significant differences from
the non-relativistic calculations in any of the considered examples. Spin-orbit coupling was not considered.


For the case of the Fenton reaction, the geometries were optimized with the DSCF and GRAD modules of Turbomole 7.0. This was done at the level of symmetry-broken unrestricted B3LYP<cit.> with def2-TZVP basis sets<cit.>, starting from the structures provided in Ref. petit:FentonReaction. The RI approximation was not employed, and solvation effects were not considered.
The characters of the starting geometry, transition state geometry, and product geometry were confirmed by computing the analytic nuclear Hessians at these points, via the AOFORCE module.<cit.>
The reaction path was computed by tracing the intrinsic reaction coordinate<cit.> in both directions,

starting at the transition state geometry.
This was done using Turbomole's DRC module.
Finally, the structures of the starting point, transition state, product, and IRC segments (of both directions) were joined, aligned, and compressed using the development version of IboView.<cit.>

Employing these geometries, the reported multi-configuration AVAS calculations were performed with PySCF, using  cc-pVTZ orbital basis sets and a non-relativistic Hamiltonian.


The used molecular geometries are supplied in the supporting information, as are selected visualizations of the obtained active spaces.
Orbital visualizations were made with IboView,<cit.> and show iso-surfaces enclosing 80% of the orbital's electron density.




§ RESULTS AND DISCUSSION


Judging the quality of an initial active space is difficult.
Here we employ two complementary criteria:




  * In the course of a CASSCF calculation, the active space orbitals are optimized.
If the overlap of the optimized active space with our initial active space guess remains high, we here regard our active space guess as being of high quality. To quantify this aspect, we compute the N_act× N_act overlap matrix

    S_change = ( C_act^final)^† S ( C_act^initial)

between the initial guess and the optimized final active orbitals, and compute its singular value decomposition.
If all singular values are close to 1.0, the active space remains mostly unchanged during the optimization. 
On the contrary, each singular value close to 0.0 indicates that an initial active orbital
had to be completely replaced by an unrelated orbital. The latter case may indicate that the initial active space lacks the capability to represent some essential features of the strongly correlated electronic structure of the given molecule (since the optimized active orbital, which presumably is needed to represent energetically important electron configurations of the system, lies completely outside the linear span of the initial active space guess); the initial active space guess may therefore need to be changed or enlarged.


  * We compare our computed excitation energies to experimental results and other high-level calculations reported in the literature.
When constructing a well-behaved series of active spaces, we should be able to see convergence or stability of the computed properties
with respect to the active space size.

We note that, strictly speaking, neither criterion can establish an initial active space as “definitely good”: While criterion 1 tests whether the CASSCF-optimized active space is close to the initial active space, there is no formal guarantee that the optimized CASSCF wave function itself is the best representative of the sought after electronic state (e.g., the optimization algorithm could be stuck in an undesirable local minimum). While criterion 2 establishes the compatibility with some specific physical properties and reference systems, it cannot guarantee that no other physical properties or systems exist which need different active spaces for a proper description.

However, both criteria could establish our active space construction to be unfit for its target applications—by showing counter-examples and negative results.






In order to investigate the basic properties of the AVAS procedure, 
we now apply these criteria in calculations on various transition-metal complexes.




 §.§ A. Ferrocene


We begin by considering the electronic structure of ferrocene, Fe(C5H5)2.
The ground state of ferrocene is dominated by a single configuration with Fe d^6. 
Both the MO analysis in Ref. fecp2_1 and our CI expansion coefficients indicate that the lowest excited states of ferrocene
have significant multiconfigurational character. 

We carried out an initial restricted  Hartree-Fock (ROHF) calculation for the singlet ground state.
We used the optimized D_5h geometry  from Ref. fecp2_geom (using the cc-pw-CVTZ basis set at the CCSD(T) level). 
In this geometry, the two cyclopentadienyl (Cp) rings are planar and the z-axis is aligned with the Cp-Fe-Cp axis.




As our target set of AOs, we first chose the five 3d orbitals of Fe.
Using a threshold of 0.1, the AVAS scheme produces a seven orbital active space: five orbitals from the occupied orbital space and two orbitals from the unoccupied orbital space are combined into a (10e,7o) active space; the orbitals are visualized in Fig. <ref>.
The five overlap eigenvalues above the threshold from the occupied space are 0.325, 0.325, 0.973, 0.973, 0.995, and the two eigenvalues from the unoccupied space are 0.675, 0.675.
There are three omitted unoccupied orbitals: two of them have only 2.7% weight in the 3d AO space and (3d_xy, 3d_x^2-y^2) character, and one has 0.6% weight in the 3d AO space, corresponding to the 3d_z^2 AO. 
The two orbitals from the occupied orbital space with eigenvalues 0.325 and 0.325, as well as the two selected MOs from the unoccupied orbital space with eigenvalues 0.675, 0.675,
both have (3d_xz, 3d_yz) character; based on this, we conclude that (3d_xz, 3d_yz) are most strongly involved in bonding with the two Cp rings, in agreement with the bonding picture established in earlier  studies <cit.>.

To construct a second, larger active space, we next included the ten p_z orbitals of the carbon atoms in the Cp rings into the target AO list. 
The AVAS construction then yields 15 orbitals above the 0.1 overlap threshold, giving an (18e,15o) active space: 9 of occupied character
and 6 of unoccupied character. 
Of the 9 occupied orbitals, four do not carry Fe 3d character, instead belonging to the π-system of the Cp rings; two more have strong mixing with the Fe (3d_xz, 3d_yz) orbitals, and the final three represent almost-pure Fe 3d_xy, 3d_x^2-y^2, and 3d_z^2 atomic orbitals (with ≥95% orbital weight on the Fe 3d AOs).
Among six unoccupied orbitals, two have no Fe 3d character, two have predominantly Fe (3d_xz, 3d_yz) character, and two have <5% character of Fe 3d_xy and 3d_x^2-y^2 orbitals, respectively.
This agrees with our observations for the smaller (10e,7o) active space; namely, the additional orbitals in the (18e,15o) space are practically non-bonding in character.
Lowering the AVAS threshold from 0.1 to 0.05 gives a (22e,17o) active space, which includes some further effectively non-bonding orbitals beyond the (18e,15o) active space.

We calculated the ground state and singlet and triplet excited states of ferrocene with these active spaces at the CASCI and CASSCF levels. 
Previous studies<cit.> indicate that there are three low-lying d → d singlet transitions (1 ^1E”_2, 1 ^1E”_1, 2 ^1E”_1) and
three low-lying d → d triplet transitions (1 ^3E”_1, 1 ^3E”_2, 2 ^3E”_1).
These d→ d transitions describe excitations from the three non-bonding orbitals (predominantly of 3d_x^2-y^2, 3d_xy and 3d_z^2 metal character,
as described above) to the two antibonding orbitals having mostly (3d_xz, 3d_yz) metal character. All
these excited states have multiconfigurational (but single-excitation) character. 
Note that the E”_1 and E”_2 states are doubly degenerate, thus there are 6 low-lying singlet and triplet excited states.
In the CASSCF calculations we therefore state-averaged over 7 roots and 6 roots in the singlet and triplet manifolds, respectively.

Table <ref> displays the excitation energies, 
compared to  experimental and theoretical numbers from the literature, 
including singly excited configuration interaction (SECI) <cit.>, symmetry adapted cluster configuration interaction (SAC-CI) <cit.>
and time-dependent density functional theory (TD-DFT) calculations <cit.>.
Table <ref>  compares the impact of using different active spaces with different methods.

The performance of CASCI for the smallest (10e,7o) active space is reasonable for most of the excited states, with the 2 ^1E”_1 and 2 ^3E”_1 states being exceptions. These two states have some Rydberg character,<cit.>
and therefore the valence CASCI overestimates these transitions by about 2 eV; this effect can also be seen in the errors of the SECI energies. 


Averaging over all the states in the CASSCF seemed to spread the error over the states, lowering all the energies.
An accurate description of the differential correlation in the 2 ^1E”_1 and 2 ^3E”_1 states thus
requires a dynamic correlation treatment. We find excellent agreement for all states at the CASSCF+NEVPT2 level, 
with a largest error of only 0.21 eV. 

Using the larger (18e,15o) active space, which includes the ligand π-orbitals, worsens the CASSCF excitation energies 
(except for the 2 ^1E”_1 and 2 ^3E”_1 states). However, incorporating dynamic 
correlation through NEVPT2 rebalances the states, improving CASCI and CASSCF energies (except for 1 ^1E”_2 and 1 ^3E”_1) and 
yielding better agreement with experiment and the (10e,7o) CASSCF+NEVPT2 excitation energies.
If we further go from the (18e,15o) active space to the (22e,17o) active space, by additionally including the truncated orbitals with eigenvalues below 0.1, the excitation energies change by less than 0.01 eV. 
Together, these observations indicate that the multiconfigurational character is already well 
converged in the smallest (10e,7o) active space.

As discussed above, a second test of the quality of the AVAS active space is provided by the SVD decomposition 
of the overlap between the CASSCF-optimized active space and the active space initial guess. In the case of ideal coincidence, the SVD eigenvalues should be equal to 1.
The smallest SVD eigenvalues we found were 0.927 for the (10e,7o) active space and 0.906 for the (18e,15o) active space, respectively.
This indicates that here the AVAS provides a stable and accurate initial guess for the CASSCF procedure.







 §.§ B. [Fe(NO)(CO)3]-


We next consider the complex anion [Fe(NO)(CO)3]-, which exhibits catalytic activity in a range of organic reactions, and has been extensively characterized both theoretically and 
experimentally (see Ref. fenoco3 and references therein).
The complex features three-center bonds along both the Fe-N-O axis and between Fe and each pair of CO ligands;
<cit.> its catalytic mode of action exhibits a highly unusual nitrosyl-ligand based oxidation in some cases,<cit.>
and response to photo-activation in other cases.<cit.>
Analysis of the ground-state CASSCF wavefunction and natural orbital occupations indicates that
it has some multiconfigurational character, and that it should be thought of as a Fe(0) species bound via two covalent 
π-bonds to the [NO^-].<cit.>




As in the previous example, we started with a RHF calculation for the singlet ground state, and 
for the simplest active space we chose five 3d orbitals of Fe as the target AOs. 
We used the geometry of Ref. fenoco3_2. 
Using an overlap threshold of 0.1 gives rise to five occupied orbitals and three unoccupied orbitals for the active space; two unoccupied orbitals with 
only 6% weight in the 3d orbital space lie below the threshold for active space inclusion.
Of the three unoccupied orbitals included in the active space, two have (3d_yz,3d_xy) and (3d_xz,3d_x^2-y^2) character, respectively, while the third one has mostly 3d_z^2 character.
The resulting active space is visualized in Fig. <ref>.
Additionally including the nitrogen 2p orbitals into the set of target AOs (with the same threshold) gives a (16e,14o) active space; this adds
three MOs with 2p character, involved in the three-center Fe-N-O bonds, and three involved in N-Fe-CO type bonds.

Unfortunately, there is no gas phase experimental excitation data for this system. However, theoretical vertical excitation energies
from state-averaged CASSCF calculations followed by MRCI+Q, using the def2-TZVPP basis set (omitting g-functions) have previously been reported,<cit.> which we can compare against.

We used CASSCF and NEVPT2 to compute the vertical transition energies
averaging over five singlet and four triplet states, as in Ref. fenoco3_2.
The smallest SVD eigenvalue for the active space overlap with the initial guess for the (10e,8o) active space is 0.806, indicating
that AVAS provides a good initial guess. For the (16e,14o) active space, the lowest SVD eigenvalue decreases to 0.652—apparently adding only the nitrogen 2p orbitals,
without also adding the carbon 2p orbitals, leads to a less balanced active space compared to the (10e,8o) AVAS initial active space.

The vertical excitation energies from CASSCF with the (10e,8o) active space are in better agreement with the CASSCF/MRCI+Q excitation 
energies than with CASSCF results from Ref. fenoco3_2. This indicates that the (10e,8o) active space constructed using AVAS 
provides a more balanced decription of electron correlation than the larger active spaces.

However, the NEVPT2 dynamical correlation treatment significantly raises the obtained excitations energies above the MRCI+Q values.
Similarly, CASSCF calculations with the larger (16e,14o) active space also yield significantly higher excitation
energies than with the (10e,8o) active space, and the CASSCF+NEVPT2 excitation energies in this larger space are also
fairly different from the values obtained from the (10e,8o) active space. This lack of stability with respect to the active space size indicates
that the excited states are not benign electronically, and that their accurate description requires a more sophisticated dynamic correlation treatment beyond 2nd order perturbation theory.
This is supported by the MRCI+Q study in Ref. fenoco3_2, where the (empirical) Q-contribution to the excitation energy is as large as 0.2 eV.

To further substantiate these claims, we also computed CASSCF+NEVPT2 results for the same manually selected (14e,9o) initial active space as described
in Ref. fenoco3_2 (in these calculations, the initial 14 active orbitals were manually selected for Fe d and NO π and π^*-character by visual inspection of KS-DFT/PBE
orbitals computed with the def2-TZVPP basis set; the CASSCF excitation energies thus obtained reproduced the results reported in the supporting information of Ref. fenoco3_2 with better than 0.01 eV accuracy).
By comparing the results of (our) NEVPT2 and (the referenced) MRCI+Q for the same active space of Ref. fenoco3_2, we can separate the effect of the dynamic correlation treatment from the quality of the active space.
We see that CASSCF+NEVPT2 calculations performed with our automatically constructed (10e,8o) active space and the manually selected (14e,9o) active space
show a fair agreement in the case of singlet excited states; however, the CASSCF+NEVPT2 method overestimates the energies of the triplet excited states with the (14e,9o) active space
by 0.5–0.75 eV more than with the (10e,8o) active space, compared to the CASSCF/MRCI+Q transition energies.
Combined, these facts strongly suggest that the approximate NEVPT2 correlation treatment is the primary cause of deviation from the
MRCI+Q reference values, rather than our automatically constructed active space, and that the smaller AVAS provides a more balanced description of this system.






 §.§ C. FeO4^2- 


As our next system, we consider the bare tetraoxoferrate (VI) ion, FeO4^2-.
We assume a tetrahedral FeO4^2- cluster with  an Fe–O distance of 1.660 Å. 

We started with a ROHF calculation for the ^3A_2 ground state. 
Including only 3d orbitals into the target AO set gives a (8e,8o) active space.
Three unoccupied MOs have 34.5% weight in the 3d orbital space, with 3d_yz, 3d_xz, and 3d_x^2-y^2 character,
but are mostly centred on the ligands; the occupied MOs have mostly metal character.
This indicates that some of the low-lying excitations are charge-transfer excitations.

An earlier study <cit.> found that the ground and excited states cannot be described by
a simple Ligand Field Theory d^2 model; rather, they contain superpositions of a large number of configurations, including  
ligand-to-metal excitations. From this it has been argued that it is insufficient
to only  consider  molecular orbitals with Fe 3d character in
the active space to describe excited states. 
Indeed, we find that CASSCF+NEVPT2 calculations with the (8e,8o) active space (generated only with the 3d orbitals in the target AO set)
significantly overestimate the excited states energies, by about ≈ 6500 cm^-1  (≈0.81 eV), compared to experiment.

For this reason, we expanded the target AO list to the five 3d orbitals of Fe and 2p orbitals of all four O atoms.
Using option <ref> to transform the alpha orbitals, our scheme with the 0.1 overlap threshold produced 14 
occupied orbitals and 3 unoccupied orbitals, resulting in a (26e,17o) active space. 

We calculated the vertical excitation energies of FeO4^2-, namely 
transition energies from the ground ^3A_2 state to the first two excited states, ^1E and ^1A_1 (see Table <ref>),
using CASCI and CASSCF (state-averaged over three singlet states and one triplet state) and CASCI+NEVPT2,
comparing to previously reported RASSCF and experimental numbers.
The CASCI calculations significantly overestimate the excitation energies; however, this is significantly improved
by optimizing the orbitals using state-averaged CASSCF.
The smallest SVD singular value for the active space overlap between the initial
and optimized active orbitals is 0.968, indicating that
our initial active orbitals provide a very good guess for the CASSCF procedure and only require
a little relaxation to yield good agrement with experiment. 

Including dynamic correlation by means of NEVPT2 on top of CASCI or CASSCF significantly improves the results.
Note that the difference between our CASSCF excitation energies and those in Ref. feo4m2  with 17 orbitals, obtained with 
the RASSCF method, reflects both the slightly different basis set as well as the truncated CI configuration space in RAS.






 §.§ D. VOCl4^2- 
	


We now consider the oxotetrachlorovanadate(IV) anion, VOCl4^2-.
We use a square pyramidal geometry for VOCl4^2-, as in Ref. Vancoillie (although 
we use a different orientation: the V atom is at the origin, the O atom is on the z axis above the x-y plane and the Cl atoms are below the x-y plane). 

In this complex, vanadium is in a d^1 configuration. 
As in the next example, here the d→ d excitations  do not have much multiconfigurational character. However, it is important
that multireference methods (and their active spaces) provide a balanced description of all states, not just multiconfigurational
ones. The vanadium complex provides a system to test this in an early transition metal (single-reference) problem.

If we choose a set of five AOs, representing the five 3d orbitals of a vanadium atom, 
we obtain five occupied and four unoccupied MOs from the ROHF reference wavefunction 
using the AO-projector option <ref>. One of the occupied MOs 
is a non-bonding 3d_xy atomic orbital, while the
other V 3d AOs strongly mix with the valence orbitals of the oxygen atom and four chlorine atoms. 
This results in four doubly occupied bonding MOs and four anti-bonding MOs which are unoccupied in the ground state. 
The two unoccupied MOs have 69.4% overlap with the 3d orbital space and carry (3d_xz,3d_yz) character,
other two have 54.0% and 70.3% overlap with the 3d orbital space and have (3d_z^2 and 3d_x^2-y^2) character, respectively;
all four unoccupied MOs have mostly metal character.

Using the (9e,9o) active space we calculated the lowest transitions, which are essentially d → d in nature.
There are four possible ligand-field transitions from the highest non-bonding 3d_xy orbital to four unoccupied MOs, thus
in CASSCF we averaged over five doublet states. Table <ref> summarizes the low-energy vertical d → d energies.
The CASSCF method with the small (9e,9o) active space gives an accurate ^2B_1 state, but strongly overestimates 
the other excited states. 
Using NEVPT2 to treat the dynamic correlation on top of CASSCF significantly improves the excited state energies,
resulting in a good agreement with the experimental values and with CASPT2 results  obtained with a (11e,10o) space <cit.>.
This (11e,11o) space is similar to ours, with the addition of the oxygen 2p shell.
 
We also construct a larger (33e,21o) active space, including the 2p orbitals of O and the 3p orbitals of Cl into the target
AO list. In this larger space, to reduce computational cost, we used CASCI+NEVPT2 rather than CASSCF+NEVPT2.
The excited states from CASCI+NEVPT2 with the (33e,21o) active space are also in excellent agreement with the experimental data.
The stability of the CASCI/CASSCF+NEVPT2 excitations with respect to expanding the active space confirms that the correlation
is well converged  by all these treatments.
 
In the CASSCF calculation with the (9e,9o) active space, the smallest SVD eigenvalue for the active space overlap between the initial and the optimized active orbitals amounts to 0.821. This implies that, in this case, AVAS provides a reasonable, but not perfect, initial guess for CASSCF. 

We also used this complex to test and compare options <ref> and <ref>, as described in Sec. <ref>, for constructing the active space with a ROHF reference determinant |Φ⟩.
The excitation energies, calculated with these two options, differ by less than 100 cm^-1 (≈0.01 eV). However, VOCl4^2-'s ground state has only one singly occupied orbital,
and it is possible that larger differences will occur for systems with ground states of higher spin.





 §.§ E. [CuCl4]^2-  
	

We finally consider the  D_4h [CuCl4]^2- complex, with a Cu–Cl bond length of  2.291 Å, as in Ref. Vancoillie.
As in the vanadium system, the d→ d transitions are single-reference in character: this complex provides a late transition metal example.

Using the 3d AOs of Cu as the target AOs and a default cutoff of 0.1,  we obtain only 5 occupied MOs and no unoccupied MOs. This result might surprise at first glance. However, in this case, the antibonding orbitals have ligand character,
and thus there are no unoccupied MOs having more than 5% 3d character. 
The lowest ligand-field transitions arise from the excitation of electrons from the doubly-occupied MOs 
with dominant 3d character to the singly occupied MOs with 3d character.

Despite the fact that the lowest transitions happen mostly within the 3d orbital space,
such a small (9e, 5o) active space is insufficient to describe them at the CASSCF level—because
the nearly filled space leaves no room for electron correlation. 
CASSCF+NEVPT2 however, provides good agreement with the experimental numbers (see Table <ref>).
To see the effect of a larger active space, we also included the
 3p AOs of Cl in the target AO list, obtaining a (33e,17o) active space. The corresponding CASSCF
excitation energies are still poor, indicating that the necessary correlation is not of valence character.
The CASSCF+NEVPT2 excitation energies in this larger space, however, remain in very good agreement with experiment.
The insensitivity to active space indicates that correlations are well converged in the CASSCF+NEVPT2 treatment.

In CASSCF calculations we averaged over five doublet states, each having four doubly occupied d orbitals and one d-orbital with a single electron.
The smallest SVD eigenvalue for the active space overlap with the initial guess in converged CASSCF calculations
is equal to 0.930 in the case of the small (9e,5o) active space and 0.985 for the (33e,17o) active space, indicating
that AVAS provides a good initial guess for CASSCF. 





 §.§ F. Fenton reaction: an example of homolytic bond dissociation

When modeling chemical reactions with multi-configuration methods, one particularly challenging problem is the selection of active spaces capable of representing all relevant configurations in a homolytic bond dissociation process.
In fact, aside from special cases, where e.g. a reasonable active space choice along the reaction path is enabled by “accidental” factors, such as molecular symmetry, or all the relevant orbitals are energetically well-separated from and unmixed with other orbitals, or the molecule is simply small enough to allow for a full-valence active space treatment, active space selection constitutes
the primary bottleneck in the real-world applications of multi-configurational methods to reaction chemistry.





However, as rationalized in Sec. <ref>, we expect the AVAS scheme to be applicable when studying homolytic bond dissociation processes. 
To illustrate this capability, we here consider one key step of a Fenton reaction. Concretely, we consider the homolytic splitting of hydrogen peroxide (H2O2) by aqueous ferrous (Fe(II)) iron at low pH, to yield OH^. radicals and an aqueous ferric (Fe(III)) species. Summarized, we examine the inner step of



Fe^2+ + H2O2 -> [Fe^2+(HOOH)] -> [Fe^3+(OH^-)] + OH^. -> Fe^3+ + OH^. + OH^-.


Variants of this reaction have been proposed as possibly relevant elementary steps in biochemical processes involving iron-oxo (Fe(IV)=O) species (including such as performed by cytochromes P450<cit.>), but this chemistry has been the subject of intense mechanistic debate since its inception in 1894<cit.>.

The history, background, and significance of this reaction, as well as the surrounding controversies, are discussed in Refs. dunford:Iron23VsHydrogenPeroxideReview, kremer:fenton1999, and petit:FentonReaction.

Here no attempt is made at quantitatively describing the reaction, or at resolving any of the controversies it involves; rather, we only examine whether AVAS is capable of producing an active space capable of qualitatively describing the reaction mechanism of a non-trivial homolytic bond dissociation along the entire IRC. In particular, no solvation effects are considered.








To model the Fenton reaction, we considered the [Fe(H2O)5(H2O2)]^2+ complex as an initial reagent. 
It represents an octahedral [Fe(H2O)6]^2+ complex, a model for the aqueous ferrous ion, 
with one water molecule substituted by the hydrogen peroxide molecule. 
Such a complex is expected to be formed when H2O2 enters the coordination sphere of [Fe(H2O)6]^2+.
The homolytic bond cleavage/dissociation of the [Fe(H2O)5(H2O2)]^2+ complex leads to the formation of two spin-coupled radical fragments: ferric [Fe(H2O)5(OH)]^2+ and hydroxyl OH radicals <cit.>. 
For this process we computed a reaction path along the intrinsic reaction coordinate, at the level of UB3LYP/def2-TZVP as described in Sec. <ref> (the geometry optimizations were started with structures from Ref. petit:FentonReaction).
The reaction model is visualized in Fig. <ref>.
The reactant complex has four unpaired electrons on iron and a total spin quantum number of S=4/2.

To build the active space, we first computed ROHF wave functions (with four unpaired electrons) for each structure along the IRC.
To ensure the convergence for the ROHF method and retain a continuous character of the ROHF solution along the entire reaction path, we used the orbitals obtained for the previous geometry as an initial guess for the next one, starting at the reactant side.
The actual active space was then formed by the AVAS procedure, in which we used the three 2p orbitals of two dissociating oxygen atoms and five 3d orbitals of Fe as target AOs.
For the open-shell treatment, we invoked option 3 from Sec. <ref>; that is, we applied the AVAS projection only on doubly-occupied and virtual orbitals from ROHF wavefunction, and the four singly-occupied 3d
orbitals of Fe were added unchanged to the active space to keep the correct spin S number for the entire system.
For the initial complex, AVAS with a 10% threshold produced 12 occupied orbitals and 1 unoccupied orbital to form the (20e,13o) active space.
To retain consistency of the active space along the reaction path, we then fixed the size of the active space to the values thus obtained; that is, for each geometry we chose the 12 occupied orbitals corresponding to the 12 largest overlap eigenvalues with the target space, and the one virtual orbital with the largest overlap with the target space.



With this active space, multi-configuration calculations along the reaction path were executed.
We wish to probe the capability of the AVAS to describe the qualitative features of the electronic structure  on both sides of the reaction, despite the fact that AVAS was constructed from a ROHF wave function, which is itself incorrect on the product side. For this reason, the orbitals were not optimized in the multiconfiguration treatment—that is, we use (approximate) CASCI with ROHF orbitals as a multi-configuration method, and not (orbital-optimized) CASSCF.
Concretely, the CASCI energies were approximated with DMRG-CI (M=1600), which should be essentially exact with this active space.

The obtained energies are displayed in Fig. <ref> (top).
During the process, the ROHF energy rises, because a ROHF determinant is incapable of qualitatively describing the electronic structure of the product complex.
In contrast, with the same ROHF orbitals, the computed CASCI energies show an energy maximum, and a decreasing energy towards the product side.
In Fig. <ref> (bottom) we also plot the CASCI-occupation numbers of two active MOs participating in this homolytic bond cleavage along the IRC.
As expected for a homolytic bond cleavage, in which one doubly occupied molecular orbital splits into two singly occupied molecular orbitals, the occupation numbers go from ≈2.0/0.0 to ≈1.0/1.0 during the process.
This confirms the formation of two radicals on the product side.
Concretely, the active orbital obtained from the virtual side of the AVAS construction at ROHF level changes its occupation
number from 0.06 to 0.96 along the IRC, at the DMRG-CI level.

As expected at this level of theory, the absolute energetics of the process are not well reproduced (compared to accurate CCSD(T)-F12 calculations<cit.>).
However, we conclude that the qualitative features of the process are correctly described with the given active space.
In particular, the AVAS contains the relevant active orbitals for describing the electronic structure of both the reactant side and product side of the reaction—even though the AVAS is constructed from a ROHF wave function which is qualitatively incorrect on the product side.
This observation is consistent with the earlier finding in DMET that relevant entanglement spaces can be constructed from Hartree-Fock wave functions even in cases where the latter are qualitatively incorrect.<cit.>
All geometries, energies and occupation numbers are available in Supporting Information.




§ SPECULATIVE EXTENSIONS FOR COMPLEX CASES

We believe that AVAS is useful for the application cases fulfilling the premises described in Sec. <ref>.
However, it is clear that various application scenarios fall outside of this range.
We here outline two straight-forward extensions of the AVAS concept, which may extend its range of applicability.

Investigations of the details and practical performance of these approaches are beyond the scope of this work and will be done elsewhere.






 §.§ Active spaces beyond the valence space

For some combinations of electronic structure methods and/or application scenarios, it may be desirable to construct active spaces capable of explicitly representing electronic degrees of freedom which cannot be represented within a (minimal) valence-AO basis.
One such case is the treatment of Rydberg excited states: describing such states on even a qualitative level will require diffuse orbitals in the active space.

The other, more common, case, is the extension of the active space in order to improve the quantitative accuracy of computations in certain cases.
For example, the accuracy of CASPT2 (and related PT2 methods) for transition metal complexes is often substantialy improved if a correlating second radial shell of d-AOs is included into the active space. 
This is known as the “double-shell effect”,<cit.> and rests on the fact that doing so, effectively, relocates the treatment of d-electron radial electron correlation from the (very approximate) “PT2 part” of the method to the (very accurate) “CAS part”.

While in this work only valence-spanned active spaces are considered, the AVAS procedure should be able to also create such extended active spaces—by including the additional desired non-valence AOs in the target AO space. This can be done by simply taking any non-minimal atomic natural orbital basis set (such as ANO-RCC,<cit.> ano-pVnZ,<cit.> ANO-VT-XZ,<cit.> or any other) as source for these target AOs—rather than the MINAO basis used here.
For example, a “double-shell” active space for a 3d transition metal atom could be created by including both 3d and 4d AOs in the “target” AO space of AVAS. One could even add an additional shell of tight f functions, to follow the idea of “correlation consistency”<cit.> of dynamic correlation—although doing so is not common practice.





 §.§ Handling complex metal/ligand interactions or multi-nuclear coordination complexes

We expect that normally the simple d-electron AVAS will be sufficient to qualitatively represent the multi-reference wave functions of ground states and (at least) d↦ d-excited states of single-metal coordination complexes. However, there are many application scenarios in which either complicated and non-transparent metal-ligand interactions, or the presence of multiple metal cores in coordination complexes, would make a straightforward active space construction via AVAS either questionable (if ligand AOs or additional metal AOs are not included in the target space) or infeasible (if such additional AOs are included, but the so-created active space becomes too large for convenient quantitative handling).

In the presence of such effects it may be warranted to combine AVAS with either one of several methods of constructing approximate FCI wave functions in large active spaces (vide infra), or with the entanglement-based active-space construction procedure of Reiher and coworkers<cit.>.

In the simplest case of suspected complex metal-ligand interactions in a single-metal complex, one might, for example, first create an AVAS by including as target AOs both the d-electrons of the transition metal, as well as all ligand AOs suspected to play an important role (e.g., all ligand valence AOs from the first coordination sphere, or all p_z-orbtials of a large π-system coordinated to the metal atom).
In practice, the active space created from this choice will frequently be too large for an accurate quantitative multi-reference calculation including dynamic correlation. 
However, there are methods capable of providing qualitative wave functions and first-order reduced density matrices at an approximate CASCI level for such active spaces.
For large active spaces without too many open shells, methods such as SHCI <cit.> or FCIQMC <cit.> can
provide an efficient description, while for active spaces with a larger number of open shells, but fewer than 50 orbitals, QC-DMRG <cit.> is a reliable approach.
Once such a qualitative wave function is present, it can be used to determine which orbitals are not needed in the active space.






Concretely, we suggest the following procedure of handling these complex cases:

   
  * Make an initial AVAS, with a target space including all AOs suspected to be possibly relevant.
   
  * Compute an approximate CASCI wave function |Ψ_approx⟩ in this active space—by using, for example, SHCI<cit.> (with high threshold ϵ) or QC-DMRG<cit.> (with low bond dimension M).
   
  * Compute the natural orbitals {ψ_i} and their occupation numbers {n_i} from the approximate CASCI wave function |Ψ_approx⟩.
   
  * Do the actual quantitative calculation with an active space {ψ_i; t≤ n_i≤ (2-t)}, where t is some numeric threshold (e.g., t=0.05). That is, use an active space composed of the natural orbitals ψ_i from step 3 for which the occupation numbers n_i differ significantly from 0.0 (empty in all configurations of |Ψ_approx⟩) and 2.0 (doubly-occupied in all configurations of |Ψ_approx⟩).

While not investigated here, we expect this combination of methods to be capable of identifying suitable active spaces in many application scenarios.




 §.§ Quantitative treatment of homolytic bond dissociation processes

The obtained data in Sec. <ref> shows that the AVAS contains the necessary orbitals for describing the electron configurations relevant in a homolytic bond dissociation process, along the entire reaction path.
However, we note that for actual quantitative calculations, particularly in simple concerted reaction processes  without hidden intermediates,<cit.> modified procedures may be beneficial.
For example, if the chemical transformation of the reaction is sufficiently characterized by a single transition state (as in Sec. <ref>), a more economical way of treating the reaction might be to compute an AVAS active space only at the transition state, use it to initialize a CASSCF calculation there, and then propagate the initial orbitals from geometry to geometry.
The main benefit of this latter approach would be that in the AVAS procedure, not only the rotated occupied and virtual orbitals with very low overlap with the target AO space could be eliminated from the active space, but also the orbitals with very high overlap (say, ≥98%) could be eliminated, as such values indicate that they stay doubly occupied/unoccupied during the reaction.
The active space size could therefore be substantially reduced in such cases, therefore allowing the application of more powerful electron correlation methods.
However, a detailed study of this is beyond the scope of this work.

We also wish to explicitly note that there are many reactions in which the chemical transformations are not characterized by only the transition state (e.g., see Ref. knizia:CurlyArrows Fig. 3, or Refs. kraka:UrvaReview,joo:MechanismOfBarrierlessReaction,kraka:StunningExampleOfComplexMechanism), and in these cases the ability to construct a consistent active space along the entire reaction path, such as afforded by the AVAS procedure, may be valuable.





























§ CONCLUSIONS


In this work, we investigated how to systematically and automatically construct
molecular active spaces solely from a single determinant wavefunction
together with a list of atomic valence orbitals. The atomic valence active space (AVAS) procedure is based on a straightforward
linear algebraic rotation of the occupied and unoccupied molecular orbital spaces which maximizes their given atomic valence character.
The method automatically detects the valence bonding partners of a given atomic valence orbital, and, by using 
a single small threshold, can also detect non-bonding orbitals without constructing
spurious partners (either occupied or unoccupied).

To test our scheme, we tested both the quality of our orbitals as initial guesses for CASSCF optimization,
as well as the accuracy and stability of the valence excitation energies calculated within our spaces. We find
high overlap of our orbitals with fully optimized CASSCF orbitals, demonstrating their high quality.
We can also obtain good CASSCF excitation energies in cases where the excitations are dominated by valence correlation.
In molecules where the excitations are not of this character, we find that the
addition of dynamic correlation (through the N-electron valence perturbation theory)
yields quantitative agreement with experiment.

No doubt it will still be necessary to experiment with active spaces in the modeling of large and very complex molecules.
Our study  provides two reasons to believe that the difficulty of
performing reliable multireference calculations in complex problems can be reduced  using the AVAS technique.
First,  the simple procedure makes it trivial to obtain not only the minimal active spaces, but 
also extended active spaces, for example, including additional ligand orbitals.
This makes it simpler to systematically explore different active spaces,  eliminating user error and subjectivity in their definition,
and allowing for convergence of properties with respect to the active space size.
Second, systematically varying active space size, while including a dynamic correlation treatment (such as NEVPT2) 
provides a straightforward way to  assess whether our active space is converged, as computed observables should become
insensitive to the active space size.
For these reasons, we believe that the AVAS construction provides a simple route to painless multireference calculations by non-experts,
particularly in complex systems involving transition metals.





§ SUPPORTING INFORMATION

Supporting information: Additional computational details; geometric structures of the studied transition metal complexes; additional figures of studied active orbital spaces; numerical and structural data of the studied reaction path of the Fenton reaction.
This material is available free of charge via the Internet at http://pubs.acs.org.



§ ACKNOWLEDGMENTS

We acknowledge the US National Science Foundation for funding this research primarily through the award NSF:CHE-1665333. Additional support for software development and to support QS was provided through NSF:CHE-1657286. We 
acknowledge additional support for GKC from the Simons Foundation through a Simons Investigatorship.




