When considering water oxidation catalysis theoretically, accounting for the transfer of protons and electrons from one catalytic intermediate to the next remains challenging: correction factors are usually employed to approximate the energetics of electron and proton transfer. Here these energetics were investigated using a closed system approach, which places the catalytic intermediate in a simulation box including proton and electron acceptors, as well as explicit solvent. As a proof of principle, the first two catalytic steps of the mononuclear ruthenium-based water oxidation catalyst [Ru(cy)(bpy)(H2O)]2+ were examined using Car-Parrinello Molecular Dynamics. This investigation shows that this approach offers added insight, not only into the free energy profile between two stable intermediates, but also into how the solvent environment impacts this dynamic evolution.
When considering water oxidation catalysis theoretically, accounting for the transfer of protons and electrons from one catalytic intermediate to the next remains challenging: correction factors are usually employed to approximate the energetics of electron and proton transfer. Here these energetics were investigated using a closed system approach, which places the catalytic intermediate in a simulation box including proton and electron acceptors, as well as explicit solvent. As a proof of principle, the first two catalytic steps of the mononuclear ruthenium-based water oxidation catalyst [Ru(cy)(bpy)(H2O)]2+ were examined using Car-Parrinello Molecular Dynamics. This investigation shows that this approach offers added insight, not only into the free energy profile between two stable intermediates, but also into how the solvent environment impacts this dynamic evolution.
Entities:
Keywords:
ab-initio molecular dynamics; density functional calculations; homogeneous catalysis; proton-coupled electron transfer; water splitting
In the optimisation of catalysts, computational techniques are increasingly employed to test proposed mechanisms and to infer catalyst design principles.1, 2, 3, 4, 5, 6, 7, 8 When evaluating proposed catalytic cycles, the free energies of the various catalytic intermediates are compared.2, 3, 4, 5, 6, 7 This comparison often uses correction factors to approximate the energetic contributions of the protons and electrons transferred during the catalytic step.3,4,6,9, 10, 11, 12, 13, 14, 15 Here the pitfalls of the commonly employed approach are highlighted, and an improved framework is suggested based on the Closed System Approach (CSA).16Consider, for example, the catalytic cycle with Proton‐Coupled Electron Transfer (PCET) reaction steps [Eq. (1)]shown in Figure 1. The PCET steps have a change in free energy [Eq. (2)]
Figure 1
The four PCET steps between the catalytic intermediates (I from I to I for water oxidation. Vertical lines denote electron transfer, horizontal lines proton transfer. Stable intermediates are shown in black. The ligand exchange I+2H2O→I+O2 is also shown.
The four PCET steps between the catalytic intermediates (I from I to I for water oxidation. Vertical lines denote electron transfer, horizontal lines proton transfer. Stable intermediates are shown in black. The ligand exchange I+2H2O→I+O2 is also shown.The free energy G has traditionally been approximated by [Eqs. (3) and (3)],following the methodology first proposed by Nørskov et al.13, 14, 15
,
and
are the enthalpy, zero‐point energy and entropic contribution respectively, as calculated for an intermediate in vacuum.
is the enthalpy calculated in the presence of a solvent, where the solvent is approximated by a continuous dielectric model. The problematic contributions in Equation (2) are the free energy changes due to proton and electron transfer. The solvation energy of a proton,
−11.45 eV,17 is significantly larger than the thermodynamic potential of water oxidation (1.23 eV). To overcome this issue, the energetic contributions of the proton and electron have often been approximated by combining the proton and electron into
H2,14,15 as the free energy of H2 is easier to compute.However, the
H2 approximation has its limitations: in the PCET reactions considered, H+ is always bonded to another ion. The transfer of H+ is a combination of the formation of one O−H covalent bond and the breaking of another. The proton transfer involves a rearrangement of the hydrogen‐bond network surrounding the catalytic site, which makes the simulation of such processes difficult.18,19 And yet, the complexity of the network rearrangement processes is ignored when using the
H2 approximation. This approximation may have been sufficient to establish trends in changes in free energy for different materials and catalysts,13,20 but as we move into the rational design of molecular catalysts, the how and why of the PCET process needs to be addressed. Two intermediates can no longer be seen as isolated from each other, instead we must also seek to optimise the processes between them.To address the details of the processes, including the thermodynamic barrier, between two intermediates, we employ the CSA. Within the CSA explicit water molecules and a metal ion (
) are included within the simulation box, to act as proton and electron acceptors respectively. In this way the charge carriers, and the processes via which they are transferred from the catalyst, can be considered explicitly. For the PCET reaction shown in Equation. (1), the equivalent equation within the CSA simulation box is Equation (5)where
denotes the solvated proton, which is part of a dynamic solvation complex.18, 19, 20 We can then decouple the PCET reaction into an electron‐ and proton‐transfer process.The electron transfer step is given by Equation (6)In the context of the CSA methodology, the energy needed to transfer an electron from the catalytic intermediate to the electron acceptor can be calculated by Equation (7)where
is the time‐averaged Kohn‐Sham (KS) energy over the simulation trajectory. One should note that
also includes the reorganisation energetic contributions resulting from the electron transfer.21 This includes contributions from internal vibrational and external solvent rearrangement. Various schemes for the calculation of redox properties based on Marcus theory have been previously implemented in the context of ab initio molecular dynamics simulations with explicit solvent.22, 23, 24 In these approaches only one redox active species is included explicitly in the simulation box, while in the CSA we include an explicit electron acceptor. Because the number of particles, charges, bonding patterns, and conformations of the reactants and products in Equation (6) remains the same, the change in entropy and zero‐point energy will be negligible, i. e.The proton‐transfer process [Eq. (8)]and the corresponding free energy change
may be investigated using constrained Car‐Parrinello Molecular Dynamics (CPMD) along the reaction coordinate of proton solvation.25, 26, 27, 28, 29 Constrained CPMD is required in order to sample relevant reaction pathways, which may be rare on the normal timescale of CPMD. Other enhanced sampling techniques, such as metadynamics,30 can also be employed to examine reaction pathways.31,32Within the CSA [Eq. (9)]where
includes energetic contributions from both the oxidation of the catalyst, as well as the reduction of the metal ion. The difference between
and
will be due to the energetic contribution of both the metal ion and solvated proton [Eq. (5)], a contribution that ideally should be a constant offset in the
throughout the catalytic cycle.Via this formalism, the way is paved for an energetic consideration of the process of a reaction step which includes both electron and proton transfer.16 Although this transcends the static consideration which uses the correction term
H2, it does introduce the extra complication of the energetic contribution due to the electron acceptor.To demonstrate this methodology the ruthenium based mononuclear molecular water oxidation catalyst (WOC) Ru‐bpy is used (see Scheme 1). Ru‐bpy provides an excellent test case as its catalytic cycle has been explored both experimentally and computationally using static methods.33 Here the first and second catalytic PCET steps of this WOC are examined within the CSA. The results of this examination are then compared to experimental data, as well as computational data obtained using static methods.
Scheme 1
Proposed catalytic cycle for water oxidation by [RuII(cy)(bpy)(H2O)]2+ (Ru‐bpy). Inset: Schematic structure of Ru‐bpy, which has a bipyridine (bpy) and a cymene (cy) moiety ligated to the Ru centre. Depicted in this way, the catalytic cycle consists of four consecutive PCET steps, i. e. along the diagonal in Figure 1, followed by the ligand exchange of O2 for H2O to regenerate the original intermediate.
Proposed catalytic cycle for water oxidation by [RuII(cy)(bpy)(H2O)]2+ (Ru‐bpy). Inset: Schematic structure of Ru‐bpy, which has a bipyridine (bpy) and a cymene (cy) moiety ligated to the Ru centre. Depicted in this way, the catalytic cycle consists of four consecutive PCET steps, i. e. along the diagonal in Figure 1, followed by the ligand exchange of O2 for H2O to regenerate the original intermediate.
Results and Discussion
The first two catalytic steps [RuII−OH2]2+→[RuIII−OH]2+ (I
1→I
2) and [RuIII−OH]2+→[RuIV=O]2+ (I
2→I
3) of the ruthenium‐based WOC were investigated within the CSA. The catalytic intermediates and electron‐ and proton‐transfer steps examined are shown in Scheme 2. Analogous to cyclic voltammetry, we first removed an electron from the catalyst (vertical arrows, Scheme 2), and then observed how the system responds. Proton transfer (horizontal arrows, Scheme 2) was investigated by means of constrained CPMD. The reaction coordinate examined during the constrained CPMD is the distance between the oxygen of an incoming water molecule and a proton on the oxygen ligated to the Ru centre (d(O→H), see also Figure 2).
Scheme 2
The electron‐ and proton‐transfer steps examined in this work. Vertical lines denote electron transfer, horizontal lines proton transfer. The regeneration step corresponds to the removal of an electron from the Mn ion, as well as the removal of the solvated proton. Intermediates observed to be stable in the simulations performed in this work are shown in black, unstable intermediates in grey. The inset shows the corresponding intermediates in Figure 1.
Figure 2
The distance constraints d(O→H), shown by the grey arrow, considered in this study for (left) [RuII−OH2]2++H2Osolv and (right) [RuIII−OH]2++H2Osolv.
The electron‐ and proton‐transfer steps examined in this work. Vertical lines denote electron transfer, horizontal lines proton transfer. The regeneration step corresponds to the removal of an electron from the Mn ion, as well as the removal of the solvated proton. Intermediates observed to be stable in the simulations performed in this work are shown in black, unstable intermediates in grey. The inset shows the corresponding intermediates in Figure 1.The distance constraints d(O→H), shown by the grey arrow, considered in this study for (left) [RuII−OH2]2++H2Osolv and (right) [RuIII−OH]2++H2Osolv.
Energetic Analysis of the PCET Step [RuII−OH2]2+→[RuIII−OH]2+
First the PCET step [RuII−OH2]2+→[RuIII−OH]2+ (cf. Scheme 2: [RuII−OH2]2++Mn3+⇌[RuIII−OH2]3++Mn2+⇌[RuIII−OH]2++H+
solv+Mn2+) is analysed. Considering the proton transfer process after electron transfer, [RuIII−OH2]3++Mn2+⇌[RuIII−OH]2++H+
solv+Mn2+, the variation of the time‐averaged constraint force
along the reaction coordinate is shown in Figure 3 (top). At d(O→H)=1.4 Å
is only slightly above zero. When d(O→H) is shortened,
increases. At d(O→H)=1.2 Å
is again very close to zero, indicating the transition state has been reached. At this point, the proton is equidistant between the two oxygen atoms: d(O→H)=1.2 Å vs. the average (1.20±0.06) Å for the unconstrained H−O distance. The dependence of
on d(O→H) is well in line with a normal activated reaction profile.
Figure 3
(top) The time‐averaged constraint force (
) as a function of d(O→H). This analysis is performed for [RuIII−OH2]3++Mn2+⇌[RuIII−OH]2++H+
solv+Mn2+. The error bars show standard deviations. The dotted line shows the fit of the calculated points. (bottom) The integral of the
fit with respect to distance. The definite integral has a value
=−0.17 eV (−3.8 kcal mol−1).
=0.05 eV (1.2 kcal mol−1)
(top) The time‐averaged constraint force (
) as a function of d(O→H). This analysis is performed for [RuIII−OH2]3++Mn2+⇌[RuIII−OH]2++H+
solv+Mn2+. The error bars show standard deviations. The dotted line shows the fit of the calculated points. (bottom) The integral of the
fit with respect to distance. The definite integral has a value
=−0.17 eV (−3.8 kcal mol−1).
=0.05 eV (1.2 kcal mol−1)At d(O→H)=1.1 Å the standard deviation of
is notably larger:
=(−1.1±0.9) eV Å−1. During the d(O→H)=1.1 Å simulation there is significant rearrangement in the solvent surrounding the reaction site, in preparation for proton diffusion into the solvent. This diffusion is facilitated by an appropriate ‘water wire’. We define a water wire to be a chain of water molecules which are hydrogen bonded together such that a proton may transfer rapidly along it. The proton can be considered to be delocalised along such water wires.34,35 The formation of a water wire has a large effect on λ. If a suitable wire has formed, the formation of the O→H bond is more favourable as the extra proton on the solvent molecule may diffuse away, and as a result λ is negative (see Figure S1). In contrast, the extra proton cannot be released from the solvent molecule if this water wire is broken; forming the O→H bond is therefore unfavourable. During these segments of the trajectory, λ is seen to be positive. After the simulation with d(O→H)=1.0 Å, the distance constraint is released and the system is allowed to evolve freely. The formed O−H bond has an equilibrium distance of 0.97 Å, and so
=0 at that distance, as shown in Figure 3 (top).The free energy profile of the complete proton transfer into solvent was obtained by integrating the fit of
and is shown in Figure 3 (bottom). This gives a
=(−0.2±0.1) eV (see Figure S2), and a transition state energy barrier
=0.05 eV (≈2 k
B
T at room temperature). Because
is less than zero we can conclude that proton transfer is thermodynamically favourable once the electron has been removed from the catalyst. An electron transfer energy
=(1.8±1.0) eV was obtained from the time‐averaged KS energy according to Equation. (7), with the distance constraint d(O→H)=1.6 Å. This distance corresponds to a typical hydrogen bond distance at equilibrium where
is around 0 (see Figure S3).
([RuII−OH2]2+→[RuIII−OH]2+)=(1.6±1.0) eV can then be calculated using Equation. (9) (see also Table 1).
Table 1
Summary of the changes in free energy
, in eV, for the first two PCET steps of the catalyst Ru‐bpy, as obtained by the CSA methodology, experiment (Exp),33 and static theoretical methods (Th).33 Δ(ΔG) is the difference between the ΔG obtained for the first two catalytic steps (top two rows of the table).The differences Δ1=CSA – Exp and Δ2=CSA – Th are also reported.
CSA
Exp
Δ1
Th
Δ2
ΔG([Ru−OH2]2+→[Ru−OH]2+)
1.6±1.0
0.67
0.93
0.87
0.73
ΔG([Ru−OH]2+→[Ru=O]2+)
2.1±0.8
1.27
0.83
1.38
0.72
Δ(ΔG)
−0.5±1.3
−0.60
0.1
−0.51
0.01
Summary of the changes in free energy
, in eV, for the first two PCET steps of the catalyst Ru‐bpy, as obtained by the CSA methodology, experiment (Exp),33 and static theoretical methods (Th).33 Δ(ΔG) is the difference between the ΔG obtained for the first two catalytic steps (top two rows of the table).The differences Δ1=CSA – Exp and Δ2=CSA – Th are also reported.
Proton Diffusion [RuIII−OH2]3+→[RuIII−OH]2+
Above it was concluded that proton transfer into the solvent is dependent on the proton at the reaction site having access to a viable water wire. During normal solvent dynamics these wires are formed and broken as water molecules move and rotate. The solvent environment during anticipated proton transfer was therefore examined, scanning to identify whether a solvent O atom was covalently bonded to a number of protons different from two i. e. whether it was an OH− or H3O+ ion. If a proton is within 1.2 Å of an O atom, it is considered to be covalently bonded.36 For the H3O+, the distance to the Ru catalytic centre, following the minimum‐image convention for periodic boundary conditions, was then plotted as a function of time (Figure 4). This analysis gives an overview of how the H+ moves through the solvent, visiting a number of different oxygens as it travels from the first to the third hydration shell.
Figure 4
The distance between Ru and H3O+, defined as an oxygen atom with 3 H within a radius of 1.2 Å, as measured for a sequence of constrained CPMD simulations. The CPMD simulations show diffusion of the single released proton for [RuIII−OH2]3++Mn2+⇌[RuIII−OH]2++H+
solv+Mn2+. The value of d(O→H) is noted in grey, and subsequent runs with decreasing d(O→H) are separated by dashed lines. According to the simulations, only one oxygen is in the H3O+ form at any time, and although the proton associates with a number of different oxygens (indicated with different colours) during the simulation, it is primarily bonded to three oxygens (blue, gold and magenta).
The distance between Ru and H3O+, defined as an oxygen atom with 3 H within a radius of 1.2 Å, as measured for a sequence of constrained CPMD simulations. The CPMD simulations show diffusion of the single released proton for [RuIII−OH2]3++Mn2+⇌[RuIII−OH]2++H+
solv+Mn2+. The value of d(O→H) is noted in grey, and subsequent runs with decreasing d(O→H) are separated by dashed lines. According to the simulations, only one oxygen is in the H3O+ form at any time, and although the proton associates with a number of different oxygens (indicated with different colours) during the simulation, it is primarily bonded to three oxygens (blue, gold and magenta).During the d(O→H)=1.2 Å simulation there is already an initial attempt at proton transfer from the first solvation shell at around 4 Å, to the second solvation shell at around 6 Å. Towards the end of the d(O→H)=1.1 Å simulation the proton resides mainly within the second solvation shell, though it is still shared with the first. At the start of the d(O→H)=1.0 Å simulation the proton has been transferred into the third solvation shell at 8 Å. After releasing the distance constraint (Unconstr, Figure 4), the proton remains in the third solvation shell for the 1 ps duration of the simulation, though it is mobile: the distance between its O atom and the ruthenium atom of the catalyst ranges between 6.5 and 8.9 Å.Of further interest here, is that the mobility of the water molecules themselves can also be observed (Figure 4). During the d(O→H)=1.1 Å simulation the proton is transferred to a water molecule in the second solvation shell (gold). During the simulation after d(O→H) has been released, at around 5 ps, the proton is again transferred to this water molecule, which has now moved into the third solvation shell.In the SI (Figure S3 and Figure S4) the [RuII−OH2]2++Mn3+⇌[RuII−OH]++H+
solv+Mn3+ proton transfer process is discussed, which would be a proton‐first pathway (e. g. I→I
as shown in the inset in Scheme 2, which is reproduced from Figure 1). This process does not yield a stable end product and is thermodynamically unfavourable.
Energetic Analysis of the PCET Step [RuIII−OH]2+→[RuIV=O]2+
In the reaction from [RuIII−OH]2+ to [RuIV=O]2+, the second PCET step (cf. Scheme 2: [RuIII−OH]2++Mn3+⇌[RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+),
remains negative, with a minimum at 1.1 Å (as shown in Figure 5, top). After the constrained dynamics with d(O→H)=1.0 Å, the constraint is again released, and the system is allowed to equilibrate. The newly formed O−H bond is stable, with an average bond distance of (0.970±0.009) Å.
Figure 5
(top) The time‐averaged constraint force (
) as a function of d(O→H) between one of the ligated water hydrogens and a solvent water oxygen for [RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+. (bottom) The integral of the interpolated
fit with respect to distance. The definite integral has a value
=−0.25 eV (−5.8 kcal mol−1).
(top) The time‐averaged constraint force (
) as a function of d(O→H) between one of the ligated waterhydrogens and a solvent wateroxygen for [RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+. (bottom) The integral of the interpolated
fit with respect to distance. The definite integral has a value
=−0.25 eV (−5.8 kcal mol−1).is integrated to give the free energy profile (Figure 5, bottom), with
=(−0.2±0.2) eV (see also Figure S2). The electron transfer energy is calculated from simulations of [RuIII−OH]2++Mn3+ and [RuIV−OH]3++Mn2+, where d(O→H)=2.04 Å. Using Equation. (7)
=(2.3±0.8) eV, and so, from Equation. (9),
([RuIII−OH]2+→[RuIV=O]2+)=(2.1±0.8) eV; this is also noted in the summarising Table 1.As in the case for the first catalytic step, for [RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+ the first attempt at proton transfer from the first solvation shell, at 4 Å, to the second solvation shell, now at around 5 Å, occurs in the d(O→H)=1.1 Å simulation (Figure 6). However, instead of moving into the third shell, it collapses back to the first solvation shell within the same d(O→H)=1.1 Å simulation. The extra proton remains shared between these two solvation shells, with brief transfers to the third shell. The water molecule in the second solvation shell receiving the proton gradually moves away from the reaction site, until, during the simulation without distance constraint (Unconstr, Figure 6), the additional proton stabilises in the second solvation shell (Figure 6, blue trace). At the end of this simulation, the additional proton is localised on a water molecule around 8.7 Å from the Ru catalytic centre.
Figure 6
The distance between Ru and the H3O+ ion associated with the proton transferred during [RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+. d(O→H) is noted in grey, and subsequent runs with decreasing d(O→H) are separated by dashed lines. The proton is primarily associated with an oxygen in the first solvation shell (brown trace), and interacts with, then transfers to, the oxygen with the dark blue trace. There are a number of momentary interactions between the proton and other oxygens; these oxygens each have a different coloured trace.
The distance between Ru and the H3O+ ion associated with the proton transferred during [RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+. d(O→H) is noted in grey, and subsequent runs with decreasing d(O→H) are separated by dashed lines. The proton is primarily associated with an oxygen in the first solvation shell (brown trace), and interacts with, then transfers to, the oxygen with the dark blue trace. There are a number of momentary interactions between the proton and other oxygens; these oxygens each have a different coloured trace.From the free energy profile shown in Figure 5 (bottom), it is clear that no thermodynamic barrier for proton transfer exists, once an electron has been transferred. Therefore, this process should occur spontaneously. To investigate this, the [RuIII−OH]2+→[RuIV=O]2+ step is also considered without an applied d(O→H).Unconstrained CPMD simulation of [RuIV−OH]3++Mn2+ shows proton transfer into the first solvation shell after 0.8 ps (Figure 7). In Figure 7 the relative positions of the O atoms bonded to a number of protons different from two are shown as the proton is released into the solvent, and the Ru−O bond contracts and stabilises. For clarity, this is decomposed into the various species present (Figure 7, bottom): H3O, and the O atom ligated to the Ru centre as first an OH, then O, ligand. After 0.8 ps, a proton is transferred from the ligated OH to the first solvation shell at 4 Å. This proton then interacts with two different water molecules in the second solvation shell. The Ru−O distance reaches equilibrium fairly quickly. This simulation encompasses the entire proton transfer into solution: [RuIV−OH]3++Mn2+⇌[RuIV=O]2++H+
solv+Mn2+.
Figure 7
(top) The distance between Ru and H3O+, OH− and O2− as measured during [RuIV−OH]3++Mn2+→[RuIV=O]2++H+
solv+Mn2+, where → indicates spontaneous proton transfer from the catalyst into the solvent. The ions are defined by the number of H within a radius of 1.2 Å from the oxygen, and are specified in (bottom): H3O – three H atoms, [Ru−OH] – one H atom, [Ru=O] – no H atoms. At first the oxygen ligated to the Ru centre has one proton bonded to it ([Ru−OH], green trace), but at around 0.8 ps a proton is transferred. The oxygen ligated to the Ru centre then has no bonded protons ([Ru=O], green trace), and the Ru – O distance is seen to shorten and stabilise. Meanwhile the excess proton is transferred to an oxygen in the first solvation shell (H3O, brown trace). It then continues to interact with three different oxygens within the solvent (brown, light green, pink).
(top) The distance between Ru and H3O+, OH− and O2− as measured during [RuIV−OH]3++Mn2+→[RuIV=O]2++H+
solv+Mn2+, where → indicates spontaneous proton transfer from the catalyst into the solvent. The ions are defined by the number of H within a radius of 1.2 Å from the oxygen, and are specified in (bottom): H3O – three H atoms, [Ru−OH] – one H atom, [Ru=O] – no H atoms. At first the oxygen ligated to the Ru centre has one proton bonded to it ([Ru−OH], green trace), but at around 0.8 ps a proton is transferred. The oxygen ligated to the Ru centre then has no bonded protons ([Ru=O], green trace), and the Ru – O distance is seen to shorten and stabilise. Meanwhile the excess proton is transferred to an oxygen in the first solvation shell (H3O, brown trace). It then continues to interact with three different oxygens within the solvent (brown, light green, pink).The difference in time‐averaged KS energies between the unconstrained calculations for [RuIV=O]2++H+
solv+Mn2+ (i. e. for the trajectory after 0.8 ps until the end of the simulation shown in Figure 7) and [RuIII−OH]2++Mn3+ is (2.7±0.7) eV. This encompasses both electron and proton transfer. However, in this trajectory bonds are formed and broken, and therefore one cannot assume that the change in zero‐point energy will be negligible. For this reason the vibrational density of states was obtained from the CPMD trajectory files,37,38 as shown in Figure S5. The zero‐point energy contribution from the changed Ru−O bond, the broken O−H bond, and the formed H+
solv/[H5O2]+ complex would amount to a correction of only −0.04 eV, which is small in comparison to the uncertainty in the change in KS energy. For the constrained case
([RuIII‐OH]2+→[RuIV=O]2+)=(2.1±0.8) eV. The two values are consistent within the large standard deviations, inherent to the thermal fluctuations for a system of this size. The differences between the two values will, in part, be due to the slightly different final configurations of the solvated proton.
Experimental Comparison and Evaluation
One of the dilemmas in comparing the values of
obtained so far to experiment is that the values of
contain a contribution from the metal ion electron acceptor and solvated proton. Taking the difference between ΔGs for two PCET steps [Eq. (10)]the energetic contributions due to the charge carriers will cancel out. This is based on the assumption that these terms are constant throughout the catalytic cycle. This provides a quantity that may be compared with experimental data, as well as data obtained from previous theoretical methodologies.33,39 The theoretical methodology used for comparison is one we have used earlier for this catalyst, which includes the use of the
H2 approximation.33 The changes in free energy for the two catalytic steps investigated here are summarised in Table 1, as are the resulting values for Δ(ΔG). For the CSA methodology Δ(ΔG)=(−0.5±1.3) eV. This is in very good agreement with experiment (see Table 1), allowing for the large standard deviation due to system size, as mentioned earlier. This is primarily due to the method of calculating
. Because of the fluctuations in KS energy during the simulations, a larger box size with more water molecules would be needed to decrease the standard deviation.The changes in free energy for each step as calculated with the CSA are compared to experimental and static theoretical values (see Table 1). Considering the static calculations, the Δ(ΔG) compares well with the CSA, where the static results were obtained using a different functional: B3LYP vs OPBE for the CSA. Using Δ(ΔG) will also lead to the partial cancellation of errors resulting from the use of the different functionals. Furthermore, there is also a good agreement between Δ(ΔG) when comparing the CSA to experimental values, within 0.1 eV. It may therefore be concluded that the initial proposition, that energetic contributions from the reduction of the metal and solvated proton should be a constant offset, does indeed hold. Future calculations can then use the calculated constant (Δ1, Table 1) to account for the use of the CSA with a Mn ion.It was mentioned in the introduction that the CSA ameliorates the need for correction factors in the case of water oxidation catalysis. Moreover, the good agreement between experiment and the CSA also leads us to conclude that it would be advantageous to apply the CSA in other contexts where relevant correction terms may not be known.
Conclusions
The energetics of the first two steps of the catalytic cycle of a ruthenium‐based WOC were examined using a closed system approach in which an electron acceptor is included in a fully solvated simulation setup. This setup allows for the independent examination of the energetic contribution of electron transfer, as well as the free energy profile and process of proton transfer. In both the steps examined, proton transfer was thermodynamically favourable after electron transfer. The first catalytic step has a small barrier of the order of k
B
T at room temperature, and the second step is barrier‐less. The closed system approach compares well with experiment, within 0.1 eV, though a large standard deviation results from the statistical uncertainty of the calculated electron transfer energies. Mechanistically, it was observed that a viable water wire is essential for proton release into the solvent, which further emphasises the importance of considering the environmental influence on a catalytic reaction.
Computational Method and Details
The CPMD program was used to examine the explicitly solvated systems.40 The solvent environment for the CPMD simulations was generated using Discovery Studio 2.5.41 The solvent was equilibrated for 0.2 ns using the CHARMM force field and CFF partial charge parameters at 300 K,42 while the catalyst was kept fixed. The volume was then adjusted using constant pressure for 0.2 ns, after which the system was further allowed to evolve with constant volume for 2 ns. Subsequently CPMD calculations were performed in the canonical NVT ensemble at 300 K, using GTH pseudopotentials for the transition metals,43 DCACP pseudopotentials for the remaining atoms,44 and the OPBE exchange‐correlation functional.45 This GGA functional has shown good performance when describing transition metal complexes.25,32,45, 46, 47 Kohn‐Sham (KS) orbitals are expanded on a plane wave basis set with an energy cut‐off of 70 Ry. A time step of 5 a.u.=0.121 fs was used. Image rendering for the CPMD output was done using VMD.48,49The general methodology for the CPMD simulations was as follows:The system was initially allowed to equilibrate with CPMD for 0.1 ps.A solvent water molecule was constrained at progressively closer distances to one of the protons of the water molecule ligated to the Ru centre (Figure 2), at 0.1 Å intervals. The system is allowed to evolve for at least 1 ps to allow the time‐averaged constraint force
to stabilise. In cases where stabilisation was difficult, simulation length was extended to a maximum of 2.5 ps.If d(O→H) was contracted to 1.0 Å, the distance constraint was then removed and the system allowed to evolve for 1 ps.If the formed bond is stable,
at this distance is set to 0.
is then fitted with a 100 pt Akima Spline,25,50 a fit based on cubic functions. Cubic‐based functions have long been used as a fitting method for
.26 This fit can then be integrated to give the free energy profile of the reaction.25, 26, 27, 28, 29
The Simulation Box
Simulations are performed on the catalytic intermediate and an Mn2+ or Mn3+ ion within a 17.52
15.78
13.65 Å3 box with 94 water molecules, total charge 5+. The water environment around the Mn ion was constrained to avoid spurious effects on
due to changes in the coordination sphere of the electron acceptor. The coordination sphere is stabilised by constraining the coordination numbers around the Mn ion. The number of oxygen atoms is constrained to four, based on unconstrained simulations where the Mn2+(H2O)4 structure was formed spontaneously (see Supplementary Information (SI), Figure S6 and S7). The coordination radius,
=2.25 Å was also determined from these unconstrained Mn simulations. At 2.25 Å the radial distribution function of O atoms around Mn was zero after the first solvation shell. The parameter
=9.4 Å−1, where
is the width of the transition region, is chosen such that
is significantly smaller than
.51 Although the four‐fold coordination is anomalous for Mn in water, it can be attributed to the relatively close location of the charged catalyst affecting the Mn ion. Mn chemistry has been shown to be very sensitive to the charge of the complex, as well as the surrounding environment.52 Each O atom was saturated with two H atoms using
=1.2 Å and
=17.6 Å−1. 1.2 Å may be considered the H−O distance at which a proton is equally shared between two oxygens,36 which makes it an appropriate
.
[RuII−OH2]2++Mn3+
In order to calculate
, a d(O→H)=1.6 Å was imposed for 2 ps for [RuII−OH2]2++Mn3+ (System 1 in Table 2). The multiplicity was then flipped to a septet to give [RuIII−OH2]3++Mn2+ (System 2 in Table 2), and the system allowed to evolve for a further 2 ps. For the calculation of
the system was allowed to evolve for at least 1 ps for each d(O→H): 1.4, 1.3, 1.2, 1.1, and 1.0 Å.
Table 2
Summary of the systems considered. qtot is the total charge of the system, 2S+1 the spin multiplicity, and Scat and SMn the spin of the catalyst and Mn ion respectively.
System
qtot
(2S+1)tot
Scat
SMn
1
[RuII−OH2]2++Mn3+
5
5
0
2
2
[RuIII−OH2]3++Mn2+
5
7
1/2
5/2
3
[RuIII−OH]2++Mn3+
5
6
1/2
2
4
[RuIV−OH]3++Mn2+
5
8
1
5/2
Summary of the systems considered. qtot is the total charge of the system, 2S+1 the spin multiplicity, and Scat and SMn the spin of the catalyst and Mn ion respectively.
[RuIII−OH]2++Mn3+
The initial geometry was obtained by removing the solvated proton once the final product had stabilised during the unconstrained
calculation for [RuIII−OH2]3++Mn2+ (System 2). The [RuIII−OH]2++Mn3+ system was equilibrated for 0.5 ps with sextet multiplicity (System 3 in Table 2). This system was allowed to evolve for 2.5 ps; then the multiplicity is switched to an octet and [RuIV−OH]3++Mn2+ (System 4 in Table 2) equilibrated for 2.5 ps. In order to calculate
the [RuIII−OH]2++Mn3+ system, with sextet multiplicity, was allowed to evolve for 2.5 ps with d(O→H)=2.04 Å, where d(O→H) was initiated from the equilibrated system. The multiplicity was then flipped to octet multiplicity [RuIV−OH]3++Mn2+ (System 4), and the system allowed to evolve for a further 2.5 ps. For the calculation of
d(O→H) was contracted from 2.04 Å: d(O→H)=1.9, 1.7, 1.5, 1.4, 1.3, 1.2, 1.1, 1.0 Å, where the initial contraction was more rapid (see Figure S8).
Conflict of interest
The authors declare no conflict of interest.As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.SupplementaryClick here for additional data file.
Authors: Ali Hassanali; Federico Giberti; Jérôme Cuny; Thomas D Kühne; Michele Parrinello Journal: Proc Natl Acad Sci U S A Date: 2013-07-18 Impact factor: 11.205
Authors: Aleksandr V Marenich; Abir Majumdar; Michelle Lenz; Christopher J Cramer; Donald G Truhlar Journal: Angew Chem Int Ed Engl Date: 2012-11-19 Impact factor: 15.336
Authors: James D Blakemore; Nathan D Schley; David Balcells; Jonathan F Hull; Gerard W Olack; Christopher D Incarvito; Odile Eisenstein; Gary W Brudvig; Robert H Crabtree Journal: J Am Chem Soc Date: 2010-10-21 Impact factor: 15.419
Authors: Roc Matheu; Mehmed Z Ertem; Jordi Benet-Buchholz; Eugenio Coronado; Victor S Batista; Xavier Sala; Antoni Llobet Journal: J Am Chem Soc Date: 2015-08-12 Impact factor: 15.419
Authors: Noam Agmon; Huib J Bakker; R Kramer Campen; Richard H Henchman; Peter Pohl; Sylvie Roke; Martin Thämer; Ali Hassanali Journal: Chem Rev Date: 2016-06-17 Impact factor: 60.622