11390
J. Phys. Chem. C 2007, 111, 11390-11396
Adsorption of As(OH)3 on the (001) Surface of FeS2 Pyrite: A Quantum-mechanical DFT Study Marc Blanchard,*,†,§ Kate Wright,‡ Julian D. Gale,‡ and C. Richard A. Catlow† Royal Institution of Great Britain, 21 Albemarle Street, London W1S 4BS, United Kingdom, and Nanochemistry Research Institute, Department of Applied Chemistry, Curtin UniVersity of Technology, P.O. Box U1987, Perth 6845, Western Australia ReceiVed: March 29, 2007; In Final Form: May 22, 2007
We have performed first-principles calculations based on density functional theory (DFT), to model interactions of arsenic with FeS2 pyrite. An understanding of pyrite reactivity on the atomic scale is needed to predict the availability and mobility of arsenic in reducing conditions. Modeling of surface processes, however, requires large systems involving a vacuum part under periodic conditions, which makes plane-wave DFT methods very expensive computationally. We take advantage here of the use of localized pseudo-atomic orbitals as implemented in the DFT code SIESTA to model bulk surface reactions of pyrite with arsenic. The latter method shows good agreement with the available data (theoretical and experimental) for the structure of FeS2 pyrite, the arsenic substitution in the bulk pyrite, as well as the energetics and relaxation of the (001) pyrite surface. This method has been employed to investigate the adsorption of As(OH)3 on the (001) surface and the potential arsenic segregation with respect to the same surface. Our calculations suggest that the best adsorption configuration is a bidentate complex with two Fe-O bonds, which is similar to the first surface complexes observed in arsenite adsorption experiments. Once the arsenic atom is incorporated at the surface by substitution with a sulfur atom, it will readily incorporate into the bulk during pyrite growth; no segregation of arsenic with respect to the (001) pyrite surface is expected.
Introduction Arsenic is recognized as one of the most serious inorganic contaminants in drinking water on a worldwide basis. It is therefore essential to know the processes involved in the mobility of this metalloid in order to understand and predict its behavior in aquifers.1 In this context, iron sulfides have an important role to play, with FeS2 pyrite as the most abundant of these minerals. Under reducing conditions, such as estuarine or deltaic environments, pyrite can delay the migration of arsenic. The arsenic released into solution by desorption from iron (hydr)oxides, and by their reductive dissolution, then adsorbs on iron sulfides. It is also known that pyrite can host up to about 10 wt % of arsenic2 and thus retard the transportation of arsenic, while in oxidizing conditions, pyrite’s breakdown can pose a potential threat to the environment. In addition to the release of the arsenic it contains, pyrite dissolution is responsible for the generation of acid rock drainage and the subsequent release of further toxic elements. In parallel to experimental studies, computational techniques can greatly contribute to our knowledge of the behavior and mobility of arsenic in groundwaters by investigating the mechanisms involved on the atomic scale. Arsenic-containing species have already been modeled,3,4 and more recent studies have focused on arsenic present in the crystal structure of pyrite and marcasite.5,6 Arsenic is known to substitute for sulfur atoms * Corresponding author. Tel.: +33 (0)1 44 27 98 22. Fax: +33 (0)1 44 27 37 85. E-mail:
[email protected]. † Royal Institution of Great Britain. ‡ Curtin University of Technology. § Present address: IMPMC, Universite ´ s Paris 6 et 7, CNRS, IPGP, 140 rue de Lourmel, 75015 Paris, France.
to form AsS dianion groups, and Monte Carlo simulations indicate that pyrite and marcasite can host up to about 6 wt % of As in solid solution before unmixing into pyrite (or marcasite) and arsenopyrite. The next step forward is to model the surface reactions since the mechanism of arsenic retention is still a matter of debate.7-9 To our knowledge, no theoretical data are available so far, which can partly be explained by the computational cost of the simulations that have to consider large and complex systems. From our previous studies,5,10 density functional theory (DFT) proved to be effective in describing the Fe-S-As system. The CASTEP code,11 in which the wavefunctions are expanded in plane waves, has been used so far. However, in surface calculations, where the simulation cells contain a vacuum region between repeating slabs due to the necessity of employing 3-D periodicity, the use of plane waves becomes a disadvantage with regard to the computational cost. Furthermore, as simulations begin to target surface adsorption, the number of atoms and geometric complexity increases. Hence, it becomes necessary to contemplate methods that have improved scaling with system size and a lower prefactor in order to make such studies tractable. An alternative approach to the use of plane waves is to use a real-space basis set, such as the numerical pseudo-atomic orbitals (PAOs) as implemented in the program SIESTA.12 In this paper, we compare the results given by the two DFT approaches (i.e., plane waves and PAOs) for the study of pure pyrite and arsenic incorporation within its bulk structure. Subsequently, the SIESTA method is used to investigate the relaxation of the (001) pyrite surface (i.e., the most stable surface of FeS2 pyrite) and arsenic segregation with respect to the same surface. Finally, we report the first calculations of the adsorption of the trivalent arsenite, As(OH)3, onto this surface. As(OH)3
10.1021/jp072468v CCC: $37.00 © 2007 American Chemical Society Published on Web 07/11/2007
Adsorption of As(OH)3
J. Phys. Chem. C, Vol. 111, No. 30, 2007 11391
is the most toxic arsenic species and occurs in reducing conditions where pyrite is stable. Methodology All calculations have been performed within the framework of nonlocal density functional theory using the generalized gradient approximation (GGA) of Perdew et al.13 The numerical solution of the Kohn-Sham equations was obtained through the use of the SIESTA methodology and software.12 Here, the core electrons and nuclei are represented by nonlocal pseudopotentials of the Kleinmann-Bylander form, and the construction was performed according to the modified scheme of Troullier and Martins14 with inclusion of the relativistic effects for iron, sulfur, and arsenic. Small core pseudopotentials were selected for the metals, with the electronic configurations used in the pseudopotential generation being S (3s2,3p3.5,3d0.5), Fe (3p6,4s2, 3d6), As (3d10,4s2,3p3), O (2s2,2p4), and H (1s1). For Fe, partial core corrections15 were also included to ensure a high degree of reproduction of the all electron results for electron promotion and ionization processes. Here, the radius was set in order that the partial core charge reproduces the maximum in density due to the 3s electrons, which demonstrate significant overlap with the 3p states. The valence electron wavefunction is expanded as a linear combination of strictly localized pseudo-atomic orbitals (i.e., these are numerical tabulations of the solutions of the pseudized atomic problem within a confinement potential). In the present work, the soft confinement scheme of Junquera et al.16 has been employed to avoid the discontinuity in the basis function first derivative at the radial cutoff. The basis set employed was of double-ζ plus polarization quality for sulfur, arsenic, oxygen, and hydrogen and double-ζ for iron (i.e., no f functions were included) except for the 4p orbitals, where only a single ζ was found to be required. Again, following the work of Junquera et al.,16 the confinement and orbital parameters were variationally optimized for relevant reference systems. For Fe and S, the parameters were optimized with respect to the pyrite phase of FeS2, while for remaining elements, the basis set was generated according to the default approach with parameters of 0.02 Ry for the radial confinement and a split-norm value of 0.15. The evaluation of Hartree and exchange-correlation potentials was conducted on an auxiliary basis set consisting of a uniform real-space grid truncated by an equivalent kinetic energy cutoff. The energy cutoff, as well as the dimension of the MonkhorstPack grid17 used for the Brillouin zone integration, have been checked to yield sufficient numerical convergence, and the chosen values are given in the Results section. The incorporation of arsenic atoms into the diamagnetic FeS2 pyrite involves the presence of unpaired electrons, and the adsorption of arsenite may possibly generate charge transfer or spin density at the surface of pyrite. Therefore, spin-polarized calculations were performed. Results Pyrite Bulk Structure. FeS2 pyrite crystallizes in cubic symmetry (space group Pa3) with four formula units per unit cell. Each iron atom is coordinated by six sulfurs in a slightly distorted octahedron, and each sulfur atom is bonded to three Fe atoms and its S2 dianion group (Figure 1). All iron and sulfur sites are equivalent. The structure is fully defined by two parameters: the unit cell length and the sulfur fractional coordinate. For these calculations, a Monkhorst-Pack k-point grid of 4 × 4 × 4 has been used, while a mesh cutoff of 200 Ry was chosen for the auxiliary basis expansion of the electron
Figure 1. Unit cell of FeS2 pyrite. Iron atoms are darker than sulfur atoms.
TABLE 1: Measured and Calculated Structural Parameters: Cell Parameter, a0, in Angstroms, and Sulfur Coordinate, u, in Fractional Unitsa a0 u
SIESTA
CASTEP
exptl.
5.449 0.382
5.406 0.384
5.418 0.385
a CASTEP results are from Blanchard et al.,5 and experimental data are from Will et al.18
density; higher cutoff values yield similar results. The optimized values of the structural parameters are listed in Table 1 and compared to the previous theoretical data and experimental measurements. Both DFT implementations (using plane waves and PAOs) provide a good description of the pyrite structure with a difference in volume of only 1% compared to the experimental structure, though there is an important qualitative difference between the calculations with the plane-wave results underestimating the volume, while the PAO approach overestimates the same quantity. The electronic structure, which is described in more detail elsewhere,10 has also been checked. The calculations reproduce the semiconducting properties of pyrite, giving an indirect band gap of 0.8 eV that compares well with the experimental value of 0.9 eV.19 Given the discrepancy between the results of calculations performed within the SIESTA methodology and those previously reported with a plane-wave basis set, it is important to explore likely sources of difference that could explain this observation. To this end, an extensive series of calculations have been performed with SIESTA to explore the influence of all computational conditions on the optimized structure of pyrite, as discussed below. The first obvious difference between the two calculations is the form of the GGA employed, with the SIESTA calculations utilizing the PBE functional while the plane-wave result corresponds to the revised form (rPBE). As can be seen from Table 2, the change from PBE to rPBE, when applied consistently to the pseudopotential and bulk calculation, leads to an increase in cell parameter by 0.0473 Å and is clearly an important factor. However, since the change to the rPBE
11392 J. Phys. Chem. C, Vol. 111, No. 30, 2007
Blanchard et al.
TABLE 2: Comparison of the Optimized Lattice Parameter for Pyrite as a Function of Calculation Parameters within the SIESTA Methodologya pseudopotential functional core size PBE rPBE rPBE rPBE rPBE rPBE rPBE rPBE rPBE rPBE rPBE rPBE
SC SC SC LC LC LC LC LC LC LC LC LC
Fe 3d core radius (Bohr) 2.0 2.0 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2
basis set OPT OPT OPT OPT DZP TZ2P QZ3P PZ4P PZ4Pf QZ3P QZ3P PZ4Pf3
TABLE 3: Calculated Unit Cell Parameter (a0) and Interatomic Distances for Pure FeS2 Pyrite and the Arsenian Pyrite Containing AsS Groups, in Angstromsa SIESTA
energy optimized shift cell (Ry) parameter (Å)
0.01 0.01 0.01 0.01 0.01 0.005 0.0025 0.0025
5.4495 5.4968 5.4891 5.4877 5.4796 5.4652 5.4619 5.4614 5.4581 5.4746 5.4820 5.4705
a The following abbreviations are used below: SC ) small core; LC ) large core; OPT ) the variationally optimized basis set used in the present study; DZP ) double-ζ polarized; TZ2P ) triple-ζ double polarized; QZ3P ) quadruple-ζ triple polarized; PZ4P ) pentuple-ζ quadruple polarized; PZ4Pf ) pentuple-ζ quadruple polarized with a single f polarization function on Fe; PZ4Pf3 ) pentuple-ζ quadruple polarized with three f polarization functions on Fe.
functional leads to an increase in cell parameter, this further exacerbates the discrepancy with the previous plane-wave results. Hereafter, we discuss all results obtained with the rPBE functional for consistency with the previously reported planewave calculations. Pseudopotentials can also be a source of error in the reliability in calculations, and so the likely errors associated with this factor have been assessed. In the case of pyrite, iron is the most challenging species for the production of a transferable pseudopotential. In the present work, both a small-core form that explicitly includes the 3p electrons and partial-core corrections are included. Through reduction of the Fe 3d core radius, to create a harder pseudopotential, a decrease of the cell parameter was achieved of approximately 0.0076 Å. When this smaller core radius was adopted for the Fe 3d projector, the difference between the results for the small- and large-core pseudopotentials was negligible when including appropriate partial-core corrections. This result is consistent with the dominant effect of the 3p electrons being through the nonlocality of the exchange-correlation potential, rather than through any change of occupancy. The remaining possible reason for the significant difference in structural predictions for pyrite between the two numerical implementations is the quality of the basis set. Within the SIESTA localized orbital method, there are two key factors that control the basis set quality, namely, the number of effective ζ’s (divisions of the pseudo-atomic orbitals into components of different radial extent) and the confinement radius chosen. Here, we have explored increasing the radial variational degrees of freedom from double-ζ to pentuple-ζ, while splitting the polarization functions to have one less ζ than the other basis functions. A systematic decrease in cell parameter is observed with this progression. However, the effect is largely saturated beyond four radial components with the change from quadrupleto pentuple-ζ being less than 0.0005 Å. The addition of f polarization functions to iron leads to a further decrease in cell parameter. In contrast, decreasing the energy shift associated with radial confinement asymptotically toward the unconfined result leads to a more significant increase in lattice parameter, as shown in Table 2. On the basis of extrapolation of the radial confinement, as well as increasing the basis set variational freedom toward the
FeS2 a0 S-S As-S Fe-S Fe-As a
5.449 2.222 2.272
supercell with 2 AsS groups 5.466 2.222-2.248 2.286 2.243-2.291 2.353
CASTEP FeS2 5.406 2.188 2.256
supercell with 2 AsS groups 5.421 2.187-2.216 2.287 2.225-2.274 2.324
CASTEP results are from Blanchard et al.5
limit, we estimate that the fully converged rPBE unit cell parameter for pyrite would be close to 5.472 Å. Given that most experimental determinations of the cell parameter are close to 5.418 Å, depending on the degree of nonstoichiometry, this would represent an overestimation of 1.0%, in line with expectations for a GGA. This value is also similar to the 5.490 Å obtained by von Oertzen et al.20 based on a PBE calculation with an all-electron Gaussian basis set. Muscat et al.21 have also noted the anomalously small cell parameter yielded by a CASTEP plane-wave calculation based on ultra-soft pseudopotentials and attribute this to a large core radius. When repeating the calculations with a norm-conserving TroullierMartins pseudopotential, the form used in the present study, they found a lattice parameter of 5.478 Å, consistent with our best estimate of the converged value at the basis set limit. Hence, having quantified the numerical errors associated with the SIESTA methodology when applied to pyrite, we can safely identify the plane-wave result as being inferior with respect to a well-controlled rPBE/PBE calculation. The use of the double-ζ polarized bulk optimized basis sets for the subsequent calculations has also shown to be a good compromise between computational cost and achieving a well-converged numerical result. Arsenic Incorporation. Our previous study5 investigated the substitution mechanisms involved in the arsenic incorporation as a solid solution into pyrite. The solution energies of incorporation reactions under both oxidizing and reducing conditions were calculated using plane-wave DFT. The results suggest that it is more energetically favorable to substitute arsenic for sulfur rather than for iron and that the formation of AsS groups is favored compared to As2 groups. A simple way to validate the latter statement consists of comparing the total energy of two 2 × 2 × 2 supercells of the same composition, Fe32S62As2, and containing either one As2 group or two AsS groups. Such a composition is relevant to natural arsenic concentrations (
