Bridging the Homogeneous-Heterogeneous Divide: Modeling Spin for Reactivity in Single Atom Catalysis

Single atom catalysts (SACs) are emergent catalytic materials that have the promise of merging the scalability of heterogeneous catalysts with the high activity and atom economy of homogeneous catalysts. Computational, first-principles modeling can provide essential insight into SAC mechanism and active site configuration, where the sub-nm-scale environment can challenge even the highest-resolution experimental spectroscopic techniques. Nevertheless, the very properties that make SACs attractive in catalysis, such as localized d electrons of the isolated transition metal center, make them challenging to study with conventional computational modeling using density functional theory (DFT). For example, Fe/N-doped graphitic SACs have exhibited spin-state dependent reactivity that remains poorly understood. However, spin-state ordering in DFT is very sensitive to the nature of the functional approximation chosen. In this work, we develop accurate benchmarks from correlated wavefunction theory (WFT) for relevant octahedral complexes. We use those benchmarks to evaluate optimal DFT functional choice for predicting spin state ordering in small octahedral complexes as well as models of pyridinic and pyrrolic nitrogen environments expected in larger SACs. Using these guidelines, we determine Fe/N-doped graphene SAC model properties and reactivity as well as their sensitivities to DFT functional choice. Finally, we conclude with broad recommendations for computational modeling of open-shell transition metal single-atom catalysts.


INTRODUCTION
Single atom catalysts (SACs) (Yang et al., 2013) are emergent catalytic materials (Yang et al., 2013;Liang et al., 2015Liang et al., , 2017 that have the promise of merging the scalability of heterogeneous catalysts with the high activity and atom economy of homogeneous catalysts, but the reactivity of SACs is poorly understood (Figure 1). Fe/N-doped graphene SACs have been demonstrated for critical transformations such as selective hydrocarbon oxidation (Liu et al., 2017), including ambient methane to methanol conversion (Cui et al., 2018), and as non-Pt oxygen reduction reaction (ORR) electrocatalysts (Li et al., 2016;Yang et al., 2018). The short-lived, variable FIGURE 1 | Iron catalysts in three classes: homogeneous tetraphenyl porphyrin (Top), N-doped graphene single atom catalysts (SACs, Middle), and heterogeneous hematite (Bottom). In all cases, carbon is shown in gray, nitrogen in blue, iron in purple, and oxygen in red.
SAC active sites that are fundamentally sub-nm-scale challenge the resolution of spectroscopic techniques (Fei et al., 2015;Wang and Zhang, 2016), making first-principles modeling essential to mechanistic study.
Although perhaps surprising in the context of heterogeneous catalysis, strong spin-state dependence in reactivity is wellknown in homogeneous catalyst (Schröder et al., 2000) analogs. Thus, it follows that paradigms that work in coordination chemistry might apply to SACs as well. The near-octahedral coordination environment around a metal center produces distinct quantum mechanical spin states (i.e., local metal magnetic moments) that are determined by the ligand-field strength as well as oxidation state and metal identity (Tsuchida, 1938). Different spin states often have distinct reaction barriers (Schröder et al., 2000;Schwarz, 2017) in a paradigm known as two-state reactivity (TSR) (Shaik et al., 1995;Schröder et al., 2000;Schwarz, 2017). TSR was first identified for Fe + ions, where oxidation of H 2 or CH 4 to H 2 O or CH 3 OH, respectively, is limited by spin inversion from a steep ground state, high spin (HS) surface to a more reactive but excited state low spin (LS) surface. For cases such as iron-oxo porphyrin systems that have nearly degenerate spin states, different pathways can indeed lead to distinct products (Kamachi and Yoshizawa, 2003;Ji et al., 2015).
In minimal model single-site catalysts, we recently demonstrated (Gani and Kulik, 2018) that bond elongation or compression has an effect similar to modulating ligand field strength, which could also alter ground state spin and reactivity in an interconnected manner. Most Fe(II)/N complexes have near degenerate HS and LS states because nitrogen ligands are of intermediate field strength, but small changes in the N-N separation of bidentate ligands that make up the octahedral complex are known to tune the experimental ground state of the material (Phan et al., 2017). In metal-doped graphene, strain has been predicted to change the ground state spin (Huang et al., 2011). Even changes in orientation of ligands (i.e., one equatorial swapped with one axial) have been experimentally observed to change the favored ground state spin of related molecular complexes (Andris et al., 2016).
The confined nature of metal d states and interactions with localized p orbitals from organic ligand atoms impart properties to open-shell SACs in a manner more closely resembling the chemical bonding of homogeneous catalysts than bulk metal counterparts. These very features, i.e., quantum size effects (Yang et al., 2013) at an open shell, high-valent metal atom, that make SACs reactive for essential catalytic transformations (Qiao et al., 2011;Yang et al., 2013;Zitolo et al., 2015Zitolo et al., , 2017Zhang et al., 2018) also make conventional computational tools used in heterogeneous catalysis (i.e., plane wave, semi-local density functional theory, or DFT) ill-suited to predictive SAC study. Well-localized electrons are disproportionately affected by selfinteraction error in approximate DFT (Cohen et al., 2011;Kulik, 2015), leading to an imbalanced treatment of differing spin and oxidation states (Ganzenmüller et al., 2005;Kulik et al., 2006;Droghetti et al., 2012;Kulik, 2015, 2017;Mortensen and Kepp, 2015;Gani and Kulik, 2017).
Despite evidence of the importance of spin in homogeneous (Abram et al., 2014;Zhu et al., 2016;Schwarz, 2017) and SAC catalysts (Liu et al., 2017), most first-principles studies of SACs (Chu et al., 2015;Ma et al., 2016;Xu et al., 2018) have avoided directly quantifying the role of metal center spin in reactivity, with few exceptions (Impeng et al., 2014;Fong et al., 2018;Sirijaraensre and Limtrakul, 2018). In most studies, the magnetic moment is calculated with a semi-local DFT functional known to produce erroneous magnetic moments (Kulik, 2015;Ioannidis and Kulik, 2017;Wilbraham et al., 2017) and the magnetization is often allowed to vary along the reaction coordinate . However, in confined metal centers, spin states are well defined and transitions between spin states occur with low probability because they are quantum mechanically forbidden. Spin state transitions can become kinetically limiting (Shaik et al., 1995;Schröder et al., 2000;Schwarz, 2017), explaining unexpected experimental reactivity (Andris et al., 2016). Within the homogeneous catalysis community (Harvey, 2014;Hernández-Ortega et al., 2015), significant effort has been made to develop tools to assess whether spin crossover is kinetically limiting, but this is not the case in SACs.
Thus, in this work, we carry out highly accurate correlated wavefunction theory calculations to develop benchmarks for transition metal complex spin state ordering with ligands that model the environment observed in single atom catalysts. Using these benchmarks, we identify trends in DFT functional performance, and then we evaluate how these observations influence predictions of the stability, reactivity, and ground state identity in models of Fe(II)/N-doped graphene SACs. Finally we provide our conclusions and outlook, including recommendations for computational modeling in this emergent space of single atom catalysis.

Octahedral Transition Metal Complexes
Initial structures of octahedral transition metal (TM) complexes with H 2 O, NH 3 , pyridine, and pyrrole ligands were built with the molSimplify toolkit (Ioannidis et al., 2016) with both ligand force-field pre-optimization and trained metal-ligand bond length features enabled. For the hexa-aqua and hexaammine complexes, M(II) Ti-Ni and M(III) V-Cu ions were studied, but the pyrrole and pyridine ligands were only studied in complex with Fe(II) or Fe(III). The formal charges assigned to the ligands were neutral in all cases except for pyrrole, which was deprotonated and given a −1 charge. High-spin (HS)-low-spin (LS) states studied in this work were defined as: triplet-singlet for d 2 Ti(II)/V(III) and d 8 Ni(II)/Cu(III), quartetdoublet for d 3 V(II)/Cr(III) and d 7 Co(II)/Ni(III), quintet-singlet for d 4 Cr(II)/Mn(III) and d 6 Fe(II)/Co(III), and sextet-doublet for d 5 Mn(II)/Fe(III). Intermediate-spin (IS) states were also studied: triplet d 4 Cr(II)/Mn(III) or d 6 Fe(II)/Co(III) and quartet d 5 Mn(II)/Fe(III).

Fe/N-Doped Graphene SAC Finite Models
Two possible Fe(II) coordination environments in finite graphitic SAC models were investigated with DFT, and in both the metal is coordinated by four nitrogen atoms substituted in the graphene structure. In both cases, we employed a hydrogen-atomterminated graphene flake to avoid increasing computational cost in accordance with prior SAC computational studies that used finite models . First, a FeN 4 C 10 compound (chemical formula: C 36 N 4 H 16 Fe) was studied in which all coordinating nitrogen atoms were in six-membered rings (i.e., pyridinic N). This structure is formed by removing two adjacent C atoms from C 42 H 16 and replacing the four C atoms surrounding the vacancy with N atoms. This active site would correspond to two adjacent point defects in graphene, as has been observed experimentally (Banhart et al., 2011). A second compound, FeN 4 C 12 (chemical formula: C 40 N 4 H 16 Fe), was also studied in which all coordinating nitrogen atoms were in fivemembered rings (i.e., pyrrolic N). This structure was formed by removing two C atoms from a C 46 H 16 structure, which contains two seven-membered rings in the center surrounded by four five-membered rings. Thus, this structure would require vacancy migration that has also been experimentally observed (Banhart et al., 2011). The two C atoms were removed from where the seven membered rings are joined, and the four inward-facing C atoms that are part of the five-membered rings were replaced with N atoms. All initial coordinates were generated by drawing the 2D structures with ChemDraw and converting the xml structures to 3D coordinates with the molSimplify (Ioannidis et al., 2016) interface to OpenBabel (O'Boyle et al., 2011) followed by force field optimization with the universal force field (Rappé et al., 1992). Singlet, triplet, and quintet spin states were studied, and all simulations had zero net charge.

Fe/N-Doped Graphene SAC Periodic Models
Periodic analogs to the flake models were studied starting from a 4 × 4 supercell of graphene at its experimental lattice parameter (Trucano and Chen, 1975). A smaller supercell than suggested (i.e., 7 × 7) in previous work (Krasheninnikov et al., 2009) was used for computational efficiency, and future work should focus on the effect of supercell size on dopant properties. The pyridinic (chemical formula: FeN 4 C 18 ) SAC model was created following the same vacancy/N-atom replacement approach as in the finite case. For the pyrrolic (chemical formula: FeN 4 C 20 ) SAC model, we started from the pyridinic case, inserting C atoms into the five-membered FeN 2 C 2 ring. Next, we adjusted the adjacent sixmembered rings into five-membered rings to create pyrrolic N atoms. In this small supercell, an eight-membered C ring was then formed next to the five membered rings. Neutral systems were studied by spin polarized, fixed magnetization periodic calculations in singlet, triplet, and quintet states.

Transition Metal Complexes
All LS, IS, and HS complexes were geometry optimized with DFT using the PBE0 (default 25% exchange) global hybrid GGA functional (Adamo and Barone, 1999) with the def2-TZVP basis set (Weigend and Ahlrichs, 2005) in ORCA v.4.0 (Neese, 2018). Singlet states were calculated in a restricted formalism, whereas all remaining calculations were open shell and required level shifting (Saunders and Hillier, 1973) in select cases to aid selfconsistent field convergence typically with a value of 1.0 eV but as large as 100.0 eV in one case ([Mn(NH 3 ) 6 ] 2+ ). The optimizations were carried out using the BFGS algorithm in redundant internal coordinates implemented to the default tolerances of 3 × 10 −4 hartree/bohr for the maximum gradient and 5 × 10 −6 hartree for the change in self-consistent field (SCF) energy between steps. All calculations at other levels of theory or with differing functional definitions were obtained as single point energies on these optimized geometries. The effect of Hartree-Fock (HF) exchange fraction choice on spin-state energetics within DFT was investigated by altering the fraction in a modified form of the PBE0 global hybrid. The HF exchange fraction was varied from as low as 0% [i.e., a pure PBE GGA (Perdew et al., 1996)] to as high as 100% HF exchange in increments of 10-20%, as indicated in the text, again using the def2-TZVP basis set. In previous work , we found tuning range-separation parameters in range-corrected hybrids to have a comparable effect on density and energetics of transition metal complexes to global exchange tuning, and therefore we focus on only global exchange tuning in this work.

Fe/N-Doped Graphene Flake Models
Geometry optimizations and single-point energy calculations were performed with ORCA v4.0. All DFT methodology was kept the same as for transition metal complexes, including geometry optimizing at PBE0 (25% exchange) and carrying out single points at modified exchange fractions in PBE0 in 10% increments in conjunction with the def2-TZVP basis set, except as noted below. All singlet, triplet, and quintet calculations were carried out in an unrestricted formalism. All calculations employed the resolution of the identity (RI) (Baerends et al., 1973;Whitten, 1973;Dunlap et al., 1979;Eichkorn et al., 1995Eichkorn et al., , 1997Kendall and Fruchtl, 1997) and the chain-ofsphere (COSX) ) approximations with the auxiliary basis set def2/J (Weigend, 2006) and def2-TZVP/C (Hellweg et al., 2007) for all atoms to accelerate the calculations while introducing marginal errors (Kossmann and Neese, 2009). Molecular structures and orbitals were visualized and plotted with VESTA (Momma and Izumi, 2011).

Periodic DFT Calculations
All systems were calculated with both the PBE (Perdew et al., 1996) semi-local GGA functional and the HSE06 (Heyd et al., 2003(Heyd et al., , 2006 local, range-separated GGA hybrid using the plane wave, periodic boundary condition Quantum-ESPRESSO (Paolo et al., 2009) code. Norm-conserving pseudopotentials for C, N, and Fe were generated with OPIUM (Rappe et al., 1990). The wavefunction and charge density cutoffs employed were 50 Ry and 200 Ry, respectively. A Monkhorst-Pack k-point grid of 8 × 8 × 1 was used for efficiency after confirming convergence of total energies with k-point mesh size. Total energies were converged to 0.1 meV and forces were converged to 1 meV/Å. Quantum-ESPRESSO post-processing tools were employed to visualize the spin density and projected density of states. Variable cell relaxation was employed to obtain final lattice parameters with PBE GGA for the pyridinic (result: 8.34 Å × 7.55 Å) and pyrrolic (result: 8.35 Å) models. A vacuum between each SAC layer of 10 Å was included along with a dipole correction to limit periodic image effects (Bengtsson, 1999). The HSE06 calculations were obtained as single point energies applied to these structures.

Correlated WFT
Complete active-space second-order perturbation theory (CASPT2) (Andersson et al., 1992) calculations were performed with OpenMolcas 1 (Aquilante et al., 2016) on M(II)/M(III) hexa-aqua and hexa-ammine octahedral complexes. Calculations were carried out with two active space definitions: the standard active space and an extended active space. For the standard active space, we followed literature recommendations for TM complexes (Pierloot, 2003;Veryazov et al., 2011) to include five orbitals with TM 3d character, two σ bonding orbitals describing covalent metal-ligand bonding, and five double-shell d orbitals for mid-row and later transition metals [i.e., Mn(II/III) and later]. For n 3d electrons, the standard active space is (n + 4, 7) (i.e., for Sc-Cr) or (n + 4, 12) (i.e., for Mn-Cu). In the extended active space, we followed additional literature recommendations (Wilbraham et al., 2017) to include the metal 3s orbital and an unoccupied counterpart, giving an active space of (n + 6, 9) (i.e., for Sc-Cr) or (n + 6, 14) (i.e., for Mn-Ni). Relativistic atomic natural orbital (ANO-rcc) basis sets (Roos et al., , 2005 contracted to [7s6p5d3f2g1h] for the metal center, [4s3p2d1f] for O and N, and [3s1p] for H were used together with the scalar relativistic Douglas-Kroll Hamiltonian (Douglas and Kroll, 1974;Hess, 1986). The 10 core orbitals were frozen in all calculations. An imaginary level shift (Forsberg and Malmqvist, 1997) of 0.1 was used, and a zeroth-order Hamiltonian empirical correction, i.e., the IPEA shift (Ghigo et al., 2004), was set as 0.5 a.u. or varied as described in the main text to identify the effect on spin state energetics. For difficult to converge complete active space self-consistent field iterations, a level shift was applied to the Hamiltonian with shift value 1.0.

Correlated WFT Results
We first conducted correlated wavefunction theory (WFT) calculations to generate reference spin-splitting energies of model first row octahedral transition metal complexes. Our focus is on weak field hexa-aqua and hexa-ammine complexes that are small (i.e., 19-25 atoms in size) and of comparable ligand field strength to the coordination environment in SACs. Although CASPT2 is often the method of choice for predicting spin state energetics for molecules that are either too large or multireference in character to be comfortably treated with CCSD(T) (Pierloot et al., 2017), a number of calculation parameters can strongly influence the CASPT2 predictions. Specifically, spin state energetics can be influenced by the active space choice and the formulation of the zeroth-order Hamiltonian, i.e., the value of the IPEA shift (Ghigo et al., 2004). Here, we investigate the effects of both of these factors and then select reference results for DFT calculations.
Several studies (Kepenekian et al., 2009;Lawson Daku et al., 2012;Vela et al., 2016;Pierloot et al., 2017) have shown that the standard IPEA shift of 0.25 a.u. in CASPT2 overstabilizes high spin states. However, there is no universal agreement about the best solution to this problem. Some (Kepenekian et al., 2009;Lawson Daku et al., 2012;Vela et al., 2016) have recommended increasing the IPEA shift to 0.5-0.7 a.u. based on comparison with experimental or MRCI results, whereas others (Pierloot et al., 2017) recommend the standard IPEA value because increased IPEA can reduce the high spin bias but only at the expense of deteriorating the CASPT2 description of valence correlation. To understand the effect of IPEA shift, we obtained CASPT2 spin-splitting energies with IPEA shifts of 0.0, 0.50, and 1.5 a.u. and standard active spaces (Figure 2 and see Computational Details). For all complexes, increasing IPEA shifts the high-spin/low-spin splitting toward more positive values, reducing high spin stabilization. The range of energetics calculated with different IPEA values for each complex provides a measure of IPEA sensitivity of each transition metal complex's spin-splitting energy (Figure 2). The IPEA sensitivity is consistently largest for cases where the high spin state has four more unpaired electrons than the low spin state (i.e., Cr 2+ through Co 3+ ).
To focus only on spin state definitions in which the high spin state has two more unpaired electrons, we compare the high spin/low spin splitting, E H−L , for early and late TMs with the high spin/intermediate spin splitting, E H−I , for midrow TMs. For these two electron differences, we observe similar sensitivities, with a ca. 10-20 kcal/mol positive shift when the IPEA is changed from 0.0 to 1.5 a.u. (Figure 2). This observation excludes only very early or late TM complexes that have smaller sensitivities. In comparison, the 4 electron-difference cases have higher sensitivities of around 30 kcal/mol or more (Figure 2). Overall, most early-and mid-row complexes remain high spin regardless of the IPEA shift (i.e., below the zero axis in Figure 2). The smaller energetic differences between states in later TM complexes [e.g., hexa-aqua Ni(III)] mean that the IPEA shift can change the ground state from high spin to low spin for very large IPEA shifts. For subsequent calculations, we selected the 0.5 a.u. IPEA shift but note that variation across the commonly employed range of 0.0-0.5 a.u. can shift the predicted E H−L ( E H−I ) spin-splitting for mid-row [e.g., Fe(II/III)] complexes by around 10 (3) kcal/mol (Figure 2). We next investigated the effect of use of a standard active space vs. a more extended active space in the CASPT2 calculations of E H−L and E H−I using the 0.5 a.u. IPEA shift (Supplementary Table 1). For most complexes, the energetic difference due to active space change is on the order of a few kcal/mol, which suggests that the calculation is converged with respect to active space size, as motivated in previous work (Wilbraham et al., 2017). The major outlier identified is [Mn(NH 3 ) 6 ] 2+ , which is strongly low-spin in the standard active space but becomes high-spin like hexa-aqua Mn(II) with the extended active space, exhibiting greater active space dependence than had been observed in Mn(II) porphyrins (Yang et al., 2016) (Supplementary Table 1). This discrepancy is likely caused by orbital rotation of some active orbitals into the inactive metal 3s/4s orbitals, as suggested in recent work (Radoń and Drabik, 2018) on aqua complexes. After removing this outlier, the mean absolute difference between the standard and extended active space results for all E H−L and E H−I combinations is 3.3 kcal/mol for the hexa-aqua and 5.9 kcal/mol for the hexa-ammine complexes. The mean signed error is near zero for the hexaaqua cases, and weakly negative (ca. −3 kcal/mol) for the hexaammines (Supplementary Table 1). Generally, discrepancies are smallest for the hexa-aqua complexes throughout and especially small for the early or late TMs (e.g., Ti 2+ -V 2+ and Ni 3+ -Cu 3+ ), typically as little as 0-2 kcal/mol (Supplementary Table 1). In mid-row cases, there is no universal preference as to whether E H−L or E H−I has more active space dependence. Given the small size of the studied octahedral complexes, we selected the extended active space calculations as reference values for comparison to DFT. Use of the larger active space changed some ground state spin assignments. When calculated with the extended active space, almost all hexa-aqua complexes are high spin, excluding only weakly low spin [Ni(H 2 O) 6 ] 3+ with a E H−L of 0.49 kcal/mol. The slightly stronger ligand field in the hexa-ammine complexes produces some additional LS late-TM complexes [e.g., Co(III) and Cu(III)] along with the analogous Ni(III) complex. Examining isoelectronic metals generally reveals that the later, more oxidized metal has only a weak high-spin-stabilizing effect that is smaller than the ligandfield effect in cases where the two metals converge to similar electronic states (Figure 2 and Supplementary Table 1).

DFT Functional Performance
Despite the high accuracy of CASPT2 for treating spin state energetics in TM complexes, the high computational cost and sensitivity to active space definition and parameters limit its application on SACs. We thus sought to identify the extent to which DFT functionals can be selected or tuned to reproduce the spin-splitting energetics obtained with CASPT2. We focused on the exchange fraction within the global hybrid PBE0 (Adamo and Barone, 1999), motivated by previous observations of comparable behavior in tuning range-separated hybrids , global hybrids with a different correlation functional (Ioannidis and Kulik, 2015), or those that incorporate meta-GGA exchange .
From all complexes studied with CASPT2, we narrowed our focus to those containing nominally 3-7 3d electrons (i.e., V 2+ to Ni 3+ ) that are most likely to be good candidates for understanding spin-state dependent single atom catalysis. We determined the effect of varying PBE0 exchange fraction on both E H−L and E H−I for the relevant subset (i.e., 4-6 3d electrons). In accordance with prior work (Droghetti et al., 2012;Ioannidis and Kulik, 2015), we anticipated the sensitivity of these quantities to exchange fraction and the optimal exchange fraction to minimize error with respect to CASPT2 to be chemistry dependent. However, we aimed to identify if these trends can be readily rationalized (Gani and Kulik, 2017) or learned  for use in SAC modeling.
As expected (Droghetti et al., 2012;Ioannidis and Kulik, 2015;Zhao and Kulik, 2018), spin-state energetics vary linearly with exchange fraction over a wide range (i.e., 0-50%) and high-spin states are stabilized with increasing exchange fraction for all transition metal complexes studied (Supplementary Tables 2-3). We quantified the exchange sensitivity of the spin-splitting energetics (Ioannidis and Kulik, 2015) by an approximate linear fit: 3d filling (i.e., varying only 2-4 kcal/(mol·HFX), see Figure 3). Conversely, hexa-ammine Mn(II), Fe(II), and Co(II) complexes have increased spin-splitting exchange sensitivity over the earlier TM complexes (Figure 3). For both ligand fields, the M(III) complexes are even more varied, with the least exchange sensitivity being observed in either late or early transition metal complexes (Figure 3). These observations are consistent with the fact that these complexes should have the least difference in electron delocalization between the two spin states, reducing exchange sensitivity (Gani and Kulik, 2017). Although M(III) complexes are more variable, isoelectronic +2/+3 complexes do have somewhat comparable exchange sensitivity (Figure 3).
In comparing E H−L for all complexes, exchange sensitivity is universally higher for these mid-row cases due to the enhanced sensitivity of the four-electron difference energetics (e.g., [M(H 2 O) 6 ] 2+ : −57 kcal/(mol·HFX) for Mn vs. −25 kcal/(mol·HFX) for V, see Supplementary Figure 1). The remaining question is whether the exchange fraction can be appropriately tuned within a modified form of PBE0 to obtain E H−L and E H−I values that match CASPT2 results. Comparing the range of spin-splitting energies obtained from 0 to 100% exchange to the CASPT2 extended active space values generally reveals that high exchange fractions (c.a. 40%) are required to reproduce CASPT2 results (Figure 4 and E H−L -only results shown in Supplementary Figure 2). With the exception of Ni(III) or V(II)/Cr(III) hexa-ammines, pure PBE GGA E values are much more positive (i.e., low-spin biased) than the WFT results (Figure 4). Increasing exchange thus in most cases improves agreement with WFT, but the optimal exchange fraction for reproducing the CASPT2 result varies significantly for the different complexes (Figure 4). For five of the cases [e.g., midrow Mn(II)(H 2 O) 6 , Fe(II)(NH 3 ) 6 , and Co(II)(H 2 O) 6 ], the optimal exchange fraction is larger than  (Figure 4). Only three complexes [i.e., V(II)(H 2 O) 6 , Cr(III)(H 2 O) 6 , and Ni(III)(H 2 O) 6 ] have an optimal exchange fraction corresponding to typically applied (Ioannidis and Kulik, 2015) values of around 10-30%.
Considering overall performance on the 20 TM complexes, incorporating any exact exchange reduces the mean absolute error (MAE) of spin-splitting energies with respect to WFT. From PBE to 20% to 40% exchange, the MAE decreases monotonically from 17 to 12 to 8 kcal/mol for hexa-ammines and comparably (i.e., 16 to 9 to 6 kcal/mol) for hexa-aquas (Figure 4). Increasing HF exchange higher than 40% increases the MAE again for both ligand fields. For these two weak ligand fields, the optimal exchange fraction is more metal dependent than ligand field dependent, producing comparable optimal exchange values for fixed metal and oxidation state (Figure 4). Recent work (Wilbraham et al., 2018) has suggested double hybrids (DH, i.e., with MP2 long-range correlation) could improve predictions of spin-state ordering. We selected the PBE0-DH, which contains 50% global HF exchange for comparison to GGA global hybrid results (Figure 4). The PBE0-DH results are comparable to those obtained with a modified PBE0 global hybrid GGA with 40% exchange (i.e., MAE of 8 kcal/mol for hexaammine complexes and 6 kcal/mol for hexa-aqua complexes, see Supplementary Table 4). Given the higher computational cost and scaling of the double hybrids, tuned GGA hybrids would remain a preferable choice for modeling larger systems.

Comparison of Nitrogen-Containing Ligands
So far, we have studied model complexes and confirmed the importance of incorporating exact exchange in DFT functionals to reproduce correlated WFT reference spin-splitting energetics. We next considered the transferability of exchange sensitivity of DFT spin-splitting energetics to octahedral transition metal complexes that contain coordination environments similar to Fe/N-doped graphene SAC models. The two ligands we used to represent these environments were pyridine (py) in which the coordinating nitrogen is in a six membered carbon-containing ring and pyrrole (pyr) in which the coordinating nitrogen is in a five membered carbon-containing ring (Figure 5 inset and see Computational Details). With a pure PBE GGA, all complexes except for hexa-pyr Fe(II) are low spin, whereas exchange fractions above ∼15% instead result in all ground states being assigned as high spin (Figure 5). The trend with exchange is again linear as in the hexa-ammine complexes, although the linearity is slightly reduced for Fe(II) vs. Fe(III) complexes ( Figure 5).
The E H−L values are within ∼7 kcal/mol for all complexes of the same oxidation state across the full range of exchange due to the relatively similar ligand field strengths of the Ncontaining ligands and comparable metal-ligand bond lengths (Supplementary Table 5). However, the hexa-ammines are generally the most low-spin-favoring, whereas the SAC-like nitrogen complexes have a slightly increased high-spin bias. In analogy to ligand field arguments for E H−L , exchange sensitivities, ∂ E H−L ∂a HF , of spin splitting are also comparable for ammonia, pyridine and pyrrole complexes (Figure 5). For , whereas pyridine has slightly larger slope of −92 kcal/(mol·HFX) in both oxidation states. These observations are consistent with the previously observed greater sensitivity to the ligand identity than to oxidation state for Fe(II)/Fe(III) (Ioannidis and Kulik, 2015). Due to the similar spin-splitting FIGURE 5 | Sensitivity to HF exchange fraction (a HF ) of high-spin/low-spin splitting ( E H−L , in kcal/mol) for Fe(II) (squares) and Fe(III) (circles) homoleptic octahedral transition metal complexes with NH 3 (red symbols), pyridine (blue symbols), and pyrrole (green symbols) ligands. The structures of pyridine and pyrrole compounds are shown in inset in ball and stick (gray carbon, blue nitrogen, white hydrogen, and orange for iron). A zero axis is shown that indicates change in favored ground state spin.
energetics and sensitivities to exchange, we expect that our observations on hexa-ammine complexes are applicable to pyridinic and pyrrolic nitrogen-containing complexes and materials as well. Thus, for larger SAC models, we recommend either typical exchange fractions for qualitative ground state spin assignment (i.e., high spin) or higher exchange fractions (ca. 40-50%) that were needed in the hexa-ammine cases to reproduce WFT results quantitatively.

Graphene Flake Models of SACs
Fe/N-doped graphitic SACs are expected on the basis of experimental spectroscopic characterization (Zitolo et al., 2015;Liu et al., 2017) to consist of Fe metal centers coordinated by pyridinic or pyrrolic nitrogen atoms. These experimental observations come from a combination of aberration corrected scanning tunneling electron microscopy to confirm well-isolated metal sites as well as numerous spectroscopic techniques (e.g., X-ray absorption spectroscopy) to confirm the metal coordination (Zitolo et al., 2015;Liu et al., 2017). Although the most reactive SAC active site remains an open question, we consider in this work two limits in finite graphene flake SAC models that contain either four pyridinic (FeN 4 C 10 ) or pyrrolic (FeN 4 C 12 ) N atoms (Figure 6). Due to the rigidity of the graphene flakes, singlet Fe(II) pyridinic (py) and pyrrolic (pyr) Fe-N bond lengths are shorter than the corresponding py or pyr octahedral complexes (py: 1.90 vs. 2.08 Å or pyr: 1.97 Å vs. 2.11 Å, see Supplementary Tables 5-6). This rigidity in the SAC models without any displacement of iron from the plane that has been observed in porphyrins (Sahoo et al., 2015) also leads to average Fe-N distances being invariant to spin state. There is a marginal (ca 0.01 Å) increase from singlet to quintet spin states for the FeN 4 C 10 SAC in comparison to large (ca. 0.16-0.20 Å) bond length increases from singlet to quintet in the Fe(II)(py) 6 complex (Supplementary Tables 5-6). The Fe-N distances are shorter in pyridinic FeN 4 C 10 than in pyrrolic FeN 4 C 12 (i.e., 1.90 Å vs. 1.96 Å) due to smaller N-N separations (FeN 4 C 10 : 2.61 Å and 2.75 Å, FeN 4 C 12 : 2.75 Å). Despite this difference in N-N separation, which has previously been noted to influence experimental spin state ordering (Phan et al., 2017), the PBE0 ground state spin is triplet in both models (Supplementary Table 7).
The two singly occupied orbitals in both triplet SAC models correspond to the d xz and d z 2 orbitals consistent with expectations for square-planar coordinated triplet Fe(II) (Figure 6). Small differences are observed in the orbital character due to the lower symmetry for the pyridinic FeN 4 C 10 flake: in this case, d xz and d yz degeneracy is broken (Figure 6). The longer N-N separation along the x-axis leads to a pure d xz orbital vs. d xz and d yz mixing for the case of pyrrolic FeN 4 C 12 (Figure 6). In both cases, weak coupling is observed between the metalcentered orbitals and p-orbitals of both the N and C atoms in the graphene flake (Figure 6).
For both SAC models at the PBE0 level of theory, singlet and quintet states reside ∼4-6 kcal/mol and 12 kcal/mol above the triplet ground state, respectively, and thus singlet-quintet E H−L is around +6-8 kcal/mol (Supplementary Table 7). These observations contrast with the octahedral models: E H−L for Fe(II)(py) 6 is −15 kcal/mol and is −30 kcal/mol for Fe(II)(pyr) 6 (Supplementary Table 4). These differences can be traced to several factors, including the coordination number (4 vs. 6) in the models as well as rigidity of the graphene flakes that compress Fe-N bonds to values more commensurate with equilibrium low-spin geometries. Finally, examining the electronic structure of the SAC models reveals distribution of spin not just on the metal but also on the flake in the high spin states, particularly for the pyrrolic FeN 4 C 12 models (Supplementary Tables 8, 9). Even in the triplet ground state this is apparent with a magnetic moment of 2.2 µ B , close to that expected (i.e., 2 µ B ) for FeN 4 C 10 , but with a larger 2.7 µ B on Fe for FeN 4 C 12 (Supplementary Tables 8, 9). The singlet FeN 4 C 12 is also open shell with a 1.4-2.0 µ B magnetic moment on Fe (Supplementary Table 9). In contrast with the molecular complexes, the quintets are particularly poorly described by a localized metal spin, with a reduced Fe moment vs. the triplet state of around 2.4 µ B on Fe for FeN 4 C 12 and a comparable one of 2.2 µ B on Fe for FeN 4 C 10 . Although spin contamination can be expected with increasing HF exchange fraction, comparison of these moments across 0-50% exchange does not ever produce a pure 4 µ B moment on Fe for the quintet FeN 4 C 12 (Supplementary Table 9). This observation could be due to lowlying unoccupied states on graphene that are populated instead of the metal states, especially in these models and at this level of theory, as is known to occur in porphyrins as well (Fujii, 2002).
Beyond PBE0 (25%) results on the graphene SAC models, we considered properties over a range that spans from typical values in periodic catalysis modeling (i.e., 0%) to larger values (40%) motivated by our octahedral complex studies (see section Spin State Ordering in Model Complexes). Over this 0-40% range of exchange fractions, singlet and triplet states become destabilized with respect to high-spin quintet states (Figure 7). The reduced dependence of spin-state ordering observed here on exchange fraction in comparison to the octahedral complexes is due to differences in coordination number and rigidity of the SAC models. Intermediate Fe spin states were observed experimentally (Liu et al., 2017) for N-doped graphitic SACs using Mössbauer spectroscopy. Although high HF exchange fractions favoring quintet states for both SAC models would suggest inconsistencies with experiment, it is important to recall that the magnetic moments of the Fe metal are intermediate in both triplet and quintet states. Therefore, high HF exchange fractions are in fact stabilizing the simultaneous presence of spin on the graphene coupled to an intermediate Fe center (Supplementary Tables 8, 9). In both FeN 4 C 12 at low exchange fractions and FeN 4 C 10 over a larger range of 0-100% exchange, higher order than typically linear sensitivities are observed to % exchange (Supplementary Figure 3). For FeN 4 C 12 , the antiferromagnetically coupled metal spin on Fe varies significantly (i.e., 0.6 µ B ) and discontinuously, leading to less smooth energetic variations (Figure 7). For the FeN 4 C 10 model where spin states are more well defined, the variations are instead linear over the expected range of HF exchange (Droghetti et al., 2012;Ioannidis and Kulik, 2015;Zhao and Kulik, 2018).
High-valent Fe(IV)=O intermediates are expected to be essential for catalytic transformations at N-doped graphitic SACs (Liu et al., 2017). Thus, we examined the spinstate-, model-, and exchange-fraction-dependence of reaction energetics for Fe(IV)=O formation. Here, we employed N 2 O as a model oxidant, but results are comparable when assuming the oxygen atom comes from triplet O 2 (Figure 8 and Supplementary Figures 4-6). Overall, pyridinic SACs form more stable oxo species across the range of HF exchange and FIGURE 7 | Spin-splitting energetics (in kcal/mol) for singlet-quintet (S-Q, red circles) or triplet-quintet (T-Q, blue circles) spin states vs. % HF exchange for pyridinic (FeN 4 C 10 , top) and pyrrolic (FeN 4 C 12 , bottom) SAC models. The 25% exchange in standard PBE0 is indicated as a vertical dashed line. spin states than pyrrolic SACs (Figure 8). Although activation energies would be needed to make firmer statements about relative active site model reactivity, the endothermic reaction energies for the intermediate spin FeN 4 C 12 above 20% exchange (ca. +10 kcal/mol at 40% exchange) suggest that the pyrrolic model could potentially be unreactive with N 2 O oxidant (Figure 8). Regardless of spin state or model, increasing exchange fraction makes formation of oxo intermediates less favorable due to the penalty for delocalization (Gani and Kulik, 2017).
Determination of ground state spin from DFT of the pristine Fe(II)/N SAC model and the Fe(IV)=O intermediate should be carried out with caution, noting that spin (ca. 0.5-1.0 µ B ) arises on the O atom in triplet and quintet states of both Fe(IV)=O SAC models but not in the singlet states (Supplementary Tables 8, 9). Although spin on the oxo species can be expected and even linked to catalytic efficiency (Liu et al., 2009;Quesne et al., 2014), these states cannot readily be described within a single Kohn-Sham determinant in DFT (Koch and Holthausen, 2015). Overall, spin state ordering of the bare Fe(II) SAC is largely preserved in the Fe(IV)=O intermediate, with singlet and triplet Fe(IV)=O states being weakly stabilized FIGURE 8 | Reaction energetics ( E rxn , in kcal/mol) for oxo formation from N 2 O oxidant vs. % HF exchange for pyridinic (FeN 4 C 10 , top) and pyrrolic (FeN 4 C 12 , bottom) SAC models. The reaction is shown in inset. In each case, singlet (red circles), triplet (blue circles), and quintet (green circles) oxo formation energies are shown. The 25% exchange in standard PBE0 is indicated as a vertical dashed line. by around 3-5 kcal/mol with respect to the bare Fe(II) case for the pyridinic SACs (Supplementary Figure 4). In the pyrrolic case, the opposite occurs, potentially due to the loss of spin on the ring in triplet pyrrolic Fe(IV)=O (Supplementary Table 9 and Figure 4). Nevertheless, further investigation of kinetic barriers is merited in future work, as close spin state ordering of both the Fe(IV)=O and pristine Fe(II) intermediates combined with comparable differences of around 5-10 kcal/mol between reaction energetics in each spin state could give rise to spin-state dependent reactivity with distinct product formation (Kamachi and Yoshizawa, 2003;Ji et al., 2015). Finally, it is noteworthy that at 0% exchange (i.e., pure PBE), singlet FeN 4 C 12 is predicted to produce a slightly more stable Fe(IV)=O than the triplet, the same ordering that is observed for FeN 4 C 10 albeit at −10 to −15 kcal/mol in the former case vs. −30 to −35 kcal/mol in the latter case (Figure 8). At the 40% exchange motivated by our careful CASPT2 characterization of Fe-N bonds (see section Comparison of Nitrogen-Containing Ligands), or even at the 25% exchange fraction motivated in stronger ligand field cases , conclusions are different. Namely, FIGURE 9 | Periodic structures for pyridinic and pyrrolic periodic SAC models in the triplet state with spin density shown. Positive spin density is shown in red and negative spin density is shown in blue, with an isosurface value of 0.03 e/Bohr 3 . All structures are shown in ball and stick representation with carbon in brown, nitrogen in light blue, and iron in orange.
at these higher exchange fractions: (i) the triplet oxo is more stable for pyrrolic SACs than the singlet, whereas the ordering remains the same for the pyridinic case, and (ii) neither form exothermically for the pyrrolic case at these exchange fractions.

Periodic Modeling of SACs
We validated our choice of finite SAC flakes by comparing to periodic models of both pyridinic and pyrrolic SAC active sites (Figure 9). We focused on the triplet intermediate spin state favored both in our finite models and in experiment (Liu et al., 2017). Shorter Fe-N 1.91 Å vs. 1.97 Å Fe-N bond lengths are observed for the pyridinic model than for the pyrrolic models, consistent with the finite models. It is more straightforward to localize the magnetic moment to the metal in these periodic systems than was observed for the molecular models (Supplementary Table 10 and Supplementary Figure 7). For the triplet cases, spin density is nearly exclusively observed on the metal center (Figure 9). In both pyridinic and pyrrolic periodic models, use of the hybrid functional leads to less electron density on the Fe center than when a GGA is employed, consistent with prior observations Zhao and Kulik, 2018) (Supplementary Table 10).
We compared the electronic structure of the pyridinic and pyrrolic SAC active site models by determining the projected density of states (PDOS) decomposed by N 2p, C 2p, and Fe 3d contributions with the HSE06 hybrid functional (Figure 10). Qualitatively, the occupied orbitals and symmetries confirm observations made on the finite graphene flake models with PBE0, which may be expected if short range effects dominate as HSE06 and PBE0 both incorporate 25% exchange in the short range mixed with pure PBE exchange. That is, the pyridinic system again has singly occupied d z 2 and d xz orbitals and no occupation of the d xy state. The pyrrolic system also differs from the pyridinic by having degenerate d xz and d yz states, consistent with the molecular flakes (Figure 10). Generally, agreement is more variable for other spin states, where convergence of the magnetic state is sensitive to starting conditions in the periodic calculation (Supplementary Figures 8-9). In the pyridinic case, FIGURE 10 | The HSE06 projected density of states (PDOS) for spin up (Top) and down (Bottom, reflected curves) triplet pyridinic (Left) and pyrrolic (Right) SAC models. The 3d Fe orbital PDOS are shown as indicated in inset legends as solid curves except for yz and xz in the pyrrolic case, which are shown as dashed lines due to their degeneracy. The average PDOS for a 2p orbital from C (brown) or N (dark blue) are shown as translucent shaded regions. All energy levels (in eV) are aligned to the Fermi level (E F ), which is shown as a vertical dashed line. Some d levels have been truncated by the y-axis range to be able to compare to the broader C 2p and N 2p features. d yz states span the Fermi level, whereas the 3d states are well separated in the pyrrolic SAC (Figure 10). The pyrrolic 3d states also mix more deeply into the C and N 2p bands, whereas in the pyridinic case, most 3d states sit at the top of the occupied C/N 2p bands (Figure 10). Overall, these observations support the use of finite models at higher levels of theory for consistent modeling, due to the unique challenges of modeling periodic systems with such methods . More analysis in larger supercells with variable graphene defects will be necessary in future work to strengthen this conclusion.

CONCLUSIONS AND OUTLOOK
We have presented an overview of the effect of computational model choice on the properties of octahedral transition metal complexes and emergent single atom catalyst (SAC) materials made from Fe centers in N-doped graphene. The octahedral transition metal complexes chosen mimic the ligand field environment observed in the SAC models but remain tractable for study with multi-reference wavefunction theory. Observations from the hexa-aqua and hexa-ammine complex studies revealed that spin state ordering of mid-row complexes could be sensitive both to the IPEA shift chosen and whether an extended active space was used, whereas late and early transition metals were far less sensitive. Then using these extended active space CASPT2 results as a benchmark, we observed that nearly all transition metal complexes benefitted from increased HF exchange. In fact, errors with respect to exchange fraction chosen monotonically decreased from 0 to 40% but with no improvement for higher exchange fractions. The 40% exchange hybrid results had comparable accuracy to the more computationally demanding double hybrid PBE0-DH.
The HF exchange tuning study confirmed the comparable behavior of Fe(II) complexes with ammonia, pyridine, and pyrrole ligands due to the overriding role of the metal, oxidation state, and ligand connecting atom in determining functional sensitivity. Comparison to CASPT2 results on the hexa-ammine system motivated us to propose higher exchange fractions (ca. 40% rather than 25% in PBE0) to be essential to counteract the low-spin bias in semi-local DFT. Using these benchmarks, we then evaluated the effect of DFT functional tuning on finite graphitic SAC models with pyridine or pyrrole nitrogen atoms. In these square planar SAC geometries, rigid structures compressed the Fe-N bond length, reducing exchange sensitivity of spin state ordering but otherwise confirming the observations in the octahedral complexes (i.e., favoring higher spin configurations and destabilizing singlet states). We observed that at the recommended higher exchange fractions, the formation of an oxo intermediate from N 2 O became unfavorable at the pyrrolic SAC active site. We also observed changes in spin state ordering of the most stable oxo intermediates. Thus, incorporation of exchange can alter predictions of reactivity at SAC active sites. Finally, we confirmed that these observations were not sensitive to choice of a finite SAC model by comparing to periodic SAC models where similar electron configurations were observed.
Overall, predictions of reactivity and spin state ordering are highly sensitive to the functional employed. We have shown that this sensitivity is broadly transferable across different ligand environments as long as the metal and direct ligating atom are kept constant. This observation can be leveraged to obtain DFT functional performance on smaller models where correlated WFT is tractable. Additionally, smaller flake models may be amenable to direct WFT calculation with methods not covered in this work. Spin state exchange sensitivity can be expected to be depressed in cases where the SAC is fully rigid and prevents expansion of the metal-ligand bond when the spin state changes. The highly ordered, symmetric cases here are expected to be the limit in this rigidity argument, and more disordered SAC models or more flexible active sites (e.g., graphitic carbon nitride) are expected to be more sensitive. Still, some outstanding challenges remain in understanding the extent to which spin arising on the graphene itself could be physical and also impart reactivity or whether it arises due to increased static correlation error for hybrid DFT. Furthermore, in periodic simulations with larger supercells than studied in this work it will become essential to employ range-separated hybrids with HF exchange only in the short range, and the comparison to CASPT2 results to range-separated hybrid tuning will be essential here. Overall, modeling in SACs will continue to benefit from this multilevel approach in assessing method accuracy and sensitivity and how method choice impacts predictions of active site geometry and reactivity.

AUTHOR CONTRIBUTIONS
HK designed the research. FL, TY, JY, EX, and AB carried out the research. FL, TY, JY, and HK wrote and revised the manuscript.

FUNDING
The authors acknowledge support by the Department of Energy under grant number DE-SC0018096 for the work on density functional theory and the support of FL and AB. The authors also acknowledge the National Science Foundation under grant number CBET-1704266 for the support of TY. HK holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund.