All DFT calculations were performed using Gaussian 16 revision C.01 65 software. First, the semi-empirical PM7 66 method was used to optimize the geometry of molecular structures solvated in PCM implicit water model. Note that for each molecule the initial atomic coordinates were prepared by Maestro, as mentioned before, and visually inspected by the authors. At this stage, the energetics of two different structures were compared in order to determine the most stable tautomers of the individual QH -, QH, and QH + forms. After the structure optimization, harmonic vibrational frequencies were assessed for free energy stability evaluation. The structure of a molecule with an imaginary frequency less than -200 cm -1 was modulated and reoptimized until the stable structure was obtained. Then, the total energy was computed for each optimized structure using M06-2X 69 together with a Def2-TZVP basis set and the more accurate SMD solvation model. This exchangecorrelation functional was found to provide the best accuracy for predicting redox potentials of organic molecules. To compute the Gibbs free energy, thermal correction to Gibbs free energy and total energy were obtained from the computationally cheaper PM7 (the first step) and high accuracy M06-2X calculations (the second step), respectively. The thermal correction includes the zero-point vibrational energy as well as vibrational enthalpy and entropic contribution to the free energy. This procedure leads to moderate computational cost and has been validated previously. For high-throughput screening, we used Random Forest Regressor (RFR) as implemented in Scikit-learn package of Python. The optimal set of hyperparameters was computed using a cross-validation score over a grid of predefined space. More specifically, the training data was split into ten groups, or folds, where nine were used to train the model and one is used to evaluate its performance. Mean squared error (MSE) was used as an evaluation metric for hyperparameter optimization. To train our RFR model, the data set was initially randomized, then 80% of the data were used for training the model, while 20% for validation.