Near-exact nuclear gradients of complete active space self-consistent field wave functions

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).⁴ 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.⁵⁷ 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 ϵ₁. 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.⁵⁶ 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.⁶⁸ 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.⁷³ 

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 ϵ₁. 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₂, 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.

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 |Ψ⟩ = ∑_ic_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.

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=⟨Da|Ĥ|Di⟩$⁠, c_i is the amplitude of the |D_i⟩, and ϵ₁ is a user-specified parameter that controls the size of the variational space, i.e., the accuracy of the variational wave function. A smaller ϵ₁ means that $V$ is closer to the FCI determinant space and in the limit that ϵ₁ 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, $∑iHkiciE0−Hkk>ϵCIPSI$⁠, 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.⁷⁷ 

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 ϵ₂ becomes smaller, and in the limit ϵ₂ → 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 ϵ₂ 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₂(ϵ₂) 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₂(ϵ₂) 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.⁷⁸ 

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.⁷⁹ 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₁₂ 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〉=∑μAOCμi|μ〉$⁠. 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 $∂Eelec∂ci=0$ the term always vanishes. Similarly, for wave functions where the molecular orbital parameters are variational, $∂Eelec∂Cμ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 $γij=∑σ⟨Ψ|aiσ†ajσ|Ψ⟩$ is the one-body reduced density matrix (RDM) and $Γijkl=12∑σσ′⟨Ψ|aiσ†ajσ′†akσ′alσ|Ψ⟩$ is the two-body RDM. Taking derivative of Eq. (6) with respect to a general nuclear coordinate a generates terms, including $∂γij∂a$ and $∂Γijkl∂a$⁠. 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.⁸⁰ 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.⁷⁹ 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 ϵ₁. 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 ϵ₁ is greater than 10 mhartree and the wave functions consisted of less than 100 determinants. Since typical HCI calculations use much smaller values of ϵ₁ and many more determinants, we do not explore these discontinuities further in this work.

For all DFT geometry optimizations, we use Gaussian16 with the def2-TZVP basis set and density fitting to reduce the computational cost.⁸¹ 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.⁸² For HCISCF calculations we use our own implementation of the HCI algorithm, called Dice,⁵⁶ in tandem with PySCF. The nuclear gradients are calculated in PySCF using the RDMs obtained from Dice and geomeTRIC⁸³ 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⁻⁶ hartree for the energy (τ_E), 6  ·  10⁻⁴ hartree/bohr for the root mean squared (rms) of the nuclear gradient (τ_grms), 9  ·  10⁻⁴ hartree/bohr for the maximum of the nuclear gradient (τ_gmax), 2.4  ·  10⁻³ Å for the rms of the atomic displacement (τ_drms), and 3.6  ·  10⁻³ Å 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.⁸⁴ We outline the procedure for selecting the active space orbitals below:

1. Run unrestricted DFT calculation using the PBE0 functional^85,86 to generate an initial guess for Step 2.

2. Run an unrestricted Hartree–Fock (UHF) calculation to generate spin orbitals without correlation from an exchange-correlation functional.

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

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

5. Use the UNOs in a “loose” HCI calculation where the ϵ₁ = 10⁻⁴ 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.

6. Use the 1-RDM from the previous step to calculate new NOs and sort based on NOONs.

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

We emphasize that Steps 1 and 2 are sensitive to the optimization strategy used for the SCF procedure, e.g., DIIS,⁸⁹ ADIIS,⁹⁰ 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⁹¹ to align geometries and report the difference in geometry as $‖Δ‖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.⁹² 

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.

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.⁸⁴ 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₂, Sc₂, trans-stilbene, and bis-(ethylene-1,2-dithiolato)nickel [Ni(edt)₂] 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₂ and Sc₂, these active space sizes correspond to the full valence space. These active spaces correspond to the full valence space for N₂ and Sc₂, the conjugated 2p_z orbitals on carbon for trans-stilbene, and the conjugated 2p_z orbitals of carbon and sulfur as well as the 3d orbitals of the nickel atom for Ni(edt)₂. Schlimgen and Mazziotti have shown that Ni(edt)₂ 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₂ 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⁹⁴ 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 ϵ₁ value, which increases the number of determinants in the HCI wave function. Using ϵ₁, we can control the accuracy of our gradients, compared to the CASSCF, and we find that we can use a relatively large value of ϵ₁, 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 ϵ₁, 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⁸¹ 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 ϵ₁ = 5  ·  10⁻⁴ 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₂ and Ni(edt)₂. For the typical range of ϵ₁ values used in research applications (ϵ₁ = 10⁻⁴–10⁻⁵ 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.

We use two geometries reported by Ortuño and Cramer, which they obtained using DFT: the first, “¹BS(1,1),” which we call A, optimized for the singlet spin symmetry, and the second, “³BS(3,1),” which we call C, optimized for the triplet.⁷³ 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.”⁷³ 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.

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 ϵ₁ = 7.5·10⁻⁵ hartree. Given the sensitivity of the MCSCF solutions to the input orbitals, we tested several values of ϵ₁ 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.⁵⁶ We report this analysis in  Appendix D. To ensure that our choice of ϵ₁ was appropriate, we tested the use of a more accurate ϵ₁ value of 5 · 10⁻⁵ hartree and found that it produced geometries that are quite similar to those produced when using ϵ₁ = 7.5·10⁻⁵. The RMSD value between the two optimized geometries is 5.6·10⁻³ Å, which is less than the loose convergence criteria used by Gaussian (6.7·10⁻³ Å). vHCISCF wave functions using ϵ₁ = 7.5·10⁻⁵ contained roughly 3.5·10⁶ determinants and those using ϵ₁ = 5·10⁻⁵ contained roughly 7.7·10⁶. 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 ϵ₁ 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 ϵ₁ to extrapolate to the FCI limit, the tightest one being 10⁻⁵, and ϵ₂ = 5·10⁻⁸ 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.

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.⁶⁸ Even though the systems studied by Park are not strongly correlated, our results for Ni(edt)₂ suggest that even strongly correlated species follow the same trend.

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⁻³ Å), 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⁻³ Å), 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.

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.⁷⁴ 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.⁷⁴ 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.

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⁷³ for reference.

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⁻² Å 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⁻³ Å) and had an extrapolated energy that agreed within 1.72·10⁻⁵ hartree, which is nearly the standard error (1.56·10⁻⁵) 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⁻¹ Å 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.

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.

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.

The authors have no conflicts to disclose.

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

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

All calculations used for extrapolation use ϵ₂ = 5·10⁻⁸ hartree and ten different values of ϵ₁ ranging from 5·10⁻⁵ to 10⁻⁵ hartree. The procedure for the extrapolation is as follows:

1. Run all SHCI calculations.

2. Collect the PT correction (E₂) (our independent variable) and total energy (E_SHCI).

3. Fit linear curves to E_SHCI as a function of E₂.

4. Estimate the error in this extrapolated energy as the standard deviation of the curve fitting parameter representing the y-intercept.

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

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.

We investigated the sensitivity of the vHCISCF orbitals to the HCI parameter, ϵ₁, analogous to the analysis from Smith et al.⁵⁶ This process involves two steps: first, converging the vHCISCF orbitals with respect to several different values of ϵ₁ 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 ϵ₁ and show the results in Fig. 9.

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.