# Prediction of Pourbaix diagrams of quinones for redox flow battery by COSMO-RS

Théophile Gaudin, Jean-Marie Aubry

► To cite this version:

Théophile Gaudin, Jean-Marie Aubry. Prediction of Pourbaix diagrams of quinones for redox flow battery by COSMO-RS. Journal of Energy Storage, 2022, 49, pp.104152. 10.1016/j.est.2022.104152. hal-04146155

HAL Id: hal-04146155

https://hal.univ-lille.fr/hal-04146155

Submitted on 29 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research Papers

# Prediction of Pourbaix diagrams of quinones for redox flow battery by COSMO-RS

Théophile Gaudin*, Jean-Marie Aubry

CNRS, UMR 8181, UCCS-Unité de Catalyse et Chimie du Solide, Centrale Lille, University of Lille, Artois University, Lille F-59000, France

<figure><img src="image_3.png" /></figure>

## ABSTRACT

Redox-flow batteries are relevant to store energy from intermittent sources such as solar panels or wind turbines, thereby smoothing their energy supply. Up to now, most redox-flow batteries are based on vanadium. Vanadium is a rare and expensive material, thus recent research has focused on redox-flow batteries based on organic compounds and, in particular, anthraquinones as electroactive materials. However, the tunability of organic chemistry poses a needle-in-haystack challenge as the structures exhibiting the most desirable electrochemical properties may be hard to pinpoint. Moreover, the low water solubility of the most readily available anthraquinones may hinder their use as battery electrolytes. To aid in such endeavor, a theoretical approach is proposed to predict Pourbaix diagrams of redox-active organic compounds, allowing in silico anticipation of their electrochemical behavior in the entire pH range. DFT/COSMO-RS predicted pKₐ and reduction potentials are in good agreement with experimental data, and the resulting calculated Pourbaix diagrams are also in agreement with 4 experimental ones from literature, proving the reliability of the method. Finally, the effect of nature and position of some functional groups on the anthraquinone backbone is discussed, illustrating the power of the method to both understand and quantify the electrochemical activity of redox active organic materials.

## 1. Introduction

Redox flow batteries (RFBs) are secondary (i.e. reusable) batteries in which the energy conversions are based on the reversible electrochemical reactions of two redox couples, which are dissolved in liquids and separated by a membrane that allows ion exchange [1]. RFBs are ideally suited for relatively large stationary applications (capacities of 1 kWh to 10 MWh) [2]. In particular, RFBs have environmental relevance as a suitable technology for energy storage from intermittent sources such as solar or wind power [3].

Water is most commonly chosen as a solvent for electrochemical reactions in RFBs, since it is widely available and provides better ion solubility and conductivity than most organic solvents [4]. In addition, as the most common and safest solvent on earth, water is the solvent of choice wherever sustainability is a key challenge. Nevertheless, it is important to note that in water, pH effects can impact which protonated form of the reduced and oxidized species participate in the redox reaction. This in turns influences their electrochemical activity as well as their water solubility, and thus taking into account pH effects is relevant to the energy storage capacity of the RFB. Because of this, pH-potential diagrams (named Pourbaix diagrams) of the electrolyte species are important to consider at early stages of the RFB's design. Thermodynamically, Pourbaix diagrams can be reconstructed based on pKₐ and standard reduction potentials of the involved species.

Currently, most RFBs are based on vanadium, taking advantage of the various oxidation states of vanadium as both the negolyte and posolyte [5], leading to easier design and operation of batteries [6]. However, vanadium is relatively rare and expensive [7]. Thus, in recent years, RFBs based entirely or partly on more sustainable organic electrolytes have attracted increasing interest [8]. For example, quinones, and in particular anthraquinones, are cheap and renewable. However, anthraquinone itself has very low solubility in water. Fortunately, quinones are a tunable class of compounds as different functional groups can be attached at their multiple substitution positions (cf. Scheme 1). Consequently, their water solubility can be increased by polar functional groups.

Nevertheless, such functionalization can change the potential of their electrochemical reaction at any pH, and therefore alter the energy storage capacity of the RFB.

Since the functional groups grafted onto the quinones, and the phenolic alcohols of their reduced forms, can be pH-sensitive in water, the reduction potential changes as a function of the pH of the solution. Pourbaix diagrams represent the evolution of reduction potential with the pH, and thus are targets of interest to anticipate the electrochemical behavior of quinones in aqueous RFB. In Fig. 1, a typical Pourbaix diagram for benzoquinone in water is sketched.

As detailed below, predictive models exist for the properties that characterize Pourbaix diagrams (reduction potentials and pKₐ). pKₐ is

<footnote>* Corresponding author.
E-mail address: gaudin.theophile@gmail.com (T. Gaudin).</footnote>
<figure><img src="image_4.png" /><figcaption>**Scheme 1.** General redox reaction of benzo-, naphtho-, and anthra-quinones (Q), giving the corresponding hydroquinones (QH₂).</figcaption></figure>

the pH at which there are equal concentrations of a species in its protonated and deprotonated form. Some species, including in particular the reduced forms of the quinones, have several labile hydrogens, and therefore several pKₐ in the experimentally accessible pH range. There is a wealth of existing methods for predicting pKa [9], including empirical methods like Quantitative Structure Property Relationships (QSPR) [10], group contributions [11], or more theoretical approaches like quantum chemistry [12] combined with explicit water molecules [13] and/or solvation models like the CONductor-Like Screening MOdel (COSMO-RS) [14]. Although errors of ± 0.2-0.4 pKₐ units are already achievable from models fitted for particular families of molecules [9,11] (with the obvious shortcoming of a smaller applicability domain), errors of ± 0.6-1.5 pKₐ units are usual for more general methods [9,10].

By definition, the reduction potential is a difference between two electric potential differences. The first potential difference is the one associated with the acquisition of one electron by a molecule, and the other one is associated with a reference reduction reaction, such as the reduction of the proton. The reference potential difference being well-known (e. g., E<sub>H+/H2</sub> = 4.42V for the reduction of the proton at ambient conditions in water [15]), when predicting the reduction potential for a given redox reaction, only the reduction of the molecule of interest remains to be modelled.

Here, two main types of reduction potentials may be the target for prediction, the electron-affinity EA, or single-electron reduction potential, and the standard electrode potential E. Both EA and E directly depend on standard Gibbs free energies of their respective reaction ΔₜG⁰, and the number of involved electrons (1 for reaction 1 and 2 for reaction 2) and Faraday's constant F Eqs. (1) and ((2)):

$$Q + \frac{1}{2}H_2 \xrightarrow{\Delta_r G_1^0} Q^{\circ-} + H^+ \Rightarrow EA = -\frac{\Delta_r G_1^0}{F} \quad (1)$$

$$Q + H_2 \xrightarrow{\Delta_2 G_2^0} QH_2 \Rightarrow E = -\frac{\Delta_2 G_2^0}{2F} \quad (2)$$

Prediction of the EA has been explored by many studies [16]. In quantum chemistry, this implies the modeling of systems with unpaired electrons (i.e. open-shell calculations), which are more complex to handle than species with paired electrons [17]. Errors of about ± 65 mV are expected for particular subclasses of compounds using tailored level of theories and solvation models, but more general and unbiased approaches lead to errors of ± 150-200 mV. The second type, E, is the reduction potential [18]. It generally implies both a reduced and an oxidized form that are thermodynamically stable, and does not involve radicals and open-shell calculations. In the context of Pourbaix diagrams, we are interested in the reduction potential when proton activity is 1 mol/L, i.e. at pH = 0, and measured with respect to a reference (generally the standard hydrogen electrode), since it is the one available from experiments [8]. This potential is also called formal electrode potential and denoted E⁰. Prediction of formal electrode potential of anthraquinones from first principles has been explored by Kim et al. [15] using DFT and an implicit solvent model with good results (errors of about 30-60 mV), but for some anthraquinones, it appeared that explicit hydrogen bonding of water could not be neglected, and Kim et al. had to use sophisticated QM/MM calculations to account for such effect, meaning that their method may fail in some cases. Additionally, researchers have been using a correlation between the standard electrode potential and the calculated heat of formation at 0 K [19]. This approach

<figure><img src="image_5.png" /><figcaption>**Scheme 2.** Methodology of prediction of Pourbaix diagrams from COSMO-RS used in this work.</figcaption></figure>

has been tested on 6 anthraquinones in a recent work, and a good correlation coefficient R² of 0.97 was obtained [20]. However, this method requires calibration of the linear fitting coefficients to a set of similar molecules before any prediction can be made. Finally, Fornari and De Silva [21] attempted prediction of pH-dependent redox potentials. However, their method is not fully from first principles as the pKa used in their method comes from machine learning. Moreover, the errors are substantial (mean absolute error of about 150 mV on pH-dependent potentials). This could be due to the use of structures optimized both in the gas phase and in solution, increasing the potential sources of error.

In this work, we present a method to predict Pourbaix diagrams based on ab initio calculations and COSMO-RS theory. Briefly, our proposed method consists of prediction of pKₐ from the COSMOtherm software [22,23], and prediction of the redox potential, either at pH = 0 by accounting for the protonation states at this standard point (labelled E⁰ in Fig. 1) or at high pH in which no proton exchange is expected (labelled E∞ in Fig. 1). Then, the Pourbaix diagram is reconstructed accounting for the relative population of each protonated state for both reduced and oxidized forms of the quinone, based on combining Nernst and Henderson-Hasselbalch equations.

## 2. Experimental

The general methodology developed in this work to predict Pourbaix diagrams (or pH-dependent reduction potentials) is presented in Scheme 2. We first establish the various structures of all oxidized and

<figure><img src="image_6.png" /><figcaption>**Scheme 3.** Redox and acid-base equilibria between benzoquinone (Q) and hydroquinone in acid (QH₂), monodeprotonated (QH⁻) or dideprotonated (Q²⁻) forms.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Scheme 4.** Thermodynamic cycle illustrating Gibbs free energies associated with reduction and protonation of a quinone.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 1.** Typical Pourbaix diagram of the redox couple benzoquinone (Q)/hydroquinone (QH₂) depicted in all their thermodynamically stable acido-basic forms in the 0–14 pH range.</figcaption></figure>

reduced forms of the molecules in all of their protonation states, for example, by drawing them on a molecule editor and saving them as SMILES strings. Then, we optimize their geometries in a perfect conductor, yielding COSMO files for each of them. Using the COSMO files, we apply COSMO-RS to get their Gibbs free energies in water which allows calculating the dissociation constant pKₐ and the formal reduction potential E⁰. Finally, using thermodynamics, we determine the predominant species as a function of pH and E in order to build their full Pourbaix diagrams.

## 2.1. Thermodynamic workflow

As we study the conversion of species between different redox forms and protonation states (cf. Scheme 3) in water, we need the levels of energy of each of the involved species in water.

To calculate this energy level, a simple thermodynamic cycle is employed (cf. Scheme 4): the Gibbs free energies of the molecules in the ideal conductor GₓCOSMO, with X = Q, QH₂, QH⁻ and Q²⁻) are first calculated (at the DFT level BP86 [24,25] with TZVP [26] basis set using the COSMOconf [27] and Turbomole [28] softwares) in the ideal conductor (COSMO [29, 30]). Then, COSMO-RS theory [30–32](which was applied using the COSMOtherm [22] software) allows to calculate

the Gibbs free energy of transfer from the ideal conductor to water. In the context of COSMO-RS theory, this Gibbs free energy corresponds to the chemical potential of the molecule in water µₓaq. By adding the two contributions together, we obtain the Gibbs free energy of the molecule in water, GₓCOSMO = GₓCOSMO + µₓaq. Finally, the relevant states are compared to yield the Gibbs free energy associated with the relevant transformation, as shown in Eqs. (3)–(6) and Scheme 4.

$$ \Delta G_{aq,0}^{\text{red,Q}} = G_{aq}^{\text{OH}_2} - G_{aq}^{\text{Q}} - 2G_{aq}^{\text{H}^+} \Rightarrow E^0 \quad (3) $$

$$ \Delta G_{aq,\infty}^{\text{red,Q}} = G_{aq}^{Q-} - G_{aq}^{Q} \Rightarrow E^{\infty} \quad (4) $$

$$ \Delta G_{aq}^{\text{diss,QH}_2} = G_{aq}^{\text{QH}^+} + G_{aq}^{\text{H}^+} - G_{aq}^{\text{QH}_2} \Rightarrow pKa_1 \quad (5) $$

$$ \Delta G_{aq}^{\text{diss,QH}^+} = G_{aq}^{Q-} + G_{aq}^{\text{H}^+} - G_{aq}^{\text{QH}^+} \Rightarrow pKa_2 \quad (6) $$

One could also formulate the thermodynamic cycle starting from the gas phase and then calculating the Gibbs free energy of solvation into the liquid phase [12,16]. However, in this approach, two calculations are required, one for the molecule in the gas phase (to get the gas phase energy) and one for the molecule in the ideal conductor (to get the COSMO file). By directly calculating the molecule in the proper reference state for COSMO-RS calculations, there is no need for any gas phase calculation, which also means that one potential source of error is avoided and the calculations are potentially more robust.

## 2.2. pKₐ of acids in aqueous solutions

The pKₐ is directly related to the Gibbs free energy of dissociation. For example, as illustrated in Scheme 4, the pKₐ of QH⁻, describing the deprotonation of QH⁻ into Q²⁻, writes as Eq. (7) [33].

$$ pK_a \cdot QH^+ = \frac{\Delta G_{aq}^{\text{diss,QH}^-}}{RT \ln 10} - \log[H_2O] \equiv a\Delta G_{aq}^{\text{diss,QH}^+} + b \quad (7) $$

$\Delta G_{water}^{\text{diss,QH}^+}$ is obtained by the procedure outlined in Scheme 4 and Eq. (6), and $G_{aq}^{\text{H}^+}$ from Eq. (6) can be either calculated or extracted from tabulated experimental data [12].

In a series of papers, Eckert et al. [23,34-36] explored the predictive power of combining DFT with COSMO-RS in order to predict the pKₐ of acids and bases in various media. They found that in water, the coefficient a and b in Eq. (7), though in the correct order of magnitude, were not exactly 1/RTln10 and log[H₂O], respectively, for acids in water [34]. Further calculations [23] have shown that this inconsistency can be fully accounted for by the occurrence of specific hydration of the ionized groups, notably for high values of pKₐ, generally associated to high concentrations of charge on the anion, which was not taken into account as $\Delta G_{water}^{\text{diss,}AH}$ was obtained without placing explicit waters around the modelled structures. The COSMOtherm software, for acids in water, automatically corrects the slope and intercept to account for this inconsistency, avoiding the rather complex [13] step of placing explicit waters at suitable positions around the relevant species, and calculates $G_{aq}^{\text{H}^+}$ from the COSMO files of water and hydronium at the appropriate temperature [23].

## 2.3. Reduction potential of quinones E⁰

The reduction potential labelled E⁰ in Fig. 1, for a given quinone Q without acid or basic side groups, is directly derived the Gibbs free energy of the reduction reaction. For the reduction of Q into Q²⁻, this Gibbs free energy of reduction can be expressed as (cf. Scheme 4 and Eq. (3)).

The standard hydrogen electrode potential, E<sub>H+/H₂</sub>, is chosen as a reference for the calculated potential, as it is typically used as a reference for experimental data available in literature [15,20]. In this work, we use the commonly chosen value [15] E<sub>H+/H₂</sub> = 4.42V. The
theoretical relationship between $\Delta G_{water}^{red,Q}$ and $E_{water}^{\infty,Q}$ for a 2-electron re-action is [37]:

$$E_{water}^{\infty,Q} = -\frac{\Delta G_{aq,\infty}^{red,Q}}{2F} - E_{H^+/H_2} \quad (8)$$

## 2.4. Prediction of the standard electrode potential $E^0$

The Nernst equation is:

$$E = E^0 - 2.303 \frac{\langle n_{H^+} \rangle_{pH} RT}{n_{e-F}} pH \quad (9)$$

where $\langle n_{H^+} \rangle_{pH}$ the number of liberated protons upon the redox reaction at a given pH and $n_e$ is the number of electrons exchanged during the redox reaction. Since experimental data are generally available for standard electrode potential, i. e., at pH = 0, $E^0$ is the property to be targeted for assessment of the predictive power of the method. Since the pH-dependent term of Eq. (9) vanishes due to the standard hydrogen electrode potential being also taken at pH=0, $E^0$ can be predicted taking advantage of the experimentally known $G_{aq,0}^H$, 272.0 kcal/mol [12]. Then, like for the deprotonated species (Eq. (8)), for a two-electron process, the standard electrode potential is computed as:

$$E_{water}^{0,Q} = -\frac{\Delta G_{aq,0}^{red,Q}}{2F} - E_{H^+/H_2} \quad (10)$$

In the more general case, standard electrode potentials should be calculated based on the predominant acidobasic form at pH = 0, which can be identified from the pKₐs of the species at hand.

## 2.5. Construction of Pourbaix diagrams from pKₐ and reduction potential

In the general case, both oxidized and reduced species on either side of the redox reaction are pH-sensitive. Furthermore, species may co-exist in different protonation states at pH close to the pKₐ of their labile functional groups. However, for clarity, we illustrate the construction of simulated Pourbaix diagrams on the simple example of benzoquinone (symbolized by Q). For general equations, the interested reader is referred to Appendix A.

First, we determine the population of each protonation form QH₂, QH⁻, and Q²⁻ of the reduced benzoquinone (hydroquinone). To do so, we combine the COSMO-RS computed dissociation constants with the material balance for all forms of hydroquinone.

With reference to Scheme 4, we have the Henderson-Hasselbalch equations for the two deprotonation reactions:

$$10^{pH-pK_a^1} = \frac{[QH^-]}{[QH_2]} \quad (11)$$

$$10^{pH-pK_a^2} = \frac{[Q^2^-]}{[QH^-]} \quad (12)$$

The following equations define the amounts (mole fractions) of each form of hydroquinone:

$$x_{Q^{2-}} = \frac{[Q^{2-}]}{C_{tot,QH_2}}; \quad x_{QH^{-}} = \frac{[QH^{-}]}{C_{tot,QH_2}}; \quad x_{QH_2} = \frac{[QH_2]}{C_{tot,QH_2}} \quad (13)$$

where $C_{tot,QH_2}$ is the total concentration of hydroquinone introduced in the system. We then introduce the material balance for hydroquinone protonation forms in solution:

$$x_{Q^{2-}} + x_{QH^{-}} + x_{QH_2} = 1 \quad (14)$$

Combining eqs. 11-14 yields expressions for the fraction of hydroquinone in each of the QH₂, QH⁻ and Q²⁻ forms at a given pH:

$$x_{QH_2} = \frac{1}{1 + 10^{pH - pK_a^1} + 10^{2pH - pK_a^1 - pK_a^2}} \quad (15)$$

$$x_{QH^{-}} = \frac{10^{pH - pK_a^1}}{1 + 10^{pH - pK_a^1} + 10^{2pH - pK_a^1 - pK_a^2}} \quad (16)$$

$$x_{Q^{2-}} = \frac{10^{2pH - pK_a^1 - pK_a^2}}{1 + 10^{pH - pK_a^1} + 10^{2pH - pK_a^1 - pK_a^2}} \quad (17)$$

With the population of each protonation form of hydroquinone at hand, the average number of exchanged protons for the reduction reaction ($\langle n_{H^+} \rangle_{pH}$ in Eq. (9)), can be calculated by summing up the contributions of each relevant form (eq. 18):

$$\langle n_{H^+} \rangle_{pH} = 2x_{QH_2} + x_{QH^{-}} \quad (18)$$

Of course, in the general case (cf. Appendix A), even the oxidized form may exhibit labile protons, and these have to be subtracted from $\langle n_{H^+} \rangle_{pH}$ according to the contribution of each oxidized protonation form. Once $\langle n_{H^+} \rangle_{pH}$ is calculated for the entire pH range, any Pourbaix diagram can be fully reconstructed from a single reduction potential and all the pKₐs of oxidized and reduced species, and the distribution of protonated forms can even be accounted for, in the vicinity of pKₐs. To the best of our knowledge, the general formulation presented in Appendix A has never been published in the literature and could find independent application in retro-fitting pKₐ values from Pourbaix diagrams. As in Eq. (9), $E^0$, by definition, does not depend on pH and, on small pH intervals, $\langle n_{H^+} \rangle_{pH}$ can also be considered pH-independent, the derivative of Eq. (8) can be evaluated numerically thanks to Eq. (19).

$$\left. \frac{dE}{dpH} \right|_{pH} = -2.303 \frac{\langle n_{H^+} \rangle_{pH} RT}{n_{e-F}} \quad (19)$$

The Nernst equation can be reformulated, in view of Eq. (19), in an integral form Eqs. (20) and (21).

$$E = E^0 + \int_{0}^{pH} \frac{dE}{dpH} dpH \quad (20)$$

$$E = E^{\infty} - \int_{pH}^{+\infty} \frac{dE}{dpH} dpH \quad (21)$$

where, for some quinone in water, $E^\infty$ can be estimated based on quantum chemistry and COSMO-RS calculations thanks to Eq. (8). Eqs. (20) and (21) can be discretized for small pH intervals (e. g. ΔpH = 0.01), cf. Eqs. (22) and (23).

$$E \approx E^0 + \sum_{pH=0}^{pH=p_{H\text{target}}} \left. \frac{dE}{dpH} \right|_{pH} \Delta pH \quad (22)$$

$$E \approx E^{\infty} - \sum_{pH=p_{H\text{target}}}^{pH=p_{H\text{max}}} \left. \frac{dE}{dpH} \right|_{pH} \Delta pH \quad (23)$$

If a reliable prediction of $E^0$ can be obtained by use of Eq. (2), then Eq. (22) is the preferred route to predict Pourbaix diagrams. Using the predicted pKₐ, values of $\langle n_{H^+} \rangle_{pH}$ can be predicted over all the relevant pH range, and numerical integration from Eq. (22) combining predicted $E^0$ with $\langle n_{H^+} \rangle_{pH}$ at each pH yields the full Pourbaix diagram. Alternatively, Eq. (23) can be used if the model underlying use of Eq. (10) has deficiencies. By this approach, predicted $E^\infty$ (Eq. (8)) can be combined with consideration that $\langle n_{H^+} \rangle_{pH} \approx 0$ at sufficiently high pH to reconstruct the Pourbaix diagram from the other end. The approach presented in this section has been implemented in an in-house Python script to reconstruct the Pourbaix diagrams from predicted reduction potentials and pKₐ.
<caption>**Table 1**</caption>

Molecules studied in this work. Bold: redox potential predictions, italic: pKa predictions, underlined: Pourbaix diagram predictions.

<table><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>2-phenylphenol</td><td>2</td><td>3-phenylphenol</td><td>3</td><td>4-indanol</td></tr><tr><td>5</td><td>1,2-benzoquinone</td><td>6</td><td>1,4-benzoquinone</td><td>7</td><td>2-hydroxyquinone</td></tr><tr><td>9</td><td>tetrahydroxyquinone</td><td>10</td><td>hydroquinone</td><td>11</td><td>tetramethylhydroquinone</td></tr><tr><td>13</td><td>2,5-dinitro-hydroquinone</td><td>14</td><td>2,5-dichloro-hydroquinone</td><td>15</td><td>1-naphtol</td></tr><tr><td>17</td><td>8-hydroxyquinoline</td><td>18</td><td>4-methyl-8-hydroxyquinoline</td><td>19</td><td>3,4-dimethyl-8-hydroxyquinoline</td></tr><tr><td>21</td><td>1,2-naphthoquinone</td><td>22</td><td>1,4-naphthoquinone</td><td>23</td><td>lawsone</td></tr><tr><td>25</td><td>1,4-dihydroxynaphтенe</td><td>26</td><td>9,10-phenanthraquinone</td><td>27</td><td>-anthraquinone</td></tr><tr><td>29</td><td>1-nitroanthraquinone</td><td>30</td><td>anthraquinone-1-carboxylic acid</td><td>28</td><td>1-hydroxyanthraquinone</td></tr><tr><td>33</td><td>2-hydroxyanthraquinone</td><td>34</td><td>2-nitroanthraquinone</td><td>31</td><td>1-aminoanthraquinone</td></tr><tr><td>37</td><td>2-aminoanthraquinone</td><td>38</td><td>anthraquinone-2-sulfonate</td><td>35</td><td>anthraquinone-2-carboxylic acid</td></tr><tr><td>41</td><td>chrysazine</td><td>42</td><td>anthraflavic acid</td><td>36</td><td>anthraquinone-2-methylcarboxylate</td></tr><tr><td>45</td><td>anthraquinone-1,8-disulfonate</td><td>46</td><td>anthraquinone-2,6-disulfonate</td><td>39</td><td>alizarine</td></tr><tr><td>49</td><td>1,4-methoxy-anthraquinone</td><td>50</td><td>1-methoxy-4-propen-oxyaanthraquinone</td><td>40</td><td>quinizarine</td></tr><tr><td>53</td><td>alizarine-3-sulfonate (red S)</td><td>51</td><td>1,4-dipropenoxy-anthraquinone</td><td>44</td><td>1-amino-2-ethyl-anthraquinone</td></tr><tr><td>57</td><td>1,5-dihydroxyanthraquinone-2,7-disulfonate</td><td>54</td><td>1-amino-2,3-dimethyl-anthraquinone</td><td>48</td><td>anthraquinone-1,8-disulfonate</td></tr><tr><td></td><td></td><td>55</td><td>1-amino-2,4-dimethyl-anthraquinone</td><td>52</td><td>3,3'-[(9,10-anthraquinone-diyl)bis(3-methylbutanoic acid)]</td></tr><tr><td></td><td></td><td>59</td><td>anthraquinol-2,7-disulfonate</td><td>56</td><td>1,8-dihydroxy-anthraquinone-2,7-disulfonate</td></tr></tbody></table>
<table><caption>**Table 2**
Experimental vs. predicted pK<sub>a</sub> for phenols, naphthols, quinones, quinols and anthraquinones and anthraquinols</caption><thead><tr><th rowspan="2">Structure</th><th colspan="2">Molecule</th><th colspan="3">Experimental pKa</th><th colspan="3">Predicted pKa</th><th rowspan="2">Ref.</th></tr><tr><th>N°</th><th>Substituents<sup>(i)</sup></th><th>1<sup>st</sup>≡NH<sup>+</sup></th><th>1<sup>st</sup>.OH</th><th>2<sup>nd</sup>.OH</th><th>1<sup>st</sup>≡NH<sup>+</sup></th><th>1<sup>st</sup>.OH</th><th>2<sup>nd</sup>.OH</th></tr></thead><tbody><tr><td rowspan="3">[Chemical Structure of 7]</td><td>7</td><td>R<sub>2</sub> = OH</td><td></td><td>4.2</td><td></td><td></td><td>3.7</td><td></td><td>[39]</td></tr><tr><td>8</td><td>R<sub>2,5</sub> = OH</td><td></td><td>2.7</td><td>5.2</td><td></td><td>3.2</td><td>4.6</td><td>[40]</td></tr><tr><td>9</td><td>R<sub>2,3,5,6</sub> = OH</td><td></td><td>4.8</td><td>6.8</td><td></td><td>4.0</td><td>4.4</td><td>[39]</td></tr><tr><td rowspan="12">[Chemical Structure of 1-4]</td><td>1</td><td>R<sub>1</sub> = Ph</td><td></td><td>10.0</td><td></td><td></td><td>10.1</td><td></td><td>[39]</td></tr><tr><td>2</td><td>R<sub>2</sub> = Ph</td><td></td><td>9.7</td><td></td><td></td><td>9.6</td><td></td><td>[39]</td></tr><tr><td>3</td><td>R<sub>1,2</sub> =</td><td></td><td>10.3</td><td></td><td></td><td>9.9</td><td></td><td>[39]</td></tr><tr><td>4</td><td>R<sub>1</sub> = NO<sub>2</sub>, R<sub>4</sub> = Ph</td><td></td><td>6.7</td><td></td><td></td><td>7.4</td><td></td><td>[39]</td></tr><tr><td>10</td><td>R<sub>3</sub> = OH (hydroquinone)</td><td></td><td>10.2<sup>(ii)</sup></td><td></td><td>11.7<sup>(ii)</sup></td><td>10.2</td><td>10.2</td><td>[39]</td></tr><tr><td>11</td><td>R<sub>1,2,4,5</sub> = Me, R<sub>3</sub> = OH</td><td></td><td>11.5</td><td></td><td></td><td>11.2</td><td></td><td>[39]</td></tr><tr><td>12</td><td>R<sub>1</sub> = NO<sub>2</sub>, R<sub>3</sub> = OH</td><td></td><td>7.6</td><td>10.1</td><td></td><td>7.4</td><td>8.9</td><td>[39]</td></tr><tr><td>13</td><td>R<sub>1,5</sub> = NO<sub>2</sub>, R<sub>3</sub> = OH</td><td></td><td>4.4</td><td>9.1</td><td></td><td>4.7</td><td>7.4</td><td>[39]</td></tr><tr><td>14</td><td>R<sub>1,5</sub> = Cl, R<sub>3</sub> = OH</td><td></td><td>7.3</td><td>10.0</td><td></td><td>8.0</td><td>10.0</td><td>[39]</td></tr><tr><td>23</td><td>R<sub>2</sub> = OH (lawsone)</td><td></td><td>4.5<sup>(ii)</sup></td><td></td><td></td><td>4.1</td><td></td><td>[39]</td></tr><tr><td>24</td><td>R<sub>5</sub> = OH (juglone)</td><td></td><td>8.9</td><td></td><td></td><td>8.4</td><td></td><td>[39]</td></tr><tr><td rowspan="4">[Chemical Structure of 5-8]</td><td>15</td><td>R<sub>1</sub> = OH</td><td></td><td>9.6<sup>(ii)</sup></td><td></td><td></td><td>9.3</td><td></td><td>[39]</td></tr><tr><td>16</td><td>R<sub>2</sub> = OH</td><td></td><td>9.8<sup>(ii)</sup></td><td></td><td></td><td>9.7</td><td></td><td>[39]</td></tr><tr><td>25</td><td>R<sub>1,4</sub> = OH</td><td></td><td>9.4</td><td>10.9</td><td></td><td>10.0</td><td>10.4</td><td>[39]</td></tr><tr><td>17</td><td>R<sub>8</sub> = OH</td><td>5.0</td><td>9.8</td><td></td><td>5.4</td><td>9.4</td><td></td><td>[39]</td></tr><tr><td rowspan="4">[Chemical Structure of 9-11]</td><td>18</td><td>R<sub>4</sub> = Me, R<sub>8</sub> = OH</td><td>5.6</td><td>10.0</td><td></td><td>6.1</td><td>9.6</td><td></td><td>[39]</td></tr><tr><td>19</td><td>R<sub>3,4</sub> = Me, R<sub>8</sub> = OH</td><td>5.8</td><td>10.1</td><td></td><td>6.2</td><td>9.6</td><td></td><td>[39]</td></tr><tr><td>20</td><td>R<sub>2,4</sub> = Me, R<sub>8</sub> = OH</td><td>6.2</td><td>10.6</td><td></td><td>6.7</td><td>9.8</td><td></td><td>[39]</td></tr><tr><td>33</td><td>R<sub>2</sub> = OH</td><td></td><td>7.1</td><td></td><td></td><td>6.5</td><td></td><td>[41]</td></tr><tr><td rowspan="7">[Chemical Structure of 12-15]</td><td>39</td><td>R<sub>1,2</sub> = OH (alizarine)</td><td></td><td>6.8</td><td>11.1</td><td></td><td>5.6</td><td>10.1</td><td>[42]</td></tr><tr><td>40</td><td>R<sub>1,4</sub> = OH (quinizarine)</td><td></td><td>9.8</td><td>12.6</td><td></td><td>8.9</td><td>10.9</td><td>[43]</td></tr><tr><td>41</td><td>R<sub>1,8</sub> = OH (chrysazine)</td><td></td><td>8.3</td><td>12.5</td><td></td><td>8.8</td><td>11.9</td><td>[44]</td></tr><tr><td>42</td><td>R<sub>2,6</sub> = OH (anthraflavic acid)</td><td></td><td>7.4</td><td>8.7</td><td></td><td>6.5</td><td>6.8</td><td>[41]</td></tr><tr><td>53</td><td>R<sub>1,2</sub> = OH, R<sub>3</sub> = SO<sub>3</sub><sup>-</sup> (red S)</td><td></td><td>5.5</td><td>11.0</td><td></td><td>3.8</td><td>8.9</td><td>[40]</td></tr><tr><td>57</td><td>R<sub>1,5</sub> = OH, R<sub>2,6</sub> = SO<sub>3</sub><sup>-</sup></td><td></td><td>(iii)</td><td>9.2</td><td></td><td></td><td>8.3</td><td>[45]</td></tr><tr><td>58</td><td>R<sub>2,7</sub> = H</td><td></td><td>7.8</td><td>11.4</td><td></td><td>8.7</td><td>10.0</td><td>[46]</td></tr><tr><td>59</td><td>R<sub>2,7</sub> = SO<sub>3</sub><sup>-</sup></td><td></td><td>7.6</td><td>10.6</td><td></td><td>6.1</td><td>9.4</td><td>[47]</td></tr></tbody></table>

(i)R<sub>x</sub> = H unless otherwise stated; (ii) values in parentheses are average of literature values; (iii) a first deprotonation is expected from molecular structure but no experimental value is available.

## 2.6. Studied molecular structures

In Table 1, we summarize the list of molecules studied herein to predict their pKa (names in italics), their redox potentials (names in bold) and / or their Pourbaix diagrams (names underlined).

# 3. Results and discussion

## 3.1. Prediction of pK<sub>a</sub>

Though it performed well for monoacids (cf. Introduction), and also appeared to perform well for amino-acid and peptides that possess multiple acido-basic functional groups [38], to the best of our

knowledge, the ability of COSMO-RS to predict multiple pK<sub>a</sub> of species containing such multiple pH-sensitive phenol functions was never assessed systematically. Thus, in this section, we check the predictive power of COSMO-RS in that respect, for various quinones and hydroquinones. Table 2 reports experimental and predicted pK<sub>a</sub> for a set of molecules close to the target molecules. As can be seen in Fig. 2 (left), a good agreement is obtained with experimental data, with a root mean square error (RMSE) of 0.9 pKa units, in line with performance of state of the art general models [9,10].

It can be noticed that pK<sub>a</sub> of 2<sup>nd</sup> deprotonations are generally underestimated by about 1 unit. This corresponds to an over-stabilization from COSMO-RS of the deprotonated form compared to the protonated form. This is ascribed to the original parameterization of
<figure><img src="image_9.png" /><figcaption>**Fig. 2.** Left: Fit plots for pKₐ predictions. Errors for the overall set of deprotonations: R² = 0.89, RMSE = 0.9 pKₐ units; for first deprotonation: R² = 0.92, RMSE = 0.7 pKₐ units; for the second deprotonation: R² = 0.91, RMSE = 1.4 pKₐ units. Right: Fit plot for prediction of standard electrode potential at pH 0. Outliers of Figure 3 are circled to show improvement obtained from use of Eq. (23).</figcaption></figure>

<table><caption>**Table 3**
Experimental and calculated standard-state potentials for a series of quinones and anthraquinones. Additional calculations based on Eq. (23) are in parentheses.</caption><thead><tr><th rowspan="2">Structure</th><th colspan="2">Molecule<sup>(i)</sup></th><th colspan="2">E<sup>0</sup> (V vs. SHE)</th><th rowspan="2">Ref</th></tr><tr><th>N°</th><th>Substituents</th><th>Exper.</th><th>Calcul.</th></tr></thead><tbody><tr><td>1,2-benzoquinone</td><td>5</td><td>-</td><td>0.83</td><td>0.81</td><td>[20]</td></tr><tr><td>1,4-benzoquinone</td><td>6</td><td>-</td><td>0.70</td><td>0.70</td><td>[20]</td></tr><tr><td>1,2-naphthoquinone</td><td>21</td><td>-</td><td>0.55</td><td>0.51</td><td>[20]</td></tr><tr><td>1,4-naphthoquinone</td><td>22</td><td>-</td><td>0.47</td><td>0.44</td><td>[20]</td></tr><tr><td>9,10-phenanthraquinone</td><td>26</td><td>-</td><td>0.44</td><td>0.37</td><td>[20]</td></tr><tr><td></td><td>27</td><td>R<sub>all</sub> = H (anthraquinone)</td><td>0.09</td><td>0.09</td><td>[20]</td></tr><tr><td></td><td>35</td><td>R<sub>2</sub> = COOH</td><td>0.21</td><td>0.17</td><td>[15]</td></tr><tr><td></td><td>36</td><td>R<sub>2</sub> = COOMe</td><td>0.22</td><td>0.18</td><td>[15]</td></tr><tr><td></td><td>31</td><td>R<sub>1</sub> = NH<sub>2</sub></td><td>0.04</td><td>0.00</td><td>[15]</td></tr><tr><td></td><td>37</td><td>R<sub>2</sub> = NH<sub>2</sub></td><td>0.08</td><td>0.01</td><td>[15]</td></tr><tr><td></td><td>54</td><td>R<sub>1</sub> = NH<sub>2</sub>, R<sub>2,3</sub> = Me</td><td>-0.03</td><td>-0.02</td><td>[15]</td></tr><tr><td></td><td>55</td><td>R<sub>1</sub> = NH<sub>2</sub>, R<sub>2,4</sub> = Me</td><td>-0.08</td><td>-0.06</td><td>[15]</td></tr><tr><td></td><td>44</td><td>R<sub>1</sub> = NH<sub>2</sub>, R<sub>2</sub> = Et</td><td>-0.04</td><td>-0.01</td><td>[15]</td></tr><tr><td></td><td>32</td><td>R<sub>1</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.20</td><td>0.23</td><td>[15]</td></tr><tr><td></td><td>38</td><td>R<sub>2</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.19</td><td>0.18</td><td>[15]</td></tr><tr><td></td><td>48</td><td>R<sub>1,5</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.24</td><td>0.33</td><td>[15]</td></tr><tr><td></td><td>45</td><td>R<sub>1,8</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.21</td><td>0.24</td><td>[15]</td></tr><tr><td></td><td>46</td><td>R<sub>2,6</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.23</td><td>0.23</td><td>[15]</td></tr><tr><td></td><td>47</td><td>R<sub>2,7</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.23</td><td>0.25</td><td>[15]</td></tr><tr><td></td><td>56</td><td>R<sub>1,8</sub> = OH, R<sub>2,7</sub> = SO<sub>3</sub><sup>-</sup></td><td>0.12</td><td>0.10</td><td>[15]</td></tr><tr><td></td><td>43</td><td>R<sub>1</sub> = OH, R<sub>4</sub> =</td><td>0.16</td><td>0.12</td><td>[15]</td></tr><tr><td></td><td></td><td>R<sub>1</sub> = OMe, R<sub>4</sub> =</td><td>0.16</td><td>0.31 (0.16)</td><td>[15]</td></tr><tr><td></td><td>50</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>51</td><td>R<sub>1,4</sub> =</td><td>0.15</td><td>0.31 (0.11)</td><td>[15]</td></tr><tr><td></td><td>49</td><td>R<sub>1,4</sub> = OMe</td><td>0.15</td><td>0.31 (0.14)</td><td>[15]</td></tr></tbody></table>

(i)Rₓ = H unless otherwise stated

COSMO-RS, based on neutral compounds, and its pKₐ model, parameterized from first protonations only. The discrepancy here could come from an over-estimation of the hydrogen-bonding term due to polarization charge densities being too high at the position of the phenolate.

## 3.2. Prediction of reduction potential E⁰

Table 3 reports the experimental and calculated standard electrode potentials at pH = 0 for various quinones and anthraquinones. Overall, a

good predictive power is obtained by the present method, with a RMSE of 0.052 eV for a range of values between -0.1 and 0.9 eV. However, three molecules substantially increase the RMSE: 1,4-dioxy(-2-propenyl) anthraquinone, 1,4-di-methoxy anthraquinone and 1-methoxy-4-oxy (-2-propenyl) anthraquinone (cf. Fig. 3). For these molecules, E⁰ is overestimated by about 0.15 V.

These 3 molecules already appeared as outliers in a previous study [15]. At that time, authors proposed that this error appeared to originate from the neglect of explicit hydrogen bonding between water molecules
<figure><img src="image_10.png" /><figcaption>Fig. 3. Oxidized and reduced forms of (a) 1,4-di-methoxy anthraquinone, (b) 1-hydroxy-4-oxy-(-2-propenyl) anthraquinone and (c) 1,4-dioxy-(-2-propenyl) anthraquinone.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 4. Predicted Pourbaix diagrams of (a) anthraquinone 27 (b) anthraquinone-2,7-disulfonate 47, (c) 1,8-dihydroxyanthraquinone (chrysazine 41) and (d) 3,3'-(9,10-anthraquinone-diyl)bis(3-methylbutanoic acid) 52 vs. experimental data [20,48–50].</figcaption></figure>

in the ether groups vicinal to the redox active carbonyls. To account for the discrepancy, Kim et al. used additional QM/MM calculations. However, as stated in the previous section, COSMO-RS parameterization for pKa calculations allows for an empirical recovery of explicit hydration effects. Thus, taking advantage of this parameterisation, we carried out calculations of E⁰ based on Eq. (23) for the three outliers. This allowed to greatly reduce the error (cf. Table 3 and Fig. 2 (right)) without requiring to use a third-party software to run QM/MM calculations. A possible explanation for these specific molecules to appear as outliers is that the sterical hindrance produced by the side groups on position 1 combined with explicit waters of hydration on the ether

moiety leads the molecule to be harder to reduce than expected, causing lower E⁰ than expected.

This suggests that this auxiliary COSMO-RS method can be used instead of a QM/MM calculation to predict E⁰ more accurately, whenever explicit hydration of a vicinal function to a reactive function of interest may impact proton or electron transfer.

### 3.3. Prediction of Pourbaix diagrams

Pragmatically, to anticipate systematic deviation pKa predictions on higher ionized units, 1 pKa unit has been added to original
<table><caption>**Table 4**
pKa and reduction potentials predicted for analyzing the impact of incremental structural variations on the Pourbaix diagram.</caption><thead><tr><th rowspan="2">N°</th><th rowspan="2">Abbreviation</th><th colspan="3">Anthraquinone</th><th colspan="2">Redox</th><th colspan="5">Anthraquinol</th></tr><tr><th>pK<sub>a</sub><sup>1</sup></th><th>pK<sub>a</sub><sup>2</sup></th><th>pK<sub>a</sub><sup>3</sup></th><th>E<sup>0</sup></th><th>E<sup>∞(a)</sup></th><th>pK<sub>a</sub><sup>1</sup></th><th>pK<sub>a</sub><sup>2</sup></th><th>pK<sub>a</sub><sup>3</sup></th><th>pK<sub>a</sub><sup>4</sup></th><th>pK<sub>a</sub><sup>5</sup></th></tr></thead><tbody><tr><td>27</td><td>anthraquinone</td><td></td><td></td><td></td><td>0.10</td><td>-0.45</td><td>8.7</td><td>11.0<sup>(b)</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>28</td><td>AQ-1-OH</td><td>9.3</td><td></td><td></td><td>0.04</td><td>-0.70</td><td>5.8</td><td>10.3<sup>(b)</sup></td><td>17.6<sup>(b)</sup></td><td></td><td></td><td></td></tr><tr><td>33</td><td>AQ-2-OH</td><td>6.5</td><td></td><td></td><td>0.07</td><td>-0.63</td><td>8.8</td><td>10.2<sup>(b)</sup></td><td>11.2<sup>(b)</sup></td><td></td><td></td><td></td></tr><tr><td>30</td><td>AQ-1-COOH</td><td>3.1</td><td></td><td></td><td>0.20</td><td>-0.46</td><td>3.7</td><td>9.0</td><td>13.3<sup>(b)</sup></td><td></td><td></td><td></td></tr><tr><td>35</td><td>AQ-2-COOH</td><td>3.7</td><td></td><td></td><td>0.19</td><td>-0.38</td><td>4.9</td><td>8.8</td><td>10.9<sup>(b)</sup></td><td></td><td></td><td></td></tr><tr><td>32</td><td>AQ-1-S</td><td>-0.5</td><td></td><td></td><td>0.28</td><td>-0.41</td><td>-2.1</td><td>8.5</td><td>14.1<sup>(b)</sup></td><td></td><td></td><td></td></tr><tr><td>38</td><td>AQ-2-S</td><td>-0.2</td><td></td><td></td><td>0.35</td><td>-0.19</td><td>-1.5</td><td>8.4</td><td>10.9<sup>(b)</sup></td><td></td><td></td><td></td></tr><tr><td>29</td><td>AQ-1-NO2</td><td></td><td></td><td></td><td>0.27</td><td>-0.14</td><td>3.9</td><td>11.8</td><td></td><td></td><td></td><td></td></tr><tr><td>34</td><td>AQ-2-NO2</td><td></td><td></td><td></td><td>0.27</td><td>-0.10</td><td>4.8</td><td>10.0<sup>(b)</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>39</td><td>alizarin</td><td>5.6</td><td>11.0<sup>(b)</sup></td><td></td><td>-0.05</td><td>-1.12</td><td>5.5</td><td>10.2<sup>(b)</sup></td><td>11.5<sup>(b)</sup></td><td>20.6<sup>(b)</sup></td><td></td><td></td></tr><tr><td>40</td><td>quinizarin</td><td>8.9</td><td>11.9<sup>(b)</sup></td><td></td><td>-0.10</td><td>-1.20</td><td>5.3</td><td>8.2<sup>(b)</sup></td><td>17.2<sup>(b)</sup></td><td>21.8<sup>(b)</sup></td><td></td><td></td></tr><tr><td>53</td><td>Red S</td><td>-3.0</td><td>4.8</td><td>9.9<sup>(b)</sup></td><td>-0.02</td><td>-1.10</td><td>-1.6</td><td>5.2</td><td>9.3<sup>(b)</sup></td><td>11.2<sup>(b)</sup></td><td></td><td></td></tr><tr><td>52</td><td>DPivOHAQ</td><td>5.1</td><td>5.8</td><td></td><td>0.08</td><td>-0.47</td><td>5.1</td><td>5.9</td><td>9.6<sup>(b)</sup></td><td>10.9<sup>(b)</sup></td><td>19.6<sup>(b)</sup></td><td></td></tr></tbody></table>

(a) E^∞ values computed by applying eq. 8 were not used to construct the displayed Pourbaix diagrams but are provided for completeness

(b) predicted pKa increased by one unit to account for the strong ionic character of the resulting anion.

COSMOtherm predictions for all protonations beyond the first. We first present the predicted pKa of anthraquinone itself using this approach in Fig. 4 (a).

To obtain the Pourbaix diagram of Fig. 4 (a), we used predicted E^0 value of 0.09 V from eq. 10 as well as predicted pKa1 = 8.7 and pKa2 = 10.0 from COSMOtherm. To anticipate systematic error at the second deprotonation, we corrected pKa2 by 1 unit (i.e. pKa2,corrected = 11.0) as suggested in Section 3.1. There is excellent agreement between experimental data [48] and the computed Pourbaix diagram.

The Pourbaix diagram of the well-known anthraquinone-2,7-disulfonate (AQ-2,7-DS) was then predicted. As sulfonic acids are strong acids, their pKa is expected to be below 0 and is not considered, and all states of protonation and oxidoreduction contain two anionic side

groups. Predicted E^0 of 0.25 V, pKa1 = 6.1 and pKa2,corr = 9.4 + 1 = 10.4 were used to predict the Pourbaix diagram in Fig. 4 (b).

Good agreement with experimental data is also obtained, though the underestimation of the first pKa value results in a small overestimation of the reduction potential of about 0.05 V at high pH. Comparing Fig. 4 (a) and (b), It can be seen that the overall increase in the value of the reduction potential due to the introduction of sulfonate groups is correctly predicted.

To further assess the predictive power of our method, we applied it to the more challenging Pourbaix diagram of 1,8-dihydroxyanthraquinone, which contains additional hydroxy groups susceptible to deprotonation. This yields a potentially more complex Pourbaix diagram with multiple slope inflections, as found by Cao et al. [49]. To construct the Pourbaix

<figure><img src="image_12.png" /><figcaption>**Fig. 5.** Predicted Pourbaix diagrams of (a) 1-hydroxy anthraquinone (AQ-1-OH) 28 and 2-hydroxy anthraquinone (AQ-2-OH) 33, (b) anthraquinone-1-carboxylic acid (AQ-1-COOH) 30 and anthraquinone-2-carboxylic acid (AQ-2-COOH) 35, (c) anthraquinone-1-sulfonate (AQ-1-S) 32 and anthraquinone-2-sulfonate (AQ-2-S) 38, (d) 1-nitro anthraquinone (AQ-1-NO2) 29 and 2-nitro anthraquinone (AQ-2-NO2) 34.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 6.** Intramolecular hydrogen bond between the hydroxy and one of the redox active oxygens, for 1-hydroxyanthraquinone (left) and 1-hydroxyanthraquinol (right).</figcaption></figure>

diagram (Fig. 4 (d)), $E^0 = -0.032\text{V}$ was predicted from Eq. (10) and pKa were predicted (and corrected for poly-anions) for both oxidized ($pKa_1 = 8.8$, $pKa_{2,corr} = 12.9$) and reduced ($pKa_1 = 4.2$, $pKa_{2,corr} = 10.8$, $pKa_3, corr = 13.0$ $pKa_{4,corr} = 19.7$) forms by COSMOtherm.

The last challenging molecule that we targeted is the 3,3'-(9,10-anthraquinone-diyl)bis(3-methylbutanoic acid) (DPivOHAQ) [50]. This molecule contains both carboxylic acid and phenolic alcohol functional groups, which can both be deprotonated. Moreover, in fact, the available experimental Pourbaix diagram was measured for a mixture of regio-isomers. So, we computed the Pourbaix diagrams for the two regio-isomers. According to our calculations, the result was identical to the eye, with no significant differences in predicted $pKa$ and reduction potentials for both oxidized and reduced species. So, we assumed that a mixing rule was not needed in this case. We obtained $E^0 = 0.079\text{V}$ by Eq. (10) and pKa (corrected for poly-anions) for oxidized ($pKa_1 = 5.1$, $pKa_{2,corr} = 5.8$) and reduced($pKa_1 = 5.1$, $pKa_{2,corr} = 5.8$, $pKa_{3,corr} = 9.6$ $pKa_4, corr = 10.9$) forms by COSMOtherm. Overall, the agreement with experimental measurements is good, though our model slightly underestimates the reduction potential over the entire pH range for this species, of about 50 mV, which is in line with our RMSE in prediction for $E^0$.

Here again, a good overall quantitative agreement is obtained between predictions and experiments, and the decrease of the redox potential due to the addition of the two hydroxyl groups is reproduced by the calculations. Furthermore, the initial slope corresponding to a single proton exchange is correctly predicted. However, it seems that the first pKa of the oxidized form is underestimated, and in fact is experimentally above the second pKa of the reduced form, resulting in a flat line on the experimental Pourbaix diagram that is not reproduced by calculations. Finally, it can be noticed that experimental Pourbaix diagrams of anthraquinones could be impacted by dimerization reactions as well as the presence of relatively stable radical intermediates [51], which were not accounted for here. Thus, it seems that the simple thermodynamics of acid dissociation and reduction, excluding radicals, already affords good predictive power for Pourbaix diagrams of anthraquinones.

### 3.4. Impact of chemical groups on Pourbaix diagrams of quinones

In this section, we predict Pourbaix diagrams of a few model anthraquinones to evaluate the impact of side groups on electrochemistry of anthraquinones (cf. Table 4 for the corresponding $pKa$ and reduction potentials). As model functional groups, we include a strong electron attractor, the $NO_2$ group, a strong electron donor, the -OH group, and we also evaluate the impact of two common ionizable groups, which are typically used to confer water solubility, $SO_3^-$ and COOH. In all predicted Pourbaix diagrams, the upper line represents the thermodynamic limit that would be imposed by the common Ferricyanide/Ferrocyanide redox posolyte on the potential difference, and the lower line represents the thermodynamic limit that is imposed by the oxidation of water itself.

First, the impact of simple addition of an OH group is discussed (cf. Fig. 5 (a)). The position of the hydroxy group introduces a noticeable difference in the diagrams. Indeed, though a OH group in the 2nd

<figure><img src="image_14.png" /><figcaption>**Fig. 7.** Structures of quinizarin, red S and simpler comparable anthraquinones.</figcaption></figure>

position decreases the reduction potential at most pHs, and especially in basic media, this is not the case when the OH groups is at the 1st position. This can be attributed to the formation of an intramolecular hydrogen-bond, which stabilizes one of the protons in higher pH (cf. Fig. 6).

The OH group at the second position, since it takes its phenate form at a lower pH, has a strong electron-donating effect on the aromatic ring, making it harder to reduce the anthraquinone, making the reaction more efficient to store energy. Thus, introducing OH moieties in places where they are not stabilized through hydrogen bonding can be a convenient approach to decrease the reduction potential and optimize the negolyte performance of the anthraquinone. In addition, OH groups introduce ionization sites that can be beneficial to increase hydro-solubility at higher pH.

Carboxylic acid moiety (cf. Fig. 5 (b)) seems to slightly increase reduction potential (by about 100 mV) up to about pH 10. As COOH groups have pKa lower than 10, it seems that even carboxylates exhibit electron withdrawing on the anthraquinones in both reduced and oxidized forms, making it less efficient as a negolyte. The effect seems to be offset at higher pH, perhaps because, as the ring becomes more electron-rich due to deprotonation of phenates on the reduced form, the carboxylate group becomes saturated in electrons and cannot withdraw more of them from the ring. As this group in carboxylate form can promote water solubility, given a suitable posolyte, the relatively weak effect on reduction potential may be tolerable.

The common sulfonic acid moiety exhibits a difference between the 1-substituted and the 2-substituted forms (cf. Fig. 5 (c)). Indeed, though the 2-substituted form maintains a noticeable increase in reduction potential through the full pH range, by more than 100 mV, again suggesting an electron-withdrawing effect by the sulfonates, and making the anthraquinone intrinsically less efficient as a negolyte, this effect is halved for the 1-substituted form, and even removed at pH 13-14. The main cause is probably sterical hindrance between the electroactive oxygen and the sulfonate for the 1-substituted form relative to the 2-substituted one, destabilizing the oxidized form more than the reduced form (in which intramolecular hydrogen-bonding can compensate the effect).

A last, interesting case is the one of the nitro-group. As can be seen in Fig. 5 (d), a single nitro-group, through its large electron-withdrawing effect, produces a large increase of 300 mV on the reduction potential, regardless of the substitution position. This indicates that polynitro-substituted anthraquinones could be good candidates as RFB posolytes.
<figure><img src="image_15.png" /><figcaption>Fig. 8. Predicted Pourbaix diagrams of well-known quinizarine 40 and Red-S 53 compared with simpler anthraquinones.</figcaption></figure>

### 3.5. Comparison of predicted Pourbaix diagrams of well-known anthraquinones

For two important and cheap anthraquinones from the dye industry, quinizarin and Red-S [52], experimental Pourbaix diagrams were not available. As can be seen in Fig. 7, these two anthraquinones are multifunctional. With the new prediction method, it is interesting to study the predicted Pourbaix diagrams for these common anthraquinones, and to compare them with simpler structures, by gradual introduction of functional groups. Based on these comparisons, molecular design guidelines can be given.

In the case of anthraquinone, it can be noticed (Fig. 8 (a)) that the introduction of a second OH group, when compared to the 1-hydroxy-anthraquinone, further decreases the reduction potential over the entire pH range. Thus, the ddp improvement observed for the addition of one -OH group can be further magnified by the addition of other -OH groups. However, this decrease becomes more moderate starting at pH ≈ 8. Indeed, while the oxidized form conserves a strongly intramolecularly H-bonded proton until high pH (as confirmed by the calculated second pKa of 11.9), the reduced form can undergo two deprotonations before it has to break a strong intramolecular hydrogen bond for further deprotonation. Thus, though it might improve the stability of the electrolyte, intramolecular H-bonded alcohols tend to decrease the achievable ddp compared to alcohols which cannot generate such bonds.

When it comes to the Red-S (Fig. 8 (b)), promisingly, the calculated reduction potential is almost identical to the one of alizarin in the entire concentration range, despite the presence of a sulfonate group that, when alone in the molecule, substantially decreases the achievable ddps. This indicates that the enhancing effect of -OH groups outweighs the inhibiting effect of sulfonate groups. Therefore, the combination of -OH groups for ddp increase and sulfonate group for solubilization can be recommended for maximum performance of negolytes in redox flow batteries.

## 4. Conclusion

In this article, a new approach to predict Pourbaix diagrams in silica, based on DFT and COSMO-RS, is proposed. Based on good-quality predictions of pKa and a phenomenological correction for higher deprotonations (with RMSE of about 0.9 units), as well as reduction potentials

(with RMSE of about 36 mV), combined with Handerson-Hasselbach and Nernst equations, the predicted Pourbaix diagrams of four known anthraquinones are in good agreement with experimental data.

Encouraged by the predictive power of the method, we provide a systematic analysis of the effect of a few functional groups at different substitution positions. This yields guidelines for the use of substituted anthraquinones as either negolyte or posolyte, with consideration of the enhancement of their water solubility.

In the future, the impact of counter ions, potentially pairing with the quinones or increasing the solubility of the corresponding salts, on Pourbaix diagrams, could be considered. Recent improvements of COSMO-RS [53] to account for species with multiple charges should eliminate the need for a phenomenological pKa correction for these species in the near future. In addition, even though in principle the effect of temperature is automatically included by our approach within the applicability range of COSMO-RS, it could be interesting to verify if the predictions are equally reliable in the entire temperature range for the liquid state of water (0 to 100°C). To the end, the inclusion of radical reactions on Pourbaix diagrams could lead to refined predictions relevant to their reactions in batteries.

### CRediT authorship contribution statement

**Théophile Gaudin:** Conceptualization, Methodology, Software, Validation, Investigation, Writing – original draft, Writing – review & editing, Visualization, Data curation. **Jean-Marie Aubry:** Supervision, Writing – review & editing, Project administration, Conceptualization, Resources.

### Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This work was financially supported by the Agence Nationale de Recherche (Grant number: ANR-19-CE05-0012). We thank T. Godet-Bar for fruitful discussions.

## Appendix A. Calculation of the number of exchanged protons in a redox reaction for a general species with pH-sensitive oxidized and reduced forms

For any pH-sensitive species i, a ratio of the protonated and de-protonated forms at equilibrium can be calculated by the Henderson-Hasselbalch equation at any pH:
$$10^{\rho H - pK_a^i} = \frac{[A^-]}{[AH]} \qquad (A.1)$$

where [A-] is the concentration of the basic form, and [AH] is the concentration of the acid form.

Let us define i as the number of deprotonated sites on a species. For example, QH⁻ corresponds to i = 1, and Q²⁻ to i = 2. For a given oxidized or reduced species, the numbers of molecules Nᵢ of all protonation states in the solution are recursively related. Rearranging eq. A.1 yields:

$$10^{\rho H - pK_a^i} = \frac{[A^-]}{[AH]} = \frac{N_{i+1}}{N_i} \qquad (A.2)$$

$$\Leftrightarrow N_{i+1} = 10^{\rho H - pK_a^i} N_i$$

Generalizing eq. A.2 by recursive substitution yields:

$$N_{i+j} = N_i \prod_{k=i}^{j} 10^{\rho H - pK_a^k} \qquad (A.3)$$

Applying eq. A.3 for i = 0 yields:

$$\frac{N_i}{N_0} = \prod_{k=0}^{j} 10^{\rho H - pK_a^k} \qquad (A.4)$$

The proportions xᵢ of each protonation state i are defined as:

$$x_i = \frac{N_i}{\sum_j N_j} \qquad (A.5)$$

Dividing eq. A.5 by N₀ on numerator and denominator and inserting eq. A.4 into it yields:

$$x_i = \frac{\prod_{k=0}^{i} 10^{\rho H - pK_a^k}}{1 + \sum_{j}^{j} \prod_{k=0}^{j} 10^{\rho H - pK_a^k}} = \frac{10^{\sum_{k=0}^{i} (\rho H - pK_a^k)}}{1 + \sum_{j}^{j} 10^{\sum_{k=0}^{j} (\rho H - pK_a^k)}} \qquad (A.6)$$

Then, the expectation value of the number of deprotonated protons on a given species at a given pH is obtained (cf. eq. A.7).

$$\langle n_{\text{dep}} \rangle_{\text{pH}} = \sum_{i} x_i (pH) \cdot i \qquad (A.7)$$

And so, the number of liberated protons upon the redox reaction is straightforwardly calculated by counting the protons on oxidized and reduced forms and calculating the difference between the counts (cf. eq. A.8)

$$\langle n_{H^+} \rangle_{\text{pH}} = (n_{\text{lab}}^{\text{red}} - \langle n_{\text{dep}} \rangle_{\text{pH}}^{\text{red}}) - (n_{\text{lab}}^{\text{ox}} - \langle n_{\text{dep}} \rangle_{\text{pH}}^{\text{ox}}) \qquad (A.8)$$

where the superscript “red” denotes the reduced species, the superscript “ox” denotes the oxidized species, and nᵢₐ₉₅ denotes the total number of labile protons of a species.

## References

[1] H. Zhang, X. Li, J. Zhang, Redox Flow Batteries: Fundamentals and Applications, CRC Press, Boca Raton, 2018.

[2] R.F. Service, Advances in flow batteries promise cheap backup power, Science 362 (6414) (2018) 508.

[3] W. Wang, V. Sprenkle, Redox flow batteries go organic, Nat. Chem. 8 (3) (2016) 204–206.

[4] Y. Huang, et al., Nonaqueous redox-flow batteries: features, challenges, and prospects, Curr. Opin. Chem. Eng. 8 (2015) 105–113.

[5] Z. Yang, et al., Electrochemical energy storage for green grid, Chem. Rev. 111 (5) (2011) 3577–3613.

[6] C. Ding, et al., Vanadium flow battery for energy storage: prospects and challenges, J. Phys. Chem. Lett. 4 (8) (2013) 1281–1294.

[7] M. Zhang, et al., Capital cost sensitivity analysis of an all-vanadium redox-flow battery, J. Electrochem. Soc. 159 (8) (2012) A1183–A1188.

[8] K. Lin, et al., Alkaline quinone flow battery, Science 349 (6255) (2015) 1529–1532.

[9] C. Nieto-Draghi, et al., A general guidebook for the theoretical prediction of physicochemical properties of chemicals for regulatory purposes, Chem. Rev. 115 (24) (2015) 13093–13164.

[10] K. Mansouri, et al., Open-source QSAR models for pKa prediction using multiple machine learning approaches, J. Chemoinform. 11 (1) (2019) 60.

[11] T. Zhou, et al., Prediction of acid dissociation constants of organic compounds using group contribution methods, Chem. Eng. Sci. 183 (2018) 95–105.

[12] J. Ho, M.L. Coote, First-principles prediction of acidities in the gas and solution phase, WIREs Comput. Mol. Sci. 1 (5) (2011) 649–660.

[13] B. Thapa, H.B. Schlegel, Improved pKa prediction of substituted alcohols, phenols, and hydroperoxides in aqueous medium using density functional theory and a cluster-continuum solvation model, J. Phys. Chem. A 121 (24) (2017) 4698–4706.

[14] M.P. Andersson, et al., Prediction of aliphatic and aromatic oil-water interfacial tension at temperatures >100°C using COSMO-RS, Fluid Phase Equilib. 476 (2018) 25–29. Part A, 25 November 2018.

[15] H. Kim, T. Goodson, P.M. Zimmerman, Achieving accurate reduction potential predictions for anthraquinones in water and aprotic solvents: effects of inter- and intramolecular h-bonding and ion pairing, J. Phys. Chem. C 120 (39) (2016) 22235–22247.

[16] A.V. Marenich, et al., Computational electrochemistry: prediction of liquid-phase reduction potentials, Phys. Chem. Chem. Phys. 16 (29) (2014) 15068–15106.

[17] C.J. Cramer, Essentials of Computational Chemistry: Theories and Models, 2nd ed., John Wiley & Sons, 2004.

[18] V.S. Bagotsky, Fundamentals of Electrochemistry, 2nd ed., Wiley, Hoboken, 2006.

[19] M.J.S. Dewar, N. Trinajstic, Ground states of conjugated molecules-XIV: redox potentials of quinones, Tetrahedron 25 (18) (1969) 4529–4534.

[20] B. Huskinson, et al., A metal-free organic-inorganic aqueous flow battery, Nature 505 (7482) (2014) 195–198.

[21] R.P. Fornari, P. de Silva, A computational protocol combining DFT and cheminformatics for prediction of pH-dependent redox potentials, Molecules 26 (13) (2021) 3978.

[22] COSMOTHERM X. 2021, Dassault systèmes.

[23] F. Eckert, M. Diedenhofen, A. Klamt, Towards a first principles prediction of pKa: COSMO-RS and the cluster-continuum approach, Mol. Phys. 108 (3-4) (2010) 229–241.

[24] A. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38 (6) (1988) 3098–3100.
[25] J.P. Perdew, Density-functional approximation for the correlation energy of the inhomogeneous electron gas, Phys. Rev. B 33 (12) (1986) 8822-8824.

[26] A. Schaefer, H. Horn, R. Ahlrichs, Fully optimized contracted Gaussian-basis sets for atoms Li to Kr, J. Chem. Phys. 97 (1992) 2571-2577.

[27] COSMOconf 4.3, COSMOlogic GmbH, KG Co, (2018). http://www.cosmologic.de.

[28] TURBOMOLE V. 7.4, https://www.turbomole.org/. 2019.

[29] A. Klamt, G. Schuurmann, COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient, J. Chem. Soc. Perkin Trans. 2 (5) (1993) 799-805.

[30] A. Klamt, The COSMO and COSMO-RS solvation models, Wiley Interdiscip. Rev. Comput. Mol. Sci 8 (1) (2018) e1338.

[31] A. Klamt, Conductor-like screening model for real solvents: a new approach to the quantitative calculation of solvation phenomena, J. Phys. Chem. 99 (7) (1995) 2224-2235.

[32] A. Klamt, et al., Refinement and parametrization of COSMO-RS, J. Phys. Chem. A 102 (26) (1998) 5074-5085.

[33] K.S. Alongi, G.C. Shields, Chapter 8-Theoretical calculations of acid dissociation constants: a review article. Annual Reports in Computational Chemistry, Elsevier, 2010, pp. 113-138. R.A. Wheeler, Editor.

[34] A. Klamt, et al., First principles calculations of aqueous pKa Values for organic and inorganic acids using COSMO-RS reveal an inconsistency in the slope of the pKa scale, J. Phys. Chem. A 107 (44) (2003) 9380-9386.

[35] F. Eckert, et al., Prediction of acidity in acetonitrile solution with COSMO-RS, J. Comput. Chem. 30 (5) (2009) 799-810.

[36] F. Eckert, A. Klamt, Accurate prediction of basicity in aqueous solution with COSMO-RS, J. Comput. Chem. 27 (1) (2006) 11-19.

[37] J. Ho, et al., Theoretical calculation of reduction potentials, in: O. Hammerich, B. Speiser (Eds.), Organic Electrochemistry, CRC Press, Boca Raton, 2016, pp. 229-259.

[38] O. Toure, C.G. Dussap, A. Lebert, Comparison of predicted pKa values for some amino-acids, dipeptides and tripeptides, using COSMO-RS, ChemAxon and ACD/Labs methods, Oil Gas Sci. Technol. Rev. IFP Energies Nouvelles 68 (2) (2013) 281-297.

[39] Z. Rappoport, CRC Handbook of tables for Organic Compound Identification, 3rd ed., CRC Press, Boca Raton, 1967.

[40] ChemicalBook, https://www.chemicalbook.com/.

[41] H. Pal, T. Mukherjee, J.P. Mittal, Pulse radiolytic one-electron reduction of 2-hydroxy- and 2,6-dihydroxy-9,10-anthraquinones, J. Chem. Soc. Faraday Trans. 90 (5) (1994) 711-716.

[42] R.W. Sabnis, Handbook of Acid Base Indicators, CRC Press, Boca Raton, 2008.

[43] S. Rouhani, S. Salimi, Optical pH sensor based on quinizarin for alkaline pH regions, Prog. Color Color. Coat. 1 (2008) 11-17.

[44] L.M. Donald, Dictionary of Analytical Reagents, Springer, New York, 1993.

[45] J.C. Macdonald, J.H. Yoe, Spectrophotometric determination of scandium with anthrarufin-2,6-disulfonic acid (disodium salt), Anal. Chim. Acta 28 (1963) 264-270.

[46] M.A.B.H. Susan, et al., Effect of pH and the extent of micellization on the redox behavior of non-ionic surfactants containing an anthraquinone group, J. Electroanal. Chem. 481 (2) (2000) 192-199.

[47] R.J. Forster, J.P. O'Kelly, Protonation reactions of anthraquinone-2,7-disulphonic acid in solution and within monolayers, J. Electroanal. Chem. 498 (1) (2001) 127-135.

[48] Y. Jing, et al., A quinone anode for lithium-ion batteries in mild aqueous electrolytes, ChemSusChem 13 (9) (2020) 2250-2255.

[49] J. Cao, et al., A highly reversible anthraquinone-based anolyte for alkaline aqueous redox flow batteries, J. Power Sources 386 (2018) 40-46.

[50] M. Wu, et al., Extremely stable anthraquinone negolytes synthesized from common precursors, Chem 6 (6) (2020) 1432-1442.

[51] A. Weiss, Spectroelectrochemical Studies of the Reductive Intermediates of Anthraquinone-2,6-Disulfonate in Aqueous Media, Loyola University Chicago, Chicago, 2012.

[52] B. Dadpou, D. Nematollahi, Electrochemical oxidation of alizarin Red-S on glassy carbon electrode: mechanistic study, surface adsorption and preferential surface orientation, J. Electrochem. Soc. 163 (7) (2016) H559-H565.

[53] S. Müller, et al., Evaluation and refinement of the novel predictive electrolyte model COSMO-RS-ES based on solid-liquid equilibria of salts and Gibbs free energies of transfer of ions, Fluid Phase Equilib 483 (2019) 165-174.