Development of a Reactive Force Field for Iron− Oxyhydroxide Systems

May 10, 2010 - We adopt a classical force field methodology, ReaxFF, which is able to ..... Journal of Computer Chemistry, Japan 2017 16 (4), 110-111 ...
0 downloads 0 Views 3MB Size
Download PDF
6298

J. Phys. Chem. A 2010, 114, 6298–6307

Development of a Reactive Force Field for Iron-Oxyhydroxide Systems Masoud Aryanpour,*,† Adri C. T. van Duin,§,‡ and James D. Kubicki Department of Mechanical Engineering, Stanford UniVersity, Stanford, California 94305 and Department of Geosciences and Department of Mechanical and Nuclear Engineering, The PennsylVania State UniVersity, UniVersity Park, PennsylVania 16802 ReceiVed: February 11, 2010; ReVised Manuscript ReceiVed: April 14, 2010

We adopt a classical force field methodology, ReaxFF, which is able to reproduce chemical reactions, and train its parameters for the thermodynamics of iron oxides as well as energetics of a few iron redox reactions. Two parametrizations are developed, and their results are compared with quantum calculations or experimental measurements. In addition to training, two test cases are considered: the lattice parameters of a selected set of iron minerals, and the molecular dynamics simulation of a model for R-FeOOH (goethite)-water interaction. Reliability and limitations of the developed force fields in predicting structure and energetics are discussed. Introduction Iron oxide and oxyhydroxide compounds are common on the Earth’s surface, in soils, and in groundwater. They play a vital role in processes that change and shape our planet. For example, in acid mine drainage, Fe hydroxide precipitates are the major sorption substrate for other metals: Cu, Zn, Co, and Mn.1 Removing toxins (e.g., arsenic) and other contaminants in water purification plants can also benefit from low-cost iron-based adsorbates such as goethite.2-4 These compounds appear in a wide range of morphologies and sizes, and hence are rich in chemical reactivity. Goethite (R-FeOOH) as the most abundant iron-oxyhydroxide compound is isomeric with akagane´ite (β-FeOOH) and lepidocrocite (γFeOOH). Lepidocrocite, often accompanied with goethite, is metastable with respect to it and forms from the oxidation of Fe2+. Goethite, on the other hand, can form from both Fe2+ or Fe3+.5 The solubility of iron compounds depends on pH and oxidation state of the iron ion. Ferrous (Fe2+) and ferric (Fe3+) exhibit fundamentally different solubility properties. Whereas Fe2+ is a soluble cation, upon oxidation to ferric ion, it precipitates and forms insoluble oxide compounds.6 However, much more quantitative information is required to answer detailed questions: for example, the thermodynamic and kinetic stability of compositionally similar phases. Navrotsky et al.7 compiled the thermodynamics of various iron oxide and oxyhydroxide minerals. They showed how the relative stability of different phasesswith and without the presence of solution (water)svaries as a function of the particle size. The hydration effect can be attributed to the variability in surface reactivity of different sizes when in contact with solution as compared to the gas phase. Precipitation of Fe-hydroxides from aqueous solutions is highly dependent on the composition of the solution and rates of crystallization. For instance, one can obtain different phases depending on the solution composition8 and different habits and surface areas of the same mineral * To whom correspondence should be addressed. E-mail: masoud@ stanfordalumni.org. † Stanford University. ‡ Department of Mechanical and Nuclear Engineering, The Pennsylvania State University. § Department of Geosciences, The Pennsylvania State University.

with different precipitation rates.9 Therefore, understanding the nucleation and growth of Fe-hydroxides requires a detailed atomic level understanding of how monomers oligomerize, how oligomers form nanoparticles, and how nanoparticles aggregate to form crystals. Characterization and understanding of mineral solid-water interfaces is a complicated problem, currently being tackled by both experiments and theoretical works.10 Focusing on the latter, classical molecular dynamics (CMD) simulations continue to be a major tool for researchers, because the time and length scales under study are normally out of the realm of ab initio molecular dynamics (AIMD). One successful CMD approach in geochemistry is CLAYFF,11 a force field in which all atoms, represented by point charges, are allowed to translate. LennardJones potentials and electrostatic interactions take into account bonding between atoms including metals and oxygen groups. The force field parameters are fitted to empirical or DFT-based structures.11,12 Larentzos et al.13 conducted an extensive AIMD and CMD (by the CLAYFF method) computational study on the structure and vibrational frequencies of the alumino-silicate talc and pyrophyllite, and compared the results with experiments whenever available. They concluded that CLAYFF consistently overestimated O-H bond length by about 0.06 Å. Also, the CLAYFF values for the Si-O-Si angle, OH surface orientation angle, and Al2OH vibration differed from measured values or by AIMD calculations. More work is thus required to capture these structural details in Al oxyhydroxides by CLAYFF. For iron-based minerals and water interactions, Rustad and co-workers have developed a classical force field based on twoand three-body interactions between Fe, O, and H atoms. Their model for a dissociable water is rooted in the work of Stillinger and David.14,15 Adopting a dissociable model for water is crucial, as failure to employ a dissociable water will not capture the Fe3+ and Fe2+ hydrolysis reactions and can lead to a wrong coordination number for Fe2+ and Fe3+.16 The Halley-RustadRahman force field for a dissociable, polarizable water considers not only Coulombic charge-charge interactions but also other electrostatic terms such as charge-dipole and dipole-dipole interactions.17 The dipole moments on O atoms are dynamically varied as a linear function of the local electric field. Furthermore, a three-body localized potential depending on the structure of

10.1021/jp101332k  2010 American Chemical Society Published on Web 05/10/2010

Reactive Force Field for Iron-Oxyhydroxide Systems H2O molecule is utilized to capture the properties of water in accord with experiment. To perform MD simulations of Fe3+ hydrolysis, Rustad et 18 al. introduced Fe-O interactions into the HRR water model. The simulations suggested that a polarizable potential scaled up from small cluster quantum mechanical calculations can do a better job than potential models fitted to, for example, bulk properties. Thus, their model could be improved by being adjusted to ab initio potential energy surfaces for Fe clusters such as [Fe(OH2)6]3+, [Fe(OH)(OH2)5]2+, and [Fe(OH)2(OH2)4]+.18 Rustad et al.19 then extended this force field to treat bulk ferric oxide and oxyhydroxides. The lattice parameters of goethite, lepidocrocite, akagane´ite, and hematite were predicted within 4% of experimental values. However, internal energetics at 0 K determined akagene´ite and lepidocrocite as more stable than the goethite polymorph, probably because of underestimation of Fe-O and H-bond strengths.19 When dealing with mineral-water interfaces, in addition to the details of bulk and surface structure, the presence of edges should not be overlooked. Rustad and Felmy20 performed largescale MD simulations on goethite (110) and (021) slabs immersed in aqueous solution. Significantly different H+ accumulation patterns were obtained on nanoparticles compared to slabs parallel to the crystal morphology. Their calculations emphasized on the inherent limitations with patching models comprised of interacting crystalline surfaces instead of an explicit account of heterogeneity in surface structures at length scales up to at least 10 nm. The interplay of size and structure-activity, as also revealed by modern experimental techniques, remains a grand challenge for theoretical works. X-ray photoelectron spectra and attenuated total reflectance Fourier transform infrared spectra (ATR-FTIR) by Cwiertny et al.21 identified differences in the formation of surface hydroxyl groups on goethite nanoparticles as compared to microparticles. Although the specific saturation uptake of oxalate in nanorod suspensions was less, the surface-area specific rate was about 4 times greater, suggesting a size-dependent surface activity. To surmount such challenges, researchers continue to improve current methodologies further, essentially to enable studying larger model systems over longer simulation times at higher levels of accuracy. As pointed out by Cygan et al,22 a combined AIMD and CMD approach seems to be promising pathway taken in the theoretical-computational field to answer pertinent questions posed at the current frontier of geochemistry mineral sciences. In this paper, we report on an effort to fill the AIMDCMD gap by developing and testing a reactive force field (ReaxFF)23 that has been especially successful in modeling bond formation/breakage in various systems. The goal is to provide a more generalized description of the Fe-water system that properly describes iron oxide and hydroxide bulk structures and can be straightforwardly extended to describe the dynamics and reactivity of Fe-ions in aqueous and bio-organic environments. To this end, we have developed a ReaxFF reactive force field for Fe/O/H interactions. ReaxFF employs a bond-order dependent potential that also contains a polarizable charge model and Coulomb-interactions between all atoms (including atoms sharing a bond). The force field in ReaxFF also includes bond-order dependent three-body (valence angles) and four-body (dihedral) interactions.23 While originally developed for hydrocarbons,23 it has been found that the combination of covalent and ionic contributions allows the ReaxFF formulation to be applied to a wider range of materials, including covalent, metallic, and ionic systems.24-28

J. Phys. Chem. A, Vol. 114, No. 21, 2010 6299 ReaxFF is orders of magnitude faster than quantum-based methods, and attempts are being made to parallelize the code, which might enable nanosecond-scale dynamics for large systems.29 A ReaxFF potential has been recently developed for water and H+ diffusion,30 whereas here we describe an extension of that potential to the Fe/O/H system. This extension involves training the Fe/O/H parameters against a QM-based database, containing structures, energies, and reactions relevant to the Fe2+/Fe3+ water system. Furthermore, to ensure the applicability of this ReaxFF potential to Fe in a wide range of environments, quantum mechanical and experimental data are applied in training the force field to describe the Fe-metal, iron oxides, and iron hydroxide condensed phases. In this paper, we report on the development of two ReaxFF parametrization focusing on major iron oxides and oxyhydroxides. The two sets of force field parameters are denoted by “full” and “oxides”. The force field parameters of the full case are based on training for calculations on single-Fe containing compounds with oxygen, OH, and H2O. Here, these parameters are trained further against the previous set of data as well as the new set on iron oxides and minerals. The new set of data include: (i) two deprotonation reactions of Fe2+ and Fe3+ in solution; (ii) dissociation of Fe-dimer in solution; (iii) heats of formation of hematite, goethite, lepidocrocite, akagane´ite, wu¨stite, and magnetite; and (iv) elastic responses of hematite and goethite. In the second set of parameters (the oxides case), only the data set presented in this work are used in the training procedure (i.e., items i through iv above). The test set for both force fields include: (i) lattice parameters of a few iron oxide minerals and (ii) a molecular dynamics simulation on goethite(010)-water interaction. The goal of this modeling approach is to address oligomerization, nucleation (nanoparticle), and growth in a more accurate way by including major redox reactions as well as bulk thermodynamics. This will allow for more realistic simulation models to investigate surface defects, complex chemistry, and more realistic concentrations (e.g., pH) to be modeled in a cost-effective manner. Force Field Method and Development The initial force field parameters for the Fe-Fe parameters were taken from an earlier force field development project on bulk-iron metal, based on DFT-calculations on antiferromagnetic BCC and FCC. The ReaxFF parameters have not been published yet, but the DFT data came from ref 31. The O/H parameters were taken from the ReaxFF bulk water description.30 The Fe/ Fe and O/H parameters were kept fixed to these initial values, whereas the Fe/O parameters were reoptimized against the quantum mechanical results presented in this manuscript. For electrostatic interactions, the charge-charge interactions in ReaxFF are truncated using a seventh order Taper-function (see Chenoweth et al.,32 Supporting Information, and de Vos Burchart33). This essentially forces the Coulomb energy and its first, second, and third derivative to go to zero at the cutoff radius, which was chosen at 10 Å for this force field. The ReaxFF polarizable charge model extends into the water-phase; it only provides polarizable point charges and no other modes, such as dipole moments. Force field parameters (given in the Supporting Information) can be used either with the standalone ReaxFF program or with the open-source LAMMPS software.34 Deprotonation Reactions of Fe2+ and Fe3+ Solvated Molecules. One of the fundamental reactions controlling the solubility of aqueous Fe species and the reactivity of Feoxyhydroxide surfaces is the deprotonation of Fe-OH 2 groups.

6300

J. Phys. Chem. A, Vol. 114, No. 21, 2010

Aryanpour et al.

Figure 1. Molecular models used in the simulation of deprotonation reactions 1 at left, and 2 at right. The highlighted O-H bonds were scanned to obtain PES for each case. The dashed lines illustrate H-bond network. Atomic coordinates can be found in the Supporting Information, Tables 3 and 4.

Below pH 5, aqueous Fe3+ is dominated by [Fe(OH2)6]3+, but as pH increases the species [Fe(OH)(OH2)5]2+, [Fe(OH)2(OH2)4]+, [Fe(OH)3(OH2)3], and [Fe(OH)4]- form with a minimum solubility for the neutral Fe(OH)3(OH2)3 species, [see ref 35, page 266]. Besides the Rustad-co-workers force field,17-19 other classical force fields have not been capable of modeling this reaction because the H2O molecules have not been dissociable. This type of reaction has been modeled previously, however, using quantum mechanical techniques.36,37 Consequently, we have used the deprotonation of [Fe(OH2)6]3+ and [Fe(OH2)6]2+ as starting points for developing a reactive force field. The ability to model this deprotonation reaction will be a first step toward modeling the oligomerization and nucleation that occurs as pH is raised within acidic, Fe-bearing aqueous solutions. Furthermore, the ability to model both the Fe3+ and Fe2+ states with the same force field is key in simulations of redox reactions in the Fe-O-H system. Two deprotonation reactions 1 and 2 were modeled by coordinating an Fe atom with six water molecules to mimic the environment around an octahedrally solvated ion. A neutral model state was simulated by adding Cl atoms in the second solvation shell: two Cl for Fe2+, and three for Fe3+. This was necessary to create a neutral model for VASP, so the same was done for Gaussian calculations. Furthermore, extra water molecules were added to each system such that each Cl atom is partially coordinated and represents a solvated anion (Figure 1).

Fe(H2O)6Cl2 · 5(H2O) f Fe(H2O)5(OH)Cl2 · [H+(H2O)5] (1) Fe(H2O)6Cl3 · 7(H2O) f Fe(H2O)5(OH)Cl3 · [H+(H2O)7] (2) In these models the minimum number of H2O molecules have been added to hydrate Fe and Cl atoms. Extra H2O molecules are needed to obtain more accurate thermodynamic values. The potential energy surfaces (PES) corresponding to these deprotonation reactions were obtained by scanning of an O-H bond belonging to the proton-donor complexes [Fe(H2O)6]2+ and

[Fe(H2O)6]3+ (both in the high-spin state as is observed for aqueous ferrous and ferric ions). Computational details were the same as for the deprotonation reactions except that in VASP one of the high-spin Fe3+ ions was a R and the other β. This was not possible in G03, so both Fe3+ ions were in the an R state, which creates a significant energy difference between the VASP and G03 calculations. ReaxFF does not explicitly treat the spin orientation of the ions, so this effect must be implicitly included in the parametrization of the force field. The potential energy surfaces (PES) corresponding to these deprotonation reactions were obtained by scanning of an O-H bond belonging to the proton-donor complexes [Fe(H2O)6]2+ and [Fe(H2O)6]3+. The rest of the bond distances and angles were relaxed to lower each sample energy point on the PES. We used the Gaussian03 (G03)38 suite of programs as well as the Vienna Ab-initio Simulation Package (VASP),39-42 to perform quantum calculations. The B3LYP level of theory was applied to the G03 calculations using basis set and pseudopotential CEP-121G for Fe, and 6-31G(d,p) for O, H, and Cl. No symmetry restriction was enforced on the structure. The results of the first round of G03 calculations were used to train ReaxFF, and the force field faced difficulties in following the quantum energy curve. Upon deprotonation reactions, the H+-acceptor H2O in solution loses one of its H atoms to other solvated water molecules. This H-hopping was avoided in the force field training through constraining two OH bonds on the water molecule, whose O atom gains another H atom. Therefore, we recalculated the energy profiles in G03 by fixing two OH bonds on the proton-acceptor H2O, while all other coordinates were optimized at each step of energy scan. That restriction on two OH bonds prevented H-hopping between the acceptor O atom and nearby H2O molecules when receiving another H from the donor H2Oads. In VASP calculations, the generalized gradient approximation PBE (GGA-PBE) pseudopotentials,43,44 with cutoff energy 500 eV (36.75 Ry), were used for all the elements. That implied an increase by 25% with respect to the default cutoff value (400 eV in PBE pseudopotentials). The whole molecules were inserted in boxes having dimensions 12 × 12 × 12, and 12 × 12 × 15 Å, respectively. A K-point sampling scheme 1 × 1 ×

Reactive Force Field for Iron-Oxyhydroxide Systems

J. Phys. Chem. A, Vol. 114, No. 21, 2010 6301

TABLE 1: Comparison of Model Structures Obtained by Gaussian and ReaxFF (Cases Full and Oxides) in Reactions 1, 2, and 3 at the Lowest Energy Pointa Fe2+, reaction 1 ReaxFF (full) ReaxFF (oxides) Fe-Ob O-Hb ∠FeOHb

2.16 1.03 115.3

2.14 1.03 110.2

Fe3+, reaction 2 G03 2.18 0.98 109.5

ReaxFF (full) ReaxFF (oxides) H2Oads 2.06 1.08 120.2

2.05 1.08 117.7

Fe-dimer, reaction 3 G03 2.01 1.01 113.9

ReaxFF (full) ReaxFF (oxides)

G03

2.08 0.99 126.9

2.06 0.98 126.4

2.09 0.98 126.3

2.14 127.2 105.6 0.0

2.06 125.9 107.9 0.0

1.99 126.4 107.0 0.0

OHads Fe-Oc ∠FeOHc ∠FeOFec dihedral FeOFeO a

Distances are in Ångstroms, and angles have units of degrees. b Averaged over Fe-H2O bonds. c Averaged over Fe-OH bonds.

1 was applied to isolate the system from its own images in other periodic cells. The VASP models were constructed to mimic those used in Gaussian calculations and to compare the results. Averaged bonds and angles of the model molecules at their lowest energy point obtained by Gaussian and ReaxFF are compared in Table 1. For the Fe2+ model, the average of the Fe-O bonds by the two ReaxFF force fields is slightly (0.02-0.04 Å) shorter than by G03, and for the Fe3+ model the Fe-O bond on average is stretched (0.04-0.05 Å). On the other hand, the force fields predict the average of the O-H bonds to be slightly longer (0.05-0.07 Å) in both Fe2+ and Fe3+. The Fe-O-H bond angles by these force fields are wider by about 1-6°. Both force fields give similar structures at the minimum point and along the reaction paths. The bond distances by ReaxFF (full) and ReaxFF (oxides) are also more or less similar, but the latter force field predicts FeOH angles closer to quantum calculations. Figure 2 illustrates the energy curves versus the O-H distance, which is elongated in the deprotonation process, from calculations by G03, VASP, and training ReaxFF. An initial minimum is observed in both reactions, corresponding to the inner minimum (O-H ≈ 1.0 Å), the ground state configuration of the molecules. Afterward, the energy increases, as mapped by stretching the O-H bond, either monotonically for Fe2+ or asymptotically for Fe3+. For the latter compound, G03 evaluates a deprotonation activation of about 45 kJ/mol, whereas VASP estimates this value to be about 35 kJ/mol. ReaxFF gives results (55-65 kJ/mol) closer to 45 kJ/mol because of being trained against the G03 set of calculations. The oxides force field is more accurate in estimating the deprotonation energy. Overall, one can conclude that although there is definitely room for improvement in ReaxFF, the general trends and values of energetics are satisfactorily reproduced. It must be emphasized here that the force field parameters are adjusted to match not only these two reactions, but also a set of various data points among which the above two reactions constitute only a part of the training set.31,45 In other words, the trained force field should be judged more reasonably by considering its performance over a complete data set, rather than by focusing on an isolated reaction profile. Note the differences between G03 and VASP results, which some would consider accurate. Dissociation of Fe-Dimer in Solution. The interaction of two Fe cations in the solution is another important aspect that must be modeled accurately enough in any reliable force field for iron oligomers and nanoparticles. Assuming two solvated Fe cations, they may make bonds if initially isolated, or dissociate from each other if initially forming a dimer. Simulation of the forward dissociation reaction 3, when

Figure 2. Deprotonation energies of corresponding to two reactions: 1 and 2.

calculating the energy curve, is easier and more straightforward than its reverse bonding reaction, because the reaction coordinate can be defined more clearly in the former case.

6302

J. Phys. Chem. A, Vol. 114, No. 21, 2010

Aryanpour et al.

[Fe2(OH)2(H2O)8]4+ · 2H2O f 2[Fe(OH)(H2O)5]2+

(3) One of the Fe-OH bond distances was chosen as the reaction coordinate. G03 calculations reveal that larger dissociation energies were obtained had the Fe-Fe distance been stretched instead of Fe-OH in the scan of PES. This is consistent with a lower energy barrier to breaking (and forming) one bond at a time rather than both simultaneously. Computational details were the same as for the deprotonation reactions. The unit cell in the VASP calculations was a cube of dimensions 13 × 13 × 13 Å. The last column of Table 1 compares bonds and angles of this molecule at the lowest energy point (in Figure 3) as obtained by Gaussian, reported in parentheses, and by ReaxFF. Good agreement is seen for Fe-O and O-H bond distances, as well as for bond angle FeOH averaged over all H2O ligands: less than 0.03 Å difference in the bond lengths, and less than 1° difference in the bond angle. ReaxFF predicts Fe-O bonds for Fe-OH ligands that are on average 0.02-0.07 Å extended with respect to the G03 values. The Fe-O bond distances for Fe2+ and Fe3+ are compatible with experimental estimations,46 2.1 and 2.02 Å, respectively. Bond angles FeOH and FeOFe differ at most by about 1-2°. The two force fields give almost a planar structure (dihedral ∼0) for the pair of Fe atoms and OH groups, in agreement with quantum calculations. Thus, both the force fields prediction of the structure for this molecule agree well with the quantum results. As Figure 3b shows, the molecule reaches its ground state around the point where Fe and O atoms are approximately 2.0 Å apart. Further stretch of the Fe-O bond raises the system energy up to 40 kJ/mol at about 2.8 Å and then starts to level off. ReaxFF (full) gives a lower activation of about 30 kJ/mol, closer to VASP results. ReaxFF (oxides) performs better in following Gaussian profile and predicts a similar (45 kJ/mol) activation energy. The observed fluctuations in the energy curves might be attributed to the rotation of end H2Oads, which can introduce ripples in energy depending on their orientation toward their adjacent H-groups. The resulting points by VASP have been obtained by relaxing structures that start from their previous configuration but shifted along the reaction coordinate. It might be possible to smooth off the local zig-zagging behavior in the VASP energy curve by performing nudge elastic band calculations. However, such further work seems unnecessary at this point, first, because these energy fluctuations are small compared to the changes in energy as the reaction proceeds. Second, the ReaxFF results here would not be affected, because they have reproduced the G03 quantum calculations reasonably accurately. The maximum discrepancy between these two quantum methods is ≈20 kJ/mol (at about 2.8 Å and beyond), and the corresponding rms value of fluctuations in the VASP curve is close to 4 kJ/mol. In summary, the two force fields for the dissociation reaction 3 have been able to faithfully regenerate the G03 quantum results in a quantitative sense over the entire PES scan, ranging from the bonding state to the dissociation state. Nevertheless, the oxides force field better reproduces the quantum-generated values. Thermodynamics of Phases and Bulk Constants. The significance of thermodynamics cannot be overemphasized when studying the stability (or metastability) of different phases.47 To simulate nucleation and growth of particles reliably, any force field must be able to predict phase stabilities at least in a relative sense. Here we account for the thermodynamic stability of goethite, lepidocrocite, akagene´ite, hematite, wu¨stite, and

Figure 3. PES scan of dissociation reaction 3: (a) molecular model where scan bond Fe-O is highlighted. Atomic coordinates can be found in Supporting Information Table 5. (b) Energy curves. ReaxFF parameters were trained against the G03 calculations.

magnetite. These are common phases of Fe-oxyhydroxides and oxides in the crust and mantle (Fe-wu¨stite). Furthermore, by using the Fe2+ mineral wu¨stite and mixed valence state magnetite (Fe3O4), we attempt to give ReaxFF the ability to model Fe redox chemistry. Lattice parameters for hematite,48 goethite,49 lepidocrocite,50 wu¨stite,51 akagane´ite,52 and magnetite53 models were obtained from experimental measurements. For a few of these minerals, thermodynamical data, bulk modulus, and its pressure derivative were collected from refs 7 and 54-59. Other than chemical composition and crystalline structure, the stability of these compounds is also affected by environmental factors, such as pH, pO2, temperature, and pressure [see ref 8, page 193], as well as hydration.7 Details of ReaxFF training set up are explained below. Heats of formation for all systems were imposed as energy cost function expressions with respect to the energy level of their constituent elements: Fe(bcc), O2, and H2 in gas. The elastic

Reactive Force Field for Iron-Oxyhydroxide Systems

J. Phys. Chem. A, Vol. 114, No. 21, 2010 6303

Figure 4. Training the force field against experimental heats of formation. Experimental values are from refs 7 and 54-59.

responses of hematite and goethite structures were included in training and were modeled by five intermediate expansion and five compression steps of their volume. Maximum compression/ expansion reached 10% of the unit cell. Training of the force field against volumetric variations was performed for hematite and goethite only. The training results for thermodynamics are given in Figure 4. For ReaxFF (full) one can see an appreciable discrepancy between the force field and empirical values in heats of formation for hematite and magnetite (∼125 kJ/mol). However, all other values by ReaxFF (full) are much closer to their corresponding experimental measurements. In particular, the absolute difference for goethite, lepidocrocite, and akagane´ite is less than 16 kJ/mol. For ReaxFF (oxides), the hematite energy is 40 kJ/mol higher, whereas the magnetite energy is 7 kJ/mol lower than the experimental levels. Energy levels for goethite, akagane´ite, and lepidocrocite are on average 80 kJ/mol less stable, whereas for wu¨stite they are 120 kJ/mol more stable than measurements. Most importantly, however, the predictions of relative stability of phases (represented by the dashed and marked line) are correct, as they closely follow the same trend of the solid line. Therefore, the relative stability of iron-oxide phases is preserved by both force fields, at least in a qualitative way. Figure 5 demonstrates the results of force field training for the elastic modulus of hematite and goethite. ReaxFF output energy curves are in good agreement with calculations based on experimentally measured quantities, showing an rms deviation of only 1-2 kJ/mol for hematite and 2-4 kJ/mol for goethite. No major difference between ReaxFF (full) and ReaxFF (oxides) is observed for this set of training data. As seen, the force fields reproduce the elastic response of hematite accurately from 10% compression up to 10% expansion. Elastic response of goethite is reproduced most accurately within about (7% volumetric change with respect to the ground state volume. Thus, elastic results generated, at least for hematite and goethite, by the developed force field will be reliable for many practical purposes. Testing Results As reliability tests on the developed force fields, we report on two sets of calculations. Lattice parameters of iron minerals

Figure 5. Training the force field against experimental elastic moduli for hematite and goethite. Experimental values are from refs 54 and 55.

and oxides, targeted in this paper, are predicted by ReaxFF, and compared to experimental measurements. The second test case consists of MD simulations performed by ReaxFF and compared to similar calculations by ab initio MD implemented in VASP. Lattice Parameters of Selected Iron Minerals. Using the two trained force fields, the lattice parameters for the selected iron minerals were predicted. 4000 steps of MD calculations (each by 0.25 fs) were performed in an NPT ensemble (T ) 5 K), which allowed the cell parameters to vary in each case. The use of low-temperature NPT calculations instead of direct relaxations of the whole computational cell is because of the involved complexities in the current version of the code. The simultaneous variation of cell parameters and optimization of ions position has not yet been fully implemented. During MD steps, the cell volume changes according to an external pressure,

6304

J. Phys. Chem. A, Vol. 114, No. 21, 2010

Aryanpour et al.

TABLE 2: Lattice Parameters Estimated by ReaxFF (full) and Compared to Experiments mineral hematite goethite lepidocrocite akagane´ite wu¨stite magnetite

supercell 3 4 3 3 6 3

× × × × × ×

3 4 1 9 6 3

× × × × × ×

3 4 3 3 6 3

atoms

a (∆a%)

b (∆b%)

c (∆c%)

V (∆V%)

240 432 144 800 500 448

4.867 (0.4) 4.494 (2.3) 3.038 (1.4) 10.679 (0.9) 5.188 (6.1) 8.032 (2.6)

4.852 (0.2) 9.902 (0.5) 14.692 (17.5) 3.020 (0.4) 2.953 (2.7) 8.032 (2.6)

13.041 (1.7) 2.963 (1.8) 3.956 (2.3) 10.40 (1.1) 2.954 (3.6) 8.032 (2.6)

266.69 (2.3) 131.86 (4.5) 176.72 (18.6) 335.34 (0.6) 37.23 (7.5) 518.06 (7.6)

and a low temperature assures that most of the kinetic energy is released from the system, therefore its effect on the cell volume is negligible. The results of cell optimizations using ReaxFF (full) are summarized in Table 2. The results by ReaxFF (oxides) are similar, and thus not presented here. Cell parameters of hematite have been calculated within 1% of experiment, with a 2.3% error in volume. Those of goethite are also close to experimental measurements. Lepidocrocite, on the other hand, shows large deviations from measurements (approximately 18%), which is a deviation that stems from an unrealistic expansion in the b-direction of the crystal. The complexity in this case evidently results from H atoms’ crucial role in the interlayer bonding of lepidocrocite. The problem is in accurate description of the Fe(OH)-OH bonds. Currently, the ReaxFF force field is parametrized only for O-HOH bonds. Because layers of the lepidocrocite structure are held together only by H-bonding (as compared with the covalent bond network of goethite where H-bonding is secondary), inaccuracies of the H-bond energy term are exposed more clearly. The polarization of the O atoms in the Fe-(OH)-O H-bonding must be greater than is present in H-(OH)-O bonds, such as those found in water that were used to parametrize the current version of ReaxFF. Consequently, the H-bond strength (i.e., potential energy vs distance) is weaker than it should be in order to hold the lepidocrocite layers closer together as is observed. Starting from an initial structure by program Cerius2,50 we performed a cell and ionic relaxation of lepidocrocite crystal using the VASP package. The intention was to detect major changes in structure and to provide further insights made possible by ab initio computations. We found that the variations in the cell parameters were less than 2% in any direction with respect to the published measurements. In the optimized position of H atoms, the crystal layers in the b-direction are connected through interlayer H-bonds between OH groups. There were convergence difficulties in DFT computations of electronic structure that might be attributed to the importance of such H-bonds in the crystalline structure of lepidocrocite. Therefore, if one wants to improve the ReaxFF force field to account for such delicate H-bonding, new energy terms and parameters have to be defined and incorporated in the current ReaxFF formulation. As another test, we compare relative stability of hematite with respect to bixbyite. Although we have no experimental heat of formation data for Fe2O3 in the bixbyite structure, it is known that hematite should be thermodynamically stable relative to bixbyite under ambient pressures. The ReaxFF (full) and (oxides) parametrizations result in bixbyite energy-minimized structures that are -169 and -253 kJ/mol lower in energy than hematite, however. Thus, ReaxFF overestimates bonding energies for the lower density bixbyite structure. The test results for akagane´ite, wu¨stite, and magnetite are again close to experiment as in the case of hematite and goethite. Mean relative error in lattice parameters of akagane´ite is about 1% with respect to empirical measurments. For wu¨stite and magnetite the average error amounts to approximately 4 and

3%, respectively. The relative errors in cell volume, indicating the collective errors in lattice parameters, is larger. However, such errors will reasonably (and hopefully) decrease as one employs a larger cell to model these crystals. The point of emphasis is that the developed force field achieves this level of accuracy not for any single mineral, but for all the selected set of minerals simultaneously. This is a step toward having a more transferable force field, a capability that allows one to apply one force field to study many systems and compare the results based on a unified approach. MD Simulation of Goethite(010)-Water Interface. The interaction of the goethite surface (010) and water was modeled by solvating a slab of four layers in a unit cell 13.7 × 9 × 19 Å. A total of 36 goethite units (FeOOH) were set up in the solid phase, and the solution consisted of 46 H2O molecules, making a system with 282 atoms (Figure 6). The initial structure was based on a model designed by Kubicki et al.60 to study the goethit-solution interface. Ab initio MD calculations were performed using VASP. Initial magnetic moments of Fe atoms were assigned as alternate positive and negative values ((5) on the crystalline layers along the Miller direction [010]. This way, an antiferromagnetic goethite slab was obtained. The onsite Coulomb interaction, the L(S)DA+U method, was applied to the Fe atoms in the sense of Dudarev et al.61 as implemented in VASP. For AIMD calculations, about 6000 steps (each at 0.5 fs) were performed using the soft potentials GGA-PBE (cutoff 400 eV). The energy profile for the last 2000 steps (1000 fs) is plotted in Figure 7 and shows equilibrium or near-equilibrium conditions.

Figure 6. Model structure of water on goethite (010) used in MD ab initio simulations by VASP and in classical simulations using the reactive force field.

Reactive Force Field for Iron-Oxyhydroxide Systems

J. Phys. Chem. A, Vol. 114, No. 21, 2010 6305

Figure 7. Energy fluctuation profiles of the last 1000 fs in the MD simulations of water on goethite (010). The zero value on each plot corresponds to the average of the system energy.

Using similar initial coordinates, classical MD (CMD) simulations were performed using the ReaxFF methodology. After an initial energy minimization, a total of 40 000 MD steps each at 0.25 fs were executed. The energy profile in Figure 7 demonstrates thermal equilibrium within 34 (40) kJ/mol rms fluctuations in the last 1000 fs for ReaxFF (full) and ReaxFF (oxides), respectively. Ab-initio energy within the same last time period shows larger fluctuations about 60 kJ/mol; perhaps because the total number of simulation steps in ab initio MD (AIMD) was far less than in classical MD, 3 ps versus 10 ps, and time step was twice as large in the AIMD (0.5 vs 0.25 fs). Figure 8a shows the radial distribution function (RDF) for all ions in the liquid phase. The presented results conform well with previous modeling and experimental studies on water structure. For example, the ratio of the second to the first peak in OH RDF by VASP/ReaxFF(full)/ReaxFF(oxides) are 0.14/ 0.15/0.14 (ca.), which are close to the value of 0.10 presented in a recent article on the properties of water by Paesani and Voth [ref 62 and references therein]. For ReaxFF we see a sharper third maximum (ratios 0.11 and 0.08 at r ) 2.1 Å) with respect to the VASP ratio ) 0.06 around r ) 2.3 Å. The ratio of the fourth maximum to the first compares similarly well. The ReaxFF (full) ratio of the fourth peak to the first is 0.15, that of ReaxFF (oxides) is 0.12, which is closer to the VASP ratio of 0.11. The RDF by the ReaxFF (oxides) force field performs better than the ReaxFF (full) force field in predicting the second and the third peak. Figure 8b illustrates the structure of O atoms in the solution. The first peak in OO bonds by VASP occurs at 2.8 Å, and ReaxFF predicts OO distances at 3.0 Å. After the first peak, the main features of the O network by ReaxFF match those in the O structure by VASP. From panels a and b of Figure 8 one can conclude that the force fields developed in this paper are

Figure 8. Structural analysis of water on goethite (010) MD simulations by VASP and ReaxFF, (a) radial distribution function between all atom pairs in solution, (b) radial distribution function between OO pairs in solution.

able to predict fairly accurately the structure of the solution phase in agreement with the DFT calculations. Figure 9a illustrates the pair distribution function for interfacial O and H atoms on the goethite surface. The three curves show their peaks close to each other: first at 1.0, second at 1.6-1.7, third at 2.3-2.5, and fourth at 3.0 Å. The ratio of the second, the third, and the fourth peaks to the first one is as follows: for VASP 0.05/0.06/0.17, for ReaxFF (full) 0.07/0.07/ 0.11, and for ReaxFF (oxides) 0.09/0.09/0.11. The fourth peak at 3.0 Å by VASP is sharper than the rest. Nonetheless, all these values are in good agreement and indicate that both force fields are able to reproduce the structure of the interfacial water quantitatively well. Finally, in Figure 9b we compare the solid structure (Fe-Fe bonds) between the three sets of simulations. The main discrepancy here is that the second peak by VASP close to 3.5 Å has been suppressed appreciably, and in ReaxFF (full) it is

6306

J. Phys. Chem. A, Vol. 114, No. 21, 2010

Figure 9. Structural analysis of water on goethite (010) MD simulations by VASP and ReaxFF, (a) radial distribution function for O and H atoms at the solid-solution interface, (b) radial distribution function between Fe atoms in the solid phase.

almost missing from the classical force field simulations. On the other hand, the two force fields have been successful in reproducing structural features at 3.1, 4.6, and 5.5 Å, very close to the ab initio results. A visual examination of the model indeed confirms the stability of the solid structure in our classical simulations. Again, discrepancies in the RDFs between VASP and ReaxFF can, at least partially, be attributed to differences in H-bonding between (H)OH-O and (Fe)OH-O. Conclusion The results of training two force fields for iron oxides/ oxyhydroxides within the ReaxFF methodology have been presented. This methodology, by going beyond the equilibrium conditions of a system, allows chemical reactions and bond formation/dissociation. Our training data set included and

Aryanpour et al. presented here: two deprotonation reactions of Fe2+ and Fe3+; dissociation of Fe-dimer in solution, and thermodynamics of a selected set of iron minerals. Two training procedures were followed: for the first force field, denoted by “full” we added a previously set of Fe compounds with OH and H2O groups to our training set. For the second force field, denoted by “oxides” we used the data set presented in this paper only. The test set consisted of estimating the lattice parameters of the selected iron minerals, plus a system of goethite and water. We showed that both force fields can reproduce the target data of the training set reasonably well. In the deprotonation reaction Fe2+, a similar increase in energetics was predicted by the full and oxides force fields in agreement with quantum calculations. In the deprotonation reaction Fe3+, ReaxFF (full) and (oxides) overpredicted the activation energy by ∼20 and ∼10 kJ/mol, respectively. In the dissociation reaction of Fedimer, ReaxFF (full) underpredicted the activation energy by ∼15 kJ/mol, whereas ReaxFF (oxides) could give almost the same activation energy as given by QC. The relative stabilities of iron-oxide mineralsshematite, goethite, lepidocrocite, akagane´ite, wu¨stite, magnetiteswere predicted correctly by our two force fields. To obtain a better accuracy for the absolute stabilities, further work seems necessary. The training results for the expansion and contraction of bulk hematite and goethite were quite satisfactory, showing the reliability of the ReaxFF force fields for such studies on other minerals. The results for the test cases are promising in bulk and surface applications of these force fields. The lattice constants for the selected set of minerals were predicted within at most 8% of the experimental measurements. The exception was a larger deviation in only one mineral (lepidocrocite) for which Hbonding plays a crucial role in between its bulk layers. For the goethite-water systems, we found ReaxFF (full) and (oxides) predict the structure of solution in a very good agreement with ab initio calculations and experimental studies. The bonding profile of the solid phase predicted by ReaxFF misses one peak out of four peaks, although the overall solid structure is well preserved. Both the ReaxFF (full) and (oxides) force fields were able to predict the structure of solid-solution interface accurately enough. To summarize the training and testing results, we have shown the promising performance of our developed reactive force fields to study iron oxide systems. Furthermore, these results are very useful in future applications, as well as in future force field developments. Substantial improvements in the Fe-O-H ReaxFF parametrization could be made by expanding the training set used to fit the force field. For example, the hydrolysis series from Fe3+ (H2O)6 to Fe(OH)4- should be modeled. Larger explicit solvation shells should also be included in order to obtain more accurate energetics.63 For Fe2+(H2O)6, only the first deprotonation reaction is generally relevant in the environment, but simulations of wu¨stite and magnetite surfaces that involve the Fe2+ oxidation state would help ensure force field accuracy for this species in mineral-water interface reactions such as reductive dissolution of goethite via aqueous Fe2+. Previously unpublished results on metallic Fe phases and corrosion products (i.e., FeOH) have been included in the full training set, but extending these calculations over a range of densities will help produce a force field that is more reliable for redox reactions among Fe0, Fe2+, and Fe3+. In addition, DFT results used for testing the force field in this current study, such as the goethite-water interface and Fe-O-H nanoparticles, could be added to the full training

Reactive Force Field for Iron-Oxyhydroxide Systems set. Finally, expanding the range of ReaxFF parameters allowed to optimize on this extended training set will make it more flexible. Acknowledgment. M.A. and J.D.K. greatly appreciate financial support from the National Science Foundation (NSF) Collaborative Research in Chemistry grant (CHE-0958664) and by the Center for Environmental Kinetics Analysis (CEKA), an NSF- and DOE-sponsored Environmental Molecular Science Institute (NSF CHE-0431328). Supporting Information Available: This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Lee, G.; Faure, G. Water Air Soil Pollut. 2007, 186 (1-4), 221– 232. (2) Mohan, D.; Pittman, C. U. J. Hazard. Mater. 2007, 142 (1-2), 1–53. (3) Cundy, A. B.; Hopkinson, L.; Whitby, R. L. D. Sci. Total EnViron. 2008, 400 (1-3), 42–51. (4) Abdallah, E. A. M.; Gagnon, G. A. Can. J. CiVil Eng. 2009, 36 (5), 881–888. (5) Sumner, M. E. Handbook of Soil Science; CRC Press: Boca Raton, Fl., 2000. (6) Houben, G. J. Ground Water 2004, 42 (1), 78–82. (7) Navrotsky, A.; Mazeina, L.; Majzlan, J. Science 2008, 319 (5870), 1635–1638. (8) Cornell, R. M. Schwertmann, U. The Iron Oxides: Structure, Properties, Reactions, Occurences and Uses, 2nd ed.; WILEY-VCH, Verlag GmbH & Co., KGaA: Weinheim, 2003. (9) Villalobos, M.; Cheney, M. A.; Alcaraz-Cienfuegos, J. J. Colloid Interface Sci. 2009, 336 (2), 412–422. (10) Brown, G. E.; Henrich, V. E.; Casey, W. H.; Clark, D. L.; Eggleston, C.; Felmy, A.; Goodman, D. W.; Gratzel, M.; Maciel, G.; Mccarthy, M. I.; Nealson, K. H.; Sverjensky, D. A.; Toney, M. F.; Zachara, J. M. Chem. ReV. 1999, 99 (1), 77–174. (11) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108 (4), 1255–1266. (12) Wang, J. W.; Kalinichev, A. G.; Kirkpatrick, R. J. Earth Planet. Sci. Lett. 2004, 222 (2), 517–527. (13) Larentzos, J. P.; Greathouse, J. A.; Cygan, R. T. J. Phys. Chem. C 2007, 111 (34), 12752–12759. (14) Stillinger, F. H.; David, C. W. J. Chem. Phys. 1978, 69 (4), 1473– 1484. (15) Stillinger, F. H.; David, C. W. J. Chem. Phys. 1980, 73 (7), 3384– 3389. (16) Curtiss, L. A.; Halley, J. W.; Hautman, J.; Rahman, A. J. Chem. Phys. 1987, 86 (4), 2319–2327. (17) Halley, J. W.; Rustad, J. R.; Rahman, A. J. Chem. Phys. 1993, 98 (5), 4110–4119. (18) Rustad, J. R.; Hay, B. P.; Halley, J. W. J. Chem. Phys. 1995, 102 (1), 427–431. (19) Rustad, J. R.; Felmy, A. R.; Hay, B. P. Geochim. Cosmochim. Acta 1996, 60 (9), 1553–1562. (20) Rustad, J. R.; Felmy, A. R. Geochim. Cosmochim. Acta 2005, 69 (6), 1405–1411. (21) Cwiertny, D. M.; Hunter, G. J.; Pettibone, J. M.; Scherer, M. M.; Grassian, V. H. J. Phys. Chem. C 2009, 113 (6), 2175–2186. (22) Cygan, R. T.; Greathouse, J. A.; Heinz, H.; Kalinichev, A. G. J. Mater. Chem. 2009, 19 (17), 2470–2481. (23) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard III, W. A. J. Phys. Chem. A 2001, 105 (41), 9396–9409. (24) Chenoweth, K.; Cheung, S.; van Duin, A. C. T.; Goddard III, W. A.; Kober, E. M. J. Am. Chem. Soc. 2005, 127 (19), 7192–7202. (25) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q. S.; Xu, X.; Goddard III, W. A. J. Phys. Chem. A 2003, 107 (19), 3803–3811. (26) Raymand, D.; van Duin, A. C. T.; Baudin, M. Hermansson, K. Surf. Sci. . accepted for publication. (27) Cheung, S.; Deng, W. Q.; van Duin, A. C. T.; Goddard III, W. A. J. Phys. Chem. A 2005, 109 (5), 851–859. (28) Buehler, M. J.; van Duin, A. C. T.; Goddard III, W. A. Phys. ReV. Lett. 2006, 96 (9), 095505. (29) Nomora, K.; Kalia, R. K.; Nakano, A.; Vashista, P.; van Duin, A. C. T.; Goddard, W. A. Phys. ReV. Lett. 2007, 99, 148303.

J. Phys. Chem. A, Vol. 114, No. 21, 2010 6307 (30) van Duin, A. C. T., Bryantsev, V.; Su, J. Goddard III, W. A. manuscript in preparation. (31) Jiang, D. E.; Carter, E. A. Phys. ReV. B 2003, 67, Article No. 214103. (32) Chenoweth, K.; van Duin, A. C. T.; Persson, P.; Cheng, M. J.; Oxgaard, J.; Goddard III, W. A. J. Phys. Chem. C 2008, 112 (37), 14645– 14654. (33) de Vos, E. Burchart Studies on zeolites: Molecular mechanics, framework stability, and crystal growth. Ph.D. Thesis, Technische University: Delft (Netherlands), 1992. (34) http://lammps.sandia.goV/. (35) Stumm, W. Morgan, J. J. Aquatic Chemistry: An Introduction Emphasizing Chemical Equilibria in Natural Waters; JohnWiley & Sons, Inc.: New York, 1981. (36) Kubicki, J. D. J. Phys. Chem. A 2001, 105 (38), 8756–8762. (37) Martin, R. L.; Hay, P. J.; Pratt, L. R. J. Phys. Chem. A 1998, 102 (20), 3565–3573. (38) Frisch, M. J. Trucks, G. W. Schlegel, H. B. Scuseria, G. E. Robb, M. A. Cheeseman, J. R. Montgomery, Jr., J. A. Vreven, T. Kudin, K. N. Burant, J. C. Millam, J. M. Iyengar, S. S. Tomasi, J. Barone, V. Mennucci, B. Cossi, M. Scalmani, G. Rega, N. Petersson, G. A. Nakatsuji, H. Hada, M. Ehara, M. Toyota, K. Fukuda, R. Hasegawa, J. Ishida, M. Nakajima, T. Honda, Y. Kitao, O. Nakai, H. Klene, M. Li, X. Knox, J. E. Hratchian, H. P. Cross, J. B. Bakken, V. Adamo, C. Jaramillo, J. Gomperts, R. Stratmann, R. E. Yazyev, O. Austin, A. J. Cammi, R. Pomelli, C. Ochterski, J. W. Ayala, P. Y. Morokuma, K. Voth, G. A. Salvador, P. Dannenberg, J. J. Zakrzewski, V. G. Dapprich, S. Daniels, A. D. Strain, M. C. Farkas, O. Malick, D. K. Rabuck, A. D. Raghavachari, K. Foresman, J. B. Ortiz, J. V. Cui, Q. Baboul, A. G. Clifford, S. Cioslowski, J. Stefanov, B. B. Liu, G. Liashenko, A. Piskorz, P. Komaromi, I. Martin, R. L. Fox, D. J. Keith, T. Al-Laham, M. A. Peng, C. Y. Nanayakkara, A. Challacombe, M. Gill, P. M. W. Johnson, B. Chen, W. Wong, M. W. Gonzalez, C. Pople. J. A. Gaussian 03, ReVision C.02.: Gaussian, Inc.,Wallingford, CT, 2004. (39) Kresse, G.; Hafner, J. Phys. ReV. B 1993, 47, 558. (40) Kresse, G.; Hafner, J. Phys. ReV. B 1994, 49, 14251. (41) Kresse, G.; Furthmller, J. Comput. Mater. Sci. 1996, 6, 15. (42) Kresse, G.; Furthmller, J. Phys. ReV. B 1996, 54, 11169. (43) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1997, 78, 1396. (45) Jiang, D. E.; Carter, E. A. Surf. Sci. 2003, 547 (1-2), 85–98. (46) Marcus, Y. Chem. ReV. 1988, 88 (8), 1475–1498. (47) Price, G. D.; Ross, N. L., Eds. The Stability of Minerals, 1st ed.; Chapman & Hall in association with the Mineralogist Society of Great Britain and Ireland, London: 1992. (48) Rozenberg, G. K.; Dubrovinsky, L. S.; Pasternak, M. P.; Naaman, O.; Le Bihan, T.; Ahuja, R. Phys. ReV. B 2002, 65 (6), 064112. (49) Huynh, T.; Tong, A. R.; Singh, B.; Kennedy, B. J. Clays Clay Miner. 2003, 51, 397. (50) Cerius2Materials Science V. 4.9; Accelrys Inc.: San Diego, CA, 2004. (51) Fjellvag, H.; Hauback, B. C.; Vogt, T.; Stolen, S. Am. Mineral. 2002, 87, 347–349. (52) Post, J. E.; Heaney, P. J.; Von Dreele, R. B.; Hanson, J. C. Am. Mineral. 2003, 88, 782–788. (53) Kuriki, A.; Moritomo, Y.; Ohishi, Y.; Kato, K.; Nishibori, E.; Takata, M.; Sakata, M.; Hamada, N.; Todo, S.; Mori, N.; Shimomura, O.; Nakamura, A. J. Phys. Soc. Jpn. 2002, 71, 3092–3093. (54) Olsen, J. S.; Cousins, C. S. G.; Gerward, L.; Jhans, H.; Sheldon, B. J. Phys. Scr. 1991, 43 (3), 327–330. (55) Gleason, A. E.; Jeanloz, R.; Kunz, M. Am. Mineral. 2008, 93 (1112), 1882–1885. (56) Sjoden, O.; Seetharaman, S.; Staffansson, L. I. Metall. Trans. B 1986, 17 (1), 179–184. (57) Stolen, S. Gronvold, F. J. Geophys. Res., [Solid Earth], 101:1153111540. (58) Faure, G. Principles and Applications of Geochemistry, 2nd ed.; Prentice Hall: Upper Saddle River, N.J., 1998. (59) Haavik, C.; Stolen, S.; Fjellvag, H.; Hanfland, M.; Hausermann, D. Am. Mineral. 2000, 85 (3-4), 514–523. (60) Kubicki, J. D.; Paul, K. W.; Sparks, D. L. Geochem. Trans. 2008, 9, 4. (61) Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Phys. ReV. B 1998, 1505. (62) Paesani, F.; Voth, G. A. J. Phys. Chem. B 2009, 113 (17), 5702– 5719. (63) Rustad, J. R.; Dixon, D. A. J. Phys. Chem. A 2009, 113 (44), 12249–12255.

JP101332K