High-Throughput Virtual Screening of Quinones for Aqueous Redox Flow Batteries: Status and Perspectives

: Quinones are one of the most promising and widely investigated classes of redox active materials for organic aqueous redox ﬂow batteries. However, quinone-based ﬂow batteries still lack the necessary performance in terms of metrics, such as speciﬁc capacity, power density, and long-term stability, to achieve mass market adoption. These performance metrics are directly related to the physicochemical properties of the quinone molecules, including their equilibrium redox potential, aqueous solubility, and chemical stability. Given the enormous chemical and conﬁgurational space of possible quinones and the high tunability of their properties, there has been a recent surge in the use of high-throughput virtual screening (HTVS) for the rational design and discovery of new high-performing molecules. In this review article, HTVS efforts for the computational design and discovery of quinones are reviewed with a special focus on the enumerated space of core quinone motif, the methods and approximations used for the estimation of performance descriptors, and the emergent structure-property relationships. The knowledge and methodological gaps in conventional HTVS efforts are discussed, and strategies for improvement are suggested


Introduction
The ever-increasing energy requirements caused by global population growth and economic growth have resulted in high consumption of fossil fuels, environmental damage, and subsequent climate change.Although the costs of renewable solar and wind energy have decreased over the last decade, their intermittent nature remains a significant issue, especially when combined with the fact that energy demand fluctuates both temporally and regionally.This has prompted the development of grid-scale energy storage systems to balance the fluctuating supply of energy and user demands.
Due to their scalability and distinctive design, which allows for the decoupling of their power and energy output [1][2][3], aqueous redox flow batteries (ARFBs) are among the most promising next-generation grid-scale energy storage devices.The chemical energy in ARFBs is provided by redox active materials dissolved in water and stored in external tanks, as shown schematically in Figure 1 below (adapted from ref. [4]), so a continuous circulation of the redox active materials using pumps is necessary to sustain the chemical reaction, similar to a fuel cell.Similar to a battery, the system can be charged and discharged because the related redox reactions are reversible.
The development of next-generation ARFBs with high cell voltage, energy density, power density, and long life depends on finding redox active materials (RAMs) with high redox potential, high aqueous solubility, high stability, and rapid redox kinetics [1][2][3].Satisfying these fundamental material requirements will form the foundation on which further cell-and pack-level optimization can be performed to bring ARFBs from the lab to the market.While vanadium ARFBs have taken the lead in commercialization among available technologies, they are still prohibitively expensive [1] and suffer from low rate kinetics and long-term stability issues [2].Inorganic alternatives based on iron [5] and zinc [6] offer potentially much lower costs; however, the increased complexity of these storage concepts (e.g., in hybrid flow batteries and metal-air systems) [7] due to the involvement of gas-phase active species and the propensity of metal deposition decreases the efficiency and operating lifetime significantly when compared to vanadium ARFBs.
Batteries 2023, 9, x FOR PEER REVIEW 2 of 28 further cell-and pack-level optimization can be performed to bring ARFBs from the lab to the market.While vanadium ARFBs have taken the lead in commercialization among available technologies, they are still prohibitively expensive [1] and suffer from low rate kinetics and long-term stability issues [2].Inorganic alternatives based on iron [5] and zinc [6] offer potentially much lower costs; however, the increased complexity of these storage concepts (e.g., in hybrid flow batteries and metal-air systems) [7] due to the involvement of gas-phase active species and the propensity of metal deposition decreases the efficiency and operating lifetime significantly when compared to vanadium ARFBs.
Figure 1.Schematic showing a typical redox flow battery with the tanks containing the dissolved catholyte (also known as posolyte) and the anolyte (also known as negolyte) materials.The electrolytes at each electrode are circulated using pumps.Reprinted with permission from Ref. [4].2022, Elsevier.
In this regard, organic molecules have emerged [8][9][10] as an attractive alternate class of RAMs for ARFBs because: (I) They are composed of only earth-abundant elements, which implies potentially much lower costs.(II) They allow a high degree of tunability of their physicochemical properties, which implies an immense potential to achieve high theoretical capacities, energy and power densities, and long-term stability.(III) They offer a huge amount of flexibility in the use of multielectron-transfer RAMs without the complications that are typically associated with inorganic multivalent ions [11].Various ARFB performance metrics, such as energy and power density, energy and charge capacity, efficiency (voltaic, coulombic), and capacity fade rate, are controlled by the physicochemical properties of organic RAMs, membranes, and electrodes.The physicochemical properties that influence these performance metrics are listed below in Table 1.A more detailed account of these relations can be found in recent review articles by Fisher et al. [8], Ding et al. [12], and Cao et al. [13].
Table 1.List of key physicochemical properties of redox active materials and the aqueous redox flow battery performance metrics that depend on them.Information obtained from Fisher et al. [8], Ding et al. [12], and Cao et al. [13].

Property
Battery Performance Metric redox potential energy and power density, cost aqueous solubility energy and power density, energy and charge capacity, cost rate kinetics power density, efficiency stability (chemical/electrochemical) efficiency, lifetime stability (calendar-life) lifetime, cost ionic conductivity power density, efficiency Schematic showing a typical redox flow battery with the tanks containing the dissolved catholyte (also known as posolyte) and the anolyte (also known as negolyte) materials.The electrolytes at each electrode are circulated using pumps.Reprinted with permission from Ref. [4].2022, Elsevier.
In this regard, organic molecules have emerged [8][9][10] as an attractive alternate class of RAMs for ARFBs because: (I) They are composed of only earth-abundant elements, which implies potentially much lower costs.(II) They allow a high degree of tunability of their physicochemical properties, which implies an immense potential to achieve high theoretical capacities, energy and power densities, and long-term stability.(III) They offer a huge amount of flexibility in the use of multielectron-transfer RAMs without the complications that are typically associated with inorganic multivalent ions [11].Various ARFB performance metrics, such as energy and power density, energy and charge capacity, efficiency (voltaic, coulombic), and capacity fade rate, are controlled by the physicochemical properties of organic RAMs, membranes, and electrodes.The physicochemical properties that influence these performance metrics are listed below in Table 1.A more detailed account of these relations can be found in recent review articles by Fisher et al. [8], Ding et al. [12], and Cao et al. [13].
Table 1.List of key physicochemical properties of redox active materials and the aqueous redox flow battery performance metrics that depend on them.Information obtained from Fisher et al. [8], Ding et al. [12], and Cao et al. [13].

Property
Battery Performance Metric redox potential energy and power density, cost aqueous solubility energy and power density, energy and charge capacity, cost rate kinetics power density, efficiency stability (chemical/electrochemical) efficiency, lifetime stability (calendar-life) lifetime, cost ionic conductivity power density, efficiency dynamic viscosity power density, efficiency permeation across the membrane energy and power density, efficiency, and lifetime toxicity environmental impact synthesizability cost The molecular structure of organic RAMs has a decisive influence on these properties, including redox potential, aqueous solubility, stability, electrochemical kinetics, synthesizability, and toxicity.This degree of freedom has been extensively utilized in both experimental and theoretical investigations of organic RAMs [4,10,11] by making systematic modifications to the molecular structures and analyzing their effects on ARFB performance.
Significant strides in organic ARFBs have been made in recent years, yet they must still undergo substantial research and development to achieve large-scale commercialization.There are several problems in employing organic RAMs in redox flow batteries, which include low energy density, low stability, insufficient rate kinetics, and high permeation through membranes, leading to irreversible capacity fade.Some of these problems are highly intrinsic to the molecular structure of organic RAMs and must be addressed at the molecular scale.According to the recent review by Fischer et al. [8], the first such problem with using organic RAMs is their comparably low energy density [8], which hardly exceeds 13 Wh/L compared to the state-of-the-art vanadium ARFBs (≈25-35 Wh/L) and zinc-bromine ARFBs (≈70 Wh/L).The reason for this has been pointed out to be the typically higher concentration of the inorganic RAMs and the battery voltage [4,8].The second problem with using organic RAMs (arguably far more critical in the long-term than the first problem) is their instability, which has several aspects ranging from chemical and electrochemical [14,15] to calendar-life [16] stability.For the majority of mechanisms that cause instability, the degradation of the organic molecule leads to irreversible capacity losses.In this regard, the reader is pointed to some excellent reviews published in the last few years, which cover a lot of ground in describing the underlying mechanisms of degradation [14 -16].
Among all the tested classes of organic molecules, quinones are one of the most widely investigated because they have shown great promise in ARFBs [33] due to the tuneability of their properties as both low-potential electrolytes (negolytes) and high-potential electrolytes (catholytes).Quinones are found in fungi, algae, bacteria, plants, and animals in various forms depending on the aromatic ring structure, size, sidechain (substituent) groups, etc.The chemical motif common to all quinones is the 6-membered ring system, in which carbonyl groups (-C=O) can be placed in para or ortho positions, depending on where the carbon atom belongs in the ring.
The chemical space of possible quinones (~10 6 ) offers a huge potential for the discovery and rational design of high-performing candidates that can eventually enable ARFB commercialization.Even the most advanced experimental screening techniques cannot screen such a large chemical space for the best candidates.However, the inherently parallel nature of the material screening problem coupled with recent advances in parallel computing and distribution techniques have given rise to the idea of high-throughput virtual screening (HTVS) [34].The central idea of HTVS is to screen out a large number of different materials enumerated in a virtual library as rapidly (and time efficient) as possible, with maximally autonomous workflows for screening the best candidates.HTVS involves the use of "funnels", such that each funnel is typically associated with a performance descriptor, which is related to one or more ARFB performance metrics.As screening proceeds, some candidate materials are eliminated at each level of the hierarchy in the funnel, depending on a descriptor screening criterion.The effectiveness of any HTVS effort in finding new material candidates that beat the state-of-the-art depends on a few crucial factors: (1) the size and diversity of the starting virtual library of candidates, (2) the accuracy and rapidity with which the performance descriptors can be calculated, and (3) the extraction of princi-ples for the rational design of new candidates that emerge from structure-property and property-property correlations.
The HTVS approach for quinones is shown schematically in Figure 2. It begins with the creation of a comprehensive virtual library of candidate molecules.This is typically performed by substituting hydrogen positions on core "backbone" rings with various substituent (or functional) groups and/or by substituting the carbon atoms in the rings with heterocyclic atoms.The candidates in the virtual library are then screened based on the values of their performance descriptors (redox potentials, aqueous solubility, chemical, and electrochemical stability, rate kinetics, synthetic accessibility, etc.), each of which is directly related to one or more ARFB performance metrics (volumetric capacity, energy and power density, capacity fade rate, etc.).The final set of screened molecules at the end of the funnel represents the "best" candidates that are worthy of further experimental validation.HTVS efforts for quinones for ARFBs [35][36][37][38][39][40][41][42][43] have not only led to the discovery of new high-performing candidates but also have elucidated several molecular design principles that are not intuitive, and which can be used to aid future rational design and discovery efforts.
performance descriptor, which is related to one or more ARFB performance metrics.As screening proceeds, some candidate materials are eliminated at each level of the hierarchy in the funnel, depending on a descriptor screening criterion.The effectiveness of any HTVS effort in finding new material candidates that beat the state-of-the-art depends on a few crucial factors: (1) the size and diversity of the starting virtual library of candidates, (2) the accuracy and rapidity with which the performance descriptors can be calculated, and (3) the extraction of principles for the rational design of new candidates that emerge from structure-property and property-property correlations.
The HTVS approach for quinones is shown schematically in Figure 2. It begins with the creation of a comprehensive virtual library of candidate molecules.This is typically performed by substituting hydrogen positions on core "backbone" rings with various substituent (or functional) groups and/or by substituting the carbon atoms in the rings with heterocyclic atoms.The candidates in the virtual library are then screened based on the values of their performance descriptors (redox potentials, aqueous solubility, chemical, and electrochemical stability, rate kinetics, synthetic accessibility, etc.), each of which is directly related to one or more ARFB performance metrics (volumetric capacity, energy and power density, capacity fade rate, etc.).The final set of screened molecules at the end of the funnel represents the "best" candidates that are worthy of further experimental validation.HTVS efforts for quinones for ARFBs [35][36][37][38][39][40][41][42][43] have not only led to the discovery of new high-performing candidates but also have elucidated several molecular design principles that are not intuitive, and which can be used to aid future rational design and discovery efforts.As is evident from the aforementioned literature references, there are already a huge number of review articles on the performance and status of organic RAMs.In this article, HTVS efforts for the computational design and discovery of quinones are reviewed with a special focus on the enumerated space of core quinone motifs, the methods used for the estimation of performance descriptors, and the emergent structure-property relationships.The main objective of this article is to inform researchers who are new to this field about the choice of available methods for quinone HTVS and their effectiveness, and to As is evident from the aforementioned literature references, there are already a huge number of review articles on the performance and status of organic RAMs.In this article, HTVS efforts for the computational design and discovery of quinones are reviewed with a special focus on the enumerated space of core quinone motifs, the methods used for the estimation of performance descriptors, and the emergent structure-property relationships.The main objective of this article is to inform researchers who are new to this field about the choice of available methods for quinone HTVS and their effectiveness, and to point out specific ways in which the overall HTVS strategy needs to be improved.Another objective is to summarize the discovered structure-property relationships and design principles for researchers who may use them to optimize the performance of quinones in experiments.This article does not address the general topic of high-throughput and machine learning-based prediction of redox potentials, aqueous solubility, or other physicochemical properties of organic compounds, which can be found elsewhere [44][45][46][47][48][49][50][51].As will be evident in the following, a majority of current HTVS efforts for quinones focus on the first problem of low nominal voltage and concentration of quinone ARFBs.In this regard, the applied methods for theoretical predictions of the redox potential and aqueous solubility are reviewed in greater detail.A list of all the abbreviations and mathematical symbols used in this article is given in Table 2.

Redox Potential
The cell voltage is a metric of utmost importance as it is directly proportional to the energy density and energy capacity of the ARFB [8,12,13].The cell voltage is related to the equilibrium redox potential of the RAM.From a theoretical standpoint, the redox potential is arguably the most extensively investigated property of quinones.This is in part due to their versatility as the RAM in both aqueous [52][53][54][55][56][57][58] and nonaqueous [59][60][61][62][63][64] flow batteries, as well as the electrode [65] in organic Li-ion batteries.The redox potential of quinones can be tuned by structural functionalization, i.e., by substituting electron-withdrawing groups (EWGs) and electron-donating groups (EDGs) into the core structure.In general, adding electron-donating groups shifts the redox potential in the negative direction, and vice versa [21].
Another decisive factor that influences the redox potential is the electrolyte's pH value during operation.Depending on the pH, the production or transport of either hydronium (H 3 O + ) or hydroxy (OH − ) ions can change their concentrations by a large amount.As a result, the equilibrium voltage might shift significantly during the redox reaction, or charge accumulation can occur, which can change the redox chemistry in an undesirable way [66].
Theoretical predictions: Redox potentials are defined in the context of reversible redox reactions at equilibrium and are related to the Gibbs free energy of change of the reaction through the Nernst Equation [67].The redox reactions of quinones proceed via the change in the formal oxidation state of oxygen atoms in the carbonyl groups that get protonated.Depending on the pH of the solution and the type of solvent, quinones can undergo either a two-e − reduction, a two-e − /one-H + (aq) reduction, or a two-e − /two-H + (aq) reduction, which is also called proton-coupled-electron-transfer (PCET) mechanism [67].For purposes of computational chemistry-based predictions, however, the PCET mechanism for quinones is widely accepted in the literature, especially in acidic aqueous environments.For the reduction of a solvated quinone molecule Q (aq) to QH upon transfer of m electrons (e − ) and n protons (H + (aq) ), the redox reaction is given as: The Gibbs free energy change in the reaction, ∆ r G, at temperature, T, can be expressed as: where the ∆ 0 f G terms represent the standard (Gibbs) free energies of formation, the terms in square brackets represent activities, F is the Faraday constant, and the free energy of the electron at the electrode, which is at equilibrium redox potential E r , has been expressed in terms of its electrochemical potential −FE r .At equilibrium, ∆ r G = 0, therefore: This general form of the redox potential does not take into account any kinetic phenomena or the order of e − /H + transfer.It assumes that the reduced and oxidized forms are thermodynamically stable, and do not involve radicals and open-shell calculations.In addition, this form assumes a single-step transfer of all e − and H + .In a more general, multi-step scheme, the second term in Equation ( 3) is often reported in terms of the acid dissociation constants (pK a ) of the reaction intermediates [37,56,57].The general e − /H + (aq) transfer scheme (adopted from [57]) is depicted in Figure 3a, which also shows a typical Pourbaix diagram (Figure 3b) of the quinone redox potential as a function of the pH.
Batteries 2023, 9, x FOR PEER REVIEW 7 of 28 Using techniques from computational chemistry and statistical thermodynamics, it is possible to theoretically estimate the standard free energy change during the redox reaction ∆ G.Many models are now available in the literature that allow for quantitative predictions of reduction potentials [67] by treating the effect of the solvent either implicitly (as a continuum dielectric) or explicitly (e.g., by including water molecules in the solvation shell).In practice, the calculated redox potentials are often reported not in absolute terms but rather with reference to a known standard reference electrode.An additional benefit of using a reference is that it enables the use of isodesmic reaction schemes.As the experimental value of the reference electrode is known, the difference between the experimental and calculated values for the reference species yields correction terms, which can be then added to the predicted value of the molecule of interest.For the case of quinones in

At standard conditions ( Q
= 1, pH = 0), this general form reduces to the well-known form of the Nernst equation: Using techniques from computational chemistry and statistical thermodynamics, it is possible to theoretically estimate the standard free energy change during the redox reaction ∆ 0 r G.Many models are now available in the literature that allow for quantitative predictions of reduction potentials [67] by treating the effect of the solvent either implicitly (as a continuum dielectric) or explicitly (e.g., by including water molecules in the solvation shell).
In practice, the calculated redox potentials are often reported not in absolute terms but rather with reference to a known standard reference electrode.An additional benefit of using a reference is that it enables the use of isodesmic reaction schemes.As the experimental value of the reference electrode is known, the difference between the experimental and calculated values for the reference species yields correction terms, which can be then added to the predicted value of the molecule of interest.For the case of quinones in ARFBs, the standard hydrogen electrode (SHE: 2H + (aq) + 2e − H 2(g) ) is the most popular choice for the standard reference electrode.The standard redox potential relative to SHE, E 0 r(SHE) , can be understood in a straightforward manner as follows: For a typical two-e − /two-H + (aq) (PCET) reduction in quinones in water, m = n = 2. Therefore, Equation ( 5) takes a form that lends itself very suitably to calculations, given by: where the explicit calculation of ∆ 0 f G H + (aq) is avoided and all the other free energies of formation can be calculated using any one of the numerous modern computational chemistry codes.Nevertheless, predicting the free energies of the formation of the molecules in aqueous media is not entirely free of error.Commonly used computational chemistry models have well-known limitations of either ignoring or not accurately accounting for contributions from explicit solvation and reorganization, entropic, zero-point, and nonelectrostatic effects.Often in practice, either the solvation (electrostatic, reorganization, cavitation, etc.) or thermal (entropic, zero-point) contributions or both are added separately as correction terms (∆ c G) to the internal energy (E DFT i ), both of which can be calculated using DFT.In this way, the Gibbs formation energy terms are approximated as This results in a less accurate yet highly effective way to rapidly estimate redox potentials given by: This formulation works rather well in practice because of the serendipitous cancellation of errors in the correction terms in Equation ( 7).As will be discussed in the following, various studies have assigned different degrees of importance to the factors influencing the correction terms.While some ignore it entirely, others take into account both the solvation and the thermal contributions.However, there exist disagreements in the literature about the relative importance of these contributions.A natural strategy that has been adopted to overcome these limitations, while maintaining predictability, is to use experimental redox potential data for fitting the parameters of a simplified linear equation: where m and c are simply the slope and intercept of the linear calibration model.Ideally, m and c should be determined by fitting experimental data obtained under consistent conditions with respect to the type of electrode, supporting electrolyte, concentration, temperature, scan rate, etc. Lastly, an even simpler set of descriptors are the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) energy levels for direct prediction of the redox potentials [68].From a mechanistic standpoint, it is clear to understand that the reduction in Q indicates e − transfer from the Fermi level at the electrode to its LUMO level, and the oxidation in QH 2 indicates e − transfer from its HOMO level to the Fermi level at the electrode.The limiting Fermi level is determined by the redox potential.Similar to Equation ( 8), the HOMO/LUMO levels can be approximately related to the redox potential as: where the constants are obtained by fitting experimental data.Irrespective of the chosen level of approximation, the inherent accuracy of predicting redox potential is significantly affected by the chosen method for calculating internal and free energies.Some theoretically rigorous approaches for predicting free energies of the reactants and products in the aqueous phase are: (a) performing molecular dynamics (MD) simulations over a long enough simulation time, followed by thermodynamic integration [52,54] to obtain ∆ 0 f G, or (b) sampling the system over a large number of microstates and performing ensemble averaging to predict the ∆ 0 f G [53].The accuracy of this approach is also affected by the choice of the interatomic potentials for energy calculations and the algorithms used for sampling.However, by design, this approach includes the contributions of explicit solvation, entropic, nonelectrostatic, and ion-pairing effects [69] without necessarily resorting to fitting or heuristics.
Liu et al. [52] predicted the redox potential of three quinone molecules using firstprinciples MD and then performed an integration of the vertical energy gap [70,71].The first-principles MD simulations were performed using the BLYP exchange-correlation functional and the computationally expensive hybrid HSE06 functional.They found that the BLYP-driven MD simulations produced much higher errors (0.2-0.7 V) compared to the HSE06-driven MD simulations, which were an order of magnitude lower (0.02-0.09V) with respect to experiments.It must be noted that their production simulations were run for only a small timespan of 6.0 ps after equilibration for 2.0 ps, which was likely due to computing time limitations associated with first-principles MD.While the appropriate MD simulation time depends on the characteristic relaxation time for a given system, 6.0 ps is not a large time window to get robust statistics, which could explain these relatively high errors despite using first-principles MD.
Kim et al. [53] created one of the earliest attempts to include explicit solvation effects in the prediction of redox potentials of quinones using quantum mechanics/molecular mechanics (QM/MM) simulations and compared their predictions to those of DFT calculations with the polarizable continuum model (PCM) of implicit solvation.They found that the DFT (PCM) model predicts the redox potentials of 16 out of 19 anthraquinones with a low mean absolute deviation (MAD) of 0.037 V.However, it is highly erroneous (MAD of 0.194 V) for three redox couples that have either group positioned for intramolecular hydrogen bonding that are not balanced with the intermolecular H-bonding of the solvent.The QM/MM simulations, which include explicit water molecules, were found to account for these imbalances, thereby reducing the MAD to 0.033 V.The best linear fit between their theoretical predictions and experiments was found to have an R 2 = 0.835.This accuracy has been improved by later prediction models that do not use computationally expensive QM/MM simulations.
A more recent study by Yu et al. [54] also used first-principles MD followed by thermodynamic integration to calculate the redox potentials of quinones and their derivatives.They used a multistep scheme in which the MD trajectories were first generated using the relatively cheaper functional BLYP, and then the configurations from the trajectories were computed with higher-level functional HSE06 and Grimme's D3 van der Waals corrections for subsequent use in thermodynamic integration.Their scheme yielded quite a low MAD (0.028 V) for 12 quinones.They also used a gas phase model akin to Equation (7) and fit the constants of the linear equation to their first-principles MD results and found an overall MAD of 0.097 V, which is over three times that of the MD results versus the experiments.Nevertheless, this deviation is within typical DFT errors (∼0.1 V), which explains very well why the linear model has been used extensively and effectively in the literature by fitting it to experimental data.Interestingly, the authors also used only 5 ps of MD sampling time, but with added van der Waals corrections, higher energy cutoffs, and better thermostats compared to that of Liu et al. [52].The first-principles MD approach has not found much traction in practice, not only because it is several orders of magnitude more computationally expensive than other approaches but also because simple linear models are already quite effective within the range of acceptable errors, as is discussed in the following.
Zhang et al. [55] performed a systematic comparison of several computational chemistry methods to point out strategies that balance the trade-offs between prediction accuracy and computational cost.They employed nine (AM1, MNDO, MNDOD, PM3, PM6, PM6-D3, PM6-D3H4X, PM7, and RM1) semiempirical quantum mechanics (SEQM) methods, two (GFN1-xTB, DFTB-D3) density functional-based tight binding (DFTB) methods, and seven (LDA, PBE, BLYP, B3LYP, PBE0, HSE06, and M08-HX) DFT functionals for performing calculations as various levels of approximations corresponding to Equations ( 7) and ( 8), as also shown schematically in Figure 4a.They performed geometry optimization and predicted formation energies at different levels of theory (SEQM, DFTB, and DFT) both with and without the effect of implicit solvation.By fitting the coefficients of Equation (8) to experimental data on 43 quinones that have been studied experimentally in ARFBs, Zhang et al. [55] inferred the following: (I) the inclusion of implicit solvation effects within the Poisson Boltzmann Formalism (PBF) reduced the root mean squared error (RMSE) across all levels of theory by ∼0.02 V, (II) while (I) is hardly surprising, it was also found that molecular geometry optimization with implicit solvation consistently leads to higher errors and a better strategy is needed to perform geometry optimization in the gas phase followed by single point formation energy calculations of ∆ 0 f G Q (aq) /∆ 0 f G QH 2(aq) with implicit solvation, (III) there is no significant effect of using zero-point and or dispersion corrections on the prediction accuracy at any DFT level of theory.Their work showed that at the DFT level, the PBE functional, which has little to no empiricism in its design, offers the best compromise between prediction accuracy (RMSE = 0.051 V) and computational cost.Notably, the computationally more demanding functional M08-HX performs best (RMSE = 0.044 V) among the selected functionals, as it is heavily parametrized to show good performance for thermochemistry.Most importantly, a comparison across all methods and approaches showed that DFTB-D3, which is at least three orders of magnitude computationally cheaper than DFT (PBE), can offer reasonable accuracy (RMSE = 0.072 V) and is most suitable for HTVS.
A recent effort by Fornari et al. [56] also performed extensive analyses of the various factors affecting the accuracy of redox potential predictions, with a special focus on the treatment of the deprotonated state of quinones at high pH.They chose B3LYP as their standard method for the calculation of formation energies along with the COSMO implicit solvation method, and their results had both similarities and dissimilarities with previous reports.Notably, they found that B3LYP performed even better than the composite G4MP2, the double hybrid rev-DOD-BLYP-D4, and the meta-hybrid M06-2X-D3 functionals by a significant amount (∼0.1 V).They compared COSMO-based predictions to the predictions from the SM12 solvation model and found only systematic differences in the order of 0.1 V, which could be removed by recalibrating the linear fit.To understand the deficiencies across implicit solvation models, they employed a cluster continuum approach with explicit water molecules to account for solvation.Using this approach, they were able to reduce the error for an outlier molecule from 0.43 V to as much as 0.15 V by adding 18 explicit water molecules.These findings are in agreement with the previous investigation by Kim et al. [53], who found that treatment of aqueous solvation using explicit water molecules leads to a significant reduction in prediction errors.In overall agreement with Zhang et al. [55], they also found that the inclusion of solvation effects lowers the prediction errors significantly and that geometry optimization with a solvation model provides only slightly better results than in the gas phase.However, they observed that it is the addition of the thermochemical contributions that causes the largest reduction in the mean absolute errors from (MAE) 0.308 V to 0.074 V.It must be noted that the prediction models from Fornari et al. [56] and Zhang et al. [55] cannot be directly compared because the latter used experimental data to fit the linear model's (Equation ( 8)) slope and intercept.Nevertheless, from the survey of these studies [52][53][54][55][56], it is quite clear that using less parametrized DFT functionals (PBE and B3LYP) can account for most of the contributions to the free energy changes during electron transfer that remain after error cancellation.The use of expensive higher-order DFT functionals offers no significant improvement in accuracy when predicting quinone redox potentials.The problem of predicting redox potentials across the entire pH scale was treated rigorously by Gaudin et al. [57] As shown in Figure 4b, they evaluated the chemical potentials of reaction intermediates (Q, Q , QH , and QH ) with different oxidation and protonation states (see also Figure 3) using the COSMO model, which yielded Gibbs free energies in an ideal conductor.These were then used as inputs in a COSMO-RS model to obtain their Gibbs free energies in water, which, summed together with the Gibbs free energy in the conductor, provided the total free energy of the reaction intermediates.These Gibbs free energies were plugged into energy balance equations akin to Equation (7), which provided the equilibrium redox potentials for Q ⇌ Q (at pH 14) and Q ⇌ QH (at pH 0).Furthermore, the limiting pK values for QH ⇌ QH and QH ⇌ Q were calculated using the same COSMO-based geometry optimization and solvation models (as for the Gibbs free energies).Furthermore, they used material balance equations, redox potentials, and pK values of the intermediates to develop a self-consistent model for finding the number of transferred protons at any given pH.With the help of this formalism, they were able to calculate Pourbaix diagrams that were in good quantitative and qualitative agreement with reported data on four different anthraquinones.Expectedly, the predictive power of this approach is dependent on the accuracy with which the Gibbs free energies of formation and pK can be predicted.The prediction of pK using first-principle simulations is a wide-open research area in itself, and some excellent reviews have already been published [72][73][74].

Aqueous Solubility
The aqueous solubility of quinones is the next important property that has been widely investigated because it is directly proportional to the energy capacity, energy To include the effect of pH on redox potential predictions, especially at higher pH values, Fornari et al. [56] used the Alberty-Legendre transform, which is similar to Equation (3).Algebraically, the use of Equation ( 3) only comes into play only when (under the assumption of simultaneous e − /H + (aq) transfers) the number of transferred e − and H + (aq) are not the same.However, at higher pH both the oxidized and the reduced form undergo deprotonation.In this case, the free energy terms in Equation ( 3) must be updated to account for this change in state.Fornari et al. [56] used the acid dissociation constants (pK a ) of the reactant and product molecules predicted from the ChemAxon cheminformatics software to account for deprotonations of both redox forms at higher pH.With this new protocol, which combined both DFT calculations and cheminformatics-based predictions, the MAE with respect to experimental potentials was found to be 0.121 V at pH = 7 and 0.144 V at pH = 13.
The problem of predicting redox potentials across the entire pH scale was treated rigorously by Gaudin et al. [57] As shown in Figure 4b, they evaluated the chemical potentials of reaction intermediates (Q, Q 2− , QH − , and QH 2 ) with different oxidation and protonation states (see also Figure 3) using the COSMO model, which yielded Gibbs free energies in an ideal conductor.These were then used as inputs in a COSMO-RS model to obtain their Gibbs free energies in water, which, summed together with the Gibbs free energy in the conductor, provided the total free energy of the reaction intermediates.These Gibbs free energies were plugged into energy balance equations akin to Equation (7), which provided the equilibrium redox potentials for Q Q 2− (at pH = 14) and Q QH 2 (at pH = 0).Furthermore, the limiting pK a values for QH 2 QH − and QH − Q 2− were calculated using the same COSMO-based geometry optimization and solvation models (as for the Gibbs free energies).Furthermore, they used material balance equations, redox potentials, and pK a values of the intermediates to develop a self-consistent model for finding the number of transferred protons at any given pH.With the help of this formalism, they were able to calculate Pourbaix diagrams that were in good quantitative and qualitative agreement with reported data on four different anthraquinones.Expectedly, the predictive power of this approach is dependent on the accuracy with which the Gibbs free energies of formation and pK a can be predicted.The prediction of pK a using first-principle simulations is a wide-open research area in itself, and some excellent reviews have already been published [72][73][74].

Aqueous Solubility
The aqueous solubility of quinones is the next important property that has been widely investigated because it is directly proportional to the energy capacity, energy density, and power density of the ARFB.The amount of quinone that can be dissolved in water is decided by the thermodynamics of aqueous solvation, which competes with the thermodynamics of crystallization.The aqueous solvation of quinones, and organic molecules in general, is a complex function of polar electrostatic interactions with water, solvent reorganization and entropic contributions, inter-and intramolecular H-bonding, and the delocalized charge density in the aromatic rings.Water is a polar solvent with a high Gutmann donor number, a high acceptor number, and a high polarizability.As a consequence, the presence of substituent groups and the size of the quinone molecule have an enormous influence on its solubility.Hydrophilic substituents, including −SO 3 H, −PO 3 H 2 , −COOH, and −OH, can increase the solubility of the core quinone by a few orders of magnitude [21,[75][76][77].In general, negatively charged substituents (−SO 3 − , −PO 3 H − , −COO − , etc.) have the most enhancing effect on solubility, followed by positively charged substituents (−C 5 H 6 N + , −C 3 H 5 N 2 + , −N(CH 3 ) 4 + , etc.) and then neutral substituents (−CN, −OH, −OCH 3 , etc.) [8].Apart from the molecule itself, the pH of the solution and the accompanying counter-ions in the solution have an enormous effect on the solubility [76,77].However, these effects have not been taken into account when making theoretical predictions of quinone solubility, as is evident in the published literature.
Theoretical predictions: Intrinsic aqueous solubility is defined as the maximum amount of a nonionized compound that can get dissolved in 1 L of water and is commonly measured in LogS (S in mol/L) units [78].Predictions of aqueous solubility are often very inaccurate, even though the solubility of a compound can be related rigorously to its other structural and physicochemical properties [79][80][81].This decades-old challenge has been tackled with broadly four types of approaches: (1) empirical methods, such as the generalized solubility equation (GSE) and its variants that use property-property relationships [82], (2) statistical methods that use quantitative structure-property relationships (QSPR) and cheminformatic descriptors [83,84], (3) physics-based methods that use Monte Carlo, MD or first-principles simulations to predict the free energy of the dissolution reaction [85,86] and (4) data-driven, machine learning, and other high-dimensional regression methods [47,87,88].
Multiple efforts have attempted to combine the strengths of various approaches, for instance, by combining first-principles calculated energies with physicochemical and cheminformatic descriptors as inputs in a data-driven machine learning model [84,89].However, it has proven very difficult to achieve error values below 0.3 LogS units on datasets of appreciable size (∼1000 points), which indicates a huge real error of 10 0.3 ≈ 2 times in solubility.Remarkably, current models do not improve even when trained on standardized data with much lower measurement RMSE values [90] (∼0.05 LogS) than usually available [91,92] (∼0.6-0.7 LogS).This suggests that there is little complementarity between the chemical information provided by the available methods, and more importantly, that an accurate method that captures the underlying thermochemistry of dissolution is still missing.Excellent reviews on these issues are already available in the published literature [79,80].
As will be evident in the next section, most HTVS studies ignore the important effect of coordinating ions [76,77,93] and pH [75] when predicting the aqueous solubility of quinones using first-principles methods.The solubility S o (in mol/L) can be predicted from the Gibbs free energy of the dissolution of the organic molecule from its crystalline form to its solvated aqueous form.At equilibrium, the free energy balance is given as: where ∆G * f Q (cry) and ∆G * f Q (aq) represent the Gibbs free energies of formation for the crystalline and aqueous phases, Q (aq) represents the activity of the molecule at equilibrium, and the * indicates the Ben-Naim standard states [86] (i.e., at 1 mol/L).The free energy balances for dissolution are depicted schematically in Figure 5.The activity Q (aq) can be expressed simply as the product of solubility, S o , and molar volume, V m .The definition of the dissolution Gibbs free energy is . Substituting these in Equation ( 10) yields: This equation is valid for the dissolution of a molecule Q in neutral water, without considering its ionic dissociation (pK a ) or acidity of the solution (pH).In practice, one uses a thermodynamic cycle for dissolution in which the crystalline phase first gets sublimated to the gas phase, followed by aqueous solvation of the gaseous molecule, as shown in Figure 5 below.

S V exp RT
This equation is valid for the dissolution of a molecule Q in neutral water, without considering its ionic dissociation (pK ) or acidity of the solution (pH).In practice, one uses a thermodynamic cycle for dissolution in which the crystalline phase first gets sublimated to the gas phase, followed by aqueous solvation of the gaseous molecule, as shown in Figure 5 below.One of the reasons for using such a cycle is that the Gibbs free energies for sublimation ( ) are more straightforward to calculate than those for dissolution.This thermodynamic cycle can take two paths, via the Ben-Naim (i.e., at 1 mol/L) or the standard (i.e., gas at 1 atm) states.As the Gibbs free energies for sublimation are typically referred to as the gas at standard conditions, one can rewrite ΔG * Q as: where ΔG Q is the standard Gibbs free energy of sublimation, p 1 atm, and T 298.15 K. Combining Equations ( 11) and ( 12) yields: As is evident from Equation ( 13), the accurate first-principles calculations of solubility necessitate full computational characterization of both the crystal and solution phases.While the prediction of the solvation free energy (ΔG * Q ) has matured significantly over One of the reasons for using such a cycle is that the Gibbs free energies for sublimation ) are more straightforward to calculate than those for dissolution.This thermodynamic cycle can take two paths, via the Ben-Naim (i.e., at 1 mol/L) or the standard (i.e., gas at 1 atm) states.As the Gibbs free energies for sublimation are typically referred to as the gas at standard conditions, one can rewrite ∆G * sub (Q) as: where ∆G o sub (Q) is the standard Gibbs free energy of sublimation, p o = 1 atm, and T = 298.15K. Combining Equations ( 11) and ( 12) yields: As is evident from Equation ( 13), the accurate first-principles calculations of solubility necessitate full computational characterization of both the crystal and solution phases.While the prediction of the solvation free energy (∆G * sol (Q)) has matured significantly over the last decade [81,86], there is ample evidence that the most erroneous contribution to the calculation of the dissolution free energy originates in the incorrect estimation of the sublimation free energy (∆G o sub (Q)), which relates to the solute's crystal structure [94,95].The algorithms for crystal structure prediction of organic crystals typically require extensive searches over space groups, unit cell dimensions, and atomic coordinates using a large number of randomized simulations in a Monte Carlo scheme [96].Such an approach for predicting the sublimation-free energy is computationally prohibitive and does not even guarantee high prediction accuracy.This is simply because of the exponential scaling in Equation ( 13) such that an error of 1 kcal/mol (which is the "heaven of chemical accuracy" for thermochemical calculations) in the prediction of the sublimation free energy implies an error of ∼6.4 times in the prediction of solubility in mol/L.Therefore, it is hardly surprising that all known HTVS efforts to date, which use firstprinciples simulations for predicting quinones' solubility, routinely ignore the sublimation free energy and use the free energy of solvation (∆G * sol ) as a descriptor for the overall solubility.Using only ∆G * sol is not entirely without merit because it is loosely correlated with experimental and cheminformatics-predicted solubility [36].Unlike the extensive work on redox potentials, there are no known studies that have systematically tried to identify rapidly calculable descriptors [97] for predicting the aqueous solubility of quinones exclusively, either using first-principles simulations or other methods.The standard method for calculation of solvation-free energies using first-principles methods is to use an implicit solvation method.Implicit solvation is a coarse-grained yet highly effective approach in atomistic simulations to account for a surrounding liquid electrolyte on the level of a continuous polarizable medium, for which the readers are pointed to an extensive recent review in ref. [98].However, this approach is not free of errors and is often corrected by using explicit water molecules as conducted in the cluster-continuum approach [56,99], or the use of QM/MM approaches [53,100].
The use of cheminformatics (e.g., in the software ChemAxon) and statistics-based methods is quite prevalent in the literature for HTVS of quinones based on solubility, as will be evident in the next section.These methods are dependent on simple and rapidly calculable descriptors and physicochemical properties of molecules (in their gaseous or solvated phases).In principle, such methods can be systematically tuned by fitting them to more available experimental data on aqueous solubility.It must be noted that the currently available amount of open-source data on aqueous solubility is neither large nor diverse enough [101] for building data-driven models that are applicable across all classes of molecules and that can predict solubility within ∼0.3 LogS units.The best models to date use careful data selection to maximize accuracy [87,89], but this strategy offers no guarantees on the accuracy of molecules "out of set".The availability of consistently produced experimental data specifically for quinones is even more scarce and scattered, with only a few tens [21,[75][76][77] of available data points.

High-Throughput Virtual Screening: Status and Challenges
For organic ARFBs, a majority of HTVS efforts have been invested in discovering new quinones as negolytes, while discovering quinones as posolytes has served as a much greater challenge to achieve.A summary of the HTVS efforts for quinone discovery from published literature is presented in Table 3.The table lists the core quinone structures that were used as backbones for functionalization along with the literature reference under the column name "Total cores/[Ref]", the substituent groups are listed under "Substituents", and the total number of redox couples enumerated in the library under "Library".It must be noted that some of these reviewed studies are simply high-throughput efforts without the screening aspect.
The HTVS of quinones begins with the creation of a virtual library in which the candidate quinones are systematically enumerated.As shown schematically in Figure 6, the enumeration can be conducted along various structural and configurational degrees of freedom, which can be broadly classified into two categories under core variations (Figure 6a) and substituent variations (Figure 6b).Core variations indicate modifications to the core molecules based on the extent and type of conjugation in the structure, the total number and position of the carbonyl groups, and the type of heteroatom substitutions.Analogously, substituent variations indicate modifications based on their type and mix of substituents, their positions relative to each other and relative to the carbonyl groups, and their total number and tendency to form intramolecular H-bonds.One of the seminal efforts to use HTVS to discover benzo-, naphtho-, and anthraquinones for ARFBs was performed by Er et al. [35].They considered two limiting conditions in the design of their virtual library for understanding the effects of the incorporated substituents on the redox potential of the quinone couples, i.e., with a single position substituted and with all positions substituted.Based on their analysis of the structure-property relationships, Er et al. [35] found that groups, which enable new intramolecular hydrogen bonds due to their positioning, strongly influence the calculated total free energies of the molecules.The quinones with carbonyl pairs in para positions were found to have a redox potential suited for the anode (i.e., close to 0 V vs SHE), whereas the groups in the ortho positions were found to be appropriate for the cathode (closer to 1.23 V vs. SHE).Irrespective of the number of aromatic rings in the core quinone structures, functionalization with EDGs (−OH, −NH 2 ) resulted in lower potential values and vice versa.Similarly, irrespective of the number of aromatic rings, incorporating the −OH, −NH 2 , −COOH, −SO 3 H, and −PO 3 H 2 groups decreased the solvation energy, thereby likely increasing their solubility.They found that substituting all positions was more effective at decreasing solvation energy than just one, however, at the likely cost of lower synthetic accessibility.
Huynh et al. [37] performed one of the most systematic analyses of the quinone redox mechanism by conducting a broad set of experiments for 25 core quinones under consistent conditions, accompanied by extensive theoretical calculations.They compared the onee − reduction potential in acetonitrile and the two-e − /two-H + (aq) redox potential in water and found no overall linear correlation.However, they did find trends between the two redox potentials within sub-groups based on the type of substituent groups, which has major implications on rational design strategies for ARFBs.To explore this aspect further, they predicted the redox potentials of several benzoquinones with numerous EDGs and EWGs using DFT calculations and found that it is the effects of the functional groups on the pK a values of the reaction intermediates that play a decisive role in tuning the two-e − /two-H + (aq) redox potentials.Their results also confirmed that using functional groups with H-bonding (−CHO, −COCH 3 , −CO 2 CH 3 , −CO 2 H) or negative charge could shift the redox potentials to more positive values, whereas using sterically bulky groups (−C(CH 3 ) 3 ) could shift the redox potentials in the negative direction.Additionally, they found that the strong electronic withdrawing inductive effect of halogens only significantly influences the one-e − redox potential and not so much the two-e − /two-H + (aq) potentials.The reason for this observation was argued to be a net "redox leveling", which implies that stabilization available for e − is canceled by H + (aq) .Nevertheless, the overall magnitude of the effect of using halogens was found to be strong enough to significantly influence the two-e − /two-H + (aq) redox potentials.The HTVS by Er et al. [35] was extended in the recent work of Schwan et al. [42], who attempted to close the gap between singly and fully substituted quinone derivatives by considering a large variational space of possible numbers and positions of functionalization.The authors also looked at the effects of explicit water, intramolecular H-bonding, steric hindrances, and inductive and resonance effects.While the effects of the EDG and EWG substituent groups were in agreement with previous reports, the effect of the number of substituents was not linear.For both benzoquinone and naphthoquinone derivatives, this effect strongly depended on the kind of substituent and its position relative to the carbonyl groups, such that larger cumulative effects on redox potentials were observed for substituents with major resonance effects than for substituents with major inductive effects.For anthraquinone derivatives, the authors found that the substituent −OCH 3 both increased and decreased the redox potential depending on the number and positioning of the functional groups owing to steric effects.They also estimated the influence of (three) explicit water molecules on the calculation of redox potential for all benzoquinone derivatives.Interestingly, this led to a lowering of the redox potential in several cases because of spontaneous Michael addition, which provides strong evidence that the experimental data on benzoquinones may be significantly affected by side reactions.Otherwise, the addition of water molecules generally led to positive shifts in the redox potentials by ∼90 meV com-pared to using just the implicit solvation model.The authors also performed a systematic analysis of the effect of intra-molecular H-bonding by varying the orientation of the hydrogens on the hydroxyl groups-an effect that is ignored routinely in the published literature (except for first-principles MD and QM/MM investigations).As shown in Figure 7, the comparison between various orientations of the hydrogens on hydroxyl groups revealed drastic differences in the relative energies of the conformers for benzoquinones (120 meV), naphthoquinones (310 meV), and anthraquinones (500 meV).
OCH3 both increased and decreased the redox potential depending on the number and positioning of the functional groups owing to steric effects.They also estimated the influence of (three) explicit water molecules on the calculation of redox potential for all benzoquinone derivatives.Interestingly, this led to a lowering of the redox potential in several cases because of spontaneous Michael addition, which provides strong evidence that the experimental data on benzoquinones may be significantly affected by side reactions.Otherwise, the addition of water molecules generally led to positive shifts in the redox potentials by ~90 meV compared to using just the implicit solvation model.The authors also performed a systematic analysis of the effect of intra-molecular H-bonding by varying the orientation of the hydrogens on the hydroxyl groups-an effect that is ignored routinely in the published literature (except for first-principles MD and QM/MM investigations).As shown in Figure 7, the comparison between various orientations of the hydrogens on hydroxyl groups revealed drastic differences in the relative energies of the conformers for benzoquinones (120 meV), naphthoquinones (310 meV), and anthraquinones (500 meV).Comparison between various orientations of the hydrogens on hydroxyl groups that result in drastic differences in the relative energies of the conformers for benzoquinones (left) and naphthoquinones (right).Reprinted with permission from Ref. [42].2022, John Wiley and Sons.
Additionally, Schwan et al. [42] also addressed the issue of stability in terms of geometric structure.Substitution of groups can lead to distortions and strain in the ideally "ring-plane" structures [102], which results in a reduced overlap between the π orbitals, reduced aromaticity, and thus, reduced stability.For all core structures considered, the authors found the distortion effect to be proportional to the number of substituent groups.It is interesting to note that the -OCH3 and -SO3 − groups were found to cause very large distortions and form unphysical structures.However, these groups also allow a very high degree of tunability of both redox potentials and solubility.The existence of such a Comparison between various orientations of the hydrogens on hydroxyl groups that result in drastic differences in the relative energies of the conformers for benzoquinones (left) and naphthoquinones (right).Reprinted with permission from Ref. [42].2022, John Wiley and Sons.
Additionally, Schwan et al. [42] also addressed the issue of stability in terms of geometric structure.Substitution of groups can lead to distortions and strain in the ideally "ring-plane" structures [102], which results in a reduced overlap between the π orbitals, reduced aromaticity, and thus, reduced stability.For all core structures considered, the authors found the distortion effect to be proportional to the number of substituent groups.It is interesting to note that the −OCH 3 and −SO 3 − groups were found to cause very large distortions and form unphysical structures.However, these groups also allow a very high degree of tunability of both redox potentials and solubility.The existence of such a correlation implies that the success in designing novel high-performance molecules by using the cumulative effect of adding more substituent groups is probably going to be limited not only by their synthetic accessibility but also by their inherent structural stability.
Flores et al. [36] considered a combinatorial library of bioinspired benzothiophenoquinones that possess structural similarities with indolequinones for HTVS.Using the data generated from DFT-calculated redox potentials, they also developed a cheminformatics model that predicts redox potentials to within ±0.09 V based solely on molecular connectivity.In contrast to Er et al. [35], they considered all possible numbers (2-4) and configurations of substitutions.Although the general trends for EDGs and EWGs were found to agree with previous reports, there were notable exceptions for trends in the opposite direction for both EDGs and EWGs when comparing single and double substitutions.They also found that, generally, the change in redox potential in either direction is higher when substitutions are made in rings opposite each other as compared to the same ring.This nonlinear dependence was not found to be true for solubility, where increasing the number of solubilizing functional groups decreased the solvation energy.They identified several new molecules both for the anode (E 0 r(SHE) < 0.25 V) and cathode (E 0 r(SHE) > 0.95 V) with solubility > 2 mol/L, that are worthy of further investigation.
Han et al.
[38] considered three different indolequinones as core structures.Their findings on the effect of substituent groups on redox potentials and solubility were quite similar to previous reports.In addition, they showed that indolequinones with carbonyl pairs in ortho positions are not suitable as negolytes and tend to have lower solubilities than the ones with carbonyl pairs at the para positions.From the 180 enumerated indolequinone derivatives, they found ten molecules suitable as negolytes (E 0 r(SHE) < 0.2 V) and nine molecules as posolytes (E 0 r(SHE) > 0.9 V).Among these 19 molecules, only three molecules were identified as having an aqueous solubility > 2 mol/L.
Tabor et al. [39] performed one of the very few HTVS investigations of a variety of quinones with a special focus on their chemical stability in water and its relationship to the redox potential.As electrophiles, quinones can undergo nucleophilic attack by water, thus, resulting in the loss of RAMs.The authors studied two mechanisms of nucleophilic addition of water that limit quinone performance in practical ARFBs, one reversible (Gem-diol formation) and one irreversible (Michael addition), as shown in Figure 8.As descriptors of stability, they considered the equilibrium constant for hydration (log K hyd ) and the change in energy for gem-diol formation and Michael addition/substitution (∆E mic ).The criteria for stability were: (1) log K hyd < 1, which implies that at equilibrium the Gem-diol form of a molecule is less than 10 times as abundant than the oxidized form of the molecule, which is not a very strict criterion because the gem-diol formation is still more favored, and (2) ∆E mic < 0.03 eV, which implies that the molecule is less susceptible to Michael addition than a known quinone alizarin red S for which ∆E mic = 0.03 eV.
data generated from DFT-calculated redox potentials, they also developed a cheminformatics model that predicts redox potentials to within ±0.09 V based solely on molecular connectivity.In contrast to Er et al. [35], they considered all possible numbers (2-4) and configurations of substitutions.Although the general trends for EDGs and EWGs were found to agree with previous reports, there were notable exceptions for trends in the opposite direction for both EDGs and EWGs when comparing single and double substitutions.They also found that, generally, the change in redox potential in either direction is higher when substitutions are made in rings opposite each other as compared to the same ring.This nonlinear dependence was not found to be true for solubility, where increasing the number of solubilizing functional groups decreased the solvation energy.They identified several new molecules both for the anode (E < 0.25 V) and cathode (E > 0.95 V) with solubility > 2 mol/L, that are worthy of further investigation.Han et al. [38] considered three different indolequinones as core structures.Their findings on the effect of substituent groups on redox potentials and solubility were quite similar to previous reports.In addition, they showed that indolequinones with carbonyl pairs in ortho positions are not suitable as negolytes and tend to have lower solubilities than the ones with carbonyl pairs at the para positions.From the 180 enumerated indolequinone derivatives, they found ten molecules suitable as negolytes (E < 0.2 V) and nine molecules as posolytes (E > 0.9 V).Among these 19 molecules, only three molecules were identified as having an aqueous solubility > 2 mol/L.
Tabor et al. [39] performed one of the very few HTVS investigations of a variety of quinones with a special focus on their chemical stability in water and its relationship to the redox potential.As electrophiles, quinones can undergo nucleophilic attack by water, thus, resulting in the loss of RAMs.The authors studied two mechanisms of nucleophilic addition of water that limit quinone performance in practical ARFBs, one reversible (Gemdiol formation) and one irreversible (Michael addition), as shown in Figure 8.As descriptors of stability, they considered the equilibrium constant for hydration (log K ) and the change in energy for gem-diol formation and Michael addition/substitution (∆E ).The criteria for stability were: (1) log K 1, which implies that at equilibrium the Gemdiol form of a molecule is less than 10 times as abundant than the oxidized form of the molecule, which is not a very strict criterion because the gem-diol formation is still more favored, and (2) ∆E 0.03 eV, which implies that the molecule is less susceptible to Michael addition than a known quinone alizarin red S for which ∆E 0.03 eV.Notably, the redox potentials and descriptors for stability were calculated at the calibrated PM7/COSMO level of theory, and the authors demonstrated excellent parity between the PM7/COSMO results and the DFT results for all of the properties of the molecules within the considered chemical space.The structure-property relationships from their work on the effect of positioning of substituent groups and carbonyls on redox potentials were in agreement with previous reports.Notably, the closer the substituent groups were to the carbonyl redox sites, the more pronounced the effect on the redox potential was.The effect was found to occur at different positions for benzo-and naphthoquinones compared to anthraquinones and phenanthrenes, owing to the number and arrangement of rings.
To understand the correlation between redox potential and water stability, the authors first divided their molecules into two groups, one with only a single pair of carbonyls and the other with two pairs of carbonyls.Expectedly, EWGs that make the molecule more electron-deficient will increase the redox potential as well as the propensity towards nucleophilic attack.This general correlation was found to hold well for molecules with a single pair of carbonyls, however, it was substantially less pronounced for molecules with two pairs of carbonyls.Generally speaking, the existence of such a correlation is desirable for negolytes because the lower the redox potential, the lower the propensity toward nucleophilic attack.However, for the same reason, this correlation is undesirable for prospective posolytes.Screening even with very benign stability criteria, the authors were able to provide reasons as to why very few quinones were successful as posolytes (>0.9 V vs. SHE) in full cell cycling tests.Based on their investigation, the authors provided recommendations for designing quinones that can overcome water stability issues.For the mitigation of the Michael addition kinetics, bulky functional groups that provide steric hindrance can be added to the molecules, presumably at the cost of synthetic accessibility.Especially for quinones with multiple carbonyl pairs, the authors suggested tuning the pK a of the reduced intermediates in such a way that the potential is "flattened out" at lower pH values, which would allow otherwise unpromising molecules at pH = 0 to become viable at higher values.This proposition alludes to the same effect at higher pH (>12) which is known from both experimental and computational investigations [57].
For organic ARFBs, there is growing interest in developing single molecules with three distinct, well-separated (>1.0 V) redox states that could potentially serve as both the negolyte and the posolyte.In principle, such a molecule can be used at both electrodes, which minimizes irreversible loss of active material via membrane crossover.These "symmetric" organic ARFBs are potentially an inexpensive, durable, and safe energy storage technology.The strategy of having multiple redox pairs by fusing quinones with aza-aromatic [103,104] motifs in the same molecule was explored by Fornari et al. [40].The authors investigated the viability of having multiple redox moieties on a compact molecular core.They performed a systematic investigation of such three-ring molecules with multiple moieties by comparing the results with a reference system composed of two molecules, each with just one redox moiety.The prototypical moieties used by them included 1,2-quinone, 1,4-quinone, pyrazine, and pyridazine.Their analysis of structure-property relations demonstrated that the electronic interaction between the two different moieties can be exploited to tune the difference between the redox potentials.Their work also confirmed that intramolecular hydrogen bonding between adjacent units has a major effect on tuning the redox potentials.Notably, their results also show that multiple moieties often improve solubility in comparison to molecules with a single moiety.While it is not clear if these molecules will be stable in water or easily react with themselves, another problem that is directly apparent is that molecules with multiple moieties are more difficult to synthesize than their single counterparts.
While many of these studies have looked at the systematic design of virtual libraries, Kristensen et al. [41] investigated a large set of entirely naturally occurring quinones, primarily derived from bacteria and fungi.Remarkably, the distribution of redox potentials spanned a wide range from −1.38 V to 1.49 V vs. SHE.While this range is well beyond the potential limits of hydrogen evolution (E 0 r(SHE) < 0 V) and oxygen evolution (E 0 r(SHE) > 1.23 V) reactions, this work offers a choice from a wide variety of core structures and unconventional functional groups that can be used to design new quinones.It is not entirely inconceivable that such quinones, which would be typically synthetically inaccessible in a chemical lab, could be mass-produced using biotechnological harvesting methods.
Inspired by the excellent performance of quinone moieties typically used as foodcoloring agents and organic dyes [105], Zhang et al. [43] enumerated the space of small quinones with carbonyl pairs on 5-membered rings (3-Cyclopentene-1,2-dione, 4-Cyclopentene-1,3-dione) and heterocyclic substitutions with sulfur, oxygen, and nitrogen atoms.In addition to the DFT calculation of redox potentials, the authors also predicted aqueous solubility using a machine learning model developed in their group [87] and explored the commercial availability of the library molecules on the ZINC database [106] by developing vendor data-mining codes.To quantify the overall effect of substituent groups on the energy density of the core structures as negolytes, they proposed C b = S(1.23 − V) as a metric for the volumetric energy density that is expressed as a product of solubility and cell voltage.It was found that the effects of functional groups on both redox potentials and solubility were in significant contrast to previous findings.For instance, except for -OH, all other functional groups (−F, −SO 3 H, −COOH, −NH 2 ) led to a decrease in solubility for molecules with a single 5-membered ring.Furthermore, heteroatom substitutions led to drastic shifts in the properties of both cyclopentenediones and benzoquinones.Except for a few cases of −OH and −SO 3 H substitution, most derivatives of these motifs were found to have worse predicted energy density C b than their respective core quinones.However, this was not found to be the case for cyclopentenediones fused with benzene rings, which led to the prediction of several substituted molecules with better energy density compared to their respective core quinones.The authors found 16 (out of 3257) molecules to be commercially available, out of which, only seven had enough solubility to allow for experimental characterization.As an important result of their HTVS, the authors identified and experimentally validated Indigo-3(SO 3 H) as a new negolyte, whose electrochemical performance and solubility were found to be competitive to that standard quinone AQDS.This molecule was recommended for further testing and characterization in symmetric ARFBs to reduce irreversible capacity losses due to membrane crossover.
In analogy to Table 3, Table 4 provides a summary of the types of approximations and methods used for calculating equilibrium redox potentials under "Redox Potential", the method/model for estimating aqueous solubility under "Solubility", and specific molecule recommendations based on screening criteria under "Candidates".A glance at this information reveals several interesting aspects about the status of HTVS in quinones.Firstly, the methods for predicting redox potentials are approaching a level of maturity that is quite sufficient for purposes of rapid screening.This is evident from the fact that all surveyed studies use simple linear models to predict the redox potentials within an error of 0.1 V vs. SHE.For an intended ARFB of cell voltage = 1 V, this represents a 10 % error and can be considered small given that the screening is typically conducted in the space of thousands of candidate molecules.It is especially encouraging to see the use of semiempirical methods, which have accuracy on par with DFT but require ∼10 −3 times the compute time.
On the contrary, a proper assessment of methods for predicting the aqueous solubility of quinones cannot be done because of the dearth of consistently generated experimental solubility data for a large and diverse set of quinones.Aqueous solubility prediction is a difficult task irrespective of the adopted approach, and the lack of data specific to quinones hinders the development of surrogate prediction models.As argued earlier, an error of 0.3 LogS units in solubility implies an error of 10 0.3 ≈2 times the actual solubility in mol/L.Given that solubility is directly proportional to the energy and power density, this error is too high even from the standpoint of HTVS.In this regard, better solubility prediction models for quinones need to be developed that are both rapidly executable and have significantly lower prediction errors.Another missing gap in the HTVS of quinones is the chemical space of quinones with a mixture of different substituent groups on the same core molecule.Enumeration of such a virtual library will be a much more demanding task, both combinatorically and computationally, as compared to enumerating a library of quinones with a single type of functional group on each core.However, this is an effort worthy of pursuing further because it can lead to the discovery of quinones that are able to overcome the performance-limiting trade-offs and correlations between various competing performance metrics.Model: Equation ( 7 The next major improvement needed in the HTVS of quinones is the inclusion of descriptors for stability.With a few exceptions, all known HTVS efforts for quinone discovery have focused exclusively on the first major problem of low cell voltage and RAM concentration.While the synthetic accessibility, viscosity, ionic conductivity, toxicity, etc., of quinones are also important to the overall ARFB performance, addressing the second major problem of low stability is of critical importance for any chance of long-term operation and high rechargeability for the ARFB.Typical mechanisms of degradation of quinones include nucleophilic substitution/addition [18][19][20][21][22], disproportionation [107], tautomerization [25], protodesulfonation [108], anthrone formation [18,107], and dimerization [109][110][111].Typically, nucleophilic substitution/addition is the dominant mechanism for quinone posolytes (e.g., benzoquinone derivatives with high redox potentials) due to the electron deficiency of their oxidized form [39].For quinone negolytes, the typical decomposition reaction involves disproportionation of the reduced form, leading to the formation of an electrochemically inactive compound, such as anthrone [18,107].
In an analogy to using experimental redox potential data for accurately predicting redox potentials, new first-principles calculable stability descriptors corresponding to each of these decomposition pathways must be identified from spectroscopic and electrochemical data, which may include information on the type of decomposition products.Understandably, the calculated stability descriptor may not be directly comparable with the available experimental data.For instance, the free energy of disproportionation is not an experimentally accessible quantity under normal circumstances.However, it can be compared directly to the capacity fade rate.The existence of a correlation between the two quantities for a set of quinones will provide a basis for employing the descriptor for HTVS.In addition, the comparison of free energy-based stability descriptors to one another may provide information on the dominance of one decomposition pathway over another.

Summary and Outlook
Driven by advancements in computing power, data analysis, and workflow automation techniques in the digital age, HTVS is a highly effective approach for accelerating the rational design and discovery of high-performance materials.HTVS efforts for quinone discovery have significantly narrowed down the space of potentially high-performing candidates of various sizes, moieties, and substituent groups.Additionally, these studies have provided broadly applicable design principles that have found solid validation in experiments.Nevertheless, wide gaps in the prediction of solubility and stability still need to be bridged to increase the "success rate" of such efforts and the chance of finding candidate quinones that exhibit disruptive performance in ARFBs.From a methodological standpoint, performing HTVS along a hierarchical pipeline, i.e., sequentially down-selecting RAMs by evaluating the redox potential, solubility, and stability, may not necessarily be the most efficient approach.The challenge of material discovery is a multi-objective optimization problem, as high cell voltage, high concentrations of active species, and high stability are all desired simultaneously for any ARFB.However, the performance descriptors are often correlated to each other in such a way that optimizing one makes the other worse, which implies the existence of unavoidable, inherent trade-offs.Unless one can design an overall 'performance rating' of candidates that simultaneously takes into account the contributions of all the performance descriptors, the conventional HTVS approach eliminates potentially high-performing candidates at each sequential step.This could happen, for instance, when a candidate scores well on two out of three performance metrics but scores rather averagely on the third, or when the method used to estimate the performance descriptor has an error large enough that the candidate fails to meet the screening criteria and is unnecessarily eliminated.To avoid such a scenario, an improved approach that goes beyond sequential HTVS is needed.
One way to improve the conventional approach would be to find the pareto-optimal hypersurface [112,113], which essentially represents the set of candidates that are equally "fit" for achieving the desired objective.As a simple example, the volumetric energy density of an ARFB is a product of the cell voltage and the RAM concentration.In this regard, the pareto-optimal boundary curve could represent all candidates for whom the product of their redox potential and solubility satisfies some minimum value, irrespective of the individual redox potential or solubility values.Such an analysis will not only make the screening more robust (i.e., less sensitive to prediction errors), but it will also help set theoretical limits on the ARFB's performance.Further, the knowledge of this boundary can be used to extract structure-property relationships, which can in turn inspire the search for new candidates that are outliers to those structure-property relationships by using powerful methods, such as reinforcement learning [114,115].

Figure 1 .
Figure 1.Schematic showing a typical redox flow battery with the tanks containing the dissolved catholyte (also known as posolyte) and the anolyte (also known as negolyte) materials.The electrolytes at each electrode are circulated using pumps.Reprinted with permission from Ref. [4].2022, Elsevier.

Figure 2 .
Figure 2. Schematic diagram of the typical implementation of an HTVS of quinones with sequential funnels associated with performance descriptors.

Figure 2 .
Figure 2. Schematic diagram of the typical implementation of an HTVS of quinones with sequential funnels associated with performance descriptors.

Figure 3 .
Figure 3. (a) Redox and acid−base equilibria between benzoquinone (Q) and hydroquinone in acid (QH ), along with their monodeprotonated (QH ) or dideprotonated (Q ) forms.(b) Typical Pourbaix diagram of the redox couple, where the regions are marked by the thermodynamically most stable acido−basic forms in the entire pH range.Figure adapted from Gaudin et al. [57] (under Creative Commons Attribution 4.0).

Figure 3 .
Figure 3. (a) Redox and acid−base equilibria between benzoquinone (Q) and hydroquinone in acid (QH 2 ), along with their monodeprotonated (QH − ) or dideprotonated (Q 2− ) forms.(b) Typical Pourbaix diagram of the redox couple, where the regions are marked by the thermodynamically most stable acido−basic forms in the entire pH range.Figure adapted from Gaudin et al. [57] (under Creative Commons Attribution 4.0).

Figure 4 .
Figure 4. (a) Schematic diagram showing the workflow for geometry optimization (OPT) and single point energy (SPE) calculations at various levels of theory used by Zhang et al. [55] (Figure adapted from Ref. [55] under Creative Commons Attribution 4.0).The hollow arrows denote the differences between gas-phase reaction energies obtained at different levels of theory.The solid black arrows denote the differences between energy computed using DFT on a fixed geometry obtained from a lower-level preceding theory and the geometry computed at that lower-level preceding theory.The dotted gray arrows represent the solvation effect from gas-phase SPE to solution-phase SPE when the implicit solvation model is considered.(b) Methodology of prediction of Pourbaix diagrams from COSMO-RS employed in the work of Gaudin et al. [57].(Figure adapted from Ref. [57] under Creative Commons Attribution 4.0)

Figure 4 .
Figure 4. (a) Schematic diagram showing the workflow for geometry optimization (OPT) and single point energy (SPE) calculations at various levels of theory used by Zhang et al. [55] (Figure adapted from Ref. [55] under Creative Commons Attribution 4.0).The hollow arrows denote the differences between gas-phase reaction energies obtained at different levels of theory.The solid black arrows denote the differences between energy computed using DFT on a fixed geometry obtained from a lower-level preceding theory and the geometry computed at that lower-level preceding theory.The dotted gray arrows represent the solvation effect from gas-phase SPE to solution-phase SPE when the implicit solvation model is considered.(b) Methodology of prediction of Pourbaix diagrams from COSMO-RS employed in the work of Gaudin et al. [57].(Figure adapted from Ref. [57] under Creative Commons Attribution 4.0).

Figure 5 .
Figure 5. Thermodynamic cycle linking Gibbs free energy of dissolution of an organic molecule to its sublimation and solvation-free energies.The linking requires a change in the standard state of the gas-phase molecule (shown in the center).

Figure 5 .
Figure 5. Thermodynamic cycle linking Gibbs free energy of dissolution of an organic molecule to its sublimation and solvation-free energies.The linking requires a change in the standard state of the gas-phase molecule (shown in the center).

28 Figure 6 .
Figure 6.Schematic showing examples of the various structural and configurational degrees of freedom that are modified when creating a virtual library: (a) core variations and (b) substituent variations.Huynh et al. [37] performed one of the most systematic analyses of the quinone redox mechanism by conducting a broad set of experiments for 25 core quinones under consistent conditions, accompanied by extensive theoretical calculations.They compared the one-e reduction potential in acetonitrile and the two-e /two-H redox potential in

Figure 6 .
Figure 6.Schematic showing examples of the various structural and configurational degrees of freedom that are modified when creating a virtual library: (a) core variations and (b) substituent variations.

Figure 7 .
Figure 7.Comparison between various orientations of the hydrogens on hydroxyl groups that result in drastic differences in the relative energies of the conformers for benzoquinones (left) and naphthoquinones (right).Reprinted with permission from Ref.[42].2022, John Wiley and Sons.

Figure 7 .
Figure 7.Comparison between various orientations of the hydrogens on hydroxyl groups that result in drastic differences in the relative energies of the conformers for benzoquinones (left) and naphthoquinones (right).Reprinted with permission from Ref.[42].2022, John Wiley and Sons.

Figure 8 .
Figure 8.The mechanisms of chemical instability considered by Tabor et al. [39].(a) Gem-diol formation, which is reversible but hampers kinetics.(b) Michael addition on the molecule in which a single vulnerable site of attack and subsequent tautomerization is shown.

Figure 8 .
Figure 8.The mechanisms of chemical instability considered by Tabor et al. [39].(a) Gem-diol formation, which is reversible but hampers kinetics.(b) Michael addition on the molecule in which a single vulnerable site of attack and subsequent tautomerization is shown.

Funding:
A.K. thanks the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for funding this research effort under the Germany's Excellence Strategy-Exzellenzcluster 2186 "The Fuel Science Center" ID: 390919832.Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.

Table 2 .
List of abbreviations and mathematical symbols and their meanings.

Table 3 .
Summary of the virtual libraries enumerated in HTVS literature for quinones.

Table 4 .
Summary of the methods used for redox potential and solubility prediction in available HTVS literature, along with the number of screened candidates and their respective criteria.