A density functional theory study of the hydrogenation and reduction of the thio-spinel Fe₃S₄{111} surface

Abstract

The mineral greigite, Fe₃S₄, shows promising electro-reduction activity, especially towards carbon dioxide conversion to small organic molecules. We have employed density functional theory calculations with correction for the long-range dispersion forces to investigate the behavior of hydrogen on the greigite{111} surface. We have studied the adsorption, diffusion, surface reduction and associative (i.e. Volmer–Tafel mechanism) and molecular desorption of hydrogen as a function of its coverage. We found that (i) the H ad-atoms adsorb on S sites far from metallic centres in the topmost surface layer; (ii) the reduction of greigite by hydrogen is energetically unfavorable at any surface coverage; and (iii) molecular hydrogen evolution has a transition state at ∼0.5 eV above the energy of the reactants on Fe₃S₄{111}, which is very similar to the barrier found experimentally on Pt{111}. We have also determined the electrode potential under room conditions at which the H₂ evolution reaction becomes energetically barrierless.

---

1. Introduction

Sulfides of iron, the most abundant transition metal element in the Earth's crust, frequently contain Fe and S in different oxidation states, yielding various types of natural iron sulfides. Among them, greigite is isomorphic with magnetite (Fe₃O₄)¹ and has been associated as a catalyst in a number of biochemical reactions, for example those associated with the “iron–sulfur world” hypothesis for the origin of life.^2–5 Thio-spinel crystal structures, such as greigite, resemble the reactive cubane structures (Fe₄S₄) found in enzymes but linked through terminal sulfides.^6–8 Although most biomimetic catalysts lack the protection and long-range interactions of the protein backbone, they show interesting catalytic activity.⁹ Fe₃S₄ also has promising applications in spintronic devices,¹⁰ electrodes for lithium-ion batteries,¹¹ drug delivery and hydrogen storage,^12–14 and it has been used as a catalyst in electro-reductive experiments.¹⁵ Understanding the redox properties of this material is therefore crucial for optimising its performance, while mitigating, for example, its corrosion. The early oxidation stages of the Fe₃S₄ surface indicate the formation of oxide shell structures (Fe₃O₄), although the inner material remains a sulfide (Fe₃S₄).^16,17 Nevertheless, the behaviour of the greigite surface under anodic potentials remains unreported. Generally, transition metal sulfides (TMS) have gained scientific interest as catalysts, in part due to their similitude to important enzymatic catalytic centres, e.g. in the H₂–CO₂ reaction couple.^18–20 TMS such as MoS₂, NiS_x, RuS₂ and spinel-type Mo₃S₄ have been tested, concluding that the TMS redox character depends on the nature of the metal atom and the sulfur coordination environment.^21–26 More sophisticated catalysts such as chalcogels enhance the evolution of H₂ but at very low turnover, even when doping the inorganic structure with other metals such as Ni, Co or Mo.²⁷

Accurate simulations of electrochemical processes, however, are a challenging task due to the presence of a solid–solvent interface, and multiple species dynamically interacting with the surface, e.g. ions, solvent, reactive species and the flow of electrons.²⁸ We have therefore followed a simple approximation in which the reaction energies are plotted as a function of the variation in the work-function under a fixed external potential.²⁹ Thus, the electronic structure calculations are related to the surface potential by the variations of the dipole.³⁰

Based on previous research and the fact that H binds on top of low-coordinated S atoms on TMS surfaces,^31–33 we have used density functional theory with a correction for the long-range dispersion forces to focus on the Fe₃S₄{111} surface, which has been shown to be exposed in Fe₃S₄ crystallites. In particular, we have analysed the H ad-atom migration, the surface reduction by either H-incorporation or by H₂S formation, and the catalytic activity towards the hydrogen evolution reaction (HER) via the Tafel reaction.³⁴

2. Computational methods

Periodic plane-wave DFT-D calculations were carried out to study the hydrogen interaction and reactivity on the greigite Fe₃S₄{111} surface. All calculations were performed using the Vienna Ab initio Simulation Package (VASP)^35,36 with the computational details that accurately describe greigite, in full agreement with recent synchrotron radiation measurements and catalytic experiments.^15,37 Hence, the ion–electron interactions were represented by the projector-augmented wave (PAW) method³⁸ within the generalized gradient approximation (GGA) with the Perdew–Wang 91 functional³⁹ and employing the spin interpolation formula of Vosko et al.⁴⁰ All calculations included the D2 long-range dispersion correction approach by Grimme (s₆ = 0.75),⁴¹ as this correction has been shown to reproduce well the interaction of the sulfides with molecules.^32,33,42,43 The PBE functional and D3 long-range dispersion showed no difference with the setup used here.⁴⁴ The Kohn–Sham valence states were expanded in a plane-wave basis set with a cut off at 600 eV for the kinetic energy.⁴⁵ This high cut off energy ensured that no Pulay stresses occurred within the cell during relaxations. Calculations were carried out with Monkhorst–Pack⁴⁶ grids of 5 × 5 × 1 K-points, which ensured electronic and ionic convergence. The geometries of all stationary points were found by the conjugate-gradient algorithm and were considered converged when the force on each ion dropped below 0.03 eV Å⁻¹, whereas the energy threshold defining self-consistency of the electron density was set to 10⁻⁵ eV. To improve the convergence of the Brillouin-zone integrations, the partial occupancies were determined using the tetrahedron method with Blöchl correction smearing, with a set width for all calculations of 0.02 eV.

The initial magnetic moment with high spin distribution describes the Fe_A(Fe_B)₂S₄ spinel structure with both types of Fe, i.e. octahedral (B) and tetrahedral (A), in a ferrimagnetic orientation, as reported previously.^37,47 For an accurate treatment of the electron correlation in the localised d-Fe orbital, we have used the U approximation,^48,49 which improves the description of local states in highly correlated systems where standard LDA and GGA functionals fail.⁵⁰ A problem with this approximation is the rather empirical character of the choice of U parameter, a feature which also appears when using hybrid functionals, since the amount of Fock exchange is system-dependent.^51–53 Here, we have used U_eff = 1 eV for greigite, based on the lattice parameters and the band gap as compared with the available experimental data.^10,47,54

2.1 Slab model

After a full relaxation of the bulk, we prepared the Fe₃S₄ surfaces by cutting the bulk structure and creating slab models using the METADISE code.⁵⁵ These models consider periodicity and provide the different atomic layer stacking resulting in a null dipole moment perpendicular to the surface plane.^32,56 The slabs contain 56 atoms (24 Fe and 32 S) per unit cell where the atoms are placed in four layers of two Fe₃S₄ units. We added a vacuum width of 12 Å between periodic slabs, which is big enough to avoid interactions between the slab and its periodic images, leading to a surface area at each end of the slab of 81.0 Ų. The slab is thick enough to relax the two uppermost layers (four Fe₃S₄ units) keeping the bottom layers frozen to model the bulk structure. To obtain the properties of the isolated molecules, we placed them in the centre of a 15 × 16 × 17 ų box to avoid lateral interactions and using the same convergence criteria as in the iron sulfide systems.

The calculations of the different minima provide the necessary information to estimate the thermodynamics of the reaction. A saddle point (TS) in the potential energy surface links these minima and determines the kinetics of the process. We have identified these particular points by means of the dimer method.^57,58 We further confirmed the TS by a vibrational frequency analysis, where only one imaginary frequency was obtained corresponding to the reaction coordinate. Afterwards, the dimer images were relaxed to the neighbouring local minima, where in a successful search, one of the images will minimise into the initial state, and the other will become the final state.

We have calculated the binding energies of the adsorbates (E_B) per number of ½ H₂ molecules on the Fe₃S₄ surface from eqn (1),

where E_n·H:Fe3S4 is the total energy of n hydrogens interacting with the slab, E_Fe3S4 is the energy of the naked slab and E_H2 is the energy of the isolated H₂ molecule in vacuo. The activation energy (E_A) of a certain process is the energy required to pass over the potential energy barrier characteristic of the transition state. We defined this barrier energy as the difference between the initial state and the transition state for the forward process. In the same way, we defined the reaction energy (E_R) as the total energy difference between the final state (products) and the initial state (reactants).

We carried out charge analyses using AIMS,^59,60 as implemented in the Henkelman algorithm.⁶¹ Within this approach, the electron (and spin) density associated with each atom is integrated over the Bader volume of the atom in question. We can use this method to monitor changes in charges, for example as an effect of surface adsorption.

3. Results and discussion

Synthetic greigite particles mainly express the {001} and {111} surfaces,^15,62 which contain two types of terminal Fe atoms derived from the tetrahedral (Fe_A) and the octahedral (Fe_B) bulk positions.⁶³ Upon relaxation of these surfaces, cations tend to move inwards into the bulk and anions outwards, i.e. the top Fe_A and Fe_B move underneath the sulfur layer.³² Here, we have focused on the {111} surface as it has shown interesting redox properties towards the reduction of CO₂.^15,42,64 This surface exposes in the uppermost atomic layer an Fe_B and two Fe_A surrounded by eight S atoms per unit cell. At the surface, both Fe_A and Fe_B atoms have a formal charge of 3+.

3.1 Hydrogen adsorption and diffusion on Fe₃S₄{111}

We have approached a single H atom on non-equivalent adsorption sites on the Fe₃S₄{111} surface to study its binding characteristics and found that H prefers to adsorb on S far from any surface Fe_B. Although the S–H distances are very similar on the different adsorption sites (d_S–H ∼ 1.36 Å), the most exothermic location (E_B = −0.4 eV) corresponds to sites 3 and 4 in Fig. 1, whereas site 1 on top of Fe_B leads to the least favourable adsorption. There, the distance Fe_B–H is 0.14 Å longer than the main d_S–H, and its E_B is endothermic by +0.7 eV. We did not find any minimum for H adsorbed on top of Fe_A. The S sites directly connected to Fe_A, i.e. 3 and 4, present a stronger basic character than the average, which is due to the distortion created by the Fe_A reconstruction upon relaxation of the naked surface that slightly modifies the local electronic structure. In fact, the H binding energy has a linear relation with the S electronic structure, see Fig. S2 (ESI†). Hence, the surface basicity is associated with the external electrons of the S-bands, and the more negatively charged is S, the weaker is the bond with the H atoms.

Fig. 1 Top-view representation of non-equivalent H adsorption sites on the Fe₃S₄{111} surface. Grey balls and sticks denote Fe cations, dark-yellow S anions and white H atoms. Numbers and arrows indicate the H position and migration direction in Table 1.

We have analysed the S–H bond by charge analysis and charge density difference flux. The latter, depicted in Fig. S3 (ESI†), shows charge density accumulation between S and H, where H couples its lone electron with one from sulfur leading to a bond with covalent character. This result agrees with the charge analysis, which indicates paired electrons with no spin density and a charge transfer of 0.6 e from the sulfur to the bond. In contrast, the low stability of the Fe_B–H interaction is characterised by a charge transfer of 0.4 e to the H, i.e. half the charge of the Fe_B.

We have investigated the migration of a single H ad-atom over the Fe₃S₄{111} surface using both thermodynamic and kinetic arguments. Although the H position on top of the S atoms is the energetically most favourable configuration, we have also considered the Fe_B site to provide a more general picture. The arrows in Fig. 1 schematically represent the H migration path over the surface for the processes summarised in Table 1. The reaction energies in this table show that H migrates easily to a position far from the Fe_B. Although at low temperatures process A is more likely than B due to its lower energy barrier (by 0.12 eV), B leads to a more stable configuration (by 0.16 eV), which may drive the diffusion process. Once on the S sites, H can migrate to another S by overcoming an average energy barrier of 0.7 eV via a transition state on top of a sub-surface Fe_B.

Table 1 Barrier (E_A) and reaction energies (E_R) for H movement over Fe₃S₄{111} between initial and final sites in Fig. 1. Note that binding energies (E_B) are for products and the energy values are not corrected by ZPE. The diffusion coefficients (D) are calculated at a temperature of 300 K

Process	Initial site		Final site	E A (eV)	E R (eV)	E B (eV)	D (cm^2 s^−1)
A	1	→	2	0.01	−0.80	−0.07	3.29· × 10^−3
B	1	→	5	0.13	−0.96	−0.23	1.41 × 10^−5
C	2	→	3	0.83	−0.32	−0.39	4.41 × 10^−17
D	2	→	7	0.39	0.14	0.07	2.58 × 10^−10
E	3	→	4	0.93	−0.02	−0.41	8.87 × 10^−19
F	3	→	7	0.78	0.46	0.07	2.16 × 10^−16
G	5	→	6	0.54	0.29	0.05	1.67 × 10^−12

---

We have derived the H diffusion coefficient (D) on the Fe₃S₄{111} surface by evaluating the changes in the vibrational contribution to the energies within the limits of validity of the transition state theory.⁶⁵ We used eqn (2) considering that a single displacement occurred when an H atom moved from one site to the next.

where ΔE_ZPE is the activation energy corrected by the vibrational contribution at the ground state and the saddle point, and nd²/2α is the Einstein relation for a random walk; n is the number of jumps to reach the next site (in this case n = 1), d is the distance between the initial and final adsorption sites, which ranges between 1.66 and 3.27 Å depending on the process, and α indicates the dimensions of the jump, which, according to the symmetry of the surface, is one.⁶⁶ The results are in line with the adsorption energies and local electronic structures and agree with the variation in acid–base character of the top S atoms on the Fe₃S₄{111}, discussed above. In summary, the H movements over the surface avoid the S-sites directly connected to Fe_B, which may be available for the adsorption of other species such as CO₂.

The resulting diffusion coefficients indicate a hindered movement of the H ad-atom in comparison to values on transition metals,⁶⁷e.g. the diffusion coefficients range from 3.35 × 10⁻⁴ to 2.2 × 10⁻³ cm² s⁻¹ on body cubic centered Fe at 233 K.⁶⁸ Nevertheless, our values are in the range of diffusion on oxide surfaces with variations in acid–base character,^69,70 showing the localization of the electron and the varying character of the sulfur sites according to their local electronic structure.

3.2 Hydrogen coverage on Fe₃S₄{111}

Actual surfaces are dynamic systems whose properties and processes depend on the concentration and types of adsorbents. The H coverage (θ_H) is, therefore, a factor to consider when developing a general picture of the surface character. We defined θ_H as the number of H atoms per S site on the surface but note that H may also adsorb on Fe_B only under strong conditions, i.e. high hydrogen chemical potentials. We increased θ_H by adding pairs of H atoms on the surface and compared the energies with a reference system: the naked slab and a number of isolated H₂ molecules in the gas phase. Fig. 2 represents the E_B per atom as a function of θ_H starting from 1 up to 8 H atoms per cell, i.e. θ_H = 0.125 to 1 ML. We tested the effect of increasing up to 10 H by placing the extra H on Fe_A and Fe_B, but this resulted in an endothermic step by 0.44 eV. A similar energy increment was found at θ_H > 1 ML on Pt{111}, which has been interpreted as a strong and repulsive interaction between H ad-atoms.⁷⁰ While neighbouring H ad-atoms are interacting by less than 0.1 eV on Fe₃S₄{111}, we related this step of 0.44 eV with the surface electronic structure. The increment of θ_H polarises the surface, shown by the average sulfur charge on the surface, see Fig. S4 (ESI†), which shows how the S sites become electron-deficient at high θ_H, while decreasing the overall stability of the system. Hence, at high θ_H, the S charge density localises along the S–H bonds generating a dipole perpendicular to the surface which drives some reconstruction.⁵⁶

Fig. 2 Binding energies per atom (E_B) as a function of H coverage (θ_H) on Fe₃S₄{111}. Solid line represents the logarithmic regression: E_B = −0.096 + 0.4623·ln(θ_H + 0.3902), R = 0.99.

As the number of H ad-atoms increases, the number of available adsorption sites decreases and the lateral interaction between hydrogens could become considerable. However, we have calculated the H–H interaction to be very weak (<0.1 eV), which is in good agreement with experimental values in another catalyst.⁷¹ In general, we have noted that the presence of –SH groups modifies the surface electronic structure in a way that enhances the basicity of nearby sulfur sites. The presence of polar adsorbates and solvents changes their orientation according to the generated dipole, i.e. creating a double layer model.^28,32,72

3.3 Hydrogen evolution on Fe₃S₄{111}

As the hydrogen coverage increases on the surface, so does the probability of their recombination into molecular hydrogen. We have studied the H₂ associative and molecular desorption processes from ad-atoms on S, Fe_A and Fe_B. In the first mode, the two H ad-atoms recombine as they leave the surface – the inverse process to dissociative adsorption. Commonly, H₂ is adsorbed dissociatively, and its desorption may follow a similar trajectory.⁷³ In the latter, the H₂ molecule is formed on the surface prior to desorption – the inverse process to molecular adsorption and dissociation. We have summarised in Table 2, and schematically represented in Fig. S5 (ESI†), the H₂ formation on the surface before desorption, i.e. molecular desorption, and the association of both H ad-atoms along the desorption process, i.e. associative desorption.

Table 2 Activation (E_A) and reaction (E_R) energies for the molecular hydrogen formation and desorption processes on Fe₃S₄{111} with respect to the naked slab and the isolated H₂ molecule. MD and AD indicate respectively molecular and associative desorption. S′ indicates a different S adsorption site

Process		Site	E A (eV)	E R (eV)
1	MD	FeB	0.09	−0.79
2	AD	FeB–S	1.36	−0.63
3	AD	FeA–S	0.17	−0.98
4	AD	S–S′	1.52	0.17
5	AD	S–S′	0.46	−0.48
6	AD	S–S′	1.41	0.63
7	AD	S–S′	2.09	0.63

---

Molecular desorption takes place only when both H atoms are adsorbed on Fe_B, e.g. process 1 in Table 2, which is very unlikely due to the instability of H adsorption on Fe_B. Examples of thermodynamically favorable associative desorption are: (i) two hydrogens on S and Fe_A and (ii) two hydrogens on neighbouring sulfurs, listed as processes 3 and 5 in Table 2, respectively. This latter process is exothermic by 0.48 eV, and its energy barrier (E_A) is only 0.5 eV, which is similar to the barrier found experimentally on Pt{111}.⁷⁴ We will focus on system 5, as it has the lowest activation barrier for the association of two H sited on the preferred S sites.

Under low H chemical potentials, i.e. low coverages, the H ad-atoms are located on top of S sites surrounding Fe_A, i.e. 3 and 4 in Fig. 1. While the hydrogen associative desorption on these sites is endothermic by 0.63 eV (process 7 in Table 2) an increase in the hydrogen coverage (θ_H) may lead to the more feasible process 5 (see Fig. S5, ESI†).

We have approximated our DFT energies (E) to the Gibbs free energies (G) by adding the zero-point energy (ZPE) and including the energy contribution from the configurational (S_conf in eqn (3)) and the vibrational (S_vib in eqn (4)) entropies. The differential configurational entropy, dS_conf/dθ_H, has been found to compare well with Monte Carlo simulations on hexagonal surfaces, such as the S site distributions,^75–77 whereas the vibrational entropy is commonly described by a Boltzmann distribution.⁷⁸ For the gas phase product (H₂), we have introduced the energy variation with temperature using the specific heat (C_p) and the rotational, translational, and vibrational energy contributions from the entropy (S) derived from the DFT partition functions. We have considered the electronic partition function as the ground state. This approach led to a relative error between the resulting H₂ free energy (G) and the thermochemistry data from ref. 79, which is smaller than 0.09% in a range of temperatures from 250 to 700 K, see Fig. S6 (ESI†).

S_conf = −k_B(θ_H·ln(θ_H) + (1 − θ_H)·ln(1 − θ_H)) (3)

As general trends, the thermodynamic saturation of the Fe₃S₄{111} takes place at the coverage of ∼0.35 ML as the adsorption of H beyond this θ_H is endothermic, see Fig. 3A. Next, we studied process 5 in Table 2 to consider the evolution of H₂ as a function of the coverage; this process initiates from hydrogens adsorbed on favorable S-sites and has the lowest activation energy. We found that the evolution of H₂ is exothermic before reaching θ_H = 0.35 ML–the H₂ evolution process is lower in energy than the H adsorption–although it is hindered by the large activation energy (ΔG^TS). According to the Tafel reaction mechanism for the evolution of H₂,⁷⁴ the activation and stability (ΔG_B) are almost unaffected by the presence of water, electric potentials, or electric fields, which are, however, important to simulate an electrochemical environment.²⁸ Hence, θ_H may go beyond 0.35 ML, e.g. as a result of electric potentials, until the activation energy is small enough to allow the H₂ evolution. It is at θ_H = 1 ML that the barrier is smaller than 0.2 eV (the difference between ΔG^TS and ΔG_B) leading to an exothermic associative desorption of H₂ by 0.88 eV (the difference between ΔG_R and ΔG_B).

Fig. 3 Free energies per H atom for the H binding (black circles), desorption (blue squares) and barrier (ΔG^TS) as a function of (A) hydrogen coverage (θ_H) and (B) electrode overpotentials (U) vs. the normal hydrogen electrode at 300 K and an H₂ pressure of 1 atm. Note in the inset that asterisks denote adsorbed and free sites.

As described in the hydrogen coverage section, an increment of θ_H polarises the surface, which affects the surface work function (ϕ)–the minimum work required to extract an electron from the top of the valence band of the uncharged surface and place it in a vacuum. We have related the changes in the work function upon evolution of H₂ (final states) to an electrode potential at zero-charge (U = ϕ − ϕ_NHE) relative to a normal hydrogen electrode (NHE) of 4.5 V,^28,30 which is within the range found in the literature, 4.4–4.9 V.⁷⁸ We have derived a Tafel transfer coefficient of 0.62, which expresses the fraction of the overpotential that helps to lower the activation energy of the homolytic H₂ formation and the potential affecting the current density. Fig. 3B indicates how the H ad-atoms on the Fe₃S₄{111} surface accumulate under cathodic potentials–solid line in Fig. 3. Although the H₂ evolution becomes thermodynamically favourable below 0.35 V (vs. NHE)–the H₂ evolution is lower in energy than adsorption–, it is kinetically hindered–big difference between ΔG^TS and ΔG_B–until U = −0.4 V (vs. NHE), where the process becomes barrierless. Indeed, bubbling has been observed at potentials of 0.6 V (vs. NHE) at room pressure and temperature.¹⁵

3.4 Fe₃S₄{111} surface reduction

Sequestration of H. The small size of the H atoms may allow them to penetrate inside the greigite structure. We have evaluated the thermodynamic viability of this process by placing one H in the centre of the cavity created between Fe_B atoms and the S in the uppermost layer. During the geometry optimisation, the intrusive H moved aside to bind a sub-surface Fe_B, see Fig. S7 (ESI†). This structure increased the surface work function in comparison with the H adsorbed on the surface. However, the charge analysis shows a negligible variation of only 0.1 e in the H atom. We have studied the same process at different surface coverages, i.e. θ_H = 0.25 ML and θ_H = 1 ML, and found that the intercalation processes remains largely endothermic compared with the H adsorbed on the surface, +0.83 and +1.62 eV respectively. These high energies indicate that H sequestration into greigite would not take place under normal conditions.

Formation of H₂S. Another way to reduce Fe₃S₄{111} is by removing sulfur from the surface, for example as H₂S. We have studied hydrogen sulfide (H₂S) formation when two H ad-atoms are adsorbed onto the same S-site. This surface configuration weakens the Fe–S bond and forms H₂S (E_R = 0.09 eV), see Fig. 4. Once formed, the molecule hovers at ∼1.6 Å above its initial position. However, the energy barrier to place both H on the same S-site is 2.35 eV under low θ_H is 2.35 eV, which prevents the process from being kinetically accessible. At a θ_H of 1 ML, the reaction energy is almost the same (E_R = 0.07 eV), but the energy barrier is even higher than at low coverage (E_A = 2.65 eV). Thus, although H₂S formation is thermodynamically possible, it is kinetically inaccessible under conditions of both low and high hydrogen coverage.

Fig. 4 Top view representation of (A) two H adsorbed and (B) H₂S on S-defective Fe₃S₄{111}. Grey balls and sticks denote Fe cations, dark-yellow S anions, and white H atoms.

4. Conclusions

We have performed a systematic study of common processes involving hydrogen on greigite. We have used DFT-based calculations with the PW91 exchange–correlation potential, including dispersion through the Grimme's D2 approach and the on-site Hubbard correction (U_eff), to evaluate H adsorption and diffusion, H₂ evolution by ad-atoms recombination, and reduction of the surface by either H-incorporation or H₂S formation on the Fe₃S₄{111}. We found that H adsorbs preferentially on surface S-sites far from the top surface Fe_B, which distorts the nearby electronic structure, making the S-sites more acidic. Thus, the H ad-atoms diffuse across the {111} surface along paths that avoid low-coordinated Fe_B. On the Fe₃S₄{111} surface, the strongest mono-atomic binding energy is −0.4 eV, which becomes weaker with coverage, following a logarithmic trend up to a θ_H of 1 ML; above this coverage, there is a trend discontinuity, not only due to the H–H interaction (<0.1 eV) but because of changes in the surface polarisation. The reduction of the Fe₃S₄{111} surface by incorporation of H into the thio-spinel structure is an endothermic process by more than 0.8 eV. H₂S formation is in thermodynamic equilibrium with two contiguous HS, but the process is kinetically prevented between 0.25 and 1 ML (E_A ∼ 2 eV). The H₂ evolution reaction takes place via associative desorption from distinctive sites, and both the barrier and the reaction energies decrease with increasing H-coverage: the reaction becomes thermodynamically favourable at coverage higher than 0.35 ML. We have also related the θ_H with the electrode potential via the work function and found that at bias of 0.4 V vs. NHE the H₂ evolution becomes thermodynamically favourable but kinetically controlled. At potentials lower than −0.45 V, the activation energy becomes negligible and the Tafel process becomes very fast and exothermic.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge the Engineering & Physical Sciences Research Council (grants EP/K035355/1 and EP/K035355/2) for funding. This work made use of the ARCHER facilities, the UK's national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc. and NAG Ltd, funded by the Office of Science and Technology through EPSRC's High-End Computing Programme and provided via our membership of the HPC Materials Chemistry Consortium (EPSRC EP/L000202). The authors also acknowledge the use of the UCL Legion High-Performance Computing Facility and the Advanced Research Computing@Cardiff (ARCCA) at Cardiff University, and associated support services, in the completion of this work. All data created during this research are openly available from the University of Cardiff Research Portal at http://doi.org/10.17035/d.2018.0053780764.

References

1. A. Roldan, D. Santos-Carballal and N. H. de Leeuw, A comparative DFT study of the mechanical and electronic properties of greigite Fe₃S₄ and magnetite Fe₃O₄, J. Chem. Phys., 2013, 138, 204712.
2. S. Mann, N. H. C. Sparks, R. B. Frankel, D. A. Bazylinski, H. W. Jannasch, D. A. Bazylinski and H. W. Jannasch, Biomineralization of ferrimagnetic greigite (Fe₃S₄) and iron pyrite (FeS₂) in a magnetotactic bacterium, Nature, 1990, 343, 258–261.
3. T. Kasama, M. Posfai, R. K. K. Chong, A. P. Finlayson, P. R. Buseck, R. B. Frankel and R. E. Dunin-Borkowski, Magnetic properties, microstructure, composition, and morphology of greigite nanocrystals in magnetotactic bacteria from electron holography and tomography, Am. Mineral., 2006, 91, 1216–1229.
4. M. Posfai, P. R. Buseck, D. A. Bazylinski and R. B. Frankel, Iron sulfides from magnetotactic bacteria: structure, composition, and phase transitions, Am. Mineral., 1998, 83, 1469–1481.
5. M. Posfai, P. R. Buseck, D. A. Bazylinski and R. B. Frankel, Reaction sequence of iron sulfide minerals in bacteria and their use as biomarkers, Science, 1998, 280, 880–883.
6. H. Dobbek, V. Svetlitchnyi, L. Gremer, R. Huber and O. Meyer, Crystal structure of a carbon monoxide dehydrogenase reveals a Ni-4Fe-5S cluster, Science, 2001, 293, 1281–1285.
7. M. Dorr, J. Kassbohrer, R. Grunert, G. Kreisel, W. A. Brand, R. A. Werner, H. Geilmann, C. Apfel, C. Robl and W. Weigand, A possible prebiotic formation of ammonia from dinitrogen on iron sulfide surfaces, Angew. Chem., Int. Ed., 2003, 42, 1540–1543.
8. U. P. Apfel and W. Weigand, Efficient activation of the greenhouse gas CO₂, Angew. Chem., Int. Ed., 2011, 50, 4262–4264.
9. B. D. Yuhas, C. Prasittichai, J. T. Hupp and M. G. Kanatzidis, Enhanced Electrocatalytic Reduction of CO₂ with Ternary Ni-Fe₄S₄ and Co-Fe₄S₄-Based Biomimetic Chalcogels, J. Am. Chem. Soc., 2011, 133, 15854–15857.
10. B. Zhang, G. A. de Wijs and R. A. de Groot, Switchable Fermi surface sheets in greigite, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 20404–20406.
11. X. Jia, Z. Chen, X. Cui, Y. Peng, X. Wang, G. Wang, F. Wei and Y. Lu, Building robust architectures of carbon and metal oxide nanocrystals toward high-performance anodes for lithium-ion batteries, ACS Nano, 2012, 6, 9911–9919.
12. F. Cao, W. Hu, L. Zhou, W. Shi, S. Song, Y. Lei, S. Wang and H. Zhang, 3D Fe₃S₄ flower-like microspheres: high-yield synthesis via a biomolecule-assisted solution approach, their electrical, magnetic and electrochemical hydrogen storage properties, Dalton Trans., 2009, 9246–9252.
13. Y.-S. Chang, S. Savitha, S. Sadhasivam, C.-K. Hsu and F.-H. Lin, Fabrication, characterization, and application of greigite nanoparticles for cancer hyperthermia, J. Colloid Interface Sci., 2011, 363, 314–319.
14. M. Feng, Y. Lu, Y. Yang, M. Zhang, Y.-J. Xu, H.-L. Gao, L. Dong, W.-P. Xu and S.-H. Yu, Bioinspired greigite magnetic nanocrystals: chemical synthesis and biomedicine applications, Sci. Rep., 2013, 3, 2994.
15. A. Roldan, N. Hollingsworth, A. Roffey, H.-U. Islam, J. B. M. Goodall, C. R. A. Catlow, J. A. Darr, W. Bras, G. Sankar, K. B. Holt, G. Hogarth and N. H. de Leeuw, Bio-inspired CO₂ conversion by iron sulfide catalysts under sustainable conditions, Chem. Commun., 2015, 51, 7501–7504.
16. S. N. A. Zakaria, N. Hollingsworth, H. U. Islam, A. Roffey, D. Santos-Carballal, A. Roldan, W. Bras, G. Sankar, G. Hogarth, K. B. Holt and N. H. de Leeuw, Insight into the Nature of Iron Sulfide Surfaces During the Electrochemical Hydrogen Evolution and CO₂ Reduction Reactions, ACS Appl. Mater. Interfaces, 2018, 10, 32078–32085.
17. D. Santos-Carballal, A. Roldan and N. H. de Leeuw, Early Oxidation Processes on the Greigite Fe₃S₄ (001) Surface by Water: A Density Functional Theory Study, J. Phys. Chem. C, 2016, 120, 8616–8629.
18. Y. Nicolet, P. Amara, J. M. Mouesca and J. C. Fontecilla-Camps, Unexpected electron transfer mechanism upon AdoMet cleavage in radical SAM proteins, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 14867–14871.
19. T. Hayashi and A. A. Stuchebrukhov, Electron tunneling in respiratory complex I, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 19157–19162.
20. H. Seino and M. Hidai, Catalytic functions of cubane-type M4S4 clusters, Chem. Sci., 2011, 2, 847–857.
21. M. Neurock and R. A. Vansanten, Theory of Carbon–Sulfur Bond Activation by Small Metal Sulfide Particles, J. Am. Chem. Soc., 1994, 116, 4427–4439.
22. F. Frechard and P. Sautet, Chemisorption of H-2 and H2S on the (100) surface of RuS₂: an ab initio theoretical study, Surf. Sci., 1997, 389, 131–146.
23. A. Travert, H. Nakamura, R. A. van Santen, S. Cristol, J. F. Paul and E. Payen, Hydrogen activation on Mo-based sulfide catalysts, a periodic DFT study, J. Am. Chem. Soc., 2002, 124, 7084–7095.
24. X. D. Wen, T. Zeng, B. T. Teng, F. Q. Zhang, Y. W. Li, H. G. Wang and H. J. Jiao, Hydrogen adsorption on a Mo(27)S(54)cluster: a density functional theory study, J. Mol. Catal. A: Chem., 2006, 249, 191–200.
25. M. V. Bollinger, K. W. Jacobsen and J. K. Nørskov, Atomic and electronic structure of MoS₂ nanoparticles, Phys. Rev. B, 2003, 67, 085410.
26. Y. D. Hou, B. L. Abrams, P. C. K. Vesborg, M. E. Bjorketun, K. Herbst, L. Bech, A. M. Setti, C. D. Damsgaard, T. Pedersen, O. Hansen, J. Rossmeisl, S. Dahl, J. K. Norskov and I. Chorkendorff, Bioinspired molecular co-catalysts bonded to a silicon photocathode for solar hydrogen evolution, Nat. Mater., 2011, 10, 434–438.
27. B. D. Yuhas, A. L. Smeigh, A. P. Douvalis, M. R. Wasielewski and M. G. Kanatzidis, Photocatalytic Hydrogen Evolution from FeMoS-Based Biomimetic Chalcogels, J. Am. Chem. Soc., 2012, 134, 10353–10356.
28. A. Roldan, Frontiers in first principles modelling of electrochemical simulations, Curr. Opin. Electrochem., 2018, 10, 1–6.
29. J. Rossmeisl, E. Skúlason, M. E. Björketun, V. Tripkovic and J. K. Nørskov, Modeling the electrified solid–liquid interface, Chem. Phys. Lett., 2008, 466, 68–71.
30. S. Trasatti, Structure of the Metal Electrolyte Solution Interface – New Data for Theory, Electrochim. Acta, 1991, 36, 1657–1658.
31. S. Haider, A. Roldan and N. H. de Leeuw, Catalytic Dissociation of Water on the (001), (011), and (111) Surfaces of Violarite, FeNi₂S₄: A DFT-D2 Study, J. Phys. Chem. C, 2014, 118, 1958–1967.
32. A. Roldan and N. H. de Leeuw, A kinetic model of water adsorption, clustering and dissociation on the Fe₃S₄{001} surface, Phys. Chem. Chem. Phys., 2017, 19, 12045–12055.
33. A. Roldan and N. H. de Leeuw, Catalytic water dissociation by greigite Fe₃S₄ surfaces: density functional theory study, Proc. R. Soc. A, 2016, 472, 20160080.
34. J. Tafel, Polarization in cathodic hydrogen evolution, Z. Phys. Chem., 1905, 50, 641–712.
35. G. Kresse and J. Furthmüller, Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15.
36. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
37. D. Santos-Carballal, A. Roldan, R. Grau-Crespo and N. H. de Leeuw, First-principles study of the inversion thermodynamics and electronic structure of FeM2X4 (thio)spinels (M = Cr, Mn, Co, Ni; X = O,S), Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 195106.
38. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
39. J. P. Perdew, J. Chevary, S. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Erratum: atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 48, 4978.
40. S. H. Vosko, L. Wilk and M. Nusair, Accurate spin-dependent electron liquid correlation energies for local spin-density calculations – A critical analysis, Can. J. Phys., 1980, 58, 1200–1211.
41. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem., 2006, 27, 1787–1799.
42. A. Roldan and N. H. de Leeuw, Selective hydrogenation of CO on Fe₃S₄{111}: a computational study, Faraday Discuss., 2017, 197, 325–336.
43. A. Roldan and N. H. de Leeuw, Methanol formation from CO₂ catalyzed by Fe₃S₄{111}: formate versus hydrocarboxyl pathways, Faraday Discuss., 2016, 188, 161–180.
44. S. Posada-Pérez, D. Santos-Carballal, U. Terranova, A. Roldan, F. Illas and N. H. de Leeuw, CO₂ interaction with violarite (FeNi₂S₄) surfaces: a dispersion-corrected DFT study, Phys. Chem. Chem. Phys., 2018, 20, 20439–20446.
45. N. D. Mermin, Thermal Properties of the Inhomogeneous Electron Gas, Phys. Rev., 1965, 137, A1441–A1443.
46. H. J. Monkhorst and J. D. Pack, Special Points for Brillouin-Zone Integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
47. A. J. Devey, R. Grau-Crespo and N. H. de Leeuw, Electronic and magnetic structure of Fe₃S₄: GGA+U investigation, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 195126.
48. V. I. Anisimov, M. A. Korotin, J. Zaanen and O. K. Andersen, Spin Bags, Polarons, and Impurity Potentials in La2-Xsrxcuo4 from 1st Principles, Phys. Rev. Lett., 1992, 68, 345–348.
49. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
50. I. D. R. Moreira, F. Illas and R. L. Martin, Effect of Fock exchange on the electronic structure and magnetic coupling in NiO, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 155102.
51. I. Ciofini, F. Illas and C. Adamo, Performance of the tau-dependent functionals in predicting the magnetic coupling of ionic antiferromagnetic insulators, J. Chem. Phys., 2004, 120, 3811–3816.
52. F. Illas and R. L. Martin, Magnetic coupling in ionic solids studied by density functional theory, J. Chem. Phys., 1998, 108, 2519–2527.
53. D. Munoz, N. M. Harrison and F. Illas, Electronic and magnetic structure of LaMnO₃ from hybrid periodic density-functional theory, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 85115–85124.
54. J. Wang, S. H. Cao, W. Wu and G. M. Zhao, The Curie temperature and magnetic exchange energy in half-metallic greigite Fe₃S₄, Phys. Scr., 2011, 83, 45702–45705.
55. G. W. Watson, E. T. Kelsey, N. H. de Leeuw, D. J. Harris and S. C. Parker, Atomistic simulation of dislocations, surfaces and interfaces in MgO, J. Chem. Soc., Faraday Trans., 1996, 92, 433.
56. P. W. Tasker, The stability of ionic crystal surfaces, J. Phys. C: Solid State Phys., 1979, 12, 4977–4984.
57. G. Henkelman and H. Jonsson, A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives, J. Chem. Phys., 1999, 111, 7010–7022.
58. A. Heyden, A. T. Bell and F. J. Keil, Efficient methods for finding transition states in chemical reactions: comparison of improved dimer method and partitioned rational function optimization method, J. Chem. Phys., 2005, 123, 224101.
59. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, Improved grid-based algorithm for Bader charge allocation, J. Comput. Chem., 2007, 28, 899–908.
60. R. F. W. Bader, Atoms in molecules, Acc. Chem. Res., 1985, 18, 9–15.
61. G. Henkelman, A. Arnaldsson and H. Jonsson, A fast and robust algorithm for Bader decomposition of charge density, Comput. Mater. Sci., 2006, 36, 354–360.
62. L. Chang, A. P. Roberts, Y. Tang, B. D. Rainford, A. R. Muxworthy and Q. Chen, Fundamental magnetic parameters from pure synthetic greigite (Fe₃S₄), J. Geophys. Res., 2008, 113, 1–16.
63. D. Santos-Carballal, A. Roldan, R. Grau-Crespo and N. H. de Leeuw, A DFT study of the structures, stabilities and redox behaviour of the major surfaces of magnetite Fe₃O₄, Phys. Chem. Chem. Phys., 2014, 16, 21082–21097.
64. A. Roldan and N. H. de Leeuw, Methanol formation from CO₂ catalyzed by Fe₃S₄{111}: formate versus hydrocarboxyl pathways, Faraday Discuss., 2016, 188, 161–180.
65. H. Eyring, The Activated Complex and the Absolute Rate of Chemical Reactions, Chem. Rev., 1935, 17, 65–77.
66. A. Einstein and R. Fürth, Investigations on the Theory of the Brownian Movement, Dover Publications, 1956.
67. L. Kristinsdóttir and E. Skúlason, A systematic DFT study of hydrogen diffusion on transition metal surfaces, Surf. Sci., 2012, 606, 1400–1404.
68. X. Shen, J. Chen, Y. M. Sun and T. Liang, Hydrogen diffusion on Fe surface and into subsurface from first principles, Surf. Sci., 2016, 654, 48–55.
69. A. B. Belonoshko, A. Rosengren, Q. Dong, G. Hultquist and C. Leygraf, First-principles study of hydrogen diffusion in α-Al₂O₃ and liquid alumina, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 024302.
70. C. Chen, H. Yu and S. Zheng, First-principles study of hydrogen diffusion mechanism in Cr₂O₃, Sci. China: Technol. Sci., 2011, 54, 88–94.
71. G. Jerkiewicz, Hydrogen sorption at/in electrodes, Prog. Surf. Sci., 1998, 57, 137–186.
72. W. Schmickler and R. Guidelli, Ionic adsorption and the surface dipole potential, J. Electroanal. Chem. Interfacial Electrochem., 1987, 235, 387–392.
73. IUPAC Compendium of Chemical Terminology, ed. M. Nič, J. Jirát, B. Košata, A. Jenkins and A. McNaught, IUPAC, Research Triangle Park, NC, 2009.
74. G. S. Karlberg, T. F. Jaramillo, E. Skúlason, J. Rossmeisl, T. Bligaard and J. K. Nørskov, Cyclic Voltammograms for H on Pt(111) and Pt(100) from First Principles, Phys. Rev. Lett., 2007, 99, 126101.
75. G. S. Karlberg, T. F. Jaramillo, E. Skúlason, J. Rossmeisl, T. Bligaard and J. K. Nørskov, Cyclic Voltammograms for H on Pt(111) and Pt(100) from First Principles, Phys. Rev. Lett., 2007, 99, 126101.
76. E. Skúlason, G. S. Karlberg, J. Rossmeisl, T. Bligaard, J. Greeley, H. Jónsson and J. K. Nørskov, Density functional theory calculations for the hydrogen evolution reaction in an electrochemical double layer on the Pt(111) electrode, Phys. Chem. Chem. Phys., 2007, 9, 3241–3250.
77. E. Skúlason, V. Tripkovic, M. E. Björketun, S. Gudmundsdóttir, G. Karlberg, J. Rossmeisl, T. Bligaard, H. Jónsson and J. K. Nørskov, Modeling the Electrochemical Hydrogen Oxidation and Evolution Reactions on the Basis of Density Functional Theory Calculations, J. Phys. Chem. C, 2010, 114, 18182–18197.
78. I. Chorkendorff and J. W. Niemantsverdriet, Concepts of Modern Catalysis and Kinetics, Wiley-VCH, Weinheim, FRG, 2003.
79. M. W. J. Chase, Journal of Physical and Chemical Reference Data. Monograph No. 9. NIST-JANAF Thermochemical Tables, 1998, p. 1951.

---

Footnote

† Electronic supplementary information (ESI) available: Details of the model used in this paper, the schematic representation of the pathways, the relation between S-band and E_B, H-coverage vs. local charge of sulfur in the pristine surface, as well as the effect on the derived potential at a neutral charge. See DOI: 10.1039/c8cp06371k

---