All geometry optimizations were carried out in ORCA 83 v.4.0.1.2 and v.4.2.1 using density functional theory (DFT) with the generalized gradient approximation (GGA) global hybrid PBE0 and the def2-TZVP basis set and employing semi-empirical D3 dispersion with Becke-Johnson damping. The resolution of the identity (RI) approximation was used, and the auxiliary basis set for def2-TZVP was generated automatically by ORCA to accelerate calculations. Consistent with prior work , the five heavy atoms of succinate and the methyl carbon atoms of His ligands were held fixed to mimic the ligand positions in the enzyme (Supporting Information Figure ). Optimizations were performed in the gas phase and in solution using the BFGS algorithm in redundant internal coordinates. The default thresholds of 5 x 10 -6 hartree for the change in self-consistent field (SCF) energy between steps and 3 x 10 -4 hartree/bohr for the maximum gradient were used. The solution calculations were carried out using the conductor-like polarizable continuum model 94 (C-PCM) solvation energies along with the conductor-like screening solvent model (COSMO) epsilon function type. A solvent dielectric value of e = 10 was used, approximately mimicking the protein environment. As observed in prior work , we find that the inclusion of environment effects through implicit solvent alters most gas-phase geometries very little (i.e., most are within 0.01-0.03 Å while a few outliers have larger changes), thus we carry out further analyses only in the gas phase (Supporting Information Table ). Low-spin (LS) singlet states were simulated in a spin-restricted formalism while non-singlet LS doublet, intermediate-spin (IS), and high-spin (HS) states were simulated in an unrestricted formalism (Supporting Information Table ). All initial and optimized structures for gas-phase and implicitly solvated optimizations are provided in the Supporting Information.