Summary

Biochemical processes and environmentally friendly applied chemistry typically take place in aqueous solution. The solvent-solute interactions, particularly between charged or polar reactants and the water environment, are strong enough to have a dramatic impact on the chemistry taking place in solution. For example, the free energy barrier of the prototype S$_\mathrm{N}$2 reaction between Cl$^-$ and CH$_3$Cl, studied in chapter 3 is changed from 9 kcal/mol in the gas phase to 27 kcal/mol in aqueous solution, leading to a difference of 13 orders of magnitude (!) in the reaction rate at a temperature of $T=300$ K. Besides the effect of solvation, the solvent molecules can be actively involved in chemical reactions, for example via their ability to transfer protons and hydrogen atoms as was shown in e.g. chapters 5 and 6. This way, the solvent environment can completely change reaction mechanisms by opening new reaction pathways. Both types of solvent effects, we have studied on a variety of chemical reactions in water on a microscopic scale using the method of ab initio molecular dynamics (AIMD).

AIMD is a very powerful simulation technique as it combines electronic structure calculations by way of the density functional theory (DFT) method with classical dynamics of the nuclei, so that both chemical reactions and solvent dynamics can be accurately modeled. The performance of AIMD to model liquid water and ionic solutions such as HCl(aq) (chapter 3) and FeCl$_2$(aq) (chapter 5) was found to be very good for structural parameters such as radial distribution functions and coordination numbers and reasonably accurate for dynamic properties, e.g. the self-diffusion coefficient of bulk water and the acidity constant of [Fe$^\mathrm{III}$(H$_2$O)$_5$]$^{3+}$(aq) (chapter 4). The accuracy of the intra- and inter-molecular interactions was found to be in general in the order of 1 kcal/mol. Although the AIMD technique is without doubt the state-of-the-art of the computer simulation methods available for these systems, calculations on "chemistry in water" are still far from trivial, simply because chemical reactions typically occur on a much longer time scale ( $\approx10^{-9}-10^{-3}$ s) than the time scale accessible to AIMD ( $\approx10^{-12}-10^{-10}$ s), so that reactions during AIMD simulations are considered rare events.

A very effective technique to obtain information on the rare but important moment of transition state barrier crossing is the method of constrained molecular dynamics. With this technique, the system can be constrained to move only in the sub-space at a certain (fixed) point on a reaction coordinate. In chapter 3, we employed the method of constrained molecular dynamics to compute the free energy barrier of the textbook identity S$_\mathrm{N}$2 reaction between Cl$^-$ and CH$_3$Cl in aqueous solution. The solvent environment not only causes a dramatic increase of the barrier height, but affects the total free energy surface, which in the gas phase has the typical double-well shape. In aqueous solution however, endothermic solvent rearrangements approximately cancel the 11 kcal/mol exothermic reactant-complex formation, so that the free energy profile becomes unimodal. Our result for the free energy barrier of the S$_\mathrm{N}$2 reaction in aqueous solution of 27 kcal/mol is in surprisingly good agreement with the experimental number of 26.6 kcal/mol, considering we encountered two important error sources for which we had to correct. In the first place, we found that DFT with the Becke-Perdew exchange-correlation functional underestimates the transition state energy barrier by as much as 8 kcal/mol, which could be attributed to an artifact in the description of symmetrical "three-center four-electron" bonds, such as the $\sigma$-bond in Cl-C-Cl by present-day GGA exchange functionals. Secondly, the free energy profile showed a hysteresis leading to an error of 3.2 kcal/mol in the barrier height, which was traced back to result from a too slowly adapting solvent environment. We found strong indications that the rearrangements in the solvation shells of the attacking and leaving Cl's (the coordination numbers of which change from five and one respectively in the reactant state to two and two respectively in the transition state) are in fact activated transformations. Although the trivial solution to the hysteresis problem would be to perform (much) longer constrained molecular dynamics simulations, one should be aware that the computed potential of mean force is a function of the chosen reaction coordinate and always will result in too low a barrier if activated environmental changes are not described by the reaction coordinate. An important lesson is to always check on the amount of hysteresis by repeating the constrained molecular dynamics exercise in the reverse reaction direction. Nevertheless, we showed that we can very accurately compute the energetic and entropic effects of the solvent environment on a prototype chemical reaction, using AIMD.

After this instructive study on the elementary S$_\mathrm{N}$2 reaction, we devoted chapters 5 till 8 to the very versatile Fenton chemistry. The term Fenton chemistry denotes a broad range of oxidation and radical reactions initiated by the application of a mixture of iron(II) ions and hydrogen peroxide in aqueous solution, also known as Fenton's reagent. Although Fenton's reagent has found many applications and has been the subject of countless studies which have led to an even further broadening of the field by incorporating other transition metals, metal ligands and solvents, understanding of this type of oxidation catalysis was still clouded by many open questions concerning, most importantly, the reaction mechanism(s) and the nature of the reactive intermediate(s). The main experimental difficulty, the detection of the very short-lived reactive intermediates (possibly hydroxyl radicals or high-valent iron-oxo species), is not a problem for theoretical methods. Modeling of Fenton chemistry is nevertheless very challenging because it requires an accurate description of a number of complicated and simultaneously occurring processes, in particular, 1) the aqueous solvation dynamics, 2) the transport of protons and hydroxyl radicals along H-bond wires in the solvent, 3) the changing oxidation state, charge and spin-state of transition metal complexes, and 4) a number of chemical reactions, among which the O-O lysis of H$_2$O$_2$, hydrogen abstraction and hydroxylation of organic substrates, hydrolysis of metal ligands and radical reactions.

In chapter 5, we set out with the computation of the energetics of the two main reaction mechanisms proposed in literature for the reaction between iron(II) ions and hydrogen peroxide, starting from the [Fe$^\mathrm{II}$(H$_2$O)$_5$(H$_2$O$_2$)]$^{2+}$ complex in vacuo. We found that the more popular Haber and Weiss mechanism, in which the very reactive hydroxyl radical is formed, is endothermic by 21 kcal/mol. The alternative Bray and Gorin mechanism, producing the much contested ferryl ion ([Fe$^\mathrm{IV}$O]$^{2+}$) via an iron(IV)dihydroxo intermediate on the other hand is exothermic by 8 kcal/mol, with two minor reaction barriers, the highest being 6 kcal/mol. These preliminary static DFT calculations already revealed the importance of solvent effect, as e.g. the micro-solvation of a single H$_2$O molecule reduces the second reaction barrier from 18 to less than 4 kcal/mol. Using AIMD, we computed reaction pathways of the reaction between an iron(II) ion and hydrogen peroxide in aqueous solution, with two different starting configurations. Starting from H$_2$O$_2$ coordinated to the hydrated Fe$^{2+}$ ion in water, we observed the formation of the ferryl via a two-step mechanism, very similar to the static DFT results in vacuo. That is, first O-O bond lysis takes place, forming [Fe$^\mathrm{III}$(H$_2$O)$_5$(OH)]$^{2+}$ and an OH. radical. The OH. is very short-lived and travels rapidly via H-abstractions along a H-bond wire through the solvent to terminate at a water ligand of an iron complex in a neighboring periodic supercell, thus forming [Fe$^\mathrm{IV}$(H$_2$O)$_4$(OH)$_2$]$^{2+}$ (and a water molecule). The acidic iron(IV)dihydroxo species is seen to be in a dynamic equilibrium with its conjugate base, via hydrolysis taking place at the water ligands. In the second step, the ferryl ion is formed as hydrolysis takes place at one of the OH ligands, donating a proton to the solvent, which was not seen to return. In the second starting configuration, we included the coordination process by starting from separated but approaching reactants, H$_2$O$_2$ and [Fe$^\mathrm{II}$(H$_2$O)$_5$]$^{2+}$ (which has a vacant coordination site) in water. Again the ferryl ion was formed spontaneously, but now via a direct mechanism, as the intermediate OH. radical produced after coordination and simultaneous O-O lysis, directly abstracts the Fe$^\mathrm{III}$-OH hydrogen.

Visualization of the computed AIMD trajectories by movies, showing at a microscopic scale and in slow-motion the coordination and subsequent chemical reactions, gives a very clear impression of the mechanisms and relative time scales of these complex processes, which chemists denote with such abstract reaction equations as: Fe$^\mathrm{II}$(aq) + H$_2$O$_2$(aq) $\rightarrow$ Fe$^\mathrm{IV}$O$^{2+}$(aq) + H$_2$O(l). A disadvantage of our limited number of illustrative reaction pathways is of course the poor statistics, which do not allow e.g. to predict which of the two observed mechanisms (i.e. two-step vs. direct mechanism) is more favorable, let alone to compute a reaction rate. Another problem is the construction of the initial configurations of the systems, which in all cases involved unphysical geometric constraints, which could make the simulations less representative. In chapter 7, we attempt to tackle these problems, by generating new reaction pathways using the transition path sampling technique. Starting from the computed trajectory in which the separated reactants coordinate and react via the direct mechanism, we generated two new reactive pathways by randomly choosing a snapshot of the old trajectory, making small random changes to the atomic momenta and integrating the equations of motion forward and backward in time. On connecting the reactant state and the product state, the new pathways were accepted and used to branch of again two new pathways and so forth. This way, we computed two sequences of both 10 reactive pathways in order to obtain relaxed pathways, which no longer have a memory of the artificial construction of the initial pathway. In both sequences of pathways, the reaction mechanism changes from the initial direct mechanism to the two-step mechanism. In the first new pathways following the two-step mechanism, the OH. radical jumps via two or three solvent H$_2$O molecules to a neighboring periodic supercell to termination, forming the iron(IV)dihydroxo intermediate as already seen when we started from the [Fe$^\mathrm{II}$(H$_2$O)$_5$(H$_2$O$_2$)]$^{2+}$(aq) complex. But in the last generated pathways, the OH. radical terminates at the same iron complex as where H$_2$O$_2$ coordinated (instead of a periodic image) via a very short H-bond wire involving only one H$_2$O molecule. This trend was rationalized in chapter 7 from the fact that, along the sequences of reactive trajectories, the solvent relaxes in particular around the reactants as it loses the memory of the artificial preparation of the starting pathway. The H-bond wires in the solvent which are formed during the solvent relaxation are followed by the OH. radical after O-O lysis, which makes the two-step mechanism more favorable than the direct mechanism.

Having established the reaction mechanism and the formed reactive species for the reaction between iron(II) and hydrogen peroxide in water, we focused on the mixture of iron(III) and hydrogen peroxide (known as the Fenton-like reagent) as iron(III) can be formed in Fenton chemistry, and the Fe$^{3+}$/H$_2$O$_2$ mixture is known to be also capable of oxidizing organic substrates, albeit less reactively than Fenton's reagent. In chapter 6, we show that O-O lysis of hydrogen peroxide coordinated to pentaaquairon(III) is energetically unfavorable in contrast to the iron(II) case. The formation of a Fe(III)-OOH species by donation of H$^+$, as proposed in literature cannot realistically be modeled by static DFT without the proper inclusion of the solvent, which screens the strong electrostatic interactions during this charge separation process. In a number of AIMD simulations of Fe$^{3+}$ and H$_2$O$_2$ in aqueous solution, we indeed observe the spontaneous formation of the iron(III)hydroperoxo species. The DFT calculations show that homolysis of the Fe-O bond and the homolysis of the O-O bond are likely candidates for a second step in the Fenton-like mechanism, in particular favoring the latter in which both a ferryl ion and a hydroxyl radical are formed, when hydrolysis of the iron complexes is taken into account. Still, the reaction energy of this O-O homolysis of [Fe$^\mathrm{III}$(H$_2$O)$_4$(OH)(OOH)]$^{+}$ is +26 kcal/mol in vacuo. In aqueous solution, we have computed the free energy barrier of the O-O homolysis to be $\Delta A^\ddagger = 21$ kcal/mol using the method of constrained molecular dynamics. Hydrolysis is indeed seen to play an important role as we observe the increasing acidity of the iron complex during the reaction. Analysis of the vibration spectra of the iron(III)hydroperoxo species in water confirms a reduced O-O bond strength and an increased Fe-O bond strength when the high-spin ($S=5/2$) iron complex is brought to the low-spin ($S=1/2$) state, as suggested by Raman spectroscopy studies. Using metal ligands that induce the low-spin state, we expect therefore a further reduction of the O-O homolysis endothermicity.

The observed spontaneous formation of the ferryl ion from Fenton's reagent alone is no prove that the ferryl ion is indeed the active intermediate capable of oxidizing organic substrates. In chapter 8, we have therefore studied the oxidation of methane to methanol by the ferryl ion in aqueous solution. We have investigated in particular two reaction mechanisms suggested in the literature: 1) the methane-coordination mechanism, in which the methane molecule first coordinates to the metal and transfers a hydrogen to the oxo ligand forming the [Fe$^\mathrm{IV}$(H$_2$O)$_4$(CH$_3$)(OH)]$^{2+}$ intermediate, after which the CH$_3^-$ ligand is transferred from the metal to the hydroxo ligand, forming methanol in the second step, and 2) the oxygen rebound mechanism, in which the ferryl ion first abstracts a hydrogen from the methane molecule forming an pentaaquairon(III)hydroxo intermediate and the CH$_3$. radical which rebounds in a second step back to the oxygen forming the methanol. Our DFT calculations strongly favor the oxygen rebound mechanism, as the reaction energy barrier for the hydrogen abstraction from methane is only 3 kcal/mol compared to 23 kcal/mol in the methane coordination mechanism. Moreover, the water-methane ligand exchange required for the methane coordination mechanism adds an endothermicity of 23 kcal/mol in vacuo. Overall, the oxidation of methane by the ferryl ion is exothermic by 47 kcal/mol. Our results show strong similarities with the methane oxidation by enzymes, such as methane mono-oxygenase and cytochrome P450. Inclusion of the full solvent environment again shows the importance of these effects when dealing with charged species. Using constrained AIMD, we estimated the free energy barrier for the H-abstraction in the oxygen rebound mechanism to be 22 kcal/mol, much higher than in the gas phase.

We have shown that correct inclusion of the solvent effects usually changes severally the results of computer simulations on chemistry in aqueous solution. Not only the chemical reaction energetics is almost always modified dramatically, but also the extra degrees of freedom provided by the solvent molecules can open new mechanistic routes, which are not found for the reaction in vacuo. These effects should not be neglected, when studying chemical reactions that typically take place in water. We found that the Car-Parrinello molecular dynamics simulation method is a most powerful technique to compute physical properties of chemicals in solution and to obtain microscopic insight in solvent structure and solvation dynamics. The most important limitations that we encountered involve the quality of the present-day DFT exchange-correlation functionals and the rather high computational demand of the method in combination with our systems. Nevertheless, first principles computer simulation of rare events in aqueous solution is no longer out of reach, not even the simulation of such complicated systems as transition metal catalyzed oxidation reactions in water.