## Abstract

Naturally occurring mackinawite (tetragonal FeS) with incorporated transition metals is an important precursor to the formation of metal sulfides in ore deposits, but experimental results have not been sufficient to establish clear trends in the structure and stability of the transition-metal-enriched mineral. Using density functional theory with dispersion corrections, we report the first systematic examination of the relationship between composition and structure for FeS incorporating bivalent transition metals. Our method was validated by successful calculations of the structures of FeS, FeSe, FeSe_{1−}* _{x}*S

*, Fe*

_{x}_{1−}

*Me*

_{x}*Se (Me = Co, Ni, Cu), and FeNi*

_{x}*Te. Two classes of transition-metal-incorporated FeS structures then were investigated: Fe*

_{x}_{1−}

*Me*

_{x}*S (metal-substituted FeS) and FeMe*

_{x}*S (FeS intercalated by a metal at either a tetrahedral or square-pyramidal interstitial site), where Me = Co, Ni, and Cu, and*

_{x}*x*≤ 0.25. We find that incorporated transition metals can both increase and decrease lattice parameters, depending on the metal and how it is incorporated into the FeS structure. As the mole fraction of substituting metal increases, the FeS unit-cell volume decreases for Co, is nearly constant for Ni, but increases for Cu. Metal substitution changes the

*c*lattice parameter, which is sensitive to interactions between the mackinawite sheets, as much as it does the

*a*and

*b*lattice parameters. Upon intercalation, the unit-cell volume and

*c*parameter increase but the

*a*and

*b*parameters decrease. Experimental structural data are consistent with our results for metal-substituted FeS. We determined the thermodynamic stability of metal-incorporated FeS by computing the free energy involved in the metal incorporation reactions as a function of chemical potential of sulfur. The thermodynamic results lead to the general conclusions that metal incorporation into mackinawite most likely occurs via substitution, which may be important to influence phase transformation pathways of mackinawite.

## Introduction

Transition-metal-enriched mackinawite (tetragonal FeS) is an excellent system to gain insight into the effects of metal incorporation on sulfide mineral phase transformations and to understand the formation processes of metal sulfide ore deposits (Blain 1978; Lennie and Vaughan 1996; Takeno and Clark 1967; Zôka et al. 1972; Zavašnik et al. 2014). Structurally, mackinawite comprises edge-sharing sheets of FeS_{4} tetrahedra arranged on a tetragonal lattice (*P*4*/nmm* symmetry, with *a* = *b* = 3.67 Å, *c* = 5.03 Å) and stacked along the *c* direction, with stabilization through van der Waals (vdW) interactions (Lennie et al. 1995; Rickard et al. 2006). The extant literature suggests that incorporated transition metals can occupy three structurally distinct sites in mackinawite (Fig. 1): Fe substitution (Sub) and intercalation between FeS_{4} tetrahedral sheets, either at a tetrahedral interstitial site (I-td) or a square-pyramidal interstitial site (I-sp), also described as an “octahedral hole” by Ward (1970).

Vaughan (1969) analyzed the structure of nickelian mackinawite (18.7 wt% Ni) from a nickel ore deposit, reporting lattice parameters slightly smaller [by 0.03 Å along the *c* axis and ≤0.006 Å along the *a* (= *b*) axis] than those of pure mackinawite. Clark (1970) inferred from micro-indentation hardness measurements that cobaltian mackinawite (16.5 wt% Co) from a copper-cobalt ore deposit had lattice parameters even smaller than those of nickelian mackinawite. The decreases in unit-cell size were attributed by Clark (1970) to substitution of Fe by Co or Ni (i.e., metal incorporation at the Sub site). In response, Vaughan (1970) cautiously proposed that metals could also intercalate between sheets (i.e., incorporation at the I-td or I-sp site). This proposal has been widely cited in the literature (Morse and Rickard 2004; Muñoz-Santiburcio et al. 2013; Mullet et al. 2002; Watson et al. 1995; Wolthers et al. 2003). Zavašnik et al. (2014) observed a significant increase in the unit-cell size (*a* = 3.76 Å, *c* = 5.43 Å) for amorphous FeS nanoparticles reacted with Cu (Cu:Fe ratio ≈ 0.1). They conjectured that Cu intercalates into mackinawite leading to the cell expansion. However, no direct evidence supporting transition metal intercalation has yet been reported.

In a systematic approach to examine the structural effects of both substitution and intercalation for three important transition metals, Co, Ni, and Cu, we studied metal-incorporated FeS using density functional theory (DFT). Quantum mechanical geometry-optimization based on DFT with periodic boundary conditions provides reliable and detailed information on solid-phase structures (Milman et al. 2000; Payne et al. 1992) and thus is a useful method to investigate mineral structures for which significant ambiguities exist in the interpretation of experimental data, or where various samples are difficult to obtain for experimental study. Conventional semilocal DFT, such as DFT under the generalized gradient approximation (GGA), does not take weak dispersion forces into account. Therefore, when vdW interactions are a significant component of the total energy of a solid, as they are for mackinawite, the errors in calculating structural parameters can also be significant (e.g., the *c* lattice parameter). In our study, we used DFT combined with a tested semi-empirical correction to account for-dispersive vdW interactions, the DFT–D method (Grimme 2006). This method offers both accuracy and a more manageable computational expense than ab initio wave function methods or quantum Monte Carlo simulations (McNellis et al. 2009; Tunega et al. 2012). We validated the DFT–D approach by successful geometry optimization of FeSe_{1−}* _{x}*S

*(anion-substituted FeSe), Fe*

_{x}_{1−}

*Me*

_{x}*Se (Me = Co, Ni, and Cu) (cation-substituted FeSe), and FeNi*

_{x}_{0.125}Te (Ni-intercalated FeTe), the structures of which are relatively well characterized by experiment. Both FeSe and FeTe are iron chalcogenides that share structural, electronic, and magnetic properties with FeS (Kwon et al. 2011; Lennie et al. 1995; Mizuguchi and Takano 2010). For a range of incorporated metal content, we examined how the structure of FeS changes when a metal is substituted at the Fe site or intercalated into the FeS interlayer. Specifically, we examined the structural effects of metal incorporation for Fe

_{1−}

*Me*

_{x}*S (metal-substituted FeS) and FeMe*

_{x}*S (metal-intercalated FeS, I-td or I-sp site), where Me = Co, Ni, and Cu and 0 ≤*

_{x}*x*≤ 0.25. Formation energies of the metal-incorporated FeS were compared to determine which incorporation process is a most favorable for a given transition metal.

## Computational details

The initial structures of transition-metal incorporated FeSe and FeS were created based on geometry optimized FeSe and FeS unit cells. For substitution calculations, one Fe in an FeSe or FeS supercell was replaced by a transition metal (Me) to give Fe_{1−}* _{x}*Me

*Se or Fe*

_{x}_{1−}

*Me*

_{x}*S, with the mole fraction*

_{x}*x*in the range 0.0625 ≤

*x*≤ 0.25. The structural formulas were specifically: Fe

_{15}Me(Se,S)

_{16}(

*x*= 0.0625), Fe

_{7}Me(Se,S)

_{8}(

*x*= 0.125) and Fe

_{3}Me(Se,S)

_{4}(

*x*= 0.25), where Me is Co, Ni, or Cu. For intercalation calculations (FeMe

*Se or FeMe*

_{x}*S), the structural formulas were: Fe*

_{x}_{16}Me(Se,S)

_{16}(

*x*= 0.0625), Fe

_{8}Me(Se, S)

_{8}(

*x*= 0.125), and Fe

_{4}Me(Se,S)

_{4}(

*x*= 0.25).

All DFT calculations were carried out using the CASTEP code (Clark et al. 2005), which implements DFT with periodic boundaries and a planewave basis set. We used ultrasoft pseudopotentials (Vanderbilt 1990) to describe the strong Coulomb potentials between atomic nuclei and core electrons. Our Fe pseudopotential treats both 3*s* and 3*p* states as the valence state, with the valence electron configuration being 3*s*^{2}3*p*^{6}3*d*^{6}4*s*^{1}^{.}^{75} because Fe pseudopotentials treating 3*p* states as core electrons do not adequately reproduce the magnetic ground state ordering energetics of many layer-type Fe chalcogenides or pnictides (Kwon et al. 2011; Mazin et al. 2008). The valence electron configurations for the S, Se, and Te pseudopotentials are 3*s*^{2}3*p*^{4}, 4*s*^{2}4*p*^{4}, and 5*s*^{2}5*p*^{4}, respectively, while the Co, Ni, and Cu pseudopotentials have 3*d*^{7}4*s*^{1.95}4*p*^{0.05}, 3*d*^{8}4*s*^{2}, and 3*d*^{7}4*s*^{0.5}4*p*^{0.001} valence electron configurations, respectively. All calculations were performed under the spin-polarized general gradient approximation for electron correlation using the Perdew, Burke, and Ernzerhof functional (Perdew et al. 1996). Because of the metallic character of these materials, the GGA+U method is not required to describe itinerant *d* electrons (Ferber et al. 2010) and, therefore, no Hubbard U was used in our calculations.

In the DFT–D approach, the total energy (*E*_{tot}) is calculated by addition of an interatomic pairwise *C*_{6}*R*^{−6} term to the DFT energy (*E*_{DFT})

where *R** _{AB}* is the separation between atoms

*A*and

*B*,

*C*

_{6}

*is the corresponding*

_{AB}*C*

_{6}coefficient, and

*R*

^{0}

*and*

_{A}*R*

^{0}

*are vdW radii (Grimme 2011). The interaction is damped at short range by the*

_{B}*f*

_{damp}function. The

*C*

_{6}

*parameters can be obtained either by DFT calculations of atomic ionization potentials and dipole polarizabilities using the London formula for dispersion [the G06 scheme (Grimme 2006)] or by time-dependent DFT calculations using the electronic-density-based atomic volume for each atom [the TS scheme (Tkatchenko and Scheffler 2009)]. We used the G06 scheme for all calculations.*

_{AB}A planewave basis set was expanded to a kinetic energy cutoff of 400 eV. The cutoff energy for the augmentation-charge density was set to 1600 eV in geometry optimizations and 6400 eV in phonon calculations. The primitive Brillouin zone was sampled with a 14 × 14 × 11 point grid in *k* space (Monkhorst and Pack 1976) for FeS or FeSe unit cells, and Gaussian broadening of 0.01 eV was applied to partially occupied bands. Proportionally reduced grids were used for supercells. Using this selection of the energy cutoffs and *k*-point grid, the atomic force converged to within 0.01 eV/Å and the total energy converged to within 0.0001 eV. The precision of our geometry optimization method was estimated to be significantly better than 0.001 Å for the *a* parameter and 0.005 Å for the *c* parameter of transition-metal-incorporated iron chalcogenides. Magnetic ordering among Fe moments was checkerboard antiferromagnetic in the sheet (i.e., each Fe was surrounded by Fe having opposite spin), with the initial magnetic moment of the transition metal being in the same direction as the Fe that was being substituted. [See Kwon et al. (2011) for details concerning Fe magnetic ordering.] The Broyden-Fletcher-Goldfarb-Shanno (BFGS) procedure was followed in the geometry optimizations with correction for any finite basis set error (Francis and Payne 1990). The residual atomic force was less than 0.01 eV**/**Å, and the root-mean-square stress was less than 0.02 GPa.

The relative stability of metal-incorporated FeS was determined by comparison of the formation energies involved in the following reactions

where *x* is the mole fraction of Me, which represents Co, Ni, or Cu: FeS = mackinawite; CoS = NiAs-type CoS (jaipurite); NiS = NiAs-type NiS; CuS = covellite. As the reference solid phases, we used metal monosulfides (MeS) rather than elemental metals to minimize incomplete error cancellation in the total energy differences between chemically dissimilar phases (Lany 2008). Takeno et al. (1982) synthesized transition-metal incorporated mackinawite by using multiple metal sulfides at high temperature. The formation energy (Δ*E** ^{f}*) for metal substituted (Sub) or intercalated (Int) FeS was calculated as:

where *F*^{phon} (*T*) is the Helmholtz free energy, and μ_{S}(*T*,*p*) is the chemical potential of sulfur gas. The relevance of *F*^{phon} (*T*) derives from the fact that phase relationships of sulfides found at the typical depths of ore deposits are relatively insensitive to pressure, and so are mostly used as geothermometers (Barton 1970). Natural mackinawite commonly occurs in contact metasomatic, pneumatolytic, and hydrothermal deposits within metamorphic rocks and ultramafic rocks (Zôka et al. 1972, and references cited therein); thus we calculated Δ*E** ^{f}* at temperatures up to 1000 K.

The Helmholtz free energy of a metal sulfide was taken as

where *E*_{tot} is the total energy at 0 K, *F*_{vib} is the purely vibrational contribution (i.e., zero point energy, vibrational energy, vibrational entropy) at temperature *T*, and *E*_{conf} is the configurational entropy contribution for metal-incorporated FeS. The value of *E*_{conf} was calculated under the assumption of ideal mixing: for Fe_{1−}* _{x}*Me

*S,*

_{x}*E*

_{conf}=

*k*

_{B}*T*[(1 −

*x*)ln(1 −

*x*) +

*x*ln

*x*], where

*k*

_{B}is the Boltzmann constant; for FeMe

*S (I-td), in which Me occupies one I-td interstitial site instead of I-sp site per formula unit,*

_{x}*E*

_{conf}= −

*k*

_{B}*T*ln2. The difference between Sub and Int in the (

*E*

_{tot}

*+F*

_{vib}) contribution to Δ

*E*

*was approximately 0.2–0.3 eV/formula unit (fu) at 1000 K, but the contribution of*

^{f}*E*

_{conf}was only about 0.03–0.06 eV/fu at 1000 K.

The vibrational contribution (*F*_{vib}) was calculated as the harmonic phonon density of states. We performed the phonon calculations using the usual finite displacement method, constructing supercells of a geometry-optimized unit cell to obtain dynamical matrix at different phonon wavevectors (Ackland et al. 1997; Parlinski et al. 1997). In a supercell, each nonequivalent atom was displaced by ±0.005 Å along the Cartesian directions and the forces on all atoms of the supercell perturbed by the atom displacement were calculated within a real space of at least 6.5 Å radius. The Fourier transform of the force constant matrix produces the dynamical matrix, whose eigenvalues (phonon frequencies) were integrated over the Brillouin zone to compute thermodynamic quantities, such as enthalpy and vibrational entropy, at different temperatures (Dove 1993). For FeS, CoS, and NiS, the supercell size was 4 × 4 × 3 (96 atoms). The supercell size for CuS, body centered cubic ferromagnetic Fe, and orthorhombic S was 4 × 4 × 1 (96 atoms), 5 × 5 × 5 (250 atoms), and 2 × 1 × 1 (256 atoms), respectively. For metal-incorporated FeS, supercells with checkerboard antiferromagnetic ordering contained 192 or 216 atoms for Sub and 204 or 243 atoms for Int.

An H_{2}S+H_{2} gas mixture is often used to control the fugacity of sulfur gas (Sack and Ebel 2006). We inferred μ_{S}(*T*,*p*) which depends on the partial pressure ratio of H_{2}S and H_{2}, under the assumption of ideal gas behavior at equilibrium with a metal sulfide (Bollinger et al. 2003; Raybaud et al. 2000)

where *h*, *s*, *ZPE*, and *E* are the enthalpy, entropy, zero point energy, and internal energy of H_{2}S_{(g)} or H_{2(g)}, respectively, and *p* is a partial pressure (*p*^{0} = 1 atm). The values of Δ*h*(*T*, *p*^{0}) and *s*(*T*, *p*^{0}) can be found in standard thermodynamic tables (e.g., http://webboook.nist.gov and http://cccbdb.nist.gov). The difference in ZPE or *E* was calculated using a 18 × 18 × 18 Å simulation box containing a single H_{2}S or H_{2} molecule.

Comparison of formation energies is meaningful only when the chemical potential varies within a bounded range (Alfonso 2010; Krishnamoorthy et al. 2013; Reuter and Scheffler 2001). The upper limit of the allowed range can be set by the chemical potential of sulfur in the solid orthorhombic state (μ_{S}^{0}). For sulfur to be stable in a metal sulfide at equilibrium, μ_{S} must be smaller than μ_{S}^{0} (i.e., μ_{S} – μ_{S}^{0} ≤ 0). The lower limit of the allowed range was determined by considering

where Δ*E*^{f}_{FeS}, *g*_{FeS}, μ^{0}_{Fe} is the formation energy of FeS, the free energy of FeS, and the chemical potential of elemental Fe, respectively. The value of μ_{Fe} should be smaller than μ^{0}_{Fe} (i.e., μ_{Fe} – μ^{0}_{Fe} ≤ 0). This condition leads to

We used Δ*E*^{f}_{FeS} = −0.68 eV computed at 0 K with reference to elemental Fe and orthorhombic elemental S. Therefore, the allowed range of μ_{S} was −0.68 eV < Δμ_{S} = μ_{S} – μ_{S}^{0} < 0 eV.

## Results and discussion

### Validation of the DFT–D method

The total energy of FeS or FeSe is compared between the DFT and the DFT–D methods in Figure 2. Without vdW dispersion corrections, a very shallow energy minimum results along with an estimated *c* parameter, which is very different from the experimental value, whereas the *c* parameter calculated using DFT–D agrees very well with experimental data. Fully geometry-optimized structural results are summarized in Table 1. Both antiferromagnetic (AFM) ordering among Fe spins and the non-magnetic (NM) state were examined, but we found that AFM ordering is more stable than the NM state for both FeS and FeSe (Kwon et al. 2011; Subedi et al. 2008). The Fe moments we calculated for FeS and FeSe were larger than experimental values, which is a common shortcoming of DFT for metallic layer type Fe chalcogenides and pnictides (Mazin and Johannes 2009). The structural parameters (i.e., lattice parameters, bond distances, and S or Se coordinate) of the AFM state differed from experiment by <1.5% for FeS and by <0.5% for FeSe, whereas those of the NM state differed by <3.8% for FeS and 3.5% for FeSe. Given these comparisons, all subsequent results described herein are based on AFM ordering.

The lattice parameters of FeSe_{1−}* _{x}*S

*calculated using DFT–D (Fig. 3) are in excellent agreement with experimental data (Mizuguchi et al. 2009), although the*

_{x}*c*parameter is slightly underestimated when the S content is large. Our DFT–D calculations also were able to reproduce two experimentally observed trends: as the S content increases in FeSe, both lattice parameters decrease, with the change in

*c*being larger than that in

*a*(=

*b*). The experimentally determined

*c*parameter at

*x*= 0.5 was larger than at

*x*= 0.4, and thus the decrease of

*c*was not linear as a function of

*x*(S content). Our calculations did not show an increase in

*c*at

*x*= 0.5, but the decreasing gradient of

*c*varied slightly depending on whether

*x*was smaller or larger than 0.5. Experimental data are not shown for 0.5 <

*x*< 1.0 in Figure 3 because only limited solid solution appears to occur between FeSe and FeS (Finck et al. 2012; Mizuguchi et al. 2009). We estimated the excess heats of FeSe

_{1−}

*S*

_{x}*formation with respect to pure FeS and FeSe (Δ*

_{x}*E*

_{ex}) from the expression:

where *E* on the right side is, respectively, the total energy of FeSe_{1−}* _{x}*S

*, FeSe, and FeS. The Δ*

_{x}*E*

_{ex}value for FeSe

_{1−}

*S*

_{x}*was nearly zero (+2 to about +6 meV/fu) for all compositions examined. Thus the existence of FeSe*

_{x}_{1−}

*S*

_{x}*in the range 0.5 <*

_{x}*x*< 1.0 cannot be ruled out based on energy considerations alone.

The effects of transition metal substitution on the structure of FeSe are more complicated than those of anion substitution. As substitution increases, *c* tends to decrease; however, *a* either increases or decreases depending the substituting metal, as shown in Figure 4. Specifically, the decrease of *c* in the case of substitution by Cu is larger than that in the case of Co or Ni; furthermore, the increase of *a* in the case of Ni is small compared to that for Cu substitution. These experimental trends for FeSe were well reproduced in our DFT–D calculations (see the solid circles in Fig. 4), although the calculated *c* parameter for Fe_{1−}* _{x}*Ni

*Se was in poor agreement with experimental data. Figure 4 also shows that metal substitution within a FeSe*

_{x}_{4}tetrahedral sheet affects not only the

*a*parameter significantly, but also the

*c*parameter, indicating that metals other than Fe that bind to Se can strongly affect the vdW forces between the sheets.

As the metal content increases, the *a* parameter decreases in the case of Co substitution, but increases in the case of Cu. In the case of Ni substitution, the lattice parameter increases very slightly as the Ni fraction varies (Fig. 4). The decrease in the lattice parameters of FeSe_{1−}* _{x}*S

*with increasing*

_{x}*x*can be understood by considering the smaller size of S relative to Se; but, for Fe

_{1−}

*Me*

_{x}*Se, comparisons of cation size alone cannot explain the trends in variation of the lattice parameters. The ionic radius of Co*

_{x}^{2+}is 0.58 Å, that of Ni

^{2+}is 0.55 Å, and that of Cu

^{2+}is 0.57 Å (Shannon 1976), all of which are smaller than that of Fe

^{2+}(0.63 Å) and so do not explain their differing structural effects. The larger metallic radius (Wells 1991) of Cu relative to Fe (1.28 vs. 1.26 Å) may explain the increase of

*a*for Fe

_{1−}

*Cu*

_{x}*Se (Vaughan 1970), but the smaller metallic radius of Ni relative to Fe (1.24 vs. 1.26 Å) cannot account for the increase in*

_{x}*a*following substitution by Ni. Thus structural trends in metal-substituted FeSe cannot be predicted by crystal chemical effects alone but require a quantitative model of bonding energetics as provided by electronic structure calculations.

We further tested our method by examining the structure of FeNi* _{x}*Te, which is a layer type Fe chalcogenide in which the interstitial sites for Ni in the interlayer are well characterized. According to a single-crystal refinement of FeNi

_{0.1}Te by Zajdel et al. (2010), Ni is located at the I-sp site rather than the I-td site. Our DFT–D calculations for FeNi

_{0.125}Te showed that the I-sp site is indeed energetically more favorable than the I-td site (see Table 2). The structural parameters also were in better agreement with experimental data when Ni was placed at the I-sp site. The

*c*lattice parameter and Z

_{Ni}for FeNi

_{0.125}Te we calculated differed from the experimental data by less than 0.2%.

### Transition metal substitution in FeS

The trend of changes in the lattice parameters of FeS following metal substitution is similar to those for FeSe, except for the case of Cu, as shown in Figure 5. The *c* parameter of Fe_{1−}* _{x}*Cu

*S increases with*

_{x}*x*, whereas that of Fe

_{1−}

*Cu*

_{x}*Se decreases with*

_{x}*x*(compare Fig. 4c and Fig. 5c). These opposing trends imply that doping the same transition metal exerts opposite pressures on the structure depending on the chalcogenide: substitution of Cu compresses the sheet structure of FeSe along the

*c*-axis, but expands the FeS sheet structure. We also examined the interlayer spacing (

*h,*the height difference of S or Se) between neighboring sheets and the thickness (

*t*) of a single tetrahedral sheet for FeSe and FeS. In Fe

_{1−}

*Cu*

_{x}*Se, the decrease in*

_{x}*c*with increasing

*x*resulted from a reduction in

*h*rather than

*t*. For example, in the case of Fe

_{0.75}Cu

_{0.25}Se,

*h*and

*t*decreased by 0.06 and 0.02 Å, respectively, as compared to pure FeSe. However, the increase in

*c*for Fe

_{1−}

*Cu*

_{x}*S we attribute to an increase in the sheet thickness (*

_{x}*t*). In Fe

_{0.75}Cu

_{0.25}S,

*h*and

*t*increased by less than 0.01 and by 0.06 Å, respectively, as compared to FeS.

When Co occupies the Fe site, the unit-cell volume of Fe_{1−}* _{x}*Co

*S monotonically decreases as a function of*

_{x}*x*(Fig. 5a), which is consistent with the observed monotonic increase of the Vickers micro-indentation hardness as a function of the Co content in FeS (Clark 1970), where it was assumed that Co substitutes for Fe and the reduction in volume can be attributed to the monotonic decrease in

*a*as a function of Co content. However, we found that the

*c*parameter did not decrease monotonically: initially it increased slightly with increasing

*x*and then decreased when

*x*> 0.125. By contrast, the volume of Fe

_{1−}

*Cu*

_{x}*S increases as a function of*

_{x}*x*because of the increase in both

*a*and

*c*, which is consistent with experimental observations for Cu-reacted amorphous FeS nanoparticles (Zavašnik et al. 2014).

In the case of substitution by Ni, the volume of the FeS unit cell remains nearly constant as a function of *x*, changing by less than 1% (Fig. 5a). This is a result of opposing trends in the variations of lattice parameters, with a greater change in *c* than in *a*. As *x* increases, *a* increases, but *c* decreases (a slight exception occurs at *x* = 0.0625). In Fe_{0.75}Ni_{0.25}S, the *a* parameter increased by 0.027 Å as compared with FeS, but *c* decreased by 0.078 Å. Vaughan (1969) also reported a much larger change in *c* (0.03 Å) than in *a* (0.006 Å) for nickelian mackinawite as compared with pure FeS.

### Transition metal intercalation into FeS

Significant expansion in the *c* parameter occurs following transition metal intercalation into the FeS interlayer, the effect being greatest with Cu and least with Ni (Fig. 6). For substitution, the change in the *a* parameter was always greater than or comparable to the change in the *c* parameter; but, for intercalation, the latter was considerably larger than the former. This trend can be characterized by the ratio *c/a*: substitution results in a smaller *c/a* ratio than intercalation. The *c/a* ratio for pure FeS was calculated as 1.35, slightly smaller than the experimental value, 5.03/3.67 = 1.37. The calculated *c/a* ratio was 1.36 for Fe_{0.75}Co_{0.25}S (substitution), and it was 1.37 (I-td site) or 1.45 (I-sp site) for FeCo_{0.25}S (intercalation). For Fe_{0.75}Ni_{0.25}S the calculated ratio was 1.32, and for FeNi_{0.25}S it was 1.39 (I-td site) or 1.46 (I-sp site). As in the Co and Ni cases, our calculated *c/a* ratio for Cu substitution was lower than the ratio for Cu intercalation: 1.34 for Fe_{0.75}Cu_{0.25}S and 1.43 (I-td site) or 1.47 (I-sp site) for FeCu_{0.25}S.

When a transition metal occupied the I-td site, the *a* parameter changed by less than 0.01 Å as compared with pure FeS (Fig. 6a). However, when the metal occupied the I-sp site, this lattice parameter decreased by as much as 0.09 Å (Fig. 6b). In the case of Co, the effect of intercalation at the I-sp site was as large as the effect of substitution (compare Fig. 5b with Fig. 6b). Because of this significant decrease in the *a* parameter, the volume of the I-sp-intercalated FeS was comparable to that of I-td-intercalated FeS, despite the much greater *c* parameter. Intercalation at the I-sp site places the metal just below an apical S, but with a slightly larger *z*-coordinate, by 0.4–0.6 Å, than the four S of the neighboring FeS sheet (see Fig. 1), resulting in a shorter distance to the apical S than for the other four S. For example, at the I-sp site in FeCu_{0.125}S, the distance between Cu and the apical S, *d*(Cu-S), was 2.14 Å, while *d*(Cu-S) was 2.68 Å for the other four S. For the I-td site, we found *d*(Cu-S) = 2.30 Å, which is similar to the same interatomic distance in chalcopyrite. The presence of a five-coordinated metal in the interlayer also resulted in a decrease in the Fe-Fe separation within the FeS sheet. In FeCu_{0.125}S, the average Fe-Fe separation was reduced from 2.61 Å in pure FeS to 2.58 Å in intercalated FeS.

### Deducing incorporation mechanisms

Our formation energy calculations indicate that substitution is a favorable mechanism when metals incorporate into mackinawite. The formation energy (Δ*E** ^{f}*) of FeMe

*S (I-td) varied linearly as a function of S chemical potential (μ*

_{x}_{S}), while Δ

*E*

*of Fe*

^{f}_{1−}

*Me*

_{x}*S (Sub) was independent of μ*

_{x}_{S}(Fig. 7). [In calculating Δ

*E*

*for intercalation, only the I-td site was considered because, in contrast to FeNi*

^{f}_{0.125}Te (Table 2), this site was more stable than the I-sp site in FeS by 50–360 meV/fu, depending on metal type and content.] We found that Sub is thermodynamically more stable than I-td for all three metals within the allowed range of μ

_{S}(−0.68 eV < Δμ

_{S}= μ

_{S}– μ

_{S}

^{0}< 0 eV). As the metal fraction increases, the stability of Sub is further enhanced, while the stability of I-td diminishes: the formation energy of Sub is more negative and of I-td, more positive. Increasing the temperature results in similar trends, but not as pronounced. The

*p*(H

_{2}S)/

*p*(H

_{2}) scales shown in Figure 7c further indicate that, at high

*p*(H

_{2}S)/

*p*(H

_{2}) (i.e., sulfidic environment), substitution is by far the more favorable mechanism of metal incorporation.

For Co or Ni incorporation, the thermodynamic stability of substitution over intercalation also agrees well with structural data, which shows that the *c* parameter for cobaltian or nickelian mackinawite is smaller than that for pure FeS (Clark 1970; Vaughan 1969). Our computations show that the *c* parameter of FeS decreases following Co or Ni substitution; however, intercalation led to a considerable increase in the *c* parameter. While Co and Ni share the trends in structural changes upon metal incorporation, the temperature dependence of relative stability differs between the metals. As temperature increases from 300 to 900 K, the value of Δ*E** ^{f}* for Co-Sub remains nearly constant or becomes slightly less negative (Figs. 7a and 7d), but that for Ni-Sub becomes noticeably more negative (Figs. 7b and 7e). This different thermodynamic behavior is due to the difference in the vibrational contribution to Δ

*E*

*(Δ*

^{f}*F*

_{vib}) between the two metals. For Fe

_{0.75}Co

_{0.25}S, Δ

*F*

_{vib}increases with temperature (≈ +50 meV/fu at 900 K) comparably to the configurational entropy contribution (

*E*

_{conf}, ≈ −40 meV/fu at 900 K). For Fe

_{0.75}Ni

_{0.25}S, however, Δ

*F*

_{vib}is negligible (<−10 meV/fu at 900 K) while

*E*

_{conf}is ≈ −40 meV/fu at 900 K.

Although Δ*E** ^{f}* for FeCu

_{0.125}S is similar to that for Fe

_{0.875}Cu

_{0.125}S near the lower limit of μ

_{S}(Fig. 7c), which corresponds to low

*p*(H

_{2}S)/

*p*(H

_{2}) (i.e., highly reducing environment), in general, Sub is more stable than I-td for Cu incorporation. This greater stability of Sub over I-td is also consistent with structural data. For Cu-reacted FeS nanoparticles, both

*a*and

*c*are larger than for pure FeS (Zavašnik et al. 2014). This increase in the lattice parameters with Cu incorporation matches only our Cu-Sub trends (Figs. 5b and 5c) and Figure 6 shows that the parameter

*a*tends to decrease when Cu is intercalated at either the I-td or I-sp site.

Metal incorporation has been postulated to increase the stability of mackinawite (Takeno and Clark 1967; Zavašnik et al. 2014). Our finding of a negative Δ*E** ^{f}* of Co-Sub and Ni-Sub supports this hypothesis, and shows further that the mechanism is substitution, not intercalation. These results also explain why, in the absence of water, it is easier to synthesize metal-incorporated mackinawite than the pure mineral (Takeno et al. 1982). On the other hand, Cu-Sub shows a negative Δ

*E*

*only above 770 K, when*

_{f}*x*= 0.125, or above 510 K, when

*x*= 0125, and it is generally less negative than for Co- or Ni-Sub. Thus we would predict that mackinawite may accommodate a smaller content of Cu than Co or Ni. Indeed, in natural mackinawite, the Cu content is typically <10 wt%, whereas the Co or Ni content is typically <20 wt% (Clark 1970; Clark and Clark 1968; Vaughan 1969; Zôka et al. 1972).

## Implications

Naturally occurring mackinawite typically contains several species of transition metal. This study is the first to examine systematically the relationships among chemical composition, structure, and thermodynamic stability of transition-metal-incorporated mackinawite. Our results show that transition metals tend to incorporate into mackinawite by substitution at the Fe sites within the FeS_{4} tetrahedral sheets accompanied by changes in bond distances. Metal substitution enhances the stability of mackinawite to a degree that depends on both the metal and temperature. In materials processing, iron sulfides are used as catalysts and catalytic activities are often controlled by doping with transition metals. Thus our findings not only help to understand the transformations of transition-metal-incorporated mackinawite into various metal sulfides in ore deposits, but may also suggest synthetic routes for developing enhanced metal-sulfide catalysts.

## Acknowledgments

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT, and Future Planning ( NRF-2013R1A1A1004657) and a 2012 Research Grant from Kangwon National University. Support for this research also was received from the University of California at Berkeley through the appointment of G.S. as a Chancellor’s Professor. Computations were performed by using resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.