H-abstraction from methane by iron(IV)oxo in water

In this section, we will examine the solvent effects on the energetics of the hydroxylation of methane by the ferryl ion in aqueous solution. We have performed constrained molecular dynamics simulations to compute the free energy barrier, using the method of thermodynamic integration (see e.g. ref. frenkel_smit). Since in vacuo the reaction barrier for the methane hydroxylation was found to be much higher in the methane coordination mechanism (22.8 kcal/mol) than in the oxygen-rebound mechanism (3.5 kcal/mol), the latter route is especially interesting to study in aqueous solution; also, because the required H$_2$O/CH$_4$ ligand exchange for the methane coordination mechanism adds an extra endothermicy of 23.4 kcal/mol in vacuo. On the other hand, one could argue that this ligand exchange reaction energy would be much less endothermic in aqueous solution, due to the hydrophobic repulsion between the methane substrate and the solvent in contrast to the hydrophilic interaction between the water molecule and the solvent. In principle, computation of the ligand exchange reaction (free) energy in aqueous solution is possible, although cumbersome and a complicated reaction coordinate would be required. Instead, using a simple bond constraint, the interaction between the methane molecule and the iron site can be computed and compared to the -6 kcal/mol attraction at an Fe-C distance of $R_{\rm {FeC}}=3.0$ Å in vacuo. We therefore performed a short constrained molecular dynamics simulation of this system (see e.g. the reactant complex structure in figure 8.1) in water, with the Fe-C bond distance constrained to 3.0 Å. The resulting average force of constraint is repulsive by 0.008 a.u., and similarly, an average force of constraint equal to 0.022 a.u. is found for a Fe-C bond distance fixed at 2.5 Å (see figure 8.4). The CH$_4$ molecule is thus repelled from the iron coordination site in aqueous solution, so that we have concluded that the methane coordination mechanism does not offer a viable route in aqueous solution.

Figure 8.4: The force of constraint (gray line) and a running average (black line) during a constrained AIMD simulation of [(H$_2$O)$_4$Fe$^{\rm {IV}}$O(CH$_4$)]$^{2+}$ in aqueous solution. Constraining the Fe-C bond length to 3.0 Å results in a repulsive force of 0.008 a.u. After moving the constraint to 2.5 Å (at $t\approx 4.25$ ps) increases the average to 0.022 a.u.

Proceeding with the oxygen-rebound mechanism, we performed 9 subsequent constrained molecular dynamics simulations at different points along the reaction coordinate of the reaction

$$[({\rm H}_2{\rm O})_5{\rm Fe}^{\rm {IV}}{\rm O}]^{2+} + {\rm ... ...})_5{\rm Fe}^{\rm {III}}{\rm OH}]^{2+} + \mbox{\.{}}{\rm CH}_3$$ (83)

in water. The cubic supercell with an edge of 10.13 Å and subject to periodic boundary conditions, contained one iron(IV)oxo molecule, one CH$_4$ molecule and a total number of 31 H$_2$O molecules. The total 2+ charge was balanced by a uniformly distributed counter charge. A Nosé thermostat maintained a temperature of $T=300$ K. The constrained reaction coordinate, $\xi$, was chosen to be a function of the O-H bond distance, $R_{{\rm OH}}$, the C-O bond distance, $R_{{\rm CO}}$, and the angle between these bonds, $\angle_{\rm HOC}$:

$$\xi = 1 - \frac{R_{\rm OH} \cos{\angle_{\rm HOC}}}{R_{\rm CO}}$$ (84)

The second term is the ratio of the projection of the O-H bond distance onto the C-O axis and the C-O bond distance, which equals $\xi =0.5$ when the hydrogen is in the middle of the oxo ligand and the carbon. For infinitely far separated reactants $\xi$ goes to zero, and for infinitely far separated products $\xi$ goes to one. Note that $\xi\geq1$ for $\angle_{\rm HOC}\geq90^\circ$, for instance when forming methanol and secondly, $\xi<0$ if $R_{\rm OH}>R_{\rm CO}$, for instance in freely rotating CH$_4$. For our 9 MD simulations, we choose $\xi=(0.30, 0.35, 0.40, 0.45, 0.50, 0.54, 0.60, 0.65, 0.70)$. Each simulation started by bringing the constrained reaction coordinate $\xi$ to the desired value during a short MD run, after which an equilibration MD run followed of at least 2.5 ps to adapt the system to the new constraint value, until there was no longer any drift in the fluctuating force of constraint. We then sampled the force of constraint, $\lambda_\xi$, during a 2 ps MD simulation. To correct $\lambda_\xi$ for the bias introduced by imposing the constraint, which limits the sampling to a constrained NVT ensemble, instead of the desired true NVT ensemble, the true force of constraint $F(\xi)$ was calculated using

$$F(\xi) = \frac{\langle\lambda_\xi Z^{-1/2}\rangle_{\xi} + \frac{1}{2}k_B T... ...{m_{i}} \nabla_i \xi \cdot \nabla_i Z \rangle_\xi}{\langle Z^{-1/2}\rangle_\xi}$$

$$Z = \sum_{i} \frac{1}{m_i} \left( \frac{\partial\xi}{\partial {\bf r}_i} \right)^2$$ (85)

where $T$ is the temperature, $k_B$ is Boltzmanns constant, $m_{i}$ and $r_{i}$ are the mass and the position of atom $i$ and $\lambda_\xi$ is the sampled Lagrange multiplier associated with the holonomic constraint $\xi-\xi^\prime({\bf r})=0$[17].

Figure: The Helmholtz free energy profile $\Delta A$ (solid thin line) for methane hydroxylation by the ferryl ion in aqueous solution, obtained from the sampled mean forces of constraint (circles and dashed line; right-hand-side axis) of 9 constrained AIMD runs. The discrepancy between the reaction energies of the hydroxylation in vacuo for the all-electron ADF calculation (crosses) and the approximate PAW calculation (squares), results in a correction to $\Delta A_\mathrm{aq, PAW}$, the corrected $\Delta A_\mathrm{aq}$ being given by the solid bold line (cf. text).

Figure 8.6: The most relevant bond distances as a function of the reaction coordinate $\xi$. The bond distances in vacuo and the (average) distances in solution are very similar, except for $R_{FeO}$. In vacuo, the Fe-O distance for [(H$_2$O)$_5$Fe$^{IV}$O]$^{2+}$ ($\xi =1$) is reached rapidly, as indicated by the horizontal dotted line, and similarly the Fe-O distance for [(H$_2$O)$_5$Fe$^{III}$OH]$^{2+}$ ($\xi =0$) at the other end. In aqueous solution the average Fe-O bond length is always significantly larger.

Figure 8.7: Solvation of the reacting species during the hydroxylation in water measured by the coordination numbers of the oxo/hydroxo ligand oxygen (upper graph) and of the methane/methyl-radical carbon and hydrogens (lower graph). Circles, squares and solid lines show the average number of solvent hydrogens with a certain radius around the atom and crosses and dashed lines show the number of oxygens with a certain radius around the atom.

The results, the forces of constraint for nine points along the reaction coordinate, are shown in figure 8.5, indicated by the open and closed circles (and the right-hand-side axis) and fitted with a quadratic spline (dashed line). Due to the limited size of the supercell, the reaction coordinate value of $\xi=0.3$ is practically the largest possible separation of the reactants in the box. The average distance between the oxo ligand and the methane carbon equals $R_{\rm {CO}}=3.44$ Å at this constraint value (see figure 8.6, for the relevant average inter-atomic distances, $R_{\rm {CO}}$, $R_{\rm {OH}}$, $R_{\rm {CH}}$, and $R_{\rm {FeO}}$ as a function of the reaction coordinate $\xi$). The average distance between the methane carbon and its hydrogen (that is to be transferred to the oxo ligand) is $R_{\rm CH}=1.114$ Å at $\xi=0.3$, which is almost equal to the other three average C-H bond distances (1.107 Å). The interaction between the reactants at $\xi=0.3$ is indeed very small, but not completely zero, as seen from the small value for the mean force of constraint at $\xi=0.3$ in figure 8.5. Moving along the reaction coordinate to larger $\xi$ results initially in an increasing repulsive force, which comes to a maximum at $\xi=0.477$ and then decreases again to cross the zero mean force of constraint at $\xi=0.580$, which marks the transition state. Initially, the reactants are just pulled towards each other, as shown by the rapid decrease in the averages $R_{\rm {CO}}$ and $R_{\rm {OH}}$ (see figure 8.6). But beyond $\xi=0.4$, also $R_{\rm {CH}}$ is seen to increase rapidly, indicating the transfer of the hydrogen from methane to the oxo ligand. The average distance between the oxo ligand and the methane carbon, $R_{\rm {CO}}$, has its minimum at the transition state position of $\xi=0.580$. During the constrained MD simulation at $\xi=0.65$ and $\xi=0.70$ (distinguished by the black filled circles in figure 8.5), the coordination number of iron dropped by one to five, as in both simulations one of the five water ligands left the coordination shell and moved into the solvent. This event did not however have a significant effect on the (running) force of constraint.

The Helmholtz free energy profile is obtained by integrating the mean force of constraint with respect to $\xi$:

$$A(\xi) = \int F(\xi)d\xi + C$$ (86)

The result is shown by the thin solid line in figure 8.5 (connected to the left-hand-side axis). We choose $C$ such that the free energy equals zero for our most left-hand-side point, $\xi=0.3$. The bold solid line shows the modified free energy profile, $\Delta A_{\rm {aq}}$, for the hydrogen abstraction from methane by the ferryl ion in water, after correcting for the discrepancy between results for the reaction energy of the hydrated complexes in vacuo obtained with the accurate all-electron ADF computation and with the approximate PAW approach, as explained in the method section 8.2. We have also plotted these reaction energies of the hydrated complexes in vacuo in figure 8.5 for comparison, using dotted lines and squares (labeled $\Delta E_{\rm {vac, ADF}}$) and crosses (labeled $\Delta E_{\rm {vac, PAW}}$). The reaction free energy barrier of 21.8 kcal/mol (with respect to the reactant complex at $\xi=0.3$) is much higher than the zero-Kelvin energy barrier found for the H-abstraction in vacuo of 3.4 kcal/mol. Including the temperature dependent internal energy, $\Delta E^\mathrm{300K}_\mathrm{int}=-4.8$ kcal/mol, and the entropy contribution, $-T \Delta S=2.5$ kcal/mol (see table 8.3) to this latter number, results in a free energy barrier of 1.1 kcal/mol in vacuo, making the difference in energetics even larger. The large difference between the free energy profile in aqueous solution and the reaction energy in vacuo is rather striking. Moving further to the product side, $\Delta A_{\rm {aq}}$ slowly decreases again by 3.3 kcal/mol, to 18.5 kcal/mol with respect to the reactant state at $\xi=0.3$.

Clearly, the aqueous solution has a strong effect on the energetics of the hydrogen abstraction reaction from methane by the ferryl ion. The hydration of the iron complex and the methyl by the surrounding dipolar water molecules plays an important role, as the charge distribution over the reactants changes along the reaction coordinate, leading to changes in the solvation energy. We have expressed the hydration of the oxo atom and the carbon and its (three) hydrogens in solvent coordination numbers, $cn$, by integrating the radial distribution functions of these atoms with respect to the solvent hydrogens and solvent oxygens up to certain radii. For the oxo/hydroxo ligand oxygen, the radial distributions, $g_{\rm {OH_s}}(r)$ and $g_{\rm {OO_s}}(r)$ showed clear peaks due to the hydration by, on average, one water molecule and the upper integration limit was therefore chosen at the minimum after the peaks at $r_{\rm {OH_s}}=2.25$ Å and $r_{\rm {OO_s}}=3.00$ Å. For the methane/methyl carbon and hydrogens, the hydrophobic solvation does not result in sharp peaks in the radial distribution functions, so that here the integration limits were chosen arbitrarily to be $r_{\rm {CH_s}}=3.0$ Å, $r_{\rm {CO_s}}=3.5$ Å, $r_{\rm {H_mH_s}}=2.5$ Å, and $r_{\rm {H_mO_s}}=3.0$ Å (the subscript s meaning solvent and m meaning methyl). Figure 8.7 shows these solvent coordination numbers for the oxo/hydroxo ligand (upper graph) and the methane/methyl substrate (lower graph).

The hydration of the hydrogen abstracting oxygen shows a clear picture: at the reactant side, the oxo ligand is coordinated by on average one water molecule and at the product side the hydroxo ligand also accepts on average one (1.2) hydrogen bond from a solvent molecule. In between, the hydration decreases to a minimum of about 0.4 hydrogen bonded solvent molecule at $\xi =0.5$. This decrease is partly explained by the approaching methane which disrupts the solvent structure around the oxo-ligand as the distance between the oxygen and the methane carbon, $R_{\rm CO}$ (figure 8.6) decreases. However, the oxo coordination is not so small at $\xi =0.55$ and $\xi=0.60$ even though $R_{\rm CO}$ is still small. Most likely, the increased charge transfer (polarization) in the transition state complex relative to the reactant complex which we observed for the hydrated complexes in vacuo (see table 8.4), also takes place during the hydrogen transfer in solution. As a result the oxygen becomes more negative in the transition state configuration (at $\xi=0.58$ in aqueous solution) and therefore more attractive for H-bond donating solvent molecules, hence the increased coordination number at $\xi =0.55$ and $\xi=0.60$, compared to $\xi=0.50$.

The coordination number of the carbon (and methyl hydrogens) initially increases, as the methane molecule becomes polarized when one C-H becomes slightly elongated (lower graph in figure 8.7). But again $cn$ drops when the distance between methane and the oxo ligand, $R_{\rm CO}$, and part of the solvation shell of methane has to make space for the oxo ligand to enter. Also here, the minimum in $cn$ is again followed by a significant increase which must be due to the increased (now positive) charge on CH$_3$, which we also noticed in vacuo (table 8.4).

Can the ferryl ion still be expected to be the active intermediate in oxidation reactions by Fenton's reagent? For the oxidation of methane in aqueous solution, the reaction free energy barrier of more than 20 kcal/mol for the rate determining step (i.e. the H-abstraction) is on the high side. However, as explained in the introduction, the reaction barrier for methane hydroxylation serves as an upper limit. The C-H bond energy in a more typical organic substrate, such as tartaric acid (as was used by H.J.H Fenton himself[128,155]) is in the order of 10 kcal/mol less than in methane. The remaining 10 kcal/mol for the barrier of the hydroxylation, is not very high at all, so that this result in combination with our previous work on the ferryl ion formation from Fenton's reagent leads us to conclude that the ferryl ion is indeed the most likely candidate for the active species in Fenton chemistry. Moreover, in most other solvents the screening effects will be less, which can now be expected to lower the barrier for the H-abstraction more.