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

The mineral greigite, Fe 3 S 4 , 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, diﬀusion, 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 B 0.5 eV above the energy of the reactants on Fe 3 S 4 {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 2 evolution reaction becomes energetically barrierless.


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.][8] Although most biomimetic catalysts lack the protection and long-range interactions of the protein backbone, they show interesting catalytic activity. 9Fe 3 S 4 also has promising applications in spintronic devices, 10 electrodes for lithium-ion batteries, 11 drug delivery and hydrogen storage, [12][13][14] and it has been used as a catalyst in electro-reductive experiments. 15Understanding 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 3 S 4 surface indicate the formation of oxide shell structures (Fe 3 O 4 ), although the inner material remains a sulfide (Fe 3 S 4 ). 16,17Nevertheless, the behaviour of the greigite surface under anodic potentials remains unreported.][23][24][25][26] More sophisticated catalysts such as chalcogels enhance the evolution of H 2 but at very low turnover, even when doping the inorganic structure with other metals such as Ni, Co or Mo. 27ccurate 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. 28We 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. 29hus, the electronic structure calculations are related to the surface potential by the variations of the dipole. 30ased on previous research and the fact that H binds on top of low-coordinated S atoms on TMS surfaces, [31][32][33] we have used density functional theory with a correction for the long-range dispersion forces to focus on the Fe 3 S 4 {111} surface, which has been shown to be exposed in Fe 3 S 4 crystallites.In particular, we have analysed the H ad-atom migration, the surface reduction by either H-incorporation or by H 2 S formation, and the catalytic activity towards the hydrogen evolution reaction (HER) via the Tafel reaction. 34

Computational methods
Periodic plane-wave DFT-D calculations were carried out to study the hydrogen interaction and reactivity on the greigite Fe 3 S 4 {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,37Hence, the ion-electron interactions were represented by the projector-augmented wave (PAW) method 38 within the generalized gradient approximation (GGA) with the Perdew-Wang 91 functional 39 and employing the spin interpolation formula of Vosko et al. 40 All calculations included the D2 long-range dispersion correction approach by Grimme (s 6 = 0.75), 41 as this correction has been shown to reproduce well the interaction of the sulfides with molecules. 32,33,42,43The PBE functional and D3 long-range dispersion showed no difference with the setup used here. 44The Kohn-Sham valence states were expanded in a plane-wave basis set with a cut off at 600 eV for the kinetic energy. 45This high cut off energy ensured that no Pulay stresses occurred within the cell during relaxations.Calculations were carried out with Monkhorst-Pack 46 grids of 5 Â 5 Â 1 K-points, which ensured electronic and ionic convergence.The geometries of all stationary points were found by the conjugategradient algorithm and were considered converged when the force on each ion dropped below 0.03 eV Å À1 , whereas the energy threshold defining self-consistency of the electron density was set to 10 À5 eV.To improve the convergence of the Brillouin-zone integrations, the partial occupancies were determined using the tetrahedron method with Blo ¨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 ) 2 S 4 spinel structure with both types of Fe, i.e. octahedral (B) and tetrahedral (A), in a ferrimagnetic orientation, as reported previously. 37,47For 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. 502][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,541 Slab model After a full relaxation of the bulk, we prepared the Fe 3 S 4 surfaces by cutting the bulk structure and creating slab models using the METADISE code.55 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 3 S 4 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 Å 2 .
The slab is thick enough to relax the two uppermost layers (four Fe 3 S 4 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 Å 3 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,58We 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 1 2 H 2 molecules on the Fe 3 S 4 surface from eqn ( 1), where E nÁH:Fe 3 S 4 is the total energy of n hydrogens interacting with the slab, E Fe 3 S 4 is the energy of the naked slab and E H 2 is the energy of the isolated H 2 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. 61Within 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.

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. 63Upon 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. 32Here, we have focused on the {111} surface as it has shown interesting redox properties towards the reduction of CO 2 . 15,42,64This 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+.
This journal is © the Owner Societies 2019

Hydrogen adsorption and diffusion on Fe 3 S 4 {111}
We have approached a single H atom on non-equivalent adsorption sites on the Fe 3 S 4 {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 B 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.
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 3 S 4 {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 subsurface Fe B .
We have derived the H diffusion coefficient (D) on the Fe 3 S 4 {111} surface by evaluating the changes in the vibrational contribution to the energies within the limits of validity of the transition state theory. 65We used eqn (2) considering that a single displacement occurred when an H atom moved from one site to the next.
where DE ZPE is the activation energy corrected by the vibrational contribution at the ground state and the saddle point, and nd 2 /2a 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 a indicates the dimensions of the jump, which, according to the symmetry of the surface, is one. 66The 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 3 S 4 {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 2 .
The resulting diffusion coefficients indicate a hindered movement of the H ad-atom in comparison to values on transition metals, 67 e.g. the diffusion coefficients range from 3.35 Â 10 À4 to 2.2 Â 10 À3 cm 2 s À1 on body cubic centered Fe at 233 K. 68 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.

Hydrogen coverage on Fe 3 S 4 {111}
Actual surfaces are dynamic systems whose properties and processes depend on the concentration and types of adsorbents.The H coverage (y H ) is, therefore, a factor to consider when developing a general picture of the surface character.We defined  1. y 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 y 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 2 molecules in the gas phase.Fig. 2 represents the E B per atom as a function of y H starting from 1 up to 8 H atoms per cell, i.e. y 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 y H 4 1 ML on Pt{111}, which has been interpreted as a strong and repulsive interaction between H ad-atoms. 70 While neighbouring H ad-atoms are interacting by less than 0.1 eV on Fe 3 S 4 {111}, we related this step of 0.44 eV with the surface electronic structure.The increment of y 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 y H , while decreasing the overall stability of the system.Hence, at high y H , the S charge density localises along the S-H bonds generating a dipole perpendicular to the surface which drives some reconstruction. 56s 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 (o0.1 eV), which is in good agreement with experimental values in another catalyst. 71In 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

Hydrogen evolution on Fe 3 S 4 {111}
As the hydrogen coverage increases on the surface, so does the probability of their recombination into molecular hydrogen.We have studied the H 2 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 2 is adsorbed dissociatively, and its desorption may follow a similar trajectory. 73In the latter, the H 2 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 2 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.
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}. 74We 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 (y 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 /dy H , has been found to compare well with Monte Carlo simulations on hexagonal surfaces, such as the S site distributions, [75][76][77] whereas the vibrational entropy is commonly described by a Boltzmann distribution. 78For the gas phase product (H 2 ), 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 2 free energy (G) and the thermochemistry data from ref. 79, which is As general trends, the thermodynamic saturation of the Fe 3 S 4 {111} takes place at the coverage of B0.35 ML as the adsorption of H beyond this y H is endothermic, see Fig. 3A.Next, we studied process 5 in Table 2 to consider the evolution of H 2 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 2 is exothermic before reaching y H = 0.35 ML-the H 2 evolution process is lower in energy than the H adsorption-although it is hindered by the large activation energy (DG TS ).According to the Tafel reaction mechanism for the evolution of H 2 , 74 the activation and stability (DG B ) are almost unaffected by the presence of water, electric potentials, or electric fields, which are, however, important to simulate an electrochemical environment. 28Hence, y 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 2 evolution.It is at y H = 1 ML that the barrier is smaller than 0.2 eV (the difference between DG TS and DG B ) leading to an exothermic associative desorption of H 2 by 0.88 eV (the difference between DG R and DG B ).
As described in the hydrogen coverage section, an increment of y H polarises the surface, which affects the surface work function (f)-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 2 (final states) to an electrode potential at zero-charge (U = f À f 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.9V. 78 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 2 formation and the potential affecting the current density.Fig. 3B indicates how the H ad-atoms on the Fe 3 S 4 {111} surface accumulate under cathodic potentials-solid line in Fig. 3.Although the H 2 evolution becomes thermodynamically favourable below 0.35 V (vs.NHE)-the H 2 evolution is lower in energy than adsorption-, it is kinetically hindered-big difference between DG TS and DG B -until U = À0.4V (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. 15

Fe 3 S 4 {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. y H = 0.25 ML and y 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 2 S. Another way to reduce Fe 3 S 4 {111} is by removing sulfur from the surface, for example as H 2 S. We have studied hydrogen sulfide (H 2 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 2 S (E R = 0.09 eV), see Fig. 4. Once formed, the molecule hovers at B1.6 Å above its initial position.However, the energy barrier to place both H on the same S-site is 2.35 eV under low y H is 2.35 eV, which prevents the process from being kinetically accessible.At a y 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 2 S formation is thermodynamically possible, it is kinetically inaccessible under conditions of both low and high hydrogen coverage.

Conclusions
We have performed a systematic study of common processes involving hydrogen on greigite.We have used DFT-based Fig. 3 Free energies per H atom for the H binding (black circles), desorption (blue squares) and barrier (DG TS ) as a function of (A) hydrogen coverage (y H ) and (B) electrode overpotentials (U) vs. the normal hydrogen electrode at 300 K and an H 2 pressure of 1 atm.Note in the inset that asterisks denote adsorbed and free sites.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 2 evolution by ad-atoms recombination, and reduction of the surface by either H-incorporation or H 2 S formation on the Fe 3 S 4 {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 3 S 4 {111} surface, the strongest mono-atomic binding energy is À0.4 eV, which becomes weaker with coverage, following a logarithmic trend up to a y H of 1 ML; above this coverage, there is a trend discontinuity, not only due to the H-H interaction (o0.1 eV) but because of changes in the surface polarisation.The reduction of the Fe 3 S 4 {111} surface by incorporation of H into the thio-spinel structure is an endothermic process by more than 0.8 eV.H 2 S formation is in thermodynamic equilibrium with two contiguous HS, but the process is kinetically prevented between 0.25 and 1 ML (E A B 2 eV).The H 2 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 y H with the electrode potential via the work function and found that at bias of 0.4 V vs. NHE the H 2 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.

Fig. 1
Fig. 1 Top-view representation of non-equivalent H adsorption sites on the Fe 3 S 4 {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 inTable1.

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

Table 1
Barrier (E A ) and reaction energies (E R ) for H movement over Fe 3 S 4 {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

Table 2
Activation (E A) and reaction (E R ) energies for the molecular hydrogen formation and desorption processes on Fe 3 S 4 {111} with respect to the naked slab and the isolated H 2 molecule.MD and AD indicate respectively molecular and associative desorption.S 0 indicates a different S adsorption site 09% in a range of temperatures from 250 to 700 K, see Fig.S6 (ESI †).S conf = Àk B (y H Áln(y H ) + (1 À y H )Áln(1 À y H ))(3)