Computational Analysis of Structure – Activity Relationships in Highly Active Homogeneous Ruthenium-based Water Oxidation Catalysts

Linear free energy scaling relationships (LFESRs) and regression analysis may predict the catalytic performance of heterogeneous and recently, homogenous water oxidation catalysts (WOCs). This study analyses twelve homogeneous Ru-based catalysts – some, the most active catalysts studied: the Ru(tpy-R)(QC) and Ru(tpy-R)(4-pic)2 catalysts, where tpy is 2,2:6,2-terpyridine, QC is 8-quinolinecarboxylate and 4-pic is 4-picoline. Typical relationships studied among heterogenous and solid-state catalysts cannot be broadly applied to homogeneous catalysts. This subset of structurally similar catalysts with impressive catalytic activity deserves closer computational and statistical analysis of energetics correlating with measured catalytic activity. We report general methods of LFESR analysis yield insufficiently robust relationships between descriptor variables. However, volcano plot-based analysis grounded in Sabatier’s principle reveals ranges of ideal relative energies of the RuIV=O and RuIV-OH intermediates and optimal changes in free energies of water nucleophilic attack on RuV=O. A narrow range of RuIV-OH to RuV=O redox potentials corresponding with the highest catalytic activities suggests facile access to the catalytically competent high-valent RuV=O state, often inaccessible from RuIV=O. Our work introduces experimental oxygen evolution rates into approaches of LFESR and Sabatier principle-based analysis, identifying a narrow yet fertile energetic landscape to bountiful oxygen-evolution activity, leading future rational design.


Introduction
The harvesting of sunlight offers a gateway into a sustainable energy future by providing a clean means to satiate the world's growing hunger for energy while providing a solution to the developing global climate change concerns [1]. The current challenge is in the harvesting and utilization of solar energy efficiently. Generation of solar fuels via artificial photosynthetic processes is an attractive method for harvesting this sunlight [2][3][4][5][6]. Numerous conceptual schemes have been proposed in which sunlight drives the flow of electrons and protons, leading to O-O bond formation and the evolution of molecular oxygen, O 2 ; the coupling of this proton and electron movement to water oxidation catalysts (WOCs) encourages the breaking of bonds in water and leads to formation of H 2 and O 2 [7,8]. The most significant bottleneck the entire catalytic process of light-driven water oxidation is oxygen-evolving reaction (OER) [9]: Current catalysts do not satisfy the need for cost efficiency, activity, and stability for this endergonic reaction; currently studied catalysts suffer from limited catalytic activity due to significant overpotentials (≥400mV). As such, research into the development and design of sufficiently stable and efficient catalysts capable of facilitating solar-driven water oxidation has grown significantly in recent years [10][11][12][13][14][15][16].
Historically, many catalysts were discovered through trial and error through synthesis and experimentation, though a more sophisticated approach adopts the mentality of rational design -development and design of novel catalysts with consideration of known catalytic theory to increase the likelihood of synthesizing an effective catalyst [17,18].
If computational theory and modelling may be employed to study and successfully characterize candidate catalysts and predict their catalytic competence, significant time can be spared on synthesizing and experimentally characterizing ineffective options. With the rapid improvement in the speed and capabilities of computational software and technologies, such an approach is becoming very alluring. Work has been done in the computational characterization of broad ranges of catalytic families and systems [19][20][21]. Popular are scaling relationships, employed to describe metal-organic frameworks [22,23], single-atom catalysts [24,25], and other heterogenous [26,27] and homogenous [28][29][30] systems. These scaling relationships relate parameters of a catalyst or its mechanism with predictors of strong catalytic activity, either oxidation rate, (low) overpotentials, turnover frequencies, or reduced energetic barriers in key mechanistic steps. The use of linear scaling relationships serves to reduce some of the systemic error [31][32][33] inherent to some computational methods, such as density functional theory (DFT), which may differ in their treatment of static correlation and electron localization. Since scaling relationships predict the relative activity of the catalysts, these systemic errors prove less significant in most cases. Despite the demonstrated need [34] for specifically chosen scaling relationships in homogenous catalysis based on ligand motifs, recent proposals still describe universal [35] scaling relationships for WOCs. Such reports are absent for the Ru(NNN)(NN)-based and similarly structured family of catalysts, Figure 1, home to some very highly active catalysts. Our computational analysis includes numerous Ru(NNN)(NN)-based catalysts, as well as Ru(bda)(N)(N) and several Ru(NNN)(QC) catalysts -some of the most active for Ru-based water oxidation catalysis. Analyzed complexes contain neutral polypyridine ligands alongside negatively charged QC (−1) and bda (−2) ligands, Figure 1.
The OER process occurs via a series of redox steps involving various reaction intermediates. The metal center coordinates water and is successively oxidized. Later, an O-O bond is formed, and molecular oxygen is ultimately released. Two primary paths have been proposed for O-O bond formation: water nucleophilic attack (WNA) and the interaction of two metal-oxo moieties (I2M) [17,36,37]; below, Figure 2 outlines a mechanism for a generic WOC, forming an O-O bond through a WNA process on the high-valent Ru V =O state.
Both methods of O-O bond formation are preceded by proton coupled oxidation, which comprises of three electron transfer events (ET) and two proton transfer events. Often, ET events are coupled to the transfer of a proton at the same time, in which case this proton-coupled electron transfer (PCET) allows for the reduction of energetic barriers during charge transfer. In the case of WNA, a metal-oxo species undergoes the attack of a solvent water molecule with a proton transfer, succeeded by further electron transfer prior to, or in concert with, release of molecular oxygen (details of these later steps were not investigated); I2M mechanisms simply involve coupling of two metal-oxo groups, resulting the O-O bond requisite for O 2 evolution. A substantial history of work on Ru-based WOCs [38][39][40][41][42] suggest that both WNA and I2M processes occur more easily upon reaching the high-valent Ru V =O state. Such studies describe the catalytic cycles leading to O 2 evolution in roughly four steps, each characterized by the removal of one electron: Ru III -OH Ru IV = O + H + + e − , Ru IV = O + H 2 O Ru III -OOH + H + + e − , and (4) Ru III -OOH Ru II + O 2 + H + + e − .
From the relative Gibbs energies in the reactions above, changes in free energies for each reaction step, ΔG i , leads to the determination of the most thermodynamically difficult step and by extension, the theoretical overpotential [17]: This implies that the ideal catalyst, with η th = 0, would have distributed overall change in free energy equally across each step in the reaction process; ΔG i =  [41,43,44]. At the same time some Ru-based WOCs boast very high rates of oxygen evolution. This suggests that the Ru V =O state might be achievable in some WOCs. Figure 2 posits an alternate pathway to this single-electron oxidation step from Ru IV =O. For Ru IV -OH species produced either by an ET step from Ru III or via protonation of the Ru IV =O species, a PCET step at a much more accessible potential (~1.4-1.6V) can be driven by sacrificial oxidants such as cerium(IV) ammonium nitrate (CAN) or electrochemically [45]. Numerous highvalent Ru V =O species of catalytically competent WOCs have been observed through use of electron paramagnetic resonance (EPR) [43,[46][47][48][49]. Ru V EPR features are rhombic and have g-tensors near (g xx ~2.1, g yy ~2.0, g zz ~1.9) [41], though one recent study report a significantly different, highly anisotropic EPR signal assigned to Ru V [50] which we suggest is better describe as a complex of Ru III with modified ligand.
Based on known mechanisms of WOC action, we tailor the generalized conventions of LFESR analysis -broadly applied to heterogeneous and solid-state catalysts -to a subset of highly active Ru-based, homogenous catalysts with Ru(NNN)(NN) or similar structures. Consideration of an additional ET/PCET event prior to O-O bond formation -necessary to yield the Ru V states is added in the analysis. These modifications reveal insights into the ideal range of relative energies of the Ru IV =O and Ru IV -OH intermediates, as well as a narrow region of optimal change in free energy of WNA on Ru V reactions: ΔG ~> −0.1eV. Finally, a narrow range of Ru IV -OH to Ru V =O redox potentials (~1.28V) corresponds to highest reported oxygen-evolution activity. These findings are based on empirical dataoxygen evolution rates and cyclic voltammetry redox couples -absent in prior Sabatierprinciple based works on similar families of catalysts.

Results
A schematic representation of the complexes considered in this study are shown below in Figure 1. Included in the set of model catalysts are those with neutral and negatively charged ligands. Each catalyst studied may be considered a deviation of another; either by substitution of an R group or similar modification can one catalyst be related to at least one other. This allows for comparisons of the precise effect of the ligand structure on catalytic activity and predictor relationships; these catalysts are chosen to enable both broad study of the Ru(NNN)(NN), Ru(NNN)(NO), and similarly structured families, as well as more narrow-scoped considerations of electronic structure effects. Table S1 shows the oxygen evolution rates as determined for each of the catalysts studied. The oxygen-evolution rate of the bda-type complexes is not considered for our analysis, as Ru(bda)(4-pic) 2 and Ru(bda) (isoq) 2 have second-order rate with catalyst concentration [51], unlike the other catalysts studied -first-order with catalyst concentration.
The twelve catalysts were modelled in each of the seven states of a typical WNA mechanism outlined in Figure 2. Each geometry was optimized, and free energies of each state were determined for each catalyst; see Methods section for computational details. The results of the DFT optimization and energy computation are reported in Tables S2-S4: the free energies of each intermediate, computed redox potentials, and changes in free energy of WNA processes. With such a large set of catalysts and species in our data set, we must validate our choice of basis set and computational methodology. This is done via comparison of computed and experimental redox potentials, Table 1. Due to systematic error inherent to DFT methods, agreement with experimental redox couples is ideally near/with ~0.2V agreement, the generally accepted margin of error for computed disagreement using similar basis set and computational methods [13,39,41,49,52]. However, due to the large number of catalysts and experimental evidence, some computed couples may have slightly larger differences; prominent examples of this include the Ru II /Ru III couples of the structurally similar Ru(EtOtpy)(bpy) and Ru(tpy)(bpy) catalysts. Both are predicted to advance to Ru III via a PCET reaction. Computation assigns a redox potential nearly 0.3V lower for these couples. However, the computed ET redox reaction for Ru(EtOtpy)(bpy) predicts a 1.21V potential for Ru II /Ru III , in closer agreement with experiment -0.98V. Computational prediction is equally dissimilar to the reported 1.04V potential for Ru(tpy)(bpy) -both are nearly 0.3V different (Tables 1 and S3). Table 1 cites experimentally reported redox potentials alongside computed redox couples for ET and PCET paths. Most of the computed redox potentials are within ~0.2V of the reported redox couples. This serves to validate the choice of basis set and the use of the generalized gradient approximation methodology, B3LYP, used in this study. Table 1 reports computational assignment of the specific redox couple at pH=0 for each catalyst's redox mechanism. Differences in the redox mechanism within catalyst of (NNN)(NN) structure suggest that specific ligand modifications sufficiently alter the electronic structure to overcome the effect of the inner coordination sphere's influence. Note that some computed redox couples, such as the Ru IV /Ru V of Ru(bda)(isoq) 2 , differ from proposed mechanism in identifying PCET/ET steps (Table 1 and S3) [49].
Having validated our choice of basis set and computational methodology, we proceeded with our adaption of conventional LFESR analysis to our set of homogenous, Ru-based catalysts outlined in Figure 1. Following prior methodologies for LFESR analysis [62], computed energies of the Ru III -OH, Ru IV =O, and Ru III -OOH intermediates (Table S5) relative to the stable Ru II species for each catalyst; these energies relative to the rest state (E RRS ) are plotted against the representative parameter of the redox potential of the PCET Ru II -H 2 O to Ru III -OH reaction in Figure 3.
Essential in the construction of volcano plots is establishing LFESRs, relating energetics of different catalytic intermediates onto a single variable [34]. In principle, different intermediates' stabilities are interrelated and cannot be significantly altered independently. Should a set of LFESRs be established for a particular intermediate, mathematical relationships derived therein may describe the entire energetic landscape in terms of a single descriptor. Plots of the (negative) change in free energy along the potential-determining step of a mechanism against a "descriptor" variable -for which the existence of LFESRs contain information on each catalyst's entire energetic landscape -indicate "ideal" catalysts; this is a volcano plot [34].
The free energy scaling relationships between the redox reaction of Ru II -H 2 O to Ru III -OH are moderately linear with the E RRS of Ru IV =O but do not have a strong linear relationship with the E RRS of the Ru III -OOH species resulting from WNA on Ru-oxo species ( Figure  3). Reports of LFESR-based analysis on catalytic systems often produce volcano plots using regressions with a high value for r 2 , usually above or near 0.90 [62,63]. Few analyses consider LFESR methods below r 2 of 0.8 [64]-which we will treat as a cutoff for sufficiently linear of a correlation. Reselection of the descriptor (x-axis variable in Figure  3) does not result in sufficiently strong correlations required for LFESR analysis; Tables S6  and S7 show the complete set of correlations of each descriptor against each E RRS . Though some descriptors, such as the Ru III -H 2 O to Ru IV -OH PCET reaction and the Ru IV =O to Ru V =O ET reaction, correlate well with early stages in the mechanism, they correlate more poorly with the WNA step. This E RRS for the Ru III -OOH intermediate does not correlate strongly with any descriptor studied (Tables S6 and S7). This step tends to have the weakest correlations with the descriptor variables, and its strongest correlation is with the Ru II -H 2 O PCET step, already dismissed above ( Figure 3). As such, these results indicate that this approach to LFESR analysis is inappropriate for the set of catalysts chosen; however, we will continue our study with Sabatier-principle and volcano-plot based analysis using empirical data of oxygen evolution rates -an approach which, to the best of our knowledge, is unpublished.
To continue our analysis, we consider alternative volcano plot / Sabatier principle-based analysis, using known oxygen-evolution rates as an experimentally derived measure of the quality of each catalyst. This approach to analysis continues in parallel with conventional LFESR-based analysis yielding volcano plots, using η th as the sole, theoretical analogue of catalytic effectiveness and quality.  Figure 4 shows the volcano-like plots produced from this approach.
Each graph in Figure 4 identifies a region of ideal catalytic activity based on a band of relative energies of intermediates in the WNA mechanism. Figure 4A explores a Sabatierbased relationship between the Ru IV =O to Ru V =O ET reaction often cited as needed for O-O bond formation -and oxygen evolution activity; easier access to the high-valent Ru V intermediate suggests high activity in O 2 evolution. Figure 4B supports the necessity of reaching Ru V =O prior to WNA reaction allowing for negative ΔG for the O-O bond formation. However, the catalysts with highest oxygen evolution rates are observed to experience a minimally downhill reaction in forming O-O bonds -ΔG ~> −0.1eV. Figure  4C suggests that ideal Ru-based catalysts have only a slightly thermodynamically uphill protonation from Ru IV =O to Ru IV -OH; unsurprisingly, this is similar to what is seen in Figure 4A. Figure 4D identifies a very narrow range of energies of the PCET Ru IV -OH to Ru V =O reaction, centered ~1.28V, where the four most active catalysts reside, surrounded by lower activity catalysts; the catalysts with lowest activity are located farthest from this region.
Having considered LFESR and Sabatier principle-based analysis for relationships with empirical oxygen evolution activity, our last approach attempts to identify linear regressions of various descriptors against oxygen evolution activity. The predictors in this approach consist of parameters specific to the species late in the WNA cycle -Ru IV , Ru V , and peroxo species. Table S8 reports the correlations of numerous computed descriptor variables against the experimentally reported oxygen-evolution activity and the theoretical overpotential of the catalytic system, Equation 6. Tables S9-S13 report the values for each of the parameters reported in Table S8. Charge and spin densities refer to the Milliken charge and spin densities as calculated by Gaussian16 [65]. Note that while computation of the Mulliken spin and density populations highly dependent upon basis set selection [66], the basis set is kept constant in this work. As such, the systematic bias inherent to basis set selection should not impact regression analysis. The four redox potentials / changes in free energies (Ru IV -OH to Ru V =O, Ru IV =O to Ru V =O, Ru IV =O to Ru IV -OH, WNA on Ru V =O) are discussed in greater depth later, as they suggest a Sabatier-like relationship with activity rather than with a linear regression-based model. The theoretical overpotential for most catalysts was taken as the very demanding ET reaction of Ru IV =O to Ru V =O proposed by most works. Thus, we consider only the strongest correlated descriptors of the catalytic mechanisms, the charge densities of the atoms most central to the Ru atom in the peroxide Ru III -OOH state. The charge densities on Ru and the adjacent oxygen atom are strongly inversely correlated with oxygen evolution activity and, suitably, by theoretical overpotential. Though some recent work [67,68] has been done on using linear regression-based methods to find suitable correlators with catalytic activity, a substantially greater amount of work uses Sabatier principle-motivated methods. Since most of the strongest correlators with activity are also discussed at length in the Sabatier-principle-based analysis, said analysis will compose the bulk of the discussion.

Discussion
Related to the employment of linear scaling relationships to find correlations and predictors of catalytic activity is the use of volcano plots, themselves rooted in Sabatier's principle. Sabatier's principle states that a catalyst should bind a substrate (water, or oxo group) neither too strongly nor too weakly. Volcano plots may predict catalytic performancevia overpotential or other thermodynamic/kinetic considerations -by consideration of a descriptor variable and linear free energy scaling relationships. Herein, this descriptor is related to estimate relative energies of other intermediates and states and, with processing, may lead to a "volcano" curve identifying predicted catalytic performance, with ideal catalyst candidates at the peak of the curve [17,26,28].
In 2018, Busch et al. [62] demonstrate the use of linear free energy scaling relationships to assess the viability of solid state and molecular OER catalysts, using a variety DFTbased functionals. The work included specifically use of generalized gradient approximation functionals (GGAs) and meta-GGA functionals. Using eight model catalysts, composed of a corrole ligand or a porphyrin equivalent, centered about a transition metal atom, Busch et al. were able to produce LFESRs among the E RRS of (* represents transition metal atom) *-OH, *=O, and *-OOH for various levels of theory for these model catalysts. From these LFESRs, Busch et al. were able to describe the thermodynamics of the WNA mechanism in relation to the overpotential of the potential-determining step; thus, they were able to determine a minimal systematic overpotential in terms of a representative parameter of these systems ---the E RRS energies of *-OH. The work did not consider mechanisms involving metal center in the formal oxidation state of (V) such as Ru V or Fe V intermediates.
We extended this LFESR-based approach to this broader group of catalysts with less ligand symmetry and, by extension, more ligand variability. The E RRS of the Ru III -OH, Ru IV =O, and Ru III -OOH intermediates are plotted against the representative parameter of the redox potential of the PCET Ru II -H 2 O to Ru III -OH reaction in Figure 3. These are based on the computationally derived energetics (Tables S3-S5). Regression against the Ru II PCET redox potential yielded r 2 of 0.78 and 0.60 for the E RRS of Ru IV =O and Ru III -OOH relative to the initial Ru II -H 2 O state. Other descriptor variables were considered to determine the existence of any sufficiently robust LFESRs (Tables S6 and S7), using the cutoff of r 2 near/above 0.80 as a lower bound of correlation. While some modest correlations were discovered, no mechanism could be constructed from pathways adequately linearly correlated against any individual predictor. As such, conventional LFESR analysis against E RRS is inappropriate for a wholly Ru-based system with even somewhat variable ligand structure. One report which employed such a technique in the past considered ligand structures exclusively a perfluoroporphyrin and corrole ligands centered on a transition metal atom [62]. This suggests that ligand differences are more highly impactful to the energetics of different intermediates; such differences are enough to change the catalytic activity by over two orders of magnitude, from Ru(tpy-Cl)(bpy) to Ru(tpy)(QC) and a 5x difference with substitution of an electron-donating EtO group in the even more structurally similar Ru(tpy)(4-pic) 2 to Ru(EtOtpy)(4-pic) 2 .
In 2019, Craig et al. [17] expand upon the groundwork established by Busch et al. prior.
Craig and coworkers consider seventeen molecular OER catalysts based on a variety of transition metals (Ru, Mn, Fe, Co, Ni, and Cu); ruthenium acts at the metal center in nine of these seventeen catalysts. They find that unmodified use of conventional scaling relations for heterogeneous systems does not accurately predict the catalytic activity of a set of homogenous catalysts; typical OER scaling relations addressing the OER for heterogenous systems neglects the ET reaction from Ru IV =O (or the PCET reaction from Ru IV -OH) yielding Ru V =O prior to WNA and O-O bond formation -Craig et al. briefly explore this by exploring the Ru IV =O to Ru V =O descriptor, discovering approximately half (five of nine) of their Ru-based catalysts -including Ru(bda)(4-pic) 2 and Ru(bda)(isoq) 2 -undergo ET for this reaction and have η th defined by this potential-limiting step. As mentioned previously, many Ru-based WOCs require access to Ru V prior to formation of an O-O bond since WNA processes on Ru IV =O tend to be thermodynamically unfavorable. Numerous catalytic reports indicate catalytic systems based on other elements, are predicted to form O-O bonds via WNA on Metal IV =O species [69][70][71]. As such, that our work focuses more on high-valent redox and WNA reactions is appropriate given our exclusive focus on the Ru metal center. Figure 4C shows the relationships between the oxygen-evolving activity and the protonation energy of Ru IV =O, often necessary for formation of the Ru IV -OH needed for PCET to Ru V . Most of the highly active catalysts show a similar trend, where a lower protonation potential to Ru IV -OH relates to higher oxygen-evolving activity. An exception to this trend is the Ru(EtO-tpy)(4-pic) 2 catalyst, found at 0.94V. This exception to the Sabatier-like trend suggests that catalytic dependence on this process is not as strong as the other relationships tested. Note that the Ru(tpy)(4-pic) 2 catalyst -most structurally like Ru(EtO-tpy)(4-pic) 2has a greater protonation activity alongside a lower oxygen evolution rate. A similar trend exists for the Ru(tpy-R)(bpy) family of catalysts as well. Taken as a whole, these data indicate that should Ru IV -OH be sufficiently accessible, oxygen evolution activity is not inhibited. Figure 4D indicates that a PCET redox potential of approximately 1.28V between Ru IV /Ru V is achievable and ideal for oxygen evolving activity. This region is surrounded by immediate and steep binding slopes lead to overly strong or weak catalyst-substrate interactions, inhibiting oxygen evolution activity. That Ru(tpy-Cl)(bpy) is weakly bound in Ru IV -OH is supported by energetics shown in Tables S3 and S5. This ideal binding region indicates that catalysts of the Ru(NNN)(NN) sand similar ligand structures may have facile access to the catalytically active Ru V species with a monotonically increasing series of redox potentials: Ru II /Ru III , Ru III /Ru IV , and Ru IV /Ru V . Figure 4B shows a similar result; oxygen-evolution activity may be inhibited by Ru V states which are too weakly bound. This characterizes most of the catalysts studied, excepting those with the highest catalytic activity, the Ru(tpy-R) (QC) family of catalysts. Ru(EtO-tpy)(4-pic) 2 offers only slightly weaker oxygen evolution activity and can be located partway along this "weakly binding slope" between the high and low activity regimes. Simply, in the absence of transition state energy considerations, the most catalytically active catalysts have minutely downhill energetics via WNA processes: ΔG ~> −0.1eV, implying a lower-energy, more accessible Ru V =O intermediate.

Materials and Methods
General Information: All chemicals and solvents were purchased from Sigma Aldrich, AK Scientific and TCI America and used as received. Aqueous solutions were prepared using ultrapure (Type 1) water (resistivity 18.2 MΩ·cm at 25 °C) from a Q-POD unit of Milli-Q integral water purification system (Millipore, Billerica, MA, USA). Solvents and chemicals were purchased from Sigma-Aldrich (St. Louis, MO, USA) and were used without further purification.

O 2 Evolution:
Oxygen evolution was measured with a Clark-type polarographic oxygen electrode with an Oxygraph System (Hansatech Instruments Ltd., King's Lynn, Norfolk, UK). Calibration was performed by measuring signal in O 2 -saturated deionized water and then adding sodium dithionite (an oxygen-depleting agent). The drop in the signal was set equal to the solubility of oxygen in water at room temperature (262 μmol/L). In Chemical catalysis the borosilicate vessel was filled with 500 μL solutions of the complex at pH = 1 (in 0.1 M nitric acid) and was constantly stirred. 20 mM. of Ce IV dissolved in nitric acid at pH = 1 was added to the chamber and oxygen concentration was recorded as a function of time.

DFT Calculation:
Density functional theory calculations were performed at the B3LYP level of theory, with the DGDZVP basis set for the Ru and Cl atoms, and the 6-31G* basis set was used for light (C, H, N, O) atoms. All molecules were modelled in water using the Conductor Polarized Continuum Model (CPCM) solvation model. All redox potentials were calculated using the DFT-calculated free energies of the products minus the reactants. From this value, 4.44 V was subtracted to account for the NHE voltage. The free energy of solvation for H+ was taken to be −11.64 eV. Geometries of intermediates are optimized, then electronic/thermal energies and vibrational frequencies are calculated as single-point calculations upon the optimized states. Optimization of some geometries required use of quadratic convergence criteria. Implicit consideration of water molecule in WNA methods with use of free energy of water molecule computed in 6-31G* at the B3LYP level of theory.

Pearson Coefficient:
The linear relationship between the two variables on a scatter plot is strongest when the points are closer to lying on a straight line. The Pearson coefficient r quantifies the strength of this linear relationship. For two variables x and y and data taken in n pairs, x 1 , y 1 , x 2 , y 2 , … x n , y n , Othe Pearson correlation coefficient is given by: Pearson correlation coefficient is given by: wherwhere x and y represent the means of the x and y variables, respectively. The strongest correlations have magnitudes of r closer to 1, and the value of the correlation coefficient is always between −1 and +1. Positive values indicate positive linear relationships; negative values represent negative linear relationships.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.