A kinetic model of water adsorption, clustering and dissociation on the Fe 3 S 4 {001} surface †

The interaction of water with catalyst surfaces is a common process which requires investigation. Here, we have employed density functional theory calculations to investigate the adsorption of up to ten water molecules on the {001} surface of greigite (Fe 3 S 4 ), which owing to its redox properties, is of increasing interest as a catalyst, e.g. in electro-catalysis. We have systematically analyzed and characterized the modes of water adsorption on the surface, where we have considered both molecular and dissociative adsorption processes. The calculations show that molecular adsorption is the predominant state on these surfaces, from both a thermodynamic and kinetic point of view. We have explored the molecular dispersion on the surface under diﬀerent coverages and found that the orientation of the molecule, and therefore the surface dipole, depends on the number of adsorbed molecules. The interactions between the water molecules become stronger with an increasing number of water molecules, following an exponential decay which tends to the interaction energy found in bulk water. We have also shown the evolution of the infra-red signals as a function of water coverage relating to the H-bond networks formed on the surface. Next we have included these results in a classical micro-kinetic model, which introduced the eﬀects of temperature in the simulations, thus helping us to derive the water cluster size on the greigite surface as a function of the initial conditions of pressure, temperature and external potential. The kinetic model concluded that water molecules agglomerate in clusters instead of wetting the surface, which agrees with the low hydrophilicity of Fe 3 S 4 . Clusters consisting of four water molecules was shown to be the most stable cluster under a wide range of temperatures and external potential.


Introduction
Understanding the properties of water at the interface with solids is of crucial importance in many areas of fundamental research and applications.Apart from the direct role of water in many (photo-)catalytic surface reactions, water-surface interactions play an intrinsic role in the description of electrochemical interfaces and in understanding wetting and corrosion, as well as fuel cell reactions which are of increasing contemporary interest due to environmental concerns. 1The properties of water manifest themselves in a variety of important applications, such as catalysis, electrochemistry, materials science, electronic devices, photocatalysis and photo-conversion, corrosion, geochemistry, adhesion, sensors, tribology, astrophysics and astro-chemistry, and even membrane science. 2major objective in studying water-surface interactions is to determine whether water is adsorbed molecularly or dissociatively on the surface.The chemical properties of the water dissociation products (OH À ,H + and O 2À ) are very different from that of the water molecule and may lead, for example, to the surface and bulk oxidation of many materials. 2,35][6] Molecular H 2 O can be used to probe sitespecific structure-reactivity relationships, especially on oxides and semiconductors, but to distinguish H 2 Of r o mO He x p e r imentally is complicated owing to the similarities in many of the properties of these two species.][17][18] Generally, the dynamic processes of H 2 O interaction with surfaces remain poorly understood, even for the extensively studied Fe 3 O 4 , due to the difficulties in the preparation of well-characterized single crystal surfaces. 19On Fe 3 O 4 , the oxide isomorph of greigite, 20 different adsorption modes have been identified upon its interaction with water: dissociative chemisorption and, at higher coverage, physisorbed H 2 O in a condensed ice conformation, 21,22 in agreement with molecular dynamics models. 17,23Studies using methods based on the density functional theory (DFT) have reported exothermic molecular and dissociative adsorption of one H 2 Oo naF e -t e r m i n a t e dF e 3 O 4 {111} surface 24,25 and a hydrogenbonded second water molecule which through its oxygen forms a hydronium ion-like structure. 26Hence, water molecules provide not only a reaction environment but also hydrogen, hydroxy and/or oxygen subspecies which can react further at the catalyst interface.Understanding processes like the synergetic adsorption and its implications on water dissociation is crucial for the development of high-performance catalysts. 27,28ater may bind to a surface through a variety of means, e.g.via electrostatics, charge transfer or hydrogen-bonding, and it has been detected on surfaces in the form of monomers, dimers, larger multimers, 2D bilayers, 3D clusters and a variety of ice structures. 1,2Clustered water is characterized by the presence of hydrogen-bonds between two or more neighboring water molecules.Owing to this ability to form hydrogen-bonding networks, water displays a variety of structures associated with clustering that span from dimers to liquid water and the diverse structures of the solid.Theoretical studies vary on the precise details of the adsorption geometry of small water clusters, although the formation of closed ring structures is generally accepted and has been verified experimentally by molecular beam electric deflection and microwave studies. 29,30At sub-monolayer coverage, water is known to diffuse on metal surfaces based on its clustering behavior. 2However, water clustering behavior is less clear on surfaces where water is strongly adsorbed, as for example on an oxide. 31,32In these cases, the minimum barrier of water diffusion from cation site to cation site, involves aligning its molecular dipole to the electronic corrugation between the two sites, but water interaction with both the substrate and other molecules dictates its orientation.Computer models are a valuable tool to identify the species on the surface, the energetics and the dipole involved in the water diffusion and clustering.In the following sections, we discuss the adsorption and dissociation of water on the Fe 3 S 4 {001} surface as a function of coverage.This surface has been identified experimentally in synthetized greigite that was used as an electro-catalyst 33 and understanding of the watersurface interface is crucial to understand the surface behavior under electro-chemical conditions.To this end, we have investigated different arrangements of the molecules, i.e. dispersed and in clusters, and implemented the results in a micro-kinetic model to close the gap between electronic structure calculations and experiment under specific conditions of temperature and water pressure.

Computational details
The greigite mineral (Fe 3 S 4 ) has an inverse spinel-structure, which is analogous to magnetite, 20 with a face-centered cubic structure and space group Fd% 3m, see Fig. S1 (ESI †).The cubic unit cell contains eight units of Fe 3 S 4 where the 32 anions are in a cubic close-packed arrangement, while eight of the tetrahedral sites and sixteen of the octahedral ones are occupied by the cations.The different cation arrangements of the thio-spinel formula can be rewritten as Fe A (Fe B ) 2 S 4 , where A and B denote tetrahedral and octahedral sites respectively. 20,34The high-spin alignment is antiparallel, leading to a half metallic character of the material, owing to the presence of Fe(II) in the octahedral sites.Hence, greigite has metallic properties along the one spin channel, whilst it presents a band gap along the other spin channel.

DFT calculations
We have carried out a systematic DFT-D2 study of the interaction of water (up to 10 molecules) with the Fe 3 S 4 {001} surface.
All calculations were performed using the Vienna Ab-initio Simulation Package (VASP). 35,36The Kohn-Sham valence states were expanded in a plane-wave basis set with a cut off at 600 eV for the kinetic energy, 37 where the ion-electron interactions were represented by the projector-augmented wave (PAW) method. 38his high value for the cut off energy ensured that no Pulay stresses occurred within the cell during relaxations.The electron exchange-correlation was calculated within the generalized gradient approximation (GGA) with the Perdew-Wang functional (PW91), 39 employing the spin interpolation formula of Vosko et al. 40 All the calculations included the long-range dispersion correction approach by Grimme, 41 which is an improvement on pure DFT when considering large polarizable atoms, 27,[42][43][44][45][46] although it may also overestimate the binding energies. 47We have used the global scaling factor parameter optimized for PBE (s 6 = 0.75).The initial magnetic moments were described with high-spin distributions in both types of Fe, i.e. octahedral (B) and tetrahedral (A) Fe in the Fe A (Fe B ) 2 S 4 spinel structure, by a ferrimagnetic orientation. 20,48e have employed Monkhorst-Pack grids of 4 Â 4 Â 1 K-points, which ensured electronic and ionic convergence. 49 eh a v eu s e dt h eH u b b a r d -l i k ea p p r o x i m a t i o n( U)f o ram o r e accurate treatment of the electron correlation in the localized d-Fe orbital of this transition metal. 50,51It improves the description of localized states in this type of systems where standard LDA and GGA functionals fail. 523][54][55] Following the approach used by Devey et al., 56 we have employed U eff =1e V which has been tested to give reliable data for catalytic processes. 33,57he geometries of all stationary points were found with the conjugate-gradient algorithm and considered converged when the force on each ion dropped below 0.03 eV Å À1 , whereas the energy threshold, which defines self-consistency of the electron density, was set to 10 À5 eV.In order 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.These smearing techniques can be considered as a form of finite-temperature DFT. 37

Slab model
The Fe 3 S 4 surfaces were prepared using the METADISE code to cut the bulk structure and create slab models. 58METADISE not only considers periodicity in the plane direction but also provides the different atomic layer stackings resulting in a null dipole moment perpendicular to the surface plane. 59The {001} surface was represented by a slab with an area of 81.0 Å 2 , containing 56 atoms (24 Fe and 32 S) per unit cell with a vacuum width of 12 Å between periodic slabs, which is big enough to avoid interactions between the slab and its images.The slabs are also thick enough to relax the four uppermost Fe 3 S 4 units until energy convergence, keeping the bottom layer frozen to model the bulk structure.To obtain the properties of an isolated H 2 O molecule, we placed it in the center of a 15 Â 16 Â 17 Å3 simulation cell to avoid lateral interactions and using the same criteria of convergence as for the iron sulfide slabs.

Characterization
We have described the atomic charges and derived magnetic moment by means of a Bader analysis, 60,61 where the electron and spin density associated with each atom is integrated over the Bader volume of the atom in question, as implemented in the Henkelman algorithm. 62Thus, due to the changes in the effective atomic radii with the oxidation state of the ion, the Bader volume is not calculated as a sphere of constant radius but is charge density dependent.Even so, the electron delocalization of the DFT method leads to an underestimation of the formal charges, although they can be used effectively in a direct comparison and to monitor changes in charges, for example as an effect of surface adsorption.
We have characterized each system by vibrational frequencies using the implemented numerical method in VASP, which takes finite displacements for every coordinate, and evaluates the second derivatives from the variation of the energy gradients in these displacements.In the harmonic approximation, these displacements have to be big enough to make a substantial variation of the energy in order to minimize the numerical errors in the calculation of the derivatives, but they have to be small enough to ensure vibrations in the harmonic regime.The eigenvalues of the diagonal Hessian matrix are the vibrational frequencies and the eigenvectors are the vibrational normal modes of the system.As strong adsorption at a surface may lead to the perturbation of the adsorbent's vibrations, we have considered the frequencies of water on the Fe 3 S 4 {001} surface decoupled from the surface phonons, in agreement with its weak interaction.The surface selection rule requires that the dynamic dipoles have a component perpendicular to the surface in order for resonant absorption to occur, the absorption intensity being related to the strength of the dynamic dipole.The adsorption symmetry of the molecule dictates whether or not a mode is infra-red (IR) active (dipole allowed) and so, for certain symmetries, even water modes with a net atomic motion parallel to the surface formally become dipole-allowed. 63The difficulty with the IR approach is that the vibrational frequencies are only a guide to the strength of the hydrogen-bond and hence also the O-O separation, and it is entirely possible that on a surface, a different local configuration leads to the same vibrational frequency. 64Frequency changes can be caused by intrinsic structural changes due to adsorption on the surface, mechanical re-normalization of the frequencies due to bonding with the surface, dipole-dipole coupling, or the coupling of vibrational and electronic states.
Dissociation of the molecule was also studied, leading to hydroxyl and hydrogen species co-adsorbed on the surface.The products are linked to the reactants by a single saddle point, the transition state (TS), which determines the kinetics of the process.Transition states were identified by means of the dimer method, 65,66 which searches for the TS by giving an initial atomic velocity towards the particular final state (product(s)).From an initial configuration, we generated the initial velocities by making two equal and opposite small finite-difference displacements in the coordinates of the reactant molecule.The method then finds a nearby saddle point by rotation and translation steps implemented by a conjugate gradient optimizer.The saddle point thus identified (TS) was confirmed by a vibrational frequency calculation, in which only one imaginary frequency is obtained, corresponding to the reaction coordinate.Afterwards, the dimer images are relaxed to the neighboring local minima.In a successful search, one of the images will minimize into the initial state and the other will become the final state.We could then define the energy barrier (E A ) for the dissociation process as the difference between the initial state (adsorbed molecule) and transition state (TS).The reaction energy (E R ) is the total energy difference between the final s t a t e( p r o d u c t s )a n dt h ei n i tial state (reactants).
We have calculated the binding energies (E B ) per molecule on the Fe 3 S 4 surfaces via eqn (1).
where Similarly, the addition interaction (E Int ) between the cluster and the surface is defined by eqn (3), where E # Fe 3 S 4 is the energy of the naked surface retaining the structure of [H 2 O] n :Fe 3 S 4 {001}.and {111} surfaces.We have focused our attention on the dominant {001} surface based on the relative surface energies. 28,33he Fe 3 S 4 {001} surface represents the (001), ( 010) and (100) terminations.It contains eight S atoms and seven Fe in the top layer of a unit cell, three Fe from a tetrahedral position in the bulk (Fe A )and four from an octahedral bulk position (Fe B ). 48 Before relaxation, it is terminated by 0.5 monolayers (ML) of 2-coordinated Fe A ions, occupying a bridge site (above two O ions) with a ffiffi ffi 2 p Â ffiffi ffi 2 p ÀÁ R45 symmetry according to Wood's notation. 67The flexibility of the S allows these Fe A to move below the uppermost layer as explained in previous work. 28The bulk structure, beneath the surface, consists of single rows in the [110] direction of 5-coordinated Fe B ions, alternating every two single rows of S ions with cubic packing. 68

Water adsorption
We first placed a H 2 O molecule on several non-equivalent positions on the surface and allowed both the surface and the molecule to relax without any restrictions.The most stable mode of adsorption, releasing an adsorption energy of 0.35 eV, takes place on one of the four Fe B present in the surface, whereas the S layer repels the lone electron pairs of oxygen, thereby blocking the attractive interaction with the Fe A .The molecule adsorbs with its oxygen at 2.335 Å from the Fe B with the H pointing towards neighboring S atoms, which stretches the H 2 O angle, see Fig. 1.The Fe B -site moves towards the molecule, rising 0.162 Å from the surface plane with respect to the clean surface, which rearranges its orbitals and modifies the Fe B spin density by 0.2 m B .There is polarization of the sulfur electron cloud (distance SÁÁÁH is 2.750 Å).
We next increased the number of H 2 O molecules up to a monolayer, see Fig. S2 (ESI †).We explored the position of a second water molecule across the surface and found the most favorable site to be close to the previously optimized H 2 Oonthe Fe 3 S 4 {001} surface, although the second H 2 O does not chemisorb on either nearby Fe B or Fe A .Upon optimization, the second molecule's oxygen interacts with the pre-adsorbed H 2 Oa sa n H-bonding acceptor, at 1.720 Å to the pre-adsorbed molecule and at B2.9 Å from the surface.This interaction enlarges the H-O* distance (we have used an asterisk to denote the chemisorbed molecules) by 0.03 Å; H 2 O* also moves closer to Fe B by 0.13 Å.
The presence of the second molecule strains the HOH angle of the pre-adsorbed molecule by 1.31.The electronic structure is mostly unaffected and only the bound Fe B site decreases its spin density by 0.3 m B compared to a single molecule adsorption state, while the sulfur polarization appears towards the second molecule H (d SÁÁÁH = 2.792 Å).
We have observed similar trends in the structures upon further increase of the number of molecules, where the formation of clusters is favorable with respect to the H 2 O dispersion on the surface, which takes place around the chemisorbed molecules.We have summarized the geometrical parameters in Table 1.While physisorbed water molecules are H-bonding donors at low coverage ([H 2 O] n n o 5), they become H-donors and acceptors at higher coverages, a sign of synergetic chemisorption.The general trend is to shrink the Fe B -H 2 O distance, while accommodating the molecules in clusters.Similar results were obtained on layered sulfide minerals, e.g.mackinawite. 69,70lusters with five or more members generate geometrical structures such as pentagons and hexagons, whereas square structures are less stable due to the geometry of the water molecule and the tension in the long-range interactions.][73] We also found that clusters with at least six molecules reverse the inward relaxation of the Fe A cations, with the cations adopting their positions prior to relaxation in a vacuum.This kind of surface reconstruction was not observed either on metals or oxides. 2 The reason for outward cation movement is because the water molecules chemisorb and thereby complete (or partially completes) the dangling bonds of the cations and reduce the dipole of the surface, thus removing the electrostatic driving force that relaxed the Fe A inwards.For coverages larger than 5 H 2 O molecules, the Fe A remains on the outside of the surface and binds up to two water molecules, see Fig. 2. When the monolayer is completed there are four water molecules bound to the surface, two on Fe B and two on the Fe A , whereas the rest of the water molecules are physisorbed and complete the H-bonding network.
The increasing number of adsorbed molecules on the surface is promoted by the partially charged bridging proton which, in a synergetic effect, compensates for the repulsive interaction between the oxygen atoms. 74The H-bond expands as new bonds are formed, converging to a basal plane ice structure similar to the one found on metals, where the water network adapts its  lateral lattice parameter to match the surface site spacing, thereby allowing the O atom of the lower water molecules to bind directly to the surface. 75Placing [H 2 O] 10 on a single Fe 3 S 4 {001} surface unit cell completely covers the surface leading to a monolayer, see Fig. 2. The full coverage leads to shorter Fe-H 2 O* bonds and decreases the distance between the surface and the hovering (physisorbed) water molecules in a comparable way to the structures found on metal surfaces. 1,2 3.2.1 Orientation of the molecules.The adsorption of water on metals leaves the molecules with one uncoordinated H atom to point away 76 from or downwards 77 onto the surface.The overall picture on the Fe 3 S 4 {001} surface, independently of the H 2 O cluster size, is that molecules orient the H down from the planar H-bond network towards the surface S-sites.This orientation does not require any hydrogen-bonds to be broken, so unless there is some additional constraint on the orientation and bonding of water (such as keeping the molecule flat to optimize bonding to the surface -as in low density water networks on the Pd(111) surface) 78 -molecule re-orientation probably requires minimal activation.
As the number of molecules increases so does the gap between the surface and the average O height, while the distance from the surface to the average H position decreases; a fact that confirms a preferred downward orientation, see side view in Fig. 2. The formation of a hydrogen-bond changes the polarization of the participating water molecules, increasing the water dipole moment and hence the ability to form further hydrogen-bonds. 79For instance, upon adsorption of the second molecule, the density flux showed substantial sulfur polarization towards the non-chemisorbed molecule. 28he water dipole, and therefore the molecule orientation, depends on its environment, i.e. adsorbate coverage, surface structure, co-adsorbates and applied potential. 2On the Fe 3 S 4 {001} surface, the orientation of the molecules and charge transfer play a role in the polarization of the system.We have calculated the variation in work-function (DF) upon water adsorption to provide information on the dipole orientation of the adsorbed cluster.DF is the difference between the Fermi energy of the different systems because the potential in the vacuum is identical for the systems compared.While the work-function varies from 4.40 to 5.49 eV, the average angles of the molecules with respect to the surface plane are in the range of À15.00 to +9.311,reachinguptoÀ83.081 and +33.041 for certain molecules.Positive and negative signs of the angles indicate the orientation of the hydrogens away and towards the surface respectively, see inset in Fig. 3.This figure shows the main angle of the plane of the molecules with respect to the surface as a function of DF.It shows that for n 4 1 the main orientation of the molecules is pointed slightly downwards to the surface.The angle and DF follow a linear trend until the surface strongly reconstructs (n 4 6), i.e. the Fe A adopts its bulk position above the S atomic layer and the surface electronic structure changes as a result.Although more points would be needed to define the limit accurately, full coverage indicates a strong orientation of the water-surface interface.In electrochemical experiments, the external potentials will modify the surface work-function, forcing the dipole orientation and thus water re-alignment. 802.2Adsorption energies.The interplay between waterwater and water-surface interactions is manifest in many of the properties of adsorbed water.In particular, the energetics and dynamics of water-water versus water-surface interactions rule the ability of water to cluster.From an energetic perspective, this relationship can be divided into three general categories: (i) where water-water interactions are comparable or stronger than those of water-surface interactions, (ii) where water-water interactions are comparable or weaker than those of water-surface interactions, and (iii) where water-water interactions are considerably weaker than water-surface interactions.The factors favoring cluster formation are rapid surface diffusion of water and weak water-surface interactions 81 compared to water-water interactions.The evaluation of these energies provides the thermodynamic basis for clustering, where strong water-surface interactions place water in an adsorption configuration in which hydrogen-bonding interactions may be structurally and thermodynamically unfavorable.Furthermore, surface diffusion is necessary in that the formation of clusters requires adsorbed water molecules to  This journal is © the Owner Societies 2017 approach each other; if diffusion is slow, bigger clusters are unlikely to form until the coverage is sufficiently high.
We have plotted in Fig. 4 the energy values corresponding to the cohesion energy (E Coh ), the adhesion energy (E Int ), and the global binding energy (E B ) per molecule as a function of the number of water molecules (summarized in Table 2).E Coh corresponds to the energy required to retain the molecules as a cluster, whilst E Int is the energy released (or required) for [H 2 O] n upon adsorption; both are calculated for water structures in the same configuration as the adsorbed one.This fixed geometry determines the number of H-bonds between molecules and the cohesion energy is therefore not directly comparable with gas phase water structures.However, the E Coh of free flat water structures show a similar trend.Compared to Mo ¨ller-Plesset and coupled-cluster computational levels, 82 the E Coh for the dimer is 0.01 eV weaker and the free structure of the trimer is 0.11 eV stronger than in our results; this value decreases with the number of water molecules.The sum of the E Coh and E Int interactions is described by E B .The slope of the E Coh trend is related to the formation of H-bonds, where a flat slope indicates thermodynamic equilibrium which is, however, still not reached with a monolayer.E Coh also decreases with the number of molecules while E Int remains at a similar energy, indicating easy exchange of molecule between layers.In general, the interaction between clusters and the surface is stronger than the cohesion between the molecules.This energy difference becomes smaller as the coverage, and number of H-bonds, increases, indicating a trend to form clusters instead of wetting of the surface.
The binding energy per water molecule (E B ) displays an asymptotic trend with the limit reaching the H 2 O-H 2 O interaction energy in an ice-Ih structure obtained with a similar computational setup. 83Full water coverage has practically the same binding energy as a molecule embedded in this hexagonal ice structure; a similar ice structure has been found on metal surfaces.Further increments in the number of monolayers would induce a partial rearrangement of the dipole due to the relocation of perpendicular H-bonds between the surface and the second monolayer.The three-dimensional H-bond mesh stabilizes the water structure by up to 0.3 eV depending on the geometry and the H-bond length. 845][26] Upon adsorption, the O-H bond of the molecule may become stretched towards the polarized sulfur on the surface until it dissociates, leading to both hydroxy and thiol groups.The dissociation energies, listed in Table 2 as E R , indicate an endothermic process by B0.8 eV which is independent of the number of water molecules; [H 2 O] 6 is stabilized by the relocation of long-range bonds without angle strains into a hexagonal shape.Although we have stretched the H-O* from one of the chemisorbed molecules towards the surface, the relaxation and optimization of the [H 2 O] n (n = 6 and 10) geometries led to hydrogen exchange within the H-network leaving the hydroxyl group adsorbed on Fe A , see Fig. 5.
The OH group subtracts on average B0.5 e À from nearby metal centers, whereas the S-H bond has covalent character.Similar charge transfers have been reported for H 2 Oo nt h ep y r i t e { 1 0 0 } surface, reflecting the greater stability of pyrite and its inability to be further oxidised. 85The dissociation kinetics are also unfavorable, with activation barriers of B1.1 eV above the adsorbed molecular state.5][26] This difference can be attributed to the lower basicity and ionic character of the sulfide compared to the oxide.The reverse process, where hydroxy and hydrogen associates to yield adsorbed H 2 O molecules, has an easily achievable transition state (B0.2 eV).

IR. Vibrational frequencies and intensities can pro-
vide qualitative information about the degree of hydrogenbonding and the geometry of water in different structures; for example, there have been several studies on oxide surfaces, which proposed monomeric water adsorption based on observations of non-hydrogen-bonded OH stretching. 1,2In this section we have characterized the simulated infra-red spectra of the different water clusters on the Fe 3 S 4 {001} surface.
The adsorption of a single molecule shows a bending vibrational mode at around 1590 cm À1 , while the symmetric and    1 and 2.
In addition to the frequency correlation, there is also a correlation between the strength of the hydrogen bond and the intensity of the IR absorption band, indicating that hydrogenbonding also increases the dynamic dipole moment, as we showed in Fig. 3. 87 As the number of water molecules increases, the stretching mode of the adsorbed molecule (n(HO*)) shifts towards lower frequency, indicating a weakening of the O-H bond.Asymmetric n a (HO*) are at higher frequencies than the symmetric stretching which, together with H-bonds at low-frequency n(OÁÁÁH) modes, 2 form the limits of bulk liquid water, 2800 and 3800 cm À1 ,dashed line in Fig. 6.The development of the H-network is detected by the slight shifts up in the scissoring vibrational mode (d(HOH)), at B1600 cm À1 , reflecting the increased stiffness of this motion.Below 1200 cm À1 , we have found the vibration modes, twisting, rocking and wagging modes, forming a wide band from B1000 cm À1 , in full agreement with the experimental IR of water.The similitudes to the bulk water IR reflect not only the quality of the model but the significant dynamics of the water clusters on the greigite surface. 71,72

Micro-kinetics
Although we have established the energetic preference for H 2 O adsorption and dissociation on the Fe 3 S 4 {001} surface, the preferred cluster dimension is still unknown.We have therefore employed a classical micro-kinetics model to reveal the relative fraction of clusters under various conditions of temperature, pressure of H 2 O and external potential, see ESI. † We derived the free energy (G) from our DFT calculations, where the internal energy is corrected by the zero point energy and its variation with the temperature by the specific heat; the entropic term is derived from the vibrational, translational and rotational partition functions of the molecules, while we kept the electronic partition function equal to the ground state. 48,88,89The energy values of the clusters containing 7, 8 and 9 molecules of water were interpolated following the linear relationships as a function of the number of water molecules which also depends on the temperature, see Fig. S3 (ESI †).We have compared the calculated free energy of H 2 O with the free energy derived from the Shomate equations 90 resulting in an average relative error of 0.15%, see Fig. S4 (ESI †).
The model considers processes of adsorption/desorption in thermodynamic equilibrium, where the cluster increases from gasphasemoleculesfollowingtheprocessesineqn( 4  The probability of the water molecules sticking on the surface was derived from the difference between the partition functions of a gas phase H 2 O with two-dimensional degrees of freedom, the third one corresponding to the reaction coordinate, and the adsorbed system. 88Within this definition, the sticking coefficient  This journal is © the Owner Societies 2017 (S 0 ) is of the order of 10 À5 , which depends on the temperature according to the expression S 0 = 20.605Â T À2.5 .It lies in the range of S 0 found on metal and semiconductor surfaces under low H 2 O coverage, 2 which also agrees with the geometric and energetic results previously discussed.
The reaction constants reveal that beyond a certain temperature it is unlikely to find water molecules scattered on the surface, because the desorption process is faster than the adsorption for most of the initial coverages, which agrees with the low hydrophilicity of greigite. 91Exceptions are the processes leading to a stable geometry configuration where the H-bond network is extended, i.e. initial y H 2 O = 0.3, 0.4, 0.8 and 0.9 ML, the last one leading to a very stable full coverage.Fig. 7 shows the natural logarithm of the adsorption/desorption reaction constants (K ads /K des ) as a function of the temperature.The behaviour of the constants agrees with the fact that at higher temperature the entropic term of the free energy becomes more important and the molecules tend to desorb from the surface as K ads /K des gets closer to one.
We have also implemented in our model a well-established approximation, extensively discussed in the literature, 80,[92][93][94][95][96] to simulate the effect of an external potential, e.g. an electrode potential at charge zero.The electrode bias modifies the free energies of the redox processes, i.e. during H 2 O dissociation there is charge transfer to the [H + OH] groups.This methodology does not include effects such as the electrical double layer parallel to the surface, which, however, is usually small due to the small dipole moment of the adsorbents on the surface and low external potentials. 97As expected, anodic potentials (oÀ0.6 V vs. NHE) enhance the dissociation of H 2 O at room temperature.Milder or positive external potentials cause the hydroxy group to recombine with hydrogen and to desorb as water, as shown in Fig. 8 and Fig. S5 (ESI †).This finding is in good agreement with the experimental potential value during the hydrogenation of CO 2 , while a larger anodic bias led to hydrogen evolution. 33,48he model allows the clusters to grow by adsorbing H 2 Ot o an already existing cluster.Thus, we have analysed the cluster distributions as a function of the temperature and found that the Fe 3 S 4 {001} surface has a hydrophobic character, in agreement with previous experimental reports. 91 We have also considered a separate kinetic model where the surface starts with pre-adsorbed water molecules.Here, the molecules move along the surface but the desorption is restricted, according to the process expressed by eqn (5).The results indicate the stability of the different clusters as a function of the initial conditions and increase in their internal energy, i.e. temperature.

Conclusions
We have studied the adsorption of water molecules up to a full monolayer on the {001} surface of the iron thio-spinel Fe 3 S 4 , which is the sulfide analogue of magnetite.All the calculations were performed by DFT+U, including the long-range dispersion forces as derived from the atomic polarizability.The sulfide allows movements perpendicularly to the surface, causing the Fe originally above the surface plane to penetrate beneath the top atomic layer, thereby decreasing the surface energy.However, the original bulk-terminated position of Fe A in the {001} surface is restored when it is surrounded by more than six waters, two of them binding Fe A and completing its tetrahedral coordination, thereby decreasing the surface dipole and the electrostatic driving force for relaxation.While chemisorbed molecules bind Fe B sites, physisorbed ones hover around the coordinated H 2 O. Hovering molecules orient themselves to form long-range interactions between themselves and with the sulfurs in the surface, producing geometric structures similar to the basal plane of hexagonal ice, whereas the water network adapts its lateral lattice parameter to match the site spacing.The average orientations of the molecules within the cluster size can be linked to the work-function through the overall dipole moment, particularly for large coverages.Thegrowthinthenumberofmoleculesincreasesthecohesion interaction between molecules by increasing the number of H-bonds, although the interaction energy with the surface remains practically invariable with increasing cluster size.This finding explains the stabilization energy by neighboring molecules whilst showing the low hydrophilicity of the sulfide.The energy balance between binding Fe B and increasing the size of existing water clusters is also shown by the micro-kinetic model, where the average size are four-and five-membered clusters (depending on the temperature).This model also showed the distribution of dissociated molecules as a function of the temperature and external potential: anodic potentials below 0.6 V vs. NHE dissociate H 2 O on a greigite electrode, which is in agreement with previous electro-catalysis work.The main interactions between adsorbates and adsorbate-surface is determined by long-range electrostatic interactions, especially in polar molecules like water, which mainly hover practically parallel to the surface, similar to those interacting with metallic surfaces.Overall, the Fe 3 S 4 {001} surface has a hydrophobic character with low desorption temperatures.In addition, wettability of the surface is very low, in agreement with experimental studies, and water clusters remain as they are and do not spill over the surface.These results show the different behavior of sulfides compared to oxides, which usually dissociate H 2 O spontaneously.
Further work is needed to ratify the accuracy of the electronic structure calculations and micro-kinetic simulations.Although in many cases they provide a reliable description of experimental features, several approximations have been considered.For instance, future higher level calculations, including long-range dispersion depending on the electron density instead of semiempirical approximations, and an explicit description of dynamical effects would help to verify the results presented here.
4 is the total energy of the cluster with n molecules on the Fe 3 S 4 {001} surface, E Fe 3 S 4 is the energy of the naked Fe 3 S 4 slab and E H 2 O is the energy of an isolated molecule in vacuum.The cohesive energy (E Coh ) of the molecules in the cluster was derived using eqn(2), where (E [H 2 O] n )isthewaterclusterenergywithoutthe slab in the same configuration as in the co-adsorbed situation.

Fig. 1
Fig. 1 Top and side view of Fe 3 S 4 {001} surface with [H 2 O] 1 molecule.The color scheme represents O in red and H in white balls and sticks while Fe and S are grey and yellow sticks respectively; black lines delimit the unit cell perimeter.All the distance values are given in Å.

Fig. 2
Fig. 2 Top and side view of [H 2 O] 10 molecules forming a monolayer on Fe 3 S 4 {001}.The color scheme represents O in red and H in white balls and sticks while Fe and S are grey and yellow sticks respectively; black lines delimit the unit cell perimeter.All the distance values are given in Å.

Fig. 3
Fig. 3 Average angle of the H with respect to the surface (see inset) as a function of the work-function variation (DF) respect to naked surface Fe 3 S 4 {001}.Number of molecules in the cluster [H 2 O] n are inset.Negative angle values indicate orientation towards the surface.Dashed line represents the interpolation for 3 Z n Z 6 following the equation angle = À5.4597À 3.3268ÁDF; R 2 = 0.99.

Fig. 4
Fig. 4 Cohesion (E Coh ), adhesion (E Int ) and binding (E B ) energies per water molecule on Fe 3 S 4 {001}.The cohesive energy is presented as a dashedred line, the adhesion energy is a dotted-blue line, and the binding energy is black-solid line.The dashed brown line indicates the DFT(PW91) H 2 O-H 2 O binding energy in ice-Ih the structure.83 asymmetric stretching peaks are practically negligible due to the modest dipole perpendicular to the surface, in agreement with the geometry in Fig.1.Upon co-adsorption, the formation of an H-bond is accompanied by shortening of the O-O distance and a slight elongation and weakening of the O-H bond.In turn, this reduces the vibrational frequency of the O-H bond (see Fig.6), which, shows a strong correlation to both the enthalpy of the hydrogen bond and the O-O distance, see Tables ),(R = 1-10).The model also includes the dissociation of a chemisorbed molecule into H + OH as a function of the coverage (R = 11-20).

Fig. 5
Fig. 5 Top view of reactants (left), transition state (middle) and products (right) of the H 2 O dissociation process on Fe 3 S 4 {001}.The color scheme represents O in red and H in white balls and sticks while Fe and S are grey and yellow sticks respectively; circles and the arrow indicate the dissociated hydrogen, the hydroxy group and the direction of proton transfer respectively.

Fig. 6
Fig. 6 Normalized infra-red spectra of [H 2 O] n (n = 1-10) clusters on Fe 3 S 4 {001}.We have shifted the y-axis position of each spectrum to enhance the visualization.The dashed line represents the IR of liquid water (Standard Reference Database 69: NIST Chemistry WebBook).
Fig. 9 shows the cluster size distribution under an applied bias of À0.6 V and open circuit (0.0 V).The process starts with the adsorption of one molecule followed by the co-adsorption of more H 2 Ot o form a five-membered cluster.The relative amount of [H 2 O] 1 depends on the initial concentration of H 2 O gas and the surface bias.Accordingly, the appearance of [H 2 O] 5 is related to the decline of a single molecule per unit cell.The adsorption and desorption processes lead to the formation of four-membered clusters [H 2 O] 4 at the expense of [H 2 O] 5 in a temperature range of 200-225 K.This four-membered cluster remains on the surface up to 275 K when the molecules desorb.The remaining molecules [H 2 O] 1 (o10 À15 ML) are in equilibrium with their dissociated form [H + OH].

Fig. 7
Fig. 7 Natural logarithm of the adsorption/desorption reaction constants as a function of the temperature.Different initial coverages are plotted (0 o y H 2 O o 0.9 ML).

Fig. 8
Fig.8Monolayers of adsorbed and dissociated water on Fe 3 S 4 {001} under different conditions of temperature and external potential two seconds after exposure of the naked surface to a pressure of 0.7 ML of gas phase water.

Fig. 9
Fig.9Normalised fraction (w) of water clusters on Fe 3 S 4 {001}, as a function of the temperature, under different applied potentials, À0.6 V and 0.0 V (left and right respectively), and for three initial pressure of H 2 O, 0.1, 0.5 and 0.9 ML from top to bottom.

Table 2
Binding (E B ), cohesion (E Coh ) and interaction (E Int ) energies per water molecule (n) for the molecular adsorption on Fe 3 S 4 {001}.The dissociation of one adsorbed H 2 O is characterized by the activation (E A ) and the reaction (E R ) energies