CO 2 interaction with violarite (FeNi 2 S 4 ) surfaces: a dispersion-corrected DFT study

The unbridled emissions of gases derived from the use of fossil fuels have led to an excessive concentration of carbon dioxide (CO 2 ) in the atmosphere with concomitant problems to the environment. It is therefore imperative that new cost-eﬀective catalysts are developed to mitigate the resulting harmful eﬀects through the activation and conversion of CO 2 molecules. In this paper, we have used calculations based on the density functional theory (DFT), including two semi-empirical approaches for the long-range dispersion interactions (-D2 and -D3), to explore the interaction of CO 2 with the surfaces of spinel-structured violarite (FeNi 2 S 4 ). This ternary sulfide contains iron ions in the highest possible oxidation state, while the nickel atoms are in the mixed 2+/3+ valence state. We found that CO 2 interaction with violarite is only moderate due to the repulsion between the oxygen lone pairs and the electronic clouds of the sulfur surface atoms. This suggests that the CO 2 activation is not dictated by the presence of nickel, as compared to the pure iron-isomorph greigite (Fe 3 S 4 ). These results diﬀer from findings in (Ni,Fe) ferredoxin enzymes, where the Ni/Fe ratio influences the redox potential, which suggests that the periodic crystal structure of violarite may hinder its redox capability. CO 2 molecule is activated upon interaction with this FeNi 2 S 4 thiospinel and the role of the Ni atom on the catalytic properties by comparison with other pure iron sulfide phases. We have employed calculations based on the density functional theory (DFT) and


Introduction
The unbridled anthropogenic emissions of carbon dioxide (CO 2 ) into the Earth's atmosphere from burning fossil fuels have caused a rise in the average global temperature, owing mostly to the greenhouse effect. 1 Current estimates suggest that CO 2 emissions will continue to increase until 2040, in line with the global energy demands, 2 with potentially irreversible consequences to life on Earth. For this reason, many efforts have been dedicated to capture and utilise CO 2 as a chemical feedstock to produce platform chemicals. 3 Commonly, the limiting rate during CO 2 conversion is the activation of the CQO bond, which typically implies the bending of the molecule, a process often induced by occupation of the lowest unoccupied molecular orbital (LUMO) of the neutral molecule, thus displaying a clear antibonding character. 4 This is, however, a challenging process because CO 2 is extremely stable and most often interacts only weakly with many surfaces. Activation may be achieved by adding promoters to an otherwise inactive surface 5 or by making use of more reactive catalytic materials. In this respect, a significant number of studies have been reported in the literature on the CO 2 adsorption, activation, and conversion on metals, 6,7 metal-oxides, 8,9 and metal organic frameworks. 10,11 Recently, it has been shown that molybdenum carbide surfaces strongly activate this very stable molecule. 12,13 In the current search for new catalysts, transition metal sulfides have gained attention as a promising alternative, due to their high abundance, low cost, and prominent catalytic features. For example, iron sulfides have been suggested as catalysts for the activation and conversion of CO 2 , 14,15 thereby also supporting the iron-sulfur hypothesis for the origin of life by Wächtershäuser and co-workers, which proposes that many of the prebiotic chemical reactions may have been catalysed by (Ni,Fe) sulfide minerals. 16,17 Previous studies have revealed that Ni is vital to modulate an adequate redox potential in (Ni,Fe) complexes, e.g. ferredoxin enzymes, and therefore to control their (bio-)catalytic activity. 18,19 However, the catalytic properties of the (Ni,Fe) sulfide environment within a periodic crystal structure, and thus whether Ni may enhance or hinder the redox capability in (Ni,Fe) sulfides compared to the pure iron sulfides, are not yet fully understood.
In this work, we have focused on determining whether the CO 2 molecule is activated upon interaction with this FeNi 2 S 4 thiospinel and the role of the Ni atom on the catalytic properties by comparison with other pure iron sulfide phases. We have employed calculations based on the density functional theory (DFT) and periodic slab models to study the CO 2 interaction with several different {001} and {111} low Miller index surface terminations of the stoichiometric ternary thiospinel violarite (FeNi 2 S 4 ). 20 In particular, we have explored the performance of a number of exchange correlation functionals and dispersion correction approximations, before comparing and rationalising our results with previous reports on greigite (Fe 3 S 4 ), and magnetite (Fe 3 O 4 ). 9,15

Computational details
The periodic DFT-based calculations described in the present work have been carried out using the Vienna ab initio simulation package (VASP). 21 The valence electron density was expanded in a plane-wave basis set and the effect of the core electrons (up to and including the 3p for Fe and Ni, 2p for S and 1s for C and O) on the valence region was described by the projector augmented wave (PAW) method by Blöchl,22 as implemented by Kresse and Joubert. 23 The kinetic energy cut-off for the plane-wave basis set was truncated at 600 eV leading to negligible Pulay stress. The threshold for the convergence of the electronic optimization was 10 À5 eV, while the relaxation of the atomic positions was allowed until the forces acting on all the atoms were smaller than 0.01 eV Å À1 . G-Centred Monkhorst-Pack grids of 5 Â 5 Â 5 and 5 Â 5 Â 1 k-points in the reciprocal space were used for the simulation of the bulk and surface slabs, respectively. 24 A single G-centred k-point was used for the simulation of the isolated CO 2 molecule in a broken symmetry cell of 9 Å Â 10 Å Â 11 Å to avoid spurious interactions.
For each of our simulations, we have optimised the geometries and evaluated the energies at two different levels of the generalized gradient approximation (GGA), i.e. using the Perdew-Wang 91 (PW91) 25,26 and the Perdew-Burke-Ernzerhof (PBE) 27,28 exchange-correlation functionals and taking into account spin polarization. The PW91 calculations were corrected with the spin interpolation formula of Vosko et al. 29 To remain consistent with previous reports, 9,14,15 long-range dispersion interactions were added to all the simulations via the D2 30 and the zero-damping D3 (zero) 31 semi-empirical methods. Note that different vdW schemes show the same trends albeit with some differences in the absolute values. 32 The global scaling parameter of s 6 = 0.75 developed for the D2 correction of the PBE functional was also used for the PW91 calculations, in line with previous studies. 9,14,33 The DFT+U formalism was used following the method suggested by Dudarev 34,35 to better describe the local character of the strongly correlated 3d electrons of the Fe and Ni atoms. Following earlier works, a U eff value of 1.0 eV was applied to the Fe 3d states when using the PW91 functional. 33 For the calculations with the PBE functional, we have chosen larger U eff values for 3.5 eV and 4.5 eV for the Fe and Ni 3d orbitals, which had been found to be required for the closely related PBEsol functional. 36 The atomic charges and spin moments were analysed using the Bader Atoms in Molecules formalism 37 according to the implementation by Henkelman et al. 38,39 which is compatible with VASP. For each surface, the most favourable adsorption mode has been analysed further by means of the charge density difference (CDD). The CDD plots have been obtained following eqn (1), where r is the CDD, r A-B is the electron density of the CO 2 on the surface, r A is that of the surface after the adsorption but without the adsorbate, and r B is that of CO 2 in the adsorption geometry without the substrate.

Results and discussion
3.1. Bulk model where U eff = 1.0 eV was applied to Fe, have revealed a completely inverse distribution of cations, 33 in agreement with powder diffraction measurements between 100 and 573 K, 40 Mössbauer data 41 and extended X-ray absorption 42 studies. However, a more recent computational study, where both Fe and Ni d orbitals were corrected with a Hubbard Hamiltonian by 3.5 and 4.5 eV, respectively, has shown that normal violarite is energetically more stable than the inverse form by 0.40 eV. 43 Santos-Carballal et al. 43 hypothesised that the synthesis of violarite produces the inverse spinel as a kinetic product, whilst its formation in the ores deep below the Earth's surface leads to the thermodynamic product with the normal distribution of cations. Calculations for the bulk phases have been carried out using the conventional cubic unit cell with an inverse spinel structure, see Fig. 1. Optimized lattice parameters of 9.402, 9.635, 9.713 and 9.709 Å have been predicted by the PW91-D2+U, PBE-D2+U, PBE-D3+U and PBE+U functionals respectively. These results compare well with the experimental value of 9.465 Å and show that, whereas changing from PW91 to PBE introduces a noticeable difference, the effect of the dispersion correction is minimal. The difference between PW91 and PBE is quite surprising and at variance with reported behaviour in other systems, e.g. the full set of transition metals. 44 Inverse violarite is a ferrimagnetic thiospinel with metallic character, where the magnetic moments of the atoms filling the tetrahedral positions are aligned antiparallel to the spin moments of the ions occupying the octahedral sites, 45,46 leading to a magnetisation of saturation of B4.00 m B per formula unit. 43 Disappointingly, the PW91+U functional is unable to describe the appropriate spinel magnetic configuration. Note that, regardless of the methodology used, we have set the initial magnetic moments parallel within each sublattice and antiparallel to the other sublattice. They were allowed to relax during our simulations. However, the calculations with the PBE+U functional predict the correct magnetic order, as shown in Table 1. In view of the close similarity of the PW91 and PBE functionals, we attribute the different descriptions of the two exchange-correlation approximations to the different U eff value used for the Fe and Ni 3d electrons.

Surface models
The different non-polar terminations of the two low Miller index surfaces studied here were created from the bulk optimised structure using the METADISE code (minimum energy techniques applied to dislocation, interface and surface energies). 47 This code analyses any dipole perpendicular to the surface following the approach pioneered by Tasker. 48 Hence, it considers the crystal structure as a stack of atomic planes parallel to the surface plane, where three possibilities arise: type 1 characterised by zero charge planes; type 2 where the planes are charged, but their periodic arrangement cancels the dipole moment; and type 3, where surface reconstructions (e.g. atom vacancies) are required to quench the dipole moment formed by the alternating charged planes. At this point it is also worth noting that the Tasker rule is not broken by surface relaxations as long as the atoms forming the stoichiometric unit do not leave their planes. We have used both possible charge arrangements for the cations, i.e. Ni 3+ (Fe 2+ Ni 3+ )S 4 2À and Ni 3+ (Fe 2.5+ Ni 2.5+ )S 4 2À to ensure the construction of all non-polar slabs. The former distribution of oxidation assumes different charges for the octahedral Ni and Fe ions, which is in agreement with experimental data for the inverse FeNi 2 S 4 , 42 whereas the latter considers that all the atoms located on the octahedral positions have the same charge (Ni +2.5 and Fe +2.5 ), in line with the itinerant electron magnetism properties of magnetite (Fe 3 O 4 ). 49 To minimise the surface stress, we relaxed different number of atomic layers until energy convergence, leading to a surface model where the top half of the slab was allowed to relax without restrictions, while the bottom half was kept frozen at the atomic bulk positions. The surface energy (g u ) of the unrelaxed slabs measures the excess energy of the pristine surface with respect to the bulk, before accounting for surface relaxation, and it is calculated as where E unrelaxed slab is the energy of the unrelaxed surface slab, E bulk is the energy per formula unit of bulk FeNi 2 S 4 , n is the number of stoichiometric units in the surface model with respect to the bulk and A is the surface area. Upon relaxation of half of the atomic layers in the surface model, the surface energy (g r ) can be obtained as in eqn (3) where E relaxed slab is the absolute energy of the relaxed slab. We express g r as a function of g u because the slab model we are using does not provide an isolated relaxed surface energy, since both the relaxed and unrelaxed sides of the slab are considered.
The stabilization by the surface energy is expressed as a percentage following eqn (4): The adsorption energy (E ads ) of CO 2 was calculated as in eqn (5) E ads = E CO 2 /Surf À (E Surf + E CO 2 ) where E CO 2 /Surf is the total energy of the CO 2 molecule interacting with the violarite surface, E Surf is the energy of the corresponding naked relaxed surface, and E CO 2 is the energy of the isolated CO 2 molecule in the gas phase.
3.2.1 Non-polar terminations of the FeNi 2 S 4 {001} surface. Table 2 collects both unrelaxed and relaxed surface energies for all the non-polar terminations of the FeNi 2 S 4 {001} surface, while their structures are shown in Fig. ESI-1 of the ESI. † Regardless of the density functional used, the most stable termination before relaxation corresponds to termination 8, from Table 1, which contains an 0.5 monolayer (ML) of Ni A atoms above the topmost bulk-like surface layer. Upon relaxation, we found that the same surface termination 8 has the lowest relaxed surface energy for the PW91-D2+U functional. Fig. 2 shows surface termination 8 optimised with the PW91-D2+U functional, where after relaxation the Ni A ad-atoms have moved towards the sub-surface layer and become octahedrally coordinated. We had expected that Fig. 1 (a) Lateral view and (b) polyhedral representation of the conventional unit cell of violarite. Balls in yellow, grey, soft blue, and dark blue represent the S, Fe, Ni A , and Ni B atoms, respectively. Ni A atoms are represented using bigger balls than Ni B ions. the lowest energy termination would be the same for each functional used. However, the optimisations with both PBE-D2+U and PBE-D3+U functionals predict that surface termination 2 is the most stable (see Table 2), a result which is again attributed to the different U eff value used in the PW91+U and PBE+U calculations. Following these findings, the CO 2 interaction has been studied on termination 2 with the PBE-(D2,D3)+U functionals, and on terminations 8 and 2 for PW91-D2+U. We have studied both surfaces with PW91-D2+U, because termination 8 is the lowest energy termination using this functional and including termination 2 allows us to compare the effect of the different functionals on the CO 2 interaction, where we disregard some surface terminations with similar 7 or even slightly lower 9 surface energies. Fig. 2 shows the top and side views of relaxed surface termination 2 for each functional used. In general terms, independently of the method used, the Ni A atoms have moved inwards during geometry optimisation and become octahedrally coordinated, despite the results from the PBE related functionals led only to a small inwards migration of the most exposed sulfur anions. The differences in extent of surface reconstruction obtained with the PW91-D2+U and PBE-(D2,D3)+U functionals are quantified by the % relaxation, see Table 2.
3.2.2 Non-polar terminations of the FeNi 2 S 4 {111} surface. Fig. ESI-2 (ESI †) displays the structure of the four {111} surface terminations before geometry optimisation. We found that the relaxed surface energies obtained using the PW91-D2+U and PBE-D3+U functionals exhibit the same trend, as is clearly seen from Table 3. The analysis of the surface energies reveals that termination 4 becomes the most stable after relaxation, although some structural differences arise due to the use of different exchange-correlation approximations. Fig. 3 illustrates the inward movement of Fe B and Ni A ions, originally sited above the S topmost layer, to (i) the subsurface layer (PW91-D2+U), or to (ii) a single plane in line with the S atoms (PBE-(D2,D3)+U).

Interaction between CO 2 and FeNi 2 S 4 surfaces
Here we consider the interaction of CO 2 with the different surface models described in the previous section. The calculations start by approaching a pre-activated CO 2 molecule to all non-equivalent sites on the FeNi 2 S 4 {001} and {111} surfaces. A pre-activated geometry was used as an initial set-up, using a bent molecule with apex angle of 1301 in three different orientations with respect to the surfaces (see Fig. ESI-3, ESI †). In all our simulations, we placed the CO 2 molecule at 1.5 Å from the surface to favour the attractive forces over the repulsive ones between the adsorbate and substrate. However, during geometry optimisation, both the surface and the molecule are free to move and change their geometries. As is also discussed in the following, the CO 2 molecule actually becomes linear during the optimisation process, independently of the site, surface termination, and functional used. . From an initial distance of 1.5 Å from the surface, the molecule moves away to about 3 Å and the pre-activated geometry relaxes to a linear CO 2 molecule. In one of the interaction modes, CO 2 hovers at a tilted angle with respect to the surface normal with the O atoms weakly coordinating two surface Fe B ions at a distance of 3.00 Å, while in the other mode, the CO 2 is only slightly tilted with respect to the surface normal, where the cations involved are Fe B and Ni B (at 43 Å). For both geometries the calculated adsorption energy is À0.21 eV. The strongest CO 2 -surface interaction using the PW91-D2+U functional (À0.26 eV) was found on surface termination 2, where the CO 2 molecule is also tilted with respect to the surface but with one of the O atoms weakly binding one Fe B surface atom (2.70 Å).
Regardless of which long-range dispersion correction method was applied, calculations with the PBE+U related functionals resulted in the CO 2 molecule recovering its linear conformation on surface termination 2 (Table ESI-2, ESI †). Using PBE-D2+U, we obtained only one geometry with a negative adsorption energy (À0.18 eV), while PBE-D3+U led to an adsorption energy of À0.68 eV due to the reconstruction experienced on the surface upon CO 2 interaction (Fig. 4a). In the latter case, the CO 2 molecule is located parallel to the surface above the top S surface atom, which migrates from the top layer towards the bulk, and it is aligned with the top S layer. Despite the rather large adsorption energy, the CO 2 molecule is not activated but recovers its linear geometry with a CQO bond distance of 1.18 Å, as in the gas phase. Furthermore, the distance between surface and the adsorbate is 3.13 Å. In general terms, the simulated wavenumbers of both the asymmetric and symmetric stretching modes remained very similar to the experimental and computational values reported for the isolated CO 2 molecule (Table ESI-3, ESI †). The weak binding, which does not depend on the functional used, can be rationalised in terms of the oxygen-sulfur repulsion between the molecule's oxygen lone pairs and the electronic clouds of the sulfur surface atoms. This repulsion is quite clear from the CDD plot in Fig. 4b for the geometry obtained with the PBE-D3+U functional, showing a large electrostatic repulsion between O and S atoms, which is the responsible for the large distance between the molecule and the surface. Table 3 Unrelaxed (g u ) and relaxed surfaces energies (g r ) calculated for all the terminations of the FeNi 2 S 4 {111} surface and the % of relaxation Surface termination g u ( J m À2 ) g r ( J m À2 ) Relaxation% Bulk charge  The Grimme dispersion effect on the CO 2 -surface interactions was quantified through comparison of the most accurate dispersion method (D3) with pure DFT (PBE) calculations. Table ESI-4 (ESI †) shows that D3 increases the CO 2 interaction by around 0.14-0.20 eV, slightly lower than has been found in previous works, where the CO 2 interaction with carbides and nitrides was studied. 50,51 3.3.2 CO 2 interaction with the FeNi 2 S 4 {111} surface. Following our predictions for the {001} surfaces, we further investigated the interaction between the CO 2 molecule and termination 4 of the FeNi 2 S 4 {111} surface. Using the PW91-D2+U functional and starting from the same three different pre-activated CO 2 geometries, we found three quasi-degenerate structures with adsorption energies between À0.18 and À0.16 eV (Table ESI-5, ESI †) and molecular geometries similar to that of the most stable interaction with the FeNi 2 S 4 {001} surface. In these structures, the oxygen atoms loosely interact with the surface Fe B atoms at average distances larger than 3 Å, which again confirms the weak interaction with CO 2 . The surface-adsorbate distances and the CO 2 stretching frequencies remain as in the gas phase. A somewhat different result is found when using the PBE-D2+U functional, leading to a moderate adsorption energy of À0.56 eV although the molecule is still practically linear (+OCO 176.31). In this adsorption configuration, the C atom is located on top of a sulfur atom in the top surface layer at a distance larger than 3.5 Å, with the two CQO bonds pointing to Ni A atoms in the top layer. In contrast to the {001} surface terminations, there are no differences in the adsorption energies between the most favourable geometries obtained with PBE-D2+U and PBE-D3+U, both of which release 0.56 eV upon CO 2 adsorption. The wavenumbers of the fundamental vibrational modes for the most favourable adsorption configurations of CO 2 are reported in Tables ESI-3 and ESI-7 (ESI †). 3.3.3 Comparison of CO 2 adsorption to other iron sulfide and oxide materials surfaces. In order to better understand the unexpected small interaction of CO 2 with the different FeNi 2 S 4 surfaces it is useful to compare to the results obtained for other similar systems. To this end, Fig. 5 summarizes adsorption energy values of CO 2 on various sulfides exhibiting different metal/sulfur ratio and also on magnetite, the spinel isomorphic oxide. In line with the data collected in this paper on violarite surfaces, previous CO 2 adsorption studies using the PW91-D2+U functional on the {001} and {111} surfaces of greigite (Fe 3 S 4 ) predicted that the interaction of the CO 2 molecule on the {111} plane was more favourable than on the {001}, although this mineral is unable to activate the CO 2 molecule. 9,15 The relatively weak interaction between CO 2 and the greigite surfaces was assigned to the large electrostatic repulsion between the lone electron pairs of the O atoms of the CO 2 molecule and the large electron clouds of the sulfur anions at the surfaces. 9 However, the present PBE-D2+U results for the interaction of CO 2 with the surfaces of the spinel FeNi 2 S 4 materials clearly differ from those corresponding to the interaction of CO 2 with the {001} surface of Fe 3 O 4 , which also features a spinel crystal structure. Using the same computational methodology, it is predicted that CO 2 is activated by the Fe 3 O 4 {001} surface with a noticeable adsorption energy.
This discussion makes it clear that the origin of the weak interaction cannot be attributed to the crystal structure alone, since violarite, greigite and magnetite all feature the spinel structure. Furthermore, although the greigite and magnetite studies were performed using two different functionals, there is enough evidence to claim that the presence of sulfur atoms at the surface hinders CO 2 adsorption and activation, as clearly summarized in Fig. 5. Nevertheless, when we compare the results for Fe 3 S 4 and FeNi 2 S 4 surfaces obtained using the same PW91-D2+U computational approach (blue bars in Fig. 5), it appears that the incorporation of Ni atoms in the spinel structure favours the interaction with CO 2 , although not sufficiently strongly to activate the molecule.

Conclusions
This work reports a systematic study of the interaction between CO 2 and different terminations and ion arrangements of the {001} and {111} surfaces of the FeNi 2 S 4 spinel violarite, using state of the art DFT-based methods, including dispersion and the on-site Hubbard correction for the 3d levels of the transition metal atoms. The accuracy of the present computational approach is established by comparing to experimental results for the bulk crystal structure. In particular, the PW91-D2+U method provides an optimized value of the lattice parameter in very good agreement with experiment, with values obtained with the PBE-(D2,D3)+U methods being only slightly inferior. The choice of the functional does not play a critical role in the structure of the {111} surface, but in the case of the {001} orientation the most stable surface termination predicted by PW91-D2+U varies from that obtained by PBE-(D2,D3)+U.
The interaction of CO 2 with both violarite surfaces is moderate, although without activation of the molecule and showing only negligible charge transfer from the surface. This weak interaction is attributed to the repulsion between the lone pair electrons of the oxygen atoms of the CO 2 molecule and the spatially extended electronic clouds of the surface sulfur atoms, in agreement with previous findings for the interaction of CO 2 with Fe 3 S 4 surfaces. The substitution of Fe atoms by Ni is found to have a strengthening effect on the binding, but not enough to activate the CO 2 molecule.

Conflicts of interest
There are no conflicts to declare.