We performed checks of the geometry 137 to determine the fidelity of the optimized geometries at different HF exchange fractions (Supporting Information Table ). We also checked for significant spin contamination, which we define as cases where the expectation value of the S 2 operator deviates from the anticipated eigenvalue (i.e., S(S + 1)) by greater than or equal to 1 µB 2 . The largest spin deviation observed in our systems was small (i.e., 0.23 µB 2 ). We also checked that the spin on the metal in the SCO complexes was close to the expected value (i.e., the number of unpaired electrons on Fe is within one of the total number of unpaired spins) to ensure that the calculation converged to the correct electronic state. For systems with spin not localized to the metal, calculations were reconverged using wavefunctions from calculations with different HF exchange fractions. After this procedure, all calculations on the 95 SCO complexes were successful, with good geometries and suitable electronic structure. All geometry optimization calculations were performed in both the gas phase and the presence of implicit solvent. Implicit solvent calculations were carried out with the conductor-like polarizable continuum 138 solvation model 139,140 model (C-PCM) as implemented in TeraChem 143 , with a dielectric constant of ε = 10.3 chosen to approximately model the effect of the crystal field. The solute cavity was built for all the elements using the defaults available in TeraChem, i.e., 1.2 x Bondi's van der Waals radii 144 . The gas-phase calculations resulted in some unreasonable geometries that failed the fidelity checks, and hence, we only consider structures obtained from geometry optimizations carried out in C-PCM (Supporting Information Table and Figures ). Following the geometry optimizations, frequency calculations were carried out in TeraChem to obtain free energy corrections for a temperature of T = 300 K using the corresponding optimized geometries and converged wave functions. We chose T = 300 K instead of the temperature of crystallization to simplify the problem due to the large number of distinct crystallization temperatures and the weak dependence of the crossover temperature on the specific temperature chosen to calculate the correction. All initial and optimized structures are provided in the Supporting Information.