Molecular simulation of the magnetite-water interface

Geochimica et Cosmochimica Acta, Vol. 67, No. 5, pp. 1001–1016, 2003 Copyright © 2003 Elsevier Science Ltd Printed in the USA. All rights reserved 0016-7037/03 $30.00 ⫹ .00

Pergamon

PII S0016-7037(02)00900-6

Molecular simulation of the magnetite-water interface JAMES R. RUSTAD,* ANDREW R. FELMY and ERIC J. BYLASKA William R. Wiley Environmental Molecular Sciences Laboratory, MS K8-96, Pacific Northwest National Laboratory, Richland, WA 99352, USA (Received December 17, 2001; accepted in revised form March 4, 2002)

Abstract—This paper reports molecular dynamics simulations of the magnetite (001)-water interface, both in pure water and in the presence of a 2.3 molal solution of NaClO4. The simulations are carried out using a potential model designed to allow the protonation states of the surface functional groups to evolve dynamically through the molecular dynamics trajectory. The primary structural quantities investigated are the populations of the surface functional groups, the distribution of electrolyte in the solution, and the surface hydrogen bonding relationships. The surface protonation states are dominated by extensive hydrolysis of interfacial ⫺ water molecules, giving rise to a dipolar surface dominated by FeOH⫹ arrangements. Triply 2 -OH2-OH coordinated, more deeply buried, surface sites are inert, probably due to the relative lack of solvent in their vicinity. The electrolyte distribution is oscillatory, arranging preferentially in layers defined by the solvating water molecules. The presence of electrolyte has a negligible effect on the protonation states of the surface functional groups. Steady-state behavior is obtained for the protonation states of the surface functional groups and hydrogen-bonding network. Although the overall structure of the electrolyte distribution is fairly well established, the electrolyte distribution has not fully equilibrated, as evidenced by the asymmetry in the distribution from the top to the bottom of the slab. Copyright © 2003 Elsevier Science Ltd et al., 1999), and controlled hydroxylation studies (Henderson et al., 1998). Previous simulations of oxide-water interfaces by our group (Wasserman et al., 1997) were run on single-processor computers on time scales of tens of picoseconds. These simulations could not make a serious attempt at establishing the extent of equilibration of either the populations of the surface functional groups or of the solvent structure and interfacial hydrogen bonding relationships. The simulations were used mainly to generate plausible hydroxylated surface structures rather than to investigate the solvation of the interface. With recent advances in parallel computing technology, it is now possible to run nanosecond time-scale molecular dynamics simulations with our potential models. We report here the results of the first simulation studies applying this model to the magnetite-water interface. The purpose of these simulations is simply to observe what happens to a hydroxylated magnetite (001) surface when placed in contact with a layer of water and water-sodium perchlorate solution. Particular questions to be addressed include the importance of the solvent in determining the populations of the surface functional groups, the structure of interfacial hydrogen bonding, and the electrolyte distribution. Also of interest are the rates of the various reactions by which the surface structure is established and, particularly, whether equilibration takes place on the time scales of these numerical experiments. We choose magnetite primarily because of the recent very high-quality experimental/modeling study of surface charging (Wesolowski et al., 2000) and also because of the interesting issues surrounding the low isoelectric point of magnetite.

1. INTRODUCTION

Current conceptual understanding of the molecular-level structure of the solvated oxide-water interface is still in a fairly primitive state. The arrangements and reactivities of individual surface functional groups, the role of double-layer structure, the specific adsorption of background electrolytes, and possibly even the usefulness of describing the surface in terms of functional groups are fundamental questions about which there still is meaningful debate (Wesolowski et al., 2000; Hiemstra et al., 1996; Rustad et al., 1996; Felmy and Rustad 1998; Fenter et al., 2000; Rudzinski et al., 2000a; Rudzinksi et al., 2000b; Boily et al., 2001). Because the key player in aqueous interfacial chemistry, the proton, is so difficult to observe, molecular simulations will play an important role in exploring some of these questions. This parallels progress on similar problems associated with proton transport in aqueous solution, for which the details of proton transfer reactions have been worked out almost entirely using molecular simulation (Geissler et al., 2001). Research in this group over the last several years has been focused on the construction and testing of a simple molecular simulation framework for investigating reactive iron oxidewater interfaces (Rustad, 2001, and references therein). The physics of the framework are based on the water model of Stillinger and David (1978) as revised by Halley et al. (1993). The construction and benchmarking phases have used information from ab initio electronic structure calculations in molecular, surface, and bulk crystalline environments (Curtiss et al., 1987; Wang et al., 1998; Rustad et al., 1999; Rustad et al., 2000; Rosso and Rustad, 2001), aqueous solutions (Rustad et al., 1995), bulk crystal structures and energies (Laberty and Navrotsky, 1998), vacuum surface terminations (Thevuthasan

* Author to whom ([email protected]).

correspondence

should

be

2. SURFACE SITES ON MAGNETITE (001)

In the bulk magnetite structure above the Verwey transition, each O2⫺ can be thought of as being coordinated to three octahedral Fe2.5⫹ ions and one tetrahedral Fe3⫹. We terminate

addressed 1001

1002

J. R. Rustad, A. R. Felmy and E. J. Bylaska

the surface along (001) using the stoichiometric autocompensated (公2⫻公2) scheme (Kim et al., 1997) (also known as the “A” termination) and then fill the coordination shells of each of the iron ions with water molecules to octahedral or tetrahedral, according to what they were in the bulk. The “B” termination (Voogt, 1998; Stanka et al., 2000) is avoided because of the nonstoichiometry, the potential structural complexity associated with the oxygen vacancies and also because, as calculated using our potentials, this termination is higher in energy than the “A” termination over a wide range of oxygen fugacities (Rustad et al., 2001). It must be stressed, of course, that our choice of the surface is only a model for what might be potentially much more complex arrangements in natural interfaces. We define the following notation for the surface functional groups: (i) octahedral Fe2.5⫹ sites are labeled Feo, tetrahedral Fe3⫹ sites are labeled Fet; (ii) the functional groups are specified by the number of Fet and Feo bound to the oxide ion, followed by the number of attached protons. Within one unitcell of the “A” terminated (001) surface, after hydration, there twelve sites: four triply coordinated FetFeo2O(H) sites, two triply coordinated Feo3O(H) sites, four singly coordinated FeoOH(2) sites, and two singly coordinated FetOH(2) sites. Also, concerning notation, we use the standard practice of referring, in a generic sense, to positive, neutral, and negative ⫺ surface functional groups as SOH⫹ 2 , SOH, SO , even when these species would not literally exist on the surface. The protonation states of these twelve sites, as determined from zero-temperature structure optimizations with the Fe-O-H interaction potentials used here, were determined previously by Rustad et al. (2001) and explored further by Rustad and Felmy (2001). In this structure, FetFeo2O(H) sites are 25% protonated, Feo3O(H) sites are 100% protonated, FeoOH(2) sites are 75% protonated, and the FetOH(2) sites are 0% protonated. Fig. 1 shows the unit cell obtained by Rustad et al. (2001). The proton affinities and gas-phase acidities of each of the functional groups were determined in Rustad and Felmy (2001) by energy minimization calculations in the absence of solvent. Independent of any ad hoc method for predicting the pKa from the gas-phase proton affinity, Rustad and Felmy (2001) found that surface functional groups having essentially the same proton binding energy may have very different pKas as predicted by the much-used multisite complexation (MUSIC) model (Hiemstra et al., 1996). For example, the proton affinity of the FeoOH site is similar to that of the FetOH site, even though the MUSIC model predicts that the pKas of these sites differ by six orders of magnitude. The overall similarity of the gas-phase reactivity of the different surface sites was due to the collective response of the hydrogen-bonding network developed between surface functional groups. 3. POTENTIAL FUNCTIONS

Potential functions for the magnetite-water system have been described previously (Rustad, 2001 and references therein). The major point to emphasize is that all oxygens, protons, and iron ions are treated on an equivalent basis. In other words, the oxygens in the crystal are the same as the oxygens in the water. The only difference between the Fe3⫹ and Fe2.5⫹ ions is the charge. For this study, we also require parameters for sodium

and perchlorate. The general functional form of our potentials is ⌽ MO共r MO; ␮ ជ O兲 ⫽ A MOe ⫺BMOrMO ⫺

C MO q Mq O ⫹ 6 r MO r MO ⫹

S MO共r兲 ⫽ 1 ⫺

ជ O 䡠 rជ MO兲 q M共 ␮ S MO共r MO兲 3 r MO 1

e

(1)

(2)

EMO共r⫺FMO兲

for an M⫹ cation of charge qM at a distance rMOfrom an oxygen of charge qO (always ⫺2) and dipole moment ␮O. SMO damps the charge dipole interactions at short range. The Na⫹-O2⫺ interaction was determined by fitting to the NaOH·H2O crystal structure (Jacobs and Metzner, 1991). Although the model does not give particularly good results on some of the higher hydrates of NaOH (these are highly delicate structures), it does give the Na-H2O radial distribution, as determined in aqueous solution, in good agreement with the first-principles calculations of White et al. (2000). The position of the first maximum is 2.4 Å, ⬇0.1 Å smaller than the value given in White et al. (2000) but in agreement with the average Na-OH2 bond length found in all the NaOH hydrate crystal structures (for a review of these structures, see Rustad et al., 2002). In particular, the model gives a nonzero value for the first minimum, indicating frequent first-shell water exchange reactions. This exchange is important in ensuring the mixing of the electrolyte in solution and minimizing the dependence of the results on the initial conditions. For this study, we have also produced a model “perchlorate” ion by creating a Cl7⫹ ion. The short range and polarization damping function parameters are chosen to fit, as closely as ⫺ possible, the structures of ClO⫺ 4 , HClO4 , and NaClO4, as well as the gas phase proton affinity and Na⫹ affinity of ClO⫺ 4 as calculated from density functional theory using the local density approximation (LDA). The LDA density functional theory was carried out with the DGauss code using the DZVP basis set (Godbout et al., 1992). Potential parameters for Na-O and Cl-O and comparisons of model results with crystallographic data and density functional calculations are listed in Tables 1, 2 , and 3, respectively. The density functional quantities are reproduced rather poorly, so, at this stage, we offer our perchlorate as only a model of a large monovalent anion with poor complexing ability. This is probably enough for our purposes here. 4. MOLECULAR DYNAMICS PROTOCOL

4.1. Run Parameters We describe our molecular dynamics calculations in the following units (molecular dynamics units, or MDU): length ⫽ Table 1. Potential parameters for Na⫹⫺O2⫺ and Cl7⫹O2⫺. All quantities reported in MDU, as defined in Section 4 in the main text. See Equations (1) and (2).

Na-O Cl-O

AMO

BMO

CMO

EMO

FMO

450.0 85

3.6 3.6

17.0 0.0

3.2 4.0

1.8 1.35

Molecular simulation of the magnetite-water interface

1003

Fig. 1. Proton distribution on the hydroxylated magnetite (001) surface and indication of the types of surface functional groups. The unit-cell is denoted by the solid lines. Large light-atoms are oxygen ions. Small dark atoms are Fe3⫹, larger grey atoms are Fe2.5⫹, and small grey atoms are protons. Two views are shown: (a) perpendicular to the slab, and (b) parallel to the slab. View (b) results from rotating (a) 90° upward.

1Å (1.8897 atomic units, AU), mass ⫽ proton mass ⫽ 1822.7 AU, energy ⫽ e2/Å (0.52918 AU). This fixes the unit of time ⫽ 110.9 AU ⫽ 2.6828⫻10⫺15 s. These MDU are commonly employed in classical molecular dynamics simulations with ionic models. All systems were run in the NVT ensemble at 300

K. We use constant-density calculations because, in simulations of pure water, the water model underestimates the density. Both dipoles and thermostats were treated using the extended Lagrangian method (Nose, 1984 [thermostats]; Saboungi et al., 1988 [dipoles]). The dipole mass was set to 0.5 MDU. To keep

1004

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Table 2. Crystallographic information for NaOH H2O (experimental data in parentheses). Crystal System

orthorhombic

Space Group Temperature a[Å] b[Å] c[Å]

Pbca 298.15 (298.15) 11.04 (11.825) 6.221 (6.213) 6.134 (6.069)

Na-O Bond Lengths [Å]

Na-O1a Na-O1b Na-O1c Na-O2a Na-O2b Na-O2c

experiment

calculated

2.31 2.68 2.88 2.34 2.34 2.45

2.25 2.68 2.78 2.30 2.31 2.45

Table 4. The populations of surface functional groups on magnetite (001) for various simulation conditions. T ⫽ 0 calculations are energy minimizations. The populations for solvated finite temperature simulations are running averages (rounded to the nearest integer value) over 106 timesteps. Average deviations for functional groups are given in parentheses. Populations with no deviations indicated have deviations less than 0.01. Group FetFeo2OH FetFeo2O Feo3OH Feo3O FeoOH2 FeoOH FetOH2 FetOH H2O OH⫺

T ⫽ 300 K T ⫽ 300 K T ⫽ 300 K (2.3 T ⫽ 0 (no solvent) (solvated, I ⫽ 0) m NaClO4) 16 48 32 0 48 16 0 32 — —

16 48 32 0 32 32 16 16 — —

16 48 32 0 61 (0.78) 3 2 (0.90) 30 346 (0.84) 14

16 48 32 0 60 (0.99) 4 3 (1.10) 29 348 (0.83) 12

a, b, and c refer to crystallographically distinct bond distances.

the dipoles “cold,” separate thermostats were used on the dipole and ion degrees of freedom (Sprik and Klein, 1988). The mass of the ion thermostat was set to 50 MDU and to a temperature of 300 K. The mass of the dipole thermostat was set to 10 MDU and to a temperature of 5 K. Because of the well-known momentum conservation problems with the two-thermostat system (Morishita and Nose, 1999), the ion momentum was reset to zero every 1000 time steps. The second order Gear algorithm was used to integrate the equations of motion (for the ions, dipoles, and thermostat) with a time step of 0.1 MDU. Simulations were carried out over 4 ⫻ 106 time steps for both sets of simulations in pure water and in the NaClO4 solution. In preliminary simulations of the solvated surface, it was found that at extended times, small numbers of defects began to appear, involving mainly the raising of an iron ion in an octahedral site to the next (empty) level 2 Å above the surface. Although these defects are certainly of interest, we cannot possibly estimate their equilibrium density from these small, relatively short simulations. Also, as stressed above, our surface is offered as only a model for what is undoubtedly a much more complex natural surface. We therefore leave the study of these defects to future work. Henceforth, we restrict our attention to the idealized surface by connecting each iron atom via a weak spring to a fixed site location. The weak spring (with spring constant 0.15 MDU) has a negligible effect on ion dynamics (characteristic frequency ⬇ 200 cm⫺1) but inhibits the iron ions from moving the 4 Å required to create the observed defects. The site location is defined by the average

position of the iron ion during a 300 K molecular dynamics simulation of the hydroxylated slab with no solvent present. All simulations were run in parallel on a 196-processor cluster of 500-MHz Pentium III processors networked with Giganet cLAN high-speed network. The program was compiled with the GNU project Fortran Compiler (g77) compiler version 0.5.24. The MPI/Pro library, version 1.5 was used for message passing. Typically, 32 to 64 processors were used at any one time in the calculations. 4.2. Initial Conditions The interfacial system was constructed from the optimal monolayer-hydroxylated magnetite (001) surface, as discussed above and shown in Fig. 1. A two-sided, eight-unit-cell slab, 24.5 Å ⫻ 24.5 Å ⫻ ⬇ 15 Å, consisting of 80 Fe3O4 units ⫹ 96 hydroxylating water molecules, served as the initial condition for the solid side of the interface. The initial condition for the solvent side is a layer of ice 11 (360 molecules), optimized with the Halley et al. (1993) water model, scaled slightly in x and y to conform to the periodicity of magnetite (001) and in z to correspond to the density of bulk water (assuming a water layer as defined by the uppermost protons on the optimized vacuumterminated slab). The system measures 24.5 Å ⫻ 24.5 Å ⫻ 33 Å and consists of 1928 atoms. The populations of the various oxide functional groups present at the surface are given in the “T ⫽ 0” column of Table 4. For the simulations with electrolyte present, these same initial conditions were used with 15 NaClO4 molecules inserted randomly into the solution. The den-

Table 3. Density functional theory and model results for ClO4, HClO4, NaClO4 molecules. model

DFT

ClO4 Cl-O HClO4 Cl-O NaClO4

1.50 Å(4) 1.66 Å (1), 1.46 Å(2), 1.47 Å (1)

1.50 Å (4) 1.71 Å (1), 1.46 Å (2), 1.45 Å (1)

Cl-O Na-O E(HClO4)-E(ClO4⫺) E(NaClO4)-E(Na⫹)-E(ClO4)

1.53 Å (2), 1.48 Å (2) 2.27 Å (2) 282 kcal/mol 136 kcal/mol

1.53 Å (2), 1.47 Å (2) 2.18 Å (2) 298.5 kcal/mol 127 kcal/mol

Molecular simulation of the magnetite-water interface

1005

Fig. 2. Potential energy as a function of time over the last 5⫻105 time steps of the simulation involving pure water.

sity was then adjusted to 1.14 g/cm3, approximately the density of a 2.3 molal aqueous solution of NaClO4 at 300 K (Lobo and Quaresma, 1981). 5. RESULTS AND DISCUSSION

5.1. Simulation of the Gas-Phase Hydroxylated Termination at 300 k As a baseline for the solvated system, the first simulation was run in the absence of solvent. Heating the system enables the exchange of protons between the FetOH and FeoOH2 functional groups. The system reverts back to the initial state upon cooling to T ⫽ 0. There is no indication that the minimum-energy state is in error; the activation barriers and energy differences between different protonation states of the FetOH and FeoOH2 functional groups are small enough so that there are many tautomeric states within kT of the minimum-energy configuration at 300 K. This system, with the proton populations in the column labeled “T ⫽ 300 K (no solvent)” in Table 4, is the reference state for comparison against the solvated systems. 5.2. Solvation of the Hydroxylated (001) Surface 5.2.1. Surface protonation In Fig. 2, we show the potential energy of the system ⌽(t)

over the final 5⫻105 time steps of the simulation, and, in Fig. 3, the populations of FetFeo2OH, Feo3OH, FetOH2, and FeoOH2 groups as a function of time over the final 5⫻105 time steps in increments of 100 time steps. Both the potential energy and the populations appear to have reached a steady state. The average protonation states of the surface functional groups in the presence of solvent are shown in Table 4. The major change in surface speciation relative to the gas phase is the large increase in the FeoOH2 population. This increase comes from two sources. Approximately 14 protons come from the deprotonation of FetOH2 sites. Interestingly, 14 additional protons come from the solvating water molecules, resulting in the production of 14 solvated hydroxide ions in solution. The “extra” proton indicated in Table 4 results from the occasional sharing of a proton between octahedral and tetrahedral sites in the configuration FetOH—H—HOFeo. This uptake occurs quite rapidly, mostly within 10 picoseconds (⬇ 35 000 time steps). This “secondary hydroxylation” is completely unexpected. It suggests a somewhat different view of the structure of surface functional groups. Instead of representing the neutral surface functional group as SOH, the molecular-level description pre⫺ dicted by our model is closer to SOH⫹ 2 —OH , or, even better, ⫹ ⫺ SOH2 -H2O-OH , since as we show later, the OH⫺ is fairly well solvated. From a thermodynamic standpoint, these are equivalent ways of viewing a neutral surface functional group.

1006

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Fig. 3. The populations of FetFeo2OH, Feo3OH, FetOH2 and FeoOH2 groups as a function of time over final 5⫻105 time steps in increments of 100 time steps.

It is merely a question of the definition of the surface. Some of the implications for this observation are discussed further in Section 6, after we have considered the hydroxide and electrolyte distribution. Also striking is the nonparticipation of the Feo3OH and FetFeo2OH functional groups in the surface protonation/deprotonation reactions. Despite the fact that the gas phase acidities and proton affinities of these sites are comparable to the singly coordinated sites, there are almost no events involving protonation/deprotonation of the triply coordinated surface functional groups. This is shown clearly in Fig. 3. One reason for the inertness of the triply coordinated groups is that they are relatively isolated from the solvent (although we shall see later that most of the Feo3OH functional groups donate a hydrogen bond to water). This isolation should inhibit deprotonation of these sites. After the ⫽SOH -⬎ SO⫺ ⫹ H⫹ ionization reaction, the deprotonated triply coordinated ⫽SO⫺ functional groups are not well solvated, making them less easily deprotonated. Of course, the inertness of these sites could be purely kinetic. Even if deprotonation of the sites were energetically favorable, the lack of solvent in communication with the triply coordinated sites could inhibit the development of conductive “proton wires” (Geissler et al., 2001) along which proton transport could take place. To address this issue, we have performed a simulation starting from a magnetite slab with no triply coordinated hydroxide ions at the surface. This slab was prepared

from the ground-state gas phase slab by introducing a harmonic potential that constrains the protons to the outermost oxygen layer, which consists entirely of singly coordinated sites. The spring constant (0.042 AU) was adjusted by trial and error to be just strong enough to prevent protonation of the triply coordinated oxide ions at the surface. This system was run for 50 ps of equilibration time before being placed in contact with a layer of 240 water molecules. This same configuration of 240 water molecules was also placed in contact with the ground-state slab having 16 FetFeo2OH and 32 Feo3OH sites. Surface functional group populations are shown in Fig. 4. Initially, there is a rapid protonation of the FetFeo2O sites (so rapid that it is not visible on Fig. 4 and a less rapid protonation of the Feo3O sites. The fast response of the FetFeo2O sites is because of hydrogen bonding to the FetOH2 sites. After the high initial load, the proportion of Feo3OH sites increases at the expense of the FetFeo2OH sites, but only slowly. As shown in Fig. 5, the ground state with 16 FetFeo2OH sites and 32 Feo3OH sites has far lower potential energy throughout the comparative simulations. This indicates that while the protonation states are not in equilibrium, the existence of protonated triply coordinated sites is an essential feature of the surface speciation. The potential energy of the sodium perchlorate system over the last 5⫻105 time steps of the simulation is shown in Fig. 6. There is a small but perceptible drift in the potential energy, which is likely related to lack of equilibration of the electrolyte

Molecular simulation of the magnetite-water interface

1007

Fig. 4. Surface protonation states on simulation started with no protons on triply coordinated surface functional groups. This system is smaller than the production calculations.

distribution as discussed below. Addition of the electrolyte has no significant effect on the surface populations. There is some indication of inner-sphere sodium binding at the surface; one of the sodium ions occupies an “octahedral” surface site (i.e., one that would be occupied by an iron ion if the bulk structure were continued further). It is tempting to consider this structure in light of the recent interest in specific electrolyte binding (Felmy and Rustad, 1998; Wesolowski et al., 2000), but, because there is only one example of this interaction, there is not enough sampling to determine whether or not the sodium is displacing protons from the surface. We leave the study of these surface complexes to future work. 5.2.2. Distribution of electrolyte and hydroxide ions in the solvent The distribution of hydroxide ions within the solvent as a function of time is shown in Figures 7 and 8 for pure water and brine, respectively. In both cases, the hydroxide ions are distributed in layers, 0˜.5 Å in thickness, largely between 8 and 14 Å from the center of the slab (or 1˜.5 to 7.5 Å from the surface of the slab). Although the bands appear to be distributed fairly uniformly throughout the water layer, they are not symmetrically occupied, as one would expect in a fully equilibrated

system. For example, the consistent band at ⫺9.5 Å is only sometimes occupied at ⫹9.5 Å, where there appears to be alternation between a hydroxide ion at 10.5 Å and one at 9.5 Å. The relatively uniform distribution, coupled with the occasional rapid transition from one layer to another, does certainly suggest a self-organized, structured arrangement. The “phase” of the layering is altered with a fairly regular shift of 0˜.5 Å in the presence of the electrolyte. The distribution of Na and ClO⫺ 4 in the electrolyte solution is shown in Figures 9 and 10 , respectively. The phase shift appears to result from the presence of a ClO⫺ 4 ion closely approaching the oxide surface at 8˜ Å from the center, or 1.5 Å from the surface, of the slab. The close approach of the perchlorate ion, at the expense of the hydroxide ion, can be explained by the weaker solvation of perchlorate. The preferential accumulation of weakly solvated species in interfacial regions has been observed for cations (James and Healy, 1972) and hydrophobic materials (Stumm, 1992). In this case, we have the additional driving force of the secondary dissociation of water molecules. An important concept is that the hydroxide ion should be displaced from the surface to some extent by more weakly solvated anions. In this case, for example, the addition of NaClO4 to the solution should have the net effect of releasing

1008

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Fig. 5. Potential energy as a function of time for simulation started with no protons on triply coordinated sites compared with the potential energy of the corresponding system in the ground state. Both systems are smaller than those used for the production calculations.

NaOH into the solution. There is considerable experimental evidence for this effect as summarized in Section 5. Three of the ClO4 layers are “doubly-occupied” by perchlorate ions, clearly indicating that the bands are not just random tracks but represent areas toward which the perchlorate ions are attracted. There is also a clear tendency for hydroxide and perchlorate to coexist. This arises from the ordering of water molecules at the interface; the positive end of the overall water dipole moment of the water molecules will tend to point toward the anion layers. 5.2.3. Solvent hydrogen-bonding The development of the hydrogen-bonding network at the interface is easily analyzed from the molecular dynamics output. In Table 5, we show how each of the surface functional groups behaves in terms of donating and accepting hydrogen bonds. The criteria for hydrogen bonding were defined using cutoff radii. The cutoff radius for an O-H bond was 1.3 Å and the cutoff radius for the H—OH bond was 2.3 Å. For the simulations with no electrolyte added, water molecules accept significant numbers of hydrogen bonds from FeoOH2, FetOH, and, perhaps surprisingly, Feo3OH functional groups. The area between two adjacent Feo3OH sites, (see Fig. 1) is favorable for adsorption of water through hydrogen bonds. Water serves as a

hydrogen bond donor, mainly to the FetOH groups. The water, as a whole, is “undersaturated” in terms of accepting hydrogen bonds, accepting on average 1˜.8 hydrogen bonds per molecule relative to the bulk value of 2.0. In terms of intrafunctional group hydrogen-bonding, FetFeo2OH groups donate significant numbers of hydrogen bonds to FetOH functional groups (with an occasional FetOH2). FetOH2 groups are rare, but donate hydrogen bonds to water molecules and FeoOH when they form. Although most of the hydrogen bonds donated by FeoOH2 groups go to water molecules, significant numbers also go to other FeoOH2 and FetOH groups. FeoOH2 groups are excellent hydrogen bond donors, donating on average 1˜.5 hydrogen bonds per group. The favorable hydrogen-bonding arrangement is probably an important contributor to the low acidity of the FeoOH2 groups. In Fig. 11 we show the time evolution of the number of hydrogen bonds donated by Feo3OH, FetOH, and FeoOH2 to water molecules over the final 106 timesteps of the simulation. From this Figure, it is reasonable to conclude that the hydrogen-bonding network developed between the surface and the solvent has equilibrated on the simulation timescales used here. The major effect of the electrolyte is, as expected, to reduce the amount of water that is hydrogen bonding with the interface. The FeoOH2 groups are affected most strongly, with the

Molecular simulation of the magnetite-water interface

1009

Fig. 6. Potential energy as a function of time over the last 5⫻105 time steps of the simulation involving sodium perchlorate solution.

number of FeoOH2—OH2 interactions decreasing from 37 to 28 as more of the water is tied up solvating the electrolyte ions. The number of Feo3OH—OH2 interactions also is decreased. In this sense, it is perhaps surprising that the presence of the electrolyte does not have a stronger effect on the surface proton populations. It might be expected that the perturbation of the solvent hydrogen-bonding network induced by the electrolyte would destabilize some of the FeoOH2 groups. The role of interfacial hydrogen bonding is of some interest, not only structurally but also because recent empirical methods of estimating surface reactivity take hydrogen bonding into account (Bleam, 1993; Hiemstra et al., 1996). In their application of the revised MUSIC model, Wesolowski et al. (2000) made the assumption that for singly coordinated sites, the sum of the number of protons and the number of accepted hydrogen bonds was equal to 2, and that for triply coordinated sites, the sum was equal to 1. From Table 5 (and populations in Table 4), it is shown that in these simulations deprotonated triply coordinated FetFeo2O sites never accept hydrogen bonds. It is also evident that deprotonated singly coordinated sites accept, on average, more than one hydrogen bond: for the zero-electrolyte case, the number of accepted hydrogen bonds for FeoOH is close to 2, and the number of accepted hydrogen bonds for FetOH groups is 1˜.5. The sensitivity of the revised MUSIC model to hydrogen-bonding results in major changes to the

predicted acidities. For example, the log KHy entries in Table 3 of Wesolowski et al. (2000) for protonation of FetFeo2O, and FeoOH sites would change from 4.9 and 12.6 to 9.2 and 8.3, respectively, if the hydrogen-bonding scheme from these simulations were used. The number of hydrogen bonds donated to deprotonated Feo3O groups could not be evaluated, as these sites never form in these simulations. 6. SUMMARY REMARKS

We have presented here the first molecular dynamics simulation of an oxide-water interface with a general interaction potential that makes no a priori assumptions about the distribution of protons in the interfacial region. The interaction potential has been tested and benchmarked in a variety of simpler systems as summarized in a previous review (Rustad, 2001). Two sets of simulations have been carried out. The first case is magnetite (001) in pure water with no electrolyte added; the second is magnetite (001) in contact with a 2.3 molal solution of sodium perchlorate. Three different facets of the interfacial structure were evaluated: the protonation states of the surface functional groups; the distribution of electrolyte; and the hydrogen bonding relationships that are established between surface functional groups and solvating water molecules.

1010

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Fig. 7. Distribution of hydroxide ions in pure water (a) above and (b) below the slab over the last 106 time steps of the simulation. The z coordinate for the oxygen of the hydroxide ion is plotted every 100 time steps.

Molecular simulation of the magnetite-water interface

Fig. 8. Distribution of hydroxide ions in 2.3 molal NaClO4 solution (a) above and (b) below the slab over the last 106 time steps of the simulation. The z coordinate for the oxygen of the hydroxide ion is plotted every 100 time steps.

1011

1012

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Fig. 9. Distribution of Na⫹, in the 2.3 molal solution (a) above and (b) below the magnetite slab over the last 106 time steps of the simulation. Coordinates are plotted every 100 time steps over the final 106 time steps.

Molecular simulation of the magnetite-water interface

6 Fig. 10. Distribution of ClO⫺ 4 , in the 2.3 molal solution (a) above and (b) below the magnetite slab over the last 10 time steps of the simulation. Coordinates are plotted every 100 time steps over the final 106 time steps.

1013

1014

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Table 5. Hydrogen bond “matrix” for surface and aqueous species. Donor relationships are listed via columns, acceptor relationships via rows. For example, H2O donates 543.6 hydrogen bonds to water molecules, 49 hydrogen bonds to OH⫺, etc; water molecules accept 543.6 hydrogen bonds from other water molecules, and 5.1 hydrogen bonds from OH⫺, etc. Blank entries indicate values less than 0.01. Values represent averages over 106 time steps. no electrolyte H2O OH⫺ FetOH2 FetOH FeoOH2 FeoOH ClO 2.3 molal NaClO4 H2O OH⫺ FetOH2 FetOH FeoOH2 FeoOH ClO

H2O

OH⫺

543.6 49.0

5.1

12.3 2.7 1.7 —

FetFeo2OH

Feo3OH

FetOH2

FetOH

FeoOH2

FeoOH

23.1 0.3

1.2 0.1

8.8 0.1

36.7 4.6

1.1

1.4

—

0.1 1.0 —

0.1 11.2 —

—

—

22.2 22.4 3.1 —

—

H2O

OH⫺

FetFeo2OH

Feo3OH

FetOH2

FetOH

FeoOH2

FeoOH

517.5 36.7

2.8 0.6

0.1

19.8

1.9

8.4

28.0 0.82

0.6

0.12 1.6

0.44

6.3 2.1 1.7 20.9

11.4 1.2

The most remarkable aspect of the surface protonation is the extensive dissociation of solvent water molecules creating the ⫺ SOH⫹ functional groups. In both the brine and in 2 -H2O-OH pure water, the OH⫺ ions are largely confined to rather regular layers 0˜.5 Å wide spaced at approximately 1 Å. The occupation within a single layer was monotonous, with interlayer jumps taking place generally only one or two times during nanosecond simulation. Longer runs or rare-events kinds of techniques would be required to evaluate the hydroxide ion distribution more quantitatively. Another important observation concerning the surface protonation is the enhanced basicity of the triply coordinated FetFeo2OH functional groups. The likely explanation is their isolation from the solvent and, hence, the limited ability of solvent polarization to stabilize the deprotonated FetFeo2O site. The presence of the electrolyte was not a significant factor in determining the populations of the surface functional groups. The electrolyte was organized in layers in a manner similar to that of the hydroxide ions. Although ions often were associated with hydroxide ions in anionic layers, the perchlorate fairly consistently tended to accumulate in the immediate vicinity of the region of free volume between the surface functional groups and the solvating water. This is attributed to the poor hydration of perchlorate relative to hydroxide. It seems likely that the neutralization of the positive surface charge of the ⫽SOH⫹ 2 groups by perchlorate would free at least some of the hydroxide ions from the ⫽SOH2(H2O)OH⫺ groups, effectively releasing sodium hydroxide into solution. Without any additional information, such as might come from simulations with thicker solution layers, the model suggests that addition of a salt of poorly solvated monovalent ions such as perchlorate or pertechnetate should increase the solution pH to an extent that depends on the water-mineral ratio of the system. It is of interest to discuss the implications of the molecular ⫺ level representations of the surface species SOH⫹ 2 — OH versus SOH — H2O for interfacial thermodynamic modeling. First, in our molecular level representation, the addition of electrolyte results in the displacement of hydroxyls from the

0.8

0.3

21.3 19.6 4.0 2.8

surface largely owing to the lower solvation energies of these anions relative to hydroxyl. These displaced hydroxyls would have a tendency to increase the solution pH, or, equivalently, require the addition of acid to maintain the same pH values. This electrolyte effect is precisely what is observed in virtually every adsorption experiment, not only for magnetite (Wesolowski et al., 2000; Blesa et al., 1984) but other iron oxides (Yates and Healy, 1975; Dzombak and Morel, 1990), aluminum oxides (Davis and Kent, 1991, and references therein), and even chromium oxides (Yates and Healy, 1975). The required additions of acid are recorded graphically as increasing positive surface charge at each pH value with increasing electrolyte, but the data are really just acid or base additions. This electrolyte effect even occurs at the PZC for iron and chromium oxides where additions of nitrate and sulfate electrolytes have been shown to alter the pH from the PZC, requiring acid additions to maintain the same pH value (Yates and Healy, 1975). Such effects of electrolyte addition are, of course, well known and have been traditionally described in site binding models (see references above) via the addition of both a proton and an anion to a neutral surface site, e.g., for Cl⫺ SOH ⫹ H⫹ ⫹ Cl⫺ ⫽ SOH2⫹ — Cl⫺ ⫹ 2

(3)

In this picture, the SOH contributes to the surface charge (␴0) and the anion is bound either in the surface plane (0 plane) or in an outer “beta” plane (triple layer model). In such models the addition of electrolyte results in the build-up of a positive surface charge that must be overcome when positively charged metal ions bind to oxide surfaces. In our molecular representation, such a surface-charge build-up does not occur from weak electrolyte addition since the surface charge is effectively zero when averaged over a small volume element near the surface. Such a representation is mathematically most closely related to the basic Stern model originally developed for Hg surfaces but applied to oxides as well (see Westall and Hohl, 1980 for a review). The two conceptualizations are not exactly equivalent since our molecular picture envisions the surface

Molecular simulation of the magnetite-water interface

1015

Fig. 11. Hydrogen bonds donated from Feo3OH, FetOH and FeoOH2 functional groups to water in the 2.3 molal NaClO4 solution as a function of time over the last 5⫻105 time steps of the simulation.

charge being zero in a small volume element near the surface whereas the basic Stern (Hg) model has a zero charge surface plane with ion binding being at the inner Helmoltz plane. Westall and Hohl (1980) showed that such models were perfectly capable of describing the effects of acid/base addition to oxides with the expected lower calculated surface charges owing to the surface bound OH⫺ ions being in the same plane as the counter ions. The implications of these issues will be discussed in a subsequent paper devoted to this topic. As far as water-surface interactions, the major hydrogenbond-donating surface functional groups SOH—OH2 are the FeOH2 and Fe3OH groups having iron in octahedral coordination. About one third of the FeOH groups with tetrahedral iron accept a hydrogen bond from water. The intrasurface hydrogen bond network is well established through FeoOH2-FeoOH2, FeoOH2-FetOH, and FetFeo2OH-FetOH bonding. The hydrogenbonding relationships are affected to a noticeable extent by the presence of the electrolyte as water molecules leave the surface to solvate the electrolyte ions. At least for the surface protonation and hydrogen bonding relationships, it appears that the simulations have equilibrated. Evidence for this includes the absence of drift apparent in the protonation states of the surface functional groups and the hydrogen-bonding analysis. However, it also is clear from simulations deliberately started from nonoptimal proton configurations, that the equilibrated arrangement of surface protons

may only be a metastable equilibrium. As evidenced by the asymmetric distribution of hydroxide and electrolyte ions above and below the mineral slab, the simulations have probably not equilibrated with respect to the hydroxide ion and electrolyte distribution in the solvent. This distribution will be the subject of future work.

Acknowledgments—The authors are grateful to Andrea Currie and Jamie Benward of the Pacific Northwest National Laboratory for editing and manuscript preparation. We also wish to thank Rebecca Sutton and two anonymous reviewers for their helpful suggestions. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Geosciences Division, Contract 18328. The research described in this manuscript was performed at the W.R. Wiley Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy’s Office of Biologic and Environmental Research located at Pacific Northwest National Laboratory. Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06 to 76RL01830. Associate editor: W. H. Casey

REFERENCES Bleam W. F. (1993) On modeling proton affinity at the oxide water interface. J. Colloid Interface Sci. 159, 312–318.

1016

J. R. Rustad, A. R. Felmy and E. J. Bylaska

Boily J. F., Lutzenkirchen J., Balmes O., Beattie J., and Sjoberg S. (2001) Modeling proton binding at the goethite (␣-FeOOH)-water interface. Coll. Surf. A: Physicochem. Eng. Aspects 179, 11–27. Curtiss L. A., Halley J. W., Hautman J., and Rahman A. (1987) Nonadditivity of ab-initio pair potentials for molecular dynamics of multivalent transition metal ions in water. J. Chem. Phys. 86, 2319 – 2327. Felmy A. R. and Rustad J. R. (1998) Molecular statics calculations of proton binding to goethite surfaces: Thermodynamic modeling of the surface charging and protonation of goethite in aqueous solution. Geochim. Cosmochim. Acta. 62, 25–31. Fenter P., Cheng L., Rihs S., Machesky M., Bedzyk M. J., and Sturchio N. C. (2000) Electrical double-layer structure at the rutile-water interface as observed in situ with small-period X-ray standing waves. J. Coll. Interface Sci. 225, 154 –165. Geissler P. L., Dellago C., Chandler D., Hutter J., and Parrinello M. (2001) Autoionization in liquid water. Science 291, 2121–2124. Godbout N., Salahub D. R., Andzelm J., and Wimmer E. (1992) Gaussian-type basis sets for local spin density functional calculations—Part I: Boron through neon, optimization technique and validation. Can. J. Chem. 70, 560 –571. Halley J. W., Rustad J. R., and Rahman A. (1993) A polarizable, dissociating molecular dynamics model for liquid water. J. Chem. Phys. 98, 4110 – 4119. Henderson M. A., Joyce S. A., and Rustad J. R. (1998) Interaction of water with the (1⫻1) and (2⫻1) surfaces of␣-Fe2O3(012). Surf. Sci. 417, 66 – 82. Hiemstra T., Venema P., and VanRiemsdijk W. H. (1996) Intrinsic proton affinity of reactive surface groups of metal (hydr)oxides: The bond valence principle. J. Coll. Interface Sci. 184, 680 – 692. Jacobs H. and Metzner U. (1991) Ungewo¨ hnliche H-Bru¨ ckenbindungen in Natriumhydroxidmonohydrat: Ro¨ ntgen- und Neutronenbeugung an NaOH•H2O bzw NaOD•D2O. Zeitschrift fu¨ r anorganische und allgemeine Chemie 597, 97–106. James R. O. and Healy T. W. (1972) Adsorption of hydrolyzable metal ions at the oxide-water interface III—A thermodynamic model of adsorption. J. Coll. Interface Sci. 40, 65– 81. Kim Y. J., Gao Y., and Chambers S. A. (1997) Selective growth and characterization of pure, epitaxial alpha-Fe2O3(0001) and Fe3O4(001) films by plasma-assisted molecular beam epitaxy. Surf. Sci. 371, 358 –370. Laberty C. and Navrotsky A. (1998) Energetics of stable and metastable low-temperature iron oxides and oxyhydroxides. Geochim. Cosmochim. Acta 62, 2905–2913. Lobo V. M. M., and Quaresma J. L. (1981) Electrolyte solutions: Literature data on thermodynamic and transport properties Vol. II. University of Coimbra, Portugal. p. 767. Morishita T. and Nose S. (1999) Momentum conservation law in the Car-Parrinello method. Phys. Rev. B. 59, 15126 –15132. Nose S. (1984) A unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 81, 511–519. Rosso K. M. and Rustad J. R. (2001) Structures and energies of AlOOH and FeOOH polymorphs from plane wave pseudopotential calculations. Am. Mineral. 86, 312–317. Rudzinski W., Panas G., Charmas R., Kallay N., Preocanin T., and Piasecki W. (2000a) (1922) A combined temperature-calorimetric study of ion adsorption at the hematite-electrolyte interface: I. Model of a homogeneous oxide surface. J. Phys. Chem. B. 104, 11912–1. Rudzinski W., Panas G., Charmas R., Kallay N., Preocanin T., and Piasecki W. (2000b) (1935) A combined temperature-calorimetric study of ion adsorption at the hematite-electrolyte interface: II. Models of a heterogeneous oxide surface. J. Phys. Chem. B. 104, 11923–1.

Rustad J. R., Hay B. P., and Halley J. W. (1995) Molecular dynamics simulation of iron(III) and its hydrolysis products in aqueous solution. J. Chem. Phys. 102, 427– 431. Rustad J. R., Felmy A. R., and Hay B. P. (1996) Molecular statics calculation of proton binding to goethite surfaces: A new approach to estimation of stability constants for multisite surface complexation models. Geochim. Cosmochim. Acta 60, 1563–1576. Rustad J. R., Dixon D. A., Rosso K. M., and Felmy A. R. (1999) Trivalent ion hydrolysis reactions: A linear free energy relationship based on density functional electronic structure calculations. J. Am. Chem. Soc. 121, 3234 –3235. Rustad J. R., Dixon D. A., and Felmy A. R. (2000) Intrinsic acidity of aluminum, chromium (III) and iron (III)␮3-hydroxo functional groups from ab initio electronic structure calculations. Geochim. Cosmochim. Acta 64, 1675–1680. Rustad J. R. Molecular models of surface relaxation, hydroxylation, and surface charging at oxide-water interfaces. In Molecular Modeling Theory: Applications in the Geosciences (eds R. T. Cygan and J. D. Kubick). 2001 Rev. Mineral. Geochem. 2, 169 –197. Rustad J. R., Wasserman E., and Felmy A. R. (2001) The magnetite (001) surface: Insights from molecular dynamics calculations. In Properties of Complex Inorganic Solids, Vol. II. (eds. A. Meike, A. Gonis, P. E. A. Turchi, and K. Rajan). Kluwer Academic/Plenum Publishers, Dordrecht, The Netherlands, pp. 499 –509. Rustad J. R. and Felmy A. R. (2001) Molecular statics calculations of acid/base reactions on magnetite (001). In Solid-Liquid Interface Theory (ed:J. W. Halley), Chap. 9, pp. 113. ACS Symposium Series 789, American Chemical Society, Washington, D.C. Rustad J. R., Felmy A. R., Rosso K. M., and Bylaska E. J. (2002) Ab initio investigation of the structures of NaOH hydrates and their Na ⫹ and OH ⫺ coordination polyhedra. Am. Mineral. (submitted). Saboungi M:L., Rahman A., Halley J. W., and Blander M. (1988) Molecular-dynamics studies of complexing in binary molten-salts with polarizable anions—MAX4. J. Chem. Phys. 88, 5818 –5823. Sprik M. and Klein M. L. (1988) A polarizable model for water using distributed charge sites. J. Chem. Phys. 89, 7556 –7560. Stanka B. I., Hebenstreit W., Diebold U., and Chambers S. A. (2000) Surface reconstruction of Fe3O4(001). Surf. Sci. 448, 49 – 63. Stillinger F. H. and David C. W. (1978) Polarization model for water and its ionic dissociation products. J. Chem. Phys. 69, 1473–1484. Stumm W. (1992) Chemistry of the Solid-Water Interface. John Wiley & Sons, New York. p. 114. Thevuthasan S., Kim Y. J., Yi S. I., Chambers S. A., Morais J., Denecke R., Fadley C. S., Liu P., Kendelewicz T., and Brown Jr. G. E. (1999) Surface structure of MBE-grown␣-Fe2O3(0001) by intermediate-energy X-ray photoelectron diffraction. Surf. Sci. 425, 276 –286. Voogt F. C. (1998) NO2-assisted molecular beam epitaxy of iron oxide films. Ph.D. Thesis, Rijkuniversiteit Gronigent, Groningen, The Netherlands. Wesolowski D. J., Machesky M. L., Palmer D. A., and Anovitz L. M. (2000) . Magnetite surface charge studies to 290 degrees C from in situ pH titrations. Chem. Geol. 167, 193–229. Wang X-G., Weiss W., Shaikhutdinov S. K., Ritter M., Petersen M., Wagner F., Schlogl R., and Scheffler M. (1998) The hematite (␣Fe2O3) (0001) surface: Evidence for domains of distinct chemistry. Phys. Rev. Lett. 81, 1038 –1041. Wasserman E., Rustad J. R., Felmy A. R., Hay B. P., and Halley J. W. (1997) Ewald methods for polarizable surfaces with application to hydroxylation and hydrogen bonding on the (012) and (001) surfaces of␣-Fe2O3. Surf. Sci. 385, 217–239. White J. A., Schwegler E., Galli G., and Gygi F. (2000) The solvation of Na ⫹ in water: First-principles simulations. J. Chem. Phys. 113, 4668 – 4673.