Automated construction of molecular active spaces from atomic valence orbitals

TL;DR

The paper presents AVAS, an automated method for constructing active orbital spaces in multireference electronic structure calculations, simplifying the process especially for transition metal chemistry and bond dissociation studies.

Contribution

Introduces AVAS, a straightforward automated technique for defining active spaces based on atomic valence orbitals, improving ease and reproducibility of multireference calculations.

Findings

Accurately calculates excitation energies for transition metal complexes

Effectively follows bond breaking in a Fenton reaction

Simplifies active space determination in MR calculations

Abstract

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 multireference (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.

Tables

Table 1. Table 1: Lowest singlet and triplet excitation energies of ferrocene: results from the (10e,7o) active space compared to other electronic structure methods. Energies are given in eV.

State	CASCI	CASSCF	CASSCF+NEVPT2	SAC-CI 80	SECI 83	TD-DFT 84	Expt.*
	(10e,7o)	(10e,7o)	(10e,7o)			PBE
$1^1{E}_2^{\prime \prime }$	2.85	2.10	2.79	2.11	2.63	2.90	2.885, 2.7 86
$1^1{E}_1^{\prime \prime }$	3.33	2.17	2.87	2.27	3.31	3.03	2.8185, 2.98 86
$2^1{E}_1^{\prime \prime }$	5.82	4.07	3.99	4.03	5.74	3.60	3.82 85, 86
$1^3{E}_1^{\prime \prime }$	1.81	0.97	1.88	1.40	1.81-1.87		1.7485
$1^3{E}_2^{\prime \prime }$	1.84	1.07	2.03	1.68	1.81-1.87		2.0585
$2^3{E}_1^{\prime \prime }$	4.26	2.25	2.50	2.60	4.56		2.29-2.3485

Table 2. Table 2: Lowest singlet and triplet excitation energies of ferrocene: comparison of active spaces. Energies are given in eV.

	CASCI		CASCI+NEVPT2		CASSCF		CASSCF+NEVPT2		Expt.*
	(10e,7o)	(18e,15o)	(10e,7o)	(18e,15o)	(10e,7o)	(18e,15o)	(10e,7o)	(18e,15o)
$1^1{E}_2^{\prime \prime }$	2.85	3.51	3.10	3.20	2.10	2.88	2.79	2.93	2.885, 2.7 86
$1^1{E}_1^{\prime \prime }$	3.33	3.90	3.20	3.21	2.17	3.45	2.87	3.11	2.8185, 2.98 86
$2^1{E}_1^{\prime \prime }$	5.82	5.81	4.24	4.00	4.07	4.61	3.99	4.06	3.82 85, 86
$1^3{E}_1^{\prime \prime }$	1.81	2.46	2.14	2.25	0.97	1.91	1.88	2.09	1.7485
$1^3{E}_2^{\prime \prime }$	1.84	2.48	2.25	2.34	1.07	2.37	2.03	2.26	2.0585
$2^3{E}_1^{\prime \prime }$	4.26	4.33	2.78	2.54	2.24	2.64	2.50	2.50	2.29-2.3485

Table 3. Table 3: Vertical excitation energies from the ground state to the lowest singlet and triplet states of \ce [Fe(NO)(CO)3]-. Energies are given in eV.

State	CASSCF		CASSCF+NEVPT2		CASSCF89	CASSCF/MRCI+Q 89	CASSCF+NEVPT2
	(10e,8o)	(16e,14o)	(10e,8o)	(16e,14o)	(14e,9o)*	(14e,9o)*	(14e,9o)*
${}_{}{}^{1}A_2^{}$	3.14	4.37	4.04	3.76	2.78	3.26	4.13
${}_{}{}^{1}E$	3.58	4.86	4.29	3.98	3.18	3.53	4.29
${}_{}{}^{1}A_1^{}$	3.54	5.17	4.33	4.25	3.22	3.64	4.40
${}_{}{}^{3}A_1^{}$	2.27	2.44	2.63	2.43	1.76	2.32	3.40
${}_{}{}^{3}E$	2.82	3.41	3.51	3.20	2.44	2.96	4.02
${}_{}{}^{3}A_2^{}$	2.89	4.04	3.80	3.82	2.56	3.15	4.39

Table 4. Table 4: Calculated and experimental low-lying excitation energies of \ce FeO4^2-. Energies are given in cm -1 .

State	CASCI	CASSCF	CASCI+NEVPT2	CASSCF+NEVPT2	RASSCF 90	Expt.*
	(26e,17o)	(26e,17o)	(26e,17o)	(26e,17o)	(26,4,0;12,5,0)
${}_{}{}^{1}E$	9882, 9899	7462	6252, 6268	6548, 6550	6300	6209, 6219 91
						6219, 6230 92
${}_{}{}^{1}A_1^{}$	14703	10710	9006	9471	9200	9119 91, 9176 92

Table 5. Table 5: Vertical excitation energies to the lowest doublet excited states from the B 2 2 superscript subscript 𝐵 2 2 {}^{2}B_{2} state ground of \ce VOCl4^2-. Energies are given in cm -1 .

State	CASSCF	CASSCF+NEVPT2	CASCI+NEVPT2	CASPT293	Exp.*
	(9e,9o)	(9e,9o)	(33e,21o)	(11e,10o)
${}_{}{}^{2}B_1^{}$	11736	13044	12123	12667	1160094
${}_{}{}^{2}E$	18650	14537	13722	13620	1370094
${}_{}{}^{2}A_1^{}$	33808	28586	26824

Table 6. Table 6: Low-lying excitation energies of \ce [CuCl4]^2-. Energies are given in cm -1 .

State	CASSCF		CASSCF+NEVPT2		CASPT2 93	Exp.*
	(9e,5o)	(33e,17o)	(9e,5o)	(33e,17o)	(11e,11o)
$1^2{B}_{2g}$	6588	6588	10675	10666	11321	10500 95
$1^2{E}_g$	8727	8728	12832	12947	13379	12800 95
$1^2{A}_{1g}$	9690	9692	14021	14135	14597

Equations

$$|\Psi ⟩=\{{\phi }_1\dots {\phi }_m\}|\mathrm{core}⟩$$

$$\hat{P}$$

$$[{\mathbf{S}}^A{]}_{ij}=⟨i|\hat{P}|j⟩,$$

$${\mathbf{S}}^A\mathbf{U}=\mathbf{U}\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_{N_{\mathrm{occ}}})$$

$$\forall i,j:\sum _k[{\mathbf{S}}^A{]}_{ik}[\mathbf{U}{]}_{kj}=[\mathbf{U}{]}_{ij}{\sigma }_j.$$

$$|i⟩\mapsto |\tilde{i}⟩=\sum _k|k⟩[\mathbf{U}{]}_{ki},$$

$$[\bar{\mathbf{S}}{]}_{ab}^A=⟨a|\hat{P}|b⟩,$$

$${\bar{\mathbf{S}}}^A\bar{\mathbf{U}}=\bar{\mathbf{U}}\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_{N_{\mathrm{vir}}})\iff \forall a,b:\sum _c[{\bar{\mathbf{S}}}^A{]}_{ac}[\bar{\mathbf{U}}{]}_{cb}=[\bar{\mathbf{U}}{]}_{ab}{\sigma }_b,$$

$$|a⟩\mapsto |\tilde{a}⟩=\sum _c|c⟩[\bar{\mathbf{U}}{]}_{ca}.$$

$$|i⟩$$

$$|a⟩$$

$$S_{ij}^A=⟨i⟩\hat{P}j=\sum _{\mu {\mu }^{\prime }}{C}_i^{\mu }{P}_{\mu {\mu }^{\prime }}{C}_j^{{\mu }^{\prime }}$$

$${\bar{S}}_{ab}^A=⟨a⟩\hat{P}b=\sum _{\mu {\mu }^{\prime }}{\bar{C}}_a^{\mu }{P}_{\mu {\mu }^{\prime }}{\bar{C}}_b^{{\mu }^{\prime }},$$

$$P_{\mu {\mu }^{\prime }}=\sum _{p{p}^{\prime }}⟨\mu |p⟩[{\sigma }^{-1}{]}_{p{p}^{\prime }}⟨p^{\prime }|{\mu }^{\prime }⟩.$$

$$|\tilde{p}⟩=\sum _{\mu }|\mu ⟩T_{\mu p},$$

$$\mathbf{P}=\mathbf{ST}({\mathbf{T}}^{†}\mathbf{ST}{)}^{-1}{\mathbf{T}}^{†}\mathbf{S}.$$

$${\mathbf{S}}_{\text{change}}=({\mathbf{C}}_{\mathrm{act}}^{\mathrm{final}}{)}^{†}\mathbf{S}({\mathbf{C}}_{\mathrm{act}}^{\mathrm{initial}})$$

Full text

1 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 1, 2, 3, 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.4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25 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. 26, 27), 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 28, 29 or Kohn-Sham calculations or from correlated calculations, such as MP2 30, 31 or approximate DMRG calculations 32, 27, 33, 34.

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 Iwata35 and Schmidt and coworkers36 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 spaces37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48.

Sec. 2 explains the motivation, background, theory, and implementation of the construction proposed here. Sec. 4 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. 6 describes conclusions and possible implication for future work.

2 Theory

2.1 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 $\{\varphi_{1}\ldots\varphi_{m}\}$. Technically, this means that the wavefunction can be written as a second quantized product

[TABLE]

where $\{\varphi_{1}\ldots\varphi_{m}\}$ denotes a general occupancy wavefunction within the active orbitals, and $|\mathrm{core}\rangle$ denotes a single determinant. It is clear that the active orbitals must span the space of our chosen specific valence atomic orbitals. However, Eq. (1) 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),49, 50, 51, 52, 53, 54, 55, 56 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 references49, 50 (but as described in the appendix of Ref. 51) 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. 2.2 we describe the isolation of entangled orbitals for the active space construction, and Sec. 2.3 will provide technical details and discuss various practical aspects relevant in the active space case. Finally, Sec. 5 outlines how the presented approach may be extended to more complex cases, such as double-shell effects or complex metal/ligand interactions.

2.2 Isolating target-overlapping orbitals for the active space

Let $A=\{\ket{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 $|\Phi\rangle$ 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 DMET49, 50, 57, 52:(a) a SCF wave function $\ket{\Phi}$ 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 wrong50), and (b) we can isolate the part of $\ket{\Phi}$ 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 $\ket{\Phi}$ 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\ldots N_{\mathrm{occ}}$ denote $N_{\mathrm{occ}}$ occupied molecular orbitals (MOs) of $|\Phi\rangle$. The projector $\hat{P}$ onto the space of atomic orbitals in $A$, is given by

[TABLE]

Here $\sigma$ denotes the $|A|\times|A|$ target AO overlap matrix with elements $[\sigma]_{pq}=\langle p|q\rangle$, and $\sigma^{-1}$ its matrix inverse.

Employing these projectors, we construct the first set of active orbitals by rotating $\ket{\Phi}$’s occupied MOs $\{\ket{i}\}$ as follows. First, we calculate the $N_{\mathrm{occ}}\times N_{\mathrm{occ}}$ overlap matrix of occupied orbitals projected onto $\operatorname{span}(A)$, the space of selected target atomic orbitals:

[TABLE]

where $i,j$ are occupied orbital indices. Next we compute the $N_{\mathrm{occ}}\times N_{\mathrm{occ}}$ (unitary) matrix of eigenvectors $[\mathbf{U}]_{ij}$ of $\mathbf{S}^{A}$, such that

[TABLE]

(where $\operatorname{diag}(\ldots)$ denotes a diagonal matrix of the given elements), or, written in component form,

[TABLE]

There are at most $|A|$ non-zero eigenvalues $\{\sigma_{i}\}$, because $\operatorname{span}(A)$ is a $|A|$-dimensional space, and $[\mathbf{S}]^{A}_{ij}$ involves a projection onto it. Furthermore, the eigenvectors $[\mathbf{U}]_{ij}$ of $S^{A}_{ij}$ define a rotation on the occupied orbitals:

[TABLE]

which clearly separates them into two groups: The at most $|A|$ rotated occupied orbitals $\ket{\tilde{i}}$ with $\sigma_{i}\neq 0$, which have non-vanishing overlap with our target atomic orbitals (and which therefore should go into the active space), and the remaining $\ket{\tilde{i}}$ with $\sigma_{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 $\{\ket{\tilde{i}}\}$ in (6) are obtained as a unitary transformation of $\ket{\Phi}$’s original occupied orbitals $\{\ket{i}\}$. Consequently, a determinantal wave function $\ket{\tilde{\Phi}}$ built from the $\{\ket{\tilde{i}}\}$ is physically equivalent (differs by at most a phase factor) from the original determinant $\ket{\Phi}$. That is, so far we have done nothing to $\ket{\Phi}$ 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 $\{\ket{a},1=1\ldots,N_{\mathrm{vir}}\}$ of $\ket{\Phi}$: We form the $N_{\mathrm{vir}}\times N_{\mathrm{vir}}$ projected overlap matrix

[TABLE]

where $a,b$ are virtual orbital indices, then find its unitary matrix of eigenvectors $\bar{\mathbf{U}}$ such that

[TABLE]

and use $\bar{\mathbf{U}}$ to rotate the virtual orbitals via

[TABLE]

Again, the at most $|A|$ of the new virtual orbitals $\{\ket{\tilde{a}}\}$ with eigenvalues $\sigma_{a}\neq 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 $|\Phi\rangle$, we can form the total active space by combining the sets of $\{\ket{\tilde{i}}\}$ and $\{\ket{\tilde{a}}\}$ with non-zero projected overlap eigenvalues. Since the combined set includes all orbitals which have non-vanishing overlap with our target space $\operatorname{span}(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.

2.3 Technical details of the construction

Sec. 2.2 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=\{\ket{p}\}$ chosen and represented? (b) How are the actual rotation matrices $\mathbf{U}$ (eq. (6)) and $\bar{\mathbf{U}}$ (eq. (9)) 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 $\ket{\Phi}$?) We will discuss these questions in the current and next subsections.

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

[TABLE]

where $\mu$ are basis functions from the (large) computational basis set $B_{1}$ (e.g., cc-pVTZ or def2-TZVPP), and $C^{\mu}_{i}=[{\mathbf{C}}_{\mathrm{occ}}]_{\mu i}$ and $\bar{C}^{\mu}_{a}=[{{\mathbf{C}}}_{\mathrm{vir}}]_{\mu a}$ are the coefficients of the basis function $\mu$ in the expansion of the occupied orbital $i$ and virtual orbital $a$, respectively. ${\mathbf{C}}_{\mathrm{occ}}$ and ${{\mathbf{C}}}_{\mathrm{vir}}$ denote the $|B_{1}|\times N_{\mathrm{occ}}$ occupied and $|B_{1}|\times N_{\mathrm{vir}}$ virtual sub-matrices of the $|B_{1}|\times|B_{1}|$ SCF orbital matrix $\mathbf{C}$ (note that $N_{\mathrm{occ}}+N_{\mathrm{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=\{\ket{p}\}$ based on a second auxiliary basis set $B_{2}$. This is a minimal basis set of tabulated free-atom AOs (here MINAO is used47; but other choices, such as subsets of ANO-RCC,58, 59 ano-pVnZ,60 or ANO-VT-XZ61, could be considered; see also Sec. 5.1). This choice leads to simple expressions for the projected overlap matrices

[TABLE]

where the matrix elements of the projector are

[TABLE]

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

• •

Let $A\subset 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 $\mathbf{\sigma}$ with elements $\sigma_{pp^{\prime}}=\langle p|p^{\prime}\rangle$, where $p,p^{\prime}\in A$, as well as $\mathbf{\sigma}$’s inverse matrix, with elements $\sigma^{pp^{\prime}}=[\sigma^{-1}]_{pp^{\prime}}$. Both matrices have dimension of $|A|\times|A|$.

• •

Form the overlap matrix $\mathbf{S}_{21}$ between the functions $p$ of $A\subset B_{2}$ and the functions $\mu$ of the large basis set $B_{1}$, with elements $[\mathbf{S}_{21}]_{p\mu}=\langle p|\mu\rangle$.

• •

Form the projector $P_{\mu\mu^{\prime}}=\sum_{pp^{\prime}}\langle\mu|p\rangle\sigma^{pp^{\prime}}\langle p^{\prime}|\mu^{\prime}\rangle$ , or $\mathbf{P}=\mathbf{S}_{21}^{\dagger}\mathbf{\sigma}^{-1}\mathbf{S}_{21}$.

• •

Form the projected overlap matrices $\mathbf{S}^{A}={\mathbf{C}}_{\mathrm{occ}}^{\dagger}\mathbf{P}\,{\mathbf{C}}_{\mathrm{occ}}$ for the occupied orbitals (eq. (11)), and $\bar{\mathbf{S}}^{A}={\mathbf{C}}_{\mathrm{vir}}^{\dagger}\mathbf{P}\,{\mathbf{C}}_{\mathrm{vir}}$ for the virtual orbitals (eq. (12)).

• •

Finally, diagonalize both projected overlap matrices to obtain the transformation matrices separating the MO sets by overlap with $\operatorname{span}(A)$. Concretely, diagonalize $\mathbf{S}^{A}$ to obtain the eigenvector matrix $\mathbf{U}$, and use it to find the transformed occupied orbital matrix ${\tilde{\mathbf{C}}}_{\mathrm{occ}}={{\mathbf{C}}}_{\mathrm{occ}}\mathbf{U}$. Then diagonalize $\bar{\mathbf{S}}^{A}$ to obtain the eigenvector matrix $\bar{\mathbf{U}}$, and use it to find the transformed virtual orbital matrix ${\tilde{\mathbf{C}}}_{\mathrm{vir}}={{\mathbf{C}}}_{\mathrm{vir}}\bar{\mathbf{U}}$. These are the expansion coefficients of $\{\ket{\tilde{i}}\}$ (eq. (6)) and $\{\ket{\tilde{a}}\}$ (eq. (9)), 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).47 If the IAOs are given as

[TABLE]

where $T_{\mu p}$ denotes the elements of the $|B_{1}|\times|B_{2}|$ IAO transformation matrix,47 this can be incorporated by updating the projection matrix in $S^{A}_{ij}$ and $S^{A}_{ab}$ in eqs. (11) and (12) as

[TABLE]

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. 5.1.

2.4 Truncating the active space

The eigenvalues $\sigma_{i}$ and $\sigma_{a}$ of the projected overlap matrices in eqs. (6) and (9) reflect the degree to which the transformed orbitals $\ket{\tilde{i}}$ and $\ket{\tilde{a}}$ overlap with the space of our target AOs. If we include every such transformed orbital with $\sigma_{a}\neq 0$ and $\sigma_{i}\neq 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 $\sigma_{i}$ and $\sigma_{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 $\operatorname{span}(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 $\ket{\Phi}$, 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. (1). 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. 52 and 57). 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.

2.5 Treatment of open-shell systems

If $\ket{\Phi}$ is a closed-shell determinant, then the active space construction algorithm can be directly used as described in Sec. 2.3. 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 62). A possible remedy for this problem would be to construct a single set of “corresponding orbitals”63 from the separate alpha- and beta-orbital sets, but this has not been tested here. 2. 2.

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 $\ket{\Phi}$ 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 $\ket{\Phi}$ have only small components on the target valence AOs, this can lead to very bad CASCI wavefunctions. 3. 3.

If a ROHF determinant $\ket{\Phi}$ is used, one may apply the construction of Sec. 2.3 exclusively to the doubly-occupied and fully unoccupied orbitals of $\ket{\Phi}$, to form the core determinant and initial part of the active space, and then include additionally all the singly-occupied orbitals of $\ket{\Phi}$ 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 2, however, we compare the different schemes in one of the systems below.

3 Computational Details

We implemented the atomic valence active space (AVAS) construction within the PySCF62 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)64, 65, 66, 67 calculations were carried out using PySCF. For active spaces with more than 16–17 orbitals we used the Block code68 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-DK69, 70 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 basis47, 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,71, 72 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 B3LYP73 with def2-TZVP basis sets74, starting from the structures provided in Ref. 75. 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.76 The reaction path was computed by tracing the intrinsic reaction coordinate77 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.78, 79 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,78, 79 and show iso-surfaces enclosing 80% of the orbital’s electron density.

4 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_{\mathrm{act}}\times N_{\mathrm{act}}$ overlap matrix

[TABLE]

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. 2. 2.

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, \ceFe(C5H5)2. The ground state of ferrocene is dominated by a single configuration with Fe $d^{6}$. Both the MO analysis in Ref. 80 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. 81 (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. 1. 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 82, 80.

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 $\pi$-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 $\geq 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 studies80 indicate that there are three low-lying $d\rightarrow d$ singlet transitions $(1^{1}E^{\prime\prime}_{2},1^{1}E^{\prime\prime}_{1},2^{1}E^{\prime\prime}_{1})$ and three low-lying $d\rightarrow d$ triplet transitions $(1^{3}E^{\prime\prime}_{1},1^{3}E^{\prime\prime}_{2},2^{3}E^{\prime\prime}_{1})$. These $d\rightarrow 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^{\prime\prime}_{1}$ and $E^{\prime\prime}_{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 1 displays the excitation energies, compared to experimental and theoretical numbers from the literature, including singly excited configuration interaction (SECI) 83, symmetry adapted cluster configuration interaction (SAC-CI) 80 and time-dependent density functional theory (TD-DFT) calculations 84. Table 2 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^{1}E^{\prime\prime}_{1}$ and $2^{3}E^{\prime\prime}_{1}$ states being exceptions. These two states have some Rydberg character,80 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^{1}E^{\prime\prime}_{1}$ and $2^{3}E^{\prime\prime}_{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 $\pi$-orbitals, worsens the CASSCF excitation energies (except for the $2^{1}E^{\prime\prime}_{1}$ and $2^{3}E^{\prime\prime}_{1}$ states). However, incorporating dynamic correlation through NEVPT2 rebalances the states, improving CASCI and CASSCF energies (except for $1^{1}E^{\prime\prime}_{2}$ and $1^{3}E^{\prime\prime}_{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. \ce[Fe(NO)(CO)3]-

We next consider the complex anion \ce[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. 87 and references therein). The complex features three-center bonds along both the \ceFe-N-O axis and between \ceFe and each pair of \ceCO ligands; 87 its catalytic mode of action exhibits a highly unusual nitrosyl-ligand based oxidation in some cases,88 and response to photo-activation in other cases.89 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 $\pi$-bonds to the [NO*-*].87

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. 89. 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. 2. 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 \ceFe-N-O bonds, and three involved in \ceN-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,89 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. 89. 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. 89. 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. 89, 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. 89 (in these calculations, the initial 14 active orbitals were manually selected for Fe $d$ and NO $\pi$ and $\pi^{*}$-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. 89 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. 89, 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. \ceFeO4^2-

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

We started with a ROHF calculation for the ${}^{3}A_{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 90 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 $\approx$ 6500 cm*-1* ($\approx$0.81 eV), compared to experiment.

For this reason, we expanded the target AO list to the five 3$d$ orbitals of Fe and 2$p$ orbitals of all four O atoms. Using option 2 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 \ceFeO4^2-, namely transition energies from the ground ${}^{3}A_{2}$ state to the first two excited states, ${}^{1}E$ and ${}^{1}A_{1}$ (see Table 4), 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. 90 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. \ceVOCl4^2-

We now consider the oxotetrachlorovanadate(IV) anion, \ceVOCl4^2-. We use a square pyramidal geometry for \ceVOCl4^2-, as in Ref. 93 (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\to 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 2. 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\rightarrow 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 5 summarizes the low-energy vertical $d\rightarrow d$ energies. The CASSCF method with the small (9e,9o) active space gives an accurate ${}^{2}B_{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 93. 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 2 and 3, as described in Sec. 2.5, for constructing the active space with a ROHF reference determinant $\ket{\Phi}$. The excitation energies, calculated with these two options, differ by less than 100 cm*-1* ($\approx$0.01 eV). However, \ceVOCl4^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. \ce[CuCl4]^2-

We finally consider the $D_{4h}$ \ce[CuCl4]^2- complex, with a Cu–Cl bond length of 2.291 Å, as in Ref. 93. As in the vanadium system, the $d\to 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 6). 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. 2.1, 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 (\ceH2O2) by aqueous ferrous (Fe(II)) iron at low pH, to yield \ceOH^. radicals and an aqueous ferric (Fe(III)) species. Summarized, we examine the inner step of

\ce

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 (\ceFe(IV)=O) species (including such as performed by cytochromes P45096), but this chemistry has been the subject of intense mechanistic debate since its inception in 189497. The history, background, and significance of this reaction, as well as the surrounding controversies, are discussed in Refs. 96, 98, and 75.

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 \ce[Fe(H2O)5(H2O2)]^2+ complex as an initial reagent. It represents an octahedral \ce[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 \ceH2O2 enters the coordination sphere of \ce[Fe(H2O)6]^2+. The homolytic bond cleavage/dissociation of the \ce[Fe(H2O)5(H2O2)]^2+ complex leads to the formation of two spin-coupled radical fragments: ferric \ce[Fe(H2O)5(OH)]^2+ and hydroxyl \ceOH radicals 75. For this process we computed a reaction path along the intrinsic reaction coordinate, at the level of UB3LYP/def2-TZVP as described in Sec. 3 (the geometry optimizations were started with structures from Ref. 75). The reaction model is visualized in Fig. 3. 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 \ceFe as target AOs. For the open-shell treatment, we invoked option 3 from Sec. 2.5; 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. 4 (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. 4 (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 $\approx$2.0/0.0 to $\approx$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 calculations75). 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.50 All geometries, energies and occupation numbers are available in Supporting Information.

5 Speculative extensions for complex cases

We believe that AVAS is useful for the application cases fulfilling the premises described in Sec. 2.1. 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.

5.1 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”,99, 100, 101 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,58, 59 ano-pVnZ,60 ANO-VT-XZ,61 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 3$d$ transition metal atom could be created by including both 3$d$ and 4$d$ 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”102, 103, 104 of dynamic correlation—although doing so is not common practice.

5.2 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\mapsto 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 coworkers26, 34.

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 $\pi$-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 105, 106 or FCIQMC 107, 108, 109, 110 can provide an efficient description, while for active spaces with a larger number of open shells, but fewer than 50 orbitals, QC-DMRG 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122 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. 2. 2.

Compute an approximate CASCI wave function $\ket{\Psi_{\mathrm{approx}}}$ in this active space—by using, for example, SHCI105 (with high threshold $\varepsilon$) or QC-DMRG121 (with low bond dimension $M$). 3. 3.

Compute the natural orbitals $\{\psi_{i}\}$ and their occupation numbers $\{n_{i}\}$ from the approximate CASCI wave function $\ket{\Psi_{\mathrm{approx}}}$. 4. 4.

Do the actual quantitative calculation with an active space $\{\psi_{i};\,t\leq n_{i}\leq(2-t)\}$, where $t$ is some numeric threshold (e.g., $t=0.05$). That is, use an active space composed of the natural orbitals $\psi_{i}$ from step 3 for which the occupation numbers $n_{i}$ differ significantly from 0.0 (empty in all configurations of $\ket{\Psi_{\mathrm{approx}}}$) and 2.0 (doubly-occupied in all configurations of $\ket{\Psi_{\mathrm{approx}}}$).

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

5.3 Quantitative treatment of homolytic bond dissociation processes

The obtained data in Sec. 4 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,123, 124, 125 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. 4), 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, $\geq 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. 79 Fig. 3, or Refs. 123, 124, 125), 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.

6 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.