Ligand additivity relationships enable efficient exploration of transition metal chemical space

To accelerate exploration of chemical space, it is necessary to identify the compounds that will provide the most additional information or value. A large-scale analysis of mononuclear octahedral transition metal complexes deposited in an experimental database confirms an under-representation of lower-symmetry complexes. From a set of around 1000 previously studied Fe(II) complexes, we show that the theoretical space of synthetically accessible complexes formed from the relatively small number of unique ligands is significantly (ca. 816k) larger. For the properties of these complexes, we validate the concept of ligand additivity by inferring heteroleptic properties from a stoichiometric combination of homoleptic complexes. An improved interpolation scheme that incorporates information about cis and trans isomer effects predicts the adiabatic spin-splitting energy to around 2 kcal/mol and the HOMO level to less than 0.2 eV. We demonstrate a multi-stage strategy to discover leads from the 816k Fe(II) complexes within a targeted property region. We carry out a coarse interpolation from homoleptic complexes that we refine over a subspace of ligands based on the likelihood of generating complexes with targeted properties. We validate our approach on 9 new binary and ternary complexes predicted to be in a targeted zone of discovery, suggesting opportunities for efficient transition metal complex discovery.


Introduction.
In recent years, virtual high-throughput screening (VHTS) 1-8 with first-principles density functional theory (DFT) and machine learning (ML) models [9][10][11][12][13][14][15][16][17] has greatly accelerated the discovery of new molecules and materials [18][19][20][21][22][23] . Nevertheless, the theoretical space of all possible compounds or materials is so large as to challenge even the most accelerated methods, with 10 30 to 10 60 theoretical drug-like molecules being enumerated [24][25][26][27][28][29] from a relatively small number of elements and atoms. 30,31 Within the space of theoretical transition metal complexes, additional variables emerge, such as the metal identity, spin, and oxidation state, as well as denticity of the ligands. 4,9 Indeed, significant analysis has been carried out by Fey and coworkers in understanding the role of privileged (e.g., phosphine) ligands in determining transition metal complex properties. [32][33][34] Jensen and coworkers have devised elegant strategies to explore the space of favored complexes, e.g., by adjusting the denticity during complex optimization carried out with efficient semi-empirical or force-field based scoring. [35][36][37][38] One key challenge for high-throughput screening with density functional theory of transition metal complex space is that the smallest nontrivial mononuclear octahedral complex consists of at least seven heavy atoms and nearly 100 electrons, also challenging the speed of conventional simulation techniques in comparison to readily computed data sets of closed-shell transition metal 39 and organic molecules 40,41 . Further compounding the challenges of exploring transition metal chemical space are potential issues with convergence success or the presence of multi-reference character. [42][43][44] Thus, it is attractive to identify the minimal set of explicit first-principles calculations that can be used to build a model of the properties of the full data set. [45][46][47] Recently, toward the goal of exploring a more diverse transition metal chemical space in comparison to complexes of frequently studied ligands, we devised a strategy for enumerating hypothetical, small (i.e., 1-2 heavy atoms per coordination site) ligands for mononuclear octahedral transition metal complexes 48 . These ligands sampled a diverse combination of coordinating atoms and their bonding environments 48 , and only a small fraction were represented in prior databases of organic molecules. 25,49,50 We showed how incorporating these molecules could improve the fidelity of artificial neural network (ANN) models 48,51 when applied to larger, realistic complexes present in the Cambridge Structural Database (CSD) 52 . That study was limited to the properties of homoleptic combinations of those ligands (i.e., all ligands equal), and therefore did not capture effects of mixing ligands that give rise to the compelling properties of many heteroleptic complexes and catalysts. Enumerating combinations of these ligands, however, would give rise to a combinatorial explosion, motivating strategies to understand which combinations are likely to be valuable or informative for a specific application.
Despite the challenge of combinatorial explosion, there are some established precedents of ligand additivity 53 that suggest that the properties of heteroleptic complexes can be inferred from combinations of homoleptic complexes 54 . For example, ligand additivity has been demonstrated in force field and DFT energetics 54 as well as DFT errors 55 . It has also been used in correction schemes, such as the DBLOC method. [56][57][58][59] We also recently exploited additivity to learn the degree of multireference character in a complex from the multireference character in its constituent ligands 60 . Additivity is also exploited heavily in fragmentation methods 61,62 and in local correlated methods 63,64 . In the present work, we carry out a survey of the symmetry classes and ligand diversity present in the CSD to confirm that the theoretical chemical space is orders of magnitude larger than the number that have been characterized. Motivated by the need to devise efficient but accurate methods for exploration of chemical space, we introduce improved interpolation schemes for heteroleptic compounds to incorporate cis and trans effects. Finally, we demonstrate how these approaches can be used for efficient but accurate discovery of transition metal complexes with targeted properties.

2a. Curation from the Cambridge Structural Database.
A set of 85,575 mononuclear octahedral transition metal complexes was curated from the Cambridge Structural Database 52 (CSD) version 5.41 (Nov 2019). This procedure employed both the Conquest graphical interface to the CSD as well as the Python application programming interface, in all cases applied to the v5.41 data set with complexes from the November 2019 data set with both the March 2020 and May 2020 data updates (Supporting Information Text S1). For the complexes identified as octahedral, equatorial planes and axial positions were assigned based on prior reported rules. 65 To identify the symmetry of the ligands, unique ligands were identified by removing the metal atom to create independent molecular graphs for each ligand. Each ligand was identified as chemically unique within a given octahedral complex if it differed from all other ligands in the complex by: 1) heavy atom chemical symbols, 2) metal-connecting-atom element or 3) more than three hydrogens (Supporting Information Text S2). The symmetry of the complex was identified by distinguishing ligand denticity overall and in the equatorial plane along with the total number of unique ligands and whether ligands that were trans to each other were identical (Supporting Information Text S2). This led to a nomenclature for 66 ligand symmetry classes (Supporting Information Text S2).
To identify the set of unique ligands in each complex, a dummy atom with identical connectivity to the metal with atomic number of 0 was introduced to preserve the connectivity of the ligands to the metal without preserving metal identity. For this ligand and dummy atom combination, the atomic-number and bond-order weighted connectivity matrix determinant was calculated as described in Ref. 65 . We also computed the determinant of the atomic-number weighted bond-order weighted connectivity matrix where the off-diagonal elements, ZiZj (i ≠j), were set to the CSD-derived bond order for each ligand. Ligands with both distinct atomic-numberweighted connectivity matrix determinants and bond-order-weighted connectivity matrix determinants were identified as distinct ligands across monometallic transition metal complexes in the CSD. A second search was carried out by requiring that oxidation states and charges be assigned by the user along with no obvious disorder or missing hydrogen atoms in the structure, leading to 17,085 unique "computation-ready" complexes. Finally, we curated a subset of 1,202 Fe(II)-containing "computation-ready" complexes, based on the oxidation state reported by the user. From the ligands identified in this Fe(II) complex set, heteroleptic calculations from CSD ligands were carried out using a previously developed procedure 60 that enabled the assignment of the per-ligand charge.

2b. Electronic Structure Calculations.
DFT geometry optimizations were performed using a development version of TeraChem v1.9. 66, 67 The B3LYP 68-70 global hybrid functional was employed with the LANL2DZ 71 effective core potential for transition metals and the 6-31G* basis 72 for all other atoms. All transition metal complexes were studied with Fe(II) centers in low-spin singlet and high-spin quintet multiplicities.
Singlet calculations were carried out in a spin-restricted formalism, while quintet calculations were unrestricted. Level shifting 73 was employed to aid self-consistent field convergence with the majority-spin and minority-spin virtual orbitals each shifted by 0.25 Ha. Geometry optimizations were carried out in the gas phase in translation rotation internal coordinates 74 using the BFGS algorithm. Default tolerances of 4.5x10 -4 hartree/bohr and 10 -6 hartree were applied in the convergence criteria for the maximum gradient and energy difference between steps, respectively. For the representative model complexes of CH3CN, H2O, CO, and NH3, initial structures were generated with molSimplify 51,75,76 , which uses OpenBabel 77, 78 as a backend. The same protocol was applied to generate homoleptic complexes of the 20 neutral ligands derived from monodentate-only, non-homoleptic Fe(II) complexes obtained from the CSD 52 after discarding one bulky ligand, OP(Ph)3, that could not form a stable homoleptic structure for steric reasons. For the 36 homoleptic Fe(II) complexes in the CSD, the structures were directly extracted for subsequent geometry optimization. Heteroleptic complexes of the 12 representative ligands were also generated with molSimplify and optimized following the same procedure. A pickle file of the curated CSD database and structures and properties of all transition metal complexes studied with DFT is provided in the Supporting Information .zip file.

3a. Symmetry Classes and Theoretical Complex Space.
The diversity present in the chemical space of transition metal complexes is derived from variability in the metal, its oxidation and spin state, as well as the chemistry of the coordinating ligands. Our experimental knowledge of this chemical space is unevenly distributed. We thus first examined the structures deposited in the CSD to uncover trends in the arrangement of ligands in previously characterized complexes (see Computational Details). From 85,575 mononuclear octahedral transition metal complexes of which 17,085 are identified as unique and computationready (e.g, have user-defined charges), the vast majority (95-98%) contain no more than three unique ligands (Figure 1 and Supporting Information Table S1). In fact, 28% of all unique computation-ready complexes are homoleptic, and a majority (76%) contain no more than two ligand types (Figure 1 and Supporting Information Table S1). Because metal identity and ligand diversity are expected to be coupled, we also evaluated statistics on a subset of 1,202 unique Fe(II) complexes and confirm that the preference for complexes with no more than two ligands is preserved and even strengthened (Figure 1 and Supporting Information Table S1).

Figure 1.
Stacked, unnormalized histogram of the number of complexes in the CSD grouped by number of unique ligand types and by the highest denticity of ligands in the complex (monodentate in brown, bidentate in gray, tridentate in orange, tetradentate in blue, pentadentate in green, and hexadentate in red, as indicated in inset legend). These counts are shown for all complexes (top), unique and computation-ready complexes (middle), and the Fe(II) unique, computation-ready subset (bottom).
To identify how the structures sampled in the CSD compare to the theoretical space of all hypothetical complexes, we enumerate the overall pool of theoretical complexes and the number in each symmetry class to compare to the most frequently characterized symmetry classes in the CSD. We applied Pólya's enumeration theorem to octahedral coordination geometries to obtain all possible symmetry classes from the cycle index (i.e., theoretical sum) of a symmetry group. [79][80][81] The cycle index for all possible octahedral complexes is where the subscript, j, on f indicates the number of unique sites (e.g., 2 for M(L1)5(L2)1) up to a theoretical maximum of 6. Each individual term in the cycle index is computed as: where N is the total number of unique ligands that can occupy the sites, and each xi corresponds to one unique ligand set choice.
For a worked example of N = 12 unique ligands, we determine the theoretical number of complexes, the cardinality, as  (Table 1).    Table 1). The less degenerate CA/TA complexes form a total of N(N-1)(N-2)/2! complexes for each of the two isomers (i.e., 1320 for N = 12, Table 1). Similarly, FA/MAC/MAT complexes can each form N(N-1)(N-2)/2! complexes in two isomers each (i.e., 3960 for N = 12, Table 1).  Returning to the diversity observed in complexes deposited in the CSD, we can qualitatively observe that relatively little of the theoretical space has been sampled. As we have shown for a representative example, for any set of unique ligands a much larger theoretical number of binary and ternary complexes can be formed in comparison to homoleptic complexes.
Nevertheless, there are far fewer binary and especially ternary complexes in the CSD (Supporting Information Table S1). Within the binary complexes, 5+1 and TS complexes are overrepresented in comparison to those for FS, MS, or CS based on our theoretical enumeration ( Figure 4 and Supporting Information Table S2). All ternary complexes are underrepresented, but those with  Table S2). If we simplify our analysis by focusing only on Fe(II) complexes, the same trends hold, although Fe(II) complexes have an even greater relative number of 5+1 and TS complexes and lower number of CS complexes ( Figure 4 and Supporting Information Table S2). An additional factor in analysis of symmetry classes is the extent to which ligand denticity plays a role. While enumeration is straightforward for monodentate ligands, higher-denticity ligands (e.g., pentadentates) may only be compatible with some of the symmetry classes. Indeed, over both the computation-ready CSD set and the Fe(II) subset, a significant number of complexes consist of higher-denticity ligands that would be incompatible with full enumeration (Figure 1 and properties from homoleptic properties are necessary to "fill in" the remainder of this unexplored space.

3b. Ligand Additivity for Interpolating Properties of Transition Metal Complexes.
We next aimed to determine the extent to which the properties of lower-symmetry heteroleptic transition metal complexes could be inferred from those of higher-symmetry complexes. We constructed complexes with two to three unique ligand types from small sets of ligands that spanned a large range of ligand field strengths: weak-field water, strong-field carbonyl, and either strong-field methylisocyanide or weak-field ammonia. For calculation of the adiabatic high-spin (e.g., quintet Fe(II)) to low-spin (e.g., singlet Fe(II)) splitting, DEH-L, it can be expected that the weak-field ligands will lead to homoleptic complexes that favor high-spin states, whereas strong-field ligands will make homoleptic complexes that favor low-spin states. Thus, heteroleptic combinations of these ligands are expected to reside between the two limits. One simple way to obtain estimates of the spin splitting of the heteroleptic complexes (e.g., with up to three unique ligand types) is to take a weighted average of the spin splitting of the parent homoleptic complexes: Indeed, we observe that heteroleptics reside between the homoleptic limits for spin splitting, and the linear averaging roughly holds for complexes with methyl isocyanide (CH3CN), H2O, and CO ligands ( Figure 5 and Supporting Information Table S7). Similar observations can be made on combinations with H2O, CO, and NH3 (Supporting Information Figure S1). Nevertheless, there are significant outliers in the interpolated versus actual DEH-L, which are particularly evident when comparing heteroleptic complexes with the same stoichiometry, and therefore the same prediction from a simple linear model, but distinct ligand symmetry ( Figure 5). For example, the CS complex Fe(II)(H2O)4(CO)2 DEH-L is predicted accurately (predicted: -7.9 kcal/mol vs calculated: -9.9 kcal/mol) from homoleptic interpolation ( Figure 5 and Supporting Information Table S8). The same prediction significantly overestimates the TS complex with the same stoichiometry (predicted: -7.9 kcal/mol vs calculated: -19.3 kcal/mol), which instead behaves much more similarly to the 5+1 complex, as observed in prior work 82 ( Figure 5). Overall, for this combination of ligands, CS or FS complexes that have one minority ligand in the axial position and another in the equatorial plane are much better predicted than the equivalent MS or TS complexes ( Figure 5).
For the case of CO, H2O, and NH3, mixing between weak-field NH3 and H2O is relatively accurately predicted from homoleptic averaging, whereas for NH3 and CO, it is the MS and TS complexes that are more accurately predicted than the FS or CS counterparts (Supporting Information Figure S2 and Tables S8-S9). Because adiabatic spin splitting involves geometry optimizations in two distinct spin states, we also evaluated HOMO levels of each singlet complex as a property that depends only on a single geometry. Overall, interpolation of HOMO energies of heteroleptic complexes of CH3CN, CO, and H2O from homoleptic complexes reproduces trends between the homoleptic limits ( Figure   6 and Supporting Information Tables S7, S10-S11). As with spin splitting, there are key differences for complexes with identical stoichiometry that cannot be captured by interpolation from homoleptic complexes alone ( Figure 6). Interestingly, for the same complexes for which the CS complex had a higher (i.e., more off-parity) DEH-L and the TS complex was more like the equivalent 5+1 complex, we find more varied results for the HOMO level, with some cases occurring where the TS and FS complexes are less accurately predicted ( Figure 6 and Supporting Information Table S10). Nevertheless, over the full range of data, the outliers in the HOMO level prediction are more modest than was observed for DEH-L. These trends also hold for the complexes with NH3, H2O, and CO (Supporting Information Figure S2 and Tables S10-S11). Given the differences between CS and TS complexes despite having identical stoichiometry, we next identified strategies for improving the interpolation. First, we employed CS and TS complexes along with homoleptic complexes to predict the spin-splitting energetics of the other binary heteroleptic complexes as follows: where E(5+1) is the interpolated energy of the M(L1)5L2 complex from the M(L1)6 HO energy and the M(L1)4(L2)2 TS energy.
Similarly, we estimate the FS and MS energies as: In practice this corresponds to computing six energies and interpolating four remaining energies for each pair of ligands for a computational savings of 40%. Indeed, we observe reduced mean absolute error (MAE) over the remaining points that are interpolated (2.4 kcal/mol vs 5.1 kcal/mol) in comparison to the HO-only interpolation ( Figure 5 and Supporting Information Tables S8-S9).
In particular, FS complexes of CO and H2O are now correctly predicted to be much more lowspin-directing than the MS complex of the same stoichiometry ( Figure 5). The HOMO levels for these and other complexes are also improved (Supporting Information Tables S10-S11). For example, the modified interpolation is able to capture the fact that MS Fe(II)(CO)4(H2O)2 has a shallower HOMO level than its FS counterpart ( Figure 6). Nevertheless, not all points are uniformly improved by this interpolation. For the case of NH3, CO, and H2O where the interpolation already performed well, some points such as the DEH-L for 5+1 Fe(II)(CO)5NH3 are slightly worsened in the modified interpolation (Supporting Information Figure S1). For the same complex, the HOMO level is equivalently predicted by both interpolation schemes, and most HOMO level estimates are improved (Supporting Information Figure S2). Overall, errors are on average significantly lower for all properties and sets of ligands considered when the modified interpolation expressions are employed.
Although the CS/TS-derived interpolation schemes greatly reduce errors in estimating the energetics of heteroleptic complexes, they still require significant computational overhead. Thus, we next aimed to identify if FS and MS complexes, which contain three ligands cis to each other or two ligands cis and two sets of ligands trans, respectively, could be used instead in the interpolation (see Figure 2). If the FS and MS complexes impart sufficient information, using them in an interpolation scheme along with homoleptic complex properties corresponds to evaluating four properties (e.g., energies) to predict six properties for a computational savings of 60%. In this interpolation scheme, we estimated the complex properties from FS and MS complexes as follows: where E(5+1) is the interpolated energy of the M(L1)5L1 complex from the M(L1)6 HO energy and the M(L1)3(L2)3 FS energy. Similarly, we obtain expressions for the CS and TS complexes as: Indeed, using this approach we achieve errors only slightly larger than that for the CS/TS averaging scheme, with the added benefit of requiring fewer energies to obtain the same fidelity of the interpolation ( Figure 5). For example, the higher spin-splitting energy of CS complexes relative to TS complexes is captured here because it can be directly derived from the strength of the high-spin directing character of the FS complex relative to that of the MS complex ( Figure 5).
Mixing of the weak-field NH3 and H2O ligands, which had slightly worsened with CS/TS interpolation, is also significantly improved in this scheme, and the other sets of complexes are of comparable accuracy (Supporting Information Figure S1 and Tables S8-S9). While for the HOMO level, CS/TS interpolation tended to underestimate the HOMO level, FS/MS interpolation slightly overestimates HOMO levels, but errors are much smaller than for the HO-only interpolation ( Figure 6 and Supporting Information Tables S10-S11). For the set of ligands including ammonia, almost all points are predicted to comparable or slightly improved values (Supporting Information Figure S2). Thus, HO-only interpolation provides a highly efficient scheme for predicting heteroleptic transition metal complexes, but the best trade-off in accuracy and computational cost for interpolation is likely achieved through estimating properties using information from FS and MS complexes.
We next investigated whether we could generalize our observations to heteroleptics with three unique ligand types ( Figure 3). This extension is motivated by the fact that 96% of all mononuclear transition metal complexes in the CSD contain no more than three ligand types, and over 99% of Fe(II) complexes contain three or fewer ligand types (Figure 1 and Supporting Information Table S1). Expansion to three ligand types introduces 29 complex energies that need computation for any set of three ligands. Homoleptic-only averaging performs poorly here for complexes that are mixtures of CH3CN, CO, and H2O, with a large difference between CA and TA complexes of H2O and CO being treated completely equivalently in this scheme ( Figure 7).
Similar differences in TA and CA HOMO levels are also missed in this averaging scheme ( Figure   8). Generally, the CA complex spin-splitting energies are better predicted by the homoleptic averaging than the TA are for both ternary complexes with CH3CN and with NH3 ( Figure 7 and Supporting Information Figure S3 and Tables S12-S13). For the HOMO level, results are more varied, with the HOMO energies of the CH3CN ternary complexes being underestimated while those of NH3 ternary complexes are overestimated ( Figure 8 and Supporting Information Figure   S4 and Tables S14-S15).  expressions for these heteroleptics are as follows: ( where the third ligand in the EA complex is the one that is trans to itself so energies involving that  S4 and Tables S14-S16). Differences in CA and TA complex properties are much better predicted by both interpolation schemes, with the CS/TS scheme performing best for the combination that includes CH3CN (Figures 5 and 6). While this may be expected as a generalization of observations on the binary heteroleptics, it is noteworthy that trends in the fully equivalent stoichiometries (i.e., all three EA isomers, DCS, or DTS complexes) are reasonably well predicted by the improved schemes, whereas they were indistinguishable with simple linear interpolation because they differ solely by cis versus trans positioning effects ( Figures 5 and 6). Overall, errors for the interpolation scheme using FS/MS ligands are sufficiently low to warrant its use in chemical space exploration (Supporting Information Table S16). We have thus demonstrated how over a set of three ligands, nine explicit calculations can be used to interpolate the properties of 18 binary and 29 ternary complexes to around 2 kcal/mol accuracy in spin-splitting energies or 0.1-0.2 eV accuracy in orbital energy levels.

3c. Interpolation of Chemical Space from Experimentally Characterized Complexes.
Given the strength of the interpolative trends we observed, we chose Fe(II) complexes from the CSD to explore the potential of interpolation from previously synthesized complexes. For these 1202 complexes, the majority (661) are homoleptics, followed closely by binary complexes (407) and ternary complexes (129), with only six having more ligand types (Figure 4 and Supporting Information Table S4).We restrict our analysis of chemical space interpolation to monodentate ligands from complexes of each symmetry type that only contain monodentate ligands to simplify the interpolation process because higher-denticity ligands impose geometric constraints that make them incompatible with certain symmetry types. Nevertheless, we note that there are a significant number of monodentate ligands that only appear in combination with higher-denticity ligands (Supporting Information Tables S4-S5). From the set of monodentate-only complexes, we identified 40 unique ligands already present in homoleptic complexes along with 48 additional unique ligands present in monodentate-only binary and ternary complexes for which we could assign a charge following the scheme introduced in Ref. 60 . Of the set of 88 ligands, only 56 are assigned a neutral charge, closed-shell electronic structure, and deemed sterically feasible for homoleptic calculations (i.e., excluding only the neutral OP(Ph)3 ligand).
We then computed the adiabatic spin-splitting energies and singlet HOMO levels of all 56 homoleptic Fe(II) complexes (see Computational Details). This set is strongly biased toward highspin structures, with only a few ligands giving rise to low-spin ground states (Figure 9 and Supporting Information Figure S5 and Table S17). The majority of spin-splitting energies are in the range of -30 to 0 kcal/mol, whereas a minority of complexes are low-spin ( Figure 9). Only one complex Fe(II)(CO)6, which was present in our study in Sec. 3b, has a HOMO level deeper than -  complexes (see Table 1). Thus, we identified a way to use the HO-only interpolation to reduce the number of unique complexes required to achieve high fidelity within a target region.
First, we select a targeted region of DEH-L values in the range of -4.0 to 4.0 kcal/mol and HOMO levels in the range of -14 to -13 eV. This region was selected because no homoleptic complexes were present in this range, but interpolation of the space predicted that heteroleptic complexes would be found there. These complexes could be of interest in chemical discovery applications that target spin-crossover candidates (i.e., with near-degenerate spin states) with good oxidative stability (i.e., deep HOMO levels). Next, we identify 12 common ligands from the parent homoleptic complexes that are predicted to give rise most frequently to binary and ternary complexes in the targeted zone on the basis of the simple HO-only model ( Figure 10 and Supporting Information Table S18 and Figure S6). These ligands consist of common ligands from our original set (i.e., CO, H2O, CH3CN, and NH3) but also introduce new chemistry (e.g, Scoordinating dimethylthioformamide, DMTF, and 2-chloropyrazine, ClPyz, Supporting Information Table S18). From this set, only 132 additional FS/MS complexes need to be studied beyond the 12 homoleptic complexes we already computed, to infer 6,776 additional HOMO level or DEH-L properties at higher fidelity than HO-only averaging (see Table 1). We carried out geometry optimizations of these 132 complexes and used their properties to evaluate the revised interpolated HOMO level and DEH-L values over this subset. Evaluating this subset also provides a validation of the accuracy of the homoleptic averaging over a larger set of ligands. Over this set, we observe that errors are comparable to the earlier tests (i.e., MAE of 4 kcal/mol), and homoleptic averaging is generally a good predictor of DEH-L (Supporting Information Figure S7 and Table   S19). For HOMO level predictions, there are more points that differ from the interpolated values, leading to higher MAEs (0.8 eV) than observed over our representative test ligands (Supporting Information Figure S7 and Table S19).  Table S20). For these same complexes, the homoleptic averaging predicted HOMO levels very well, within around 0.1 eV (i.e., lower error than for DEH-L, Supporting Information Table S20).
Returning to the properties predicted from FS/MS-interpolation, we identified a total of three additional binary complexes predicted to be in the targeted zone from FS/MS-interpolation but not in the targeted zone according to HO-only interpolation and we computed their properties ( Figure 10). These complexes are 5+1 Fe(II)(MeCN)5(CO) and CS/TS Fe(II)(ClPyz)4(CO)2 ( Figure 11 and Supporting Information predicted to fall within the target zone, the difference between the HO-only to FS/MS interpolation schemes causes some regions of property space to be predicted to be enriched while others are predicted to be depleted ( Figure 10). Figure 11. Properties of 8 binary (red or green) and 6 ternary (orange or blue) complexes in validation set for HO-interpolation (circles) and FS/MS-augmented interpolation (squares) for DEH-L (in kcal/mol, left) and HOMO level (in eV, right). The targeted zone for each quantity is shown as a turquoise square, and three representative complexes are shown in inset with their symmetry class and the associated points are indicated with gray arrows. Structures are colored as follows: brown for Fe, gray for C, blue for N, white for H, green for Cl, red for O, and yellow for S. A dotted parity line is also shown.
As a control, we also identified two complexes predicted to be in zone by the HO-only interpolation but out of zone by FS/MS-interpolation, CS/TS Fe(II)(MeOH)4(CH3CN)2 (Supporting Information Table S20). The homoleptic averaging predicts the DEH-L value for both complexes to be -2.6 kcal/mol, while calculated DEH-L values are significantly lower (-9.1 and -13.9 kcal/mol). The FS/MS interpolation captures well the relative CS/TS energetics in both this case and the in-zone example and only slightly overestimates the DEH-L value (Supporting Information Table S20). Overall, FS/MS-interpolation MAEs for these binary complexes is low (2.5 kcal/mol for DEH-L and 0.1 eV for HOMO levels) and less than half that observed from HOonly averaging (Supporting Information Table S20).
As a more stringent test of our FS/MS interpolation scheme, we also selected six representative ternary complexes predicted to be in the targeted zone from FS/MS interpolation but out of the targeted zone when estimated with homoleptic interpolation (Figure 11 and Supporting Information Table S21). All complexes were generally predicted to be more high-spin favoring by homoleptic averaging than predicted from FS/MS interpolation, and several were also predicted to have deeper HOMO levels than the targeted region (Supporting Information Table   S21). Over this set, explicit DFT calculations show that five of the complexes have DEH-L values in the targeted zone and four have HOMO levels in the targeted zone (Supporting Information   Table S21). These include FA Fe(II)(ClPyz)3(CO)2(DMTF), which was predicted to be strongly  These results demonstrate that coarse interpolation of large spaces with homoleptic complexes can be followed up by improved interpolation using selected FS/MS compounds to identify the most promising binary or ternary complexes for explicit calculation. Overall, the demonstrated approach represents a data-efficient strategy to infer properties across large compound spaces with systematically improvable fidelity.

Conclusions
A large-scale analysis of the mononuclear octahedral complexes deposited in the Cambridge Structural Database revealed a propensity towards specific, higher-symmetry classes.
In addition, few complexes contained more than three unique ligand types. To assess the relative diversity of these complexes compared to the enumerated chemical space, we obtained expressions for the theoretical number of complexes of the five binary and eight ternary symmetry classes for octahedral complexes. We showed that even for a relatively small number of neutral, monodentate ligands present in Fe(II) complexes, the total theoretical space of 816,200 binary and ternary complexes far exceeded those that had been characterized in the CSD. Page S9 Table S11 MAEs for HOMO binary complex estimates from interpolation schemes Page S10 Figure S3 Interpolated vs calc'd DEH-L for TMCs with three ligands: NH3/H2O/CO Page S10 Table S12 Interpolated and actual DEH-L for ternary H2O/CO/CH3CN complexes Page S11 Table S13 Interpolated and actual DEH-L for ternary H2O/CO/NH3 complexes Page S12 Figure S4 Interpolated vs calc'd HOMO for TMCs with three ligands: NH3/H2O/CO Page S13 Table S14 Interpolated and actual HOMO for ternary H2O/CO/CH3CN complexes Page S14 Table S15 Interpolated and actual HOMO for ternary H2O/CO/NH3 complexes Page S15 Text S3 CS/TS-derived interpolation expressions for ternary complexes Page S16 Table S16 MAEs of HOMO and DEH-L for ternary complexes Page S16 Table S17 Properties of 52 homoleptic Fe(II) complexes Page S17 Figure S5 Scatter of HOMO and DEH-L for 52 homoleptic Fe(II) complexes Page S18 Table S18 12 most common ligands in target zone for Fe(II) complexes Page S18 Figure S6 Homoleptic-only interpolated target zone for Fe(II) complexes Page S19 Figure S7 Parity plot for predicted vs calculated 132 FS/MS Fe(II) complexes Page S19 Text S1. Description of mononuclear complex curation. The Conquest interface was used to query for structures containing a transition metal atom (here, groups 3-12, periods 4-6, excluding La) that forms at least 4 bonds with either p-block (here, the first five rows of groups 13−17), noble gas, or group 1 or 2 elements. The only constraint imposed on the search was that 3D coordinates were determined. This initial query produced 462,012 Refcodes for unique structures. The CSD Python API was then used to iterate through every "component" molecule in each x-ray crystallographic structure to identify all non-polymeric molecules with distinct molecular weight (within the same crystal structure) and only one transition metal atom. For each of the resulting 240,117 mononuclear transition metal complexes the CSD mol2 files were saved after adding any missing hydrogen atoms, preserving both the connectivity and bond orders present in the CSD. From this diverse set, 6-coordinate metal centers without edge-bound or sandwich-type ligands were evaluated for deviation of the angles (as defined in Ref. 1) between the metal-connecting atoms from the VSEPR ideal angles for octahedral, pentagonal pyramidal, or trigonal prismatic structures. Structures with lowest angular deviation from ideal octahedral molecular structure over other geometries were identified as octahedral.
Text S2. Description of complex symmetry identification.
Within each complex, a 3-hydrogen difference criterion was introduced to heuristically capture ligands that only differed in protonation state. If none of the equivalence conditions were met the ligands were identified as chemically identical. Ligand symmetry was identified by binning the octahedral structures by: (1) Ligand denticities (e.g., 5-denticity + 1-denticity) (2) maximum number of connecting atoms from a single ligand in the equatorial plane (e.g., 3denticity + 3-denticity has at most 2 connecting atoms from the same ligand in the equatorial plane -> fac binding) (3) total number of unique ligands (e.g., for a monodentate-only complex with 5 unique ligands, two must be duplicates) and (4) for lower denticities, counts of pairs of trans ligands that are identical (e.g., for a monodentate-only complex with 3 pairs of identical ligands, if only 1 pair of trans ligands are identical -> other unique pairs must be cis). A nomenclature for ligand symmetry was created to describe the observed symmetries following these rules (1) Denticities are denoted by numbers as the base of the ligand symmetry representation with higher denticities listed first. (2) If ligands are identical they are surrounded by "||" markers, with ligands with the highest number of identical copies of the same denticity listed first. (3) Tridentate ligands are identified as either fac-or mer-by f/m labels, and corresponding sets of three identical monodentate ligands are given equivalent f/m labels. (4) Pairs of identical monodentate ligands are identified as either cis-or trans-with c/t labels, pairs of non-identical monodentate ligands are given c/t labels only when there is one pair of non-identical monodentate ligands, and the ligand symmetry of identical pairs is ambiguous. Examples of ligand symmetry using this nomenclature are "|111111|" for a homoleptic monodentate-only complex, "4|11|t" for an equatorial tetradentate with trans identical monodentate ligands, and "2|11|c11" for a complex with a bidentate and 2 identical monodentate ligands in the equatorial plane with 2 chemically-distinct monodentate ligands in the axial positions. Table S1. Number of complexes grouped by number of unique ligands for all (all) mononuclear octahedral transition metal complexes, all unique complexes with user-defined charges and nondisordered structures ("computation-ready"), and the subset of unique complexes with Fe(II) centers (Fe(II)). The percent of each category is also shown.  Table S2. Number of complexes grouped by symmetry class for all mononuclear octahedral transition metal complexes, all complexes with user-defined charges and non-disordered structures ("computation-ready"), and the subset of complexes with Fe(II) centers. The percent of each category is also shown.  Table S3. Number of unique ligands of each denticity grouped by symmetry class for all complexes with user-defined charges and non-disordered structures.           Figure S4. Calculated vs. linearly interpolated HOMO level (eV) for singlet Fe(II) complexes with at least one each of three ligands: NH3, H2O, and CO. From left to right: interpolation from homoleptic complexes (HO only), interpolation using homoleptic complexes as well as cis symmetric and trans symmetric complex energies derived from pairs of ligands (CS+TS), or interpolation using homoleptic complexes as well as fac and mer symmetric complex energies derived from pairs of ligands (FS+MS). Points are colored according to the ligand with the highest stoichiometric coefficient: H2O (red circles), CO (gray squares), and NH3 (blue down triangles), or equal weight of all ligands (green up triangles), as indicated in inset legend. Key isomers are annotated. In all panes, a black dotted parity line is shown.
The interpolated values are shown as a 2D histogram colored from low (purple) to high (yellow) density, and the same data is shown as a kernel density estimate on the histogram panes.   Table S21. Ternary complexes with their chemical formula and symmetry along with DE H-L (in kcal/mol) and HOMO level (in eV) obtained through three approaches: homoleptic-only interpolation, FS/MS-augmented interpolation, and explicit calculation. The resulting calculated properties that are in the targeted zone are shown in bold.