Active Learning Exploration of Transition Metal Complexes to Discover Method-Insensitive and Synthetically Accessible Chromophores

Transition metal chromophores with earth-abundant transition metals are an important design target for their applications in lighting and non-toxic bioimaging, but their design is challenged by the scarcity of complexes that simultaneously have optimal target absorption energies in the visible region as well as well-defined ground states. Machine learning (ML) accelerated discovery could overcome such challenges by enabling screening of a larger space, but is limited by the fidelity of the data used in ML model training, which is typically from a single approximate density functional. To address this limitation, we search for consensus in predictions among 23 density functional approximations across multiple rungs of Jacobs ladder. To accelerate the discovery of complexes with absorption energies in the visible region while minimizing MR character, we use 2D efficient global optimization to sample candidate low-spin chromophores from multi-million complex spaces. Despite the scarcity (i.e., approx. 0.01\%) of potential chromophores in this large chemical space, we identify candidates with high likelihood (i.e.,>10\%) of computational validation as the ML models improve during active learning, representing a 1,000-fold acceleration in discovery. Absorption spectra of promising chromophores from time-dependent density functional theory verify that 2/3 of candidates have the desired excited state properties. The observation that constituent ligands from our leads have demonstrated interesting optical properties in the literature exemplifies the effectiveness of our construction of a realistic design space and active learning approach.


Introduction
Transition metal chromophores are an important design target because they play an important role in many chemical and biological processes ranging from natural light harvesting 1-3 , light-emitting technologies 4 to photocatalysis 5,6 . Due to the delicate interplay 7,8 required to tune complex properties, it is challenging to use a standard Edisonian approach 9 to simultaneously alter metal-ligand interactions, ligand field strength, electron donating/withdrawing effects, and the relative energetic positioning between the ground and excited state potential energy surfaces. A notable exception is a recent work 10 that utilized high-throughput experiments to identify heteroleptic Ir(III)-based chromophores. Nevertheless, to facilitate scalable materials design, transition metal chromophores made with 3d earth-abundant metals with d 6 electron configuration are preferred relative to their state-of-the-art 4d and 5d metal (e.g. Ru(II) and Ir(III)) analogs 7,8 .
Despite the success of this active learning approach in many applications, there remain significant challenges that prevent experimental realization of the predictions yielded by computational workflows. First, the DFT outcome depends on the density functional approximation (DFA) choice. A DFA that works well on certain systems may fail prominently on other systems due to the approximations made in the exchange-correlation functional 14,45,46 . When a single-DFA approach is used in VHTS, the DFA choice can lead to large biases in the data sets generated, which in turn biases the candidates the ML models recommend 47 . For transition metal complexes (TMCs) in particular, the electronic structure of a TMC is sometimes dominated by static correlation 48 that would make DFT error-prone, and predictions can be highly sensitive to DFA choice. Additionally, it is difficult to guarantee that the predicted lead molecules are synthesizable, despite the ability to add explicit constraints to ML models [49][50][51] . For TMCS the synthesizability problem becomes multiplicative 52 (i.e., all ligands comprising a TMC need to be synthesizable and compatible with complex formation).
In this work, we apply efficient global optimization 53 (EGO) to discover 3d 6 Fe(II)/Co(III) transition metal chromophores in a design space with 32.5M TMCs. We address the outstanding challenge of synthesizability of candidate chromophores by carefully crafting the design space with constraints using synthetically accessible fragments and ligand symmetries in Cambridge Structural Database (CSD). We avoid bias from DFA choice by applying a DFA-consensus approach 47 that considers property evaluation as an ensemble of predictions from 23 DFAs that span multiple rungs of "Jacob's ladder" 54 . Our active learning approach successfully identifies promising transition metal chromophores and is estimated to exhibit a 1,000-fold acceleration compared to the random sampling. We reveal that Co(III) complexes with large, strong-field ligands with more saturated bonds are preferred as candidate transition metal chromophores. By invoking Hammett tuning effects 55 and introducing functionalization on compounds, we further enriched the number of potential transition metal chromophores and verified our most promising candidates with time-dependent DFT (TDDFT) calculations.

Results and Discussion
Design space. We construct and explore a hypothetical design space of TMCs where all the constituent fragments (i.e., metal ions and ligands) are synthetically accessible (Figure 1). We further constrain the TMCs in the space to contain three bidentate ligands (e.g. Fe(II)(bpy)3) and restrict ourselves to d 6 Fe(II) or Co(III) metal centers, based on the precedent of this coordination type supporting efficient chromophores. We limit the number of unique ligands in a complex to two in order to promote the likelihood of synthesizability. We started with 5,173 CSD ligands that we curated in previous work 56 , including by assigning the charge states of the ligands. From this set, we selected bidentate ligands that contain common elements (i.e., H, B, C, N, O, F, Si, P, S, Cl, Br, and I) and ≤ 25 heavy (i.e., non-hydrogen) atoms, leaving a set of 812 ligands (Supporting Information). Combined with either Fe(II) or Co(III), the constraint of forming a complex with up to two unique bidentate ligands in an octahedral complex with three bidentate ligands produces 2´812=1,624 homoleptic and 2´812´811=1,317,064 heteroleptic TMCs. We refer to these 1.3M TMCs as the "base complexes". Hammett tuning is a commonly adopted strategy in experiments to fine tune the electronic properties of a complex by adding electron-donating or electronwithdrawing functional groups on conjugated rings. Here, we consider three distinct functionalization positions and ten functional groups, expanding the final design space to 32.5M functionalized TMCs, which we refer to as the "functionalized complex space" (Figure 1, see details in Exploring the functionalized design space).
Active learning procedure and design criteria. In transition metal chromophores, the photoexcited state should have a lifetime that is sufficiently long, such that the resulting chemical potential can be redirected before it is lost to unproductive competing pathways. Correspondingly, it is advantageous to have the photo-excited electron populate a long-lived metal-ligand charge transfer (MLCT) state and avoid low-lying metal-centered (MC) states that deactivate electron transfer from the expected photo-excited state. Therefore, we target complexes with low-spin (LS) ground states to increase the likelihood of MLCT states and to destabilize MC states 8 . We also desire target complexes to have weak multi-reference (MR) character. Avoiding high MR character has the benefit of avoiding complexes for which even a consensus DFT approach is likely to be inaccurate. It is possible to efficiently estimate MR character from fractional occupation number DFT as the contribution from nondynamical correlation 57,58 (i.e., rND 59 , see Methods).
When this value is low, we also anticipate a lack of deleterious low-lying electronic states. In addition, the absorption energy should fall within the wavelengths of the visible spectrum, ranging from 1.5 eV (825 nm) and 3.5 eV (350 nm). The absorption energy is estimated from D-SCF calculations 60 , which are more robust to DFA choice than the HOMO-LUMO gap from orbital energies (see Methods).
We use EGO with a 2D probability of improvement (P[I]) criterion to sample TMCs with LS ground states in a target zone of [1.5 eV, 3.5 eV] for D-SCF gap and [0, 0.307] for rND as candidate transition metal chromophores (Figure 1). The 2D P[I] is used to estimate the overall model probability (i.e., total area from the prediction and its model uncertainty) of residing within the target region and lacks the term present in some active learning algorithms that favor data points with higher uncertainty beyond a defined Pareto front. At odds with Pareto front-focused optimization, the 2D P[I] score employed here is amenable to the design goal of discovering a range of equally valid D-SCF gaps with modest rND values. The cutoff of 0.307 for rND was chosen based on our previous work 61 as a distinguishing cutoff for TMCs with weak vs. strong MR character. For EGO, we largely followed the established protocols from our previous work 43,44 : complexes in each new generation were selected by k-medoids sampling over the full design space.
We then used DFT to evaluate the D-SCF gap, rND, and ground spin state of these complexes. After combining the new data with data from previous generations, we retrained our ML models. Lastly, we used the updated ML models to evaluate the ground spin state and 2D P[I] of the complexes.
We selected the top 2,000 TMCs for k-medoids sampling to obtain the 200 complexes for DFT simulation in the next generation. Importantly, because both the ground state assignment and D-SCF gap are sensitive to DFA choice, we adopted a DFA consensus procedure 47 , where we considered an ensemble of 23 DFAs that cover the broad spectrum of "Jacob's ladder" of functionals to increase the robustness of our lead candidate chromophores (Supporting Information Figures S1-S3). Specifically, we only retained complexes when a majority of DFAs (i.e., 70% or > 16 of 23) predict the complex to have a LS ground state (Supporting Information Table S1). During the evaluation of 2D P[I], we consider 23 ML models separately trained on each DFA, from which the D-SCF gap and its corresponding model uncertainty (i.e., from a calibrated distance in latent space 62 ) are estimated. The rND is computed only from an ML model trained on a single DFA (i.e., B3LYP) because trends in rND values have been shown to be insensitive to the functional once calibrated 63 . On the contrary, the D-SCF gap is averaged over the 23 models due to its relatively high DFA sensitivity (Supporting Information Figure S3  All ligands in the CSD are first filtered to retain bidentate ligands with common elements and ≤ 25 heavy atoms, and known charge. The resulting 812 bidentate ligands are paired with either Fe(II) or Co(III), under the constraint that each complex has ≤ 2 unique ligand types, to form a design space of 1.3M base complexes. Lastly, these complexes are expanded to the full design space of 32.5M complexes with functionalization on the coordinating rings with a series of electron-donating or electron-withdrawing functional groups. (right) Active learning for discovering DFA-consensus-designed transition metal chromophores. DFT simulations are performed with 23 DFAs that span multiple rungs of "Jacob's ladder", which are used to iteratively train ML models. These ML models are applied to predict the ground spin state and D-SCF gap of complexes for 23 DFAs, and the MR character and their corresponding uncertainties. These quantities are used to select complexes with LS ground states and to evaluate the 2D P[I] of the design space to sample candidate complexes to compute in the next generation. The inset is an illustration of the ML prediction (solid dot), uncertainty (shaded area), and effective 2D P[I] area (solid outline) for multiple DFAs (blue and green) with respect to a target zone (rectangle with dashed lines).
Active learning on the 1.3M base complexes. We observe a strong negative linear correlation between D-SCF gap and rND for the 2,000 TMCs sampled in the initial generation (gen-0), which introduces difficulties for identifying candidates with simultaneously low D-SCF gap and rND ( Figure 2). This negative linear correlation exists because a small D-SCF gap generally suggests the existence of low-lying excited states, which would lead to high rND as MR character arises from near-degenerate occupied and virtual orbitals. In addition, we find that LS complexes often have stronger MR character, and thus higher rND relative to their HS counterparts because they can access more configuration state functions 61 (Supporting Information Figure S4). The nature of this multiple objective search for transition metal chromophores suggests that TMCs that can fulfill all our design requirements will be scarce. Indeed, for the 2,000 TMCs in gen-0, no complex with an LS ground state matches our target criteria with desirably low D-SCF gap and rND ( Figure 2). To put this in context, the lack of compounds in gen-0 suggests an extremely low probability, p, of a TMC residing in the target: when p is as low as 0.030%, there would still be 1/3 chance of finding at least one target complex in 2,000 random trials. Our ML models trained on gen-0 data also give a similarly conservative estimate that only 0.018% of TMCs have a 2D P[I] > 1/3, i.e., have a onethird chance of simultaneously fulfilling the two design criteria ( Figure 2).
Despite the initial absence of promising transition metal chromophores, we use active learning to discover lead TMCs in the target zone. The distributions of the sampled points in gen-1 to gen-3 shift towards the target zone, due to identification of compounds that overcome the trade-off of the negative linear correlation between D-SCF gap and rND present in gen-0 ( Figure   2). Although only 200 complexes are sampled at each subsequent generation, we discover numerous TMCs in the target zone once their DFT properties are explicitly computed. This enrichment is greatest in gen-2, where we identify 14 new TMCs that fulfill the design criteria, leading to a rather high (i.e., 7%) lead conversion rate (i.e., number of leads over number of samples). A conservative estimate using the binomial distribution shows that one would need to sample 200,000 TMCs randomly in the base complex space to produce 14 lead complexes, demonstrating a 1,000-fold acceleration of our active learning approach compared to random sampling. In addition, the ML models improve systematically as active learning proceeds from gen-0 to gen-3, as exemplified by reduction in relative MAEs of predicting the DFT results for the set-aside test data with increasing model generation (  The number of complexes in the target zone (black) and the percentage of the design space that has a 2D P[I] > 1/3 (brown) at each generation (top). The 2D P[I] at gen-0 is not available as the ML models have only been trained after gen-0. MAE for rND (middle) and D-SCF gap (bottom) from gen-0 to gen-3. At each generation, the ML models are trained on the combined training set of all previous generations and are tested on the set-aside test set of each generation separately. For example, the gen-2 model (blue bars) was trained on the combined training set of gen-0, gen-1, and gen-2 data. Generations are colored as follows throughout: green for gen-0, blue for gen-1, red for gen-2, and orange for gen-3. Gen-0 represents a k-medoids sampling of the 1.3M base TMC space.
After three generations of active learning, both the number of TMCs landing in the target zone (i.e., after computed by DFT) and the percentage of high 2D P[I] complexes decrease, indicating that most candidate base TMCs have likely been identified as well as the leveling off of ML model performance on the 1.3M base complexes. Therefore, we used the gen-3 models to screen through the base complex design space to reveal chemical trends for the 2,432 TMCs that have a reasonable probability of residing in the target zone (i.e., 2D P[I] > 1/6). Here, we use a smaller cutoff for 2D P[I] (i.e., 1/6 compared to 1/3) to retain a reasonable number of complexes for statistical analysis. From this set, we find that complexes with Co(III) and strong-field ligands (e.g., coordinating atom combinations of: CC, CN, NP, and PP) are significantly enriched ( Figure   3). This likely occurs due to our requirement of a LS ground state during the screening procedure ( Figure 1). Because we prefer a small D-SCF gap, the complexes that are favored by 2D P[I] tend to have large ligands, consistent with our previous observation that D-SCF gap has a negative linear correlation with complex size 61 ( Figure 3). Lastly, we find that complexes with reasonable (i.e., > 1/6) 2D P[I] tend to consist of ligands that are more saturated, as measured by their increased inverse average bond order 56 ( Figure 3). This trend can be understood by the fact that unsaturated ligands tend to contain higher MR character, and complex properties correlate (i.e., are additive) with those of their constituent ligands 56 (Supporting Information Figure S6). In general, we learn from our ML models that a complex with Co(III) and large, strong-field, and relatively saturated ligands would have an increased chance of being a transition metal chromophore with the desired properties. Exploring the functionalized design space. Hammett tuning, i.e., functionalization on conjugated rings, is a common procedure applied in experiments to fine tune the electronic properties of a TMC without dramatically scarifying its synthesizability. 64, 65 We considered three possible functionalizable positions, categorized by the bond depth (d) of the H atom on a ring with respect to the metal (Figure 4). For a six-membered ring, d=3, 4, or 5 corresponds to the ortho, meta, and para positions, respectively. For a five-membered ring, only d=3 or 4 is feasible. We consider a wide range of 10 electron-donating or electron-withdrawing functional groups ( Figure   4). To retain good likelihood of synthesizability, we constrain the in silico functionalization procedure to consist of one unique functionalizable position and one unique functional group for a TMC and disallow any combinations with multiple functionalizable positions or functional groups. Despite this constraint, the base design space will roughly be expanded by a factor of 25 after accounting for rings that are not functionalizable, leading to 32.5M TMCs (Figure 1).
However, the effect of Hammett tuning is not expected to be large enough to move all of the base complexes into the target zone, so those far enough from the target zone can be immediately discarded. For a representative Co(III) complex that resides in the target zone, we find that Hammett tuning can roughly tune the D-SCF gap by 1.0 eV and rND by 0.04 ( Figure 4). Therefore, we use the gen-3 ML models to screen through the 1.3M base complexes and only keep TMCs with a predicted D-SCF gap < 4.5 eV and rND < 0.35 as candidates for Hammett tuning. These  Since we have created a new design space with functionalized complexes that our ML models have not seen before, we expected that the 2D P[I] computed based on the trained model predictions and uncertainties would not be able to directly guide the exploration of candidate chromophores. Therefore, we repeated the k-medoids sampling that we performed in gen-0 but this time limited to the new 0.8M functionalized design space, selecting a set of 200 complexes.
Indeed, we find that the gen-3 models have significantly higher errors on the predictions of the kmedoids sampled gen-4 data ( Figure 5). However, the ML models improve quickly after retraining on the functionalized complexes in gen-4. Therefore, we expect the 2D P[I] to regain its predictive power and undertook two generations of active learning using 2D P[I] criteria. During these two generations, the ML models achieve MAEs that are comparable to those on the base complexes  Table S2). This observation showcases the effectiveness of our strategy for using Hammett tuning to further enrich a pool of candidate chromophores and improve their electronic properties.  To verify that complexes discovered through the active learning process have the desired excited state properties of transition metal chromophores, we computed the excited state properties with TDDFT of lead chromophores with the lowest rND for representative D-SCF gap values in the target zone from all seven generations (see Methods). Using this approach, we verified that six out of nine of our lead complexes have desired transition energy < 3.5 eV, where the photoexcitation process involves a MLCT state ( Figure 6). This observation is surprising because we do not explicitly set design objectives that involve explicit excited state calculations of our complexes during the active learning. Still, we achieve high likelihood (i.e., 67%) for obtaining lead complexes with promising excited state properties by carefully crafting the ground-state design objectives and using DFA consensus. The remaining three complexes have excited states with energy < 3.5 eV but are "dark" states (i.e., with small oscillator strengths), which cannot be foreseen by the ground state calculations that we performed during the active learning procedure (Supporting Information Table S3). If we had not required a consensus low D-SCF gap (i.e., < 3.5 eV) and LS ground state, the likelihood for obtaining leads with promising excited state properties would have been much lower, at around 15% (Supporting Information Table S4).
Since all the bidentate ligands considered are synthetically accessible, we propose these complexes as promising candidates for experimental verification. The quality of the leads selected by the consensus approach is expected to be better than from a standard single-functional screen because the functionals from different rungs agree, and previous work has shown DFT consensus more frequently produces leads similar to those that are experimentally validated. 47 Moreover, although complexes A-F are not in the CSD, they all contain ligands present in other synthesized compounds that have demonstrated photo-induced properties (Supporting Information Table S5).
For example, despite the fact that complex F has not been characterized experimentally, one of its constituent ligands (i.e., 4',5'-diaza-9'-[4,5-bis(methylthio)-1,3-dithiol-2-ylidene]fluorene) has been studied 66 for its interesting nonlinear optical properties in Co, Cu, and Cd complexes ( Figure   6). These observations showcase the power of our strategy for identifying experimentally relevant candidate transition metal chromophores that were potentially missed by previous experimental exploration due to the combinatorial challenges in chemical discovery.

Conclusions
We applied EGO with a 2D P[I] criteria to discover potential 3d 6 Fe(II)/Co(III) transition metal chromophores in a design space of 32.5M compounds that simultaneously fulfill three design objectives: LS ground state, a D-SCF gap corresponding to an electronic transition in the visible region of the electromagnetic spectrum, and weak MR character (i.e., low rND). We avoid common biases that arise from DFA choice in VHTS and ML-accelerated chemical discovery by applying a DFA-consensus approach that considers the property evaluations from 23 DFAs that span multiple rungs of "Jacob's ladder". We also addressed the challenge of synthesizability for computationally designed functional molecules by constraining the design space construction to ligands that are synthetically accessible and symmetry classes that are easy to access in experiments. Compounds discovered through this active learning workflow therefore have a higher likelihood to be synthesizable and higher fidelity (i.e., more robust to changes in DFA choice).
Despite the scarcity of potential transition metal chromophores in our design space, judged by the fact that no compounds in the initial 2,000 samples land in our target objective zone, our active learning process gradually shifts the distributions of sampled compounds towards the target zone and successfully identifies many leads. A conservative estimate suggests that our active learning approach achieves a 1,000-fold acceleration relative to the random sampling.
Interrogation of our ML models revealed that Co(III) complexes with large, strong-field, and relatively saturated ligands are preferred as candidate transition metal chromophores. To fine-tune their electronic properties, we further introduced Hammett tuning by functionalization of the base complexes, which further increased the number of complexes that satisfied the design criteria.
Lastly, we performed TDDFT calculations on the nine most promising leads. We found that six of the nine compounds demonstrated desired excited state properties with MLCT states and contain ligands that have been previously studied experimentally due to their interesting optical properties.
We expect our strategy for design space construction and DFA-consensus enhanced active learning workflow to be broadly useful in discovering candidate molecules and materials that are more synthesizable and computationally robust in transition metal chemical space.

DFT calculation details.
All initial geometries were generated using molSimplify 67,68 , where initial ligand geometries were derived from the crystal structures of transition metal complexes containing the ligands (Supporting Information). DFT geometry optimizations were carried out using TeraChem 69 , as automated by molSimplify 67,68 with a 24 h wall time per run with up to five resubmissions. These calculations used the B3LYP 70-72 global hybrid functional with the LACVP* basis set, which corresponds to the LANL2DZ 73 effective core potential for transition metals (i.e., Fe, Co) and heavier elements (i.e., I or Br) and the 6-31G* basis for all remaining elements. These geometries were optimized using the L-BFGS algorithm in translation rotation internal coordinates (TRIC) 74 to the default tolerances of 4.5 ´ 10 -4 hartree/bohr for the maximum gradient and 10 -6 hartree for the energy change between steps. All HS (i.e., quintet) states were calculated with an unrestricted formalism and LS (i.e., singlet) states with a restricted formalism. In all calculations, level shifting of 0.25 Ha was employed between the occupied and virtual spin orbitals. Geometry checks were applied to eliminate optimized structures that deviated from the expected octahedral shape following previously established metrics without modification. 75 Open-shell structures were also removed from the data set following established protocols if the expectation value of the S 2 operator deviated from its expected value 75 of S(S + 1) by >1 µ 2 B (Supporting Information Table   S6).
For optimized TMCs, we followed our established protocol 47  DFAs as in our previous work 47 that were chosen to be evenly distributed among the rungs of "Jacob's ladder" 77 (Supporting Information Table S7).
We evaluated the rND diagnostic 57,58,78 by performing finite-temperature DFT 79 calculations using TeraChem 69 . Specifically, we followed a literature recommendation 57,58,78 to use a temperature of 9,000 K for B3LYP. Here, we evaluated fractional occupation numbers (FON) from a broadened distribution (i.e., with Fermi-Dirac statistics).
We performed linear response TDDFT calculations with the Tamm-Dancoff approximation using wB97X-D/def2-TZVP following by the recommendation of a recent benchmark study 80 in Psi4 1.4 76 . We used a polarizable continuum implicit solvent model with water (e = 80) as the solvent. Because we focus most on the lowest few excitations, only the first 30 states were computed. We broadened simulated spectra using Lorentzian functions, and we considered only excited states with significant oscillator strength (i.e., fosc > 0.01 a.u.).

ML models.
As in our prior work, we use extended revised autocorrelations 81,82 (eRACs) as descriptors for all our machine learning models. The eRAC features are sums of products and differences of six atom-wise heuristic properties (i.e., topology, identity, electronegativity, covalent radius, nuclear charge, and group number in the periodic table) on the 2D molecular graph.
As motivated previously on large TMCs 43 , we applied the maximum bond depth of four and eliminated RACs that were invariant over the mononuclear octahedral transition metal complexes.
We used metal oxidation state and total ligand charge of a complex as two additional features.
Because we would like to discover transition metal chromophores with a LS ground state with certain ranges of Δ-SCF gap and rND, we built ML models to predict these three properties.
Specifically, we built i) a classification model to predict whether a complex fulfills the consensus LS condition (i.e., > 16 DFAs categorize the ground state to be LS), ii) a regression model to predict rND of a LS complex, and iii) 23 separate models to predict the Δ-SCF gap of a LS complex for each DFA (Supporting Information Table S1). In our workflow, we first used the ground state classification model to filter out complexes that do not satisfy the consensus LS condition. We used the energy from both the HS and LS optimization of a complex as training data for the model to determine its ground spin state. On the contrary, only the LS calculation was used for building the ML models that predict Δ-SCF and rND. For the 23 Δ-SCF gap models, we adopted our established workflow 47 to fine-tune the 22 non-B3LYP models initialized by the weights of the B3LYP model to avoid randomness in the weight initialization and to increase the consistency between ANN models trained with DFT data derived from different DFAs.
During each generation of the active learning, we partitioned the data using a random 80%/20% train/test split and used 20% of the training data (i.e., 16% overall) as the validation set.
As in our prior work 43 , all ANN models were trained using Keras 83 with a Tensorflow 84 backend and Hyperopt 85 for hyperparameter selection for gen-0 data (Supporting Information Table S8).
For all other generations, the models were only fine-tuned with a reduced learning rate (i.e., 10 -5 ) on the combined training set of all previous generations. All ANN models were trained with the Adam optimizer up to 2,000 epochs, and dropout, batch normalization, and early stopping were applied to avoid overfitting.

Supporting Information Statement
Histogram for average ground state spin; D-SCF gap computed at different DFAs and spin states; Page S3 Table S1 Ground state labeling with DFA consensus Page S3 Figure S4 rND of optimized structures at LS and HS state Page S3 Figure S5 Model AUC for ground spin state classification Page S4 Figure S6 Ligands' rND of the 2,432 selected and 1.3M base complexes Page S4           Table S5. Summary of ligands in complexes A-F that have parent complexes in CSD that demonstrate interesting photo-induced properties and the associated references. The complex ID corresponds to the index of the six complexes in Figure 6 of the main text. The ligand ID is our internal database ID, for which the corresponding geometry can be found in the Supporting Information Data. For each ligand, the refcode for the CSD complex that contains this ligand and demonstrates interesting photo-induced properties is also shown.  Table S6. Summary of the filtering statistics for the TMCs at each generation, including the number attempted, number of complexes with valid initial geometry (i.e., no ligand clashing), and the number of complexes with good geometry and <S 2 > (judged by our established protocols 6-7 ) after converged geometry optimizations. Generations that explored the base complexes space are shaded in blue and those that explored the functionalized complexes are shaded in green.