- Research article
- Open Access
Periodic density functional theory calculations of bulk and the (010) surface of goethite
Geochemical Transactions volume 9, Article number: 4 (2008)
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).
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.
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.
Goethite (α-FeOOH) is a common and reactive mineral in the environment . α-FeOOH is the most thermodynamically-stable form of the Fe-oxyhydroxides found in soils, groundwater, and acid mine drainage precipitates . α-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-H2O interfaces. For example, Rosso and Rustad  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.  used a force field that allowed for the dissociation of H2O, in order to perform classical molecular dynamics (MD) simulations of the (110) α-FeOOH-H2O interface. Rakovan et al.  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  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.  employed a polarizable version of the SPC/E water model, along with a α-FeOOH-H2O force field , to describe the structure of electrolyte solutions near the (010) α-FeOOH surface. Most recently, Aquino et al.  performed DFT and Møller-Plesset calculations to investigate the H-bonding of H2O 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) α-FeOOH-water 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. . 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 H2O at the α-FeOOH-H2O interface.
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-H2O interface.
Bulk α-FeOOH and the (010) surface models were built using the Crystal Builder and Surface Builder modules of Cerius2 4.9 , respectively. The α-FeOOH lattice parameters and initial atomic positions were based upon the previously published experimental X-ray diffraction measurements of Szytula et al. . It is imperative to note that the Cerius2 4.9 program uses the Pnma space group for α-FeOOH. In the Pbnm space group the (010) surface, which is commonly used for α-FeOOH , 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. . The 2-D periodic slab generated by this process was stoichiometric, neutral, and symmetric. The bare surface consisted of bridging OH groups (Fe2OH or μ-hydroxyl), bridging oxo groups (Fe3O or μ3-oxo), and 5-fold coordinated Fe atoms. Thus, to hydrate the surface and keep it neutral, H2O 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-OH2 functional groups did not indicate a positively-charged surface of protonated Fe-OH functional groups, as is commonly assumed for metal oxide surfaces. The question of whether or not H2O 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 Fe2OH2 + site), so if the (Fe)-OH2 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., Fe3O + H+ → Fe3OH+). 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 H2O molecule and this could dissociate to form Fe-OH and protonated O1 surface atoms (i.e., FeOH2 + FeO → 2FeOH). However, since Rakovan et al.  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; ). The projector-augmented wave (PAW) method [16, 17] was used in combination with the Perdew, Burke, and Ernzerhof (PBE)  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 d5 electronic configuration), with alternating (010) planes in positive and negative spin directions (Fig. 1; ). This ensured that the model mimicked the anti-ferromagnetic 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  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  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 . 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 al. . See Rollmann et al.  for details on the SP-GGA+U method applied to α-Fe2O3.
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 H2O molecule per 5-coordinate surface Fe atom with a bond distance of approximately 2.1 Å (2 H2O molecules in the 1 × 1 and 18 H2O molecules in the 3 × 3 models). The H2O molecules were then relaxed using the COMPASS force field  while all of the goethite slab atoms were fixed. For the hydrated models, these were filled with H2O molecules to approximate a density of 1 g/cm3 as closely as possible. This gave a total of 8 and 54 H2O 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  for experimental verification of these layers on TiO2 and SnO2 surfaces). The H2O molecules were initially positioned to optimize their H-bonding with the surface. An orientation that allowed the H2O to act as an H-bond donor to the O2 surface atoms and H-bond acceptor from the surface H2O 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-OH2 surface groups should be more acidic than the H+ ions on the O2 surface groups ; hence the Fe-OH2 sites should be better H-bond donors and the O2 atoms better H-bond acceptors ( 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  and the observations of Szytula et al. . For example, the a, b and c lattice parameters in  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 , 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–34]. The vibrational frequencies predicted by molecular clusters are generally in good agreement with observation for a range of compounds and materials . 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 . 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 stretching frequencies by approximately 200 cm-1 compared to the SP-GGA method and brings the calculated values closer to observed frequencies (Table 2). Whether or not this apparent agreement is the result of improved H-bond modeling is not clear because we do not know the percent anharmonicity in these O-H vibrations. However, if we consider the difference between calculated and observed values to be solely due to anharmonicity, the SP-GGA results predict an anharmonicity of approximately 8% whereas the SP-GGA+U method predicts approximately 3% anharmonicity. The latter value is more consistent with previous DFT results comparing observed and predicted vibrational frequency anharmonicities , so we believe the SP-GGA+U method is more accurate.
VASP 4.5 does not calculate the IR intensities of the vibrational modes. Thus, it is difficult to assign calculated frequencies unequivocally to observed vibrational modes. Because there are a large number of calculated vibrational frequencies, correlations between model and observed frequencies may be fortuitously accurate. However, given these caveats, the DFT-calculated displacements in Table 2 reasonably correspond to experimentally-measured vibrational modes both from IR spectroscopy and inelastic incoherent neutron scattering . 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
After cleaving the (010) α-FeOOH surface, H2O molecules were added to each 5-fold coordinated Fe atom (Fig. 4). Unless metal oxides are cleaved under ultra-high vacuum, gaseous H2O will adsorb at the surface to satisfy under-coordinated metal atoms. The layer of coordinated H2O molecules in this study was analogous to the L1 layer observed by Mamontov et al.  for TiO2 and SnO2. We term this state a "hydrated" surface to distinguish it from a bare surface without adsorbed H2O and the "solvated" surface with at least two H2O layers H-bonded to the hydrated surface. The structures of the hydrated surface models are displayed in 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(H2) bond length was significantly longer (Fig. 4) than the Fe-O bond it replaced in the bulk structure (2.42 versus 1.87 Å; Fig. 2b, Fig. 3a or 2.38 versus 1.93 for the SP-GGA+U results, Fig. 3c versus 4c). Furthermore, the surface Fe-O bonds shortened from approximately 2.13 to 2.02 Å to compensate for the added H2O molecules, as would be predicted by Pauling's rules. A small difference was predicted between the O-H bond lengths of the adsorbed H2O molecules and the OH functional groups that occupy the (010) α-FeOOH surface. The former were 0.98 Å and the latter were 1.00 Å due to differences in H-bonding. For example, the adsorbed H2O was only weakly H-bonded to the O atom of the Fe-(OH)-Fe linkage (H---O = 2.25 Å and O-H---O = 140°), while the α-FeOOH OH functional groups formed moderately strong H-bonds (H---O = 1.82 Å and O-H---O = 175°). It should be noted that a shorter H---O bond distance and more linear O-H---O angle both correspond to a stronger H-bond . Although this difference may seem insignificant, it has an influence on how solvating H2O molecules interacted with the surface, as discussed below.
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 H2O 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.
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 solvated (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 H2O is made.
The most important difference between the hydrated and solvated (010) α-FeOOH surfaces was the Fe-O(H2) bond distance (Fig. 5; Additional file 7). The DFT-calculated Fe-O(H2) bond distance changed from 2.42 to 2.04 Å with the addition of two layers of H2O molecules H-bonded to the L1 layer (2.38 to 2.02 Å for SP-GGA+U; see Additional file 8). This result indicates that if bond distances predicted via quantum chemical calculations are used to model the mineral-water interface (e.g. ), solvation beyond the L1 layer should be considered. For example, metal-oxygen bond distances used in surface complexation models such as the MUSIC model  (see Discussion below) will be more accurate if solvent H2O molecules are included beyond the L1 layer. Shortening of the Fe-O(H2) bond also resulted in lengthening of the remaining Fe-O bonds around the surface Fe atoms, in closer agreement with bulk values (e.g. the 1.85 Å bond in Fig. 4a increased to 1.91 Å in Fig. 5).
H-bonding between the L1 and L2 layers significantly affected the Fe-O bond lengths in the models used in this study. One H2O molecule acts as a proton donor and acceptor with the surface forming a 1.60 Å H-bond between (Fe)-OH2 functional groups and a 1.55 Å H-bond to the O atom of the Fe-(OH)-Fe surface linkage (Fig. 5). H-bonds can transfer electron density from the O-H bond which allows the O atom in the adsorbed H2O to form stronger bonds with Fe(III) on the surface.
The 1 × 1 slab was highly constrained by symmetry, as evidenced by the changes observed for the energy-minimized 3 × 3 calculation (Fig. 6; Additional file 9). The Fe-O(H2) bonds were longer (2.29 Å) for the 3 × 3 slab model, in comparison to the 1 × 1 slab model. However, this bond distance was still not as long as that predicted for the hydrated surface calculations (2.44 Å; 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)OH2---OH2 and 2.87 Å for OH2---(OH)Fe2 (Fig. 4b). Consequently, a H-bond network between (Fe)-OH2 functional groups and Fe-(OH)-Fe linkages was generated, but steric constraints kept these H-bonds fairly weak with H---O distances of 2.3 to 2.5 Å and O-H---O angles of approximately 110°. The only relatively strong H-bonds on the surface were between the Fe-(OH)-Fe hydroxyl groups and the underlying O atoms as discussed for the hydrated surface above.
As mentioned above in the Hydrated surface section, the OH functional groups of α-FeOOH were H-bonded to the surface, preventing them from behaving as proton donors to solvent H2O molecules. Because the Fe-OH2 functional groups only acted as proton donors in H-bonds to solvent H2O molecules, a pattern of H2O molecules formed in the L2 layer whereby each H2O molecule formed a donating H-bond to the O of the Fe-(OH)-Fe group and an accepting H-bond to the Fe-OH2 functional group (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; ). However, this vibrational mode is very broad and the frequencies are strongly affected by H-bonding. With the addition of H2O 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 H2O from the atmosphere that affect the structure and vibrational spectra of metal-oxides.
The (010) α-FeOOH surface can form different configurations with respect to H+ positions. The physisorbed H2O surface has been discussed above with an Fe2OH site from the original bulk structure, an FeOH2 terminal group (molecular adsorbed H2O) and no H+ on the Fe3O sites. Dissociated H2O configurations are also possible with H2O adding a H+ to the Fe3O site to form an Fe3OH + Fe2OH + FeOH surface (Fig. 5). Theoretically, a surface comprised of Fe3OH + Fe2O + FeOH2 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 . (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 Fe3OH + Fe2OH + FeOH, as suggested by Rakovan et al.  (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 Fe3OH---HOFe2 and 1.58 Å between the O atom of the FeOH group and the closest H2O molecule. However, the calculated potential energy of this configuration was -213.2039 eV. The Fe3OH + Fe2O + FeOH2 surface (Fig. 8) did not result in a stable potential energy minimum because the Fe3OH proton was transferred to form an Fe2OH group, resulting in a configuration similar to the original associated surface (i.e., Fe3O + Fe2OH + FeOH2 – 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 .
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-OH2 functional groups. This is similar to the result of Rustad et al.  where Fe-OH2 functional groups were generated from their classical MD simulations. However, the initial structure in  was a Fe-OH-terminated surface and the Fe-OH2 functional groups were created as solvent H2O molecules transferred protons to the surface and OH- ions H-bonded to the surface.
A similar discrepancy exists between the current study and the (010) α-FeOOH surface model of Rakovan et al. . In , a Fe-OH-terminated surface was the final configuration. Rakovan et al.  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 H2O 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 H2O molecules, resulting in a Fe-OH-terminated surface charge-balanced by H3O+ ions or whether a protonated surface and aqueous OH- forms 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 H2O layers would need to re-arrange significantly to adjust to this H+-transfer.
In agreement with a study performed by Kerisit et al. , who found ordering of 4 to 5 H2O layers, the H2O molecules in our solvated model were ordered to at least 3 layers of H2O 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-H2O interface.
The studies discussed above did not report Fe-O(H2) bond distances, nor H-bond distances between solvent H2O molecules and the α-FeOOH surface, so we cannot compare our DFT-calculated (010) α-FeOOH-H2O interface structures to these studies quantitatively. However, Aquino et al.  did report O-H bond distances and H-bond distances from their DFT calculations, using cluster models of the (110) α-FeOOH surface. The (110) surface in the study of Aquino et al. , 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 Fe2OH groups (μ-hydroxo), and Fe3OH (μ3-hydroxo). These values were similar to those predicted in our periodic DFT calculations, but H-bonding played a role in lengthening the Fe2OH hydroxyl OH bonds of the hydrated surface compared to terminal Fe-OH2 groups in the current work. This effect was not observed for the solvated models.
The range of H-bond distances, 1.6 to 2.3 Å, calculated by Aquino et al.  overlaps the range predicted for the solvated 1 × 1 slab of this study. However, as noted above, the H-bond distances were significantly lengthened when the model was expanded to a 3 × 3 surface slab. Consequently, strong DFT-calculated H-bonds were probably an artifact of the system size and the neglect of H-bonding to H2O layers above the L2 layer. Qualitatively, this was confirmed by the observation that up to three H-bonds were formed between the α-FeOOH surface and one H2O molecule in the Aquino et al.  study. On the other hand, such extensive H-bonding was not observed in this study. The closest result obtained in this work to the H-bonding arrangement in  is in Fig. 9 where the FeOH2 forms both donor and acceptor H-bonds to an H2O molecule. The third H-bond is to the O atom of the adjacent Fe2OH group, however. Recent work by Nangia et al.  has demonstrated the importance of surface H-bonding networks on influencing the structure of the mineral-water interface, so this factor should be considered in models of mineral surfaces.
The MUlti-SIte Complexation (MUSIC) model developed by Hiemstra et al.  is useful for predicting the PZC, as well as macroscopic adsorption behavior of mineral surfaces. Consequently, this type of thermodynamic model has great value in geochemistry and environmental chemistry. MUSIC requires charges, metal-oxygen bond lengths and the number of H-bonds to surface sites in order to calculate the pKa values of these sites that control surface charging and adsorption. Fitts et al.  used second-harmonic generation spectroscopy to show that surface bond lengths and H-bonding arrangements calculated from periodic DFT and classical MD simulations significantly improve the MUSIC-predicted PZC, in contrast to using bulk bond lengths and assumed H-bonding structures. Thus, the structural predictions reported in this study should be useful as input to MUSIC model calculations on charging and adsorption onto α-FeOOH (010) surfaces. The most reliable values are likely to be those determined on the solvated 3 × 3 slab because this model allows for the greatest relaxation and accounts for surface-H2O H-bond interactions that will be present in bulk potentiometric experiments.
Periodic DFT calculations can reasonably reproduce the structures and vibrational frequencies of bulk α-FeOOH. These methods were also applied to predict the surface structures of the (010) α-FeOOH surface. Model predictions of the surface metal-oxygen bond distances and H-bonding configurations should use a simulation cell comprised of multiple surface unit cells to allow for relaxation beyond what is achievable with single unit cell models. Additional layers of H2O molecules beyond the L1 layer at metal-oxide surfaces are necessary to reproduce what would be found at the bulk mineral-water interface. The α-FeOOH surface was predicted to have a configuration consisting of FeOH2, Fe2OH and Fe3O groups (i.e., an associated water adsorption model).
Brown GE, Henrich VE, Casey WH, et al: Metal oxide surfaces and their interactions with aqueous solutions and microbial organisms. Chemical Reviews. 1999, 99: 77-174. 10.1021/cr980011z.
Majzlan J, Grevel KD, Navrotsky A: Thermodynamics of Fe oxides: Part II. Enthalpies of formation and relative stability of goethite (α-FeOOH), lepidocrocite (γ-FeOOH), and maghemite (γ-Fe2O3). American Mineralogist. 2003, 88: 855-859.
Violante A, Del Gaudio S, Pigna M, et al: Coprecipitation of arsenate with metal oxides. 2. Nature, mineralogy, and reactivity of iron(III) precipitates. Environmental Science & Technology. 2007, 41: 8275-8280. 10.1021/es070382+.
Zhang JS, Stanforth R, Pehkonen SO: Irreversible adsorption of methyl arsenic, arsenate, and phosphate onto goethite in arsenic and phosphate binary systems. Journal of Colloid and Interface Science. 2008, 317: 35-43. 10.1016/j.jcis.2007.09.024.
Rakovan J, Becker U, Hochella MF: Aspects of goethite surface microtopography, structure, chemistry, and reactivity. American Mineralogist. 1999, 84: 884-894.
Kerisit S, Ilton ES, Parker SC: Molecular dynamics simulations of electrolyte solutions at the (100) goethite surface. Journal of Physical Chemistry B. 2006, 110: 20491-20501. 10.1021/jp0636569.
Rosso KM, Rustad JR: Structures and energies of AlOOH and FeOOH polymorphs from plane wave pseudopotential calculation. American Mineralogist. 2001, 86: 312-317.
Rustad JR, Felmy AR, Hay BP: Molecular statics calculations for iron oxide and oxyhydroxide minerals: Toward a flexible model of the reactive mineral-water interface. Geochimica et Cosmochimica Acta. 1996, 60: 1563-1576. 10.1016/0016-7037(96)00035-X.
Shroll RM, Straatsma TP: Molecular dynamics simulations of the goethite-water interface. Molecular Simulation. 2003, 29: 1-11. 10.1080/0892702031000065683.
de Leeuw NH, Parker SC: Molecular-dynamics simulation of MgO surfaces in liquid water using a shell-model potential for water. Physical Review B. 1998, 58: 13901-13908. 10.1103/PhysRevB.58.13901.
Aquino AJA, Tunega D, Haberhauer G, Gerzabek MH, Lischka H: Quantum chemical adsorption studies on the (110) surface of the mineral goethite. Journal of Physical Chemistry C. 2007, 111: 877-885. 10.1021/jp0649192.
Cerius2 Materials Science v. 4.9. 2004, Accelrys Inc. San Diego, CA
Szytula A, Burewicz A, Dmitrijevic Z: Neutron Diffraction Studies of α-FeOOH. Physica Status Solidi. 1968, 26: 429-434. 10.1002/pssb.19680260205.
Sheals J, Sjoberg S, Persson P: Adsorption of glyphosate on goethite: Molecular characterization of surface complexes. Environmental Science & Technology. 2002, 36: 3090-3095. 10.1021/es010295w.
Kendall TA, Hochella MF, Becker U: Computational modeling of azotobactin-goethite/diaspore interactions: applications to molecular force measurements and siderophore-mineral reactivity. Chemical Geology. 2005, 216: 17-35. 10.1016/j.chemgeo.2004.10.007.
Kresse G, Furthmüller J: Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical Review B. 1996, 54: 11169-11186. 10.1103/PhysRevB.54.11169.
Blöchl PE: Projector Augmented-Wave Method. Physical Review B. 1994, 50: 17953-17979. 10.1103/PhysRevB.50.17953.
Kresse G, Joubert D: From ultrasoft pseudopotentials to the projector augmented-wave method. Physical Review B. 1999, 59: 1758-1775. 10.1103/PhysRevB.59.1758.
Perdew JP, Burke K, Ernzerhof M: Generalized gradient approximation made simple. Physical Review Letters. 1996, 77: 3865-3868. 10.1103/PhysRevLett.77.3865.
Coey JMD, Barry A, Brotto J-M, et al: Spin flop in goethite. Journal of Physics: Condensed Matter. 1995, 7: 759-768. 10.1088/0953-8984/7/4/006.
Monkhorst HJ, Pack JD: Special Points for Brillouin-Zone Integrations. Physical Review B. 1976, 13: 5188-5192. 10.1103/PhysRevB.13.5188.
Bandura AV, Sofo J, Kubicki JD: Derivation of force field parameters for SnO2-H2O surface systems from DFT-PW calculations. Journal of Physical Chemistry B. 2006, 110: 8386-8397. 10.1021/jp0523423.
Bandura AV, Sykes D, Shapovalov V, Truong TN, Kubicki JD: Adsorption of water on the TiO2 (rutile)  surface: A DFT study. Journal of Physical Chemistry B. 2004, 108: 7844-7853. 10.1021/jp037141i.
Methfessel M, Paxton AT: High-precision sampling for Brillouin-zone integration in metals. Physical Review B. 1989, 40: 3616-3621. 10.1103/PhysRevB.40.3616.
Dudarev S, Botton GA, Savrasov SY, Humphreys CJ, Sutton AP: Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study. Physical Review B. 1998, 57: 1505-1509. 10.1103/PhysRevB.57.1505.
Rollmann G, Rohrbach A, Entel P, Hafner J: First-principles calculation of the structure and magnetic phases of hematite. Physical Review B. 2004, 69: 165107-10.1103/PhysRevB.69.165107.
Sun H: COMPASS: An ab initio force-field optimized for condensed-phase applications – Overview with details on alkane and benzene compounds. Journal of Physical Chemistry B. 1998, 102: 7338-7364. 10.1021/jp980939v.
Mamontov E, Vlcek L, Wesolowski DJ, et al: Dynamics and structure of hydration water on rutile and cassiterite nanopowders studied by quasielastic neutron scattering and molecular dynamics simulations. Journal of Physical Chemistry C. 2007, 111: 4328-4341. 10.1021/jp067242r.
Venema P, Hiemstra T, Weidler PG, van Riemsdijk WH: Intrinsic proton affinity of reactive surface groups of metal (hydr)oxides: Application to iron (hydr)oxides. Journal of Colloid and Interface Science. 1998, 198: 282-295. 10.1006/jcis.1997.5245.
Onfroy T, Clet G, Houalla M: Quantitative IR characterization of the acidity of various oxide catalysts. Microporous and Mesoporous Materials. 2005, 82: 99-104. 10.1016/j.micromeso.2005.02.020.
Novak A: Hydrogen Bonding in Solids. Structure and Bonding 18: Large Molecules. Edited by: Dunitz JD, Hemmerich P, Holm RH, Ibers JA, Jorgensen CK, Neilands JB, Reinen D, Williams RJP. 1974, Springer, New York, 177-216.
Rustad JR, Rosso KM, Felmy AR: Molecular dynamics investigation of ferrous-ferric electron transfer in a hydrolyzing aqueous solution: Calculation of the pH dependence of the diabatic transfer barrier and the potential of mean force. Journal of Chemical Physics. 2004, 120: 7607-7615. 10.1063/1.1687318.
Bargar JR, Kubicki JD, Reitmeyer RL, Davis JA: ATR-FTIR characterization of inner-sphere and outer-sphere carbonate surface complexes on hematite. Geochimica et Cosmochimica Acta. 2005, 69: 1527-1542. 10.1016/j.gca.2004.08.002.
Larentzos JP, Greathouse JA, Cygan RT: An ab initio and classical molecular dynamics investigation of the structural and vibrational properties of talc and pyrophyllite. Journal of Physical Chemistry C. 2007, 111: 12752-12759. 10.1021/jp072959f.
Paul KW, Kubicki JD, Sparks DL: Sulphate adsorption at the Fe (hydr)oxide-H2O interface: Comparison of cluster and periodic slab DFT predictions. European Journal of Soil Science. 2007, 58: 978-988. 10.1111/j.1365-2389.2007.00936.x.
Kubicki JD: Interpretation of Vibrational Spectra using Molecular Orbital Theory Calculations. Molecular Modeling: Theory: Applications in the Geosciences, Reviews in Mineralogy and Geochemistry. Edited by: Cygan RT, Kubicki JD. 2001, Geochemical Society, 42: 459-483.
McMillan PF: Vibrational spectroscopy in the mineral sciences. Microscopic to Macroscopic: Atomic Environments to Mineral Thermodynamics. Edited by: Kieffer SW, Navrotsky A. 1985, MSA, Washington DC, 14: 9-63.
Scott AP, Radom L: Harmonic vibrational frequencies: An evaluation of Hartree-Fock, Møller-Plesset, quadratic configuration interaction, density functional theory, and semiempirical scale factors. Journal of Physical Chemistry. 1996, 100: 16502-16513. 10.1021/jp960976r.
Hall PG, Clarke NS, Maynard SCP: Inelastic neutron scattering (TFXA) study of hydrogen modes in α-FeOOH (goethite) and γ-FeOOH (lepidocrocite). Journal of Physical Chemistry. 1995, 99: 5666-5673. 10.1021/j100015a059.
Fitts JP, Machesky ML, Wesolowski DJ, Shang X, Kubicki JD, Flynn GW, Heinz TF, Eisenthal KB: Second-harmonic generation and theoretical studies of protonation at the water/α-TiO2 (110) interface. Chemical Physics Letters. 2005, 411: 399-403. 10.1016/j.cplett.2005.03.152.
Hiemstra T, Venema P, Van Riemsdijk WH: Intrinsic proton affinity of reactive surface groups of metal (hydr)oxides: The bond valence principle. Journal of Colloid and Interface Science. 1996, 184: 680-692. 10.1006/jcis.1996.0666.
Cole KE, Paik Y, Reeder RJ, Schoonen MAA, Grey : CP 2H MAS NMR studies of deuterated goethite (α-FeOOD). Journal of Physical Chemistry B. 2004, 108: 6938-6940. 10.1021/jp0486090.
Bunting JW, Stefanidis D: A systematic entropy relationship for the general-base catalysis of the deprotonation of a carbon acid – A quantitative probe of transition-state solvation. Journal of the American Chemical Society. 1990, 112: 779-786. 10.1021/ja00158a043.
Nangia S, Washton NM, Mueller KT, et al: Study of a family of 40 hydroxylated beta-cristobalite surfaces using empirical potential energy functions. Journal of Physical Chemistry C. 2007, 111: 5169-5177. 10.1021/jp0678608.
Szytula A, Burewicz A, Dmitrijevic Z: Neutron Diffraction Studies of α-FeOOH. Physica Status Solidi. 1968, 26: 429-434. 10.1002/pssb.19680260205.
Liese HC: Selected terrestrial minerals and their infrared absorption spectral data (4000-300 cm-1). Infrared and Raman Spectroscopy of Lunar and Terrestrial Minerals. Edited by: Karr C Jr. 1975, Academic Press, New York, 197-229.
Rochester CH, Topham SA: Infrared Study of Surface Hydroxyl-Groups on Goethite. Journal of the Chemical Society – Faraday Transactions I. 1979, 75: 591-602. 10.1039/f19797500591.
Parfitt RL, Russell JD, Farmer VC: Confirmation of Surface-Structures of Goethite (α-FeOOH) and Phosphated Goethite by Infrared Spectroscopy. Journal of the Chemical Society, Faraday I. 1976, 72: 1082-1087. 10.1039/f19767201082.
Russell JD, Parfitt RL, Fraser AR, Farmer VC: Surface-Structures of Gibbsite Goethite and Phosphated Goethite. Nature. 1974, 248: 220-221. 10.1038/248220a0.
Clarke NS: Ph.D Thesis, University of Exeter. 1984
Atkinson RJ, Parfitt RL, Smart RStC: Infrared Study of Phosphate Adsorption on Goethite. Journal of the Chemical Society, Faraday Transactions I. 1974, 70: 1472-1479. 10.1039/f19747001472.
Hartert E, Glemser OZ: Ultrarotspektroskopische Bestimmung der Metall-Sauerstoff-Abstande in Hydroxyden, Basischen Salzen und Salzhydraten. Zeitschrift fur Elektrochemie. 1956, 60: 746-751.
We are grateful to U. Becker for providing the structure of the O2-terminated (010) goethite surface. K. W. Paul appreciates financial support from a University of Delaware fellowship. J. D. Kubicki acknowledges the support of NSF grants CHE-0221934 and CHE-0431328. VASP calculations were performed at the Center for Environmental Kinetics Analysis, a National Science Foundation and Department of Energy Environmental Molecular Sciences Institute at The Pennsylvania State University (Grant No. CHE-0431328). We appreciate advice and guidance from Andrei Bandura and Jorge Sofo regarding VASP utilization, and the constructive criticisms of three anonymous reviewers who greatly increased the quality of this study.
JDK performed most of the DFT calculations and was the primary author of the manuscript. KWP performed DFT calculations on the 3 × 3 × 1 supercell of the hydrated surface and determined how to perform the SP-GGA+U calculations. KWP was secondary author of the manuscript, especially with regard to the introduction and methods sections. DLS supervised and reviewed the manuscript throughout its development. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Kubicki, J.D., Paul, K.W. & Sparks, D.L. Periodic density functional theory calculations of bulk and the (010) surface of goethite. Geochem Trans 9, 4 (2008). https://doi.org/10.1186/1467-4866-9-4
- Surface Complexation Model
- Hydrated Surface
- Classical Molecular Dynamic Simulation
- Pnma Space Group