Evaluation of Density Functional Theory Methods for Studying

Theoretical studies of arsenite adsorption and its oxidation mechanism on a perfect TiO2 anatase (101) surface. Applied Surface Science 2011, 258, 119...
0 downloads 0 Views 337KB Size
Download PDF
Environ. Sci. Technol. 2005, 39, 4816-4822

Evaluation of Density Functional Theory Methods for Studying Chemisorption of Arsenite on Ferric Hydroxides NIANLIU ZHANG, PAUL BLOWERS, AND JAMES FARRELL* Department of Chemical and Environmental Engineering, University of Arizona, Tucson, Arizona 85721

Understanding adsorption of arsenic on ferric hydroxide surfaces is important for predicting the fate of arsenic in the environment and in designing treatment systems for removing arsenic from potable water. This research investigated the binding of arsenite to ferric hydroxide clusters using several density functional theory methods. Comparison of calculated and experimentally measured As-O and As-Fe bond distances indicated that As(III) forms both bidentate and monodentante corner-sharing complexes with Fe(III) octahedra. Edge-sharing As(III) complexes were less energetically favorable and had As-O and As-Fe distances that deviated more from experimentally measured values than corner-sharing complexes. The hydrated bidentate complex was the most energetically favorable in the vacuum phase, while the monodentate complex was most favored in the aqueous phase. Structures optimized using the Harris and Perdew-Wang local functionals were close to both experimental data and structures optimized using the nonlocal Becke-Lee-YangParr (BLYP) functional. Binding energies calculated with the gradient-corrected BLYP functional were only weakly dependent on the method used for geometry optimization. The approach of using low-level structures coupled with higher level single-point energies was found to reduce computational time by 75% with no loss in accuracy of the computed binding energies.

Introduction Toxic arsenic compounds are often present in water supplies due to the weathering of rocks, industrial waste discharges, and agricultural use of arsenical herbicides and pesticides (1). Under the pH and redox conditions of most groundwaters and surface waters, dissolved arsenic exists as the As(V) (arsenate) species H2AsO4- and HAsO42- and the As(III) (arsenite) species H3AsO30 and H2AsO3- (2). Both As(V) and As(III) are strongly absorbed by ferric oxides and hydroxides, such as hematite, goethite, and lepidocrocite (3-16). Adsorption to ferric hydroxides is often used for removing arsenic from potable water, and also affects the fate and transport of arsenic in soil and groundwater (17, 18). Although As(III) species are more toxic and water soluble than As(V) species, most prior studies investigating arsenic binding to ferric hydroxides have focused on As(V) (19), and little is known about the mechanisms of As(III) adsorption with ferric hydroxides. * Corresponding author phone: (520) 621-2465; (520) 621-6048; e-mail: [email protected]. 4816

9

fax:

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 13, 2005

Previous extended X-ray absorption fine structure (EXAFS) studies suggest that As(III) adsorbs on ferric hydroxides by forming inner-sphere surface complexes by ligand exchange with OH- or OH2 groups (3, 8). However, the type of surface complex could not be positively identified from the EXAFS results. Adsorption of As(III) on goethite was first observed by Manning et al. (3), who proposed that As(III) formed bidentate, binuclear surface complexes on ferric hydroxides. Additionally, Raven et al. (6) proposed that a monodenate bonding mechanism might play an important role during As(III) adsorption at high pH values. Computational chemistry methods, such as density functional theory, are ideally suited for evaluating the different binding modes of arsenite with ferric hydroxides. Density Functional Theory. Density functional theory (DFT) is a quantum mechanical modeling approach based on expressing the total energy of a system as a functional of the electron density (20-23). The accuracy of results from DFT is comparable to that of correlated Hartree-Fock methods, but requires substantially less computational effort (24). Because of its computational advantages, DFT has evolved as the most important quantum mechanical approach in solid-state physics (25, 24), and can handle large systems that are intractable by correlated Hartee-Fock methods (26, 27). DFT methods are also particularly useful for transition-metal compounds where molecular orbital calculations often give poor results (28, 29). DFT has proven to be useful for studying metal-ligand interactions in organometallic and coordination complexes (30, 31), and for studying adsorption processes on metal (32) and metal oxide (33, 34) surfaces. DFT is based on the fact that the ground-state properties of a molecular system are uniquely defined by the electron density (F(r)) (22). In DFT the energy functional (E[F]) is written as a sum of four terms (35):

ZA

M

E[F(r)] ) -

∑∫|r - R |F(r) dr +

( )

A)1 N

∑∫ψ (r) i

i)1

32 2

A

ψi(r) dr +

1

F(r1)F(r2)

∫∫ |r 2

1

- r 2|

dr1 dr2 + Exc[F(r)] (1)

where ZA is the charge on nucleus A, r is a position vector, RA is the position vector of nucleus A, M is the number of nuclei, N is the number of electrons, ψi is the spatial orbital occupied by the ith electron, and Exc[F(r)] is the exchangecorrelation energy functional. The first term on the righthand side of eq 1 corresponds to the potential energy due to nucleus-electron interactions, the second term is the kinetic energy for N noninteracting electrons, and the third term represents the energy for electron-electron repulsions. Internuclear repulsions are constant and are accounted for separately. Equation 1 serves to define Exc[F(r)] as the energy needed to be added to the other terms to reproduce the exact energy of the molecular system. Thus, Exc[F(r)] contains contributions associated with electron exchange and correlation, and also contains a contribution arising from the difference between the true kinetic energy and the kinetic energy of N noninteracting electrons. The functional approximation used for Exc[F(r)] is the distinguishing feature of different DFT methods. The simplest approximation to Exc[F(r)] is the local density approximation (LDA), which assumes that Exc[F(r)] for the molecular system is the same as that for a homogeneous 10.1021/es050271f CCC: $30.25

 2005 American Chemical Society Published on Web 05/26/2005

electron gas (23). This assumption makes Exc[F(r)] dependent on only the local value of the electron density. Because of including the exchange-correlation term, LDA approximations are more accurate than Hartree-Fock approximations with similar computational cost (24). Because the electron density in a real molecule varies with position, Exc[F(r)] is dependent not only on the local electron density, but also on the electron density at other points in the system. To get a more accurate approximation of the exchange-correlation energy, Exc[F(r)] functionals that include both the electron density and the gradient of the electron density have been developed. Methods using these functionals are called generalized gradient approximations (GGAs). Experience has shown that the LDA can be used to determine accurate molecular geometries (36, 37) and harmonic frequencies (38, 39). However, the energies calculated by LDA functionals may contain large errors, even when the particular functional yields accurate molecular structures (40). Previous studies for systems composed of main group elements have shown that GGA calculations generally yield more accurate energies than LDA calculations, but do not always produce more accurate molecular geometries than LDA methods (41). However, for the limited range of transition-metal systems that have been studied, such as carbonyl complexes, GGA methods normally yield more accurate approximations than LDA methods for both geometries and energies (29). In addition to LDA and GGA approximations for Exc[F(r)], other methods have been developed that can be applied to systems that are too large to be treated using most LDA or GGA functionals. One such method that greatly reduces computational time was developed by Harris (42). In the Harris approximation, atoms in the molecular system are treated as isolated fragments and a correction is included to account for this approximation. The Harris method greatly reduces the amount of self-consistency cycling that is necessary to solve the Kohn-Sham equations (42). Although the Harris approximation was developed for weakly interacting systems, it has been found to give useful results for strong covalent bonds and even ionic systems (43-47). Unlike ab initio quantum chemical methods, DFT is not based purely on fundamental theory, which makes it necessary to test and evaluate each DFT method for a particular application (41). Although a hybrid molecular orbital (MO)/DFT approach has recently been used for comparing bidentate As(III) binding to aluminum and iron hydroxides (48), there have been no prior studies using DFT methods for determining the different binding modes of As(III) on ferric hydroxides. The objectives of this study were to determine the most favorable binding structures for different As(III) species with a ferric hydroxide cluster, and to compare the accuracies and computational costs of several DFT methods for studying these systems. Modeling adsorption complexes of As(III) will be useful in designing experiments, for interpreting experimental and spectroscopic data, and for predicting the fate of As(III) under various environmental conditions.

Methods DFT calculations were performed using the Dmol3 (49) package in the Accelrys Materials Studio modeling suite (50). The calculations were performed on a personal computer operating with a 2.80 GHz Pentium II processor. All simulations used double-numeric with polarization (DNP) basis sets (49, 51), and the nuclei and core electrons were described with DFT semilocal pseudopotentials (52). Geometry optimizations and energies were calculated using the Harris approximation (42), Perdew-Wang (PWC) LDA functional (53), and Becke-Lee-Yang-Parr (BLYP) GGA functional (54, 55). The effects of spin polarization were investigated by

FIGURE 1. Fully hydrated cluster used to simulate ferric hydroxide (Fe2(OH)2(H2O)84+).

TABLE 1. AsH3 Optimized Geometries and Vibrational Frequencies Expressed as Wavenumbers for Each Functional Compared with Experimental Data parameter

Harris

LDA

GGA

experiment

As-H length (Å) H-As-H angle (deg) wavenumber (cm-1)

1.496 89.6 880 1001 2101 2149

1.534 90.53 878 959 2124 2146

1.543 91.3 917 981 2114 2131

1.523 (59) 92 (59) 906 (60) 1003 (60) 2116 (60) 2123 (60)

comparing calculations using spin-restricted and unrestricted methods. In the unrestricted calculations, the spin state was determined by comparing the energies of optimized structures with different spin states. All geometries were optimized in the gas phase without any constraints, and aqueous solvation was incorporated using the conductor-like screening model (COSMO) (56). Single-point energies for geometries optimized using Harris, LDA, and restricted GGA methods were also calculated using an unrestricted GGA approach. Vibrational frequencies were calculated using the harmonic approximation and corrected for anharmonicity using scaling factors (57). Ferric hydroxide was simulated using a cluster consisting of two iron atoms in octahedral coordination with 10 oxygen atoms, as illustrated in Figure 1. This cluster is identical to the one used by Sherman and Randall (34) for simulations of As(V) binding to ferric hydroxides. Three different binding geometries were investigated for reactions of different possible As(III) species with ferric hydroxide. Binding geometries for bidentate corner sharing, bidentate edge sharing, and monodentate corner sharing were studied. The As(III) species used in the calculations were hydrated (H3AsO3) and unhydrated (HAsO2) arsenious acid, metaarsenite (AsO2-), hydrogen arsenite (HAsO32-), and dihydrogen arsenite (H2AsO3-). These species were selected for their relevance under environmental conditions. The pKa of hydrated arsenious acid is 9.2, and the ratio of hydrated to unhydrated arsenious acid is approximately 9:1 under environmental conditions (58).

Results and Discussion Prior to calculations of As(III) species binding to ferric hydroxide being performed, the three different DFT functionals, the DNP basis set, and core pseudopotentials were validated for arsenic by calculating molecular geometries and vibrational spectra of AsH3. This species was selected to validate each DFT implementation because it was the only arsenic species for which both the geometry and vibrational spectra could be obtained. The optimized geometries and vibrational frequencies of AsH3 calculated using spinrestricted Harris, LDA, and GGA functionals are shown in Table 1. In all cases, the DFT calculations are in good agreement with the experimental results for both the geometry and VOL. 39, NO. 13, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

4817

FIGURE 2. As(III) surface complexes with the three lowest energies: bidentate complex A, AsO2dFe2(OH)2(H2O)63+; bidentate complex B, HAsO3dFe2(OH)2(H2O)62+; monodentate complex C, H2AsO3-Fe2(OH)2(H2O)73+. vibrational frequency. The Harris functional underestimated the As-H bond length by 0.027 Å (1.77%) and underestimated the H-As-H angle by 2.4° (2.61%) as compared to the experimental results. In contrast, the LDA and GGA functionals overestimated the As-H bond length by about 0.02 Å (1.31%) and overestimated the H-As-H bond angle by 1.5° (1.63%). GGA and LDA methods gave better geometry approximations than the Harris approach, but the GGA results did not show an obvious improvement as compared to the LDA method. The scaled GGA frequencies are closer to experimental results than the LDA and Harris frequencies. These results indicate that the DFT implementations used here can reasonably describe the properties of arsenic species. The accuracy of each method for calculating the binding geometry of three different surface complexes using the Harris, LDA, and GGA methods was determined by comparing As-O and Fe-O bond lengths and As-Fe distances with experimental data. EXAFS data for As(III) adsorbed on goethite (R-FeOOH) suggest that there is a monodentate complex with an Fe-As distance of 3.46 Å, and a bidentate complex of As surrounded by three oxygen atoms at 1.78 Å and two Fe atoms at 3.34 Å (3, 8). The complexes for which the calculated As-O and As-Fe distances were closest to the experimental data are shown in Figure 2. In all three complexes As(III) is bound at corner sites on the iron octahedra. Edge-sharing As(III) complexes were less energetically favorable and had As-O and Fe-As distances that deviated more from experimentally measured values than those of corner-sharing complexes. Therefore, all subsequent calculations focused on these three complexes. Bond lengths in complexes A-C calculated using different methods are listed in Table 2 and compared in Figure 3. The As-O bond lengths calculated by all methods were within 3.0% of the experimental value of 1.78 Å, with the exception of the Harris calculation for complex C. The GGA approach consistently predicted longer As-O bonds than the LDA approach, but both methods had errors of similar magnitude as compared to experiment. For complexes A and B, the Harris approximation gave As-O bond lengths closer to experiment than the LDA and GGA methods. However, for complex C, the Harris approach had a much larger error than the other methods. For the Fe-As bond lengths, the GGA approach was consistently better than the LDA and Harris methods, which had similar errors. The LDA calculations consistently underestimated As-Fe bond lengths by 4-10% (0.33-0.12 Å). This result is consistent with prior studies showing LDA methods generally underestimate metal-ligand bond lengths (41). The GGA functional improved over the LDA results for As-Fe bond lengths with errors ranging from 1% to 4% (0.13-0.03 Å). Although there are no experimental results for Fe-O bond lengths for As(III) complexes with ferric iron that can be compared to the calculated results, Fe-O bond lengths in 4818

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 13, 2005

TABLE 2. Bond Lengths Calculated with the Harris, LDA Spin-Restricted (LDA-re) and Unrestricted (LDA-un), and GGA Spin-Restricted (GGA-re) and Unrestricted (GGA-un) Methodsa bond lengths (Å) structure

Harris

LDA-re

As-Fe As-O Fe-O

3.262 1.761 1.745

3.216 1.728 1.816

As-Fe As-O Fe-O

3.183 1.796 1.732

As-Fe As-O Fe-O

3.598 2.012 1.842

LDA-un

GGA-re

GGA-un

EXAFS

Complex A 3.214 3.317 1.734 1.752 1.822 1.882

3.307 1.75 1.872

3.34 1.78 1.88b

3.014 1.761 1.785

Complex B 3.103 3.21 1.782 1.794 1.775 1.834

3.299 1.822 1.814

3.34 1.78 1.88b

3.342 1.764 1.887

Complex C 3.344 3.422 1.781 1.84 1.885 1.879

3.423 1.829 1.905

3.46 1.78 1.88b

a Also shown are bond lengths for As(III)-ferric hydroxide complexes measured using EXAFS (8). b Measured Fe-O bond length in Fe3O4 (18).

other iron oxides have been measured at 1.88 Å (18). The GGA-calculated Fe-O bond lengths most closely approximate this value. For the Fe-O bond lengths, there is a general trend with the GGA method yielding the largest value and the Harris method yielding the shortest. In all cases, the calculated values were close to 1.88 Å, with the largest errors less than 8%. Surface complexes formed with ferric oxides may be openshell systems due to iron’s unsaturated electronic configuration (41). This means that calculation methods allowing unpaired electrons may be required because spin-restricted calculations may lead to mixing configurations with high spin states (61). For all three complexes, the unrestricted LDA calculations produced the lowest energy for a spin state of 2, while the GGA calculations produced the lowest energy for a spin state of 0. Figure 3 compares the bond lengths calculated with each method and indicates that differences between spin-restricted and unrestricted calculations were small, and thus restricted calculations can be used for geometry optimizations for this system. This will save substantial computational time since unrestricted calculations may take more than twice as long as spin-restricted methods. The results in Figure 3 also show that the Harris approximation can yield reasonably good optimized geometries for this system. The ability to use the Harris approximation will be important in modeling clusters larger than those used in this work, such as those needed to investigate interactions between adsorbed arsenic species. In eight of the nine bond length calculations, the Harris method, with a mean error of 3.8%, gave results that were similar in accuracy to those of the LDA methods, which had a mean error of 3.4%. However,

FIGURE 3. Comparison of calculated As-O, As-Fe, and Fe-O bond lengths in surface complexes A-C using different methods. the error in the Harris approximation for the complex C As-O bond length of 13% was the largest error in any of the calculations. This error is larger than those previously reported for the Harris approximation, which were found to be less than 9% for diatomic molecules (44). Thus, although the average accuracy of the Harris method is on par with that of the LDA/PWC method, it may sometimes yield results with large errors. Energies. Relative energies of the three complexes calculated using the Harris, LDA, and GGA methods are compared in Figure 4. The unrestricted GGA calculations gave the lowest total energy for all three complexes. This is consistent with past results showing that GGA functionals usually deliver much better energy approximations than LDA methods (41). However, in past work with transition-metal clusters, GGA functionals have been found to yield binding energies with errors up to 50% as compared to experimental values (62). Figure 4 shows the relative energy error for each method using the unrestricted GGA calculation as the reference value. Also shown in Figure 4 are the relative errors associated with the unrestricted GGA single-point energies for each cluster calculated using the structures determined with lower level methods. Relative errors in energy calculated with both the Harris and LDA methods were approximately 0.8% for all three clusters. However, GGA single-point energies calculated for the Harris and LDA structures gave much better approximations to the energy calculated using the unrestricted GGA geometry. This indicates that high-level singlepoint energy calculations can be paired with lower level

FIGURE 4. Energy errors for each method relative to the unrestricted GGA energy for each cluster. GGA-un/GGA-un energies for each cluster are (A) -3.61 × 106 kJ/mol, (B) -3.81 × 106 kJ/mol, and (C) -4.01 × 106 kJ/mol. Also shown are the relative errors associated with single-point unrestricted GGA energies calculated using geometries optimized by each method.

TABLE 3. Computational Time (min) for Complex A Geometry Optimizations Using Different Methods and Unrestricted GGA Single-Point Energies for Each Structure Harris LDA-re LDA-un GGA-re GGA-un geometry optimization 62.8 single-point energy (GGA-un) 19.6

204 16.5

528 18.9

627 20.3

874 18.2

geometry optimizations to yield energies similar to those obtained with high-level geometry optimizations. This should result in considerable savings in computational time since geometry optimizations take much longer than single-point energy calculations. The computational times for complex A using different methods are shown in Table 3. Cluster geometries optimized using the Harris, LDA, and GGA methods each required approximately 20 min to calculate single-point energies using the unrestricted GGA method. However, the time for geometry optimization varied by a factor of 14 between the different methods. Time differences between the unrestricted GGA optimization and other methods would increase as the size of the cluster increases. Single-point energy calculations with the unrestricted GGA method for the LDA geometries VOL. 39, NO. 13, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

4819

lowered the average energy error from 0.84% to 0.0021% for the LDA spin-restricted geometries, and from 0.84% to 0.0019% for unrestricted LDA geometries. Single-point energy calculations with the unrestricted GGA functional for the Harris geometries reduced the average energy error from 0.78% to 0.0091%. Another area where computational time can be saved is by performing restricted spin calculations. Unrestricted and spin-restricted calculations yielded energies that differed by less than 0.0008%, which was less than 34 kJ/mol for all clusters. This indicates that spin-restricted calculations can give reasonable energies for these open-shell transition-metal complexes. Binding energies can be used to determine the most energetically favorable surface complex, and for determining the relative affinity of different species for adsorption sites on the iron oxide. The binding energy for a surface complex is defined as the energy difference between the reactants and products for the adsorption reaction. The reactions for the two bidentate complexes are

complex A: HAsO20 + Fe2(OH)2(H2O)84+ T

AsO2-Fe2(OH)2(H2O)63+ + H2O + H3O+ (2)

complex B: H3AsO30 + Fe2(OH)2(H2O)84+ T

HAsO3-Fe2(OH)2(H2O)62+ + 2H3O+ (3)

and for the monodentate complex the reaction is

complex C: H3AsO30 + Fe2(OH)2(H2O)84+ T

H2AsO3-Fe2(OH)2(H2O)73+ + H3O+ (4)

The As(III) binding energies for each complex were calculated by subtracting the total energy of the reactant species from the total energy of the product species in both the vacuum and aqueous phases. Caution must be used in interpreting the binding energies for the aqueous-phase reactions. Past research shows that polarizable continuum models such as COSMO can yield large errors when the solvation energies of ions or hydrogen-bonding species are calculated (63). Determining accurate hydration energies for ionic and hydrogen-bonding species requires explicit treatment of one to several solvation shells around each solute (64). Explicit inclusion of a sufficient number of water molecules for solvating the large complexes in this study is not feasible with the available computational resources. Although there may be some cancellation of hydration energy errors because the charged iron oxides appear on both sides of each reaction, hydration energy errors for H3AsO3 or HAsO2 and H3O+ are not likely to cancel, and may contribute significant error to the calculated aqueous-phase binding energies. For example, the hydration free energy calculated for H3O+ of -385 kJ/mol is 51 kJ/mol less exergonic than the experimentally measured value of -435 kJ/mol (65). However, hydration energy errors for the uncharged H3AsO3(aq) and HAsO2(aq) species are likely to be much less than that for H3O+(aq). Systematic errors in the calculations notwithstanding, the gas-phase binding energies (∆E) indicate that bidentate complex B is the most energetically favorable. This result agrees with EXAFS data for goethite suggesting that As(III) forms bindentate, binuclear complexes with octahedrally coordinated Fe atoms (3). However, the aqueous-phase ∆E values, standard enthalpies (∆H°), and Gibbs free energies (∆G°) for the binding reactions suggest that monodentate complex C is the most energetically favorable. The ∆G° values for the bidentate complexes in Table 4 are considerably more exergonic that those calculated in a 4820

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 13, 2005

TABLE 4. Binding Energies, Standard Enthalpies, and Free Energies for As(III) Complexation Reactionsa in vacuo

in solution

reaction

complex

∆E

∆E

∆H°(298 K)

∆G°(298 K)

2 3 4

A B C

-814 -1208 -919

-629 -349 -727

-629 -337 -712

-667 -366 -699

a Thermal corrections to ∆E values were based on gas- and aqueousphase vibrational spectra and the change in mole numbers for the reactions.

TABLE 5. Binding Energy (∆E) Errors for Each Method Using the GGA-un/GGA-un Binding Energy as the Reference Valuea energy error (kJ/mol)/(%) method (geometry/energy)

complex A

complex B

complex C

average

Harris/Harris Harris/GGA-un LDA-re/LDA-re LDA-re/GGA-un LDA-un/LDA-un LDA-un/GGA-un GGA-re/GGA-re GGA-re/GGA-un GGA-un/GGA-un

292/+282 -3/-3.2 -30/-29.0 -2/-1.6 -17/-16.3 0.6/+0.6 8/+7.7 -0.7/-0.69 0

435/+206 39/+18 -54/-25.5 6/+2.96 -45/-21.1 5/+2.5 18/+8.7 6/+2.7 0

444/+199 234/+105 -106/-47.8 -19/-8.66 -113/-50.6 -14/-6.5 2/+1.1 21/+9.5 0

389/229 92/42 63/34 9/4.4 58/29 7/3.2 9/5.8 9/4.3 0

a

Also shown are the average of the absolute errors for each method.

previous MO/DFT investigation which ranged from +14 to -59 kJ/mol (48). Part of this disparity can be attributed to differences in computational techniques and to modeling different complexation reactions. Experimental evidence suggests that the binding free energies for As(III) with iron oxides are on the order of several hundred kilojoules per mole. Dixit and Hering (11) reported monondentate binding constants for arsenite on several iron oxides of ∼1039, which corresponds to a ∆G° value of -220 kJ/mol. Thus, considering that different reactions were modeled, and that DFT/GGA binding energies for transition-metal compounds have been found to be in error by up to 50% (62), the binding energy values in this study are in reasonable agreement with experimental data. Table 5 lists the vacuum binding energy errors for each method assuming that unrestricted GGA binding energies calculated for complexes with geometries optimized using the unrestricted GGA method (GGA-un/GGA-un) are the correct values. Harris binding energies calculated for Harrisoptimized structures had errors ranging from 290 to 440 kJ/mol. This is consistent with past speculation that the energies calculated with the Harris method are not reliable for strong covalent bonds (42). LDA binding energies calculated for LDA-optimized structures consistently yielded binding energies that were 16-50% more exergonic than the unrestricted GGA binding energies. This is consistent with past observations that LDA methods predict stronger binding energies than GGA methods (41). Unrestricted GGA singlepoint energies calculated for structures optimized at lower levels greatly reduced the binding energy errors. With the exception of the Harris method, the maximum binding energy error arising from the lower level geometry optimizations was 21 kJ/mol (9.5%), with average errors ranging from 3.2% to 4.4%. For both the LDA and GGA methods, the unrestricted calculations improved both the geometry optimizations and the single-point energies. Comparison of errors between the LDA-re/GGA-un and LDA-un/GGA-un methods indicates that relaxing the spin restriction improved the geometries

slightly and reduced the average error in the single-point binding energy calculations from 4.4% to 3.2%. For the LDA methods, the improvement in the binding energy errors was 5 kJ/mol or less for unrestricted versus restricted geometry optimizations. Comparison of errors between the GGA-re/ GGA-re and GGA-re/GGA-un methods indicates that relaxing the spin restriction improved the single-point binding energies for complexes A and B, but increased the binding energy error for complex C. Results from this study indicate that a variety of DFT methods can be used to investigate As(III) complexes with ferric hydroxides. In all but one calculation, all methods yielded bond lengths that were within 10% of experimental values. This indicates that the Harris method may be useful for geometry optimization for large clusters where higher level methods may not be able to reach convergence in the self-consistent field iterations, or may be too time consuming. For geometry optimizations on large clusters, the Harris method may be able to reduce computational times by more than 2 orders of magnitude. Single-point energies calculated with the unrestricted GGA method yielded gas-phase binding energies for LDA geometries that were all within 21 kJ/mol of the binding energies calculated with the unrestricted GGA geometry optimization. This indicates that LDA methods can be used for the time-consuming geometry optimizations and unrestricted GGA methods are required only for single-point energy calculations. This approach will balance the requirements of calculation accuracy with the computational affordability, and may increase the tractable cluster size for which calculations can be performed. The structures investigated in this study represent only a small fraction of the possible binding modes of As(III) with ferric hydroxides. There may be other energetically favorable binding modes for As(III) on different faces of crystalline iron oxides. Additionally, the effects of structural defects and substitutions in the iron oxide, surface charge, and interactions between As(III) and other adsorbed species need to be considered to determine the full range of environmentally relevant As(III) complexes. However, the limited results from this study showing the presence of three energetically favorable binding modes for As(III) should be considered when experimental data are interpreted, or when surface complexation models for studying the adsorption behavior of As(III) compounds in complex mixtures are employed.

Acknowledgments Thanks are due to James Kubicki for helpful discussions and Eric Case for computer support. This project was made possible by Grant Number 2P42ES04940-11 from the National Institute of Environmental Health Sciences of the National Institutes of Health, with funds from the U.S. Environmental Protection Agency.

Literature Cited (1) Bhumbla, D. K.; Keefer, R. F. Arsenic mobilization and bioavalilability in soils. In Arsenic in the environment. Part I: Cycling and characterization; Nriagu, J. O., Ed.; John Wiley & Sons: New York, 1994; pp 51-82. (2) Cullen, W. R.; Reimer, K. J. Arsenic speciation in the environment. Chem. Rev. 1989, 89, 713-764. (3) Manning, B. A.; Fendorf, S. E.; Goldberg, S. Surface structures and stability of arsenic(III) on goethite: Spectroscopic evidence for inner-sphere complexes. Environ. Sci. Technol. 1998, 32, 2383-2388. (4) Farrell, J.; Wang, J.; O’Day, P.; Conklin, M. Electrochemical and spectroscopic study of arsenate removal from water using zerovalent iron media. Environ. Sci. Technol. 2001, 35, 20262032. (5) Goldberg, S.; Johnston, C. T. Mechanisms of arsenic adsorption on amorphous oxides evaluated using macroscopic measurements, vibrational spectroscopy, and surface complexation modeling. J. Colloid Interface Sci. 2001, 234, 204-216.

(6) Raven, K. P.; Jain, A.; Loeppert, R. H. Arsenite and arsenate adsorption on ferrihydrite: Surface charge reduction and net OH- release stoichiometry. Environ. Sci. Technol. 1999, 33, 1179-1184. (7) Roberts, L. C.; Hug, S. J.; Ruettimann, T.; Billah, M. M.; Khan, A. W.; Rahman, M. T. Arsenic removal with iron(II) and iron(III) in water with high silicate and phosphate concentrations. Environ. Sci. Technol. 2004, 38, 307-315. (8) Manning, B. A.; Hunt, M. L.; Amrhein, C.; Yarmoff, J. A. Arsenic (III) and arsenic(V) reactions with zerovalent iron corrosion products. Environ. Sci. Technol. 2002, 36, 5455-5461. (9) Appelo, C. A.; Weiden, M. J. J. V. D.; Tournassat, C.; Charlet, L. Surface complexation of ferrous iron and carbonate on ferrihydrite and the mobilization of arsenic. Environ. Sci. Technol. 2002, 36, 3096-3103. (10) Su, C.; Puls, R. W. In situ remediation of arsenic in simulated groundwater using zerovalent iron: Laboratory column tests on combined effects of phosphate and silicate. Environ. Sci. Technol. 2003, 37, 2582-2587. (11) Dixit, S.; Hering, J. G. Comparison of arsenic(V) and arsenic(III) sorption onto iron oxide minerals: Implications for arsenic mobility. Environ. Sci. Technol. 2003, 37, 4182-4189. (12) Raven, K. P.; Jain, A.; Loeppert, R. H. Arsenite and arsenate adsorption on ferrihydrite: Kinetics, equilibrium, and adsorption envelopes. Environ. Sci. Technol. 1998, 32, 344-349. (13) Waychunas, G. A.; Fuller, C. C.; Rea, B. A.; Davis, J. A. Wideangle X-ray scattering (WAXS) study of “two-line” ferrihydrite structure: Effect of arsenate sorption and counterion variation and comparison with EXAFS results. Geochim. Cosmochim. Acta 1996, 60, 1765-1781. (14) Waychunas, G. A.; Davis, J. A.; Fuller, C. C. Geometry of sorbed arsenate on ferrihydrite and crystalline FeOOH: re-evaluation of EXAFS results and topological factors in predicting sorbate geometry, evidence for monodentate complexes. Geochim. Cosmochim. Acta 1995, 59, 3655-3661. (15) Waychunas, G. A.; Rea, B. A.; Fuller, C. C.; Davis, J. A. Surfacechemistry of ferrihydrite. 1. EXAFS studies of the geometry of coprecipitated and adsorbed arsenate. Geochim. Cosmochim. Acta 1993, 57, 2251-2269. (16) Fendorf, S.; Eick, M. J.; Grossl, P.; Sparks, D. L. Arsenate and chromate retention mechanisms on goethite: 1. Surface structure. Environ. Sci. Technol. 1997, 31, 315-320. (17) Su, C.; Puls, R. W. Significance of iron (II, III) hydroxycarbonate green rust in arsenic remediation using zerovalent iron in laboratory column test. Environ. Sci. Technol. 2004, 38, 52245231. (18) Melitas, N.; Wang, J.; Conklin M.; O’Day, P.; Farrell, J. Understanding soluble arsenate removal kinetics by zerovalent iron media. Environ. Sci. Technol. 2002, 36, 2074-2081. (19) Manceau, A. The mechanism of anion adsorption on iron oxides: evidence for the bonding of arsenate tetrahedra on free Fe(O,OH)6 edges. Geochim. Cosmochim. Acta 1995, 59, 36473653. (20) Fermi, E. A. Statistical method for the determination of some properties of atoms. II. Application to the periodic system of the elements. Z. Phys. 1928, 48, 73-79. (21) Thomas, L. H. The calculation of atomic fields. Proc. Cambridge Philos. Soc. 1926, 23, 542-548. (22) Hohenberg P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. B 1964, 136, 864-871. (23) Kohn, W.; Sham, L. Self-consistent equations including exchange and correlation effects. Phys. Rev. A 1965, 140, 1133-1138. (24) Wimmer, E. Computational materials design and processing: perspectives for atomistic approaches. Mater. Sci. Eng., B 1996, 37, 72-82. (25) Politzer, P.; Seminario, J. M.; Concha, M. C.; Murry, J. S. Some applications of local density functional theory to the calculation of reaction energetics. Theor. Chim. Acta 1993, 85, 127-136. (26) Dahl, J. P., Avery, J., Eds. Loacal Density Approximations in Quantum Chemsity and Solid State Physics; Plenum Press: New York, 1984. (27) Labonowski, J. K., Andzelm, J. W., Eds. Density Functional Methods in Chemistry; Springer-Verlag: New York, 1991. (28) Bauschlicher, C. W.; Ricca, A.; Partidge, H.; Langhoff, S. R. In Recent Advances in Density Funtional Methods. Part II; Chong, D. P., Ed.; World Scientific: Singapore, 1995. (29) Jonas, V.; Thiel, W. Theoretical study of the vibrational spectra on the transition metal carbonyls M(CO)6 [M)Cr, Mo, W], M(CO)5 [M)Fe, Ru, Os], and M(CO)4 [M)Ni, Pd, Pt]. J. Chem. Phys. 1995, 102, 8474-8484. VOL. 39, NO. 13, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

4821

(30) Durate, H. A.; Carvalho, S.; Paniago, E. B.; Almeida, W. B. Interaction of N-hydroxyacetamide with vanadate. A density functional study. J. Inorg. Biochem. 1998, 72, 71-77. (31) Torrent, M.; Deng, L.; Sola, M.; Ziegler, T. A density functional study of [2+3] versus [2+2] addition of ethylene to chromiumoxygen bonds in chromyl chloride. Inorg. Chem. 1998, 37, 13071314. (32) Neyman, K. M.; Rosch, N. Structural and vibrational properties of Ni(111)/NO adsorption complexes-ALCGTO-LDF model cluster investigation. Surf. Sci. 1994, 309, 1193-1199. (33) Matveev, A. V.; Neyman, K. M.; Pacchioni, G.; Rosch, N. Density functional study of M-4 clusters (M) Cu, Ag, Ni, Pd) deposited on the regular MgO(001) surface. Chem. Phys. Lett. 1999, 299, 603-612. (34) Sherman, D. M.; Randall, S. R. Surface complexation of arsenic(V) to iron(III) hydroxides: Structural mechanism from Ab Initio molecular geometries and EXAFS spectroscopy. Geochim. Cosmochim. Acta 2003, 67, 4223-4230. (35) Leach, A. R. Molecular Modeling Principles and Applications, 2nd ed.; Prentice Hall: New York, 2001. (36) Dickson, R. M.; Becke, A. D. Basis-set-free local densityfunctional calculations of geometries of polyatomic molecules. J. Chem. Phys. 1993, 99, 3898-3905. (37) Andzelm, J.; Wimmer, E. Density functional Gaussian-typeorbital approach to molecular geometries, vibrations, and reaction energies. J. Chem. Phys. 1992, 96, 1280-1303. (38) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. The performance of a family of density functional methods. J. Chem. Phys. 1993, 98, 5612-5626. (39) Zhou, M.; Wheeless, C. J. M.; Liu, R. Density fuctional theory study of vibrational spectra. 1. Performance of several density functional methods in predicting vibrational frequencies. Vib. Spectrosc. 1996, 12, 53-63. (40) Scheier, A. C.; Baker, J.; Andzelm, J. W. Molecular energies and properties from density functional theory: Exploring basis set dependence of Kohn-Sham equation using several density functionals. J. Comput. Chem. 1997, 18, 775-795. (41) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory, 2nd ed.; Wiley-VCH: Weinheim, Germany, 2000. (42) Harris, J. Simplified method for calculating the energy of weakly interacting fragments. Phys. Rev. B 1984, 31, 1770-1779. (43) Averill, F. W.; Painter, G. S. Harris functional and related methods for calculating total energies in density-functional theory. Phys. Rev. B 1990, 41, 10344-10353. (44) Kobayashi, K.; Kurita, N.; Kumahora, H.; Tago, K. Molecularbond-energy calculations based on the harris-functional approximation coupled with the generalized-gradient approximation. Phys. Rev. B 1992, 45, 11299-11304. (45) Yang, S. H. Ab Initio local-orbital density-functional method for transition metals and semiconductors. Phys. Rev. B 1998, 58, 1832-1838. (46) Sankey, O. F.; Demkov, A. A.; Windl, W.; Fritsch, J. H.; Lewis, J. P.; Fuentes-Cabrera, M. The application of approximate density functionals to complex systems. Int. J. Quantum Chem. 1998, 69, 327-340. (47) Bankhead, M.; Watson, G. W.; Hutching, G. J.; Scott, J.; Willock, D. J. Calculation of the energy profile for the fluorination of dichloromethane over an R-alumina catalyst. Appl. Catal., A 2000, 200, 263-274.

4822

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 13, 2005

(48) Kubicki, J. D. Comparison of As(III) and As(V) complexation onto Al- and Fe-hydroxides. In Arsenic in the Environment; O’Day, P., Vlassopoulos, D., Benning, L., Eds.; American Chemical Society Symposium Series; American Chemical Society: Washington, DC, 2005; in press. (49) Delley, B. An all-electron numerical-method for solving the local density functional for polyatomic-molecules. J. Chem. Phys. 1990, 92, 508-517. (50) Accelrys, Inc. Material Studio 2.2, San Deigo, CA, 2001. (51) Delley, B. Fast calculation of electronstatics in crystals and large molecules. J. Phys. Chem. 1996, 100, 6107-6110. (52) Delley, B. Hardness conserving semilocal pseudopotentials. Phys. Rev. B 2002, 66, Art. No. 155125(1-9). (53) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron gas correlation energy. Phys. Rev. B 1992, 45, 13244-13249. (54) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 30983100. (55) Lee, C.; Yang, W.; Parr, R. G. Development of the colle-salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785-789. (56) Klampt, A. Conductor-like screening model for real solvents: A new approach to quantitative calculation of solvation phenomena. J. Phys. Chem. 1995, 99, 2224-2235. (57) Lewars, E. Computational Chemistry: Introduction to the Theory and Applications of Molecular and Quantum Mechanics; Kluwer Academic Publishers: Norwell, MA, 2003. (58) Pourbaix, M. Atlas of Electrochemical Equilibria in Aqueous Solutions, 2nd ed.; National Association of Corrosion Engineers: Houston, TX, 1974. (59) Loomis, C. C.; Strandberg, M. W. P. Microwave spectrum of phosphine, arsine, and stibine. Phys. Rev. 1951, 81, 798-807. (60) Shimanouchi, T. Tables of Molecular Vibrational Frequencies Consolidated, Volume I; National Bureau of Standards: Washington, DC, 1972; pp 1-160. (61) Sherman, D. M. Electronic structures of Fe3+ coordination sites in iron oxides: Applications to spectra, bonding and magnetism. Phys. Chem. Miner. 1985, 12, 161-175. (62) Ziegler, T.; Li, J. Bond energies for cationic bare metal hydrides of the first transition series: A challenge to density functional theory. Can. J. Chem. 1994, 72, 783-789. (63) Mejias, J. A.; Lago, S. Calculation of the absolute hydration enthalpy and free energy of H+ and OH-. J. Chem. Phys. 2000, 113, 7306-7316. (64) Tossel, J. A. Theoretical studies on arsenic oxide and hydroxide species in minerals and aqueous solution. Geochim. Cosmochim. Acta 1997, 61, 1613-1623. (65) Barone, V.; Cossi, M. A new definition of cavities for the computation of solvation free energies by the polarizable continuum model. J. Chem. Phys. 1997, 107, 3210-3221.

Received for review February 10, 2005. Revised manuscript received April 20, 2005. Accepted April 22, 2005. ES050271F