Exploiting Ligand Additivity for Transferable Machine Learning of Multireference Character Across Known Transition Metal Complex Ligands

Accurate virtual high-throughput screening (VHTS) of transition metal complexes (TMCs) remains challenging due to the possibility of high multi-reference (MR) character that complicates property evaluation. We compute MR diagnostics for over 5,000 ligands present in previously synthesized transition metal complexes in the Cambridge Structural Database (CSD). To accomplish this task, we introduce an iterative approach for consistent ligand charge assignment for ligands in the CSD. Across this set, we observe that MR character correlates linearly with the inverse value of the averaged bond order over all bonds in the molecule. We then demonstrate that ligand additivity of MR character holds in TMCs, which suggests that the TMC MR character can be inferred from the sum of the MR character of the ligands. Encouraged by this observation, we leverage ligand additivity and develop a ligand-derived machine learning representation to train neural networks to predict the MR character of TMCs from properties of the constituent ligands. This approach yields models with excellent performance and superior transferability to unseen ligand chemistry and compositions.

Of these MR diagnostics, those derived from wavefunction theory (WFT) are generally more predictive than those derived from DFT. 53 However, WFT calculations scale at least N 4 with the system size and thus are impractical in VHTS of transition metal chemistry. Moreover, certain MR diagnostics have been observed to be more transferable over different chemical spaces (i.e. organic molecules vs. TMCs) than others. [54][55][56] Therefore, data-driven approaches [57][58][59][60][61] have been developed to bridge the gap between DFT-and WFT-based diagnostics 53 , making system-specific MR/SR decisions for method selection 39 , and automating the active space selection for MR WFT calculations 57,60,62 .
Strong MR character in TMCs still poses challenges for VHTS. Most studies that investigate the MR character of TMCs focus on model systems with a few ligands. This is likely due to the prohibitive computational cost of studying a large data set of realistic, synthetically accessible TMCs with multiple WFT calculations. In addition, despite the fact that the MR character of both TMCs and organic molecules have been studied in the literature, the relationship between the MR character of a TMC and that of its constituent ligands remains unknown.
In this work, we study TMCs and thousands of their constituent ligands that are known to be synthetically accessible by virtue of their presence in the Cambridge Structural Database (CSD) 63 . We devise a strategy to quantify the total MR character of TMCs based on ligandderived measures of MR character. This approach makes tractable the diagnosis of high MR character in TMCs. To accomplish this task over a large set of known ligands, we first introduce an iterative approach for ligand charge assignment on CSD ligands based on all charge and oxidation state information available within the CSD. We then perform high-throughput calculations to obtain MR diagnostics for 5,163 synthetically accessible ligands. For these ligands, we observe that ligand MR character linearly correlates with the inverse value of the averaged bond order over all bonds in the ligand. We demonstrate that the concept of ligand additivity holds for MR character, which suggests that the MR character of a TMC can be inferred from the sum of the MR character of the ligands. Leveraging this ligand additivity, we develop a new ligand-derived ML representation to train ANNs to predict TMC MR character.
These models demonstrate excellent performance and superior transferability to unseen ligand chemistry and compositions in comparison to alternative approaches.

Data Sets
We searched the CSD 63 version 5.41 (Nov. 2019) + 3 data updates to find octahedral mononuclear complexes. This query returned 28,006 complexes that have user-defined charge and oxidation states. We used molSimplify 64,65 to deconstruct these complexes into metals and their corresponding ligands, which resulted in 9,279 unique ligands. Duplicate ligands were eliminated using our previous molecular graph determinant 66 approach. We applied our iterative CSD-consistent charge assignment approach on these complexes and successfully assigned the charge of 8,805 ligands (see Results and Discussions 4.1). We excluded ligands that have more than 25 heavy atoms and those that contain uncommon chemical elements, as judged by their low abundance (< 0.1%) in the ligands of mononuclear octahedral complexes in CSD. After removing these ligands from the overall set, we were left with 6,131 ligands in total (Supporting Information Table S1).
We compared our CSD-consistent charge with charge assignment using the octet rule 67 .
Only 5,905 ligands can be assigned a charge using the octet rule, because the remaining examples lack bond order information. We identified 5,713 cases where our assigned charge matches an octet rule charge and 192 cases where the two charges differed. We manually inspected all 192 cases where the charges differed and observed that the most common reason for a mismatch was due to incorrect treatment of radicals or ligands containing multiple atoms that have variable oxidation states (e.g., B and P) in the octet rule assignment (Supporting Information Figure S1). We confidently identified 57 of these 192 cases as radicals or diradicals by inspecting the corresponding literature. A complete list of these ligands is provided in a .csv file in the Supporting Information. We then computed the MR diagnostics (see MR diagnostics calculations) of the 5,713 closed-shell ligands and 57 radicals or diradicals following our established procedure 53 . If we could not successfully compute all 14 MR diagnostics (e.g., due to self-consistent field convergence failure), we removed that ligand from the data set (Supporting   Information Table S2). Our final data set contains 5,163 ligands from octahedral mononuclear CSD complexes for which we have ascertained charges and computed MR diagnostics, which we refer to as OctLig.
To investigate the effect of mixing ligands with distinct MR character, we also selected ligands from the OctLig set that have the strongest or weakest MR character for their respective denticities of 1, 2, 5, or 6, which corresponds to a set of 8 ligands in total (Supporting Information Figure S2). We combined these ligands with Fe centers to make octahedral complexes. We considered both Fe(II) and Fe(III) in high-spin (HS) (i.e., quintet or sextet) or low-spin (LS) (i.e., singlet or doublet) states and all possible ligand symmetries, which produces a set of 172 transition metal complexes that we refer to as extTMC (Supporting Information   Table S3).
We also assembled 1,000 octahedral complexes from monodentate ligands in the OctLig set under the constraints that the complex consists of up to 2 unique ligands, that it was no more than 40 atoms in size, and had a net charge of -4 to +2. Each complex contains a mid-row transition metal (Cr, Mn, Fe, or Co) in oxidation state II or III and in its HS or LS state (Supporting Information Table S4). We performed DFT geometry optimizations on these 1,000 TMCs, and 855 of them converged to suitable octahedral geometries (Supporting Information   Table S5). We refer to this set of 855 complexes from diverse CSD ligands as divTMC.

DFT Geometry Optimizations.
Density functional theory (DFT) geometry optimizations with the B3LYP 68-70 global hybrid functional were carried out using a developer version of the graphical processing unit (GPU)-accelerated electronic structure code TeraChem [71][72][73] . The LANL2DZ effective core potential 74 basis set was used for metals, Br, and I, and the 6-31G* basis for all other atoms.
Singlet spin states were calculated with the spin-restricted formalism while all other calculations were carried out in a spin-unrestricted formalism. All DFT calculations employed level shifting 75 of 0.25 Ha to both majority-and minority-spin virtual orbitals to aid self-consistent field (SCF) convergence. All geometry optimizations were performed with the L-BFGS algorithm in translation rotation internal coordinates (TRIC) 76 to the default tolerances of 4.5 × 10 −4 Hartree/bohr for the maximum gradient and 1 × 10 −6 Hartree for the energy change between steps.
For the OctLig set, we selected the initial ligand geometry as the structure that produced the lowest B3LYP single-point energy among the multiple copies of the ligand present in the CSD. In the subsequent geometry optimization of that structure, we only optimized the location of H atoms, which are not necessarily resolved in the crystal structure, while we froze all heavy atoms. For the extTMC and divTMC sets, we assembled transition metal complex initial geometries using molSimplify 64,65 to combine the optimized ligand structures from the OctLig set with transition metals in octahedral geometries. Geometry checks 77,78 were applied to eliminate optimized TMC structures that deviated from the expected octahedral shape following previously established metrics without modification 77,78 (Supporting Information Table S5).

MR Diagnostic Calculations.
Following our prior studies 39 Information Tables S6-S7). These DFT-based and WFTbased diagnostics are all explicitly calculated, unlike some previous work where we trained ML models to predict the higher cost WFT-based diagnostics from a combination of descriptors and DFT-based diagnostics. 79 We used our established workflow 39,54,79 to ensure convergence to a consistent electronic state when evaluating different MR diagnostics (Supporting Information Text S1 and Figure S3). We used the restricted formalism for closed-shell species and restricted open-shell formalism otherwise. We chose the percentage of correlation energy recovered by  Information Table   S8).

ML Model Training.
We first partitioned the divTMC set with a random 80%/20% train/test split, with 20% of the training data (i.e., 16% overall) used as the validation set for hyperparameter selection using HyperOpt 81 with 200 evaluations (Supporting Information Table S9). To examine the transferability of our newly developed ligand diagnostic revised autocorrelation functions (RACs), we further applied two non-random partition schemes on the divTMC set based on ligand identity to ensure distinct unseen chemistry in the test set. First, 50 ligands (out of 344 ligands in total) were chosen, and all complexes that contain any of those 50 ligands were partitioned into a test set (Supporting Information Table S10). Second, all complexes that contained any ligand with a 3p (i.e., period 3) connecting atom (e.g., Si, S, P, or Cl) were partitioned into the test set (Supporting Information Table S10). For these two partitioning schemes, we adopted the hyperparameters from hyperparameter optimization for the random partition (Supporting Information Table S9). All ANN models were trained up to 2,000 epochs with the Adam optimizer using Keras 82 with Tensorflow 83 as a backend, and dropout, batch normalization, and early stopping to avoid over-fitting.

Consistent Ligand Charge Assignment Using Collective Information in the CSD.
There are 28,006 octahedral mononuclear transition metal complexes deposited in the Cambridge Structural Database (CSD) 63 (version 5.41 (Nov. 2019) + 3 data updates) that have user-labeled complex charges and metal oxidation states. This information facilitates computational studies of their properties (e.g., HOMO-LUMO gap) that can be obtained from DFT. 84 Similar by-ligand charge information, however, is absent for the ligands that comprise these CSD complexes, preventing large-scale computational screening of individual synthetically accessible ligands. Thus, to assess a ligand's MR character and how it correlates to the MR character of TMCs that contain it, we need a method that assigns charges of individual ligands.
We thus developed an iterative approach to obtain the ligand charges that is consistent over multiple copies of the ligand in the CSD (Supporting Information Text S2). This strategy applies the constraint that the total complex charge should be the sum of the metal oxidation state and the sum of all constituent ligand charges (Supporting Information Text S2). We start with homoleptic complexes (i.e., with only one unique ligand) and directly derive the charge of the ligand involved. Next, we search for complexes that only have one ligand with unknown charge while others have their charges assigned from a previous iteration of our algorithm. This process continues until we cannot obtain charges for any new ligands ( Figure 1). This iterative approach has several advantages. First, the deduced ligand charges are consistent across the CSD and are thus likely to be correct. Second, the number of ligands with known charge will increase as more structures are deposited in the CSD (Supporting Information Figure S4). Lastly, we can use the same logic to obtain either the metal oxidation state or the total complex charge for complexes when this information is not available in the deposited CSD structure. Ligands are categorized into three groups: successfully assigned (blue), failed assignment due to conflicting information in CSD entries (red), and unknown (gray). An illustration of the iterative CSD-consistent charge assignment for the example of two complexes, [Fe(II)(H2O)6] 2+ and neutral Fe(II)(C8S8)(H2O)4, is also shown (top). (right) Examples where CSD-consistent charge assignment will fail: a homoleptic Fe(II) complex (refcode: DIXDIX) with two identical tridentate ligands that leads to a non-integer charge assignment for each tridentate ligand (i.e., 0.5, top) and a Ni(II) complex (refcode: APZPQN) that is comprised of two ligands, both of which rare in the CSD. Atoms are colored as follows: orange for Fe, green for Ni, yellow for S, gray for C, blue for N, red for O, and white for H. With this iterative approach, we successfully obtain CSD-consistent charges for 96.3% of ligands seen in octahedrally coordinated complexes in the CSD ( Figure 1). The remaining 3.7% of ligand failures are mostly due to mislabeling or lack of information in the CSD. For 0.7% of the ligands, we failed to assign the charge due to potential mislabeling of the total complex charge and/or the metal oxidation state by the user who deposited the crystal structure. For example, an Fe(II) complex (refcode: DIXDIX) is reported with a total complex charge +3 but contains two identical tridentate ligands, leading to a fractional ligand charge assignment (i.e., 0.5) for each tridentate ligand ( Figure 1). The remaining 3.0% of the ligands are rarely present in the CSD and thus their charges are challenging to determine through inference. For example, a Ni(II) complex (refcode: APZPQN) contains two ligands that only appear once in all octahedral complexes in the CSD, leading difficulties for their CSD-consistent charge assignment ( Figure 1).
For the ligands where we are able to assign a CSD-consistent charge, we obtain good agreement (96.7%) with charges derived from the octet rule. 67 When the CSD-consistent charge disagrees with the octet rule charge, the ligands are either radicals, diradicals, or contain multiple atoms that have variable oxidation states (e.g., B and P), where the octet rule is expected to fail (Supporting Information Figure S1).  Due to the good linear correlation between %Ecorr[(T)] and the inverse average BO, ligands with higher denticities have weaker MR character, because they tend to be more saturated than ligands with lower denticities (Figure 3). This observation suggests that single-site catalysts with high-denticity ligands are preferable targets for VHTS not only due to their ease of their synthesis and introduction of steric effects that can be beneficial "knobs" to tune in catalyst design 85,86 but also due to their simpler electronic structure. Over the 5,163 unique ligands from the OctLig set, our algorithm identifies 57 radical ligands, which we confirm by literature analysis (details are provided in a .csv file in the Supporting Information .zip file). Surprisingly, these radical ligands do not have increased MR character relative to the closed-shell singlets when their inverse average bond orders are comparable ( Figure 2).

Incorporating Additivity of MR character for TMCs Improves ML Transferability.
We next investigated the relationship between the MR character of individual ligands and the MR character of the transition metal complexes (TMCs) that they comprise. We started with a controlled data set of Fe(II) complexes, extTMCs, in which all TMCs are constructed from ligands that have either strong (e.g., CS) or weak (e.g., H2O) MR character (see Data Sets). By ligands that have individually higher MR character, we conclude that the nature of rND (i.e., the ratio of non-dynamical to total correlation) intrinsically fails to capture additivity effects in MR character. The rND diagnostic also sometimes shows larger sensitivities to isomer choice than to the character of the ligands, which diverges from the other computed diagnostics and figure of merit (Supporting Information Figure S10). ( Figure 6). However, we find that this eRACs-based ANN model is not sufficiently transferable.
The eRACs/ANN model performance is worsened if tested on ligands that were held out in a group split (R 2 =0.81) or chemical elements (i.e., unseen 3p elements when only 2p elements were in the training set, R 2 =0.51) that were not present in the training set ( Figure 6 and see Data Sets). We thus next investigated whether additivity could be exploited to devise new descriptors from a ligand-focused perspective that would increase model transferability to unseen chemistry.  To improve model transferability, we designed a new set of descriptors to exploit the ligand additivity effects in MR character in TMCs, which we refer to as ligMR-RACs.
Specifically, we adapt RACs to encode the MR character of a TMC through the MR diagnostics of its fragments (i.e., metal center and ligands). Because most of the MR diagnostics cannot be localized to individual atoms, we first convert the 2D molecular graph to a fragment-based graph, where the MR diagnostics of each fragment (i.e., ligand) have been pre-computed (Supporting Information Figure S13). To compute ligMR-RACs for a TMC, we follow the same procedure as in RACs, but we instead operate on the ligand-based graph and replace the  Figure   S14). We next considered whether this combined feature set improves transferability where the eRACs/ANN model performance had been poor.
To test our model transferability, we returned to the two non-random partition schemes on the divTMC set based on ligand identity to ensure distinct unseen chemistry in the test set (see . For a grouped split where specific ligands were excluded from the training set and partitioned into the test set, adding ligMR-RACs significantly improves model performance in comparison to using only eRACs (R 2 =0.81 to R 2 =0.96, Figure 6, Supporting Information Figure   S14). Notably, adding ligMR-RACs greatly improves model accuracy when the model is only trained on the TMCs containing 2p elements but tested on those containing unseen 3p elements (R 2 =0.51 to R 2 =0.93, Figure 6, Supporting Information Table S12). For models with both eRACs and ligMR-RACs as inputs, their performance is thus robust against each of these data partitioning schemes (R 2 ranges from 0.93 to 0.96). In addition, the fraction of training data in the three partition schemes are comparable (i.e., 73%-82%), suggesting that the inclusion of ligMR-RACs to our models allows us to learn DLPNO-CCSD %Ecorr[(T)] comparably with either in-distribution (i.e., random partition) or out-of-distribution (i.e., excluding ligands and 2p/3p partition) test data (Supporting Information Table S10). This observation highlights the transferability facilitated by ligMR-RACs in predicting DLPNO-CCSD %Ecorr[(T)] for TMCs in diverse chemical spaces. As a comparison, we also trained ANN models using only ligMR-RACs. Despite the slightly worse performance in comparison to the models trained on the combined feature set, the ligMR-RAC-only models maintain good transferability we identified with ligMR-RACs. Specifically, these models have rather uniform performance in all three partitioning schemes (R 2 ranges from 0.89 to 0.93, Supporting Information Figure S14 and Table   S12). This comparison emphasizes the synergistic benefit of using both eRACs and ligMR-RACs in chemical exploration.

Conclusions
The accurate screening of transition metal complexes is challenged by the fact that they may suffer from varying degrees of strong correlation. Given the large size of transition metal complexes and their constituent ligands, it is nearly impossible to exhaustively screen these complexes with multi-reference wavefunction theory methods, motivating low-cost detection of Page S9 Text S2 Pseudocode for iterative CSD-consistent charge assignment Page S10 Figure S4 Percentage of ligands vs. percentage of CSD complexes Page S11 Figure S5 %Ecorr[(T)] vs. inverse average bond order of molecules at different sizes Page S11 Table S11 Detailed information for select ligands with triple bonds Page S12 Figure S6 Unsigned Pearson's r matrix of MR diagnostics and %Ecorr[(T)] Page S12 Figure S7 Spearman's r matrix of MR diagnostics and %Ecorr[(T)] Page S13 Figure S8 2D chemical structure of phenylmethylidyne and pyridine-N-oxide Page S13 Figure S9 MR diagnostics vs. the number of CS ligands in LS Fe(II)(CS)n(H2O)6-n Page S14 Figure S10 MR diagnostics vs. the number of CS ligands in HS Fe(II)(CS)n(H2O)6-n Page S15 Figure S11 MR diagnostics for LS "5+1" complexes Page S16 Figure S12 MR diagnostics for HS "5+1" complexes Page S17 Figure S13 2D molecular graph for RACs and ligMR-RACs Page S18 Text S3 Extended description of RACs and ligMR-RACs Page S18 Figure S14 Scaled MAE for ML predicted DLPNO-CCSD %Ecorr[(T)] Page S20  Figure S1. Examples where the octet rule gives incorrect ligand charge for a radical 1 (left, from a TMC with refcode WEQZOI) and a small cluster-containing ligand 2 (right, from a TMC with refcode EHUCUG). The CSD-consistent charge was verified by reviewing the original publication associated with the structure. Atoms are colored as the following, gray for C, blue for N, red for O, orange for P, white for H. For each ligand, a placeholder for metal to which the coordinating atoms connect is also shown as white.   Table S3. Ligand symmetries considered in the extTMC set. For each symmetry, there are four choices for the metal core (i.e., singlet Fe(II), quintet Fe(II), doublet Fe(III), and sextet Fe(III)), resulting in 43*4=172 complexes in total in the extTMC set. We used a short-hand notation, MR denticity for the ligand. For example, a strong MR character bidentate ligand is S 2 . Thus, a ligand symmetry made of two bidentate ligands of strong MR character and two monodentate ligands of weak MR character is (S 2 )2(W 1 )2. We use the bracket (i.e., []) to group isomers together. For symmetries made of two bidentate and two monodentate ligands, we require the two bidentate ligands reside in the equatorial plane.
description count symmetries Six monodentate 10 One hexadentate 2 (S 6 )1, (W 6 )1  Table S5. Geometry cutoffs used in this work, which are the same as our prior work 3 . These geometry checks include: coordination number must be 6, mean and max. deviation of connecting-atom-metal-connecting-atom angles, max. deviation of any metal-ligand bond length, max. deviation of metal-ligand bond lengths for equatorial ligands, max. root mean square deviation (RMSD) of the ligand from its initial structure, and mean and max. deviation of an expected linear ligand from a 180° angle (which are only computed for linear ligands).  Table S6. Summary of MR diagnostics grouped by type and method used. The only change from prior work 4 is that D1 and D2 diagnostics were removed considering their high linear correlations with some of the diagnostics (D1 with T1 and D2 with C0 2 ) and max(t1) was added as in our recent work 31 . nHOMO[MP2] 6,9 MP2 occupations MP2 highest occupied natural orbital occupation nLUMO[MP2] 6,9 MP2 occupations MP2 lowest unoccupied natural orbital occupation T1 10 CCSD excitations Frobenius norm of the single-excitation amplitude vector normalized by the square root of the number of electrons in CCSD max(t1) 11 CCSD excitations The largest eigenvalue of the matrix derived from the single-excitation amplitudes. %TAE[(T)] 12 CCSD(T) TAE Percent difference in TAE from CCSD vs. CCSD(T) C0 2 [14] 10, 13 CASSCF occupations CASSCF leading coefficient CSF at an active space of 14 orbitals nHOMO [14] 6,14 CASSCF occupations CASSCF highest occupied natural orbital occupation at an active space of 14 orbitals nLUMO [14] 6,14 CASSCF occupations CASSCF lowest unoccupied natural orbital occupation at an active space of 14 orbitals  Figure S3. Workflow for computing 14 MR diagnostics following our previous established procedure [16][17][18] . All calculations were performed with ORCA 4.2.1 15 .
Text S1. Workflow for the calculations of 14 MR diagnostics. We converged a B3LYP calculation and used it to initialize both DFT calculations with other density functional approximations (i.e., BLYP, B1LYP, PBE, and PBE0) as well as HF calculations. This ensured we converged a consistent electronic state across multiple calculations and simultaneously saved computational time. The converged HF wavefunction was then used for MP2 and CCSD(T) calculations following an established procedure in our previous work 17 ( Figure S3). All MP2 natural orbitals with occupations between 0.02 and 1.98 were initially selected for the CASSCF active space. If these thresholds led to an active space that would produce greater than 1 × 10 7 configuration state functions, the thresholds for occupation-based inclusion were tightened, limiting the maximum number of orbitals in the active space to 19 across the OctLig set. During the computation of total atomization energy (TAE)-based diagnostics, we assumed homolytic dissociation for the atoms in the ligands. We chose the percentage of correlation energy recovered by CCSD compared to CCSD(T) (i.e., %Ecorr  (Table S8).              Circles are used to represent atoms and squares for fragments (i.e., a metal center or ligand) in TMCs. Atoms are colored as follows: orange for Fe, gray for C, red for O, and blue for N.

Text S3. Extended description of RACs and ligMR-RACs
We have introduced a systematic approach to featurize molecular inorganic complexes that blends metal-centric and whole-complex topological properties in a feature set referred to as revised autocorrelation functions (RACs). 21 These RACs, variants of graph autocorrelations (ACs), 22 are sums of products and differences of atomic properties, i.e., electronegativity ( ), nuclear charge (Z), topology (T), covalent radius (S), and identity (I). Recently, we incorporated the group number (G) as an additional atomic property for better generalization across rows in the periodic table. We name these new sets of descriptors eRACs 23 . Standard ACs are defined as where Pd is the AC for property P at depth d, δ is the Dirac delta function, and dij is the bondwise path distance between atoms i and j.
In our approach, we have five types of RACs: • : standard ACs start on the full molecule (f) and have all atoms in the scope (all).
• and : restricted-scope ACs that start on the full molecule (f) and separately evaluate axial or equatorial ligand properties where nax/eq is the number of atoms in the corresponding axial or equatorial ligand and properties are averaged within the ligand subtype.
• : restricted-scope, metal-centered (mc) descriptors that start on the metal center (mc) and have all atoms in the scope (all), in which one of the atoms, i, in the i,j pair is a metal center: • :and : restricted-scope, metal-proximal ACs that start on a ligand-centered (lc) and separately evaluate axial or equatorial ligand properties, in which one of the atoms, i, in the i,j pair is the metal-coordinating atom of the ligand: We also modify the AC definition, Pʹ, to property differences rather than products for a minimum depth, d, of 1 (as the depth-0 differences are always zero): where scope can be axial, equatorial, or all ligands.
Here, we adapt RACs to encode the MR character of a TMC through the MR diagnostics of its fragments (i.e., metal center and ligands). Since most of the MR diagnostics cannot be localized to individual atoms, we first convert the 2D molecular graph to a fragment-based graph, where the MR diagnostics of each fragment have been pre-computed ( Figure S13). To compute ligMR-RACs for a TMC, we follow the same routine in RACs, but simply operating on this fragmentbased graph and replacing the atomic property in RACs by ligand MR diagnostics. Since the maximum bond depth (d) is 2 in fragment-based graph for all TMCs, we truncate ligMR-RACs at depth 2 as well. We remove four MR diagnostics (i.e., IND[B3LYP], rND[B3LYP], nHOMO [MP2], and nHOMO [14]) that have Pearson's r coefficients larger than 0.95 with at least one of the remaining diagnostics, leaving 10 diagnostics and %Ecorr[(T)] ( Figure S6). Due to the constraints on our fragment-based molecular graph,