Computational details

The microscopic reaction pathways presented in this paper were computed using the Car-Parrinello method.[49] This method combines classical molecular dynamics (MD) with a quantum mechanical computation of the electronic structure. The forces on the nuclei are obtained from the electronic ground-state energy, rather than from an empirical forcefield, which is why Car-Parrinello MD often is referred to as ab initio molecular dynamics (AIMD). The Car-Parrinello technique differs from other AIMD methods by the dynamical optimization method, known as simulated annealing, for the electronic wavefunction degrees of freedom (i.e. the basis set expansion coefficients) which can be treated simultaneous with, and on the same footing as, the Newtonian nuclear dynamics. The parameter of inertia, usually called ``fictitious electronic mass'', controlling the response of the basis set coefficients to the potential in the simulated annealing approach, is chosen much smaller than the atomic masses. This way, the wavefunction adapts instantaneously to the moving nuclei, keeping the electrons sufficiently close to the correct ground-state. The Verlet algorithm is used to integrate the equations of motion for the nuclear dynamics and for the coefficient optimization. In the present work, the fictitious electron mass was chosen 500 a.u. ( $\approx 4.555\cdot 10^{-28}$ kg), which limits the time step to 0.145 fs.

The electronic structure is calculated with density functional theory (DFT).[46] We used the Becke-88[36] and Perdew-86 [38] gradient corrected functionals for exchange and electron correlation respectively. The Car-Parrinello simulations were performed using the CP-PAW program developed by Blöchl,[24] who integrated the projector augmented wave (PAW) method with the ab initio molecular dynamics. The PAW method uses an augmented plane-wave basis for the electronic valence wavefunctions and, in the current implementation, frozen atomic wavefunctions for the core states. The plane-wave basis expansion was cut off after functions with a kinetic energy of 30 Ry. The 1s electrons of oxygen and up to the 3p electrons for iron were kept frozen. For the augmentation for H and O, one projector function per angular momentum quantum number was used for $s$- and $p$-angular momenta. For Fe, one projector function was used for $s$- and $p$- and two for $d$-angular momenta.

The advantage of PAW over the more commonly used pseudopotential approach is that transferability problems should be largely avoided. However, as with all methods in which the core is represented approximately, there will be some loss of accuracy. Extensive tests are therefore required and we have previously shown that bond energies and geometries of small molecules and complexes computed with CP-PAW agree very well with highly accurate all-electron DFT results obtained with the ADF[70] program[144]. However, we have in the course of the present work noted that reaction energies of chemical reactions involving a change in the formal oxidation state of iron sometimes exhibit relatively large (several kcal/mol) discrepancies with accurate (large basis set, all-electron) ADF calculations. Table 7.1 shows four reaction energies of gas phase reactions involving simple aqua iron complexes, calculated with CP-PAW and with ADF. The same exchange-correlation functional was used in both calculations. The Kohn-Sham orbitals were expanded in the ADF calculation a large even-tempered all-electron Slater-type basis set containing: 4 s, 2 p, and 1 d functions for hydrogen; 6 s, 4 p, 2 d, and 1 f functions for oxygen; and 11 s, 7 p, 5 d, and 1 f functions for iron[175]. The CP-PAW calculation used plane-waves with a cut-off of 30 Rydberg. Reaction A and B, show the typical small differences in the order of 1 kcal/mol due to the differences in basis set and the frozen core approximation used with PAW. However, reactions C and D, which involve a change in the iron oxidation state, show larger discrepancies, up to 6.3 kcal/mol for the reaction in which iron(II) is oxidized to iron(III). The error does not seem to be due to the plane-wave cutoff of 30 Ry (it is not reduced when going to 50 Ry) and can be attributed to the partial waves for the inner region of the valence electrons and the projector functions used in the PAW calculations. The purpose of the present work is the application of transition path sampling in order to obtain reaction pathways which have lost the memory of their initial artificial construction. This memory effect, as we will see, is related to the solvent motion and the formation and breaking of H-bond networks in the solvent, which are well represented by the present type of CPMD simulations, as demonstrated earlier [144]. We therefore feel the present CPMD method can be applied to this type of study, but we have to take the limitations due to the limited accuracy of the reaction energetics in mind; obviously one cannot realistically calculate a reaction rate with this type of CP-PAW for a reaction involving change of oxidation state of Fe.

Table 7.1: Bond dissociation energies and rearrangement energies in kcal/mol for isolated (``gas phase'') complexes, calculated with the Amsterdam Density Functional program (ADF) and the CP-PAW program using the Becke-Perdew exchange correlation functional.

	Gas phase reaction			CP-PAW	ADF	$\Delta \Delta$E
A	Fe $^{\mathrm{II}}$(H$_2$O)$_6$	$\rightarrow$	Fe $^{\mathrm{II}}$(H$_2$O)$_5$ + H$_2$O	21.7	22.1	-0.4
B	Fe $^{\mathrm{II}}$(H$_2$O)$_6$	$\rightarrow$	Fe $^{\mathrm{II}}$(H$_2$O)$_5$-H$_2$O	-2.3	-3.5	1.2
C	Fe$^{\rm {II}}$(H$_2$O)$_5$ + H$_2$O$_2$	$\rightarrow$	Fe$^{\rm {III}}$(H$_2$O)$_5$OH + OH.	-8.4	-2.1	-6.3
D	Fe$^{\rm {III}}$(H$_2$O)$_5$OH + OH.	$\rightarrow$	Fe$^{\rm {IV}}$(H$_2$O)$_5$O + H$_2$O	-30.8	-28.7	-2.1

For the AIMD simulations of the Fenton reaction in aqueous solution in this paper, we applied periodic boundary conditions to a cubic box with an edge of 9.900 Å, containing one iron ion, one hydrogen peroxide molecule and 31 water molecules. The positive 2+ charge of this system was compensated by a uniformly distributed counter charge. In the reaction path sampling technique, we quenched the fictitious basis set coefficient dynamics after every fourth reaction pathway generation, to avoid deviations from the Born-Oppenheimer surface. A Nosé thermostat[69] maintained a constant temperature of $T=300$ K.