Modeling the Catalyst Activation Step in a Metal–Ligand Radical Mechanism Based Water Oxidation System

Abstract

Designing catalysts for water oxidation (WOCs) that operate at low overpotentials plays an important role in developing sustainable energy conversion schemes. Recently, a mononuclear ruthenium WOC that operates via metal–ligand radical coupling pathway was reported, with a very low barrier for O–O bond formation, that is usually the rate-determining step in most WOCs. A detailed mechanistic understanding of this mechanism is crucial to design highly active oxygen evolution catalysts. Here, we use density functional theory based molecular dynamics (DFT-MD) with an explicit description of the solvent to investigate the catalyst activation step for the [Ru(bpy)${}_2$(bpy–NO)]${}^{2+}$ complex, that is considered to be the rate-limiting step in the metal–ligand radical coupling pathway. We find that a realistic description of the solvent environment, including explicit solvent molecules and thermal motion, is crucial for an accurate description of the catalyst activation step, and for the estimation of the activation barriers.

1. Introduction

With a long-term increase in the demand for energy coupled with the wide spread use of fossil fuels, there is an urgent need to transition to efficient carbon-neutral fuels [1]. Sunlight-driven water splitting to produce molecular hydrogen is an attractive option in this regard. The overall water splitting reaction consists of the oxygen evolution reaction (OER) and the proton reduction reaction to form H${}_2$. The OER in particular suffers from a significant overpotential, thereby limiting the overall efficiency of this reaction [2]. A number of previous studies have proposed that the origin of this substantial overpotential is due to the non-ideal scaling relationship between *OH and *OOH, two intermediates of the OER [3,4,5]. While the ideal energetic separation between these intermediates should be 2.46 eV, it is found to be ∼3.2 eV for most metal oxides and a number of molecular catalysts [3,6]. Therefore significant efforts have been directed in the past few decades towards identifying highly efficient OER catalysts [2,7].

Among molecular catalysts, a number of mononuclear ruthenium based complexes have been reported to be highly active towards OER [8,9,10,11,12]. They can be broadly classified based on the operating mechanism during the crucial O–O bond formation step: the water nucleophilic attack (WNA) or radical-oxo coupling (I2M) as shown in Figure 1a,b. WNA involves the formation of a peroxo species (*OOH) resulting from the nucleophilic attack of a H${}_2$O molecule on a high-valent ruthenium oxo moiety (Ru${}^{(IV/V)}$=O). The I2M mechanism involves radical coupling of two high-valent oxo species, followed by the release of O${}_2$. Since this mechanism essentially involves the coupling of two radicals, it is expected to proceed with minimal activation barriers [13]. Additionally, it avoids the formation of the *OOH intermediate, rendering the above mentioned scaling relationship irrelevant [4]. Therefore, while it is highly desirable to have OER catalysts that operate via the I2M mechanism, there are only a few examples of ruthenium complexes that are known to operate via this mechanism, with the Ru(bda)L${}_2$ family being the most prominent examples [14,15] while other Ru based WOCs have also been reported [16]. The I2M mechanism involves the dimerization of the high-valent oxo species, and at the low concentration limit the reaction is essentially diffusion limited. Therefore, a crucial factor that decides the activity of these complexes is their affinity towards dimerization. Since this process involves a face-to-face approach of these moieties, the steric hindrance between them must be minimal in order to allow for easy coupling of the oxo moieties [17]. In addition, the solvent can also contribute to the barrier during this approach [18]. Naturally, a favorable interaction between the axial/equatorial ligands can help reduce the overall steric hindrance, thereby facilitating interaction of the ruthenium–oxos with minimal reorganization [14,15].

In a recent work, Pushkar and co-workers reported a mononuclear Ru based complex [Ru${}^{II}$(bpy)${}_2$(bpy–NO)]${}^{2+}$ (A) that is proposed to operate via a metal-ligand radical coupling pathway (MLC) during the O–O bond formation step [19]. This involves an intramolecular coupling between the metal and ligand generated radicals (Figure 1c), thereby eliminating the need for the dimerization of the oxo species in the right orientation and overcoming steric and solvation barriers, that is the main limitation of the I2M based complexes. They use isotope labeling studies and EPR measurements to show the involvement of the ligand in providing an oxygen during the O–O bond formation step. In addition, computational studies done in gas phase with implicit solvent corrections using the CPCM model show a negligible barrier for this step (<0.1 eV), via a ligand recoil mechanism involving an Ru${}^{IV}$=O species and the [ligand–NO]${}^{+}$ radical. However, the catalyst activation step that involves the dissociation of the ligand NO moiety and the binding of a water molecule on the metal center is most likely the rate-limiting step in this system, and only 10% of the reactant complexes ([Ru${}^{II}$(bpy)${}_2$(bpy–NO)]${}^{2+}$) enter the catalytic cycle as detected by EPR. In this regard, it is important to have an accurate understanding of the catalyst activation pathway, as it can provide insights on future design of highly active metal-ligand radical coupling WOCs.

In this work, we use an ab-initio molecular dynamics approach (DFT-MD) with an explicit description of the aqueous solvent to model the catalyst activation step for this system. The importance of incorporating thermal fluctuations and explicit solvent in order to account for the role of the solvent environment has been shown in a number of previous studies on organometallic complexes [20,21,22,23,24] including those on water oxidation systems [25,26,27,28,29,30,31] We find that the solvent actively participates along the catalyst activation pathway via hydrogen bonding interactions resulting in a higher barrier for this pathway, when compared to a gas phase model that lacks such interactions.

2. Results and Discussion

The catalyst activation step involves the dissociation of the NO ligand moiety that is bound to the metal center, with the adsorption of a (solvent) water molecule on the vacant metal site. To model this reaction in aqueous solution, we set up a system consisting of a reactant complex A and 156 water molecules in a cubic periodic box (L = 15.63 Å). The free energy profile of the activation step is determined using a biased sampling method. Here, we employed constrained molecular dynamics (CMD), with the distance Q between the metal and the oxygen of the ligand NO moiety (Ru–NO${}_{ligand}$) controlling the reaction pathway. In the reactant state, the Ru–O distance is ∼2.15 Å since the ligand NO moiety is bound to the metal center, while this distance is ∼3.75 Å close to the product state. Therefore, Q was varied between 2.15 Å and 3.75 Å. The average constrained forces and the resulting free energy profile in solution for this pathway is shown in Figure 2a,b. The computed activation barrier for this process is 29 kcal mol${}^{-1}$, with an estimated statistical sampling error of a few kcal mol${}^{-1}$. The barrier height indicates that the catalyst activation is a difficult step, consistent with the experimental observation that only 10% of the starting complexes enter the catalytic cycle. In order to quantify the effect of solvent, we also performed static DFT calculations of the corresponding catalyst activation pathway without explicit solvent molecules (gas-phase), and instead using an implicit solvent model for the gas-phase model. The activation barrier and the relative stability of the product state obtained is shown in Figure 2b (gas-phase). Our estimate for the gas-phase model shows a quantitatively different free energy profile with the barrier for the forward reaction significantly lower (19 kcal mol${}^{-1}$), and a higher stability of the product state as shown in Figure 2b ($\Delta$G = 12.5 kcal mol${}^{-1}$). Pushkar and co-workers report a value for $\Delta$G for the catalyst activation step in gas-phase (17.3 kcal mol${}^{-1}$), using a hybrid (B3LYP) functional with an implicit solvent model (CPCM) for water [19].

Representative snapshots of the initial, transition and product states during the catalyst activation pathway are shown in Figure 3. For clarity, the solvent molecules involved in hydrogen bonding interactions are highlighted. Initially, the ligand NO moiety forms one hydrogen bond on average with the solvent water molecules, as shown in Figure 3a. Close to the transition state corresponding to Q = 3.37 Å, a solvent water molecule (H${}_2$O${}_{substrate}$) approaches the metal center and binds to the Ru metal center in the process, as shown in Figure 3b. This provides an indication that as soon as the metal–NO${}_{ligand}$ bond is weakened, there is spontaneous adsorption of the water molecule on the (vacant) metal site. This also results in the metal–NO${}_{ligand}$ moiety forming upto three hydrogen bonds with the (solvent) water molecules (Figure 3b). At the product state (Q = 3.75 Å), the water molecule is strongly adsorbed on the metal center, and the ligand NO moiety is completely detached from the metal center, again forming upto three hydrogen bonds with the aqueous solvent, as shown in Figure 3c.

All these observations indicate a substantial involvement of solvent along the entire catalyst activation pathway. There is a significant rearrangement of the hydrogen-bond network near the reacting species, resulting in the product state showing very strong hydrogen bonding interactions of the ligand NO moiety with three water molecules, including the H${}_2$O${}_{substrate}$ that is involved in the reaction. This is reflected in the observed qualitative difference between the calculated free energy profiles of the gas phase and explicit solvent model, the most important characteristic being a substantially larger activation barrier ( 10 kcal mol${}^{-1}$) for the explicit solvent model. These observations also suggest that for a realistic description of catalytic water oxidation systems, simulation studies should account for the role of explicit solvent molecules.

In addition to the solvent effects discussed above, it is also evident that the NO${}_{ligand}$ moiety has a strong interaction with the Ru metal center in reactant complex A, due to a high barrier involved in the weakening of this bond in the catalyst activation pathway. In this regard, selective functionalization of the NO${}_{ligand}$ moiety, for instance, with electron-donating groups that can destabilize the metal–oxygen bond can potentially lead to a reduced barrier for the catalyst activation step, resulting in greater amounts of the starting catalytic complexes entering the catalytic cycle. This can be a topic of future high-throughput screening studies, [32] to help in the rational design of metal–ligand radical coupling based WOCs.

3. Materials and Methods

The reactions incorporating explicit water solvent were studied using DFT-based molecular dynamics (DFT-MD) with the Born-Oppenheimer approach, as implemented in CP2K [33]. The BLYP functional supplemented by Grimme’s D3 dispersion corrections was used [34,35,36] with GTH pseudopotentials for the nonvalance electrons (DZVP-GTH for Ru and TZVP-GTH for all other atom types) [37,38]. The auxiliary plane waves were expanded up to 280 Ry, and the system consisted of the Ru complex with 108 water molecules in a periodic cubic box (L = 15.63 Å). The temperature was controlled by a CSVR thermostat and set at T = 330 K [39] and a time step of 0.5 fs was used. The electronic structure was converged to an accuracy so that the energy is conserved within 0.6 kcal mol${}^{-1}$, for trajectories up to 15 ps. We used the constrained molecular dynamics (CMD) method to simulate the catalyst activation pathway and determine the associated free energy profile [40,41]. In this method, using a chosen reaction coordinate, Q, simulations are performed at several fixed values of this coordinate, ranging from Q${}_1$ to Q${}_2$. The free-energy difference ($\Delta$G) upon changing from Q${}_1$ to Q${}_2$ is then computed according to Equation (1). Here, is the average constraint force measured for each value of Q during the DFT-MD simulation. For each value of Q, a 3 ps equilibration run, followed by a 7–8 ps production run was performed. The data from production runs was used to compute the average forces and the resulting free energies, with the corresponding error bars.

$$\Delta G=-{\int }_{Q_1}^{Q2}〈F\left(Q\right)〉dQ$$

(1)

Gas phase model calculations of the barrier for the catalyst activation pathway were performed using the BLYP functional, with a TZVP basis set for all atoms, supplemented with D3 dispersion correction and implicit solvent corrections for water (COSMO) as implemented in the Orca package (Version 3.0.3) [42]. Hessian matrix calculations were performed to characterize all minima (no imaginary frequencies) and transition states (one imaginary frequency). Zero point energy and gas-phase thermal corrections (298 K, 1 bar, entropy and enthalpy) were obtained from these analyses.

4. Conclusions

In summary, we performed modeling studies on a mononuclear ruthenium based water oxidation catalyst that is proposed to operate via a metal–ligand radical coupling (MLC) pathway during the O–O bond formation step. Using a first-principles modeling approach, the importance of incorporating explicit solvent and thermal motion for an accurate description of the catalyst activation pathway was demonstrated. Possible directions for the design of active MLC based water oxidation systems was proposed.