Accurate multi-objective design in a space of millions of transition metal complexes with neuralnetwork-driven efficient global optimization
The accelerated discovery of materials for real world applications requires the achievement of multiple design objectives. Satisfaction of such constraints requires exploration of multi-million compound libraries over which even density-functional theory (DFT) screening is intractable. Machine learning (ML, e.g., artificial neural network, ANN, or Gaussian process, GP) models for this task are limited by training data availability and predictive uncertainty quantification (UQ). We overcome such limitations by using efficient global optimization (EGO) with the multi-dimensional expected improvement (EI) criterion. EGO balances exploitation of a trained model with acquisition of new DFT data at the Pareto front, the region of chemical space that contains the optimal trade-off between multiple design criteria. We demonstrate this approach for the simultaneous optimization of redox potential and solubility in candidate M(II)/M(III) redox couples for redox flow batteries from a space of 2.8M transition metal complexes designed for stability in practical RFB applications. We employ latent-distance-based UQ with a multi-task ANN to enable model generalization that surpasses that of a GP. With this approach, ANN prediction and EI scoring of the full 2.8M complex space is achieved in minutes. Starting from ca. 100 representative points, EGO improves both properties by 3-4 standard deviations in only five generations. Analysis of lookahead errors confirms rapid ANN model improvement during the EGO process, achieving suitable accuracy for predictive design in the space of transition metal complexes. The ANN-driven EI approach achieves at least 500-fold acceleration over random search, identifying a Pareto-optimal design in around five weeks instead of fifty years.
Introduction<!>Results and Discussion<!>Figure 1.<!>Conclusions<!>Methods
<p>The discovery of new functional materials is essential to address outstanding challenges in energy and resource utilization. Computational screening, e.g., with first-principles density functional theory (DFT), is a powerful complement to experimental efforts due to the ease and speed with which properties can be evaluated. Materials design often necessitates that multiple objectives be satisfied, such as meeting an energetic criterion while also exhibiting stability or low cost. 1 However, satisfying these multiple criteria will necessitate exploration of a wide space (ca. millions) of candidates, making even efficient computational screening intractable. In these cases, it is beneficial to optimize properties along a Pareto front, the region of chemical space that contains the optimal tradeoffs between desired properties, as has occasionally been carried out in chemistry [2][3][4][5][6] , e.g., for catalyst design 2 .</p><p>In recent years, machine learning (ML) has emerged as a promising alternative to direct, physics-based simulation by enabling property prediction in seconds that would normally take days or weeks for full DFT-based evaluation 1,7 . Although ML models were first demonstrated to achieve near-DFT accuracy in small molecule organic chemistry 8 , recent developments have focused on challenging materials spaces, such as transition metal chemistry 7,[9][10][11][12][13][14][15][16] and inorganic materials [17][18][19][20][21][22] . Challenges remain to deliver on the promise of ML-accelerated materials discovery. Prior to property prediction, ML models still require substantial training data, typically on the order of thousands [23][24] to millions [25][26] of DFT-evaluated properties. Even when such models are carefully trained, performance on a set aside test set can differ 10,14,[27][28] from the true generalization performance to out-of-sample compounds that are the typical target of discovery efforts.</p><p>Active learning has been recognized 5,[29][30][31][32][33][34][35][36][37] as an attractive paradigm for balancing between data acquisition in ML model training (i.e., exploration) and ML-model-based prediction (i.e., exploitation). Within chemical discovery 5,[30][31][32][33][34][35][36][37] , Bayesian optimization [38][39] approaches, such as efficient global optimization (EGO) with the expected improvement (EI) criterion [33][34][35][36][37] , have been employed to balance exploration with exploitation. In cases where a suitable metric for model uncertainty is available, these approaches have been fruitfully applied to accelerate ML model training 29 and single-objective discovery efforts [33][34][35][36][37] .</p><p>In this work, we develop an approach to multi-objective transition metal complex design with application to redox flow batteries (RFBs). RFBs [40][41] decouple power delivery and cell capacity by storing the redox-active species separately and delivering or replacing them on demand, making RFBs ideal for economical, large-scale energy storage [42][43] . Transition metal ions (e.g., vanadium) or complexes have been widely used in RFBs due to their ability to cycle through a range of oxidation states. However, these smaller species with high redox potentials tend to diffuse across the membrane in a phenomenon known as crossover. Complexation with bulkier organic ligands can both tune redox properties and prevent diffusion 44 , but these same ligands can limit solubility 45 in the polar, organic solvents (e.g., acetonitrile) favored for RFB electrolytes.</p><p>Computational optimization of transition metal complex redox couples for RFBs should therefore address crossover, stability, solubility, and redox potential. We address the first two challenges by constructing a 2.8M compound space of 700k candidate bulky ligands that are expected to form stable complexes with four open shell transition metal ions. To optimize solubility and redox potential within this space, we then employ two-dimensional EI to balance exploration and exploitation at the Pareto front. We show that latent-distance-based UQ metrics 14 with a multi-task ANN generalize better than a Gaussian process, achieving good predictive accuracy from a combination of data sets from prior work 7,[9][10]13 and representative points from the design space. This active learning approach rapidly improves properties and model prediction errors in only a handful of generations, revealing design rules for RFB redox couples satisfying multiple objectives. The ANN-driven EI approach achieves at least 500-fold acceleration over random search, identifying a Pareto-optimal design in weeks instead of decades.</p><!><p>We begin by reviewing the multiple concerns we must address in the design of transitionmetal complex-based redox couples for redox flow batteries (RFBs). To address the challenge of crossover and stability, we created a smoothly varying but diverse complex space of nearly three million unique candidate redox couples with bulkier ligands that should avoid the problem of crossover (Figure 1). This space was created through a series of hierarchical steps that started with core heterocyclic ligands inspired by experiment [46][47] that we subsequently modified, optionally joined, and functionalized. By using experimentally-inspired design rules for heterocycle construction, we also anticipate and then validate with DFT calculations that these ligands will form stable complexes. Nitrogen or oxygen atoms always act as the metalcoordinating atoms in eight distinct five-or six-membered heterocycles, and further modification (e.g., insertion of NH, S, O, or a C-C double bond) of the non-metal-coordinating ring structure produces 38 non-redundant heterocycles (Figure 1 and Figure S1). Along with the monodentate core ligands (including pyridine, furan, and oxazole), we fuse all possible heterocycle combinations into 741 bidentate ligands by forming a carbon-carbon bond between sites adjacent to the metal-coordinating atom (Figure 1 and Figure S1).</p><!><p>Approach for hierarchical assembly of a 2.8M homoleptic complex design space for multi-objective RFB design. The sequential stages that lead to each combinatorial step are indicated in the flowchart at left starting with selection and modification of rings, optionally fusing them to form 779 unique monodentate and bidentate ligands, functionalization of ligands with 897 unique functional groups, and finally complexation with four metals. At right, example skeleton structure components and their assembly or modification coinciding with each step is shown to form a bidentate oxygen-coordinating (coordinating atoms shown with red or blue highlight) ligand. The functional group assembly step is indicated in the black inset with resulting structure shown. Finally a ball and stick structure is shown for an assembled complex.</p><p>To ensure a smoothly varying space for property fine-tuning, core ligands were hierarchically functionalized 3 bonds away from the metal-coordinating atom (i.e., para in the six-membered rings, Text S1). Methylene or phenyl groups are sequentially added up to two times, optionally functionalized (i.e., for methylene only) with up to two functional groups (e.g., OH, NH 2 , or CH 3 ), and finally capped with a terminal functional group (e.g., Cl, OH, or NH 2 ) for a total of 897 functionalizations (Figure 1 and Figure S2 and Text S1). These functionalizations are added symmetrically to all ligands, producing approximately 700k distinct ligands that form 2.8M unique homoleptic complexes with four transition metals (Figure 1 and Text S1).</p><p>The performance of an RFB is determined by the cell voltage, E cell , which is affected by both the redox potential of the redox couple as well as its concentration in the RFB electrolyte:</p><p>where ΔG solv is the Gibbs free energy of the oxidation or reduction process, n is the number of electrons ionized, F is Faraday's constant, and C is the concentration of the redox active species in the electrolyte.</p><p>We thus focus on optimizing in our complex space the remaining two design objectives: i) maximizing redox potential and ii) maximizing solubility. As in prior work 9,11 , the ionization process we optimize is computed from the ground state of the solvated M(II) (M = Cr, Mn, Fe, or Co) complex to the lowest energy M(III) complex that differs by loss of a single spin-up or spindown electron (see Methods and Table S1). The ΔG solv of each redox couple is modeled as a oneelectron ionization of the M(II) complex:</p><p>evaluated as the gas phase adiabatic ionization potential combined with single point solvation free energy corrections (see Methods). To estimate solubility 48 in polar solvents, we compute the standard hydrophilicity or logP (i.e., partition coefficient between octanol and water) from the ground state M(II) complex as:</p><p>where the numerator and denominator both correspond to solvation free energies (see Methods).</p><p>To ensure computational efficiency, we model solvent effects implicitly with a thermodynamic cycle approach, neglect vibrational corrections, and combine hybrid DFT with a modest, polarized double-ζ basis set (Figure S3, see Methods). These choices also allow us to leverage data sets of transition metal complexes from prior work 7,[9][10]13 .</p><p>To identify potential lead compounds in the design space, we use the expected improvement (EI) criterion 49 , widely-employed [33][34][35][36][37] for optimization with probabilistic surrogate models. EI ranks data points based on how much they are likely to improve on the current best characterized point, emphasizing uncertain points based on the model used to make these predictions. For one dimensional optimization, the simplest acquisition function (i.e., a way to quantify the priority with which we should evaluate a point, x) corresponds to the expected decrease in the value of the objective, or improvement, I:</p><p>at a point x with estimated value ŷ(x) (i.e., from a model prediction) relative to the current best design, x * , with known value y(x * ). The expected improvement acquisition function is evaluated as:</p><p>where p x ( ŷ) is the distribution of model predictions. When the distribution is Gaussian with mean value µ(x) with variance σ 2 (x), one can obtain an analytical integral of eqn. ( 5) in an approach known as efficient global optimization 49 (EGO).</p><p>Within EGO, the EI acquisition function at every new point x is then:</p><p>where Φ and φ are the cumulative and distribution functions of the standard normal distribution.</p><p>The Φ term, also referred to as the probability of improvement (i.e., P[I]), encourages exploitation of the model, whereas the second term favors exploration of high uncertainty points.</p><p>While 1D EGO has proven useful in accelerating chemical design [33][34] , RFB redox couple design considerations are inherently multi-faceted. In this case, we must define a Pareto front of the best possible tradeoffs 50 between high redox potential (i.e., high ΔG solv ) and solubility (i.e., lower logP) in the design space (Figure 2). We therefore extend the two-dimensional EI [51][52] approach to chemical design (Algorithm S1). In 2D EI, the improvement is now defined as the Euclidean distance a point x lies beyond the current Pareto front, and the total probability mass beyond the front defines E[I] (Figure 2). For this approach, we need a ML-model predicted distribution of the ΔG solv and logP values for a new complex</p><p>where μΔG solv and μlogP are the predicted mean values and σ ΔG solv  To use 2D EI for accelerated exploration, we require a surrogate model capable of rapidly predicting properties that generalizes well across a 2.8 M compound design space as well as a quantitative uncertainty metric. A Gaussian process (GP) provides inherent uncertainty estimates from the structure of the learned posterior distribution 53 that are suitable for EGO, but ANNs do not. We recently developed a latent space distance uncertainty metric for ANNs 14 that both demonstrated superior performance to commonly employed (i.e., ensemble variance) estimates 29,54 in active learning experiments and also does not require the cost of training multiple models.</p><p>To convert the latent space distance, d, to the units of property prediction, we modeled the generalization error with a conditionally Gaussian distribution:</p><p>where σ 1 and σ 2 are parameters obtained from maximum likelihood estimation on a small (10%) set of out-of-sample complexes (Figure S4). Thus, our error distribution can be treated as Gaussian for both ANNs and GPs, reducing the choice of ML model in EGO to the one with the best generalization performance.</p><p>We therefore trained and compared GPs and single-task (ST) ANN models to independently predict DFT values of ΔG solv and logP as well as a multitask (MT) ANN trained to predict both ΔG solv and logP. All ML models are trained on revised autocorrelations (RACs) 9 , a molecular graph-based representation (see Methods). RACs are a representation tailored for open-shell transition metal chemistry 9,[11][12]16 , yielding modest errors (ca. 0.2-0.3 eV) in combination with kernel or neural network models on small data sets (ca. 200-300 points) of redox/ionization potentials 9,12 . Notably, this representation incorporates no geometric information, which enables rapid (ca. 3 minutes) ML-model assessment and uncertainty quantification of the 2.8M complex design space but precludes prediction of property conformational dependence (Text S2).</p><p>To compare the baseline prediction accuracy of GP and ANNs, we assembled a 'hot start' data set of 235 predominantly heteroleptic transition metal (M = Cr, Mn, Fe, or Co) complexes from prior work 7,[9][10]13 and assessed train/test errors from a 90/10 split (see Methods and Figure S5 and Table S2). For ΔG solv , training mean absolute errors (MAEs) are low across all models (ST: 0.11, MT: 0.15, and GP: 0.13 eV). The ST ANN has fewer outliers, as evidenced by a lower training root-mean square error (RMSE) than the other two models (Table S3). As expected, all test MAEs are higher, with the MT ANN performing the best (0.32 eV MAE), consistent with some prior observations on MT ANNs 55 (Figure 3 and Table S3). For logP, training MAEs are lowest for the GP (ST: 1.45x10 -4 , MT: 1.45x10 -4 , and GP: 1.05x10 -4 ), the ST ANN produces the lowest test error (1.49x10 -4 ), and all train/test MAEs are relatively modest (Tables S3  As a test of model generalization, we evaluate the performance of the models trained on 'hot start' data to predict properties of design space complexes. It is evident from principal component analysis in the RAC representation that the 'hot start' data and design space compounds have limited overlap (Figure 3). The 'hot start' data consists of smaller, lower symmetry complexes with low-denticity ligands as compared to the homoleptic complex space, which give them relatively high redox potentials but make them unsuitable for applications in RFBs (Figure 3 and Figure S8). To identify the model with the best generalization characteristics, we carried out diversity-oriented clustering on the complete design space and obtained DFT properties of a representative set of 107 transition metal complexes (see Methods).</p><p>As expected, MAEs are significantly higher for all models (2-3x for ΔG solv and 1-2 orders of magnitude for logP) on this representative set than the corresponding test MAEs (Figure 3 and Table S3 and Figure S6). Among the model choices, the MT ANN distinguishes itself by exhibiting the lowest generalization errors for both ΔG solv (MT: 0.71 eV, ST: 0.80 eV, GP: 1.47 eV) and logP (MT: 3.04x10 -3 , ST: 5.92x10 -3 , GP: 1.31x10 -2 ) properties (Figure 3 and Tables S3-</p><p>and Figure S6). While the representative data set properties differ significantly from the 'hot start' data for both ΔG solv and logP, poorer generalization of the logP models can likely be attributed to the more substantial differences in the logP distributions of the two data sets than ΔG solv (Figure 3). Both GP and ANN UQ metrics can be used to isolate higher error points, and all metrics identify higher prediction uncertainty for the representative members of the design space than for the 'hot start' data test partition (Figure S9). Due to its superior generalization, we select the MT ANN for multi-objective expected improvement optimization of redox couples.</p><p>At the start of the expected improvement optimization, we use the 107-complex representative data set to define our initial Pareto front (Figure 4). This Pareto front, i.e., the set of known compounds that represents the current best combination of ΔG solv and logP, consists of six points formed from three of the metals (i.e., Fe, Co, or Mn) complexed with bidentate ligands (Figure 4). An 82-atom, low-spin Fe complex with a nitrogen-coordinating bidentate ligand has the best (i.e., lowest) logP of -5.31x10 -2 (ΔG solv = 6.36 eV), likely due to the ligand's functionalization with a large number of strongly polar groups (i.e., chloride and carbonyl) (A in Figure 4). A 145-atom, high-spin Mn complex with an oxygen-coordinating bidentate ligand has the highest ΔG solv of 7.37 eV (logP = -4.38x10 -2 ) in the initial data set (B in Figure 4). The other four Pareto complexes naturally represent a tradeoff of the two quantities between these limits. We then retrain the MT ANN and calibrate its UQ metric on both 'hot start' and the representative data set and refer to this as the generation 1 ANN model (Figure 2 and Figure S10, also see Methods). Using the generation 1 model, we evaluate the probability that each candidate in the 2.8M complex design space will improve on the Pareto front (i.e., P[I]), as well as the associated expected improvement (E[I]). While the P[I] uniformly favors any point predicted by the ANN to be beyond the front, the E[I] focuses on the region of the front that will lead to the largest possible property improvement (Figure 4). In conventional EGO only the top-scoring candidate would be selected, but due to the parallel nature of our calculations we cluster the topscoring complexes to obtain representative points for characterization with DFT in order to generate additional data for model retraining (Figures 2 and 4, also see Methods). After only a single iteration, the top-scoring complexes making up the Pareto front indeed are primarily comprised of these new generation 2 complexes (Table S6).</p><p>Carrying out this optimization requires balancing computationally demanding DFT calculations at each generation (i.e., over a week of wall time for 100 complexes) with the much more rapid ANN-derived assessment of the full design space (i.e., 2.8M complexes in ca. 3 minutes, Text S2 and Table S7). While there is no universal standard for stopping EGO, we choose to stop after four generations, supported in part by the rapid decline of mean E[I] values of clustered points from the top (i.e., 10k of 2.8M) lead complexes from generation 1 (mean E[I] = 0.45) to 4 (mean E[I] = 0.04, Figure S11). We can also assess if model predictive power is improving during this process, given the small but growing size of the DFT training data with each generation (Table S8). To quantify model improvement, we calculated the 'lookahead' errors of models trained in one generation on data from later generations. Although the model trained only on 'hot start' data (i.e., generation 0) exhibits the poorest MAEs for both ΔG solv and logP on later generation (i.e., 3-5) data sets, model prediction accuracy for both quantities improves most significantly (> 50% reduction) with incorporation of generation 1 data (Figure 5 and Tables S9-S10 and Figure S12). These lookahead MAEs generally do not improve  Over each generation of EI, new points are added to extend beyond the previous generation's Pareto front. The final Pareto front consists of eight complexes, all from the last three generations of the EGO run (Figure 6 and Table S6). The direction of improvement evolves during EGO, with improvements to ΔG solv primarily occurring in generations 2 and 3 (Figure 6).</p><p>Conversely, the later generations 3-5 capture more logP improvements as ΔG solv increases diminish, with the six best (i.e., lowest) logP values occurring in these three generations (Figure 6). Although the DFT-evaluated Pareto front moves out at each generation, the improvements are indeed diminishing, consistent with a decrease in model-predicted E[I] values (Figure 6 and Figure S11). Broad chemical trends are evident in the evolution of the Pareto front during EGO. While the ligands in our design space support both high and low spin redox couples, only high-spin and predominantly Mn complexes are favored in later generations, consistent with some prior observations on smaller-scale, redox-focused computational screening 56 . One confounding effect in this process is that both time limits and calculation quality criteria prevent all attempted DFT calculations suggested by the model from succeeding and being incorporated into subsequent generations (Tables S7-S8). However, we note there is no metal-dependence in the success rates, meaning that the model acquires information both from 'hot start' data and from the generation 1 diversity-oriented points that inform the model's preference for high-spin Mn complexes (Figure S13). In fact, a high-spin iron complex with a bidentate ligand consisting of an oxathiin fused with 2,5-dihydrofuran and functionalized with halide capped phenyl groups also appears in the Pareto front at generation 5 because it has a lower logP (-5.60x10 -2 ) than any prior sampled complex (C in Figure 6). The highest ΔG solv complex consists of fused bidentate oxygencoordinating heterocycles, which each also contain a sulfur modification (i.e., 4H-1,3-oxathiin and 1,3-oxathiolan, structure A in Figure 6). This complex also exhibits a better logP (-4.46x10 -</p><p>2 ) than the highest ΔG solv complex in prior generations, likely due to the strongly polar hydroxymethyl functionalization (Figure 6). Assuming equal importance of improving both properties, one high-spin Mn complex with an unsaturated fused pyran and oxathiole structure along with ketone functional groups presents the highest redox potential (ΔG solv = 8.05 eV) while preserving good solubility (logP = -5.52x10 -2 ) at the Pareto front (complex B in Figure 6).</p><p>Commonalities are evident in the ligand structures of the Pareto front high spin Mn complexes with the highest ΔG solv values ranging from 8.05-8.51 eV (Figure 7). All contain oxygen-coordinating bidentate ligands with a combination of a six-and a five-membered ring, with at least one ring also incorporating a sulfur heteroatom (Figure 7 and Figures S14-S15). The EI process generally avoids selecting nitrogen-coordinating ligands or in-ring modification with nitrogen because these complexes are typically superseded by oxygen-coordinating structures with sulfur in-ring modifications (Figure 7 and Figure S1 and Table S6). A Pareto set complex after generation 3 with high ΔG solv (8.24 eV) that does have in-ring nitrogen-modification but also significantly lower logP (-4.61x10 -2 ) than other complexes (Table S6). As with metaldependence, all ligand types are sampled in early generations but the optimization quickly converges toward a few select heterocycles in later generations (Figure 7 and Figures S16-S17). The EI process most frequently selects Mn complexes that are uniformly smaller (ca. 90 atoms) than the initial diversity-oriented Pareto front complexes (ca. 120 atoms) in size (Table S6). This smaller complex size is consistent both with prior observations 9, 11 that smaller complexes have higher redox potentials and the fact that the highest logP iron complex noted earlier is both the largest (184 atoms) and lowest ΔG solv of all complexes in the final Pareto front.</p><p>Since high solubility appears to be heavily influenced by the presence and nature of polar functional groups, the EI scoring is likely favoring atom economy to maintain high redox potential while introducing groups that can improve solubility (Figure 7). The polar carbonyl oxygen and hydroxyl groups are most favored in the optimization along with chloride termination, especially on methylene group backbones, but a wide range of functionalizations are sampled during optimization (Figure 7 and Tables S6 and S11). Comparing two Pareto front complexes with the same core ligand, we observe that adding more polar functional groups (e.g., amine, carbonyl, and chloride vs. only carbonyl addition) does sometimes lead to significant logP enhancement (-4.49x10 -2 vs. -4.96x10 -2 ), albeit with a reduction in the ΔG solv (ca. 0.3 eV, Table S12). Overall, EGO using the EI criterion appears to have converged on key design elements to maximize redox potential while not limiting solubility.</p><p>To evaluate the efficiency and accuracy of the ANN-based EGO in the design space, we compare to observations made on a uniform random sample of the space (i.e., distinct from our initial diversity-oriented clustering, Figure 5, Table S8, and see Methods). Within EGO using the EI criterion, surrogate model errors could lead to poor prediction accuracy with simultaneous over-confidence causing specific areas of the design space to never be sampled with DFT calculations, as the surrogate model becomes myopically focused on the first initially promising region of chemical space. To assess this, we evaluated errors of each generation's ANN model on this 122-complex random set (Figure 5). After rapid reductions for both properties from the incorporation of generation 1 data (ΔG solv gen. 0: 0.70 eV, gen. 1: 0.41 eV, gen. 5: 0.38 eV), random set MAEs remained above test MAEs and comparable to or lower than generation 5 lookahead errors (Figure 5 and Tables S9-S10). Since the random set more evenly covers the space while the EI criterion has emphasized acquisition of points in one region of the design space, the low errors on the random set are encouraging evidence of the model maintaining similar accuracy across the full space to that obtained at the Pareto front. However, because this optimization procedure relies on a surrogate model, it is not guaranteed that the model-identified compounds will be optimal when evaluated with DFT.</p><p>To quantify optimization efficiency, we compare the distribution of properties observed in the 122-complex random set with the 8 complexes that form the Pareto front after five generations of EI (Figure 6 and Figures S18-S21, Tables S13-S15, and Text S3). The average of random set complex redox properties is significantly lower than the Pareto complexes, with a ΔG solv of 5.80 eV (Figure 6). The mean logP for the random set of -4.45x10 -2 also indicates poorer solubility characteristics. Of the 84 complexes characterized through E[I]-based selection, nearly one third ( 26) have a higher ΔG solv than the highest value observed in the random set, four have a lower logP than the lowest random set value, and none of the random set would lie on the generation 5 Pareto front. Design principles evident from the EI generations, such as large polargroup-functionalized phenyl rings for the lowest logP or smaller (ca. 100-atom) complexes for the highest ΔG solv , are also evident in the random set (Table S15).</p><p>Overall, the generation 5 Pareto set complexes correspond to over 4 standard deviations (std. dev.) beyond the mean of the random set for improved solubility characteristics (logP random set std. dev. = 2.75x10 -3 ) and 3 std. dev. above the mean for random set redox properties (ΔG solv random set std. dev. = 0.89 eV). Since the ΔG solv and logP random set properties are not normally distributed and are more concentrated at the mean, any distribution-based analysis will provide a conservative estimate of the EI approach's efficiency and performance (Figure S19).</p><p>Assuming the two properties are independent, this analysis suggests at least a 500-fold acceleration of EGO with 2D-EI over a random search, consistent with or exceeding prior applications of EGO in chemical discovery 30,34 (Figure S19). That is, a fully random search would be expected to require 500 times more DFT calculations to identify the same number of top performing compounds obtained with the EI approach. Given that EGO required around 4-6 weeks of parallel GPU-accelerated computation to acquire all of the DFT properties on these large transition metal complexes, this time favorably compares to approximately 50 years for a fully random search (Text S2).</p><p>Interestingly, one complex from the random set featuring a bidentate-nitrogen coordinating ligand in complex with high-spin Co had a better ΔG solv value (7.9 eV) than the best Co/N complex sampled during EGO (Tables S6 and S15). This suggests improvements to initial sampling of the space could improve the optimization strategy. We could also leverage recently developed ML models 57 that can predict when EI-requested complex calculations will succeed, increasing yield. Overall, EGO with 2D-EI provides a targeted method of orders-of-magnitude acceleration of the discovery of optimal complexes and design rules from a compound space intractable for direct simulation.</p><!><p>We have demonstrated a multi-objective materials design strategy applied to the accelerated discovery of candidate M(II)/M(III) redox couples for redox flow batteries that are both soluble in polar solvents and have high redox potentials. We enumerated a large 2.8M compound space comprised of bulky ligands designed to prevent crossover while forming stable transition metal complexes to identify redox couples that possessed both desirable solubility and redox characteristics. Since direct calculation of the full space with DFT would have been predictive design in the full compound space. Analysis of the property distribution of a subsequently collected random test set suggests that the ANN-driven EI approach achieved at least 500-fold acceleration over random search. In practice, this acceleration corresponds to EGO identifying a Pareto-optimal design in approximately five weeks as compared to a random search that would have instead required at least fifty years. A natural extension to this work is application to the design of earth-abundant, stable, and selective catalysts, for which multiobjective design is expected to be essential.</p><!><p>Data sets and calculation parameters. Initial structures selected from the design space for simulation were generated with molSimplify 58 , which assembles complexes with ANN-predicted metal-ligand bond lengths 7,11 and uses OpenBabel 59 as a backend for ligand structure generation from SMILES strings. All geometry optimizations and single point energy calculations were performed with a developer version of the TeraChem 60 code. For all DFT calculations, the B3LYP hybrid generalized-gradient approximation exchange-correlation functional was employed, which in TeraChem uses the VWN1-RPA form for the LDA component of the correlation energy. The LANL2DZ effective core potential was employed for transition metals and the 6-31G* basis set was used for the remaining atoms. We previously showed 11 this choice to introduce small errors in relative redox potential prediction with respect to larger basis sets.</p><p>Geometry optimizations were carried out with the TRIC 61 optimizer using default tolerances of 4.5x10 -4 hartree/bohr for the maximum gradient and 1x10 -6 hartree for the change in selfconsistent field (SCF) energy between steps. Level-shifting 62 values of 0.25 Ha for both virtual and occupied orbitals were applied to all non-singlet complexes, whereas singlets were computed in a restricted formalism.</p><p>The quantities selected for training and evaluation of each complex, ΔG solv and logP, were determined first by using the implicit-aqueous-solvent-corrected energies of the M(II) complex to identify its ground state as either low-or high-spin (LS or HS). LS and HS multiplicities for M(II)/M(III) complexes grouped by nominal electron configuration were as follows: quintet-singlet for both d 4 Mn(III)/Cr(II) and d 6 Co(III)/Fe(II), sextet-doublet for d 5</p><p>Fe(III)/Mn(II), and quartet-doublet for both d 3 Cr(III) and d 7 Co(II) (Table S1). Solvation energies for both logP (i.e., ΔG s,water and ΔG s,octanol ) and ΔG solv (i.e., ΔG s,water (M(II)) and ΔG s,water (M(III))) were obtained from single point energies with the conductor-like solvation model (COSMO) 63 implicit solvent model, as implemented 64 in TeraChem, for both octanol (ε = 10.3) and water (ε = 78.39) on all completed gas phase geometry optimizations. The solute cavity was built using defaults available for non-metals in TeraChem (i.e., 1.2x Bondi's van der Waals radii 65 ), and we provided standard van der Waals radii 66 for metals, which were also scaled by 1.2x. Although the target application in practical RFBs involves acetonitrile (ε = 38.8),</p><p>we choose water as our model polar solvent, given the similarity in the dielectric constants and the ability to use prior data.</p><p>Once the LS or HS ground state was identified, the oxidation process was computed as a one electron removal to a M(III) LS or HS complex. From this spin state, the adiabatic, gas phase ionization potential (IP), ΔE III-II , was computed:</p><p>from the electronic energies of the gas phase, geometry optimized complexes in the LS or HS ground state of the M(II) complex and the corresponding LS or HS M(III) complex. To obtain the approximate redox potential, ΔG solv , solvation free energies of the two end states are then added as described above, but we neglected vibrational corrections that we incorporated in prior work 9 because they are small in magnitude while significantly increasing the computational cost of screening (Figure S3). Aqueous and non-aqueous contributions to logP were also obtained as single point energies on gas-phase optimized ground state M(II) complexes.</p><p>Calculations were automatically submitted and monitored using molSimplify Automatic Design (mAD) [12][13] . The mAD workflow flags and excludes calculations that remain unconverged after the default maximum five 24 hour job resubmissions 12 . It also excludes unrestricted calculations with deviations of the expectation of the <S 2 > operator of more than 1.0 µ B from the expected value and checks structures for quality and preservation of expected connectivity 12 . For the bidentate ligands studied in this work, we have loosened some of the prior default 12 angular and root mean square deviation (RMSD) thresholds (Table S16). To obtain the target properties of interest, ΔG solv and logP, a minimum of three gas phase geometry optimizations 11 must pass all of these checks, i.e., M(II) complex optimizations in both spin states as well as the corresponding M(III) complex in the ground state spin, reducing overall success rates (Tables S7-S8 and Figures S20-S21).</p><p>The 'hot start' data 7,[9][10]13 was generated following the same workflow and spin/oxidation state definitions with only two differences: i) level shifting values for virtual and occupied states were 1.0 and 0.1 Ha, respectively, and ii) the L-BFGS implementation in Cartesian coordinates with the DL-FIND 67 was sometimes employed (Table S2). All optimized geometries and total energies are provided in the Supporting Information.</p><p>Machine learning (ML) models. We trained ML models on geometry-free, graph-based RACs 9 .</p><p>RACs are products and differences of heuristic properties (i.e., nuclear charge, topology, identity, Pauling electronegativity, and covalent radius) for atoms separated by a specific number of bond paths over an entire complex. Additional RACs are obtained by fixing the scope, i.e., by always incorporating the metal or ligand-coordinating atom. The maximum distance of 3 bonds for properties to be correlated has been previously motivated, leading to a total of 155 features for octahedral transition metal complexes, 153 from RACs 9 along with oxidation state and fraction of Hartree-Fock exchange, but the latter two are not relevant to this work.</p><p>ANNs were trained with hyperparameters selected using Hyperopt followed by manual fine-tuning in Keras with the Tensorflow backend. We hold out a uniform-random 10% of training data at each stage as a test set in order to calibrate the probabilistic ANN uncertainty model 14 (Figure S4). We independently optimized hyperparameters for 1500 iterations with 10% validation data for the three ANNs (two ST and one MT) for predicting ΔG solv and logP (Figure S22 and Table S17). Inputs (i.e., RAC feature vectors) and outputs (i.e., DFT properties) were normalized to make the scales of ΔG solv and logP consistent. The MT and ΔG solv ST ANN architectures were similar (i.e., 2-3 layers of 100 tanh nodes), whereas the ST logP ANN optimized to a smaller residual-type ANN with skip connections and two layers of 50 ReLu nodes (Figure S22 and Table S17). For ANN retraining, architectures were held fixed, and models were initialized with previously converged weights using a more aggressive learning rate decay and a larger batch size for smoother optimization (Table S17).</p><p>We trained GP regression models to independently predict ΔG solv and logP with a singleparameter isotopic Gaussian covariance kernel models using kernlab in R v.3.6.1.</p><p>Hyperparameters (i.e., kernel width and regulation strength) were optimized with a 2500-point logspaced grid search (10 -8 to 10 -1 ) using 10-fold CV, as implemented in CVST and reoptimized with data addition (Table S18 and Figure S23).</p><p>Diverse points for DFT simulation were selected with k-medoids clustering, as implemented in the Cluster package in R. All model weights are provided in the Supporting Information. In the first clustering of the complete design space, we obtained 300 medoids, which we confirmed to be representative of the space (Text S4, Figures S24-S25, and Table S19). Other clustering was carried out to obtain 100 medoids from top 10k 2D EI leads and employed the WeightedCluster package to weight points by expected improvement scores. These medoids of the top 10k leads are also representative: in the first generation, the selected 100 compounds have a comparable E[I] value (i.e., 0.46) to the top 10k candidates (i.e., 0.45) and much higher than the complete design space average (i.e., 0.04) assessed with the generation 1 model. In practice, only a subset of medoids succeed, either due to the time limits on the calculations or other checks (Tables S7-S8).</p>