Periodic density functional theory calculations of bulk and the (010) surface of goethite

Background Goethite is a common and reactive mineral in the environment. The transport of contaminants and anaerobic respiration of microbes are significantly affected by adsorption and reduction reactions involving goethite. An understanding of the mineral-water interface of goethite is critical for determining the molecular-scale mechanisms of adsorption and reduction reactions. In this study, periodic density functional theory (DFT) calculations were performed on the mineral goethite and its (010) surface, using the Vienna Ab Initio Simulation Package (VASP). Results Calculations of the bulk mineral structure accurately reproduced the observed crystal structure and vibrational frequencies, suggesting that this computational methodology was suitable for modeling the goethite-water interface. Energy-minimized structures of bare, hydrated (one H2O layer) and solvated (three H2O layers) (010) surfaces were calculated for 1 × 1 and 3 × 3 unit cell slabs. A good correlation between the calculated and observed vibrational frequencies was found for the 1 × 1 solvated surface. However, differences between the 1 × 1 and 3 × 3 slab calculations indicated that larger models may be necessary to simulate the relaxation of water at the interface. Comparison of two hydrated surfaces with molecularly and dissociatively adsorbed H2O showed a significantly lower potential energy for the former. Conclusion Surface Fe-O and (Fe)O-H bond lengths are reported that may be useful in surface complexation models (SCM) of the goethite (010) surface. These bond lengths were found to change significantly as a function of solvation (i.e., addition of two extra H2O layers above the surface), indicating that this parameter should be carefully considered in future SCM studies of metal oxide-water interfaces.


Introduction
Goethite (α-FeOOH) is a common and reactive mineral in the environment [1]. α-FeOOH is the most thermodynamically-stable form of the Fe-oxyhydroxides found in soils, groundwater, and acid mine drainage precipitates [2]. α-FeOOH is an excellent adsorbent of contaminants (e.g., arsenate) and nutrients (e.g., phosphate) and is an electron receptor for anaerobic bacterial respiration under anoxic conditions ( [3,4] and references therein). Consequently, the surface chemistry of α-FeOOH is important for understanding many environmental processes. One of the most stable and well-studied low-index surfaces of α-FeOOH is the (010) surface (Pbnm space group = (100) in the Pnma space group; [5,6]). As a result, this study will focus on the (010) α-FeOOH surface.
Several molecular modeling studies have been published related to α-FeOOH and goethite-H 2 O interfaces. For example, Rosso and Rustad [7] published local density approximation (LDA) and generalized gradient approximation (GGA) density functional theory (DFT) calculations on the structures of diaspore, goethite, boehmite, lepidocrocite, akaganeite, guyanaite and grimaldiite that reproduced observed structures of these minerals to within 3%. Rustad et al. [8] used a force field that allowed for the dissociation of H 2 O, in order to perform classical molecular dynamics (MD) simulations of the (110) α-FeOOH-H 2 O interface. Rakovan et al. [5] combined atomic force microscopy, X-ray photoelectron spectroscopy (XPS), low-energy electron diffraction, and periodic molecular orbital calculations, and proposed that the most stable termination of the (010) α-FeOOH surface was the O2-termination. In the O2-termination, the surface is cleaved between protonated O atoms rather than the O1 non-protonated O atoms of the bulk structure ( Fig. 1). Shroll and Stratsmaa [9] performed classical MD simulations with an AMBER-type force field to examine goethite-water interactions and found a mineral influence on the water structure out to 5.5 Å. Kerisit et al. [6] employed a polarizable version of the SPC/E water model, along with a α-FeOOH-H 2 O force field [10], to describe the structure of electrolyte solutions near the (010) α-FeOOH surface. Most recently, Aquino et al. [11] performed DFT and Møller-Plesset calculations to investigate the H-bonding of H 2 O to clusters of the (110) α-FeOOH surface.
The calculations reported here are different from these previous studies in that they are periodic DFT energy minimizations of bulk goethite and the (010) α-FeOOHwater interface. The use of quantum chemistry to model the surface distinguishes this work from previous classical simulations. Furthermore, the use of periodic boundary conditions separates these results from those of Aquino et al. [11]. Energy minimizations performed using periodic DFT can be used to verify the accuracy of the computational approach before expending dramatically more computational resources to perform DFT-based MD simulations of the α-FeOOH-water interface.
The purposes of this study were to test the ability of periodic DFT to: (1) reproduce the structures and vibrational frequencies of bulk α-FeOOH, (2) predict the surface bond distances for use in surface complexation models, and (3) investigate the H-bonding of H 2 O at the α-FeOOH- Fe-bearing minerals present a significant challenge for DFT calculations because the electronic ground-state of Fe is typically in a high-spin state. Furthermore, α-FeOOH is anti-ferromagnetic and the spin states of individual Fe atoms within a model must be carefully controlled. Once the DFT methodology, as implemented in VASP, can be shown to reliably reproduce static observables, such as crystal structures and vibrational frequencies, then similar computational methods can be employed to perform quantum MD simulations of reactions at the mineral-H 2 O interface.

Methods
Bulk α-FeOOH and the (010) surface models were built using the Crystal Builder and Surface Builder modules of Cerius 2 4.9 [12], respectively. The α-FeOOH lattice parameters and initial atomic positions were based upon the previously published experimental X-ray diffraction measurements of Szytula et al. [13]. It is imperative to note that the Cerius 2 4.9 program uses the Pnma space group for α-FeOOH. In the Pbnm space group the (010) surface, which is commonly used for α-FeOOH [13], is the (100) surface in the Pnma space group (P.J. Heaney, pers. comm.).
The (010) α-FeOOH surface employed in this study was cleaved from the (100) plane of α-FeOOH in the Pnma space group with the O2-termination (Fig. 1). This surface structure is consistent with that of Kendall et al. [14]. The 2-D periodic slab generated by this process was stoichiometric, neutral, and symmetric. The bare surface consisted of bridging OH groups (Fe 2 OH or μ-hydroxyl), bridging oxo groups (Fe 3 O or μ 3 -oxo), and 5-fold coordinated Fe atoms. Thus, to hydrate the surface and keep it neutral, H 2 O molecules were coordinated to the surface Fe atoms to fulfill octahedral coordination. In this configuration, the hydrated model was neutral and the presence of Fe-OH 2 functional groups did not indicate a positivelycharged surface of protonated Fe-OH functional groups, as is commonly assumed for metal oxide surfaces. The question of whether or not H 2 O favorably dissociates on this surface was also examined in this study. The O2 atoms are unlikely to accept a second proton (i.e., form an Fe 2 OH 2 + site), so if the (Fe)-OH 2 groups were to donate a proton it would be transferred to the O atoms bonded to three Fe atoms just below the α-FeOOH-water interface (i.e., Fe 3 O + H + → Fe 3 OH + ). The relative potential energies of these two configurations were compared. If the O1-termination were used, then dissociation would be more probable because the 5-coordinate surface Fe atoms could adsorb an H 2 O molecule and this could dissociate to form Fe-OH and protonated O1 surface atoms (i.e., FeOH 2 + FeO → 2FeOH). However, since Rakovan et al. [5] determined that this termination was less stable, we did not model it in our study.
All calculations were performed with the Vienna Ab Initio Simulation Package (VASP; [15]). The projector-augmented wave (PAW) method [16,17] was used in combination with the Perdew, Burke, and Ernzerhof (PBE) [18] exchange-correlation functional. The O_h and H_h pseudopotentials, as implemented in VASP, were used for the O and H atoms, respectively. The Fe_pv pseudopotential, which includes 14 electrons in the valence shell and treats the 3-p electrons explicitly, was used for the Fe atoms. The plane-wave kinetic energy cutoff was set to 700 eV. (Note: Softer pseudopotentials and a lower energy cutoff were tested and resulted in much poorer reproduction of the observed goethite crystal structure.) Energy differences of 1 × 10 -4 eV/atom and energy difference gradients of -0.02 eV/Å were used as convergence criteria. The number of unpaired electrons for the Fe atoms were set to five (i.e., high-spin d 5 electronic configuration), with alternating (010) planes in positive and negative spin directions ( Fig.  1; [19]). This ensured that the model mimicked the antiferromagnetic ground-state of α-FeOOH and that the overall magnetization remained close to zero.
Two bulk α-FeOOH simulation cells were used, namely a 1 × 1 × 1 unit cell model and a 1 × 3 × 2 supercell model. P1 symmetry was applied in both cases. The Monkhorst-Pack [20] scheme was used to generate the k-point sampling grids within the Brillouin zone. For the bulk calculations, a 2 × 6 × 4 grid was used for the 1 × 1 × 1 unit cell model because the unit cell was approximately 9.9 × 3.0 × 4.6 Å. For the 1 × 3 × 2 supercell model, a 1 × 1 × 1 grid (i.e., gamma-point) was used because the simulation cell was nearly cubic with dimensions of 9.9 × 9.0 × 9.2 Å. This combination of direct and reciprocal lattices provided nearly isotropic accuracy with respect to the energy calculations [21,22]. Any difference in the accuracy of the energy between small and large simulation cells is not significant because we only compared calculated structures and not the energies between different models. The "Accurate" precision level (which uses the energy cut-off input as 700 eV and doubles the number of grid points in the fast Fourier transforms used to describe the charge density), as implemented in VASP, was used in all cases. Each k-point grid for sampling the Brillouin zone utilized a first-order Gaussian smearing function [23] of width sigma = 0.1 eV (i.e., ISMEAR = 0 and SIGMA = 0.1 in the VASP input file) in each calculation.
Variable-cell energy minimizations were performed for the calculations of bulk α-FeOOH, but only the atomic positions were allowed to relax for the energy minimizations of the (010) α-FeOOH surfaces. Vibrational frequencies were calculated for the energy-minimized bulk and surface slab models, using a finite difference method (i.e., each atom is displaced individually around its equilibrium position by approximately 0.1 Å) and numerical solution of the Hessian matrix. This solution assumes a harmonic oscillator, so any anharmonicity that is present in the actual vibrational modes will not be accounted for. This is especially important for O-H stretching modes with stronger H-bonds. Energy minimizations and frequency calculations were performed with both the SP-GGA and SP-GGA+U methods in most cases (only the extraordinarily long frequency calculations for the larger models were not repeated with SP-GGA+U). In both cases, the GGA functional corresponded to the SP-PBE exchange-correlation functional. Dudarev's rotationally invariant approach to the SP-GGA+U method was used here [24]. The effective on-site Coulomb and exchange interaction parameters for each Fe atom were set to 4 eV and 1 eV, respectively, as recommended by Rollmann et

-termination defined between non-protonated O atoms and the O2-termination defined between the protonated O atoms.
Rakovan et al. [5] have determined that the O2-termination is the stable termination, so this was the termination selected for this study.
Periodic slabs with surface areas of 13.86 Å 2 (1 × 1) and 124.76 Å 2 (3 × 3) were used for the surface model calculations to test the effect of increasing the surface area on the predicted structures. Slab thickness was 8.88 Å in both cases. Vacuum gaps with z-dimension of 10 Å were used to separate the slabs. Solvated models were built by adding one H 2 O molecule per 5-coordinate surface Fe atom with a bond distance of approximately 2.1 Å (2 H 2 O molecules in the 1 × 1 and 18 H 2 O molecules in the 3 × 3 models). The H 2 O molecules were then relaxed using the COMPASS force field [26] while all of the goethite slab atoms were fixed. For the hydrated models, these were filled with H 2 O molecules to approximate a density of 1 g/cm 3 as closely as possible. This gave a total of 8 and 54 H 2 O molecules in the 1 × 1 and 3 × 3 models, respectively. This represents three layers of hydration (L1 -directly bound to the surface Fe atoms, L2 -H-bonding to the L1 layer, and L3 -H-bonding to the L2 layer; see [27] for experimental verification of these layers on TiO 2 and SnO 2 surfaces). The H 2 O molecules were initially positioned to optimize their H-bonding with the surface. An orientation that allowed the H 2 O to act as an H-bond donor to the O2 surface atoms and H-bond acceptor from the surface H 2 O molecules was selected. This initial configuration has not been verified experimentally and will influence the final results of the energy minimizations. Molecular dynamics simulations should be performed to test the reliability of this initial configuration. However, we justify this current decision based on the acidity of the surface groups present. The H + ions in the Fe-OH 2 surface groups should be more acidic than the H + ions on the O2 surface groups [28]; hence the Fe-OH 2 sites should be better H-bond donors and the O2 atoms better H-bond acceptors ( [29] and references therein).

Bulk structure and vibrational frequencies
The first test of the computational methodology was to compare the observed and calculated crystal structures. The results in Table 1 show that both the 1 × 1 × 1 unit cell (Additional file 1) and the 1 × 3 × 2 supercell (Additional file 2) models reproduced the observed lattice parameters and fractional coordinates, within a small percentage error. The current results decrease the already small discrepancy between the calculated values of Rosso and Rustad [7] and the observations of Szytula et al. [13]. For example, the a, b and c lattice parameters in [7] are 9.80, 3.00 and 4.49 Å compared to the 9.95, 3.00 and 4.60 Å computed in this study. Hence, the maximum error is decreased from 3% to 0.5% via use of the pseudopotentials employed in this study. Compared to the unit cell, the supercell methodology provided a more stringent test for predicting the observed crystal structure because there was more flexibility for the ions to relax from their initial positions. Figures 2 and 3 graphically illustrate the comparison between the observed and calculated crystal structures. The agreement was good and did not significantly depend on whether the 1 × 1 × 1 unit cell or the 1 × 3 × 2 supercell was used for the DFT calculation. Differences were on the order of a few percent in each case. The SP-GGA+U method (Additional file 3) gave the same lattice parameters as calculations without the +U correction although the fractional coordinates were slightly different between the two methods (Table 1). These results suggest that the structures predicted with the methodology described above are accurate and independent of the size of the model system for the bulk α-FeOOH mineral.
One significant discrepancy between the observed and calculated structural parameters was the (Fe)O-H bond length. Experimentally, the (Fe)O-H bond length is 0.88 Å, whereas the DFT-calculated (Fe)O-H bond length was approximately 0.99 Å. Given that the IR frequencies and intensities of Fe-OH vibrational modes strongly depend upon O-H bond lengths [30], a difference of 0.1 Å is expected to be significant. Other calculations, such as molecular orbital theory calculations, also predict O-H bond lengths in the range of 0.96 to 1.00 Å in minerals and Fe-hydroxide molecular clusters [31][32][33][34]. The vibrational frequencies predicted by molecular clusters are generally in good agreement with observation for a range of compounds and materials [35]. Furthermore, X-ray diffraction methods are insensitive to the positions of H atoms in crystalline materials and thus their atomic positions must be inferred. Consequently, it is possible that the calculated (Fe)O-H bond length is more accurate than the experimentally estimated (Fe)O-H bond length.
The reproduction of crystal structures can test the ability of a particular method to determine a system's minimum energy position on its potential energy surface. Vibrational frequencies probe the second derivatives of the potential energy surface (i.e., the Hessian matrix) around this minimum and are related to the force constants of bonds. Atomic vibrations can also be used to calculate entropies and other thermodynamic properties of crystals [36]. Consequently, vibrational frequencies are useful observables to calculate and serve to further validate a computational methodology. Table 2 lists the observed and calculated vibrational frequencies for α-FeOOH. In general, there is good correspondence between the calculated and observed vibrational frequencies, especially considering that the DFT-calculated frequencies are determined numerically within the harmonic oscillator approximation. Use of the SP-GGA+U method decreases the calculated O-H stretch-    Table 2 reasonably correspond to experimentally-measured vibrational modes both from IR spectroscopy and inelastic incoherent neutron scattering [38]. Therefore, we can be confident that the calculations are reasonably modeling the potential energy surface around the minimum energy structure.

Surface structure, H-bonding and vibrational frequencies Hydrated surface
After cleaving the (010) α-FeOOH surface, H 2 O molecules were added to each 5-fold coordinated Fe atom (Fig.  4). Unless metal oxides are cleaved under ultra-high vacuum, gaseous H 2 O will adsorb at the surface to satisfy    Fig. 3a). System size effects were minimal for this model system. (C) Shows the 1 × 1 slab structure with SP-GGA+U included, and the bond lengths change by a few hundredths of an Angstrom compared to Fig. 4a Figure 4 (Additional files 4-6).
Regardless of whether the 1 × 1 or 3 × 3 slabs were used in the DFT calculations, the DFT-calculated bond lengths were similar to within a few hundredths of an Angstrom (Fig. 4a versus Fig. 4b). However, changes in the interatomic distances can be on the order of 0.1 Å when this surface was modeled with the SP-GGA+U method compared to SP-GGA (e.g., compare Fig. 4a and 4c). The SP-GGA+U method appears to predict O-H and H-bonding more realistically based on the bulk goethite results, so values in Fig. 4c are likely to be more accurate. Because it was impractical to perform a frequency calculation for the 3 × 3 slab using our currently available computational resources, we focus on the 1 × 1 slab model (Fig. 4a).
The calculated Fe-O(H 2 ) bond length was significantly longer (Fig. 4)  The DFT-calculated IR vibrational frequencies for the hydrated (010) α-FeOOH surface generally agreed well with observation ( Table 3). The agreement is noteworthy because significant error is expected due to the fact that harmonic frequencies are calculated whereas anharmonic frequencies are observed. One exception was the observed vibrational mode near 800 cm -1 , which had its closest corresponding DFT-calculated frequency at 727 or 758 cm -1 for the SP-GGA and SP-GGA+U results, respectively (Table  3). In addition, the O-H stretching frequencies of the model were predicted, in general, to have higher energies than the observed vibrational mode at 3160 cm -1 (although this error is less in the SP-GGA+U calculation, Table 3). Because the frequencies of the O-H vibrational modes are highly dependent upon the H-bonds they form, this difference was probably due to the presence of more than one layer of adsorbed H 2 O on the (010) α-FeOOH surface. The model O-H functional group with the strongest H-bond, for example, had a corresponding DFT-calculated vibrational frequency of 3210 cm -1 -50 cm -1 or less than 2% deviation from experiment. This suggestion is tested below in discussing the solvated α-FeOOH surface.

Solvated surface
In this section, we first compare the solvated (010) α-FeOOH surface to the hydrated (010) α-FeOOH surface for the 1 × 1 slab discussed above. Second, a comparison of the solvated (010) α-FeOOH surfaces generated via energy minimization of the 1 × 1 and 3 × 3 slabs are compared. Third, the vibrational frequencies of the 1 × 1 sol- Frequencies less than 600 cm -1 have not been included, but the output file is available as supplemental information. a = 1 × 1 unit cell b = 1 × 1 with SP-GGA+U c = average of two closely spaced frequencies d. [46] e. [47] f. [48] g. [49] h. [50] i. [51] vated (010) α-FeOOH slab are examined in comparison to experimental vibrational spectra. Last, a comparison of the potential energies of the solvated surfaces with associated and dissociated H 2 O is made.
H-bonding between the L1 and L2 layers significantly affected the Fe-O bond lengths in the models used in this study. One H 2 O molecule acts as a proton donor and acceptor with the surface forming a 1.60 Å H-bond between (Fe)-OH 2 functional groups and a 1.55 Å H-bond to the O atom of the Fe-(OH)-Fe surface linkage (Fig. 5 Fig. 4b). In the 3 × 3 slab calculation, the L2 layer preferred to create a strong H-bonding network within the L2 layer, rather than to the (010) α-FeOOH surface. H-bonds between the L2 layer and the surface did exist, but these were weak with H---O bond distances of 2.43 Å for (Fe)OH 2 ---OH 2 and 2.87 Å for OH 2 ---(OH)Fe 2 (Fig. 4b). Consequently, a H-bond network between (Fe)-OH 2 functional groups and Fe-(OH)-Fe linkages was generated, but steric constraints kept these H-bonds  (Fig. 5). The DFT energy minimizations exaggerated this pattern, so MD simulations should be performed to test how stable this structure is at finite temperatures. However, the underlying structure of the mineral should influence the H-bonding pattern of the L2 structure.
As discussed above for the hydrated (010) α-FeOOH surface, the DFT-calculated frequencies in the O-H stretching region had significantly higher energies than the observed vibrational mode near 3200 cm -1 for α-FeOOH (Table 4; [41]). However, this vibrational mode is very broad and the frequencies are strongly affected by H-bonding. With the addition of H 2 O to form the solvated (010) α-FeOOH surface model, the calculated O-H stretching frequencies decreased significantly and ranged from 2370 to 3680 cm -1 -a similar range to the observed IR vibrational mode (Table 4). Deconvoluting this broad vibrational mode and relating frequencies to individual hydroxyl groups was not practical. However, the agreement between observation and model suggests that the calculations can reasonably predict H-bonding at the surface. Furthermore, the DFT results are consistent with the idea that oxides typically adsorb at least two to three layers of H 2 O from the atmosphere that affect the structure and vibrational spectra of metal-oxides.  (Fig. 5). Theoretically, a surface comprised of Fe 3 OH + Fe 2 O + FeOH 2 is another combination. The relative energies of these configurations were investigated using the SP-GGA+U method for the solvated surface because this method is thought to result in more accurate energies [25]. (Note: this was not done for the hydrated surface because H-bonding to the L2 likely plays a significant role in stabilizing the protonation state of the surface.) It is also worth noting here that the interatomic distances predicted using the SP-GGA+U method are similar to those obtained for the SP-GGA method.
A minimum potential energy configuration was found for the dissociated surface with Fe 3 OH + Fe 2 OH + FeOH, as suggested by Rakovan et al. [5] (Fig. 7; Additional file 10). The H-bonding network between the (010) α-FeOOH surface and water, as well as at the surface, shows that all atoms have formed H-bonds. Many of these H-bonds were fairly strong with H---O distances of 1.51 Å between the Fe 3 OH---HOFe 2 and 1.58 Å between the O atom of the FeOH group and the closest H 2 O molecule. However, the calculated potential energy of this configuration was -213.2039 eV. The Fe 3 OH + Fe 2 O + FeOH 2 surface (Fig. 8) did not result in a stable potential energy minimum because the Fe 3 OH proton was transferred to form an Fe 2 OH group, resulting in a configuration similar to the original associated surface (i.e., Fe 3 O + Fe 2 OH + FeOH 2 - Fig. 9; Additional file 11). The calculated potential energy of this configuration was -214.3466 eV or approximately 75 kJ/mol (where the "mol" indicates a mole of H + on the α-FeOOH (010) surface) lower in energy than the dissociated surface in Fig. 7. Consequently, the current results suggest that the associated surface is likely to be more stable at 25°C unless the entropy of the dissociated surface is greater than 250 J mol -1 K -1 (75 kJ mol -1 /298K) which is large for a deprotonation reaction [42].

Discussion
Although the above results represented energy-minimized structures on a neutral surface, and further testing must be performed under various surface charge states with MD simulations, it is instructive to compare the results to previous studies of the goethite-water interface. In our model, we began and ended with a (010) α-FeOOH surface terminated by Fe-OH 2 functional groups. This is similar to the result of Rustad et al. [8] where Fe-OH 2 functional groups were generated from their classical MD simulations. However, the initial structure in [8] was a Fe-OH-terminated surface and the Fe-OH 2 functional groups were created as solvent H 2 O molecules transferred protons to the surface and OHions H-bonded to the surface.
A similar discrepancy exists between the current study and the (010) α-FeOOH surface model of Rakovan et al. [5]. In [5], a Fe-OH-terminated surface was the final configuration. Rakovan et al. [5] did not use a solvation layer above their (010) α-FeOOH surface and interpreted their XPS data to conclude that 5-fold coordinated surface Fe atoms should be coordinated with OH functional groups. In our model, such a termination was metastable, but the associated H 2 O surface was lower in potential energy. This is based on a neutral surface that should represent α-FeOOH (010) at the point of zero-charge (PZC). However, the water layers in this model do not contain H + or OH -. Consequently, DFT-based MD simulations would be useful to test whether H + transfers occur from the surface to solvent H 2 O molecules, resulting in a Fe-OH-terminated surface charge-balanced by H 3 O + ions or whether a protonated surface and aqueous OHforms as should occur for this surface, which is below its PZC at pH 7. Simply moving H + in these energy minimizations could not test this because the H 2 O layers would need to re-arrange significantly to adjust to this H + -transfer.
In agreement with a study performed by Kerisit et al. [6], who found ordering of 4 to 5 H 2 O layers, the H 2 O molecules in our solvated model were ordered to at least 3 layers of H 2 O molecules. However, the large degree of ordering in our model structures is partly an artifact of the highly ordered initial state and the fact that energy minimizations at 0 K were performed, not MD simulations at finite temperature. Again, DFT-based MD simulations are necessary to examine the configuration space of the α-FeOOH-H 2 O interface.
The studies discussed above did not report Fe-O(H 2 ) bond distances, nor H-bond distances between solvent H 2 O molecules and the α-FeOOH surface, so we cannot compare our DFT-calculated (010) α-FeOOH-H 2 O interface structures to these studies quantitatively. However, Aquino et al. [11] did report O-H bond distances and Hbond distances from their DFT calculations, using cluster models of the (110) α-FeOOH surface. The (110) surface in the study of Aquino et al. [11], however, was terminated by OH functional groups, which complicates comparison. Their O-H bond distances ranged from 0.97 to 1.00 with no systematic variation among the terminal Fe-OH (hydroxo), bridging Fe 2 OH groups (μ-hydroxo), and Fe 3 OH (μ 3 -hydroxo). These values were similar to those predicted in our periodic DFT calculations, but H-bonding played a role in lengthening the Fe 2 OH hydroxyl OH bonds of the hydrated surface compared to terminal Fe-OH 2 groups in the current work. This effect was not observed for the solvated models.