Methodology

All electronic structure calculations were performed using the density functional theory (DFT) method (see e.g. ref. DFT). We used the Becke-88 gradient corrected exchange functional[67] and the Perdew-86 gradient corrected correlation functional.[38] The ab initio (DFT) molecular dynamics calculations of the systems including the solvent environment were done with the Car-Parrinello (CP) method[49] as implemented in the CP-PAW code developed by Blöchl[24]. The one-electron valence wave functions were expanded in an augmented plane wave basis up to a kinetic energy cutoff of 30 Ry. The frozen core approximation was applied for the 1s electrons of O, and up to 3p for Fe. 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 characteristic feature of the Car-Parrinello approach is that the electronic wave function, i.e. the coefficients of the plane wave basis set expansion, are dynamically optimized to be consistent with the changing positions of the atomic nuclei. The mass for the wave function coefficient dynamics was $\mu_e$=1000 a.u., which limits the MD timestep to $\delta t=0.19$ fs. To maintain a constant temperature of $T=300$K, a Nosé thermostat [69] was applied with a period of 100 femto-second. Periodic boundary conditions were applied to the cubic systems containing one iron ion, one hydrogen peroxide molecule and 31 water molecules. The size of the cubic box was 9.900 Å. The positive charge of the systems was compensated by a uniformly distributed counter charge.