Introduction

Although the strength of molecular dynamics (MD) simulations lies in its use as a statistical tool to compute physical properties as ensemble averages, it is also used to simulate and analyze only single (or very few) illustrative events, for instance to show the possibility of certain chemical reaction pathways (see for some diverse examples refs. MoMaPaZi,AaMeBu98,MaHa01,YoShKaYa00,speybroeck1). Chemical reactions are commonly activated processes, i.e. the reactants have to pass a free energy barrier (transition state) before transforming into the products. If this reaction barrier is high compared to the thermal energy of the system, the probability becomes very small to find the reactants in the transition state and thus observe a reaction event on the timescale of the thermal molecular vibrations. In molecular dynamics simulations, which is a good technique to sample the thermal molecular vibrations, the chemical reaction is then a rare event. To simulate a reactive pathway, one therefore has to manipulate the system and enforce a reactive encounter by using for instance a geometric constraint or umbrella potential or by introducing kinetic energy in some translational, rotational or vibrational mode. Obviously, this manipulation makes the dynamics of the illustrative reaction pathway less realistic. Whether the found pathway is indeed a representative one can be tested using the technique to generate new pathways from earlier found pathways known as the transition path sampling method developed by Bolhuis, Dellago and Chandler[212,213]. This amounts to making small random changes to the atomic momenta of a randomly chosen configuration along an existing pathway and integrating the equations of motion backwards and forwards in time from this new point in phase space. If this leads again from reactants to products, it is accepted as a new pathway, which can then again be used to generate new pathways, and so forth.

We have recently studied the dissociation of hydrogen peroxide by iron(II) ions in aqueous solution, using ab initio (DFT) molecular dynamics simulations (AIMD)[171,146]. The mixture of Fe$^{2+}$ and H$_2$O$_2$, also known as Fenton's reagent[155], is known to oxidize organic substrates either via OH. radicals[156] or via a high-valent iron oxo species[157] (ferryl ion) as the reactive intermediate. In this previous work, we have generated reactive pathways of the reaction between iron(II) and hydrogen peroxide, starting from two different initial conditions: a) starting from hydrogen peroxide coordinated to pentaaqua iron(II) in water, we found a pathway which resulted in the formation of a dihydroxo iron(IV) moiety, which transformed into the ferryl ion in a second step a short time (3 ps) later; b) starting from hydrogen peroxide and pentaaqua iron(II) separated from each other an arbitrary distance ($\approx$ 4 Å), we found a pathway which lead to the ferryl ion via a more direct rebound mechanism. In both cases, we used bond distance constraints in the initial MD steps to prepare a configuration that would lead to a reactive pathway. Unfortunately, the resulting pathways are likely to contain a memory of the preparation using the unphysical constraints.

In this work, we will apply the transition path sampling technique to generate new reaction paths from the earlier found pathway of the hydrogen peroxide coordinating to and reacting with pentaaqua iron(II) in water. In principle, transition path sampling can be (and, in classical MD, has been) used to generate thousands of reaction paths, in which case again the statistical strength of MD is put to good use and quantities such as the reaction rate and rates of energy dissipation can be obtained. Unfortunately at the present time, the computational cost of ab initio molecular dynamics, which is needed to accurately describe the bond-breaking and making and the changing oxidation state of iron during the course of our reaction in water, does not allow for the simulation of more than a few tens of reaction pathways for systems as large as ours (see for another example ref. GeDeChHu). This should however be enough to obtain reaction paths which have lost the memory of our initial artificial system preparation and show that either or both of the earlier found pathways are indeed realistic illustrative reaction paths or perhaps even reveal a new, third, reaction mechanism.

This article is structured as follows. We will first give a short recapitulation of our previous results on the Fenton reaction and our starting reaction pathway. Then, we briefly summarize the computational details of our Car-Parrinello MD simulations in section 7.3, followed by an introduction of the transition path sampling technique in section 7.4. The results are presented in section 7.5, starting with the estimation of the transition state of the initial path, following a new procedure. After that, the results of the pathway relaxation are presented in section 7.5-7.5.2. These results, interestingly showing a change in the reaction mechanism, and in fact confirming previous predictions, are rationalized in the discussion part (section 7.6), which is followed by the conclusions.