The iron(III) complexes [Fe(H2O)n(OH)m]3−m (n + m = 5, 6, m ≤ 3) and corresponding proton transfer reactions are studied with total energy calculations, the nudged elastic band (NEB) method, and molecular dynamics (MD) simulations using ab initio and a modification of reactive force field potentials, the ReaxFF-AQ potentials, based on the implementation according to Böhm et al. [J. Phys. Chem. C 120, 10849–10856 (2016)]. Applying ab initio potentials, the energies for the reactions [Fe(H2O)n(OH)m]3−m + H2O → [Fe(H2O)n−1(OH)m+1]2−m + H3O+ in a gaseous environment are in good agreement with comparable theoretical results. In an aqueous (aq) or alkaline environment, with the aid of NEB computations, respective minimum energy paths with energy barriers of up to 14.6 kcal/mol and a collective transfer of protons are modeled. Within MD simulations at room temperature, a permanent transfer of protons around the iron(III) ion is observed. The information gained concerning the geometrical and energetic properties of water and the [Fe(H2O)n(OH)m]3−m complexes from the ab initio computations has been used as reference data to optimize parameters for the O–H–Fe interaction within the ReaxFF-AQ approach. For the optimized ReaxFF-AQ parameter set, the statistical properties of the basic water model, such as the radial distribution functions and the proton hopping functions, are evaluated. For the [Fe(H2O)n(OH)m]3−m complexes, it was found that while geometrical and energetic properties are in good agreement with the ab initio data for gaseous environment, the statistical properties as obtained from the MD simulations are only partly in accordance with the ab initio results for the iron(III) complexes in aqueous or alkaline environments.

The [Fe(H2O)n(OH)m]3−m (n + m = 5, 6, m ≤ 3) complexes in an aqueous environment are of importance in industrial processes, such as in the electrochemical machining (ECM) process. Considering the modeling of an ECM process on the atomistic scale,2 the aquoferric layer is one essential subsystem of the full iron-oxide/electrolyte system.3 Experimentals4–9 and theoretical10–16 investigations have been performed for iron(III) in an aqueous environment. Fe–O distances between 1.99 and 2.05 Å could be determined from measurements (cf. Table 2 of Ref. 5). In the theoretical studies by Martin et al.,10 De Abreu et al.11, and Ottonello and Zuccolini12 performed at the level of (hybrid) density functional theory (DFT), typically at the first step, the geometrical and energetic properties of [Fe(H2O)n(OH)m]3−m complexes in a gaseous environment are evaluated and afterward corrected for the solute–solvent interaction with continuum models.17,18 In this way, accurate deprotonation energies and, partly, changes of geometries have been predicted. Furthermore, with molecular dynamics simulations13–16 for Fe3+ in a pure aqueous environment, important statistical properties, such as radial distribution functions, and the coordination number distribution have been investigated.

Taking into consideration that with material removal rates of mm/min ≃ 10−7 Å/ps, the ECM process is a slow process at the atomistic level and thus, to observe relevant events in a molecular dynamics (MD) simulation, run time and system size have to be in range of ns and ∼1M atoms. Rather than ab initio based approaches, the reactive force field20–22 (ReaxFF) method might be a more appropriate method for such a scale.19 The ReaxFF method has been used to describe solids, liquids, and interfaces with an adequate compromise between computation speed and accuracy in many cases.23 

In the present study, we focused on the modified reactive force field version implemented as reported by Böhm et al.1 [denoted as ReaxFF-AQ (RxAQ)]. Within the ReaxFF-AQ approach, the most energy contributions from the original ReaxFF method are evaluated, but, in particular, modifications for the Coulomb energy and the charge energy are applied. Additionally, within the reactive force field implementation, according to Ref. 1, a parameter optimization procedure was included, which has been partly used as a basis for the present study.

In the described context, the aim of the present study is the development of a parameter set for the RxAQ method to be able to describe the properties of iron(III) in aqueous and alkaline environments with RxAQ in accordance with the predictions at the level of DFT. To this end, first step ab initio methods are applied to investigate the system properties and provide adequate reference data [or sampling points of the corresponding ab initio potential energy surfaces (AI-PESs)] to adjust the RxAQ parameters for the O–H–Fe interaction. After obtaining an optimized parameter set for the RxAQ approach, the predictions at the level of empirical potentials are compared with the ab initio results.

This article is organized as follows: First, the methods for the ab initio calculations are presented. RxAQ parameter optimization and calculations are described in Sec. II. The results are shown in Sec. III. Finally, the study is summarized in Sec. IV.

The gaseous environment of the [Fe(H2O)n(OH)m]3−m complexes is realized by the application of an otherwise empty 30 × 30 × 30 Å3 supercell. In an aqueous environment, the [Fe(H2O)n(OH)m]3−m (n + m = 5, 6, m ≤ 3) complexes are modelled by a periodic supercell approach containing one water molecule per 3.1 × 3.1 × 3.1 Å3 (corresponding to the equilibrium water density of 1 g/cm3) and with cubic supercells containing 27 or 64 water molecules (i.e., with lattice constants of 9.1 or 12.4 Å). With the replacement of up to three water molecules by a hydroxyl group OH around the iron ion, an alkaline environment with pOH ≤ 0.05 is configured.

The calculations are performed at the level of spin-polarized DFT with the GGA-PBE functional26 as implemented in SIESTA.27,28 The electron–electron interaction has been computed by applying localized basis sets at the level of double-ζ plus polarization with a real space sampling cutoff energy of 150 Ry. Norm-conserving non-local pseudopotentials of the Troullier–Martin29 type have been applied for all elements. To ensure charge neutrality for the applied supercells, a compensating background charge has been applied in the calculations.

To investigate the minimum energy path (MEP) of the proton transfer reactions
[Fe(H2O)n(OH)m]3m+H2O(aq)[Fe(H2O)n1(OH)m+1]2m+H3O+(aq),
(1)
the nudged elastic band method, as implemented in the Atomic Simulation Environment (ASE),30 is used. In detail, the forces used to apply the improved tangent method are according to Ref. 31. After convergence of the minimum energy path, the climbing image method is used to compute an accurate value of the energy barrier (cf. Ref. 32). A convergence threshold of 0.05 eV/Å for the force has been applied for the nudged elastic band (NEB) calculations.

Molecular dynamics simulation is performed at room temperature by applying the NVT Nosé thermostat33 as implemented in SIESTA. A nose mass of 100 Ryfs2 and a time step of 1.0 fs are used in the simulations. The radial distribution function gAB (RDF) between atom species A and B and the distance dependent coordination numbers (as the radial integral of the RDF, e.g., Ref. 34) are computed for every time step (after an equilibration time) and averaged over the time. Furthermore, during the run time of the NVT MD simulations, the first hydration shell around the Fe3+ ion (Fe–O cutoff distance 2.3 Å) is checked for an appearance of H2O molecules and OH groups. A cutoff distance of dOH*=1.2 Å in the HO–H* complex is used to distinguish between a water molecule and a hydroxyl group with an additional H* in the vicinity.

Furthermore, the calculations not only for the iron(III) structures but also for general water structures (i.e., the sampling points of the considered AI-PESs) are used as reference data for the RxAQ parameter optimization (as presented below).

In the present work, the standalone implementation of the reactive force field approach from Ref. 1 was initially used and has been modified. The potential energy within the RxAQ method (consisting of one-, two-, and three-body energy contributions) is given as follows:
ERxAQ=Ebond+Elp+Eover+Eunder+Eval+EH-bond+EvdW+ECoul+Ech+Ecorr,
(2)
where the bond energy (Ebond), the lonepair energy (Elp), the over co-ordination and under co-ordination energies (Eover and Eunder), the valency energy (Eval), and the hydrogen bond energy (EH-bond) depend on the bond order between two atoms, as in the original ReaxFF method. Here, EvdW is the van der Waals energy, and Ecorr has been implemented within this work and is described in  Appendix A. Additional parameters to extend the description of the charge energy (Ech) and the Coulomb energy (ECoul) are described in Ref. 1. As a result, the Coulomb energy and the charge energy also exhibit a dependence on the bond order within the RxAQ formalism.

The Coulomb energy and the charge energy depend generally on the atomic charges. The charges on the atoms are computed within a minimization procedure for the sum of the charge energy and the Coulomb energy. A charge equilibration procedure in the spirit to the extended charge equilibration procedure35 (EQeq) has been implemented, but the constraint for the total charge Q = ∑iqi is incorporated by substitution into the quadratic minimization problem, and a compensating contribution is added to the atomic charges to ensure charge neutrality (see  Appendix B for details).

The RxAQ empirical parameters are optimized by minimizing the deviation ϵ, depending on the difference between the empirical and DFT reaction energies (EIRxAQ and EIDFT) as well as charges [q(rj)RxAQ and q(rj)DFT] for the chosen geometries I = (ijk…) and j,
ϵIωIL(EIRxAQEIDFT)+ωqjrjL(q(rj)RxAQq(rj)DFT),
(3)
where ωI and ωq are the weights and L is a loss function. The DFT reference geometries related to the properties of the water generated within this work are based on Refs. 22, 34, and 3638. The reference data for the [Fe(H2O)n(OH)m]3−m complexes were generated from basic calculations, NEB calculations, and NVT MD simulations. All together about 4000 structures, generated with DFT, were used as reference geometries. The reference charges of the atoms for selected geometries have been computed according to the Bader method,39–42 due to the stability of this method in combination with the calculations based on localized basis sets.43 

The RxAQ potential energy can depend on up to 200 empirical parameters. To obtain an optimal parameter set (i.e., the parameter set with the smallest deviation to the reference data), initial parameter sets were taken from the ReaxFF-related works of Zhang and van Duin34 and Aryanpour et al.44 The deviation ϵ was minimized by applying a combination of the genetic algorithm,45 a conjugate gradient, and a generation of additional auxiliary DFT structures.2 

To perform MD simulations within the RxAQ framework, the RxAQ package has been interconnected with the LAMMPS code.46 Transferring the calculated positions of the atoms and the forces on the atoms between the two packages opens up the possibility of performing MD simulations in NVE and NVT ensembles. For NVT MD simulations, the Nosé–Hoover thermostat47 with a temperature damping parameter (Tdamp) of 100 time steps was applied. A time step of 0.25 fs has been used in the MD simulations with RxAQ potentials.

As for ab initio potentials, NEB simulations have been performed applying the ASE. For RxAQ potentials, however, the convergence threshold for the forces was set at 0.13 eV/Å, and the climbing image procedure was not applied since convergence of the MEP was hardly achieved. The computed reaction paths are to be considered approximate.

Before the characteristics of iron(III) in aqueous and alkaline environments are considered with RxAQ potentials, the description of the properties of water with hydronium and hydroxyl ions at room temperature is examined. Apart from verifying the radial distribution function for water, the respective Grotthuss hopping functions of hydronium and hydroxyl ions are computed. The Grotthuss hopping mechanism mainly contributes to the mobility of protons compared to oscillatory shuttling. The Grotthuss hopping function h(t) can be defined by34,48
h(t+Δt)=h(t)+Δh(Δt),
(4)
with h(0) = 0 and the time step being Δt. In addition, Δht) is 0, 1, and −1 for no proton transfer between water molecules, for proton hopping to a new water molecule, and proton hopping to the previous water molecule within the time Δt. To detect H3O+ and OH complexes, the oxygen atom of the corresponding complex as the reference point is identified. For H3O+, the proton H+ is defined as the hydrogen atom with the largest distance to the central oxygen atom. Non-vanishing contributions to the hopping function are related to changes of the reference oxygen atom.

The DFT energies of the considered [Fe(H2O)n(OH)m]3−m complexes (n + m = 5, 6 and m ⇐ 3) and related reaction energies E = EFEI (where EI and EF are the energies of the initial and final geometries, respectively) in a gaseous environment are summarized in Table I. The different geometries for the [Fe(H2O)4(OH)2]+ complex (cis, trans, and outer) are denoted in accordance with Ref. 10. The proton transfer reactions [Fe(H2O)n(OH)m]3−m + H2O → [Fe(H2O)n−1(OH)m+1]2−m + H3O+ are exothermic in all the considered cases with values between −193.2 and −37.9 kcal/mol. The computed energies are in agreement with the results from De Abreu et al.11 by applying DFT-PBE-LCAO potentials. Max. deviations of 40% from the computations of Martin et al.10 (using DFT-B3LYP potentials49) are observed. In particular, the absolute values for the energies of the reactions [Fe(H2O)5OH]2+ + H2O → [Fe(H2O)4(OH)2]+ + H3O+ are overestimated as compared to Ref. 10 in the present calculations. The deviations are explained with the dependence of the values on the description of the exchange and correlation (XC) energy (PBE in the present work) as well as on the basis set. The scattering of the reaction energies with the basis set and the XC functional is also presented in Table 2 of Ref. 11. In particular, for the configuration PBE/DZVP, a value of −37.6 kcal/mol can be computed from the work of De Abreu et al. for the reaction [Fe(H2O)5OH]2+ + H2O → [Fe(H2O)4(OH)2]+ + H3O+ to be closer to the present results.

TABLE I.

Energies E (in kcal/mol) of proton and H2O transfer reactions in a gaseous environment calculated with ab initio (AI) and RxAQ (RxAQ) potentials compared with the theoretical results (Theo.).

ReactionsE(AI)E(RxAQ)Theo.
(Fe(H2O)6)3+ + H2O → (Fe(H2O)5(OH))2+ + H3O+ −147.5 −196.5 −143.7a, −150.0b 
(Fe(H2O)5OH)2+ + H2O → cis-(Fe(H2O)4(OH)2)+ + H3O+ −45.1 −26.8 −27.5a, ∼−37.6b 
(Fe(H2O)5OH)2+ + H2O → trans-(Fe(H2O)4(OH)2)+ + H3O+ −37.6 −31.2 −25.5a 
(Fe(H2O)5OH)2+ + H2O → outer-(Fe(H2O)4(OH)2)+ + H3O+ −43.2 −26.1 −26.6a 
(Fe(H2O)6)3+ + 2H2O → cis-(Fe(H2O)4(OH))+ + 2H3O+ −193.6 −222.2 −187.6b 
Outer-(Fe(H2O)4(OH)2)+ → (Fe(H2O)3(OH)2)+ + H218.4 11.5 9.6a 
(Fe(H2O)6)3+ + 2H2O → (Fe(H2O)2(OH)3) + 3H3O+ −108.1 −75.1 −112.4b 
ReactionsE(AI)E(RxAQ)Theo.
(Fe(H2O)6)3+ + H2O → (Fe(H2O)5(OH))2+ + H3O+ −147.5 −196.5 −143.7a, −150.0b 
(Fe(H2O)5OH)2+ + H2O → cis-(Fe(H2O)4(OH)2)+ + H3O+ −45.1 −26.8 −27.5a, ∼−37.6b 
(Fe(H2O)5OH)2+ + H2O → trans-(Fe(H2O)4(OH)2)+ + H3O+ −37.6 −31.2 −25.5a 
(Fe(H2O)5OH)2+ + H2O → outer-(Fe(H2O)4(OH)2)+ + H3O+ −43.2 −26.1 −26.6a 
(Fe(H2O)6)3+ + 2H2O → cis-(Fe(H2O)4(OH))+ + 2H3O+ −193.6 −222.2 −187.6b 
Outer-(Fe(H2O)4(OH)2)+ → (Fe(H2O)3(OH)2)+ + H218.4 11.5 9.6a 
(Fe(H2O)6)3+ + 2H2O → (Fe(H2O)2(OH)3) + 3H3O+ −108.1 −75.1 −112.4b 
a

Martin et al.,10 Table 3.

b

De Abreu et al.,11 Table 2 (PBE/DZVP values).

Geometrical and electronic properties of the [Fe(H2O)n(OH)m]3−m complexes in a gaseous environment are summarized in Table II. Generally, the Fe–OH2O distances (for oxygen bound in water) are larger than the Fe–OOH distances (for oxygen bound in the OH-group). Furthermore, increasing Fe–O distances with decreasing total charge of the [Fe(H2O)n(OH)m]3−m complex are observed. For the [Fe(H2O)6]3+ complex, the computed average Fe–OH2O distance of 2.09 Å shows a slight deviation from the B3LYP results of Moin et al.15 with 2.05 Å. The tendency of increasing Fe–O distances with the number of present hydroxyl groups is compatible with the results of De Abreu et al.11 and Martin et al.10, both reporting minimal and maximal Fe-O distances. Considering the Bader charges on the atoms, it can be observed that with increasing number of hydroxyl groups (for m = 0, 1, 2), a charge accumulation on oxygen and iron atoms occurs. For m = 2 → 3, the trend is reversed for iron but still holds for oxygen. The charge on hydrogen atoms is continuously decreasing with an increasing number of hydroxyl groups.

TABLE II.

Average Fe–O, Fe–H distances (d̄ in Å), and charges (q̄ in |e|) on atoms as computed with the Bader approach based on the ab initio electronic density (AI) and within the RxAQ approach (RxAQ) for [Fe(H2O)nOHm]3−m (m = 0…3) in a gaseous environment. OH2O and HH2O (OOH and HOH) denote oxygen and hydrogen atoms, respectively, bound in a water molecule H2O (in a hydroxyl group OH).

Property (AI/RxAQ) (Fe(H2O)6)3+ (Fe(H2O)5OH)2+ cis-(Fe(H2O)4OH2)1+ (Fe(H2O)2OH3
d̄(Fe–OH2O2.09/2.08 2.17/2.12 2.23/2.21 2.29/2.28 
d̄(Fe–OOH⋯ 1.78/1.79 1.84/1.83 1.90/1.90 
q̄(Fe) 1.34/1.69 1.54/1.75 1.66/1.68 1.63/1.50 
q̄(OH2O−0.61/−0.71 −0.70/−0.75 −0.72/−0.78 −0.76/−0.81 
q̄(OOH⋯ −0.79/−0.81 −0.87/−0.85 −0.89/−0.91 
q̄(HH2O0.44/0.46 0.42/0.43 0.40/0.40 0.40/0.40 
q̄(HOH⋯ 0.52/0.52 0.38/0.45 0.32/0.42 
Property (AI/RxAQ) (Fe(H2O)6)3+ (Fe(H2O)5OH)2+ cis-(Fe(H2O)4OH2)1+ (Fe(H2O)2OH3
d̄(Fe–OH2O2.09/2.08 2.17/2.12 2.23/2.21 2.29/2.28 
d̄(Fe–OOH⋯ 1.78/1.79 1.84/1.83 1.90/1.90 
q̄(Fe) 1.34/1.69 1.54/1.75 1.66/1.68 1.63/1.50 
q̄(OH2O−0.61/−0.71 −0.70/−0.75 −0.72/−0.78 −0.76/−0.81 
q̄(OOH⋯ −0.79/−0.81 −0.87/−0.85 −0.89/−0.91 
q̄(HH2O0.44/0.46 0.42/0.43 0.40/0.40 0.40/0.40 
q̄(HOH⋯ 0.52/0.52 0.38/0.45 0.32/0.42 

In aqueous and alkaline environments, the energies of the proton transfer reactions are computed from randomly chosen local geometrical minima in the corresponding supercell. In the present investigations, local energetic minima for geometries [Fe(H2O)n(OH)m]3−m + H3O+ (aq) have been found only by considering a separation layer of water between the [Fe(H2O)n(OH)m]3−m complex and the H3O+ ion. The energies on the MEP (computed with NEB) for the proton transfer reactions are shown in Fig. 1. The corresponding reaction energies and activation energies are summarized in Table III. The reaction energies are between 3.7 and 10.8 kcal/mol and clearly different in sign and absolute values to the results for the gaseous environment. The computed reaction energies are in the same order of magnitude as those of the results from Martin et al.,10 De Abreu et al.11 (both applying a continuum model to account for solute–solvent interactions) and experimental values from Flynn24 and Meagher.25 For the activation energies of the reactions EA = Emax≠FEI (where Emax≠F is the energy of the transition state) up to 14.6 kcal/mol have been determined. No energy barrier has been found for the proton transfer reaction [Fe(H2O)6]3+ + H2O (aq) → [Fe(H2O)5(OH)]2+ + H3O+ (aq).

FIG. 1.

Energies on the minimum energy path of the proton transfer reactions [Fe(H2O)n(OH)m]3−m + H2O ↔ [Fe(H2O)n−1(OH)m+1]2−m + H3O+ (m = 0…2) in an aqueous/alkaline environment as computed with the NEB method using ab initio (AI) and RxAQ (RxAQ) potentials (also presented in Table III). The notation of the reactions R0…R2 is in accordance with Table III.

FIG. 1.

Energies on the minimum energy path of the proton transfer reactions [Fe(H2O)n(OH)m]3−m + H2O ↔ [Fe(H2O)n−1(OH)m+1]2−m + H3O+ (m = 0…2) in an aqueous/alkaline environment as computed with the NEB method using ab initio (AI) and RxAQ (RxAQ) potentials (also presented in Table III). The notation of the reactions R0…R2 is in accordance with Table III.

Close modal
TABLE III.

Energies E and activation energies EA (in kcal/mol) of proton transfer reactions [Fe(H2O)n(OH)m]3−m ↔ [Fe(H2O)n−1(OH)m+1]2−m + H+ (m = 0…2, labeled with Rm) in an aqueous or alkaline environment as calculated by applying the NEB method using ab initio (AI) and RxAQ (RxAQ) potentials compared with the further theoretical results (Theo.).

ReactionsLabelE(AI/RxAQ)Theo.EA(AI/RxAQ)
(Fe(H2O)6)3+ + H2O (aq) → (Fe(H2O)5(OH))2+ + H3O+ (aq) R0 3.7/9.1 1.6a, −1.3b, 3c -/- 
(Fe(H2O)5OH)2+ + H2O (aq) → cis-(Fe(H2O)4(OH)2)+ + H3O+ (aq) R1-i 7.4/12.8 18.3a, 5d 8.3/- 
(Fe(H2O)5OH)2+ + H2O (aq) → trans-(Fe(H2O)4(OH)2)+ + H3O+ (aq) R1-ii, 9.5/12.1, 16.1a 12.2/16.7 
R1-iii 10.8/3.6 14.6/9.2 
(Fe(H2O)3(OH)2)1+ + H2O (aq) → (Fe(H2O)2(OH)3) + H3O+ (aq) R2 3.8/7.7  6.8/11.6 
ReactionsLabelE(AI/RxAQ)Theo.EA(AI/RxAQ)
(Fe(H2O)6)3+ + H2O (aq) → (Fe(H2O)5(OH))2+ + H3O+ (aq) R0 3.7/9.1 1.6a, −1.3b, 3c -/- 
(Fe(H2O)5OH)2+ + H2O (aq) → cis-(Fe(H2O)4(OH)2)+ + H3O+ (aq) R1-i 7.4/12.8 18.3a, 5d 8.3/- 
(Fe(H2O)5OH)2+ + H2O (aq) → trans-(Fe(H2O)4(OH)2)+ + H3O+ (aq) R1-ii, 9.5/12.1, 16.1a 12.2/16.7 
R1-iii 10.8/3.6 14.6/9.2 
(Fe(H2O)3(OH)2)1+ + H2O (aq) → (Fe(H2O)2(OH)3) + H3O+ (aq) R2 3.8/7.7  6.8/11.6 
a

Martin et al.,10 Table 4.

b

De Abreu et al.,11 Table 2 (PBE/DZVP values).

c

Flynn.24 

d

Meagher.25 

Since a separation layer of water is present between the (Fe(H2O)n(OH)m)3−m complex and the H3O+ ion in an aqueous/alkaline environment, a collective proton movement with at least two involved protons occur within the proton transfer reactions. The initial, intermediate, and final geometries for the reactions with m = 0 and 1 are shown in Figs. 2 and 3.

FIG. 2.

Initial, intermediate, and final geometries for the reaction (Fe(H2O)6)3+ + H2O (aq) → (Fe(H2O)5(OH))2+ + H3O+ (aq) (R0). The white, red, and yellow spheres denote hydrogen, oxygen, and iron atoms. Surrounding water molecules are illustrated as sticks. The transition of protons is indicated by the arrows.

FIG. 2.

Initial, intermediate, and final geometries for the reaction (Fe(H2O)6)3+ + H2O (aq) → (Fe(H2O)5(OH))2+ + H3O+ (aq) (R0). The white, red, and yellow spheres denote hydrogen, oxygen, and iron atoms. Surrounding water molecules are illustrated as sticks. The transition of protons is indicated by the arrows.

Close modal
FIG. 3.

Initial, intermediate, and final geometries for the reaction (Fe(H2O)5OH)2+ + H2O (aq) → cis-(Fe(H2O)4(OH)2)+ + H3O+ (aq) (R1-i). The colors and symbols are in accordance with Fig. 2.

FIG. 3.

Initial, intermediate, and final geometries for the reaction (Fe(H2O)5OH)2+ + H2O (aq) → cis-(Fe(H2O)4(OH)2)+ + H3O+ (aq) (R1-i). The colors and symbols are in accordance with Fig. 2.

Close modal

NVT MD simulations have been performed at room temperature for up to 25 ps for the [Fe(H2O)n(OH)m]3−m complexes in aqueous/alkaline environments [denoted by [Fe(OH)m]3−m(aq) in the following]. The statistical properties have been evaluated for the last ∼15 ps. To check for proton transfer, the number of OH groups around the central Fe3+ ion has been monitored with the NVT MD run time. The results are shown in Fig. 4. Indeed, proton transfer was found in all [Fe(OH)m]3−m(aq) systems. The proton hopping is most frequently occurring in the [Fe(OH)]2+(aq) system. The appearance of H3O+ ions or Zundel ions in the NVT MD trajectories is also demonstrated in the NVT MD snapshot geometries, which are shown in Fig. 5 (at times, shown in Fig. 4).

FIG. 4.

Monitoring of the number of OH groups around the Fe3+ ion with the NVT MD run time (for every 1 fs) using ab initio potentials for the systems [Fe(OH)m]3−m (aq) (m = 0…3).

FIG. 4.

Monitoring of the number of OH groups around the Fe3+ ion with the NVT MD run time (for every 1 fs) using ab initio potentials for the systems [Fe(OH)m]3−m (aq) (m = 0…3).

Close modal
FIG. 5.

Snapshots of the ab initio NVT MD trajectory of [Fe]3+(aq) for times T0 (a), T1 (b), and T2 (c) as indicated in Fig. 4 showing zero OH groups, one OH group (and a Zundel cation), and two OH groups (and a H3O+ group and a Zundel cation) around the [Fe]3+ ion. The colors and symbols are in accordance with Fig. 2.

FIG. 5.

Snapshots of the ab initio NVT MD trajectory of [Fe]3+(aq) for times T0 (a), T1 (b), and T2 (c) as indicated in Fig. 4 showing zero OH groups, one OH group (and a Zundel cation), and two OH groups (and a H3O+ group and a Zundel cation) around the [Fe]3+ ion. The colors and symbols are in accordance with Fig. 2.

Close modal

The respective radial distribution functions for Fe–O and (distance-dependent) coordination numbers for the [Fe(OH)m]3−m(aq) complexes are shown in Fig. 6. The main peak (at a binning of 0.05 Å) in the RDF is shifted toward lower values with a decreasing system charge. One reason is the more frequent existence of OH groups and its lower average Fe–OOH− distance compared to the Fe–OH2O value in aqueous/alkaline environments (see Table IV). The difference between the average Fe–OH2O and Fe–OOH distance is a maximum of 0.2 Å. With experiments, Fe–O distances in the range of 1.99 to 2.05 Å for Fe3+(aq) could be determined (cf. Table 2 of Ref. 5). Consequently, the Fe–O distances are slightly overestimated with the present calculations.

FIG. 6.

Fe–O (a) and Fe–H (b) radial distribution functions (RDF) and accumulated Fe–O coordination numbers (CN, c) as calculated from the NVT MD simulations with ab initio potentials for the systems [Fe(OH)m]3−m (aq) (m = 0…3).

FIG. 6.

Fe–O (a) and Fe–H (b) radial distribution functions (RDF) and accumulated Fe–O coordination numbers (CN, c) as calculated from the NVT MD simulations with ab initio potentials for the systems [Fe(OH)m]3−m (aq) (m = 0…3).

Close modal
TABLE IV.

Average Fe–O distances (d̄ in Å) for [FeOHm]3−m (m = 0…3) in aqueous and alkaline environments as computed from NVT MD simulations with ab initio (AI) and RxAQ (RxAQ) potentials. Notation of oxygen atoms is in accordance with Table II.

System[Fe]3+ (aq)[Fe(OH)]2+ (aq)[Fe(OH)2]1+ (aq)[Fe(OH)3] (aq)
d̄(Fe–OH2O) (AI/RxAQ) 2.09/2.13 2.06/2.15 2.10/2.16 2.00/2.17 
d̄(Fe–OOH) (AI/RxAQ) 1.99/1.85 1.91/1.86 1.91/1.86 1.90/1.87 
System[Fe]3+ (aq)[Fe(OH)]2+ (aq)[Fe(OH)2]1+ (aq)[Fe(OH)3] (aq)
d̄(Fe–OH2O) (AI/RxAQ) 2.09/2.13 2.06/2.15 2.10/2.16 2.00/2.17 
d̄(Fe–OOH) (AI/RxAQ) 1.99/1.85 1.91/1.86 1.91/1.86 1.90/1.87 

For the Fe3+(aq) system, the coordination number of the first hydration shell is 6, while for the systems [Fe(OH)m]3−m(aq) (m > 0), the coordination numbers for the first hydration shell are 5, 5, and 4 for m = 1, 2, and 3, respectively. The coordination number of 6 for the first hydration shell for the Fe3+(aq) system is in agreement with the MD simulations of Mandal et al.14 (with AI potentials), Rustad et al.16 (with Born–Meyer type potentials), and Moin et al.15 (applying the QMCF-MD approach). We are not aware of ab initio MD simulations for [Fe(OH)m]3−m(aq) (m > 0) systems, but a coordination number of 6 for the first hydration shell of the [Fe(OH)]2+(aq) system has been obtained by Rustad et al.16 (by applying Born–Mayer type potentials) in disagreement with the present results.

Before the properties of iron(III) in aqueous and alkaline environments with RxAQ potential are evaluated, the basic modell for water is examined. A RxAQ parameter set for the H–O interaction has been optimized, resulting in the following statistical properties for water at room temperature. Performing NVT MD simulations for up to 70 ps, one obtains the RDFs for O–H, O–O, and H–H as shown in Figs. 7(a)7(c). The main O–H peaks located at 2.63 and ∼4.6 Å are slightly shifted toward lower values compared to the ReaxFF results of Zhang and van Duin34 (with 2.67 and ∼4.4 Å) or experimental results of Soper50 (with 2.75 and ∼4.5 Å). Deviations for the O–O and H–H peaks are in the same range of magnitude.

FIG. 7.

O–H (a), O–O (b), and H–H (c) radial distribution functions (RDF) for water at a density of 1 g/cm3 and at a temperature of 300 K, as calculated from NVT MD simulations with RxAQ potentials, compared with experimental results of Soper50 (Exp) and ReaxFF calculations of Zhang and van Duin34 (Theo.).

FIG. 7.

O–H (a), O–O (b), and H–H (c) radial distribution functions (RDF) for water at a density of 1 g/cm3 and at a temperature of 300 K, as calculated from NVT MD simulations with RxAQ potentials, compared with experimental results of Soper50 (Exp) and ReaxFF calculations of Zhang and van Duin34 (Theo.).

Close modal

A NVE MD simulation has been performed for water. The behavior of the potential, kinetic, and total energies as well as the temperature with the run time is shown in Fig. 8. A drift of the total energy of 0.2 kcal/(mol ps) due to truncation errors is visible. A possibility for the reduction of the drift is pointed out by Furman and Wales.52,53

FIG. 8.

Monitoring of the total (Etot), potential (Epot), kinetic (Ekin) energies (in kcal/mol) and the temperature with the NVE MD run time for 125 H2O molecules. E0 = Etot(0) = −2.51 × 104 kcal/mol is the total energy of the NVE MD simulation for the time t = 0.

FIG. 8.

Monitoring of the total (Etot), potential (Epot), kinetic (Ekin) energies (in kcal/mol) and the temperature with the NVE MD run time for 125 H2O molecules. E0 = Etot(0) = −2.51 × 104 kcal/mol is the total energy of the NVE MD simulation for the time t = 0.

Close modal

Furthermore, the properties of aqueous systems with an excess or deficiency of H+ ions have been considered. The hopping functions for the H3O+ ion and the OH ion are shown in Fig. 9(a). The mean squared displacements (MSDs) for H2O, O of OH, and H+ of H3O+ are shown in Fig. 9(b). As a result of the applied computation procedure,51 the hopping functions and the MSDs show only small deviations from a linear behavior with time.

FIG. 9.

Grotthuss hopping function of H3O+ and OH for one trajectory (trj) and an average over 15 trajectories (avg) (a), MSDs for H+(H3O+), O(OH), and H2O (b) as well as approximate MSDs (aMSD) for H+(H3O+) and O(OH) and MSDs for H3O+(aq) and OH(aq) (c), as obtained from NVT MD simulations with RxAQ potentials.51 

FIG. 9.

Grotthuss hopping function of H3O+ and OH for one trajectory (trj) and an average over 15 trajectories (avg) (a), MSDs for H+(H3O+), O(OH), and H2O (b) as well as approximate MSDs (aMSD) for H+(H3O+) and O(OH) and MSDs for H3O+(aq) and OH(aq) (c), as obtained from NVT MD simulations with RxAQ potentials.51 

Close modal
Assuming linearity in time for the MSD as well as for the hopping functions and a decoupling of the true movement of potentially hopping atoms r̃A from the drift due to the virtual hopping, one can compute the approximate MSD (aMSD)
aMSD=[r̃A(t)r̃A(0)]2+h̄̇tdhp2,
(5)
where h̄̇ and dhp are the average hopping rate and the hopping distance, respectively, and ⟨·⟩ denotes an ensemble average. The derivation of Eq. (5) is presented in  Appendix C. The hopping distance, for both H3O+ and OH, is dhp = dOO ≃ 2.6 Å (i.e., the position of the first peak in the O–O RDF for water). Average hopping rates for H3O+ and OH can be computed from the corresponding hopping functions shown in Fig. 9. The first contribution [r̃A(t)r̃A(0)]2 in the aMSD is approximated by the MSD for the H3O+(aq) and OH(aq) systems.

For H+(H3O+) and O(OH), maximal deviations of 15% are found between MSD and aMSD for t > 10 ps. Furthermore, on comparing MSD and aMSD, one finds that the ratio between contributions from virtual hopping and true movement of atoms to the MSD for H+(H3O+) and O(OH) is about 2:1.

The hopping rates of h̄̇(H3O+)=0.42 ps−1 and h̄̇(OH)=0.30 ps−1 are of the correct tendency (i.e., h̄̇(H3O+)>h̄̇(OH)) and close to the results of Zhang and van Duin34 [with h̄̇(H3O+)=0.53ps1 and h̄̇(OH)=0.35ps1] or Wu et al.48 [with h̄̇(H*)=0.16ps1].

The diffusion coefficient D can be computed from the Einstein relation as
D=16ddtMSD(t).
(6)
In particular, for a system obeying the diffusion equation, the exponent of the power law between the moments54,55 |r(t)r(0)|ν (ν = 1, 2, 3, 4) and Dt (i.e., time × diffusion coefficient) is given as the slope of the respective log–log plot. It is shown ν = 2 (i.e., for the MSD) in the inset of Fig. 9(b). The deviations from a unit slope are visible only for times ≲1 ps.

The average diffusion coefficients of D(H2O) = 0.19 Å2/ps, D(O(OH)) = 0.52 Å2/ps, and D(H+(H3O+)) = 0.70 Å2/ps can be extracted. The values are also of the correct tendency [i.e., D(H2O) <D(O(OH)) <D(H+(H3O+))] and not far from the computations of Zhang and van Duin34(with 0.27, 0.78, and 1.04 Å2/ps by applying the ReaxFF method), Lee and Rasaiah56 (with 0.23, 0.46, and 0.76 Å2/ps by applying the OSS2 modell), and Wu et al.48 [with D(H2O) = 0.23…0.35 Å2/ps, and D(H+) = 0.29…0.50 Å2/ps using the MS-EVB3 modell]. The experimental values57 computed for infinite dilution34 are 0.24, 0.56, and 0.93 Å2/ps.

In summary, since the main trends for the water characteristics are correctly described with the present RxAQ parameter set for the H–O interaction, the corresponding water model is considered as acceptable to incorporate the H–O–Fe interactions in order to study the properties of iron(III) in aqueous and alkaline environments with RxAQ potentials.

The parameter set for the H–O–Fe interaction has been obtained by optimizing parameters for H–Fe and O–Fe two-body contributions and H–O–Fe three-body contributions of the RxAQ potentials (without modifying the H–O interaction). The RxAQ energies for proton transfer reactions in a gaseous environment are presented in Table I. Maximal deviations of up to 41% from the DFT reference data are found for the reactions [Fe(H2O)5OH]2+ + H2O → [Fe(H2O)4(OH)2]+ + H3O+. The energetic ordering of the final states also deviate from the ab initio results as trans-[Fe(H2O)4(OH)2]+ is the most stable geometry within the RxAQ approach. The average Fe–OH2O and Fe–OOH distances, presented in Table II, for the geometries [Fe(H2O)n(OH)m]3−m are close to the DFT values with a maximal deviation of 0.04 Å. The increasing Fe–O distances for the decreasing total charge are also in agreement with the DFT results. For the RxAQ charges on the atoms, one observes the same trends as for DFT–Bader charges for oxygen and hydrogen. The charges on iron are overestimated within the RxAQ approach by a max. of 26% and with the maximal charge on iron for the [Fe(H2O)5(OH)]2+ complex (with 1.75 |e|) in a part disagreement with the ab initio results.

In aqueous and alkaline environments, the energies of the reactions R0…R2 are higher than the corresponding ab initio values in 4 of 5 cases (cf. Table III and Fig. 1). No activation energy has been found for the reaction R1-i in disagreement with the DFT results. For the reactions R1-ii and R2, the RxAQ activation energies (with 16.7 and 11.6 kcal/mol) are larger than the corresponding DFT values (with 11.6 and 6.8 kcal/mol). Once again, it should be mentioned that since the convergence criteria have been relaxed for the RxAQ approach, only approximate values and paths are reported.

NVT MD simulations with RxAQ potentials have been performed for 100 ps. To compare the RxAQ results with the ab initio data, a similar supercell containing 64 water molecules has been used. The Fe–O radial distribution functions and coordination numbers for the [Fe(OHm)](3−m) systems resulting from the RxAQ NVT MD simulations are shown in Fig. 10. Differences from the DFT results can be detected. For Fe3+ (aq), one finds a single peak at 2.13 Å but the first shell coordination number is 5. For [Fe(OHm)](3−m) (m > 0), the first shell coordination numbers are in agreement with the DFT results, but two peaks between 1.83 and 2.18 Å (roughly corresponding to the average Fe–OOH and Fe–OH2O distances) are visible. The reason is the larger separation of the average Fe–OOH and Fe–OH2O distances by about 0.3 Å compared to DFT values of 0.1…0.2 Å (cf. Table IV).

FIG. 10.

Fe–O (a) and Fe–H (b) radial distribution functions (RDF) and accumulated Fe–O coordination numbers (CN) (c) as calculated from the NVT MD simulations with RxAQ potentials for the systems [Fe(OH)m]3−m (aq) (m = 0…3).

FIG. 10.

Fe–O (a) and Fe–H (b) radial distribution functions (RDF) and accumulated Fe–O coordination numbers (CN) (c) as calculated from the NVT MD simulations with RxAQ potentials for the systems [Fe(OH)m]3−m (aq) (m = 0…3).

Close modal
As for ab initio NVT MD, the number of OH groups around the Fe3+ ion has been monitored for the RxAQ NVT MD. The results for the NVT MD simulations with temperatures of 300 and 350 K are shown in Fig. 11. Compared to the ab initio NVT MD, far less proton transitions are visible when applying RxAQ potentials. A slight increase in proton hopping events is observed when the temperature is increased to 350 K in qualitative agreement with the Arrhenius law (or transition state theory), assuming that energy barriers are easier to overcome with higher temperatures. The lower proton transfer frequencies applying RxAQ potentials might be related to the mostly larger reaction energies and activation barriers for the reaction R0…R2 within the RxAQ approach. In particular, the NVT MD simulation for [Fe(OH)2]+(aq) and the reaction R2 can be tentatively related due to equal first shell coordination numbers. The minimal difference in activation energies for reaction R2 within DFT and RxAQ can be estimated from NVT MD simulations. The reaction constant kDFT for the reaction m → m + 1 can be computed by the average time the [Fe(OH)2]+(aq) system remains with m = 2 until a transition to m = 2 + 1 occurs, and with the assumption that the time the system remains with m ≠ 2 is negligible. One finds
kAI=#(mm+1)/TMD,
(7)
where #(mm + 1) ≈ 30 is the number of transitions m → m + 1, and TMD = 14 ps is the corresponding MD simulation time. For RxAQ, it is estimated that kRxAQ > 1/100 ps−1 for T = 300 K, as shown in Fig. 11, (i.e., m → m + 1 transitions have not been observed during the NVT MD simulation). Using the Arrhenius law (and assuming a similar pre-exponential factor for DFT and RxAQ), one obtains
ΔEA=EA(RxAQ)EA(AI)>ln(kAI/kRxAQ)/kBT8.9kcalmol.
(8)
FIG. 11.

Monitoring of the number of OH groups around the Fe3+ ion with the NVT MD run time (for every 50 fs) using RxAQ potentials for the [Fe(OH)m]3−m (aq) (m = 0…3) systems at temperatures of 300 and 350 K.

FIG. 11.

Monitoring of the number of OH groups around the Fe3+ ion with the NVT MD run time (for every 50 fs) using RxAQ potentials for the [Fe(OH)m]3−m (aq) (m = 0…3) systems at temperatures of 300 and 350 K.

Close modal

With the estimated value, a lower threshold for ΔEA is given. The result is slightly larger than the value of 4.8 kcal/mol, which can be extracted from the corresponding NEB calculations (cf. Fig. 1 and Table III). The comparison presented here can be considered as a check for the orders of magnitude.

The iron(III) complexes [Fe(H2O)n(OH)m]3−m (n + m = 5, 6, m ⇐ 3) and related proton transfer reactions have been studied using ab initio methods in the first part. The total energy calculations for the reactions [Fe(H2O)n(OH)m]3−m + H2O → [Fe(H2O)n−1(OH)m+1]2−m + H3O+ revealed that all reactions are exothermic in a gaseous environment. Evaluating the geometrical properties and the atomic Bader charges for the iron(III) complexes, it is observed that Fe–O distances are enhanced and mostly a charge separation between iron and oxygen occurs with an increasing number of OH groups. In an aqueous or alkaline environment, the proton transfer reactions are endothermic with energies between 3.7 and 10.8 kcal/mol. With the application of the nudge elastic band method, activation energies of up to 14.6 kcal/mol and a collective proton transfer with at least two involved protons have been found. NVT MD simulations have been performed for [Fe(OH)m]3−m(aq) systems (m = 0, …, 3). By monitoring the number of hydroxyl groups around the iron(III) ion in the trajectory, it can be concluded that a permanent proton transfer in the first hydration shell occurs. Fe–O and Fe–H radial distribution functions and distance-dependent coordination numbers have been evaluated. One finds a coordination number of six for the first hydration shell of Fe3+(aq). The first hydration shell coordination number is decreasing to 5, 5, and 4 for the systems [Fe(OH)]2+(aq), [Fe(OH)2]1+(aq), and [Fe(OH)3](aq) [i.e., with the number of hydroxyl groups around the iron(III) ion], respectively.

Similar investigation methods combined with ReaxFF-AQ potentials have been applied to the iron(III) complexes in the second part. The corresponding O–H–Fe interaction within the ReaxFF-AQ approach has been adjusted by mainly using the present ab initio data. Comparing the ab initio and the ReaxFF-AQ results, one finds quantitative differences, but many trends remain conserved. The latter applies to the Fe–O distance changes and atomic charge accumulations with an increasing number of hydroxyl groups around the iron(III) ion as well as to the energies of proton transfer reactions in gaseous and aqueous/alkaline environments. Larger deviations between ab initio and ReaxFF-AQ predictions can be found evaluating the statistical properties of the molecular dynamics simulations for iron(III) in aqueous/alkaline environments. The first hydration shell coordination number of Fe3+(aq) is five within ReaxFF-AQ NVT MD (in disagreement with the ab initio predictions). Furthermore, the proton transfer frequency in the first hydration shell is clearly lower for ReaxFF-AQ potentials than for ab initio potentials. Nevertheless, with this study, it can be demonstrated that it is partly possible to describe the characteristics of iron(III) in aqueous/alkaline environments with ReaxFF-AQ potentials or that the AI-PESs of the [Fe(OH)m]3−m(aq) (m = 0, …, 3) systems are partly correctly modeled with ReaxFF-AQ.

Three possible reasons for the deviations between the ab initio and the ReaxFF-AQ results can be given: (I) an insufficient or inappropriate sampling of the AI-PESs could have been applied. (II) The optimized parameter set for ReaxFF-AQ might be not accurate enough to describe the considered properties of AI-PESs. (III) The ReaxFF-AQ model could, in general, not be accurate enough to describe the multidimensional AI-PESs for [Fe(OH)m]3−m(aq) (m = 0, …, 3).

The methodological extensions for the ReaxFF-AQ approach can be conducted from the reasons for the deviations: (I) certain aspects concerning systematic techniques for sampling and interpolation of AI-PESs58,59 can be adopted. (II) The parameter optimization procedure can be extended or modern methods can be explored. (III) In addition, four body contributions and further degrees of freedom can be considered in the ReaxFF-AQ model.

From a broader perspective, particularly proton transfer nuclear quantum effects (NQEs), such as tunneling and zero-point energy, may play an important role.60,61 Investigations performed by Ganeshan et al.62 applying path integral MD in combination with ReaxFF, density functional tight binding and NequIP neural network potentials showed that the diffusion coefficient of water is influenced as well as proton transfer barriers in bulk and under confinement in Ti3C2 MXenes are reduced by NQEs.

This project was funded by the Federal Ministry for Economic Affairs and Climate Action (Grant Reference Nos. ZF4792901DF9 and ZF4105821DF9), following a decision of the German Bundestag.

The authors have no conflicts to disclose.

Arthur Riefer: Conceptualization (equal); Data curation (lead); Investigation (lead); Methodology (equal); Software (lead); Writing – original draft (lead). Matthias Hackert-Oschätzchen: Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (supporting). Philipp Plänitz: Conceptualization (equal); Methodology (equal); Project administration (supporting); Software (supporting); Writing – review & editing (equal). Gunnar Meichsner: Project administration (supporting); Resources (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The correction energy Ecorr within RxAQ is defined as
Ecorr=Ebc+E3BII.
(A1)
The semi-parabolic two-body correction energy Ebc is defined as
Ebc=I[AI(BIBI0)]2BI2CI2+1Θ(BIBI0)+Ebc,I0,
(A2)
where I = (ij) and Θ is an analytical approximation of the step function. The additional parameters AI, BI0, CI, and Ebc,I0 per atom pair (ij) can be optimized. The parabolic three-body correction energy E3BII is defined as
E3BII=TAI1Θ(BI10BI1)+AI1+Θ(BI1BI10)BI12CI12+1(BI1BI10)2+AI2Θ(BI20BI2)+AI2+Θ(BI2BI20)BI22CI22+1(BI2BI20)2+cos2(θTθT0)ATBI1BI2BI1BI2+CT2+E3BII,T0,
(A3)
with T = (I1, I2), I1 = (ij), and I2 = (jk). The additional parameters AI1, AI1+, BI10, CI1, AI2, AI2+, BI20, CI2, θT0, AT, CT, and E3BII,T0 can be optimized.
The charges on the atoms within RxAQ are determined by solving the following optimization problem:
E({qĩ})Ech({qĩ})+Ecoulomb({qĩ})min.,
(B1)
with the effective charges q̃i=qiQN and the constraint ∑iqi = Q. The charge modifications qiqiQN are applied for structures with atom numbers N > 10.
The approximate MSD can be computed from the following assumptions: (I) linearity in time for the MSD for the true movement r̃A of potentially hopping atoms A, (II) linearity in time for the hopping function h(t), and (III) a decoupling of the true movement r̃A from the virtual drift δ(t)eδ(t) [with the drift direction eδ(t)] due to hopping,
MSD(t)aMSD(t)=
(C1)
[r̃A(t)r̃A(0)+δ(t)eδ(t)]2.
(C2)
Using the Taylor expansion up to the first order in Δt in all functions [i.e., f(t + Δt) ≈ f(t) + tf(t) · Δt], one finds for the difference quotient,
1Δt[aMSD(t+Δt)aMSD(t)]=const.=2(r̃A(t)r̃A(0))dr̃Adt(t)+2ddt[(r̃A(t)r̃A(0))eδ(t)δ(t)]+2δ(t)dδdt(t).
(C3)
The middle contribution is proportional to |(r̃A(t)r̃A(0))|ϕ(t) (or to its time derivative) with the phase
ϕ(t)=r̃A(t)r̃A(0)|(r̃A(t)r̃A(0))|eδ(t).
(C4)
Contributions of that kind are neglected in the present derivation.
The first contribution on the right-hand-side (rhs) of Eq. (C3) is related to the MSD of the true movement of atoms and is a constant. The consequence is that the third term on the rhs of Eq. (C3) must be constant, i.e.,
2δ(t)dδdt(t)=C=const.
(C5)
The solution of the ordinary differential equation is
δ(t)2=Ct+C0.
(C6)
The two constants C0 and C can be determined from the initial condition δ(0) = 0 and from consideration of the average time for one hopping event 1/h̄̇ and the hopping distance of dhp. The result is Eq. (5) after an expansion of Eq. (C2).
1.
O.
Böhm
,
S.
Pfadenhauer
,
R.
Leitsmann
,
P.
Plänitz
,
E.
Schreiner
, and
M.
Schreiber
, “
ReaxFF+—A new reactive force field method for the accurate description of ionic systems and its application to the hydrolyzation of aluminosilicates
,”
J. Phys. Chem. C
120
,
10849
10856
(
2016
).
2.
A.
Riefer
,
M.
Hackert-Oschätzchen
,
P.
Plänitz
, and
G.
Meichsner
, “
Derivation of parameter sets for the ReaxFF+ method for modeling an electrochemical machining process
,”
Procedia CIRP
117
,
231
236
(
2023
).
3.
C.
Rosenkranz
,
M.
Lohrengel
, and
J.
Schultze
, “
The surface structure during pulsed ECM of iron in NaNO3
,”
Electrochim. Acta
50
,
2009
2016
(
2005
).
4.
B. S.
Brunschwig
,
C.
Creutz
,
D. H.
Macartney
,
T.
Sham
, and
N.
Sutin
, “
The role of inner-sphere configuration changes in electron-exchange reactions of metal complexes
,”
Faraday Discuss. Chem. Soc.
74
,
113
127
(
1982
).
5.
P.
D’Angelo
and
M.
Benfatto
, “
Effect of multielectronic configurations on the XAFS analysis at the Fe K edge
,”
J. Phys. Chem. A
108
,
4505
4514
(
2004
).
6.
M.
Benfatto
,
J. A.
Solera
,
J.
García Ruiz
, and
J.
Chaboy
, “
Double-channel excitation in the x-ray absorption spectrum of Fe3+ water solutions
,”
Chem. Phys.
282
,
441
450
(
2002
).
7.
J. E.
Enderby
, “
Ion solvation via neutron scattering
,”
Chem. Soc. Rev.
24
,
159
168
(
1995
).
8.
G.
Herdman
and
G.
Neilson
, “
Ferrous Fe(II) hydration in a 1 molal heavy water solution of iron chloride
,”
J. Phys.: Condens. Matter
4
,
649
(
1992
).
9.
D.
Powell
and
G.
Neilson
, “
The concentration dependence of the Ni2+ hydration geometry in aqueous solution
,”
J. Phys.: Condens. Matter
2
,
3871
(
1990
).
10.
R. L.
Martin
,
P. J.
Hay
, and
L. R.
Pratt
, “
Hydrolysis of ferric ion in water and conformational equilibrium
,”
J. Phys. Chem. A
102
,
3565
3573
(
1998
).
11.
H. A.
De Abreu
,
L.
Guimarães
, and
H. A.
Duarte
, “
Density-functional theory study of iron(III) hydrolysis in aqueous solution
,”
J. Phys. Chem. A
110
,
7713
7718
(
2006
).
12.
G.
Ottonello
and
M. V.
Zuccolini
, “
Ab-initio structure, energy and stable Fe isotope equilibrium fractionation of some geochemically relevant H–O–Fe complexes
,”
Geochim. Cosmochim. Acta
73
,
6447
6469
(
2009
).
13.
S. A.
Bogatko
,
E. J.
Bylaska
, and
J. H.
Weare
, “
First principles simulation of the bonding, vibrational, and electronic properties of the hydration shells of the high-spin Fe3+ ion in aqueous solutions
,”
J. Phys. Chem. A
114
,
2189
2200
(
2010
).
14.
S.
Mandal
,
R.
Kar
,
B.
Meyer
, and
N. N.
Nair
, “
Hybrid functional and plane waves based ab initio molecular dynamics study of the aqueous Fe2+/Fe3+ redox reaction
,”
ChemPhysChem
24
,
e202200617
(
2023
).
15.
S. T.
Moin
,
T. S.
Hofer
,
A. B.
Pribil
,
B. R.
Randolf
, and
B. M.
Rode
, “
A quantum mechanical charge field molecular dynamics study of Fe2+ and Fe3+ ions in aqueous solutions
,”
Inorg. Chem.
49
,
5101
5106
(
2010
).
16.
J. R.
Rustad
,
B. P.
Hay
, and
J.
Halley
, “
Molecular dynamics simulation of iron(III) and its hydrolysis products in aqueous solution
,”
J. Chem. Phys.
102
,
427
431
(
1995
).
17.
J.
Tomasi
and
M.
Persico
, “
Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent
,”
Chem. Rev.
94
,
2027
2094
(
1994
).
18.
L. R.
Pratt
,
G.
Hummer
, and
A. E.
Garcia´
, “
Ion pair potentials-of-mean-force in water
,”
Biophys. Chem.
51
,
147
165
(
1994
).
19.
L. E.
Ratcliff
,
S.
Mohr
,
G.
Huhs
,
T.
Deutsch
,
M.
Masella
, and
L.
Genovese
, “
Challenges in large scale quantum mechanical calculations
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
7
,
e1290
(
2017
).
20.
A. C.
van Duin
,
S.
Dasgupta
,
F.
Lorant
, and
W. A.
Goddard
, “
ReaxFF: A reactive force field for hydrocarbons
,”
J. Phys. Chem. A
105
,
9396
9409
(
2001
).
21.
A. C.
van Duin
,
A.
Strachan
,
S.
Stewman
,
Q.
Zhang
,
X.
Xu
, and
W. A.
Goddard
, “
ReaxFFSiO reactive force field for silicon and silicon oxide systems
,”
J. Phys. Chem. A
107
,
3803
3811
(
2003
).
22.
A. C.
van Duin
,
C.
Zou
,
K.
Joshi
,
V.
Bryantsev
, and
W. A.
Goddard
, “
A Reaxff reactive force-field for proton transfer reactions in bulk water and its applications to heterogeneous catalysis
,”
Comput. Catal.
14
,
223
(
2013
).
23.
T. P.
Senftle
,
S.
Hong
,
M. M.
Islam
,
S. B.
Kylasa
,
Y.
Zheng
,
Y. K.
Shin
,
C.
Junkermeier
,
R.
Engel-Herbert
,
M. J.
Janik
,
H. M.
Aktulga
et al, “
The ReaxFF reactive force-field: Development, applications and future directions
,”
npj Comput. Mater.
2
,
15011
(
2016
).
24.
C. M.
Flynn
, Jr.
, “
Hydrolysis of inorganic iron(III) salts
,”
Chem. Rev.
84
,
31
41
(
1984
).
25.
A.
Meagher
, “
Ferric hydrolysis in water: An iron-57 mössbauer study using iron-exchanged nafion
,”
Inorg. Chim. Acta
146
,
19
23
(
1988
).
26.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
27.
A.
García
,
N.
Papior
,
A.
Akhtar
,
E.
Artacho
,
V.
Blum
,
E.
Bosoni
,
P.
Brandimarte
,
M.
Brandbyge
,
J. I.
Cerdá
,
F.
Corsetti
,
R.
Cuadrado
,
V.
Dikan
,
J.
Ferrer
,
J.
Gale
,
P.
García-Fernández
,
V. M.
García-Suárez
,
S.
García
,
G.
Huhs
,
S.
Illera
,
R.
Korytár
,
P.
Koval
,
I.
Lebedeva
,
L.
Lin
,
P.
López-Tarifa
,
S. G.
Mayo
,
S.
Mohr
,
P.
Ordejón
,
A.
Postnikov
,
Y.
Pouillon
,
M.
Pruneda
,
R.
Robles
,
D.
Sánchez-Portal
,
J. M.
Soler
,
R.
Ullah
,
V. W.-z.
Yu
, and
J.
Junquera
, “
Siesta: Recent developments and applications
,”
J. Chem. Phys.
152
,
204108
(
2020
).
28.
J. M.
Soler
,
E.
Artacho
,
J. D.
Gale
,
A.
García
,
J.
Junquera
,
P.
Ordejón
, and
D.
Sánchez-Portal
, “
The SIESTA method for ab initio order-N materials simulation
,”
J. Phys.: Condens. Matter
14
,
2745
2779
(
2002
).
29.
N.
Troullier
and
J. L.
Martins
, “
Efficient pseudopotentials for plane-wave calculations
,”
Phys. Rev. B
43
,
1993
2006
(
1991
).
30.
A.
Hjorth Larsen
,
J.
Jørgen Mortensen
,
J.
Blomqvist
,
I. E.
Castelli
,
R.
Christensen
,
M.
Dułak
,
J.
Friis
,
M. N.
Groves
,
B.
Hammer
,
C.
Hargus
,
E. D.
Hermes
,
P. C.
Jennings
,
P.
Bjerre Jensen
,
J.
Kermode
,
J. R.
Kitchin
,
E.
Leonhard Kolsbjerg
,
J.
Kubal
,
K.
Kaasbjerg
,
S.
Lysgaard
,
J.
Bergmann Maronsson
,
T.
Maxson
,
T.
Olsen
,
L.
Pastewka
,
A.
Peterson
,
C.
Rostgaard
,
J.
Schiøtz
,
O.
Schütt
,
M.
Strange
,
K. S.
Thygesen
,
T.
Vegge
,
L.
Vilhelmsen
,
M.
Walter
,
Z.
Zeng
, and
K. W.
Jacobsen
, “
The atomic simulation environment—A Python library for working with atoms
,”
J. Phys.: Condens. Matter
29
,
273002
(
2017
).
31.
G.
Henkelman
and
H.
Jónsson
, “
Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points
,”
J. Chem. Phys.
113
,
9978
9985
(
2000
).
32.
G.
Henkelman
,
B. P.
Uberuaga
, and
H.
Jónsson
, “
A climbing image nudged elastic band method for finding saddle points and minimum energy paths
,”
J. Chem. Phys.
113
,
9901
9904
(
2000
).
33.
S.
Nosé
, “
A unified formulation of the constant temperature molecular dynamics methods
,”
J. Chem. Phys.
81
,
511
519
(
1984
), https://pubs.aip.org/aip/jcp/article-pdf/81/1/511/9722086/511_1_online.pdf.
34.
W.
Zhang
and
A. C.
van Duin
, “
Second-generation ReaxFF water force field: Improvements in the description of water density and OH-anion diffusion
,”
J. Phys. Chem. B
121
,
6021
6032
(
2017
).
35.
C. E.
Wilmer
,
K. C.
Kim
, and
R. Q.
Snurr
, “
An extended charge equilibration method
,”
J. Phys. Chem. Lett.
3
,
2506
2511
(
2012
).
36.
U. W.
Schmitt
and
G. A.
Voth
, “
The computer simulation of proton transport in water
,”
J. Chem. Phys.
111
,
9361
9381
(
1999
).
37.
E. E.
Dahlke
,
M. A.
Orthmeyer
, and
D. G.
Truhlar
, “
Assessment of multicoefficient correlation methods, second-order Møller–Plesset perturbation theory, and density functional theory for H3O+(H2O)n (n = 1–5) and OH(H2O)n (n = 1–4)
,”
J. Phys. Chem. B
112
,
2372
2381
(
2008
).
38.
M. E.
Tuckerman
,
D.
Marx
, and
M.
Parrinello
, “
The nature and transport mechanism of hydrated hydroxide ions in aqueous solution
,”
Nature
417
,
925
929
(
2002
).
39.
W.
Tang
,
E.
Sanville
, and
G.
Henkelman
, “
A grid-based Bader analysis algorithm without lattice bias
,”
J. Phys.: Condens. Matter
21
,
084204
(
2009
).
40.
E.
Sanville
,
S. D.
Kenny
,
R.
Smith
, and
G.
Henkelman
, “
Improved grid-based algorithm for Bader charge allocation
,”
J. Comput. Chem.
28
,
899
908
(
2007
).
41.
G.
Henkelman
,
A.
Arnaldsson
, and
H.
Jónsson
, “
A fast and robust algorithm for Bader decomposition of charge density
,”
Comput. Mater. Sci.
36
,
354
360
(
2006
).
42.
M.
Yu
and
D. R.
Trinkle
, “
Accurate and efficient algorithm for Bader charge integration
,”
J. Chem. Phys.
134
,
064111
(
2011
).
43.
C.
Fonseca Guerra
,
J.-W.
Handgraaf
,
E. J.
Baerends
, and
F. M.
Bickelhaupt
, “
Voronoi deformation density (VDD) charges: Assessment of the Mulliken, Bader, Hirshfeld, Weinhold, and VDD methods for charge analysis
,”
J. Comput. Chem.
25
,
189
210
(
2004
).
44.
M.
Aryanpour
,
A. C. T.
van Duin
, and
J. D.
Kubicki
, “
Development of a reactive force field for iron-oxyhydroxide systems
,”
J. Phys. Chem. A
114
,
6298
6307
(
2010
).
45.
M.
Dittner
,
J.
Müller
,
H. M.
Aktulga
, and
B.
Hartke
, “
Efficient global optimization of reactive force-field parameters
,”
J. Comput. Chem.
36
,
1550
1561
(
2015
).
46.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in ’t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
,
S. J.
Plimpton
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
47.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
48.
Y.
Wu
,
H.
Chen
,
F.
Wang
,
F.
Paesani
, and
G. A.
Voth
, “
An improved multistate empirical valence bond model for aqueous proton solvation and transport
,”
J. Phys. Chem. B
112
,
467
482
(
2008
).
49.
A. D.
Becke
, “
Density-functional thermochemistry. III. The role of exact exchange
,”
J. Chem. Phys.
98
,
5648
5652
(
1993
), https://pubs.aip.org/aip/jcp/article-pdf/98/7/5648/11091662/5648_1_online.pdf.
50.
A.
Soper
, “
The radial distribution functions of water and ice from 220 to 673 K and at pressures up to 400 MPa
,”
Chem. Phys.
258
,
121
137
(
2000
).
51.

The hopping functions and the MSDs have been computed with the following procedure. For water, H3O+(aq) and OH(aq) supercells with 64H2O, H3O+/63H2O and OH/63H2O in the NVT MD and corresponding geometry snapshots for every 50 fs have been considered. For water, H3O+(aq) and OH(aq) average properties for 1, 15 and 15 trajectories have been computed. Furthermore, the MSDs have been computed from unfolded NVT MD trajectories for a MD run time of 160 ps considering a maximal time lag of 80 ps.

52.
D.
Furman
and
D. J.
Wales
, “
Transforming the accuracy and numerical stability of ReaxFF reactive force fields
,”
J. Phys. Chem. Lett.
10
,
7215
7223
(
2019
).
53.
D.
Furman
and
D. J.
Wales
, “
A well-behaved theoretical framework for ReaxFF reactive force fields
,”
J. Chem. Phys.
153
,
021102
(
2020
).
54.
P.
Heitjans
and
J.
Kärger
,
Diffusion in Condensed Matter: Methods, Materials, Models
(
Springer
,
Berlin, Heidelberg
,
2006
), ISBN: 9783540309703, SpringerLink: Springer e-Books.
55.
S.
Fritzsche
,
R.
Haberlandt
,
J.
Kärger
,
H.
Pfeifer
, and
K.
Heinzinger
, “
An MD simulation on the applicability of the diffusion equation for molecules adsorbed in a zeolite
,”
Chem. Phys. Lett.
198
,
283
287
(
1992
).
56.
S. H.
Lee
and
J. C.
Rasaiah
, “
Proton transfer and the mobilities of the H+ and OH ions from studies of a dissociating model for water
,”
J. Chem. Phys.
135
,
124505
(
2011
).
57.
T. S.
Light
,
S.
Licht
,
A. C.
Bevilacqua
, and
K. R.
Morash
, “
The fundamental conductivity and resistivity of water
,”
Electrochem. Solid-State Lett.
8
,
E16
(
2005
).
58.
M. A.
Collins
, “
Molecular potential-energy surfaces for chemical reaction dynamics
,”
Theor. Chem. Acc.
108
,
313
324
(
2002
).
59.
H.-Y.
Kwon
,
Z.
Morrow
,
C.
Kelley
, and
E.
Jakubikova
, “
Interpolation methods for molecular potential energy surface construction
,”
J. Phys. Chem. A
125
,
9725
9735
(
2021
).
60.
T.
Miyazaki
,
Atom Tunneling Phenomena in Physics, Chemistry and Biology
(
Springer Science & Business Media
,
2004
), Vol.
36
.
61.
T. E.
Markland
and
M.
Ceriotti
, “
Nuclear quantum effects enter the mainstream
,”
Nat. Rev. Chem.
2
,
0109
(
2018
).
62.
K.
Ganeshan
,
R.
Khanal
,
M. G.
Muraleedharan
,
M.
Hellström
,
P. R.
Kent
,
S.
Irle
, and
A. C.
van Duin
, “
Importance of nuclear quantum effects on aqueous electrolyte transport under confinement in Ti3C2 MXenes
,”
J. Chem. Theory Comput.
18
,
6920
6931
(
2022
).