In this paper, we study the nuclear gradients of heat bath configuration interaction self-consistent field (HCISCF) wave functions and use them to optimize molecular geometries for various molecules. We show that HCISCF nuclear gradients are fairly insensitive to the size of the “selected” variational space, which allows us to reduce the computational cost without introducing significant errors. The ability of the HCISCF to treat larger active spaces combined with the flexibility for users to control the computational cost makes the method very attractive for studying strongly correlated systems, which require a larger active space than possible with a complete active space self-consistent field. Finally, we study the realistic catalyst, Fe(PDI), and highlight some of the challenges this system poses for density functional theory (DFT). We demonstrate how HCISCF can clarify the energetic stability of geometries obtained from DFT when the results are strongly dependent on the functional. We also use the HCISCF gradients to optimize geometries for this species and study the adiabatic singlet–triplet gap. During geometry optimization, we find that multiple near-degenerate local minima exist on the triplet potential energy surface.

## I. INTRODUCTION

Systems that contain transition metals, covalent bond breaking, and electronically excited states are often challenging for theoretical methods because their electronic structure can be dominated by more than one electronic configuration. To study such systems, we require multireference approaches such as the complete active space self-consistent field (CASSCF) method.^{1–3} In the CASSCF, we perform the full configuration interaction (FCI) procedure on a subset of the molecular orbitals termed the active space. The FCI expansion enumerates all possible configurations for a given set of orbitals and scales combinatorially with the size of the active space. This steep cost limits the size of active spaces that can be treated to roughly 22 electrons in 22 orbitals, which we abbreviate as (22e,22o).^{4} A substantial amount of research in quantum chemistry is devoted to reducing the cost of configuration interaction (CI) methods through a variety of approximations: restricted active space,^{3,5} generalized active space,^{2,6,7} density matrix renormalization,^{8–29} selected configuration interaction (SCI),^{30–49} and FCI quantum Monte Carlo.^{50–53} The heat bath configuration interaction (HCI) is a particularly efficient implementation of the SCI method.^{54–56} These methods can correlate more than 40 electrons in 40 orbitals in routine calculations, making the study of larger and more complex systems accessible.

With the ability to treat larger active spaces, extending these algorithms to calculate molecular properties other than single-point energy is vital for connecting theoretical and experimental research. One such property of interest is the nuclear gradients of the electronic energy, which are crucial for chemical applications because they are used to obtain minimum energy geometries, transition states, and reaction pathways, and are used in *ab initio* molecular dynamics.^{57} Many of the approximate CI/CASSCF schemes have already been extended to calculate the nuclear gradients,^{58–68} and in this work, we extend the family of HCI methods.

The recently developed HCI algorithm has proved to be an efficient approximation of FCI where the size of the CI expansion is controlled with a single user parameter *ϵ*_{1}. This parameter allows one to continuously grow the wave function from a single determinant, e.g., Hartree–Fock (HF), to all configurations in the FCI expansion. The energy of this wave function can be corrected with a perturbative correction, which is often necessary to reduce the error, relative to FCI, to an acceptable level. In recent work, two of the authors extended the HCI algorithm to be used in CASSCF-like calculations and implemented this method in the open source quantum chemistry package PySCF.^{69–72} This work showed that the multireference orbitals can be obtained from HCISCF with relatively loose thresholds and as long as a final step in which the optimized orbitals are used to perform a single tight threshold HCI calculation, the final results are in good agreement with the full CASSCF calculation.^{56} In practice, this insensitivity of the final results to the optimized orbitals is advantageous because it means that accurate multireference orbitals can be obtained from relatively cheap HCI calculations. In general, other properties may also be insensitive to the accuracy of the HCI wave function, which motivates this work and the study of nuclear gradients.

Previous work on approximate CI/CASSCF schemes has shown that for most systems, properties such as gradients and relaxed geometries are often quite accurate even when the energy is not.^{58–68} In particular, Park has studied the gradients of a different selected CI scheme, adaptive sampling CI (ASCI), and demonstrated that CASSCF gradients and geometries were well approximated by this method.^{67,68} Park showed that geometries converged faster than energies with respect to the number of determinants in the wave function and that the convergence of both energy and relaxed geometries could be improved by a second order perturbation theory correction (PT2) and by extrapolating these corrections.^{68} As a result, important properties of these molecules such as singlet–triplet (S-T) gaps were also improvable with these corrections, but we note that for most systems, the difference after the correction was typically around 1 kcal/mol. However, Park’s work focused largely on aromatic systems, which are not particularly strongly correlated, while in this work, we focus on the transition metal complex Fe(PDI), which is highly multireference.^{73}

The rest of this paper is organized as follows: First, we briefly review HCI theory and its self-consistent variant along with analytical gradients of CI wave functions. Next, we discuss the gradients of CASSCF-like wave functions where HCI is used as the CI solver and verify that they converge smoothly with the HCI parameter *ϵ*_{1}. Then, we highlight several strategies to improve the quality of the gradients when increasing the size of the CI expansion is not feasible, e.g., when memory is limited. Finally, we demonstrate the use of these gradients to find equilibrium geometries and adiabatic singlet–triplet gaps for the model catalytic system (PDI)Fe–N_{2}, where PDI is 2,6-bis[1-(2,6-dimethylphenyl-imino)ethyl]pyridine. For simplicity, we refer to this complex as Fe(PDI) and show it in Fig. 1. For this challenging system, we also examine the important practical consideration of selecting the active space in a robust numerical manner.

## II. THEORY

### A. Heat-bath configuration interaction algorithm

HCI is a two-step procedure consisting of a variational calculation and a perturbative correction. In the first stage of the HCI algorithm, a set of *important* determinants is iteratively built and used to expand the variational wave function, then a second-order correction to the energy is computed using multireference Epstein–Nesbet perturbation theory (PT).^{75,76} Like other selected CI (SCI) + PT methods, it derives from FCI and linearly parametrizes the wave function |Ψ⟩ = *∑*_{i}*c*_{i}|*D*_{i}⟩. Unlike FCI, where the expansion is a set of all possible determinants, SCI methods express the wave function as a set of *important* determinants. This allows SCI+PT methods to correlate much larger systems than possible with FCI.

#### 1. Variational stage

In the first step, the HCI algorithm identifies a set of *important* determinants that approximate the full many-body wave function and we refer to this set as the variational space $(V)$. At the beginning of the variational step, the user specifies a small number of determinants to add to $V$, and usually, this includes only the Hartree–Fock determinant. Then, the algorithm adds all connected determinants |*D*_{a}⟩, which satisfy the following importance criteria:

where |*D*_{i}⟩ is a determinant already in the variational space, $Hai=\u27e8Da|H\u0302|Di\u27e9$, *c*_{i} is the amplitude of the |*D*_{i}⟩, and *ϵ*_{1} is a user-specified parameter that controls the size of the variational space, i.e., the accuracy of the variational wave function. A smaller *ϵ*_{1} means that $V$ is closer to the FCI determinant space and in the limit that *ϵ*_{1} goes to 0, we recover our FCI wave function and energy. We repeat this selection of important determinants and grow the variational space until the HCI energy converges.

The CIPSI selection criteria, $\u2211iHkiciE0\u2212Hkk>\u03f5CIPSI$, is more complicated and uses first order perturbation theory to determine whether a determinant should be added to the variational space. Despite its similarity to the CIPSI selection criteria, the HCI importance function offers several key advantages over the *original* CIPSI algorithm. Since the matrix elements between determinants connected by double-excitations are just the two-electron integral elements (and a parity factor of +1 or −1), evaluating the HCI selection criteria requires no calculations. While avoiding calculations is advantageous, the real gain in performance comes from the fact that the two electron integrals are stored and sorted *once* before the variational stage. Selecting *important* determinants only requires iterating through the sorted two-electrons integrals and the algorithm stops once those values fall below a certain threshold. As a result, this means that the algorithm never even “touches” the unimportant connected determinants, which often outnumber the *important* ones by many orders of magnitude. We note that HCI was inspired by challenges in the original formulation of CIPSI and more modern implementations have shown substantial improvements compared to their predecessors.^{77}

#### 2. Perturbative stage

We then correct the variational energy using a second order correction from Epstein–Nesbet perturbation theory, evaluated as

where the inner summation is “screened” in a similar manner to the variational stage. Analogously, the perturbative expansion grows as *ϵ*_{2} becomes smaller, and in the limit *ϵ*_{2} → 0, the perturbative correction becomes exact.

While screening the inner sum in Eq. (2) mitigates the large memory requirement for calculating a perturbative correction, it does not eliminate it entirely. Calculating an accurate correction requires a small *ϵ*_{2} value and for large systems, evaluating Eq. (2) deterministically becomes prohibitively expensive. In these cases, we can reduce the memory requirement further by calculating *E*_{2}(*ϵ*_{2}) semi-stochastically. By that, we mean that we calculate as much of the correction deterministically as computational resources allow and then calculate the remainder stochastically. This procedure is advantageous compared to a purely stochastic one because the portion of *E*_{2}(*ϵ*_{2}) calculated stochastically is smaller, and it takes fewer iterations to converge the stochastic error. We refer to calculations that use this procedure as semi-stochastic HCI (SHCI) to distinguish them from HCI calculations that only use the variational step. There is another even more efficient scheme, which was recently introduced by one of the authors of this paper, but we do not use that scheme in this work.^{78}

### B. Gradients of HCISCF wave functions

We briefly review analytical energy gradients of variational wave functions. For a detailed discussion on the subject, we refer to the reader to the review by Yamaguchi and Schaefer.^{79} We start with the electronic Hamiltonian in the second quantization,

We define the one- and two-body integrals in the usual way,

where *ϕ* are orbitals, *Z*_{I} are the nuclear charges, *r*_{I} the electron–nuclear separations, and *r*_{12} the electron electron separation. We use *I*, *j*, *k*, and *l* as orthogonal orbital indices and *σ*, *σ*′ as spin indices. In most cases, these are molecular orbitals (MOs), which are linear combinations of atomic orbitals (AOs) $|i\u3009=\u2211\mu AOC\mu i|\mu \u3009$. For the remainder of this work, we use the following convention for orbital indices: Greek letters, e.g., *μ* and *ν* refer to AOs, while *i*, *j* refer to MOs.

The active space Hamiltonian is a function of the basis set, the MO coefficients **C**, and the nuclear coordinates {*a*_{i}}, while the electronic wave function |Ψ⟩ is a function of the parameters {*c*_{i}} making the electronic energy a function of all four. Expanding the nuclear gradients in terms of these variables, we get

The first term is the contribution from the basis set, i.e., the atomic orbitals, the second term is from the MO coefficients, and the final term is from the wave function parameters. Since the CI coefficients are always determined variationally so $\u2202Eelec\u2202ci=0$ the term always vanishes. Similarly, for wave functions where the molecular orbital parameters are variational, $\u2202Eelec\u2202C\mu i=0$ and the second term on the right-hand side of Eq. (5) also vanishes. To derive an explicit form for the gradients, we start by expressing the electronic energy in the second-quantized form:

where $\gamma ij=\u2211\sigma \u27e8\Psi |ai\sigma \u2020aj\sigma |\Psi \u27e9$ is the one-body reduced density matrix (RDM) and $\Gamma ijkl=12\u2211\sigma \sigma \u2032\u27e8\Psi |ai\sigma \u2020aj\sigma \u2032\u2020ak\sigma \u2032al\sigma |\Psi \u27e9$ is the two-body RDM. Taking derivative of Eq. (6) with respect to a general nuclear coordinate *a* generates terms, including $\u2202\gamma ij\u2202a$ and $\u2202\Gamma ijkl\u2202a$. Since the RDMs are just functions of the CI coefficients, we know that we can ignore their derivatives because the final term in Eq. (5) vanishes. As a result, we arrive at a concise expression for then the nuclear gradients of the energy,

Expanding this equation further, we can separate the terms relating to the MO response,

where $hija$, $gijkla$, and $Sija$ are the “skeleton” derivatives of the one-electron, two-electron, and overlap integrals defined as

$Uija$ is called the orbital response and is defined by

It can be determined by solving the z-vector equations.^{80} Finally, *X*_{ij} is the Lagrangian from the Generalized Brillouin Theorem, defined as

For the CASSCF, *X*_{ij} is symmetric, i.e., *X*_{ij} = *X*_{ji}, and we can simplify Eq. (8) to

and avoid solving the z-vector equations. We direct the reader to the derivation of Eq. (133) in the excellent review by Yamaguchi and Schaefer.^{79} All three terms in Eq. (12) arise from the response of the basis set to changes in nuclear coordinates, and no terms from the response of MO or CI coefficients are present.

This derivation hinges on the fact that *X*_{ij} is symmetric and for approximate CASSCF schemes such as HCISCF and DMRGSCF, *X*_{ij} is generally not symmetric and since we ignore the third term in Eq. (12), it is not exact. In practice, one can solve include rotation among the active space orbitals, which we refer to as active–active (AA) rotations during the orbital optimization to eliminate the error of using Eq. (12) with approximate schemes such as DMRGSCF and HCISCF. Despite the formal lack of exactness, we highlight that if these approximate CI methods are sufficiently converged, the expression for the gradients is a justified approximation. If these methods are not converged to the FCI limit, Eq. (12) must include a response term due to the AA rotations to be formally exact.

Finally, we note that screening in Eqs. (1) and (2) is dependent on the geometry, meaning that the potential energy surface can be discontinuous for a fixed value of *ϵ*_{1}. During the finite-difference calculations used to verify the analytic gradients, we only observed discontinuities for extremely inaccurate HCI wave functions, where the HCI parameter *ϵ*_{1} is greater than 10 mhartree and the wave functions consisted of less than 100 determinants. Since typical HCI calculations use much smaller values of *ϵ*_{1} and many more determinants, we do not explore these discontinuities further in this work.

## III. COMPUTATIONAL DETAILS

For all DFT geometry optimizations, we use Gaussian16 with the def2-TZVP basis set and density fitting to reduce the computational cost.^{81} We use the default convergence thresholds in the Gaussian16 package. In addition to finding optimal geometries, we used harmonic frequency analysis to confirm that the reported geometries are all genuine local minima, i.e., contain no imaginary frequencies. We perform all other quantum chemical calculations using PySCF^{69–72} and checked several of the mean-field calculations using Q-Chem.^{82} For HCISCF calculations we use our own implementation of the HCI algorithm, called Dice,^{56} in tandem with PySCF. The nuclear gradients are calculated in PySCF using the RDMs obtained from Dice and geomeTRIC^{83} is used to perform the geometry updates. Due to the challenging nature of the Fe(PDI) problem, we relaxed the geometry convergence thresholds for all optimizations of the Fe(PDI) species, which are each double the default thresholds used by both Gaussian16 and geomeTRIC. Our convergence thresholds are 2 · 10^{−6} hartree for the energy (*τ*_{E}), 6 · 10^{−4} hartree/bohr for the root mean squared (rms) of the nuclear gradient (*τ*_{grms}), 9 · 10^{−4} hartree/bohr for the maximum of the nuclear gradient (*τ*_{gmax}), 2.4 · 10^{−3} Å for the rms of the atomic displacement (*τ*_{drms}), and 3.6 · 10^{−3} Å for the maximum atomic displacement (*τ*_{dmax}).

For the Fe(PDI) species, we use a numerical procedure to select the active space to reduce the heuristics which typically plague active space calculations. All multireference calculations use the cc-pVDZ basis set.^{84} We outline the procedure for selecting the active space orbitals below:

Run unrestricted DFT calculation using the PBE0 functional

^{85,86}to generate an initial guess for Step 2.Run an unrestricted Hartree–Fock (UHF) calculation to generate spin orbitals without correlation from an exchange-correlation functional.

Use the procedure outlined in Ref. 87 to generate unrestricted natural orbitals (UNOs).

^{88}Use the UHF natural orbital occupation numbers (NOONs) along with the spin density information from the UHF calculation to ensure that the targeted state shows a strong correlation and qualitatively correct spin density.

Use the UNOs in a “loose” HCI calculation where the

*ϵ*_{1}= 10^{−4}hartree and the active space is (100e,100o). We select the orbitals based on their NOON so the 50 orbitals below the highest occupied natural orbital (HONO) and 50 above the lowest unoccupied natural orbital (LUNO) are part of the active space.Use the 1-RDM from the previous step to calculate new NOs and sort based on NOONs.

Use these orbitals in all subsequent HCISCF calculations, selecting the orbitals solely based on their NOONs analogous to the way used in Step 5. We select the orbitals in pairs (i.e., one occupied and one virtual) around the HONO-LUNO gap. For all Fe(PDI) calculations, we use an active space of (40e,40o), which includes more than those recommended by the criteria UNO methods, i.e., 0.02 < NOON < 1.98.

^{88}

We emphasize that Steps 1 and 2 are sensitive to the optimization strategy used for the SCF procedure, e.g., DIIS,^{89} ADIIS,^{90} or Newton’s Method. We surveyed a range of choices for these strategies in both steps and select the best pair of strategies for each species. By “best” here, we mean that the SCF solution is the lowest in energy, stable with respect to its variational parameters, and has a non-negligible spin-density on the Fe atom. All results from this survey are shown and described in Appendix C. While we have tried to search thoroughly for the lowest energy UHF solution, we note that it is possible that even lower energy solutions exist in this complicated energy landscape. From our experience, these different SCF solutions could yield different multireference orbitals, i.e., distinct multiconfigurationalself-consistent field (MCSCF) wave functions, and as a result different optimized molecular structures.

For all geometry comparisons, we use Kabsch’s algorithm^{91} to align geometries and report the difference in geometry as $\Vert \Delta \Vert 2na$ where Δ is the difference in aligned nuclear coordinates and *n*_{a} is the number of atoms. For brevity, we refer to this quantity as RMSD for the remainder of the paper. To align the geometries we use the open source package rmsd.^{92}

The input scripts necessary to reproduce **all** of the results from this work are publicly available on GitHub at https://github.com/jamesETsmith/hci_nuc_gradients.

## IV. RESULTS

### A. HCISCF gradients

In this subsection, we discuss the errors that arise when using Eq. (12) for approximate CASSCF wave functions. In Fig. 2, we compare different variants of HCISCF gradients to CASSCF gradients, all using the cc-pVDZ basis set.^{84} As mentioned earlier, the HCISCF gradients approximate CASSCF gradients because HCI is an approximation of FCI and as a result, the energy is not invariant for AA orbital rotation. In our discussions of gradients, we note that all gradients discussed suffer from basis set incompleteness errors, but since that is constant throughout the calculations, we do not discuss it further.

Without solving the response term, there are two ways that we can eliminate the error of this approximation: we can include AA rotations in the MCSCF optimization or we can converge the HCI wave function to the FCI limit. Including the AA rotations in the MCSCF optimization removes the need for the z-vector equation completely and ensures that Eq. (12) is exact. However, it can significantly increase the number of MCSCF iterations required to reach convergence. AA rotations make the MCSCF procedure more challenging because the corresponding parameters are strongly coupled to the CI coefficients. Due to these challenges, optimizing AA rotations is often the most expensive way to improve the gradients. In the second approach, we can either increase the size of the variational space or include a perturbative correction and increase the size of the perturbative space to better converge an HCI wave function. By better converging to the FCI limit, we eliminate any effect of AA rotations on the energy, and as a result, any need for the z-vector equation, again ensuring that Eq. (12) is exact. It is more cost-efficient to increase the number of configurations treated perturbatively, and if we are memory limited, then this may be the only option. When performing CASSCF calculations using only the variational step of HCI, we refer to these as vHCISCF wave functions and when we use both the variational and the perturbative, we refer to them as HCISCF. In Fig. 2, we compare different flavors of HCI-based gradients for four systems: N_{2}, Sc_{2}, *trans*-stilbene, and bis-(ethylene-1,2-dithiolato)nickel [Ni(edt)_{2}] using actives spaces of (10e,8o), (6e,18o), (14e,14o), and (18e,13o), respectively. The latter two systesms are shown in Fig. 3. For N_{2} and Sc_{2}, these active space sizes correspond to the full valence space. These active spaces correspond to the full valence space for N_{2} and Sc_{2}, the conjugated 2*p*_{z} orbitals on carbon for *trans*-stilbene, and the conjugated 2*p*_{z} orbitals of carbon and sulfur as well as the 3*d* orbitals of the nickel atom for Ni(edt)_{2}. Schlimgen and Mazziotti have shown that Ni(edt)_{2} exhibits a strong multireference character making it an ideal candidate to study gradients of approximate CASSCF wave functions.^{63,93} The geometries for all systems studied are in the publicly available Github repository at https://github.com/jamesETsmith/hci_nuc_gradients. We start with an intentionally inaccurate vHCISCF wave function, to leave considerable room for improvement, and show that adding AA rotations or PT correction can decrease the gradient error. As expected, we find that AA rotations and the PT correction improve the quality of gradients and with the exception of N_{2} provide the largest reduction of error when used together.

Although AA rotations and the PT correction in our MCSCF optimization can reduce the error in our nuclear gradients, they do so at a cost that is often difficult to predict. Both options increase the number of MCSCF iterations required for convergence and with AA rotations that the calculations can require more than 50 iterations and may not converge at all. This behavior is consistent with the previous DRMG-SCF^{94} and ASCI-SCF^{67,68} work but could likely be improved by using more sophisticated optimization like the ones described by Yao and Umrigar or Kreplin *et al.*^{94,95}

Another strategy is to use a smaller *ϵ*_{1} value, which increases the number of determinants in the HCI wave function. Using *ϵ*_{1}, we can control the accuracy of our gradients, compared to the CASSCF, and we find that we can use a relatively large value of *ϵ*_{1}, i.e., a smaller number of determinants in our wave function, and still produce accurate gradients. For the same two systems as above, we examine the convergence of the analytic gradients with respect to *ϵ*_{1}, i.e., the size of the variational space. Again, we use the cc-pVDZ basis set, and Fig. 4 shows the results for both systems. To contextualize the results, we show the loose, default, and tight gradient convergence tolerances used by the Gaussian16 package^{81} during geometry optimizations. In almost all cases, we find that AA rotations decrease the error relative to the CASSCF gradients, but, in several instances, such as *ϵ*_{1} = 5 · 10^{−4} hartree in *trans*-stilbene, it does not. Since AA rotations are not guaranteed to reduce the vHCISCF or HCISCF gradient error relative to the CASSCF, it is not surprising that some cases exist where it increases the error slightly. For all wave functions, we observe a steady convergence of the gradient to the exact (CASSCF) value. The rate of convergence varies but shows similar behavior for all four systems, which is encouraging given the multireference nature of Sc_{2} and Ni(edt)_{2}. For the typical range of *ϵ*_{1} values used in research applications (*ϵ*_{1} = 10^{−4}–10^{−5} hartree), the HCISCF gradients agree well with their CASSCF counterparts. These encouraging results allow us to reduce the computational time required to calculate nuclear gradients of vHCISCF without significantly reducing their accuracy. For the remainder of this paper, all gradients are calculated with vHCISCF and no AA rotations, unless otherwise specified.

### B. Geometry optimization of Fe(PDI)

We use two geometries reported by Ortuño and Cramer, which they obtained using DFT: the first, “^{1}BS(1,1),” which we call **A**, optimized for the singlet spin symmetry, and the second, “^{3}BS(3,1),” which we call **C**, optimized for the triplet.^{73} Both structures were calculated using unrestricted M06-L with the def2-TZVP basis and density fitting. In this work, we use the shorthand of (M,**X**) to describe a state with multiplicity M at geometry **X**.

Ortuño and Cramer note in their paper that “although we report DFT energies, we note that they are likely quantitatively unreliable owing to significant spin contamination and unpredictable sensitivity to [multireference] character.”^{73} Given this warning, we first studied Fe(PDI) with several DFT functionals to assess how sensitive properties such as energy and spin density (on the Fe atom) were to these choices. We chose several of the Minnesota functionals (M06-2X, M06-L, and MN15)^{97–98} because of their success with challenging transition metal complexes. We also optimized geometries to test how sensitive the relaxed structures were to the choice of these functionals.

For all optimizations, we used **A** as the starting geometry for singlets and **C** as the starting geometry for triplets. We found that for each multiplicity, the three functionals led to distinct geometries, i.e., six geometries in all. We compare the geometries in detail in Appendixes A–E, see Table VI. In Fig. 5, we illustrate the difference in relative energy and Fe spin density at the relaxed geometries. The left panel demonstrates that the spin symmetry of the ground state depends on the functional, and the right panel highlights that spin density on the iron atom ranges from two to four unpaired electrons. We also found that the NOONs at each of these relaxed structures in Table I vary noticeably for different functionals. This qualitative difference makes it challenging to determine whether singlet and triplet geometries should start from **A** or **C** in subsequent multireference calculations. Throughout the study of this species, we have found that geometry optimization use of DFT and vHCISCF wave functions can depend on the initial geometry, which prompted us to select our initial geometries with care.

Functional . | Mult. . | Geom. . | HONO-1 . | HONO . | LUNO . | LUNO+1 . |
---|---|---|---|---|---|---|

M06-2X | 1 | A | 1.325 | 1.000 | 1.000 | 0.675 |

M06-2X | 3 | C | 1.188 | 1.000 | 1.000 | 0.812 |

M06-L | 1 | A | 1.921 | 1.000 | 1.000 | 0.079 |

M06-L | 3 | C | 1.566 | 1.000 | 1.000 | 0.434 |

MN15 | 1 | A | 1.617 | 1.000 | 1.000 | 0.383 |

MN15 | 3 | C | 1.686 | 1.000 | 1.000 | 0.314 |

Functional . | Mult. . | Geom. . | HONO-1 . | HONO . | LUNO . | LUNO+1 . |
---|---|---|---|---|---|---|

M06-2X | 1 | A | 1.325 | 1.000 | 1.000 | 0.675 |

M06-2X | 3 | C | 1.188 | 1.000 | 1.000 | 0.812 |

M06-L | 1 | A | 1.921 | 1.000 | 1.000 | 0.079 |

M06-L | 3 | C | 1.566 | 1.000 | 1.000 | 0.434 |

MN15 | 1 | A | 1.617 | 1.000 | 1.000 | 0.383 |

MN15 | 3 | C | 1.686 | 1.000 | 1.000 | 0.314 |

#### 1. Selecting initial geometries

Since we hope to calculate accurate S-T gaps, which we can qualitatively compare to experiments, it is challenging to “decide,” which geometries and functionals that we should use for multireference calculations. To address this, we can use SHCI energies to determine which DFT geometry is the most stable in an unbiased manner and use these geometries in the subsequent study of the S-T gaps. For each geometry (**A** and **C**), we used the procedure outlined above to select the active space of size (40e,40o) and then ran single point vHCISCF calculations followed by a single tight SHCI calculation. We deliberately correlate a larger number of orbitals than would be selected by the NOON numerical thresholds mentioned earlier in order to confidently include all strongly correlated orbitals in the active space. For the vHCISCF calculations, we used a value of *ϵ*_{1} = 7.5·10^{−5} hartree. Given the sensitivity of the MCSCF solutions to the input orbitals, we tested several values of *ϵ*_{1} and found that the MCSCF orbitals converge faster than the energy with respect to this parameter agreeing with the previous work by several of the authors.^{56} We report this analysis in Appendix D. To ensure that our choice of *ϵ*_{1} was appropriate, we tested the use of a more accurate *ϵ*_{1} value of 5 · 10^{−5} hartree and found that it produced geometries that are quite similar to those produced when using *ϵ*_{1} = 7.5·10^{−5}. The RMSD value between the two optimized geometries is 5.6·10^{−3} Å, which is less than the loose convergence criteria used by Gaussian (6.7·10^{−3} Å). vHCISCF wave functions using *ϵ*_{1} = 7.5·10^{−5} contained roughly 3.5·10^{6} determinants and those using *ϵ*_{1} = 5·10^{−5} contained roughly 7.7·10^{6}. The ability to use a larger value of epsilon reduced the time necessary to solve the CI problem from more than 1600 s to roughly 700 s. We point out that the savings would increase when compared to smaller values of *ϵ*_{1} and this speed-up is a rough lower bound. All timings reported here were recorded on a single AMD Rome node with 128 cores. For the final SHCI calculation, we used several values of *ϵ*_{1} to extrapolate to the FCI limit, the tightest one being 10^{−5}, and *ϵ*_{2} = 5·10^{−8} hartree. We detail this extrapolation procedure in Appendix A. Extrapolations for singlet and triplet states at all three geometries are shown in Figs. 6 and 7). The numerical data are shown in Table II.

Mult. . | Geom. . | Initial (hartree) . | Final (hartree) . | Relaxation (kcal/mol) . |
---|---|---|---|---|

1 | A | −2497.697 79 | −2497.735 52 | −23.67 |

1 | B | −2497.696 78 | −2497.735 50 | −24.29 |

1 | C | −2497.671 40 | −2497.719 16 | −29.97 |

3 | A | −2497.644 96 | −2497.676 56 | −19.83 |

3 | B | −2497.680 06 | −2497.709 22 | −18.30 |

3 | C | −2497.683 59 | −2497.710 49 | −16.88 |

Mult. . | Geom. . | Initial (hartree) . | Final (hartree) . | Relaxation (kcal/mol) . |
---|---|---|---|---|

1 | A | −2497.697 79 | −2497.735 52 | −23.67 |

1 | B | −2497.696 78 | −2497.735 50 | −24.29 |

1 | C | −2497.671 40 | −2497.719 16 | −29.97 |

3 | A | −2497.644 96 | −2497.676 56 | −19.83 |

3 | B | −2497.680 06 | −2497.709 22 | −18.30 |

3 | C | −2497.683 59 | −2497.710 49 | −16.88 |

Figure 6 indicates that **A** is the best structure to start geometry optimization for the singlet spin symmetry and **C** for the triplet, agreeing with the DFT results using M06-L and MN15. For completeness, we optimized the Fe(PDI) structure starting from both multiplicities at all initial geometries for a total of four vHCISCF geometry optimizations. For all vHCISCF gradients used in geometry optimizations, we use Eq. (12) without including any AA rotations or response terms mentioned in Sec. II making the gradients approximate. This approach is consistent with several previous studies where gradients of approximate CASSCF methods were calculated.^{59–61,67,68,99,100} In recent work, Park found little difference between ASCI-SCF geometries compared to the CASSCF, so we believe that this approximation is justified.^{68} Even though the systems studied by Park are not strongly correlated, our results for Ni(edt)_{2} suggest that even strongly correlated species follow the same trend.

#### 2. Comparing vHCISCF optimized geometries

In the following discussion of optimized geometries, we consider differences in two geometries to be small when they are roughly equal to the tolerance that we use for convergence, *τ*_{drms} (2.4 · 10^{−3} Å), and large when they are more than one order of magnitude larger. The RMSD between initial and final geometries for all vHCISCF optimization is between 0.1 and 0.4 Å, indicating that the structures changed in a substantial manner, and the detailed values are shown in Table III. Since this difference is orders of magnitude larger than our threshold for convergence (3.6 · 10^{−3} Å), the final geometries are all distinct from their starting structures. When comparing optimized geometries starting from different initial coordinates, the structural difference shrink for all singlet states after vHCISCF geometry optimization, while triplet states all diverge from each other. We show a numerical comparison of the geometries in Table IV.

Geometries . | Singlet RMSD . | Triplet RMSD . |
---|---|---|

A | 1.7·10^{−1} | 1.3·10^{−1} |

C | 3.9·10^{−1} | 3.4·10^{−1} |

Geometries . | Singlet RMSD . | Triplet RMSD . |
---|---|---|

A | 1.7·10^{−1} | 1.3·10^{−1} |

C | 3.9·10^{−1} | 3.4·10^{−1} |

Geometries . | Initial RMSD . | Optimized singlet RMSD . | Optimized triplet RMSD . |
---|---|---|---|

A,C | 7.8·10^{−1} | 5.2·10^{−1} | 7.9·10^{−1} |

Geometries . | Initial RMSD . | Optimized singlet RMSD . | Optimized triplet RMSD . |
---|---|---|---|

A,C | 7.8·10^{−1} | 5.2·10^{−1} | 7.9·10^{−1} |

Next, we examine the differences between the experimental and theoretical structures of Fe(PDI). For consistency, we use the same numbering scheme as Stieber *et al.*^{74} Figure 8 shows errors in the initial and final geometries of state (1,**A**). To re-iterate, all initial geometry was optimized with symmetry-broken DFT, while the final geometries were obtained with the vHCISCF. We compare selected bond lengths from both geometries to experimental ones reported in Ref. 74, which were obtained from crystallographic experiments. We only report a subset of the bond lengths because (1,**A**) is only a model system for the experimental complex and they differ slightly in composition. Despite using multireference methods, we found that the optimized geometries did not outperform the M06-L, MN15, and M06-2X. However, the DFT geometry optimizations use def-TZVP while the vHCISCF calculations use cc-pVDZ. The def-TZVP basis is considerably larger than cc-pVDZ and should *a priori* provide a better description of geometry. Due to computational limitations, we were not able to run vHCISCF calculations using triple zeta basis sets. We found that the time to calculate the nuclear gradients increased by nearly 50 fold while the time per MCSCF iteration increased by less than a factor of two, making the gradients the clear bottleneck to using triple zeta basis sets. We note that agreement between the DFT-optimized and experimental geometries may be fortuitous since the Fe(PDI) complex contains more alkyl groups and the experimental evidence suggests that the electronic structure is sensitive to these changes.^{74} However, a further investigation is necessary to test this claim. The other geometry optimizations starting from other geometries and/or spin symmetry showed similar behavior, and we report them in Figs. 10–12 in Appendix E.

#### 3. Singlet-triplet gaps

Reference 74 argued that the ground state is a singlet with a thermally accessible triplet state based on the strong temperature dependence of NMR data for species similar to Fe(PDI). However, the set of compounds that they report in Table II of their paper differ by the alkyl (methyl, ethyl, or isopropyl) groups attached to C(2) and C(8) compared to Fe(PDI). The type of functional group noticeably impacts the extent of the temperature dependence, and it is not clear if those same arguments apply to our model species. In this work, we define thermally accessible as roughly 1 kcal/mol, i.e., “chemical accuracy.” Reference 73 also suggested that the S-T gap of Fe(PDI) small and reported vertical excitation values of 2.2 kcal/mol using RASSCF and 3.4 kcal/mol using RASPT2. Their RAS calculations correlated 22 electrons in 22 orbitals and use a mixed basis where a triple-zeta basis is used for Fe, a double-zeta basis is used for C and N, and a minimal basis is used for H.

Our calculations agree with previous experimental and theoretical work that the ground state has singlet spin symmetry but do not indicate that there is a low-lying triplet state. Our vertical excitation energies at **A** are roughly an order of magnitude larger than those reported by Ortuño and Cramer, while the initial adiabatic gap is more comparable. During optimization, the singlet species relax between 24 and 30 kcal/mol, while the triplet species relax by only 17–20 kcal/mol, leading to larger S-T gaps in contrast to Ortuño and Cramer’s expectations. We report the S-T gaps in Table V from extrapolated SHCI calculations along with DFT and RAS calculations from Ortuño and Cramer^{73} for reference.

Method . | States . | Initial . | Final . |
---|---|---|---|

Ortuño and Cramer^{73} | |||

M06-L | (3,C)–(1,A) | 9.5 | ⋯ |

RASSCF (22e,22o) | (3,A)–(1,A) | 2.2 | ⋯ |

RASPT2 (22e,22o) | (3,A)–(1,A) | 3.4 | ⋯ |

This work | |||

SHCI (40e,40o) | (3,A)–(1,A) | 33.2 | 37.0 |

SHCI (40e,40o) | (3,C)–(1,A) | 8.9 | 15.7 |

Method . | States . | Initial . | Final . |
---|---|---|---|

Ortuño and Cramer^{73} | |||

M06-L | (3,C)–(1,A) | 9.5 | ⋯ |

RASSCF (22e,22o) | (3,A)–(1,A) | 2.2 | ⋯ |

RASPT2 (22e,22o) | (3,A)–(1,A) | 3.4 | ⋯ |

This work | |||

SHCI (40e,40o) | (3,A)–(1,A) | 33.2 | 37.0 |

SHCI (40e,40o) | (3,C)–(1,A) | 8.9 | 15.7 |

#### 4. Multiple local minima

During initial DFT geometries optimizations, we found another triplet geometry optimization with M06-L, which we call **B**. We found that this geometry was more similar to **A** than **C** and when compared to **A** had an RMSD value of 2.8·10^{−2} Å and differed in energy by roughly 1 kcal/mol. After running both singlet and triplet vHCISCF geometry optimizations starting from **B**, we found that the singlet geometry became much closer to **A** (8.0·10^{−3} Å) and had an extrapolated energy that agreed within 1.72·10^{−5} hartree, which is nearly the standard error (1.56·10^{−5}) hartree, leading us to believe that these are the same structure. On the other hand, after geometry optimization for triplet spin symmetry, we found that the geometry differed from the relaxed **C** geometry by 7.9·10^{−1} Å but had an extrapolated energy that differs by less than 1 kcal/mol implying that there are at least two distinct stable triplet structures and that initial geometry choices are important for this species.

## V. CONCLUSION

We have benchmarked nuclear gradients of HCI-type wave functions and demonstrated that like orbitals, nuclear gradients are relatively insensitive to the empirical parameter in HCI. While the nuclear gradients are formally inexact, this error tends to lead to only minimal differences in the final optimized geometries. If more accurate gradients or optimized structures are required, we note that the gradients are systematically improvable if more accurate gradients are required, they are systematically improvable. As a result, accurate molecular structures can be obtained with less computational cost and if accurate energies are desired, one can run a final more accurate SHCI calculation at this geometry.

We studied the challenging catalytic model, Fe(PDI), and used HCI-type wave functions to assess geometries obtained from DFT and to optimize new structures. Given the sensitivity of geometry optimization to the starting coordinates, HCI-based methods can unambiguously compare different geometries and determine the lowest energy structure when DFT cannot. For multireference systems, where single-reference methods, such as DFT and coupled-cluster, will give ambiguous results or even fail, HCI-type wave functions are a useful alternative. Our results agree with previous experimental and theoretical work that the ground state of the Fe(PDI) complex is a singlet. However, we found that the S-T was larger than previously predicted and was *not* thermally accessible. As alluded to earlier, there are several factors that could explain this difference such as the missing functional groups in our model system as well as different correlated methods and basis sets. When compared to the experimentally obtained structure, the vHCISCF performs similarly, albeit slightly worse, than the three DFT functionals studied. Fe(PDI) is a challenging system even for accurate multireference methods, and we encourage others in the community to use it as a benchmark system going forward.

## ACKNOWLEDGMENTS

This work utilized resources from the University of Colorado Boulder Research Computing Group, which was supported by the National Science Foundation (Award Nos. ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University. S.S. acknowledges funding from NSF Award No. CHE1800584 and the Sloan research fellowship. J.E.T.S. gratefully acknowledges support from a fellowship through The Molecular Sciences Software Institute under NSF Grant No. ACI-1547580. J.E.T.S. was also supported by the Flatiron Institute, which is a division of the Simons Foundation. J.L. thanks David Reichman for encouragement and support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**James E. T. Smith**: Investigation (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (lead). **Joonho Lee**: Investigation (supporting); Supervision (supporting). **Sandeep Sharma**: Supervision (lead); Writing – original draft (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are openly available at https://github.com/jamesETsmith/hci_nuc_gradients.

### APPENDIX A: EXTRAPOLATION PROCEDURE

All calculations used for extrapolation use *ϵ*_{2} = 5·10^{−8} hartree and ten different values of *ϵ*_{1} ranging from 5·10^{−5} to 10^{−5} hartree. The procedure for the extrapolation is as follows:

Run all SHCI calculations.

Collect the PT correction (

*E*_{2}) (our independent variable) and total energy (*E*_{SHCI}).Fit linear curves to

*E*_{SHCI}as a function of*E*_{2}.Estimate the error in this extrapolated energy as the standard deviation of the curve fitting parameter representing the y-intercept.

### APPENDIX B: DFT GEOMETRY OPTIMIZATION

Here, we provide the numerical data for Fig. 5 along with ⟨*S*^{2}⟩ and the RMSD compared to the singlet and triplet geometries from Ortuño and Cramer in Table VI.

Functional . | Multiplicity . | Starting geometry . | Energy (hartree) . | Fe spin density . | $\u27e8S\u03022\u27e9$ . | RMSD (Å) . |
---|---|---|---|---|---|---|

M06-2X | 1 | A | −2506.345 333 | 1.949 442 | 1.9345 | 1.1·10^{−1} |

M06-2X | 3 | C | −2506.364 616 | 3.721 854 | 3.9329 | 3.3·10^{−1} |

M06-L | 1 | A | −2506.768 403 | 1.454 747 | 1.2123 | 6.6·10^{−4} |

M06-L | 3 | C | −2506.753 225 | 3.120 475 | 2.9871 | 4.4·10^{−3} |

MN15 | 1 | A | −2505.529 508 | 1.767 367 | 1.6642 | 3.3·10^{−2} |

MN15 | 3 | C | −2505.526 531 | 2.040 751 | 2.5562 | 7.8·10^{−1} |

Functional . | Multiplicity . | Starting geometry . | Energy (hartree) . | Fe spin density . | $\u27e8S\u03022\u27e9$ . | RMSD (Å) . |
---|---|---|---|---|---|---|

M06-2X | 1 | A | −2506.345 333 | 1.949 442 | 1.9345 | 1.1·10^{−1} |

M06-2X | 3 | C | −2506.364 616 | 3.721 854 | 3.9329 | 3.3·10^{−1} |

M06-L | 1 | A | −2506.768 403 | 1.454 747 | 1.2123 | 6.6·10^{−4} |

M06-L | 3 | C | −2506.753 225 | 3.120 475 | 2.9871 | 4.4·10^{−3} |

MN15 | 1 | A | −2505.529 508 | 1.767 367 | 1.6642 | 3.3·10^{−2} |

MN15 | 3 | C | −2505.526 531 | 2.040 751 | 2.5562 | 7.8·10^{−1} |

### APPENDIX C: STABLE UHF SOLUTIONS

Since the UDFT and UHF calculations used to prepare orbitals for the multireference calculations were so sensitive to their optimization settings, we surveyed a large number of these settings to find the lowest energy state possible. We emphasize that such care is needed for this species because MCSCF optimizations are sensitive to the initial orbitals for this species. We examined six different optimization options for the DFT stage and three for the UHF stage leading to a total of 108 ways to prepare the orbitals. In Table VII, we show a subset of these calculations where the final UHF produced a stable wave function with respect to the orbital rotation parameters.

Multiplicity . | Geometry . | DFT opt. strategy . | UHF opt. strategy . | Fe spin density . | ⟨S^{2}⟩
. | Energy (hartree) . |
---|---|---|---|---|---|---|

1 | A | Fast Newton+ADIIS | NEWTON | −2.50 | 3.16 | −2497.225 037 |

1 | A | ADIIS | NEWTON | −2.50 | 3.16 | −2497.225 036 |

1 | A | NEWTON | NEWTON | 2.50 | 3.19 | −2497.224 119 |

1 | A | Fast Newton+ADIIS | ADIIS | −2.01 | 2.69 | −2497.219 975 |

1 | A | ADIIS | ADIIS | −2.01 | 2.70 | −2497.219 970 |

1 | A | Fast Newton+DIIS | NEWTON | 7.50 × 10^{−2} | 2.75 | −2497.188 196 |

1 | A | DIIS | NEWTON | −9.79 × 10^{−2} | 2.75 | −2497.188 062 |

1 | A | Fast Newton+DIIS | ADIIS | −3.90 × 10^{−4} | 2.63 | −2497.187 794 |

1 | A | DIIS | ADIIS | 9.10 × 10^{−4} | 2.61 | −2497.175 954 |

1 | C | ADIIS | NEWTON | −1.98 | 3.60 | −2497.207 645 |

1 | C | NEWTON | NEWTON | −1.78 | 4.16 | −2497.200 086 |

1 | C | Fast Newton+ADIIS | NEWTON | 1.77 | 4.13 | −2497.199 389 |

1 | C | Fast Newton+DIIS | NEWTON | −1.77 | 4.13 | −2497.199 389 |

1 | C | Fast Newton+DIIS | ADIIS | −1.76 | 4.12 | −2497.198 453 |

1 | C | NEWTON | ADIIS | −1.80 | 3.89 | −2497.194 584 |

1 | C | Fast Newton+ADIIS | ADIIS | 1.80 | 3.88 | −2497.194 527 |

1 | C | DIIS | NEWTON | −1.33 × 10^{−2} | 3.54 | −2497.176 719 |

1 | C | DIIS | ADIIS | 7.18 × 10^{−3} | 3.34 | −2497.175 636 |

3 | A | ADIIS | NEWTON | 2.06 | 3.81 | −2497.206 466 |

3 | A | ADIIS | ADIIS | 2.05 | 3.58 | −2497.205 565 |

3 | A | DIIS | NEWTON | 1.25 × 10^{−1} | 4.10 | −2497.204 954 |

3 | A | Fast Newton+DIIS | NEWTON | 1.88 | 4.61 | −2497.196 780 |

3 | A | NEWTON | NEWTON | 1.88 | 4.61 | −2497.196 780 |

3 | C | Fast Newton+DIIS | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | NEWTON | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | Fast Newton+ADIIS | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | ADIIS | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | DIIS | NEWTON | 1.84 | 5.11 | −2497.200 823 |

3 | C | DIIS | ADIIS | 1.99 | 4.33 | −2497.196 357 |

Multiplicity . | Geometry . | DFT opt. strategy . | UHF opt. strategy . | Fe spin density . | ⟨S^{2}⟩
. | Energy (hartree) . |
---|---|---|---|---|---|---|

1 | A | Fast Newton+ADIIS | NEWTON | −2.50 | 3.16 | −2497.225 037 |

1 | A | ADIIS | NEWTON | −2.50 | 3.16 | −2497.225 036 |

1 | A | NEWTON | NEWTON | 2.50 | 3.19 | −2497.224 119 |

1 | A | Fast Newton+ADIIS | ADIIS | −2.01 | 2.69 | −2497.219 975 |

1 | A | ADIIS | ADIIS | −2.01 | 2.70 | −2497.219 970 |

1 | A | Fast Newton+DIIS | NEWTON | 7.50 × 10^{−2} | 2.75 | −2497.188 196 |

1 | A | DIIS | NEWTON | −9.79 × 10^{−2} | 2.75 | −2497.188 062 |

1 | A | Fast Newton+DIIS | ADIIS | −3.90 × 10^{−4} | 2.63 | −2497.187 794 |

1 | A | DIIS | ADIIS | 9.10 × 10^{−4} | 2.61 | −2497.175 954 |

1 | C | ADIIS | NEWTON | −1.98 | 3.60 | −2497.207 645 |

1 | C | NEWTON | NEWTON | −1.78 | 4.16 | −2497.200 086 |

1 | C | Fast Newton+ADIIS | NEWTON | 1.77 | 4.13 | −2497.199 389 |

1 | C | Fast Newton+DIIS | NEWTON | −1.77 | 4.13 | −2497.199 389 |

1 | C | Fast Newton+DIIS | ADIIS | −1.76 | 4.12 | −2497.198 453 |

1 | C | NEWTON | ADIIS | −1.80 | 3.89 | −2497.194 584 |

1 | C | Fast Newton+ADIIS | ADIIS | 1.80 | 3.88 | −2497.194 527 |

1 | C | DIIS | NEWTON | −1.33 × 10^{−2} | 3.54 | −2497.176 719 |

1 | C | DIIS | ADIIS | 7.18 × 10^{−3} | 3.34 | −2497.175 636 |

3 | A | ADIIS | NEWTON | 2.06 | 3.81 | −2497.206 466 |

3 | A | ADIIS | ADIIS | 2.05 | 3.58 | −2497.205 565 |

3 | A | DIIS | NEWTON | 1.25 × 10^{−1} | 4.10 | −2497.204 954 |

3 | A | Fast Newton+DIIS | NEWTON | 1.88 | 4.61 | −2497.196 780 |

3 | A | NEWTON | NEWTON | 1.88 | 4.61 | −2497.196 780 |

3 | C | Fast Newton+DIIS | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | NEWTON | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | Fast Newton+ADIIS | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | ADIIS | NEWTON | 3.76 | 5.42 | −2497.264 950 |

3 | C | DIIS | NEWTON | 1.84 | 5.11 | −2497.200 823 |

3 | C | DIIS | ADIIS | 1.99 | 4.33 | −2497.196 357 |

### APPENDIX D: MCSCF SENSITIVITY TO *ϵ*_{1}

We investigated the sensitivity of the vHCISCF orbitals to the HCI parameter, *ϵ*_{1}, analogous to the analysis from Smith *et al.*^{56} This process involves two steps: first, converging the vHCISCF orbitals with respect to several different values of *ϵ*_{1} and second, running an SHCI calculation using these orbitals with identical settings. Despite the fact that Fe(PDI) is more multireference than the Fe(porphyrin) system studies in Ref. 56, we find that the orbitals are relatively insensitive to *ϵ*_{1} and show the results in Fig. 9.

### APPENDIX E: GEOMETRY COMPARISONS

We show the error between several theoretical results and the experimentally determined bond lengths for a subset of bond lengths near the central Fe atom (Figs. 10–12). For all states, we compare vHCISCF to M06-L and for (1,**A**), (3,**A**), and (3,**C**) and we compare to MN15 and M06-2X as well, since those calculations were already performed during our earlier analysis. In general, vHCISCF performs similarly to the ensemble of DFT results but does not typically provide the most accurate bond lengths.

## REFERENCES

_{4}CaO

_{5}cluster in photosystem II