Low-cost machine learning prediction of excited state properties of iridium-centered phosphors

Prediction of the excited state properties of photoactive iridium complexes challenges ab initio methods such as time-dependent density functional theory (TDDFT) both from the perspective of accuracy and of computational cost, complicating high-throughput virtual screening (HTVS). We instead leverage low-cost machine learning (ML) models and experimental data for 1380 iridium complexes to perform these prediction tasks. We find the best-performing and most transferable models to be those trained on electronic structure features from low-cost density functional tight binding calculations. Using artificial neural network (ANN) models, we predict the mean emission energy of phosphorescence, the excited state lifetime, and the emission spectral integral for iridium complexes with accuracy competitive with or superseding that of TDDFT. We conduct feature importance analysis to determine that high cyclometalating ligand ionization potential correlates to high mean emission energy, while high ancillary ligand ionization potential correlates to low lifetime and low spectral integral. 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 use uncertainty-controlled predictions to identify promising ligands for the design of new phosphors while retaining confidence in the quality of the ANN predictions.


Contents Table S1
Cyclometalating ligand structures and labels Page S3 Table S2 Ancillary ligand structures and labels Page S7 Figure S1 Histograms of the three target properties Page S9 Text S1 Electronic structure method and CSD details Page S10 Text S2 Feature set details Page S12 Page S17 Table S8 Features in the electronic structure feature sets Page S18 Figure S2 Correlation between xTB features Page S19 Table S9 Performance of different charge features Page S20 Figure S3 Feature set performance on lifetime Page S20 Figure S4 Feature set performance on spectral integral Page S21  Figure S5 Uncertainty quantification for the lifetime xTB ANN Page S26 Figure S6 Uncertainty quantification for the spectral integral xTB ANN Page S27 Table S20 Comparison of different ML models Page S27 Figure S7 Correlation between xTB features and target properties Page S28 Figure S8 Histograms of xTB features Page S29 Figure S9 Ligand substitution effect on lifetime Page S30 Figure S10 Ligand substitution effect on spectral integral Page S30 Figure S11 Sixteen CSD ligands that lead to extreme hypothetical complexes Page S31   Figure S16 Phosphor geometry differences for singlet vs triplet state Page S42 Figure S17 TDDFT-calculated energy with CAM-B3LYP Page S43 Figure S18 TDDFT-calculated energy with ωB97X-D3BJ Page S43

References
Page S44  Figure S1. Histograms of the three target properties across the 1,380 complexes reported in the experimental study of DiLuzio et al., 1 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 main text Computational Details).
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. 2,3 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
Neutral complexes are generated when NN40, NN41, or NN42 is the ancillary ligand. These three ancillary ligands contain tetrazole or pyrazole moieties that deprotonate upon metal coordination. 1 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 S32).
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 identical HLS CN ligands and an outof-HLS CSD NN ligand, two identical out-of-HLS CSD CN ligands and an HLS NN ligand, or two identical out-of-HLS CSD CN ligands and an out-of-HLS CSD NN ligand. From the final complexes in Table S32, 153 unique out-of-HLS CN ligands and 269 unique out-of-HLS NN ligands were identified. Consequently, 60,816 hypothetical complexes were considered. Application of uncertainty quantification cutoffs left 70 unique out-of-HLS CN ligands and 42 unique out-of-HLS NN ligands spread out over 3,598 hypothetical complexes, which we analyzed with our ANNs.
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, which is also sometimes referred to as the eRAC set when the latter two features are included. 4,5 T is topology, I is identity, c 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. 6 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. 7 We also considered the popular Tanimoto 7 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 property for atom , is the Dirac delta function, and "# is the bond-wise path distance between atoms and .
For CD-RACs, 8 there is an adjustment for the spatial distance between two atoms. CD-RACs are defined as where !,&' is the Coulomb-decay graph autocorrelation for property at depth , is the number of atoms, ' is the property for atom , '( is the spatial distance between atoms and , is the Dirac delta function, and "# is the bond-wise path distance between atoms and . Difference CD-RACs are analogously defined.

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 total.
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. 9 The CD-RAC feature set consists of 222 features.  (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 eRAC 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  4_NN  Table S7. The 222 features in the CD-RAC feature set. Each category of features contains twentyeight features to start (i.e., seven atomic properties at four depths, 0, 1, 2, or 3). Invariant features that are removed are listed below along with a final feature count. The notation of features and property type is the same as in RACs but incorporates a distance-dependent term as outlined in ref. 8 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. 11. 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. 11. The UQ metric is normalized such that the largest UQ metric is scaled to 1.     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 sixletter identifiers are CSD refcodes of each complex from which the ligand was extracted. 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 secondhighest 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 thirdhighest 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 secondhighest predicted lifetime.  Figure S12. 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 (S0) geometries using the B3LYP functional. Lifetime was calculated using excitation energy and transition dipole moment output from the TDDFT calculation (see main text Computational Details). The dotted line is included as a reference and corresponds to perfect agreement between prediction and experiment.    True  False  True  True  True  True  res  True  True  True  True  False  False   Table S32. 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.