Low-cost machine learning approach to the prediction of transition metal phosphor excited state properties

Photoactive iridium complexes are of broad interest due to their applications ranging from lighting to photocatalysis. However, the excited state property prediction of these complexes challenges ab initio methods such as time-dependent density functional theory (TDDFT) both from an accuracy and a computational cost perspective, complicating high throughput virtual screening (HTVS). We instead leverage low-cost machine learning (ML) models to predict the excited state properties of photoactive iridium complexes. We use experimental data of 1,380 iridium complexes to train and evaluate the ML models and identify the best-performing and most transferable models to be those trained on electronic structure features from low-cost density functional theory tight binding calculations. Using these models, we predict the three excited state properties considered, mean emission energy of phosphorescence, excited state lifetime, and emission spectral integral, with accuracy competitive with or superseding TDDFT. We conduct feature importance analysis to identify which iridium complex attributes govern excited state properties and we validate these trends with explicit examples. As a demonstration of how our ML models can be used for HTVS and the acceleration of chemical discovery, we curate a set of novel hypothetical iridium complexes and identify promising ligands for the design of new phosphors.


Introduction.
The interactions between light and matter underpin phenomena ranging from photovoltaics 1 to photosynthesis 2 to bioluminescence. 3 The design of new functional materials that can leverage light-matter interactions has led to significant technological advancements. [4][5][6] Exemplary of these advancements are iridium photoactive complexes that have been investigated extensively due to their applications in lighting and display technology 7-10 , photocatalysis 11-13 , and bioimaging. 14,15 The spin-orbit coupling (SOC) produced by iridium causes these complexes to efficiently convert excitons into light or chemical energy. 16 Simultaneously, iridium uniquely limits nonradiative decay rates by destabilizing a metal-centered ( 3 MC) triplet excited state due to strong metal-ligand bonding 17 , further improving complex efficiency. In these iridium-centered complexes, the judicious selection of ligands allows for the modulation of phosphorescence color (i.e., emission energy) and efficiency/brightness (i.e., via modulating excited state lifetime and photoluminescence quantum yield). The accurate prediction of these three excited state properties is important for the discovery of novel iridium complexes to enable vibrant display technologies and green photocatalysis.
To screen a large number of compounds, computational modeling with time-dependent density functional theory (TDDFT) can enable rapid transition metal complex property prediction.
While these methods are commonly employed to estimate emission energies 18-26 , the calculation of lifetime and quantum yield is more challenging both from an accuracy and a computational cost perspective. The calculation of lifetime 21,27-31 requires the inclusion of SOC in TDDFT to estimate the transition dipole moment between the T1 sublevels and S0. The calculation of photoluminescence quantum yield further requires the calculation of nonradiative rates, which entails the use of the thermal vibration correlation function rate theory [32][33][34] and excited state geometry optimization. 34,35 Thus, while ab initio computational methods have provided valuable insight into the properties of iridium complexes, they are computation-intensive, requiring around one day of computation time per complex for the least demanding calculations, and may not reach the accuracy required to enable rational design.
Supervised machine learning (ML) has emerged as a powerful complement to ab initio methods in recent years due to its capacity to reproduce ab initio results at significantly lower cost [36][37][38][39][40] , enabling screening vast regions of chemical space. 41 Furthermore, ML models can be trained on experimental data, enabling the prediction of properties that challenge ab initio methods, such as material stability. 42 With regard to excited state properties, ML models have been successfully applied for the prediction of phosphorescence energies 43 , fluorescence rates 44 , and fluorescence energies and quantum yields 45 after training on ab initio or experimental data.
However, to our knowledge ML has yet to be applied to predict the excited state properties of transition metal (e.g., iridium) photoactive complexes, which are workhorses of photocatalysis and OLED technology. Thus, the need to predict for and simultaneously optimize multiple properties challenging to predict with TDDFT in iridium complexes motivates the application of ML.
In this work, we use ML to predict three key properties of iridium complexes: Em50/50 (mean emission energy), excited state lifetime, and emission spectral integral (brightness). We train and evaluate artificial neural networks (ANNs) on a recent experimental dataset 46 of 1,380 iridium(III) phosphors and their properties. We show that features generated with density functional tight binding lead to the most predictive ANN performance and generalization on outof-sample data. Using these features, we identify trends in phosphor properties, and we extend our models to a new set of hypothetical iridium phosphors. These experimentally-informed ANNs enable fast, accurate prediction of iridium phosphor properties for the rapid exploration of chemical space.

2a. Dataset
We constructed the structures of bidentate ligands used in the experimental study of DiLuzio et al. 46 on Ir(III) complexes of the form [Ir(CN)2(NN)] + (Figure 1). We assigned each ligand as either cyclometalating (CN) or ancillary (NN), as determined by the two iridiumcoordinating atom identities. We studied the same 60 CN ligands and 23 NN ligands from the prior experimental study 46 , excluding only the monodentate DMSO ligand in the prior work, giving rise to a combinatorial set of 1,380 [Ir(CN)2(NN)] + phosphor complexes. This set of ligands will be referred to as the high-throughput ligand set (HLS), and we use the same labeling as in the prior study when referring to individual ligands (ESI Tables S1 and S2). We used experimental data from the prior study 46 on the three target properties, Em50/50, excited state lifetime, and emission spectral integral. The experimental values for these properties were reported for each of the 1,380 iridium phosphors, and were used for ML model training and performance assessment (ESI Figure   S1). CN ligands were generated in their neutral form (i.e., with a proton added) when featurized.
Because of this, all ligands are neutral with the exception of three NN ligands (ESI Text S1). After ligand construction using the draw tool in Avogadro v1.1.2 47,48 and force field (i.e., UFF) optimization, we used these ligands to generate the structures of all possible iridium complexes with one distinct type of CN ligand and NN ligand using molSimplify v1.6.0 49,50 and force field optimized again. Atoms are colored as follows: white for hydrogen, gray for carbon, blue for nitrogen, red for oxygen, light blue for fluorine, and green for chlorine.

2b. Feature Sets
We evaluated eight representations of the iridium complexes to identify the most suitable set of features for training ML models to predict Ir phosphor properties. We compared eight feature sets: those based on similarity (i.e., Morgan, Dice), those based on graph descriptors (i.e., ligandonly revised autocorrelation, RACs 51 , whole-complex RACs, and Coulomb-decay RACs 52 , CD-RACs), and those based on electronic structure calculations (i.e., xTB, wPBEh, and B3LYP).
For the similarity feature sets, we generated Morgan fingerprints 53,54 , which have been used previously in machine learning chemistry applications [55][56][57][58][59][60] , by one-hot encoding of groups of atoms in a structure. We computed these with a radius of three and 2,048 bits on the isolated CN and NN ligands to capture the presence and absence of chemical substructures. We also generated Dice similarity coefficients 53 of ligand Morgan fingerprints. In this approach, we compare the CN and NN ligands of each new iridium complex to all HLS CN or NN ligands through the Dice similarity metric (ESI Text S2 and Table S3). Dice similarity was selected after determining it outperformed the commonly employed Tanimoto similarity (ESI Table S4).
For the graph-based descriptor sets, we generated RACs 61,62 for both the isolated ligands and the full iridium complex structures (ESI Text S2). RACs are connectivity-based representations that have shown good performance for transition metal complex (TMC) property prediction. 41,61,63 For RACs, a TMC is represented as a molecular graph, with vertices for atoms and unweighted (i.e., no bond length or order information) edges for bonds. Each RAC feature is the sum of products or the sum of differences of heuristic atomic properties at depth d on a TMC molecular graph, where d indicates the number of edges separating the starting and ending atoms (ESI Text S2). The RACs include those that span over the entire complex as well as weighted averages over the equatorial ligands and axial ligands. We use the largest set of heuristic properties described in previous work, including both group number 64 and number of bonds 52 , leading to a final RAC feature set that contains 196 features (ESI Table S5). For the ligand-only RAC feature set, we generated full-scope product RACs (i.e., all atoms are used as starting atoms) on isolated CN and NN ligands for each TMC. We concatenate individual feature vectors for the CN ligands and NN ligands (ESI Text S2). The ligand-only RAC feature set contains significantly fewer (i.e., only 70) features than the RAC feature set (ESI Table S6). We also generated Coulomb-decay RACs 65 on the iridium complex structures that were optimized with UFF (see Dataset). CD-RACs are a variant of RACs that also encode internuclear distance information between the atoms in the RAC feature (ESI Text S2). The CD-RAC feature set contains Coulomb-decay versions of the features in the RAC feature set but is of higher dimension (i.e., 222 features) due to the added information from the geometry (ESI Table S7).
Finally, we computed descriptors obtained from electronic structure theory. Specifically, we selected electronic properties of the isolated ligands due to the lower computational cost in comparison to whole-complex properties. These ligand-based decriptors include the HOMO and LUMO of each ligand type, the ionization potential (IP) and electron affinity (EA) of each ligand type, and the partial charges (i.e., Mulliken) of each of the metal-coordinating atoms (ESI Table   S8). For the xTB feature set, we utilized a specially reparametrized 66,67 version of GFN1-xTB, a low-cost, semi-empirical tight binding method that has parameters for most elements in the periodic table. 68 We calculated ligand-only xTB features on UFF-optimized CN and NN ligands.
We confirmed the suitability of ligand-only feature calculation due to the lower computational cost. The xTB features consist of contain electronic structure information on the CN and NN ligands of a phosphor, and are in some cases correlated (ESI Table S8 and Figure S2). For the B3LYP and wPBEh DFT feature sets, we performed density functional theory (DFT) calculations on isolated CN and NN ligands using the B3LYP 69-71 or wPBEh 72 exchange correlation functionals respectively (see Methods). Mulliken charges were used after they were found to outperform natural charges on the grouped split test set (ESI Table S9).

3a. Features for ANN Models and Performance.
The representation for the phosphor is a crucial piece in determining which machine learning (i.e., ANN) models are most likely to both predict experimental properties accurately and generalize to unseen complexes. A model and feature set that performs well on one property may perform poorly on another. We thus trained models with each of the eight feature sets and the three target properties (Em50/50, excited state lifetime, and emission spectral integral) and assessed their prediction performance on both a random split and grouped split. For the first random split, the Dice, Morgan, and xTB feature sets lead to the lowest model test set errors across all three target properties. Conversely, the B3LYP DFT and RAC feature sets lead to the largest model test set errors, while the remaining feature sets CD-RAC, ligand-only RAC, and wPBEh DFT exhibit intermediate performance (Figure 2 and ESI Figures S3 and S4 and Tables S10-S13). The composition-based Dice and Morgan feature sets lead to the best predictions, judged on the basis of scaled MAEs of 0.03 to 0.05 for the three target properties (ESI Tables S11-S13). There is a substantial difference in performance: the percent difference in random split mean absolute error (MAE) between using the optimal feature set and the worst feature set for mean emission energy is 80% and a similar performance erosion is observed for the other two properties (i.e., 59% for spectral integral and 40% for phosphorescence lifetime). We rationalize the relative performance of each feature set in the random split assessment by considering what aspects of ligands and complexes the different feature sets capture. We attribute the predictive power of the Dice feature set to the fact that phosphorescent properties are very ligand-directed, and in the random split, each ligand is represented in both the train and test sets. We attribute the predictive power of the Morgan feature set in the random split to the identification of substructures in ligands that affect phosphor excited state properties by tuning energy levels and ligand rigidity. Thus, features encoding ligand similarity to previously observed ligands and substructures present in ligands lead to the best performance on the random split. On the other end of the spectrum, the random split ANNs trained on the RAC feature set (scaled MAE: 0.05 to 0.08) likely perform poorest because they encode significant metal-local information that does not vary across this set, since all complexes have an iridium center and an identical first coordination shell (ESI Tables S5 and S11-S13).
With regard to the ligand-only electronic structure feature sets, we surprisingly observe improved random split ANN performance with the xTB feature set relative to the two DFT feature sets. While xTB is expected to be faster than DFT for feature generation, electronic structure properties from xTB alone should not necessarily be more accurate. The Em 50/50 xTB model error is 30% lower than that of the corresponding B3LYP DFT model (i.e., 0.021 eV vs 0.029 eV).
Nevertheless, most (i.e., eight of twelve) xTB features have high (>0.5) linear correlation with their B3LYP DFT and wPBEh DFT counterparts, as determined by Pearson correlation coefficients (ESI Table S14). Thus, the GFN1-xTB method provides reliable electronic structure information from which our ANNs can generate accurate predictions, even if this electronic structure information is not of the same accuracy as that generated by DFT.
To assess the utility of these models for discovery of out-of-distribution complexes, we   next used the feature sets to train new ANNs on a grouped split, where five ligands are present   only in the test set and are consequently unseen by the ANNs during training (see Methods). We used the same grouped split across each property prediction task. We find that the test accuracy of the grouped split ANNs is worse than that of the corresponding random split ANNs for all features, due to the presence of unseen ligands in the grouped split test sets (ESI Tables S15-S17). However, this worsened performance is more significant for the spectral integral and Em 50/50 target properties.
Overall, the MAE averaged over all feature sets is significantly worse for these two properties (157% or 164% worse on average for spectral integral and Em 50/50 , respectively) than for phosphorescence lifetime (28% worse, ESI Table S18).
With the grouped split, the predictive power of the xTB feature set is further improved relative to the other feature sets (Figure 2 and ESI Figures S3 and S4). For Em 50/50 prediction, the xTB feature set improves from the third best feature set to the best feature set as a result of its scaled MAE increasing less than the random split best performers (i.e., 0.04 to 0.078 for xTB versus Dice 0.031 to 0.138, ESI Tables S18 and S19). In practice, this means that the Em 50/50 xTB MAE doubles from 0.021 eV to 0.041 eV, while the Dice MAE nearly quadruples from 0.016 eV to 0.072 eV. We attribute this particularly worsened performance of the Dice feature set to the reduction of the total number of features in the Dice feature set in the grouped split (ESI Table   S3). The poor generalizability of the Dice feature set can also be attributed to the pseudo one-hot encoding of ligands via the similarity scores. The xTB features, in contrast, convey physical information that extrapolates beyond the ligands seen in the training data. We ultimately chose the xTB feature set for further analysis in evaluating hypothetical complexes, because the xTB feature set has favorable performance on the grouped split for all three properties, indicating that the ANNs using the xTB feature set generalize well.
To quantify ANN uncertainty in predictions for new phosphors outside of our initial training set, we use the latent space distance as a measure of how similar a new phosphor is to the phosphors used to train the model. 73 To confirm that this is a good measure of similarity that quantifies uncertainty for the current prediction task, we assessed the influence of latent space distance on test set prediction accuracy of the random split ANNs that use xTB inputs. Following prior work 73 , we computed the average distance to ten nearest neighbors in the latent space formed  The uncertainty quantification (UQ) cutoff versus test set mean absolute error (in eV) of the random split ANN model trained on the xTB feature set and predicting Em50/50. The data fraction is the number of test set complexes under the corresponding UQ cutoff, and the MAE is calculated on this subset of complexes. The UQ metric used is the average latent space distance to the ten nearest neighbors in the training set following the protocol introduced in Ref. 74 . The UQ metric is normalized such that the largest UQ metric is scaled to 1.

3b. Feature Importances and Trends
Given the high accuracy of the xTB-trained ANN models, we next sought to determine if simpler and more interpretable linear and random forest models trained on xTB features could attain similar accuracy. These models allow us to more transparently gain insight into which features most heavily influence phosphor property prediction. We trained random forest regression models that use xTB features to predict each of the three target properties using scikit-learn. 75 Using the same random split train/test partition, these random forest models have comparable performance to the ANNs, indicating that they can be analyzed to determine which features inform phosphor properties (ESI Table S20). Random forest models significantly outperform linear ridge regression models, which preclude their use for feature analysis but could have been anticipated from the low linear correlation of individual xTB features to target properties (ESI Table S20 and Figure S7). From random forest models, we extracted impurity-based feature importances to interpret the ML model predictions (Figure 4). We find that xTB features of the CN ligand are more important than those of the NN ligand in predicting Em 50  Focusing more on Em 50/50 , we find that ionization potential (IP) and electron affinity (EA) are important for model predictions, as are the charges of metal-coordinating atoms ( Figure 4).
Specifically, the top three xTB features for predicting Em 50/50 are the IP of the CN ligand and two of the nitrogen-coordinating charges (i.e., N charge (CN) and N1 charge (NN)). This suggests that the information provided by the charge features allows the Em 50/50 ML model to make better predictions. The information provided by IP (CN) for the cyclometalating ligand also is important for model performance, and this conforms to prior observations that ligand energy levels affect emission energy. 76,77 We also emphasize that these three xTB features vary significantly over the experimental dataset. The IP (CN) varies by nearly 1.5 eV (i.e., from 7.56 eV to 9.03 eV), and the partial charges have a 0.1 a.u. range (i.e., N charge (CN) from -0.35 a.u. to -0.24 a.u. and N1 charge (NN) from -0.37 a.u. to -0.28 a.u) ( Figure 5 and ESI Figure S8). Thus, tuning these three features in a coordinated fashion should enable tuning of Ir phosphor complex Em 50/50 .

3c. New Compound Exploration
We applied the three random split ANNs that use xTB features as inputs to screen hypothetical iridium complexes generated from CSD ligands (see Methods and ESI Text S1). We For this set of curated hypothetical complexes, we evaluated which ligands are present in the complexes with the highest and lowest ANN-predicted properties (ESI Figure S11). We find that specific ancillary ligands tend to be well-represented in complexes with extreme properties, indicating that phosphor properties are tuned by these ancillary ligands (ESI Table S21). For example, the ligand that appears most often in hypothetical complexes with high predicted lifetime is the ancillary ligand from the CSD structure with refcode RASGAV. This conjugated ligand has a relatively low IP (NN) of 7.79 eV, which contributes to a longer lifetime following the previously identified trend ( Table S21). We also note clear xTB feature trends in predictions for spectral integral and Em 50/50 . The low IP (NN) of the ancillary ligand from RASGAV leads to a hypothetical complex with one of the highest predicted spectral integrals ( Figure 4 and 8 and ESI Table S22). With regard to Em 50/50 , the fluorinated cyclometalating ligand from the CSD structure with refcode RADTEZ has a high ionization potential (9.24 eV). The high IP (CN) feature appears to contribute to a high emission energy, as the RADTEZ CN ligand is present in the three highest predicted Em 50/50 hypothetical complexes ( Figure 4 and Figure 8 and ESI hypothetical complexes ( Figure 4 and Figure 8 and ESI Table S22). Thus, we find that many ligands that lead to extreme hypothetical phosphor predicted properties follow our identified xTB feature trends from the experimental data. This lends interpretability to our model predictions and indicates that these predictions are derived from the electronic structure properties of the ligands.

Figure 8.
Ligands mined from the CSD that lead to very high or very low phosphor properties predicted by the ANNs along with their percentile rank of the relevant property in the context of the experimental complexes. Coordinated nitrogen and carbon atoms are indicated with blue and gray circles respectively. Atoms are colored as follows: white for hydrogen, gray for carbon, blue for nitrogen, red for oxygen, light blue for fluorine, and yellow for sulfur.
To further validate random split ANN predictions, we obtained TDDFT excited state energy and lifetime predictions and compared them to ANN predictions over complexes in both the experimental dataset and the hypothetical dataset. We found that using singlet geometries instead of triplet geometries as inputs to TDDFT leads to better agreement with experiment, although geometries do not differ greatly from a RMSD perspective (ESI Table S23 and Figure   S12). Emission energies calculated with B3LYP were found to correlate with experiment better than those calculated with the range-separated hybrid functional CAM-B3LYP, motivating our use of B3LYP for TDDFT (ESI Figure S13).
Over a group of 26 representative test set complexes from the experimental dataset, we find that TDDFT overestimates experimental emission energy by 0.3 eV on average, and further find that TDDFT predictions correlate with experiment less well than the Em 50/50 ANN predictions ( Figure 9 and ESI Tables S24 and S25). These results show that the Em 50/50 ANN achieves excellent performance. Even after applying a rigid downward shift to TDDFT energy predictions, they exhibit a larger spread around the experimental values than the predictions of our Em 50/50 ANN.
Over the same 26 complexes, TDDFT lifetime predictions trend with experiment and ANN predictions; however, unlike the case of Em 50/50 , TDDFT predictions outperform our lifetime ANN for complexes with long lifetimes (ESI Figure S14). This shortcoming of the lifetime ANN can be rationalized by the lower amount of phosphors with long lifetimes in the experimental dataset used for model training (ESI Figure S1 and  In the insets, atoms are colored as follows: white for hydrogen, gray for carbon, blue for nitrogen, and dark blue for iridium. The dotted line is included as a reference and corresponds to perfect agreement between prediction and experiment.

Conclusions.
While ab initio methods like TDDFT are useful tools for studying iridium phosphor excited states, they are both computation-intensive and can also have eroded accuracy, motivating the use of machine learning. Using experimental data on 1,380 iridium phosphors, we trained ANNs to predict three experimental properties: Em50/50, excited state lifetime, and emission spectral integral.
We found that features calculated with xTB led to the best overall performance across the three properties on out-of-sample complexes, outperforming the standard Morgan fingerprint features and features based on RACs. We then used random forest regression models to determine which xTB features most influence phosphor properties, deriving chemical insight. We found that high cyclometalating ligand ionization potential is indicative of high Em50/50, while high ancillary ligand ionization potential correlates to low lifetime and low spectral integral. These observations illustrate how phosphor properties can be altered through judicious ligand selection.
We next demonstrated how our ANNs can be applied to chemical exploration by considering hypothetical iridium phosphors derived from ligands found in the CSD. We identified cyclometalating and ancillary ligands that lead to edge-of-distribution properties, such as an ancillary ligand predicted to result in both high lifetime and spectral integral phosphors. We confirmed the validity of these predictions by comparing to TDDFT, showing that for Em50/50 the ANN significantly outperforms TDDFT, while for lifetime the corresponding ANN performs well only in regimes of sufficient training data. The ANN models for iridium phosphor property prediction that we present here are promising tools for chemical screening and the acceleration of chemical discovery, as they can be used to quickly evaluate thousands of hypothetical iridium phosphors to identify promising candidates for follow-up synthesis.

5a. Feature and Structure Generation
We generated feature sets to represent the 1,380 Ir phosphor complexes as inputs to ML models (see Feature Sets). We generated all ligands using the draw tool in Avogadro v1.1.2 47,48 and subsequently optimized them with UFF. 78 We used molSimplify v1.6.0 for the generation of RAC feature sets on either ligands or complexes. 49,50 The xTB features were generated using xTB 6.4.0 68 , and DFT-based features were generated using the B3LYP [69][70][71]  hartree/bohr for the maximum gradient and 1 × 10 -6 hartree for the change in energy between steps.
We then performed single point energy calculations on the optimized neutral ligand geometries at two different charges: +1 and -1 (ESI Text S1). We used the LACVP* basis set, which for the HLS ligands corresponds to the LANL2DZ 83 effective core potential for Br and the 6-31G* basis set for all remaining elements. We calculated all non-singlet states with an unrestricted formalism and singlet states with a restricted formalism. Level shifting of 0.25 Ha was employed on both virtual and occupied orbitals to facilitate self-consistent field convergence. We used the hybrid DIIS 84 /A-DIIS 85 scheme for the self-consistent field procedure. We used TeraChem dynamic precision and a grid with approximately 3,000 points per atom. Like the xTB feature set, DFT-generated features encode electronic structure information (ESI Table S8).

5b. ML Models
We trained multiple artificial neural networks (ANNs) with for each of the seven feature sets to predict three target properties: Em50/50, excited state lifetime, and emission spectral integral.
For all random split ANNs, we used a random 70%/15%/15% train/validation/test split of the 1,380 complexes from the prior study 46 and refer to this hereafter as the random split. To assess the generalizability of our best-performing ANNs, we also carried out grouped splits where we excluded from the training and validation data any complex containing a ligand from a select

5c. Out-of-distribution Complexes
We identified hypothetical out-of-distribution iridium complexes which we enumerated combinatorially using CN and NN ligands not in the HLS. We selected these ligands by screening the CSD v5.42 + 2 updates, released in November 2020, for iridium complexes with two CN ligands and one NN ligand by specifying the first coordination sphere around iridium in a ConQuest 2021.1.0 search. Complexes selected by the screening were then examined by hand, and those that were not fit for analysis were eliminated (ESI Table S31). The molSimplify code was used to identify unique ligands from the remaining complexes on the basis of their atom-weighted molecular graph determinants 95 , and we used any ligands not already in the HLS in combination with the HLS to generate new hypothetical [Ir(CN)2(NN)] + complexes (ESI Text S1).

5d. TDDFT Calculations
For ab initio validation of predictions using TDDFT, iridium phosphors were first geometry optimized, and TDDFT was then run on the optimized geometries using the ORCA Due to relativistic SOC caused by iridium, the T1 manifold is split into three sublevels (zero-field splitting). For the calculation of ab initio energy, the energies of these three lowest triplet sublevels from the TDDFT calculation were averaged for each complex. For ab initio lifetime, radiative rate and radiative lifetime were calculated as in prior work 27-29,31,103 using output from TDDFT calculations. The radiative rate ! from a triplet sublevel is given by: where ! is the radiative lifetime of sublevel , " = ∆ ! is the excitation energy in atomic units from the ground state to the sublevel , and $ ! is the -axis projection of the transition dipole moment in atomic units between the ground state and the sublevel .
The overall radiative lifetime from the three triplet sublevels is calculated as a Boltzmann average of radiative rates that depends on the energy differences between triplet sublevels.
∆ =,, is the energy difference between sublevels 1 and 2, and ∆ =,# is the energy difference between sublevels 1 and 3. T = 300K was used. This equation for lifetime does not take into account nonradiative decay, which can be significant in some cases. In order to account for the DMSO solvent, the ab initio lifetime was divided by the square of the refractive index of DMSO

DATA AVAILABILITY
The datasets supporting this article have been uploaded as part of the ESI.

CONFLICTS OF INTEREST
The authors declare no competing financial interest.   Figure S1 Histograms of the three target properties Page S7 Figure S2 Correlation between xTB features Page S8 Figure S3 Feature set performance on lifetime Page S9 Figure S4 Feature set performance on spectral integral Page S9 Figure S5 Uncertainty quantification for the lifetime xTB ANN Page S10 Figure S6 Uncertainty quantification for the spectral integral xTB ANN Page S11 Figure S7 Correlation between xTB features and target properties Page S12 Figure S8 Histograms of xTB features Page S13 Figure S9 Ligand substitution effect on lifetime Page S14 Figure S10 Ligand substitution effect on spectral integral Page S14 Figure S11 Sixteen CSD ligands that lead to extreme complexes Page S15 Figure S12 Phosphor geometry differences for singlet vs triplet state Page S17 Figure S13 TDDFT-calculated energy with CAM-B3LYP Page S18 Figure S14 Calculated lifetime vs experiment Page S18 Figure S15 ANN Em50/50 vs TDDFT energy Page S19 Figure S16 ANN lifetime vs TDDFT lifetime Page S19 Figure S17 Confusion matrices for 2 μs lifetime cutoff Page S20  Text S1. Additional details of calculations and CSD search.

Additional DFT calculation details
As mentioned in the main text, we conducted single point energy calculations on the optimized neutral ligand geometries at two different charges: +1 and -1. In combination with information from the neutral geometry optimization, this allowed us to calculate vertical IP and EA values analogous to those generated by GFN1-xTB. 1,2 The remaining information (HOMO, LUMO, and Mulliken charges of coordinating atoms) was extracted solely from the neutral geometry optimization calculation for a given ligand. We specified ligands to be in a singlet state when neutral and in a doublet state otherwise.

Description of NN40, NN41, and NN42
These three ancillary ligands contain tetrazole or pyrazole moieties that deprotonate upon metal coordination. 3 Consequently, we generated these ligands without a hydrogen on the pyrazole/tetrazole nitrogen coordinating atom. As a result, while all other HLS ligands are neutral, these ligands have a charge of -1. We specified this charge for GFN1-xTB calculations, and also adjusted our DFT workflow to account for it as follows: We performed geometry optimization on the NN40, NN41, and NN42 ligand geometries generated with Avogadro while specifying a charge of -1. We then conducted single point energy calculations on the optimized anionic ligand geometries at two different charges: neutral and -2. We specified these ligands to be in a singlet state when at -1 charge and in a doublet state otherwise.

Additional CSD details
In a ConQuest search, we searched for structures containing an iridium atom bonded to two carbon atoms and four nitrogen atoms. We set bond types to "Any" and the cyclicity of the carbon and nitrogen atoms to "Cyclic." We searched with 3D coordinates determined and an R factor <= 0.05. Complexes in the hitlist were exported as mol2 files with the "Export largest molecule only" and "One file per entry" options, where the former option removes solvent molecules. From 700 hits, those with multiple iridium atoms and no iridium atoms were removed from consideration, as were hits that still had solvent or counterions, hits that did not have three bidentate ligands, duplicate hits as determined by CSD refcodes, and hits that had a CC ligand (Table S31).
Classification of CSD ligands as CN or NN was accomplished by analyzing the exported mol2 files using molSimplify, which identifies coordinating atoms in addition to identifying ligands. Hydrogen atoms were added to coordinating carbons of CSD CN ligands using molSimplify. The presence or absence of a CSD ligand in the HLS was determined through atomweighted molecular graph determinants, computed using molSimplify. CSD ligands absent from the HLS will be referred to as out-of-HLS CSD ligands. Atom-weighted molecular graph determinants were also used to ensure each out-of-HLS CSD ligand structure was only considered once in hypothetical complex enumeration, since multiple hit complexes might have CN or NN ligands in common. Hypothetical complexes either had two HLS CN ligands and an out-of-HLS CSD NN ligand, two out-of-HLS CSD CN ligands and an HLS NN ligand, or two out-of-HLS CSD CN ligands and an out-of-HLS CSD NN ligand. The two CN ligands were always the same in any given complex. From the final complexes in Table S31, 153 unique out-of-HLS CN ligands and 269 unique out-of-HLS NN ligands were identified. Consequently, 60,816 hypothetical Text S2. Extended description of feature sets. Explanation of feature notation.
For a given iridium complex, features in the ligand-only RAC set, the xTB set, the B3LYP DFT set, the wPBEh DFT set, and the Dice set require only the molecular geometry of the CN and NN ligand of the complex.
For RAC-style feature sets, the atomic properties considered were topology, identity, electronegativity, covalent radius, nuclear charge, group number, and number of bonds by the octet rule. 4 T is topology, I is identity, chi is electronegativity, S is covalent radius, Z is nuclear charge, Gval is group number, and NumB is number of bonds. For the RAC and CD-RAC feature sets, mc indicates metal-centered, lc indicates ligand-centered, D indicates difference (and its absence indicates product), and depth is indicated by the number in the feature name. The terms all, ax, and eq refer to the extent of a RAC, i.e. whether it spans over the whole complex or only over the axial or equatorial ligands. The four miscellaneous features that describe charge or denticity of a ligand were all removed. For the ligand-only RAC feature set, atom-wise properties, depth, and ligand type are indicated by the feature name.
For the electronic structure feature sets, IP stands for ionization potential and EA stands for electron affinity. The first coordinating nitrogen of the NN ligand (N1) is chosen such that the number of nitrogen atoms in its ring is less than or equal to the number of nitrogen atoms in the ring of the second coordinating nitrogen.

Invariant features
Invariant features, i.e. features that are the same across all complexes in the experimental training data, were removed during pre-processing. Only the Morgan, RAC, and CD-RAC feature sets had invariant features.

Morgan and Dice feature sets
A Morgan fingerprint indicates the presence of substructures in a molecule by hashing any given substructure into a X-bit integer, and effectively storing these integers as the indices of bits set to 1 in a 2 X -size bitset. 5 The Morgan feature set initially contains 4,096 features, with 2,048 bits (X=11) allocated to both the CN ligand and NN ligand; however, over 75% of these features are invariant over the training data and are dropped from the set. The Dice coefficient between two Morgan fingerprints A and B is defined as !" #$% , where a is the number of bits set to 1 in A, b is the number of bits set to 1 in B, and c is the number of bits set to 1 in both A and B 6 . We also considered the popular Tanimoto 6 coefficients ( " #$%&" ) but found these to hold less predictive power than Dice coefficients in the present application.

Revised autocorrelation functions details
RAC features can be full-scope, metal-centered, or ligand-centered; these distinctions indicate the positions of the starting atoms used to generate the RAC feature. For a full-scope RAC feature, every atom in the complex can be used as the starting atom. For a metal-centered RAC feature, the metal center serves as the starting atom. For a ligand-centered RAC feature, the coordinating atoms on the ligands serve as starting atoms.
Mathematically, autocorrelations are defined as ) where ! is the graph autocorrelation for property at depth , ! $ is the analogous difference graph autocorrelation, is the Dirac delta function, and "# is the bond-wise path distance between atoms and .

RAC-style feature sets details
For the RAC feature set, we allowed d to range from zero -corresponding to the correlation of an atom property to itself -to three. The RAC feature set consists of 196 features.
For the ligand-only RAC feature set, we allowed d to range from zero to four. We used a larger maximum depth than that of the RAC feature set due to the comparatively low number of features in the ligand-only case. This feature set was motivated by our anticipation that separate treatment of CN and NN ligands would improve predictive power. Furthermore, the absence of metal-centered information in this feature set is acceptable because it was anticipated that the metal-centered RACs in the RAC feature set would be uninformative -all complexes in this study have an iridium metal center and an identical first coordination sphere of two carbon atoms and four nitrogen atoms. The ligand-only RAC feature set consists of 70 features.
For the CD-RAC feature set, we allowed d to range from zero to three. Both the RAC feature set and the ligand-only RAC feature sets can be generated from structures that are not geometry optimized. This is because RACs are connectivity-dependent features but do not depend on geometry information such as bond lengths. In contrast, CD-RACs are affected by geometry optimization. CD-RACs were calculated on molSimplify-generated structures optimized with UFF. This feature set was motivated by the impact of bond distances on Ir phosphor properties like quantum yield. 7 The CD-RAC feature set consists of 222 features. Figure S1. Histograms of the three target properties across the 1,380 complexes reported in the experimental study of DiLuzio et al. 3 , excluding the baseline solvato complexes that contain a DMSO ligand. The range of phosphorescence lifetimes is restricted to omit eleven outliers with long lifetimes ranging from 14 to 75 μs. Complexes with low spectral integral (less than 1 × 10 5 photon counts) are considered dim in this work and are excluded from further Em50/50 and lifetime analysis (see Methods). Figure S2. The signed Pearson correlation coefficients (-1 in blue to +1 in red with gray for 0 as indicated in the colorbar on the right) of the xTB features with each other across the original dataset in the experimental study of DiLuzio et al. 3 , excluding the baseline solvato complexes that contain a DMSO ligand and complexes with a spectral integral below 10 5 photon counts. Figure S3. The test set performance of ANNs trained on different feature sets in predicting lifetime (in units of μs) for both random (red bars) and grouped splits (blue bars). The Morgan feature set leads to the best performance on the random split, and the RAC feature set leads to the best performance on the grouped split. Figure S4. The test set performance of ANNs trained on different feature sets in predicting spectral integral (in units of photon counts) for both random (red bars) and grouped splits (blue bars). The Dice feature set leads to the best performance on the random split, and the Morgan feature set leads to the best performance on the grouped split. Figure S5. The uncertainty quantification (UQ) cutoff versus test set mean absolute error (in μs) and data fraction of the random split ANN model trained on the xTB feature set and predicting lifetime. The data fraction is the number of test set complexes under the corresponding UQ cutoff, and the MAE is calculated on this subset of complexes. The UQ metric used is the average latent space distance to the ten nearest neighbors in the training set following the protocol introduced in Ref. 8 . The UQ metric is normalized such that the largest UQ metric is scaled to 1. Figure S6. The uncertainty quantification (UQ) cutoff versus test set mean absolute error (in photon counts) of the random split ANN model trained on the xTB feature set and predicting spectral integral. The data fraction is the number of test set complexes under the corresponding UQ cutoff, and the MAE is calculated on this subset of complexes. The UQ metric used is the average latent space distance to the ten nearest neighbors in the training set following the protocol introduced in Ref. 8 . The UQ metric is normalized such that the largest UQ metric is scaled to 1. Figure S7. The signed Pearson correlation coefficients (-0.5 in blue to +0.5 in red with gray for 0 as indicated in the colorbar on the right) of the xTB features with the phosphor properties across the original dataset in the experimental study of DiLuzio et al. 3 , excluding the baseline solvato complexes that contain a DMSO ligand and complexes with a spectral integral below 10 5 photon counts.   white for hydrogen, gray for carbon, blue for nitrogen, red for oxygen, and green for chlorine. Figure S11. Sixteen CSD ligands that lead to extreme predicted phosphor properties. The six letter identifiers are CSD refcodes. Coordinated nitrogen (carbon) atoms are indicated with blue (gray) circles. ACUROI NN appears frequently in hypothetical complexes with high predicted spectral integral. CIDDAX NN is present in the hypothetical complex with the lowest predicted lifetime. CIGKIP CN is present in the hypothetical complex with the highest predicted spectral integral. FEQSEB NN appears frequently in hypothetical complexes with high predicted Em50/50 and lifetime. KAFPEO NN appears frequently in hypothetical complexes with low predicted spectral integral. LEZJAD NN is present in two of the hypothetical complexes with the lowest predicted Em50/50. LISMIK NN appears frequently in hypothetical complexes with low predicted Em50/50. MAXWIS CN appears frequently in hypothetical complexes with low predicted lifetime and is present in the hypothetical complex with the third lowest predicted lifetime. MIMYEO NN appears frequently in hypothetical complexes with high predicted lifetime. OVALEE NN appears frequently in hypothetical complexes with low predicted lifetime. QEQVOA NN appears frequently in hypothetical complexes with low predicted spectral integral. RADTEZ CN is present in the three hypothetical complexes with the highest predicted Em50/50 and in the hypothetical complex with the second highest predicted spectral integral. RASGAV NN appears frequently in hypothetical complexes with high predicted Em50/50 and lifetime and is present in the hypothetical complex with the third highest predicted spectral integral and the hypothetical complex with the third highest predicted lifetime. TOTPAW NN appears frequently in hypothetical complexes with low predicted Em50/50 and lifetime and is present in the hypothetical complex with the highest predicted lifetime and the hypothetical complex with the second lowest predicted lifetime. TUZHEE NN appears frequently in hypothetical complexes with low predicted Em50/50 and is present in the hypothetical complex with the second lowest predicted Em50/50. YUWWOD NN appears frequently in hypothetical complexes with high predicted spectral integral and Em50/50 and is present in the hypothetical complex with the third highest predicted Em50/50 and the hypothetical complex with the second highest predicted lifetime. Triplet geometries are shown in brown, while singlet geometries are colored normally with white for hydrogen, gray for carbon, blue for nitrogen, red for oxygen, light blue for fluorine, yellow for sulfur, green for chlorine, and dark blue for iridium. The RMSD values between these singlettriplet structure pairs are 0.065 Angstroms, 0.188 Angstroms, and 0.164 Angstroms, respectively.  Figure S14. Comparison of random split ANN and TDDFT lifetime prediction to experiment (in μs) across 26 test set iridium complexes in the experimental dataset. These complexes were chosen to span the range of emission energies and lifetimes of the full set. TDDFT was carried out on optimized singlet geometries using the B3LYP functional. Lifetime was calculated using excitation energy and transition dipole moment output from the TDDFT calculation (see Methods). The dotted line is included as a reference and corresponds to perfect agreement between prediction and experiment. Figure S15. Comparison of random split ANN and TDDFT Em50/50 predictions (in eV) across 21 hypothetical iridium complexes. TDDFT was carried out on optimized singlet geometries using the B3LYP functional. The energy of the three lowest triplet sublevels was averaged for the TDDFT energy, which is used to approximate Em50/50. Figure S16. Comparison of random split ANN and TDDFT lifetime predictions (in μs) across 21 hypothetical iridium complexes. TDDFT was carried out on optimized singlet geometries using the B3LYP functional. Lifetime was calculated using excitation energy and transition dipole moment output from the TDDFT calculation (see Methods). Figure S17. Confusion matrices indicating the agreement between different lifetime labels. 0 corresponds to a lifetime τ <= 2 μs and 1 corresponds to τ > 2 μs. Over the 26 representative test set experimental complexes, we show experimentally measured lifetimes versus random split ANN predictions (top) and TDDFT predictions (middle). Over the 21 hypothetical complexes, we show TDDFT predictions versus random split ANN predictions (bottom). The TDDFT results were generated with singlet geometries and the B3LYP functional.  (3) lc-I-0-eq, lc-I-1-eq, lc-T-0-eq 25 difference (10) D_lc-chi-0-eq, D_lc-Z-0-eq, D_lc-I-0-eq, D_lc-I-1-eq, D_lc-I-2-eq, D_lc-I-3-eq, D_lc-T-0-eq, D_lc-S-0-eq, D_lc-Gval-0-eq, D_lc-NumB-0-eq Table S6. The seventy features in the ligand-only RAC feature set. Each combination of ligand and atomic property yields five features due to depths ranging from 0 to 4. Features and notation are described in detail in Text S2.  ligand  atomic property  feature list  CN  topology  T-0_CN, T-1_CN, T-2_CN,  T-3_CN, T-4_CN  identity  I-0_CN, I-1_CN, I-2_CN, I-3_CN, I-4_CN  electronegativity  chi-0_CN, 4_NN  Table S7. The 222 features in the CD-RAC feature set. Invariant features are removed from the original twenty-eight features of that category (seven atomic properties and four depths, i.e. 0, 1, 2, and 3). The notation of features and property type is the same as in RACs but incorporates a distance-dependent term as outlined in ref. 10  difference (10) D_lc-chi-0-eq, D_lc-Z-0-eq, D_lc-I-0-eq, D_lc-I-1-eq, D_lc-I-2-eq, D_lc-I-3-eq, D_lc-T-0-eq, D_lc-S-0-eq, D_lc-Gval-0-eq, D_lc-NumB-0-eq    Table S31. Statistics on structure attrition for hit CSD complexes used to identify CN and NN ligands outside of the HLS. There were six reasons why any CSD complex could be eliminated from consideration: the presence of multiple iridium atoms; the absence of any iridium atoms due to the presence of another molecule larger than the iridium complex in the CSD entry (combined with the "Export largest molecule only" option, as described in Text S1); the incomplete removal of solvent or counterions leading to index errors; the presence of non-bidentate ligands such that a complex was not 2-2-2, i.e. did not have three bidentate ligands in octahedral geometry; refcode duplicates such as HULVEQ and HULVEQ01; and the presence of a CC ligand.