Method

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] For the calculations of energies, frequencies and optimizations of geometries of the molecular structures in vacuo, we used the Slater type orbital based ADF package.[70] The Kohn-Sham orbitals were expanded in 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 carbon and oxygen; and 11 s, 7 p, 5 d, and 1 f functions for iron[175]. Finite-temperature reaction enthalpies at $T=300$ K and the entropies were estimated using

$$\Delta H^\mathrm{300K} = \Delta E_0 + \Delta E_\mathrm{ZPE} + \Delta E^\mathrm{300K}_\mathrm{int}$$ (80)

$$\Delta E^\mathrm{300K}_{{\rm int}} = \Delta E^\mathrm{v}_T + \Delta E^\mathrm{t} + \Delta E^\mathrm{r}$$ (81)

$$\Delta S = R \ln(Q^{\mathrm{t}}Q^{\mathrm{r}}Q^{\mathrm{v}})$$ (82)

with $E_0$ the sum of the electronic energy in a static nuclear field (Born-Oppenheimer approximation) and the nuclear electrostatic repulsion. The zero-point vibrational energy $E_\mathrm{ZPE}$ and the temperature dependent vibrational energy $E^\mathrm{v}_T$ were calculated from the unscaled DFT-BP frequencies, within the harmonic approximation. The change in translational energy $\Delta E^\mathrm{t}$ and rotational energy $\Delta E^\mathrm{r}$ were obtained using the ideal gas law, associating $\frac{1}{2}k_BT$ to each degree of freedom. The partition function $Q$ is the product of translational, rotational and vibrational contributions (see e.g. chapter 20 in reference [78]).

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 time step to $\delta t=0.19$ fs.

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 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 8.1 shows four reaction energies of gas phase reactions involving simple aqua iron complexes, calculated with CP-PAW and with ADF.

Table 8.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$]$^{2+}$ $\rightarrow$ [Fe $^{\mathrm{II}}$(H$_2$O)$_5$]$^{2+}$ + H$_2$O	21.7	22.1	-0.4
B	$[$Fe $^{\mathrm{II}}$(H$_2$O)$_6$]$^{2+}$ $\rightarrow$ [Fe $^{\mathrm{II}}$(H$_2$O)$_5$]$^{2+}$-H$_2$O	-2.3	-3.5	1.2
C	$[$Fe$^{\rm {II}}$(H$_2$O)$_5]^{2+}$ + H$_2$O$_2$ $\rightarrow$
	$[$Fe$^{\rm {III}}$(H$_2$O)$_5$OH$]^{2+}$ + OH.	-8.4	-2.1	-6.3
D	$[$Fe$^{\rm {III}}$(H$_2$O)$_5$OH$]^{2+}$ + OH. $\rightarrow$
	$[$Fe$^{\rm {IV}}$(H$_2$O)$_5$O$]^{2+}$ + H$_2$O	-30.8	-28.7	-2.1

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. Also in this work, we will simulate a reaction involving a changing iron oxidation state (from iron(IV) to iron(III)) in aqueous solution. We will therefore only discuss the solvent effects on this reaction, by comparing the results with the same reaction in vacuo, rather than attempting to compute the reaction rate for this reaction in aqueous solution.