Catalytic formation of oxalic acid on the partially oxidised greigite Fe₃S₄(001) surface

Abstract

Greigite (Fe₃S₄), with its ferredoxin-like 4Fe–4S redox centres, is a naturally occurring mineral capable of acting as a catalyst in the conversion of carbon dioxide (CO₂) into low molecular-weight organic acids (LMWOAs), which are of paramount significance in several soil and plant processes as well as in the chemical industry. In this paper, we report the reaction between CO₂ and water (H₂O) to form oxalic acid (H₂C₂O₄) on the partially oxidised greigite Fe₃S₄(001) surface by means of spin-polarised density functional theory calculations with on-site Coulomb corrections and long-range dispersion interactions (DFT+U−D2). We have calculated the bulk phase of Fe₃S₄ and the two reconstructed Tasker type 3 terminations of its (001) surface, whose properties are in good agreement with available experimental data. We have obtained the relevant phase diagram, showing that the Fe₃S₄(001) surface becomes 62.5% partially oxidised, by replacing S by O atoms, in the presence of water at the typical conditions of calcination [Mitchell et al. Faraday Discuss. 2021, 230, 30–51]. The adsorption and co-adsorption of the reactants on the partially oxidised Fe₃S₄(001) surface are exothermic processes. We have considered three mechanistic pathways to explain the formation of H₂C₂O₄, showing that the coupling of the C–C bond and second protonation are the elementary steps with the largest energy penalty. Our calculations suggest that the partially oxidised Fe₃S₄(001) surface is a mineral phase that can catalyse the formation of H₂C₂O₄ under favourable conditions, which has important implications for natural ecosystems and is a process that can be harnessed for the industrial manufacture of this organic acid.

---

1. Introduction

Oxalic acid (H₂C₂O₄), which is the smallest dicarboxylic low molecular weight organic acid (LMWOA),¹ has traditionally been regarded as a major player in the chemical industry^2,3 and natural systems.⁴ Oxalic acid is widely used for the treatment of metallic surfaces^5–8 and textiles,^9–11 as a bleaching agent,^12–15 for the preparation of esters,^16–18 as a reagent in chemical synthesis^19–21 and for separation of rare-earth elements.^22–25 The most sophisticated methods currently used in the industrial production of oxalic acid involve the oxidation of carbohydrates,^26–29 ethylene glycol^30,31 and carbon monoxide.³² Moreover, oxalic acid is the metabolic product of plants,^33–35 soil fungi^36–38 and bacteria,^39–41 and to a lesser extent of animals,^42–44 playing a key role in modulating the interactions between these organisms.^4,45,46 Three main pathways have been suggested for the biosynthesis of oxalic acid, including the cytoplasmatic pathway,⁴⁷ the tricarboxylic acid pathway⁴⁸ and the glyoxylate pathway.^36,49,50

From a catalytic point of view, ferredoxins have gained renewed significance owing to their biological role in electron transfers in several metabolic reactions.^51,52 Ferredoxins are iron–sulfur proteins, with one of the most notable families containing cubane Fe₄–S₄ catalytic centres.⁵³ The catalytic centre of ferredoxins can be either exposed or hidden to solvent, leading to low reduction potentials with a minimum of −600 mV or high potentials with a maximum of +350 mV, respectively.^53,54 Structurally, ferredoxins contain Fe ions in the oxidation states +2 and +3, which are forming tetrahedrally coordinated FeS₄ units, whereas each S atom is bonded to three Fe ions.⁵⁴ Ferredoxin oxidoreductases are responsible for the coupling of the C–C from CO₂ at cellular level in the reductive tricarboxylic acid cycle.⁵⁵ Attempts to synthesise biomimetic Fe₄–S₄ clusters have been successful in controlling the binding of substrate through remote steric effects, but their preparation requires carefully controlled conditions, toxic chemicals and complex steps.⁵⁶

Greigite (Fe₃S₄) is a naturally occurring mineral that contains Fe₄–S₄ clusters with the cubane structure,⁵⁷ very similar to the catalytic centre in ferredoxins. Unsurprisingly, this mineral phase has been found to be catalytically active towards the activation and conversion of CO₂ into small organic molecules, including formic acid, acetic acid, pyruvic acid and methanol.^58–60 Recent studies have suggested that the surfaces of Fe₃S₄ become partially oxidised easily in air or wet conditions, forming core–shell iron sulfide–iron oxide nanoparticles,^61,62 which enhances the catalytic activity.

The important biocatalytic role of ferredoxins and the structural similarities with Fe₃S₄ have led us to study the catalytic conversion of CO₂ and H₂O into H₂C₂O₄ over the partially oxidised Fe₃S₄(001) surface by means of spin polarised density functional theory (DFT) calculations. We have modelled the bulk phase of Fe₃S₄, as well as the two reconstructed Tasker type 3 terminations of its (001) surface, and discuss the changes in several structural, electronic and magnetic properties. We have also calculated the effect of the ratio of the partial pressures of H₂O, i.e. the reactant, and H₂S, i.e. the product of the oxidation, and temperature on the oxygen coverage of the Fe₃S₄(001) surface. Next, we have considered the adsorption of the single CO₂ and H₂O molecules at different sites on the partially oxidised Fe₃S₄(001) surface and assessed a number of descriptors of the chemical activation. We have calculated the structures and energies of possible surface-bound intermediates and saddle points and constructed the reaction profiles for three mechanisms, i.e. via (i) oxalate, (ii) carboxylate and (iii) bicarbonate pathways.

2. Computational methods

2.1. Calculation details

The bulk phase, pristine and partially oxidised (001) surfaces of Fe₃S₄, as well as the reaction mechanisms were calculated using unrestricted DFT calculations as implemented in the Vienna ab initio simulation package (VASP).^63–66 The strongly constrained and appropriately normed (SCAN) functional, which meets all 17 known exact constraints within the meta-generalised gradient approximation (meta-GGA), was used for the exchange-correlation energy functional.⁶⁷ Meta-GGA functionals, which comprise up to the second order derivative of the electron density, i.e. the non-interacting kinetic energy density term, provide energies and structures that are in excellent agreement with experiments.^68–72 The projector augmented wave (PAW) formalism was used to treat the core states of Fe: [Ar], O: [He], C: [He] and S: [Ne], their kinetic energies and the interactions with the valence orbitals.^73,74 The aspherical components of the density gradient were also considered within the PAW one-centre terms, which is required by meta-GGA functionals to describe properly the strongly directional Fe 3d states. The D2 semiempirical method of Grimme was included in our calculations to correct the long-range dispersion interactions,⁷⁵ which is particularly important for an appropriate description of materials and interface properties.^{59,68,76–85} A periodic plane-wave basis set with an upper kinetic energy threshold fixed at 400 eV was employed to expand the Kohn-Sham valence states. The electronic ground state was determined using a preconditioned conjugate gradients minimisation algorithm, which optimises completely and self-consistently the total (free) energy, which is the variational quantity within this methodology.^86–88 This stable and robust electronic minimisation technique is preferred for modelling surface slabs with vacuum gaps when using meta-GGA functionals, as it applies an all band simultaneous update of orbitals. The time step needed to carry out the line optimisation of the energy along the preconditioned conjugate gradient was scaled automatically. The optimisation of the electronic density was terminated when the energy difference between two consecutive self-consistent loop steps was smaller than 10⁻⁵ eV. The internal coordinates were relaxed to their equilibrium positions using a Newton line optimiser, an efficient conjugate gradients technique,^89,90 which was stopped when the Hellmann–Feynman forces were below 0.01 eV Å⁻¹.

We have also used a mean-field Hubbard model, i.e. the so-called DFT+U method, to enhance the calculation of the delocalised Fe 3d states.^91–93 We have applied a simplified rotationally invariant strong intra-atomic Coulomb interaction⁹⁴ and developed the value for the effective parameter U_eff = 5.26 eV following the linear response approach of Cococcioni et al.,⁹⁵ which has been used successfully for the simulation of (Ni,Fe) oxyhydroxides (Ni_1−xFe_xOOH)⁹⁶ and defective ceria CeO₂.⁹⁷ We have carried out a series of static calculations to determine the U_eff value, where we used the conventional cubic unit cell of Fe₃S₄ with ferrimagnetic ordering with the experimental lattice parameter and internal positions.⁵⁷ We broke the symmetry of the tetrahedral (Fe_A) and octahedral Fe (Fe_B) sublattices one at a time, by treating one of these atoms differently to the rest. We calculated the number of d electrons (N_d) for (i) a standard DFT, (ii) a non-self-consistent response function (NSCRF) and (iii) a self-consistent response function (SCRF) to various spherical potentials (V) acting on both channels of the spin of that cation site. The NSCRF and SCRF were found by fitting linearly N_dvs. V, see Fig. SI1 (ESI†). The optimal U_eff value was obtained by calculating the weighted arithmetic mean of the difference between the reciprocals of the gradients of the linear fittings.

2.2 Surface and molecular models

The (001) surface was investigated using the periodic p(1 × 1) supercell, which was created from the bulk Fe₃S₄ using the METADISE code.⁹⁸ A vacuum of 18 Å was added in the direction perpendicular to the surface to prevent interaction with the periodically repeated images and to ensure that the surface can accommodate, in any conformation, the H₂C₂O₄ molecule, which has a length of 4.7 Å. We considered four formula unit (f.u.) layers for the simulation of the (001) surface slab, with an area of 102.517 Ų. The atoms in the two topmost layers were allowed to relax during geometry optimisations, to mimic the surface, while those atoms in the remaining two bottom layers were kept at their optimised bulk positions. The Brillouin zone of the surface was sampled using a Γ-centred Monkhorst-Pack (MP) mesh of 5 × 5 × 1 k-points.⁹⁹ For the geometry optimisations, the partial occupancies were determined using the finite temperature Gaussian smearing method⁸⁷ with 0.01 eV for the width of the smearing, which was tested to ensure that the electronic entropy (TS) vanished. The free energy (F), which is the variational quantity in the finite temperature Gaussian smearing approach differs from E(σ) by TS.^87,100 However, we employed the tetrahedron method with Blöchl corrections¹⁰¹ for the calculation of accurate energies as well as the electronic and magnetic properties in static simulations of the optimised structures. Convergence of energy differences within 1 meV per atom was further tested for the simulation supercell using different vacuum thicknesses, cut off energy values as well as various numbers of total and relaxed atomic layers. Dipole corrections were also included perpendicular to the surface plane to compensate for any dipole created by the chemical species added onto the relaxed side of the slab and to enhance the electronic convergence.^102,103

The isolated CO₂ and H₂O molecules were modelled in a periodic box of 14 × 15 × 16 ų to ensure negligible interaction with their images in the neighbouring cells. The Gaussian scheme was used to determine the electronic partial occupancies during geometry optimisations and energy calculations, sampling only the Γ point of the Brillouin zone.

The effective Bader atomic charges were obtained using an improved grid of charge density values without lattice bias,^104–106 whereas the magnetic moments were integrated within the same atomic basins. The work function (Φ), which is the minimum thermodynamic work required for an electron to leave the Fe₃S₄(001) surface, was calculated as Φ = E_vac − E_F, where E_vac is the potential at the vacuum level and E_F is the Fermi level of the slab. The STM images were calculated using the Tersoff-Hamann method¹⁰⁷ in the implementation of Vanpoucke and Brocks.¹⁰⁸ The wavenumbers of the fundamental harmonic vibrational modes were determined using the central finite differences approach, by allowing each atom to move by a small displacement in the three Cartesian directions. The atomic displacements are used to calculate the elements of the Hessian matrix, which are computed by numerical differentiation of the analytical gradients of the energy with respect to the energy. The saddle points and minimum energy pathways (MEP) between initial and final states were calculated using the climbing image nudged elastic band (cNEB) method.^109,110 The cNEB algorithm used for modelling the MEPs has been employed successfully before to obtain the saddle point structures and energies of several processes in heterogeneous catalysis.^{59,61,68,78,111} We have used five images to model the MEP, which were optimised globally by means of the limited-memory Broyden–Fletcher–Goldfarb–Shannon (LBFGS) method.^89,112 The saddle points were characterised by a single imaginary frequency along the reaction coordinate, which were determined using the central finite differences method.

2.3 Surface energy diagrams

The surface energies for the slab before (γ_u) and after relaxation (γ_r) were evaluated as,where E_u, E_r and E_b are the energies of the unrelaxed slab, the half-relaxed slab and one Fe₃S₄ f.u. in the bulk, respectively, whereas n_b and A are the number of Fe₃S₄ f.u. in the supercell and the surface area, respectively. The degree of relaxation (R) was calculated as R = 100(γ_u − γ_r)/γ_u.

The temperature-corrected energy [E_M(T,p°)] for the adsorbate molecule (M) was calculated according to,

E_M(T,p°) = E_M − TS(T,p°) (3)

where E_M is the calculated DFT energy for the isolated adsorbate, T is the temperature and S(T,p°) is the experimental entropy in the standard state, extracted from thermodynamic tables.¹¹³ The average adsorption energy for the adsorbates [E_ads(T,p°)] was obtained via,where E_c is the energy of the slab with the adsorbed species and n_M is the number of surface-bound adsorbates.

The process of partial oxidation of the Fe₃S₄(001) surface can be represented by the following chemical reaction,

Fe₂₄S₃₂(001) + nH₂O(g) = Fe₂₄S_32−nO_n(001) + nH₂S(g) (5)

where n is the number of O atoms oxidising the surface.

For the partially oxidised surfaces, the effect of temperature on the surface free energy [σ(T,p)] was introduced as follows,^85,114–116

where C is the coverage of oxygen and R is the ideal gas constant. The last term represents the change in free energy of the H₂O and H₂S gases (assuming ideal gas behaviour) in equilibrium with the surface at constant temperature when the ratio of their partial pressures is allowed to change from an arbitrarily small to an arbitrarily large value.

The activation free energy (F_ai) was calculated as the difference between the free energy of the saddle point and the free energy of the reactants, where the index i refers to each of the three elemental steps that we calculated. The free energy of the elemental step (ΔF_i) was obtained as the difference of the free energy of the products and the free energy of the reactants. The energy diagrams were prepared using the most stable binding modes of the reactants, intermediates and final products, which were connected by saddle points and MEPs. The Helmholtz free energy (ΔF) was calculated as ΔF = ΔE − TΔS_vib where S_vib is the vibrational entropy, which was obtained as follows for the adsorbed species,

where ℏ is the reduced Planck constant, ω_i are the vibrational frequencies and k_B is the Boltzmann constant. We have used the entropies reported in thermodynamic tables for the isolated CO₂, H₂O, H₂C₂O₄ and O₂, to avoid calculating the significant translational and rotational components of the entropy for these gas phase molecules.¹¹⁷

3. Results and discussion

3.1. Bulk phase of greigite (Fe₃S₄)

We first investigated the optimised bulk structure of Fe₃S₄, which is later used to create the (001) surface for the calculation of the phase diagrams, molecular adsorptions and reaction mechanisms. Our starting point for Fe₃S₄ is the bulk structure characterised by Li and collaborators using room temperature powder X-ray diffraction (XRD).⁵⁷ Fe₃S₄ has the typical face-centred cubic (fcc) lattice of the spinel structure and crystallises in the space group Fd3̄m (no. 227). Fig. 1 depicts the conventional cubic unit cell of Fe₃S₄ containing 8 f.u. and the periodically repeated images of those ions located at the edges and corners of the cell. The sulfide anions are formally in the 2− oxidation state and exhibit a nearly regular cubic close-packed (ccp) distribution along the [111] direction. The S atoms occupy the tetrahedrally distorted Wyckoff 32e positions with coordinates (u, u, u) and are surrounded by 12 other equidistant nearest-neighbour anions. The difference of the sulfur parameter u from ¼, in direct coordinates, represents the shift of the anions in the [111] direction to accommodate cations of different coordination number and radii in the interstitial sites. Fe₃S₄ is an inverse 2–3 type spinel, as it contains two 3+ and one 2+ cations per formula unit. The tetrahedral Fe_A atoms are in the highest 3+ oxidation state and occupy the 8a crystallographic positions at (1/8, 1/8, 1/8). The octahedral Fe_B cations, which have two equal populations in the 2+ and 3+ valence states each, fill the 16d octahedral holes with coordinates (1/2, 1/2, 1/2). However, only 1/8 of the total number of tetrahedral positions and half of the total number of octahedral holes are filled by the cations. The cations form two different types of alternate layers in the [111] direction, one composed only by Fe_B ions and the other by both Fe_A and Fe_B.

Fig. 1 Polyhedral model of the face-centred cubic (fcc) conventional unit cell containing 8 formula units (f.u.) of Fe₃S₄. Crystallographic directions are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta and S atoms are in light yellow.

Table 1 summarises the calculated and previously reported structural, electronic and magnetic properties for the Fe₃S₄ bulk material. Our calculations reveal that the lattice parameter for the cubic unit cell of Fe₃S₄ was overestimated by 0.25 Å with respect to experiment,⁵⁷ but is in excellent agreement with the value reported by Devey et al. when using PW91+U where U_eff = 5 eV.¹¹⁸ The atomic positions were allowed to relax fully, but the cations remained in their ideal crystallographic sites. However, the value obtained for the sulfur u parameter is just slightly larger than in experiment, indicating that the anions suffer a relatively larger displacement in the simulation cell. The anion–cation bond distances, which are mathematically related to the sulfur parameter by and d(Fe_B–S) = a(1/2 − u), were overestimated by 0.08 and 0.05 Å, respectively, compared to a previous study¹¹⁹ but in agreement with a larger lattice parameter. Our computational setup provides better atomic Bader charges (q) and magnetic moments (m_s) than calculations using small values of U_eff,^118,120 and our results are in better agreement for comparable Hubbard parameters to the one used in this study.¹¹⁸ Although the atomic charges of 1.371e⁻ for Fe_A, 1.189e⁻ for Fe_B and −0.937e⁻ for S are still underestimated with regards to their oxidation states, this is a common finding with Bader charges, and they do provide evidence that our Fe₃S₄ is indeed an inverse spinel as q(Fe_A) > q(Fe_B), which has been observed before.^57,121 The larger magnetic moment of −4.177μ_B atom⁻¹ for Fe_A refers to the high-spin electronic distribution , whereas we can assign the configuration to the antiparallel aligned Fe_B, implying that the conduction mechanism involves both cation sublattices. Our total spin magnetisation of saturation (M_S) of 2.00μ_B f.u.⁻¹ is in excellent agreement with the early works of Morrish and collaborators using Mössbauer spectroscopy.^121,122 Note that a wide range of different values have been reported for Fe₃S₄, corresponding to other electronic and magnetic properties.^{119,123–127}Fig. 2 displays the electronic density of states (DOS) for the bulk phase of Fe₃S₄. The occupied e and t₂ states of the Fe_A ions appear strongly hybridised as a narrow single peak in the minority channel of spins at around −8.5 eV. The Fe_A cations also have a small e band crossing the Fermi level in the α channel, illustrative of the itinerant electron provided by the t_2g valence band of the Fe_B sublattice in the minority channel of the spins. The t_2g and e_g valence bands of the Fe_B cations are located between −7.0 and −3.0 eV in the majority channel of the spins. The virtual e and t₂ states of Fe_A and the t_2g and e_g conduction bands of Fe_B appear above 1.5 eV. Importantly, the hybridisation of the S p levels with the t_2g states of the Fe_B cations observed in the minority spin channel of the DOS of the bulk of Fe₃S₄ explains the small magnetic moment of −0.263μ_B calculated for the anion.

Table 1 Unit cell lattice parameter (a), sulfur parameter (u), interatomic distances (d), atomic charges (q), atomic magnetic moments (m_s) and total spin magnetisation (M_S) for the thiospinel greigite Fe₃S₄ from this work and previous reports. Negative atomic magnetic moments indicate antiparallel alignment

Property	This work	Previous reports	Ref.
a (Å)	10.1251	9.8719	57
u	0.2586	0.2546	57
d(FeA–S) (Å)	2.29	2.21	119
d(FeB–S) (Å)	2.47	2.42	119
q(FeA) (e^− atom^−1)	1.371	1.1	120
q(FeB) (e^− atom^−1)	1.189	1.0	120
q(S) (e^− atom^−1)	−0.937	−0.8	120
m s(FeA) (μB atom^−1)	−4.177	−2.8	120
m s(FeB) (μB atom^−1)	3.614	3.0	120
m s(S) (μB atom^−1)	−0.263	0.1	120
M S(Fe3S4) (μB f.u.^−1)	2.000	2.00	121

---

Fig. 2 Atomic projections of the spin decomposed total density of states (PDOS) for the bulk phase of Fe₃S₄. α and β stand for the majority and minority channel of the spins, respectively. Fe_A and Fe_B contributions are from the 3d bands, whereas S contributions are from the 3p orbitals. Fe_A projections are in black, Fe_B projections are in red and S projections are in green.

3.2. Pristine Fe₃S₄(001) surface

We have also simulated the two non-polar stoichiometric and symmetric terminations of the pristine Fe₃S₄(001) surface, where we have kept the cell borders fixed at the parameters calculated for the bulk, represented in Fig. 3. The stacking of the atomic planes is (Fe_A)–(Fe_B4–S₈)–(Fe_A) and (Fe_B2–S₄)–(Fe_A2)–(Fe_B2–S₄) for the terminations Fe_A and Fe_B–S, respectively, where the atoms within parenthesis are located approximately within the same layers. We have employed single surface unit cells to simulate terminations Fe_A and Fe_B–S of the pristine Fe₃S₄(001) surface, which are the smallest periodically repeating surface unit cells that we can use to model the stacking sequences of these systems. Both terminations Fe_A and Fe_B–S are reconstructed Tasker type 3 surfaces,¹²⁸ where the dipole moment was removed before geometry optimisation by displacing half of the ions from the exposed stoichiometric stacking sequence at the top of the slab to the bottom of the slab. Termination Fe_A has a bulk-like structure comprising a full monolayer (ML) of S atoms, where 25% are 4-fold and the remaining 75% have one dangling bond. The topmost plane of termination Fe_A also contains 0.5 ML of 5-fold Fe_B cations, which form alternating rows in the [110] direction, on top of which are 0.125 ML of 2-fold Fe_A ions with the symmetry, according to the vectorial notation of Wood.¹²⁹ Termination Fe_B–S presents 0.5 ML of 2-fold S atoms, unsurprisingly also forming alternating rows in the [110] direction, alongside 0.25 ML of Fe_B with 3 dangling bonds, whereas the layer beneath contains 0.25 ML of 2-fold Fe_A ions.

Fig. 3 Side (top panels) and top (bottom panels) views of terminations (a) Fe_A and (b) Fe_B–S of the Fe₃S₄(001) surface after geometry optimisation. Layers containing atoms with dangling bonds are highlighted. Crystallographic directions and stacking sequence of the atomic layers before geometry optimisation are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta and S atoms are in light yellow.

Fig. 4 shows the relaxation of the interplanar distances, which was calculated as where d_ij is the separation between the relaxed i and j = i + 1 layers and is the separation for the equivalent planes in the bulk.¹³⁰ Following the geometry optimisation of both terminations, we found a relaxation behaviour that was increasingly attenuated towards the bulk. Each pair of layers Fe_A-i and Fe_B–S-j that were allowed to relax showed a preference to move together in alternating directions. For termination Fe_A, the protruding twofold Fe_A-1 ions suffered a strong inward relaxation of 141%, recovering their bulk tetrahedral coordination, when they merged with the subsurface layer Fe_A-3 where 0.25 ML of fourfold Fe_A are located, see Fig. 3(a) and 4. The negatively charged subsurface plane Fe_B–S-2 experienced only a moderate 13% shift towards the bulk, as it was pulled electrostatically by the merged layer beneath containing a larger number of Fe_A ions and therefore more positive charge carriers than in the bulk. The outward displacement of 16% of the subsurface layer Fe_A-3 can also be rationalised in terms of the Coulomb attraction exerted by the plane Fe_B–S-2 above, which lies closer as a result of its inward relaxation. The subsurface layer Fe_B–S-4 shifted by 9% towards the surface, to compensate for the expansion of the plane Fe_A-3 above, whereas Fe_A-5 preferred to reduce its interplanar distance to the closer Fe_B–S-6 layer that was not allowed to relax. The relaxation pattern of termination Fe_A for the Fe₃S₄(001) surface is in excellent agreement with the relaxation reported for the same termination of the magnetite Fe₃O₄(001) surface.¹¹⁶ Moreover, the contractions and expansions of the interplanar distances Fe_B–S-2 and Fe_A-3 are considerably larger for termination Fe_B–S that than for termination Fe_A, see Fig. 4. However, the inward shift of the topmost layer Fe_B–S-2 of termination Fe_B–S is only half of the contraction of the exposed layer Fe_A-1 of termination Fe_A.

Fig. 4 Relaxation (Δ_ij) of the interplanar distance for the terminations (a) Fe_A and (b) Fe_B–S of the Fe₃S₄(001) surface in black and red, respectively.

From our calculations of the surface energies before and after geometry optimisation, we have found that termination Fe_A is the most stable of the Fe₃S₄(001) planes, see Table 2. The largest relaxation of R = 50.68% obtained for termination Fe_A is a testament of the strong inward displacement of the protruding Fe_A ion to the subsurface layer. We have calculated the atomic Bader charges, which are lower for the exposed Fe_A and S ions than for their bulk counterparts. The (001) surface cations, Fe_A and Fe_B, are more ionic at termination Fe_A than at termination Fe_B–S, whereas we observed the opposite trend for the S anion. We also found an increment by 0.068e⁻ atom⁻¹ in the positive charge of the Fe_B ions in termination Fe_A compared to their bulk values. The magnetic moments calculated for all the under-coordinated atoms are larger in termination Fe_A than in termination Fe_B–S. Our calculations only predict larger magnetic moments than in the bulk for the Fe_B ions with dangling bonds in both terminations, Fe_A and Fe_B–S, and for the exposed S atoms in the termination Fe_A. In contrast to the bulk, we found that the charges and magnetic moments are larger for the exposed Fe_B cations than for the Fe_A ions for the two terminations considered in this study. We also noted that the small magnetic moments of the S atoms, which are aligned parallel to the Fe_A ions in the bulk, prefer to align parallel to the Fe_B cations in both terminations of the Fe₃S₄(001) surface. The work function indicates that termination Fe_A is marginally less reactive, by 0.055 eV, than termination Fe_B–S, which can provide more easily the loosest held electron required for the catalytic formation of H₂C₂O₄. From here, we do not consider the Fe_B–S termination for further analysis due to its large relaxed surface energy, which makes it unlikely to be expressed in the crystal morphology of Fe₃S₄.

Table 2 Surface energies before (γ_u) and after relaxation (γ_r), percentage of relaxation (R), average atomic charges (q), average magnetic moments (m_s) and work function (Φ) for terminations Fe_A and Fe_B–S of the pristine Fe₃S₄(001) surface. Negative atomic magnetic moments indicate antiparallel alignment

Termination	FeA	FeB–S
γ u (meV Å^−2)	66.3	118.7
γ r (meV Å^−2)	32.7	75.2
R (%)	50.68	36.69
q(FeA) (e^− atom^−1)	1.187	1.004
q(FeB) (e^− atom^−1)	1.257	1.022
q(S) (e^− atom^−1)	−0.845	−0.919
m s(FeA) (μB atom^−1)	−3.858	−3.679
m s(FeB) (μB atom^−1)	4.040	3.713
m s(S) (μB atom^−1)	0.318	0.139
Φ (eV)	5.613	5.558

---

3.3. Partially oxidised Fe₃S₄(001) surface

We next analysed the relative stabilities of termination Fe_A of the Fe₃S₄(001) surface containing different coverages of oxygen as a function of the ratio of the partial pressures of H₂O and H₂S and temperature. We have successfully used these ab initio thermodynamic techniques previously to model the oxidation¹¹⁵ and hydrogenation⁷⁸ of cobalt, the hydration of platinum⁸⁵ and the redox behaviour of Fe₃O₄.¹¹⁶ We have calculated several coverages, by replacing in turn each of the S atoms with dangling bonds by O atoms. There is strong experimental evidence, backed up by calculations, that different iron and nickel sulfide phases, including Fe₃S₄,^61,62 troilite (hexagonal FeS),¹³¹ pyrrhotite (hexagonal Fe_1−xS),^68,132 polydymite (Ni₃S₄)⁶⁸ and violarite (FeNi₂S₄)⁶⁸ oxidise partially and easily when they are in contact with water, forming catalytically active core–shell sulfide–oxide nanoparticles. The dimension of our computational slab allows the substitution of up to eight S atoms from termination Fe_A of the Fe₃S₄(001) surface. We have considered several distributions of the S and O atoms for each O coverage and used eqn (6) to calculate their surface free energies. The atomic configurations with the lowest surface free energies for each ratio of the two anions were selected to build the surface phase diagram. The surface free energies for each O coverage are represented by non-planar potential energy surfaces that intersect each other at specific temperatures and ratios of partial pressures of H₂O and H₂S, as shown in a bird's eye view in Fig. SI2 (ESI†). The curves formed where two potential energy surfaces cross correspond to the conditions required to modify the extent of the partial oxidation of the Fe₃S₄(001) surface. However, the most convenient way to depict and analyse this type of diagram is by making a bi-dimensional projection along the surface free energy axis onto the plane formed by the temperature and ratio of partial pressures of H₂O and H₂S, see Fig. 5(a). The surface phase diagram shows that termination Fe_A of the Fe₃S₄(001) surface requires temperatures above 230 K and a ratio of partial pressures of H₂O and H₂S below 100 to remain pristine. The lowest possible oxygen coverage of C = 0.98 O nm⁻² is achieved for a maximum temperature of 300 K and a minimum ratio of partial pressures of H₂O and H₂S of 10^2.2. The coverages 2.93 and 3.90 O nm⁻², with very limited ranges of stabilities and equivalent to three and four O atoms, respectively, appear at higher temperatures and higher ratio of partial pressures of H₂O and H₂S than the 0.98 O nm⁻² coverage. The surface is capable of accommodating five O atoms, i.e. a coverage of 4.88 O nm⁻² if the temperature is increased up to 530 K for a ratio of partial pressures of H₂O and H₂S smaller than 10^5.5. Fig. 5(b) displays this coverage of O atoms, which is consistent with a 62.5% partial oxidation, and is in good agreement with our previous work on FeNi₂S₄, where we found that 75% of the surface S could be replaced by O.⁶⁸ Our results for the Fe₃S₄(001) surface suggest that a coverage of C = 4.88 O nm⁻² has a relatively wide range of stabilities between 300 and 530 K when H₂O is the predominant gas. For comparison, FeNi₂S₄ suffers spontaneous surface oxidation when it is calcined at 473 K in air or moisture.⁶⁸ The coverage is expanded to 6.83 O nm⁻², or seven O atoms, if the ratio of partial pressures of H₂O and H₂S is increased by a further 10^0.5 with respect to C = 4.88 O nm⁻² and the temperature does not exceed 600 K. The full coverage of C = 7.80 O nm⁻² only becomes thermodynamically favourable for a temperature and ratio of partial pressures of H₂O and H₂S larger than the values required for the incorporation of seven O atoms in our surface. Coverages C = 1.95 and 5.85 O nm⁻², corresponding to two and six O atoms, respectively, are not accessible, as they are the least energetically stable partially oxidised Fe₃S₄(001) surface slab models within the conditions of pressure and temperature of our phase diagram.

Fig. 5 (a) Phase diagram for the partially oxidised termination Fe_A of the Fe₃S₄(001) surface as a function of the ratio of the partial pressures of H₂O and H₂S and temperature. The O coverages (C) are indicated for each region. (b) Top view of the Fe₃S₄(001) surface with an oxygen coverage of C = 4.88 O nm⁻². Layers containing atoms with dangling bonds are highlighted. Crystallographic directions and the adsorption sites Fe_A (■), Fe_B (▲), O (×) and S (●) are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta, S atoms are in light yellow and O atoms are in red.

The STM images of both the pristine termination Fe_A of the Fe₃S₄(001) surface and with an oxygen coverage of C = 4.88 O nm⁻² are represented in Fig. 6. The negative sample bias of −2.5 eV used to generate these STM images indicates that the electrons hop from the valence band to the probe tip, as expected for the half-metal properties predicted in the DOS for the bulk of Fe₃S₄. The S atoms appear as the brightest spots in the pristine and oxidised surfaces, whereas the cations are not well resolved. The O atoms can be identified in the oxidised surface, but they are less bright and smaller than the S anions, in agreement with their smaller ionic radius. The pristine surface displays a pattern of well-ordered arrays of S atoms along the [11̄0] direction, which is lost after oxidation.

Fig. 6 Scanning tunnelling microscopy (STM) images of termination Fe_A of the Fe₃S₄(001) surface, where (a) represents the pristine surface and (b) the surface with an oxygen coverage of C = 4.88 O nm⁻². The images were generated using a bias of V = −2.5 eV, a tip distance of d = 2.00 Å, and a density of (a) ρ = 0.0101 e Å⁻³ and (b) ρ = 0.0120 e Å⁻³. Crystallographic directions are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta, S atoms are in light yellow and O atoms are in red.

To provide additional insight into the Fe₃S₄(001) surface, we have also analysed its thermodynamic stability, electronic and magnetic properties, as well as the surface reactivity as a function of the oxygen coverage. The surface free energy (σ) at 0 K, which quantifies the work needed to create the surface of a material and indicates its stability, increases gradually with the oxygen coverage, indicating that partially oxidising the Fe₃S₄(001) facet is an energetically unfavourable process, see Fig. 7(a). Our calculations suggest that the surface free energy increases relatively rapidly for coverages from 0.00 O nm⁻² < C < 1.95 O nm⁻² and especially for 4.88 O nm⁻² < C < 5.85 O nm⁻² and 6.83 O nm⁻² < C < 7.80 O nm⁻². The surface free energy barely changes for the coverages 2.93 and 3.90 O nm⁻², explaining the narrow range of conditions which need to be controlled precisely to develop these particular partial degrees of oxidation. Interestingly, the partially oxidised surface with a coverage of 6.83 O nm⁻² is approximately 13 meV Å⁻² more stable than the facet with C = 5.85 O nm⁻². The atomic charges and magnetic moments do not change noticeably and remain fairly constant for any oxygen coverage, with the exception of the smaller negative charge of the S atom for C = 6.83 O nm⁻² as shown in Fig. 7(b) and (c). The work function shows two minima at the coverages C = 1.95 and 5.85 O nm⁻², suggesting that these partial degrees of oxidation make the Fe₃S₄(001) surface considerably more reactive, see Fig. 7(d).

Fig. 7 (a) Surface free energies (σ) at 0 K, (b) average atomic charges (q), (c) average magnetic moments (m_s) and (d) work function (Φ) for different oxygen coverages of the Fe₃S₄(001) surface. Negative atomic magnetic moments indicate antiparallel alignment.

3.4. Molecular adsorptions on the partially oxidised Fe₃S₄(001) surface

We have investigated the interaction of the single molecules CO₂ and H₂O on the Fe₃S₄(001) surface with a coverage 4.88 O nm⁻², which from hereon is the surface model of our catalyst. In order to simplify the notation used in this paper, we will refer to the Fe₃S₄(001) surface with a 62.5% degree of partial oxidation as “62.5O–Fe₃S₄(001).” We have assessed the effect of these molecular adsorptions on the surface free energy, charge transfer, atomic magnetic moment and work function of the 62.5O–Fe₃S₄(001) surface. We also discuss the geometries of the most favourable adsorption modes and the shifts of the fundamental scaled vibrational modes of the adsorbates. Fig. 5(b) shows the four adsorption sites investigated in this study, including atop the (i) Fe_A, (ii) Fe_B, (iii) S and (iv) O atoms.

Table 3 lists the calculated bond distances (d), bond angles (∠) as well as the harmonic (unscaled) and scaled asymmetric stretching (ν_asym), symmetric stretching (ν_sym) and bending (δ) fundamental vibrational modes for the CO₂ and H₂O adsorbates, which we have used to benchmark our computational SCAN-D2 setup. Our calculations indicate that we have obtained, up to the second decimal, the correct equilibrium bond distances for the CO₂ and H₂O molecules compared to the experimental values.^133,134 The bond angle of 180° predicted for the CO₂ molecule is in perfect agreement with the value inferred from electronic spectroscopy,¹³³ whereas the bond angle calculated for the H₂O molecule is only slightly overestimated with respect to experiment.¹³⁴ Despite their different point groups, i.e. C_2v for H₂O and D_∞h for CO₂, both molecules are characterised by only three fundamental vibrational modes. Compared to previous reports, our calculated asymmetric stretching, symmetric stretching and bending harmonic vibrational modes are red-shifted for CO₂¹³³ and H₂O,^135,136 although the computational description is better for the former adsorbate than for the latter. The smallest deviation from the experimental value was calculated for the bending mode of CO₂ and the largest difference was calculated for the asymmetric stretching mode of H₂O, which were overestimated by 2 and 167 cm⁻¹, respectively. The discrepancies with experiment are larger for H₂O than for CO₂ due to the larger anharmonicity of the vibrations for the hydrogen-containing bonds with respect to the C=O bonds. Our harmonic vibrational frequencies are calculated as the second derivative of the potential energy with respect to the atomic positions in the vicinity of the minimum of the well. This section of the potential energy surface shows quadratic behaviour, leading to harmonic vibrational energy levels that are equally separated. However, experimental fundamental modes represent the transition between the ground state and the first state vibrational energy levels of an anharmonic potential energy surface. We have also determined the empirical scaling factors needed to convert our harmonic frequencies calculated for the isolated molecules into anharmonic vibrational modes, which can be compared directly with experiment.^137–141 The scaling factor (c) was estimated as

where ω_exp are the experimental wavenumbers and ω_calc are the calculated harmonic vibrational wavenumbers. The uncertainty (u) of the scaling factors is calculated asWe obtained the values of c = 0.9655 ± 0.0001 for CO₂ and 0.96089 ± 0.00004 for H₂O, which supports the excellent performance of the meta-GGA SCAN functional for the calculation of vibrational frequencies. After scaling the wavenumbers to account for anharmonicity effects, we found that the asymmetric stretching for both molecules and the symmetric stretching of H₂O were overestimated, whereas the remaining vibrational frequencies were underestimated with respect to the experimental values. Unsurprisingly, the largest difference of 35 cm⁻¹ between scaled and experimental wavenumbers was obtained for the bending mode of H₂O and the smallest shift of only 1 cm⁻¹ was calculated for the symmetric stretching mode of H₂O. In Sections 3.4.1 and 3.4.2 only the scaled vibrational numbers are discussed to facilitate comparison with future experiments.

Table 3 Calculated and previously reported intramolecular bond distances [d(C=O) and d(H–O)], bond angles [∠(O=C=O) and ∠(H–O–H)] and wavenumbers for the fundamental vibrational modes for the isolated CO₂ and H₂O molecules, respectively. The presented harmonic vibrational modes, both unscaled and scaled, are the asymmetric stretching (ν_asym), symmetric stretching (ν_sym) and bending (δ) modes for the O=C=O and H–O–H groups in the CO₂ and H₂O molecules, respectively

		CO2	H2O
d (Å)	Calculated	1.165	0.964
	Experimental	1.162 (ref. 133)	0.958 (ref. 134)
∠ (°)	Calculated	180.00	105.09
	Experimental	180.00 (ref. 133)	104.48 (ref. 134)
ν asym (cm^−1)	Unscaled	2447	3923
	Scaled	2363	3770
	Experimental	2349 (ref. 133)	3756 (ref. 135)
ν sym (cm^−1)	Unscaled	1366	3807
	Scaled	1319	3658
	Experimental	1333 (ref. 133)	3657 (ref. 136)
δ (cm^−1)	Unscaled	669	1623
	Scaled	646	1560
	Experimental	667 (ref. 133)	1595 (ref. 136)

---

3.4.1. CO₂ interaction with the partially oxidised Fe₃S₄(001) surface. We have probed the adsorption of a single CO₂ molecule in a pre-activated configuration, where we reduced its point group from D_∞h to C_2v, by bending the apex angle to ∠(O=C=O) = 130°.^{58,59,76,111,142–146} To characterise the adsorption configurations of our CO₂ molecule, we have borrowed the notation widely used in organometallics to describe the bonding and structural types of metal–CO₂ complexes, which is based on the hapticity of the ligand (μ_n) and the number of surface atoms (m) it can bridge (η^m).^147–149 We approached the bent CO₂ molecule in the μ₃–η³ configuration, with both the C₂ rotational axis and the vertical mirror plane σ_v perpendicular to the surface, i.e. with the C atom coordinating the exposed O or S atoms of the catalyst and the molecular O atoms bonding the Fe cations with dangling bonds. We have also considered the initial η¹-O adsorption modes, where the symmetry elements C₂ and σ_v of the adsorbate lie parallel and perpendicular, respectively, to the surface and the molecular O atom is interacting with the under-coordinated Fe_A or Fe_B ions. We placed the CO₂ molecule at 1.8 Å from the partially oxidised surface of the thiospinel before allowing the atomic positions of the interface to relax.

Our calculations indicate that the CO₂ molecule energetically prefers to interact molecularly with the O ion of the 62.5O–Fe₃S₄(001) surface in the μ₃–η³ configuration, where it releases 0.746 eV, see Table 4 and Fig. 8(a). We also found evidence that the CO₂ molecule binds exothermically to the Fe_A ion in the end-on η¹–O configuration, with an adsorption energy E_ads = −0.188 eV, which is less favourable than at the O site. Despite starting in the μ₃–η³ configuration, the optimised structure of the adsorption of CO₂ on the S site is the side-on μ₂–η²–(C,O) mode, which is a slightly endothermic process requiring 0.016 eV. The η¹–O interaction configuration with the exposed Fe_B atoms of the 62.5O–Fe₃S₄(001) surface displays the largest positive energy at E_ads = 0.589 eV. The trend of binding energies on the anion sites can be rationalised in terms of the number of interactions that the CO₂ molecule forms with the surface, i.e. three in the μ₃–η³ configuration on the most exothermic O position, and two in the μ₂–η²–(C,O) mode on the S atom. However, the trend of the adsorption strength on the cation sites depends on their total number of dangling bonds, i.e. two and one for the Fe_A and Fe_B atoms, respectively. Quesne and collaborators reported that chemisorption of the activated CO₂ molecule on the surfaces of transition metal carbides occurs either directly, spontaneously and barrierless, or via a linearly physisorbed intermediate with a small energy barrier.¹⁵⁰ Both the chemisorption and physisorption processes of CO₂ on the surfaces of transition metal carbides are exothermic, with the former releasing more energy than the latter.¹⁵⁰ Thus, regardless of the adsorption mechanism of CO₂ on the partially oxidised 62.5O–Fe₃S₄(001) surface, we expect this process to take place on the O site, given its relatively large exothermic adsorption energy of −0.746 eV.

Table 4 Structural type, adsorption energies (E_ads) at 0 K, average interatomic distances (d), bond angle (∠), scaled wavenumbers for the fundamental vibrational modes and charge transfers (Δq) calculated for a single CO₂ molecule interacting with the different adsorption sites considered for the partially oxidised 62.5O–Fe₃S₄(001) surface. The presented vibrational modes are the asymmetric stretching (ν_asym), symmetric stretching (ν_sym) and average bending (δ) modes. The surface free energies (σ) at 0 K and work functions (Φ) are also reported. Negative values of Δq indicate charge transfer from the surface to the adsorbate

Adsorption site	O	S	FeA	FeB
Structural type	μ 3–η^3	μ 2–η^2–(C,O)	η ^1–O	η ^1–O
E ads (eV)	−0.746	0.016	−0.188	0.589
d(Omol–Fe) (Å)	2.047	1.997	2.308	2.157
d(C–Osurf) (Å)	1.353	—	—	—
d(C–Ssurf) (Å)	—	1.899	—	—
d(C=Omol) (Å)	1.277	1.255	1.171	1.177
∠(O=C=O) (°)	127.98	130.47	177.46	172.67
ν asym(O=C = O) (cm^−1)	1520	1690	2329	2272
ν sym(O=C=O) (cm^−1)	1225	1088	1297	1278
δ(O=C=O) (cm^−1)	795	671	581	594
Δq (e^−)	−0.209	−0.537	0.003	0.003
Φ (eV)	5.172	5.687	4.662	4.221
σ (meV Å^−2)	65.1	72.5	70.5	78.1

---

Fig. 8 Molecular adsorption of (a) CO₂ and (b) H₂O on the O and Fe_A, Fe_A sites, respectively, of the partially oxidised 62.5O–Fe₃S₄(001) surface. Side (top panels) and top (bottom panels) views are displayed. Layers containing atoms with dangling bonds are highlighted. Crystallographic directions are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta, S atoms are in light yellow, O atoms are in red, H atoms are in white and C atoms are in black.

The interactions have a noticeable impact on the geometry of the interface, the vibrational properties of the CO₂ molecule, as well as on the electronic properties and thermodynamic stability of the partially oxidised surface of our catalyst. The largest interfacial O–Fe distances, with the values d = 2.308 and 2.157 Å were calculated for the η¹–O structures on the Fe_A and Fe_B sites, respectively, indicating that the only coordination bond for these adsorption modes is weak. The CO₂ molecule is more strongly attached to the anion sites, and particularly to the surface O atom, since the distances for both the interfacial secondary O–Fe and main C–S (C–O) coordination bonds are the smallest overall. The intramolecular C=O bond distances, which experience minor elongations for the η¹–O interactions with the two types of Fe cation sites, are stretched by 0.090 and 0.112 Å for the μ₃–η³ adsorption on the O and μ₂–η²–(C,O) configuration on the S positions, respectively. Despite introducing a pre-activated bent CO₂ molecule, which remained in this geometry on the anion sites during geometry optimisation, we observed that the adsorbed molecule became approximately linear on the cation positions, similar to the molecule in its isolated state. This has also been found in other iron sulfide and oxide catalytic materials, where a pre-activated molecule is introduced but allowed to relax unconstrainedly during geometry optimisation.^{58,60,111,146}

The changes observed in the intramolecular C=O bond distances are directly proportional to the red-shift of the asymmetric stretching vibrational modes compared to the isolated adsorbate, with the smallest value of 1520 cm⁻¹ obtained for the ground state μ₃–η³ binding configuration on the surface O site. Our results suggest that the variations in the symmetric stretching and average bending vibrational modes cannot be rationalised in terms of the reduction of the intramolecular C=O bond length. As expected, the symmetric stretching vibrational modes for the CO₂ molecule adsorbed on the anion sites are smaller than on the cation positions, but the largest wavenumbers at each type of site were calculated for the most stable adsorption modes on the O and Fe_A positions. The interaction of CO₂ with the surface O and S atoms leads to a blue-shift of the average bending vibrational mode with regard to the isolated adsorbate, whereas the η¹–O coordination to the Fe sites causes a red-shift. Table SI1 (ESI†) summarises both the unscaled and scaled vibrational frequencies calculated for the adsorbed CO₂ molecule.

The charge transfers denote that the cation sites, which are deficient in electrons, receive a minor charge Δq = 0.03e⁻ from the adsorbate. However, the electron-rich O and S positions donate −0.209 and −0.537e⁻, respectively, to the σ* antibonding molecular orbital of the CO₂ molecule, which weakens the intramolecular C=O bond, forcing a bent configuration. The adsorption of the CO₂ molecule reduces the work function, which makes the 62.5O–Fe₃S₄(001) surface of our catalyst more reactive, except for the μ₂–η²–(C,O) interaction with the S site. The lowest work function values Φ = 4.662 and 4.221 eV required to remove the loosest held electron were predicted when the adsorbate binds the Fe_A and Fe_B sites, respectively, as the surfaces have increased marginally their electron density after adsorption of CO₂. The exothermic interactions of a single CO₂ molecule enhance the thermodynamic stability of our 62.5O–Fe₃S₄(001) surface, with the lowest value of the surface free energy σ = 65.087 meV Å⁻² computed for the μ₃–η³ adsorption on the O site.

Previous studies have shown that the (001) and (111) surfaces of the thiospinels Fe₃S₄^58,111 and FeNi₂S₄⁷⁶ are unable to activate the CO₂ molecule, which only remains physisorbed. Among the pure sulfides, the (111) surface of mackinawite (tetragonal FeS) is the only catalyst able to dissociate a chemically activated CO₂ molecule into stable fragments.^111,146 The Fe₃O₄(001) and (111) surfaces are also capable of strongly chemisorbing the adsorbate, but the dissociation reaction has large activation energies, while the CO and O fragments are at least 1.5 eV less stable than the CO₂ molecule.¹¹¹ To date, the catalytic activity towards the reduction of bicarbonate has been reported for the major surfaces of the partially oxidised sulfides Fe₃S₄,^58,62 hexagonal Fe_1−xS,^68,132 Ni₃S₄⁶⁸ and FeNi₂S₄,⁶⁸ with the latter described as the most efficacious system due to the synergism between Fe and Ni and the O content. Our DFT calculations show that the activation of the adsorbed CO₂ molecule can be enhanced via partial oxidation of the least reactive Fe₃S₄(001) surface.

3.4.2. H₂O interaction with the partially oxidised Fe₃S₄(001) surface. In this section, we focus on discussing the molecular and dissociative adsorption of a single H₂O molecule on the partially oxidised 62.5O–Fe₃S₄(001) surface, which are competing processes that have been reported on the surfaces of oxides,^151–159 sulfides^160–164 and metals.⁸⁵ For the molecular binding modes, we initially placed the adsorbate in the μ₃–η³ configuration, with both symmetry elements C₂ and σ_v perpendicular to the surface of our catalyst, i.e. with the molecular O atom coordinating the under-coordinated cations and the H atoms interacting with the exposed anions. We also explored the possibility that the C₂ axis and σ_v plane of the H₂O molecule are oriented parallel and perpendicular to the surface, respectively, with the H atoms forming hydrogen-bonds with the surface O anions in the η¹–O structure. For the dissociative adsorption configurations, we introduced one OH group coordinating to the exposed Fe_A and Fe_B ions and forming hydrogen-bonds with the surface anions in the μ₂–η²–(O, H) configuration and with the dissociated H binding a nearby O atom. The O (H) interacting atoms from the adsorbate were located at 1.8 (1.00) Å from the surface of our catalyst before carrying out a full geometry optimisation. However, we have not investigated the η¹–O adsorption configurations via H-bonds, as the H₂O rotated during geometry optimisation leading to the same binding modes considered for the surface Fe ions.

Table 5 lists the exothermic adsorption energies calculated for H₂O, which are clearly more favourable than for the interaction of the CO₂ molecule with the 62.5O–Fe₃S₄(001) surface. The dissociative adsorption modes release 0.507 and 0.219 eV less than the molecular binding configurations at the (Fe_A, Fe_A) and Fe_B sites, respectively, suggesting that H₂O prefers to remain undissociated at the surface of our catalyst. The order of decreasing stabilities are E_ads (Fe_A, Fe_A) < E_ads (Fe_A) ≪ E_ads (Fe_B) and E_ads (Fe_A, Fe_A) ≪ E_ads (Fe_B) for the molecular and dissociative binding modes, respectively, showing that the adsorption energies have a strong dependence on the type of Fe site. The trend of adsorption energies can be explained in terms of the total number of dangling bonds of the cation sites. The most stable Fe_A, Fe_A binding site with the μ₄–η⁴ configuration has a total of three dangling bonds, see Fig. 8(b), whereas the Fe_A and Fe_B sites with the μ₃–η³ structure have two and one dangling bonds, respectively.

Table 5 Structural type, adsorption energies (E_ads) at 0 K, average interatomic distances (d), bond angle (∠), scaled wavenumbers for the fundamental vibrational modes and charge transfers (Δq) calculated for a single H₂O molecule interacting with the different adsorption sites considered for the partially oxidised 62.5O–Fe₃S₄(001) surface. The presented vibrational modes are the asymmetric stretching (ν_asym), symmetric stretching (ν_sym) and average bending (δ) modes. The surface free energies (σ) at 0 K and work functions (Φ) are also reported. Negative values of Δq indicate charge transfer from the surface to the adsorbate

Adsorption site	FeA	FeA, FeA		FeB
Type	Molecular	Molecular	Dissociative	Molecular	Dissociative
Structural type	μ 3–η^3	μ 4–η^4	μ 3–η^3	μ 3–η^3	μ 3–η^3
E ads (eV)	−1.094	−1.318	−0.811	−0.510	−0.291
d(Omol–Fe) (Å)	2.151	2.301	2.060	2.089	1.909
d(H1–Osurf) (Å)	2.436	2.800	2.704	—	—
d(H2–Osurf) (Å)	—	—	0.980	1.466	1.023
d(H1–S) (Å)	—	—	—	2.599	2.753
d(H2–S) (Å)	2.242	1.963	2.650	—	—
d(H1–Omol) (Å)	0.981	0.977	0.973	0.979	0.974
d(H2–Omol) (Å)	1.001	1.044	2.728	1.069	1.639
∠(H–O–H) (°)	107.73	106.65	175.77	106.77	110.15
ν asym(H–O–H) (cm^−1)	3461	3549	3602	3511	3578
ν sym(H–O–H) (cm^−1)	3053	2339	3356	2003	2585
δ(H–O–H) (cm^−1)	1514	1468	807	1514	904
Δq (e^−)	−0.005	−0.032	0.020	0.010	0.041
Φ (eV)	5.660	4.996	5.504	5.623	5.405
σ (meV Å^−2)	61.7	59.5	64.5	67.4	69.5

---

The average interfacial O–Fe binding distances for H₂O, which are marginally smaller than in the interaction of CO₂, are inversely proportional to the adsorption energies. The smallest distances d(O_mol–Fe) = 1.909 and 2.060 Å were calculated for the thermodynamically least stable dissociative adsorptions on the Fe_B and (Fe_A, Fe_A) sites, respectively, as the OH group is a better ligand than H₂O. Our calculations indicate that the molecular and dissociated H₂O molecules are able to form hydrogen-bonds with the surface O and S atoms. The interfacial hydrogen-bonds with S are stronger than with O for the interactions with the Fe_A and (Fe_A, Fe_A) sites, as evidenced by the H–S distances, which are smaller than the H–O_surf lengths. We also found that the hydrogen-bond lengths with the surface O atoms are directly proportional to the adsorption energies for the molecularly adsorbed H₂O molecule, with the largest value d(H–O) = 2.800 Å calculated for the ground state μ₄–η⁴ adsorption mode. However, the opposite effect of the binding energies was observed for the hydrogen-bond lengths with the S atoms, with the smallest value d(H–S) = 1.963 Å obtained for the most stable μ₄–η⁴ binding configuration. The H–O_surf distance of 2.704 Å suggests that the OH group adsorbed on the (Fe_A, Fe_A) site in the μ₃–η³ structure is able to form a hydrogen-bond with the surface, which is not observed for the dissociative adsorption mode on the Fe_B site. Furthermore, the dissociated H atom prefers to sit slightly closer at 0.043 Å to the surface O atom for the interaction with the (Fe_A, Fe_A) site than for the Fe_B site, in agreement with their relative adsorption energies. The adsorption of H₂O induces minor changes on the intramolecular H–O distances, which are between 0.009 and 0.105 Å larger for the two types of adsorption modes than in the isolated adsorbate. Our calculations show that the dissociated H atom diffuses 1.089 Å further away from the O of the OH group for the interaction with the most stable (Fe_A, Fe_A) site than with the least favourable Fe_B position. Note that we did not find evidence of a dissociated H atom binding to a nearby surface S atom. The intramolecular bond angle is less sensitive to adsorption of H₂O than CO₂, with the largest value ∠ = 107.73° calculated for the molecular μ₃–η³ adsorption on the Fe_A site. However, this geometrical descriptor can be used to discriminate the dissociative from the molecular adsorption modes of H₂O, as the intramolecular bond angle values are typically larger for the former than for the latter.

The three fundamental vibrational modes become smaller as a result of the molecular and dissociative adsorption of H₂O. The asymmetric stretching vibrational mode is below 3549 cm⁻¹ for the molecular adsorption configurations, whereas it lies above 3578 cm⁻¹ for the dissociated OH groups. The symmetric stretching vibrational modes, which appear in a wide range of values between 3356 and 2003 cm⁻¹, show no clear correlation with any of the properties discussed. The bending vibrational modes were estimated between 1514 and 1468 cm⁻¹ for the molecular adsorption modes and the calculated values were smaller than 904 cm⁻¹ for the dissociative binding configurations. Table SI2 (ESI†) displays both the harmonic and scaled vibrational modes observed for the adsorbed H₂O molecule. We calculated only minor charge transfers, with the H₂O molecule receiving no more than 0.032e⁻ for the most stable molecular adsorptions on the Fe_A and (Fe_A, Fe_A) sites. Our results indicate that, unlike the linear CO₂ molecule, the activation and adsorption of the angular H₂O species does not require charge density from the surface. The work functions are Φ ∼ 5.5 eV for all the adsorption modes, except for the molecular μ₄–η⁴ binding configuration on the (Fe_A, Fe_A) site, which is around 0.5 eV smaller than the value for the 62.5O–Fe₃S₄(001) surface. The large reduction of the surface free energy after adsorption of the H₂O molecule demonstrates that the 62.5O–Fe₃S₄(001) surface becomes more stable, even compared to the catalyst interacting with CO₂. In line with our results, previous calculations have also shown that the molecular adsorption of water is energetically preferred over the dissociative mode for CaO,¹⁵⁶ CaF₂^156,165 and the Fe₃S₄(001) surface.¹⁶¹ Moreover, the partially oxidised 62.5O–Fe₃S₄(001) surface also displays an improved activity towards the adsorption of H₂O than the (001) facet of pure Fe₃S₄.

3.5 H₂C₂O₄ formation on the partially oxidised Fe₃S₄(001) surface

We have examined three pathways for the catalytic conversion of CO₂ and H₂O into H₂C₂O₄ on the 62.5O–Fe₃S₄(001) surface, see Fig. 9. Our selection of mechanisms, which describe the combination of the reactants into key intermediates that are transformed into the desired H₂C₂O₄ and by-product O₂, are discussed in terms of their structural, local electronic and vibrational properties. The Helmholtz free energies (ΔF), which are plotted at the representative temperature of 300 K, are referenced in our energy diagrams compared to the isolated reactants and the heterogeneous catalyst, i.e. two CO₂ molecules, one H₂O molecule and the 62.5O–Fe₃S₄(001) surface slab. Note that the vibrational frequencies of the adsorbed species were scaled by 0.9632, which is the average of the scaling factors calculated for the isolated CO₂ and H₂O molecules, before estimating the vibrational entropy of the adsorbed species, see Tables SI3 and SI4 (ESI†). For the purpose of comparison, the Helmholtz free energy profiles at 0 and 600 K are provided in Fig. SI3 and SI4 (ESI†). The first step, which is common to the three pathways, is the co-adsorption of the reactants to the surface of the catalyst. Note that the adsorbed species are indicated with the symbol *.

Fig. 9 Minimum energy pathways (MEPs) for the conversion of CO₂ and H₂O into H₂C₂O₄ on the partially oxidised 62.5O–Fe₃S₄(001) surface at 300 K. Minimum states are denoted by bold lines and saddle point are represented by narrow lines linked by dashed lines. Energies are referenced to the isolated 2CO₂(g) + H₂O(g) molecules, including the surface slab. Adsorbed species are presented followed by the symbol *.

3.5.1. Oxalate pathway. In line with our observations for the μ₃–η³ adsorption of the single molecules, the initial co-adsorption configuration comprises CO₂ interacting with two nearby surface O atoms and four Fe cations, whereas the H₂O species was introduced coordinating the (Fe_A, Fe_A) site and forming two hydrogen-bonds in the μ₄–η⁴ structure. We had assumed that the (co-)adsorption of the reactants is barrierless, but found that this is actually an exothermic process that releases ΔF = 417 meV and is therefore 2.39 eV less favourable than the sum of the adsorption energies of the single molecules, as shown in Fig. 9. The co-adsorbed CO₂ molecules interact with the O sites at 1.340 and 1.399 Å, and with the C₂ axis remaining perpendicular to the surface only for the closest adsorbed species, see Fig. 10(a). Both, the perpendicular and tilted CO₂ molecules are able to coordinate to the surface Fe_A and Fe_B cations at an average distance of ∼2.1 Å, but the former uses its two O atoms, while the latter employs only one O atom, which explains their different relative orientations with respect to the surface. Our calculations suggest that both CO₂ molecules become activated upon adsorption, with the apex angle ∠ at < 126.34°, which is already smaller than the value calculated for the most stable μ₃–η³ adsorption mode of the single molecule. The intramolecular bonds in the CO₂ molecules are elongated to an average value of d(C=O) ∼ 1.28 Å, in agreement with the single molecule adsorption. The H₂O species sits on the (Fe_A, Fe_A) site at approximately the same distance calculated for the most stable single molecule μ₄–η⁴ adsorption, with the intramolecular distances and bond angle remaining relatively undisturbed. However, the surface S and O atoms forming the hydrogen-bonds are 0.120 Å further away and 0.262 Å closer to the H₂O molecule, respectively, than in the ground state single molecule adsorption configuration. The charge analysis of the co-adsorption state indicates that the perpendicularly adsorbed CO₂ molecule received 0.027e⁻ less charge than in the case of the single molecule adsorption, which was compensated for by the identical electron density gained by the tilted species. On the other hand, the H₂O adsorbate received only a negligible charge of Δq = 0.002e⁻, which is considerably less than in the case of the single molecule μ₄–η⁴ adsorption. We attribute the smaller charges received by the H₂O and perpendicular CO₂ molecules to the fact that they share the same surface Fe_A ion.

Fig. 10 Side views of the adsorption configurations of (a) reactants 2CO₂* and H₂O*, intermediates (b) C₂O₄* and H₂O* and (c) HC₂O₄* and OH*, as well as (d) final products H₂C₂O₄* and O* of the oxalate pathway on the partially oxidised 62.5O–Fe₃S₄(001) surface. Interatomic distances are provided in Ångstroms; crystallographic directions are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta, S atoms are in light yellow, O atoms are in red, H atoms are in white and C atoms are in black.

The oxalate pathway leads to the combination of the two CO₂* molecules into the C₂O₄* intermediate, which is an endothermic process with a reaction free energy of 2.100 eV and a saddle point A at 3.296 eV, as shown in Fig. 9. The large free energy barrier of this process is required to break the strong covalent bonds within the adsorbed carbonate-like groups, leading to the C atoms moving to at least 3.0 Å away from the nearest surface O atoms, see Fig. 10(b). The product of this elementary step, which is the least stable intermediate, is a staggered C₂O₄* group with the dihedral angle ∠(O–C–C–O) = 66° and the point group of reduced symmetry D_2d. The C₂O₄* species is anchored to the exposed Fe cations of the surface through three O atoms, with d(O–Fe) ∼ 2.0 Å, which is slightly smaller than for the co-adsorption state. The C–C bond distance is 1.539 Å for the C₂O₄* intermediate, whereas the intramolecular C=O distances and O=C=O angles do not change noticeably with respect to the previous state. The calculated Bader charges show that the C₂O₄* species draws −1.367e⁻ from the surface, which is considerably larger than for two co-adsorbed CO₂ molecules, explaining the low stability of this intermediate.

The C₂O₄* species accepts one H atom from the H₂O molecule to form HC₂O₄* in state 4 of our energy diagram, see Fig. 9. The free energy released during this elementary step is 0.846 eV, but it requires 1.679 eV to cross the saddle point B. The protonated O atom, with the typical distance of 0.98 Å for the O–H bond, is the closest to the donor H₂O molecule, belonging to the bidentate CO₂ moiety in the previous state, as shown in Fig. 10(c). The dihedral angle becomes more orthogonal with ∠(O–C–C–O) = 81.82° after protonation, whereas the largest intramolecular distances of d(C–O) ∼1.31 Å were calculated in each molecular moiety for the protonated O and the O coordinating the surface Fe ion. The OH* fragment moves 0.442 Å closer to the Fe_A ion that is shared with the HC₂O₄*, as the basicity of the former increases and the latter decreases, compared to H₂O* and C₂O₄* respectively. However, the shift of the OH* species happens at the cost of an increase in the H⋯S hydrogen bond-distance by 0.383 Å. The protonation reduces almost by half the charge of the HC₂O₄* intermediate to q = −0.704e⁻ with respect to C₂O₄*, whereas the OH* is able to retain −0.650e⁻ from the proton that it ceded.

Fig. 9 illustrates that the HC₂O₄* intermediate is capable of admitting a second proton from the OH* group, which is also an endothermic process with a reaction free energy of 2.202 eV and an activation free energy C of 4.027 eV. The protonation increases the symmetry of H₂C₂O₄* to point group C_2v, implying that the molecule becomes essentially flat with the dihedral angle ∠(O–C–C–O) = 164°, see Fig. 10(d). The intramolecular bond distances d(C–C) = 1.548 Å, d(C = O) ∼1.22 Å, d(H–O) ∼1.00 Å and d(C–OH) ∼ 1.341 Å and angle ∠(O–C–O) = 125.8° of H₂C₂O₄* are very similar to the values calculated for the C₂O₄*. Our calculations suggest that H₂C₂O₄* is still coordinating one Fe_A and two Fe_B, with the smallest Fe–O distance of 2.248 Å for the unprotonated O and the largest distances of 2.92 Å for the protonated O atoms. The O* species moves approximately 0.2 Å towards the bulk after losing its H atom. The transfer of the second proton quenches the negative charge of H₂C₂O₄* to −0.037e⁻, whereas O* has −1.108e⁻, almost duplicating the electron density of OH* in the previous state.

In the final step, the surface releases one H₂C₂O₄ and 0.5 O₂ molecules, leaving the system ready for the next catalytic cycle. These desorptions are endothermic processes that require 0.477 eV. The final state is the overall highest in free energy in our energy diagram, lying 3.518 eV above the reactants. To calculate these energies, we modelled the isolated O₂ molecule in the triplet state and the isolated H₂C₂O₄ in the point group C_2h, which are the electronic and structural ground states, respectively, of these species.

3.5.2. Carboxylate pathway. The protonation of one of the CO₂* species in the carboxylate pathway leads to the most stable intermediate in this study, which is 1.157 eV below the reactants in our free energy diagram, see Fig. 9. This exothermic process has a reaction free energy of −0.740 eV and a transition state D at 1.324 eV. The O–H bond distance is 1.01 Å for the bicarbonate-like group, which also forms a hydrogen-bond of 1.727 Å to the neighbouring carbonate-like species, as shown in Fig. 11(a). The intermolecular hydrogen-bond is enabled by the rotation of the two intermediates around the C₂ axis perpendicular to the surface, which also reduces their interfacial C–O distance by an average of 0.05 Å and their apex O–C–O angles by at least 5°. The protonation weakens the C–OH bond, as its distance increases by 0.1 Å, which implies that part of the electron density is shared with the newly added H atom. We found that the S atom that formed the hydrogen-bond to the H₂O molecule, moved to 1.751 Å away from the CO₂* group during the proton transfer. The Bader analysis assigns 0.445e⁻ to the CO₂H*, which receives a large Coulomb attraction from the negatively charged surface O atom, explaining the stability of these intermediates. Despite the CO₂* species remaining bent, it donates electron density back to the surface S atom and is only able to keep −0.032e⁻. Our calculations show that the OH* species displays very similar Fe–O and O–H bond distances and electron density charge as the OH* group that coexists with HC₂O₄* in the oxalate pathway.

Fig. 11 Side views of the adsorption configurations of the intermediates (a) CO₂H*, CO₂* and OH* and (b) 2CO₂H* and O* of the carboxylate pathway, as well as intermediates (c) CO₃H*, CO₂* and H* and (d) CO₃H* and CO₂H* of the bicarbonate pathway on the partially oxidised 62.5O–Fe₃S₄(001) surface. Interatomic distances are provided in Ångströms; crystallographic directions are indicated. Fe_A atoms are in dark yellow, Fe_B atoms are in magenta, S atoms are in light yellow, O atoms are in red, H atoms are in white and C atoms are in black.

Fig. 9 depicts the protonation of the second CO₂* species, which is an endothermic process with a reaction free energy of 1.547 eV and a saddle point E at 3.066 eV. We found that the structural properties of the two CO₂H* species, such as the H–O and C–O bond distances and the O–C–O angle are very similar to the protonated group of the previous state, as shown in Fig. 11(b). However, our calculations show that each CO₂H* group can only coordinate a single Fe_B cation via the unprotonated O atom at the slightly different distances of 2.058 and 2.245 Å. The two CO₂H* species have a positive charge of ∼0.5e⁻, whereas the O* gains approximately half of one electron.

The combination of the two carboxylate groups to form the adsorbed H₂C₂O₄* molecule is an endothermic process with a reaction free energy of 2.650 eV, which is 0.448 eV larger than the free energy required to transfer the second proton to the HC₂O₄* in the oxalate pathway, see Fig. 9. The free energy of the saddle point F is 3.998 eV for the formation of the C–C bond in the carboxylate pathway, whereas it is ∼0.7 eV smaller for the saddle point A in the oxalate pathway.

3.5.3. Bicarbonate pathway. In the bicarbonate pathway, the H₂O* molecule dissociates and the H* atom migrates to a nearby surface O atom, whereas the OH* group diffuses to one of the CO₂* species. This is an endothermic process with a reaction free energy of 0.699 eV, i.e. the lowest in this study, and a transition state G at 2.318 eV, shown in Fig. 9. The CO₃H* species moves 1.921 Å away from the surface upon hydroxylation, remaining only bound to the undercoordinated Fe ions at the average distance of 2.0 Å in the μ₂–η²–(O,O) configuration, see Fig. 11(c). The C–OH distance is 1.330 Å in the CO₃H* species, which tilts to allow the formation of a short hydrogen-bond of 1.764 Å with the CO₂* molecule. The interfacial Fe_B–O and C–O bonds between the CO₂* molecule, which prefers to stay in the side-on μ₂–η²–(C,O) structure, and the surface are 2.000 and 1.413 Å, respectively. The intramolecular C=O bond distance is 1.28 Å for the CO₃H* species and 1.300 and 1.233 Å for the O end coordinating the surface Fe_B and the free O end, respectively, of the CO₂* molecule. The dissociated H* atom, which stays coordinated to the S atom that was part of the hydrogen-bond, diffuses to the subsurface layer to increase its stability. The incorporation of the OH group into the CO₃H* species increases its negative charge to −0.672e⁻ with respect to the value of −0.213e⁻ for the CO₂* molecule, whereas the subsurface H* atom becomes almost neutral with 0.053e⁻.

Fig. 9 depicts the protonation of the second CO₂* molecule in our free energy profile. This elementary step releases −0.384 eV, which is the smallest exothermic free energy value reported in this study, and has a saddle point of 2.120 eV. The newly formed CO₂H* species remains bound to a surface O anion at 1.295 Å and to a Fe_B ion at 1.755 Å through the molecular C and OH group, respectively, as shown in Fig. 11(d). Our calculations suggest that the intramolecular OH distance is 0.976 Å. However, we found that the proton diffusion only leads to negligible changes in the structure of the CO₃H* species, including the length of the Fe–O and hydrogen-bonds with the surface and CO₂H* intermediate, respectively. The protonation provides the CO₂*H intermediate with a positive charge of 0.382e⁻, whereas CO₃H* increases its negative charge by 0.038e⁻.

Finally, the coupling of the CO₃H* and CO₂H* species is an endothermic elementary step with a reaction free energy ΔF = 3.143 eV and the largest activation energy value of 5.177 eV calculated in this work, see Fig. 9.

4. Conclusions

We have used DFT methods to model the catalytic conversion of CO₂ and H₂O into H₂C₂O₄ on the partially oxidised Fe₃S₄(001) surface. First, we have modelled the bulk phase of Fe₃S₄ and found that the structural, electronic and magnetic properties are in good agreement with previous reports. We have elucidated the relaxation patterns of the interplanar distances, as well as the values for the surface energies, atomic charges, atomic magnetic moments and work functions for the two reconstructed Tasker type 3 terminations of the Fe₃S₄(001) surface. The thermodynamically most stable facet, which terminates in two-fold Fe_A ions with symmetry, has the largest relaxation, work function, ionic character and magnetisation. We have predicted the phase diagram as a function of the ratio of partial pressures of H₂O and H₂S and temperature and found that the processes of replacing sequentially each S atom with dangling bonds by an O atom are endothermic. Only selected coverages of O are allowed in the Fe₃S₄(001) surface, with 62.5% of partial oxidation becoming prominent at the typical experimental conditions in which samples of the catalyst are calcined.

We have also reported the interaction between single molecules of CO₂ and H₂O and the Fe₃S₄(001) surface with 62.5% of partial oxidation. Adsorption at the O site, which is energetically preferred, bends and activates the CO₂ molecule and the σ* antibonding molecular orbital receives electronic density from the partially oxidised Fe₃S₄(001) surface. The H₂O molecule releases the largest adsorption energy when it coordinates two surface Fe_A cations and forms hydrogen-bonds with the exposed S and O anions, but no charge transfers were identified. The co-adsorption of two CO₂ and one H₂O molecule at nearby surface sites is also an exothermic process, although 2390 meV less favourable than the interaction of the single species.

We have calculated three minimum energy pathways for the formation of H₂C₂O₄, where we have considered the formation of an oxalate, carboxylate or bicarbonate intermediate in the first elementary step. The energy profiles show that C₂O₄* and HC₂O₄* are particularly unstable intermediates, whose formation requires crossing saddle points of large energy barriers. However, the carboxylate CO₂H* intermediate is the most stable species when it is co-adsorbed with CO₂*, and to a lesser extent, with the bicarbonate CO₃H*. The saddle points with the largest energies were observed for the coupling of the C–C bond, whereas the transfer of protons are the elementary steps with the smallest activation barriers. Our reaction mechanisms suggest that the energy released during the co-adsorption of the reactants is not enough to cross all the saddle points or to reach the final state. We would therefore argue that the feasibility of the conversion of CO₂ into H₂C₂O₄ depends strongly on the source of hydrogen and the reactivity of the surface. The high energy saddle points and intermediate species would become accessible if the catalytic process over the partially oxidised Fe₃S₄(001) surface is carried out within a continuous electrochemical cell with an appropriate voltage.

Future work will be focused on calculating the catalytic formation of pyruvic acid, which is an important intermediate in several metabolic pathways, on the partially oxidised surfaces of Fe₃S₄, including not only the stable (001) surface, but also more reactive surfaces such as the (111) plane. The reaction profiles for the conversion of CO₂ into pyruvic acid will allow us to compare this process and interconnect it into the catalytic formation of acetic acid, formic acid and oxalic acid.

Conflicts of interest

The authors have no conflicts of interest to declare.

Acknowledgements

Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202/1 and EP/R029431/1), this work made use of the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk). This research was undertaken using the supercomputing facilities at Cardiff University operated by Advanced Research Computing at Cardiff (ARCCA) on behalf of the Cardiff Supercomputing Facility and the HPC Wales and Supercomputing Wales (SCW) projects. We acknowledge the support of the latter, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government. Further work was undertaken on ARC4, part of the High-Performance Computing facilities at the University of Leeds, UK. All data are provided in full in the Results and Discussion section of this paper.

References

1. B. W. Strobel, Influence of vegetation on low-molecular-weight carboxylic acids in soil solution—a review, Geoderma, 2001, 99, 169–198.
2. W. Riemenschneider and M. Tanifuji, Ullmann's Encyclopedia of Industrial Chemistry, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2011.
3. H. Sawada and T. Murakami, Kirk-Othmer Encyclopedia of Chemical Technology, Wiley, 2000.
4. F. Palmieri, A. Estoppey, G. L. House, A. Lohberger, S. Bindschedler, P. S. G. Chain and P. Junier, Advances in Applied Microbiology, Elsevier Inc., 1st edn, 2019, vol. 106, pp. 49–77.
5. D. Panias, M. Taxiarchou, I. Paspaliaris and A. Kontopoulos, Mechanisms of dissolution of iron oxides in aqueous oxalic acid solutions, Hydrometallurgy, 1996, 42, 257–265.
6. W. W. Frenier and F. B. Growcock, Mechanism of Iron Oxide Dissolution—A Review of Recent Literature, Corrosion, 1984, 40, 663–668.
7. E. Baumgartner, M. A. Blesa, H. Marinovich and A. J. G. Maroto, Heterogeneous electron transfer as a pathway in the dissolution of magnetite in oxalic acid solutions, Inorg. Chem., 1983, 22, 2224–2226.
8. S. O. Lee, T. Tran, B. H. Jung, S. J. Kim and M. J. Kim, Dissolution of iron oxide using oxalic acid, Hydrometallurgy, 2007, 87, 91–99.
9. B. Hong, G. Xue, X. Guo and L. Weng, Kinetic study of oxalic acid pretreatment of moso bamboo for textile fiber, Cellulose, 2013, 20, 645–653.
10. B. Hong, L. Chen, G. Xue, Q. Xie and F. Chen, Optimization of oxalic acid pretreatment of moso bamboo for textile fiber using response surface methodology, Cellulose, 2014, 21, 2157–2166.
11. D. Marković, J. Ašanin, T. Nunney, Ž. Radovanović, M. Radoičić, M. Mitrić, D. Mišić and M. Radetić, Broad Spectrum of Antimicrobial Activity of Cotton Fabric Modified with Oxalic Acid and CuO/Cu₂O Nanoparticles, Fibers Polym., 2019, 20, 2317–2325.
12. J. A. Krasowski and J. Marton, The Formation of Oxalic Acid During Bleaching of Kraft Pulp, J. Wood Chem. Technol., 1983, 3, 445–458.
13. A. Brolin, J. Gierer and Y. Zhang, On the selectivity of ozone delignification of softwood kraft pulps, Wood Sci. Technol., 1993, 27, 115–129.
14. M. B. Roncero, J. F. Colom and T. Vidal, Why oxalic acid protects cellulose during ozone treatments?, Carbohydr. Polym., 2003, 52, 411–422.
15. M. B. Roncero, J. F. Colom, T. Vidal and M. A. Queral, The Role of Oxalic Acid in the Ozone Bleaching Kinetics of an Xo-Kraft Pulp, J. Wood Chem. Technol., 2000, 20, 147–167.
16. G. Gardea-Hernández, R. Ibarra-Gómez, S. G. Flores-Gallardo, C. A. Hernández-Escobar, P. Pérez-Romo and E. A. Zaragoza-Contreras, Fast wood fiber esterification. I. Reaction with oxalic acid and cetyl alcohol, Carbohydr. Polym., 2008, 71, 1–8.
17. L. Chen, J. Y. Zhu, C. Baez, P. Kitin and T. Elder, Highly thermal-stable and functional cellulose nanocrystals and nanofibrils produced using fully recyclable organic acids, Green Chem., 2016, 18, 3835–3843.
18. D. Li, J. Henschen and M. Ek, Esterification and hydrolysis of cellulose using oxalic acid dihydrate in a solvent-free reaction suitable for preparation of surface-functionalised cellulose nanocrystals with high yield, Green Chem., 2017, 19, 5564–5567.
19. R. Bharti, C. B. Reddy and P. Das, Oxalic Acid as Sustainable CO Source for Pyrrolone-Fused Benzosuberenes Synthesis through Palladium Catalyzed Carbonylative Cyclization, ChemistrySelect, 2017, 2, 4626–4629.
20. N. R. Guha, V. Thakur, D. Bhattacherjee, R. Bharti and P. Das, Supported Rhodium Nanoparticle-Catalyzed Intermolecular Regioselective Carbonylative Cyclization of Terminal Alkynes using Oxalic Acid as Sustainable C₁ Source, Adv. Synth. Catal., 2016, 358, 3743–3747.
21. P. Langer, J. Wuckelt, M. Döring and R. Beckert, Reaction of Ambident Dianions with Oxalic Acid Dielectrophiles – Effect of the Heteroatoms of the Dinucleophile on the Regiochemistry of Cyclization, Eur. J. Org. Chem., 1998, 1467–1470.
22. E. Jorjani and M. Shahbazi, The production of rare earth elements group via tributyl phosphate extraction and precipitation stripping using oxalic acid, Arab. J. Chem., 2016, 9, S1532–S1539.
23. W. Zhang, A. Noble, B. Ji and Q. Li, Effects of contaminant metal ions on precipitation recovery of rare earth elements using oxalic acid, J. Rare Earths, 2022, 40, 482–490.
24. T. V. Hoogerstraete, B. Blanpain, T. Van Gerven and K. Binnemans, From NdFeB magnets towards the rare-earth oxides: a recycling process consuming only oxalic acid, RSC Adv., 2014, 4, 64099–64111.
25. R. Chi and Z. Xu, A solution chemistry approach to the study of rare earth element precipitation by oxalic acid, Metall. Mater. Trans. B, 1999, 30, 189–195.
26. M. Gürü, A. Y. Bilgesü and V. Pamuk, Production of oxalic acid from sugar beet molasses by formed nitrogen oxides, Bioresour. Technol., 2001, 77, 81–86.
27. S. D. Deshpande and S. N. Vyas, Oxidation of Sugar to Oxalic Acid and Absorption of Oxides of Nitrogen to Sodium Nitrite, Ind. Eng. Chem. Prod. Res. Dev., 1979, 18, 69–71.
28. K. E. Richardson and N. E. Tolbert, Oxidation of Glyoxylic Acid to Oxalic Acid by Glycolic Acid Oxidase, J. Biol. Chem., 1961, 236, 1280–1284.
29. S. K. Dube, P. Vasudevan and B. L. Khandelwal, Oxalic acid manufacture, J. Chem. Technol. Biotechnol., 2007, 32, 909–919.
30. T. Matsumoto, M. Sadakiyo, M. L. Ooi, S. Kitano, T. Yamamoto, S. Matsumura, K. Kato, T. Takeguchi and M. Yamauchi, CO₂-Free Power Generation on an Iron Group Nanoalloy Catalyst via Selective Oxidation of Ethylene Glycol to Oxalic Acid in Alkaline Media, Sci. Rep., 2015, 4, 5620.
31. T. Matsumoto, M. Sadakiyo, M. L. Ooi, T. Yamamoto, S. Matsumura, K. Kato, T. Takeguchi, N. Ozawa, M. Kubo and M. Yamauchi, Atomically mixed Fe-group nanoalloys: catalyst design for the selective electrooxidation of ethylene glycol to oxalic acid, Phys. Chem. Chem. Phys., 2015, 17, 11359–11366.
32. F. Goodridge and G. Presland, The electrolytic reduction of carbon dioxide and monoxide for the production of carboxylic acids, J. Appl. Electrochem., 1984, 14, 791–796.
33. P. A. Nakata, Advances in our understanding of calcium oxalate crystal formation and function in plants, Plant Sci., 2003, 164, 901–909.
34. M. A. Webb, Cell-Mediated Crystallization of Calcium Oxalate in Plants, Plant Cell, 1999, 11, 751–761.
35. J. C. Yang and F. A. Loewus, Metabolic Conversion of l-Ascorbic Acid to Oxalic Acid in Oxalate-accumulating Plants, Plant Physiol., 1975, 56, 283–285.
36. M. V. Dutton and C. S. Evans, Oxalate production by fungi: its role in pathogenicity and ecology in the soil environment, Can. J. Microbiol., 1996, 42, 881–895.
37. J. A. L. Kan, M. W. Shaw and R. T. Grant-Downton, Botrytis species: relentless necrotrophic thugs or endophytes gone rogue?, Mol. Plant Pathol., 2014, 15, 957–961.
38. J.-C. Arvieu, F. Leprince and C. Plassard, Release of oxalate and protons by ectomycorrhizal fungi in response to P-deficiency and calcium carbonate in nutrient solution, Ann. For. Sci., 2003, 60, 815–821.
39. R. Hamel, R. Levasseur and V. D. Appanna, Oxalic acid production and aluminum tolerance in Pseudomonas fluorescens, J. Inorg. Biochem., 1999, 76, 99–104.
40. P. A. Nakata, The oxalic acid biosynthetic activity of Burkholderia mallei is encoded by a single locus, Microbiol. Res., 2011, 166, 531–538.
41. P. A. Nakata and C. He, Oxalic acid biosynthesis is encoded by an operon in Burkholderia glumae, FEMS Microbiol. Lett., 2010, 304, 177–182.
42. D. Markovich and P. S. Aronson, Specificity and Regulation of Renal Sulfate Transporters, Annu. Rev. Physiol., 2007, 69, 361–375.
43. T. Wang, A. L. Egbert, T. Abbiati, P. S. Aronson and G. Giebisch, Mechanisms of stimulation of proximal tubule chloride transport by formate and oxalate, Am. J. Physiol. Physiol., 1996, 271, F446–F450.
44. S. Albrecht, H. Brandl and C. Schönfels, Human Oxalate—Really Just an End-Product of Metabolism?, Angew. Chem., Int. Ed. Engl., 1994, 33, 1780–1781.
45. M. B. Rudnick, J. A. van Veen and W. de Boer, Oxalic acid: a signal molecule for fungus-feeding bacteria of the genus C ollimonas?, Environ. Microbiol. Rep., 2015, 7, 709–714.
46. L.-M. Zeng, J. Zhang, Y.-C. Han, L. Yang, M. Wu, D.-H. Jiang, W. Chen and G.-Q. Li, Degradation of oxalic acid by the mycoparasite Coniothyrium minitans plays an important role in interacting with Sclerotinia sclerotiorum, Environ. Microbiol., 2014, 16, 2591–2610.
47. Y. Han, H.-J. Joosten, W. Niu, Z. Zhao, P. S. Mariano, M. McCalman, J. van Kan, P. J. Schaap and D. Dunaway-Mariano, Oxaloacetate Hydrolase, the C–C Bond Lyase of Oxalate Secreting Fungi, J. Biol. Chem., 2007, 282, 9581–9590.
48. A. Mattevi, G. Obmolova, E. Schulze, K. Kalk, A. Westphal, A. de Kok and W. Hol, Atomic structure of the cubic core of the pyruvate dehydrogenase multienzyme complex, Science, 1992, 255, 1544–1550.
49. M. R. Mäkelä, K. Hildén and T. K. Lundell, Oxalate decarboxylase: Biotechnological update and prevalence of the enzyme in filamentous fungi, Appl. Microbiol. Biotechnol., 2010, 87, 801–814.
50. C. Plassard and P. Fransson, Regulation of low-molecular weight organic acid production in fungi, Fungal Biol. Rev., 2009, 23, 30–39.
51. J. M. Schuller, J. A. Birrell, H. Tanaka, T. Konuma, H. Wulfhorst, N. Cox, S. K. Schuller, J. Thiemann, W. Lubitz, P. Sétif, T. Ikegami, B. D. Engel, G. Kurisu and M. M. Nowaczyk, Structural adaptations of photosynthetic complex I enable ferredoxin-dependent electron transfer, Science, 2019, 363, 257–260.
52. G. Kurisu, M. Kusunoki, E. Katoh, T. Yamazaki, K. Teshima, Y. Onda, Y. Kimata-Ariga and T. Hase, Structure of the electron transfer complex between ferredoxin and ferredoxin-NAPD+ reductase, Nat. Struct. Biol., 2001, 8, 117–121.
53. A. Dey, T. Glaser, M. M. J. Couture, L. D. Eltis, R. H. Holm, B. Hedman, K. O. Hodgson and E. I. Solomon, Ligand K-Edge X-ray Absorption Spectroscopy of [Fe₄S₄] 1+, 2+, 3+ Clusters: Changes in Bonding and Electronic Relaxation upon Redox, J. Am. Chem. Soc., 2004, 126, 8320–8328.
54. R. H. Holm, P. Kennepohl and E. I. Solomon, Structural and Functional Aspects of Metal Sites in Biology, Chem. Rev., 1996, 96, 2239–2314.
55. P. Y.-T. Chen, B. Li, C. L. Drennan and S. J. Elliott, A Reverse TCA Cycle 2-Oxoacid:Ferredoxin Oxidoreductase that Makes C–C Bonds from CO₂, Joule, 2019, 3, 595–611.
56. A. C. Brown and D. L. M. Suess, Controlling Substrate Binding to Fe₄S₄ Clusters through Remote Steric Effects, Inorg. Chem., 2019, 58, 5273–5280.
57. G. Li, B. Zhang, F. Yu, A. A. Novakova, M. S. Krivenkov, T. Y. Kiseleva, L. Chang, J. Rao, A. O. Polyakov, G. R. Blake, R. A. de Groot and T. T. M. Palstra, High-Purity Fe₃S₄ Greigite Microcrystals for Magnetic and Electrochemical Performance, Chem. Mater., 2014, 26, 5821–5829.
58. A. Roldan, N. Hollingsworth, A. Roffey, H.-U. Islam, J. B. M. Goodall, C. R. A. Catlow, J. A. Darr, W. Bras, G. Sankar, K. B. Holt, G. Hogarth and N. H. De Leeuw, Bio-inspired CO₂ conversion by iron sulfide catalysts under sustainable conditions, Chem. Commun., 2015, 51, 7501–7504.
59. D. Santos-Carballal, A. Roldan and N. H. de Leeuw, CO₂ reduction to acetic acid on the greigite Fe₃S₄{111} surface, Faraday Discuss., 2021, 229, 35–49.
60. A. Roldan and N. H. de Leeuw, Methanol formation from CO₂ catalyzed by Fe₃S₄{111}: formate versus hydrocarboxyl pathways, Faraday Discuss., 2016, 188, 161–180.
61. D. Santos-Carballal, A. Roldan and N. H. de Leeuw, Early Oxidation Processes on the Greigite Fe₃S₄(001) Surface by Water: A Density Functional Theory Study, J. Phys. Chem. C, 2016, 120, 8616–8629.
62. S. N. A. Zakaria, N. Hollingsworth, H. U. Islam, A. Roffey, D. Santos-Carballal, A. Roldan, W. Bras, G. Sankar, G. Hogarth, K. B. Holt and N. H. de Leeuw, Insight into the Nature of Iron Sulfide Surfaces During the Electrochemical Hydrogen Evolution and CO₂ Reduction Reactions, ACS Appl. Mater. Interfaces, 2018, 10, 32078–32085.
63. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
64. G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
65. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
66. G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15–50.
67. J. Sun, A. Ruzsinszky and J. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Phys. Rev. Lett., 2015, 115, 1–6.
68. C. E. Mitchell, D. Santos-Carballal, A. M. Beale, W. Jones, D. J. Morgan, M. Sankar and N. H. de Leeuw, The role of surface oxidation and Fe–Ni synergy in Fe–Ni–S catalysts for CO₂ hydrogenation, Faraday Discuss., 2021, 230, 30–51.
69. D. Mejía-Rodríguez and S. B. Trickey, Meta-GGA performance in solids at almost GGA cost, Phys. Rev. B, 2020, 102, 1–4.
70. J. H. Yang, D. A. Kitchaev and G. Ceder, Rationalizing accurate structure prediction in the meta-GGA SCAN functional, Phys. Rev. B, 2019, 100, 1–10.
71. Y. Zhang, J. Sun, J. P. Perdew and X. Wu, Comparative first-principles studies of prototypical ferroelectric materials by LDA, GGA, and SCAN meta-GGA, Phys. Rev. B, 2017, 96, 035143.
72. J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, M. L. Klein and J. P. Perdew, Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional, Nat. Chem., 2016, 8, 831–836.
73. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
74. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
75. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem., 2006, 27, 1787–1799.
76. S. Posada-Pérez, D. Santos-Carballal, U. Terranova, A. Roldan, F. Illas and N. H. de Leeuw, CO₂ Interaction with Violarite (FeNi₂S₄) Surfaces: A Dispersion-Corrected DFT Study, Phys. Chem. Chem. Phys., 2018, 20, 20439–20446.
77. O. Lupan, D. Santos-Carballal, N. Ababii, N. Magariu, S. Hansen, A. Vahl, L. Zimoch, M. Hoppe, T. Pauporté, V. Galstyan, V. Sontea, L. Chow, F. Faupel, R. Adelung, N. H. de Leeuw and E. Comini, TiO₂/Cu₂O/CuO Multi-Nanolayers as Sensors for H₂ and Volatile Organic Compounds: An Experimental and Theoretical Investigation, ACS Appl. Mater. Interfaces, 2021, 13, 32363–32380.
78. D. Santos-Carballal, A. Cadi-Essadek and N. H. de Leeuw, Catalytic Conversion of CO and H₂ into Hydrocarbons on the Cobalt Co(111) Surface: Implications for the Fischer–Tropsch Process, J. Phys. Chem. C, 2021, 125, 11891–11903.
79. O. Lupan, N. Ababii, D. Santos-Carballal, M. I. Terasa, N. Magariu, D. Zappa, E. Comini, T. Pauporté, L. Siebert, F. Faupel, A. Vahl, S. Hansen, N. H. de Leeuw and R. Adelung, Tailoring the selectivity of ultralow-power heterojunction gas sensors by noble metal nanoparticle functionalization, Nano Energy, 2021, 88, 106241.
80. X. Liang, X. Lin, G. Wei, L. Ma, H. He, D. Santos-Carballal, J. Zhu, R. Zhu and N. H. De, Leeuw, Competitive adsorption geometries for the arsenate As(V) and phosphate P(V) oxyanions on magnetite surfaces: Experiments and theory, Am. Mineral., 2021, 106, 374–388.
81. L. M. Botha, D. Santos-Carballal, U. Terranova, M. G. Quesne, M. J. Ungerer, C. G. C. E. van Sittert and N. H. de Leeuw, Mixing thermodynamics and electronic structure of the Pt_1−xNi_x (0 ≤ x ≤ 1) bimetallic alloy, RSC Adv., 2019, 9, 16948–16954.
82. B. Ramogayana, D. Santos-Carballal, P. A. Aparicio, M. G. Quesne, K. P. Maenetja, P. E. Ngoepe and N. H. de Leeuw, Ethylene carbonate adsorption on the major surfaces of lithium manganese oxide Li_1−xMn₂O₄ spinel (0.000 < x < 0.375): a DFT+U-D3 study, Phys. Chem. Chem. Phys., 2020, 22, 6763–6771.
83. L. Reguera, N. L. López, J. Rodríguez-Hernández, M. González, C. E. Hernandez-Tamargo, D. Santos-Carballal, N. H. de Leeuw and E. Reguera, Synthesis, Crystal Structures, and Properties of Zeolite-Like [T₃(H₃O)₂M(CN)₆]₂·uH₂O (T = Co, Zn; M = Ru, Os), Eur. J. Inorg. Chem., 2017, 2980–2989.
84. A. Vahl, O. Lupan, D. Santos-Carballal, V. Postica, S. Hansen, H. Cavers, N. Wolff, M.-I. Terasa, M. Hoppe, A. Cadi-Essadek, T. Dankwort, L. Kienle, N. H. de Leeuw, R. Adelung and F. Faupel, Surface functionalization of ZnO:Ag columnar thin films with AgAu and AgPt bimetallic alloy nanoparticles as an efficient pathway for highly sensitive gas discrimination and early hazard detection in batteries, J. Mater. Chem. A, 2020, 8, 16246–16264.
85. M. J. Ungerer, D. Santos-Carballal, A. Cadi-Essadek, C. G. C. E. van Sittert and N. H. de Leeuw, Interaction of H₂O with the Platinum Pt(001), (011), and (111) Surfaces: A Density Functional Theory Study with Long-Range Dispersion Corrections, J. Phys. Chem. C, 2019, 123, 27465–27476.
86. T. A. Arias, M. C. Payne and J. D. Joannopoulos, Ab initio molecular dynamics: Analytically continued energy functionals and insights into iterative solutions, Phys. Rev. Lett., 1992, 69, 1077–1080.
87. M. J. Gillan, Calculation of the vacancy formation energy in aluminium, J. Phys.: Condens. Matter, 1989, 1, 689–711.
88. I. Štich, R. Car, M. Parrinello and S. Baroni, Conjugate gradient minimization of the energy functional: a new method for electronic structure calculation, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 4997–5004.
89. D. Sheppard, R. Terrell and G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys., 2008, 128, 134106.
90. M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand., 1952, 49, 409–436.
91. V. I. Anisimov, J. Zaanen and O. K. Andersen, Band theory and Mott insulators: Hubbard U instead of Stoner I, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 943–954.
92. A. Rohrbach, J. Hafner and G. Kresse, Electronic correlation effects in transition-metal sulfides, J. Phys.: Condens. Matter, 2003, 15, 979–996.
93. O. Bengone, M. Alouani, P. Blöchl and J. Hugel, Implementation of the projector augmented-wave LDA+U method: application to the electronic structure of NiO, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 16392–16401.
94. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
95. M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA+U method, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 035105.
96. D. Friebel, M. W. Louie, M. Bajdich, K. E. Sanwald, Y. Cai, A. M. Wise, M.-J. Cheng, D. Sokaras, T.-C. Weng, R. Alonso-Mori, R. C. Davis, J. R. Bargar, J. K. Nørskov, A. Nilsson and A. T. Bell, Identification of Highly Active Fe Sites in (Ni,Fe)OOH for Electrocatalytic Water Splitting, J. Am. Chem. Soc., 2015, 137, 1305–1313.
97. S. Fabris, S. de Gironcoli, S. Baroni, G. Vicario and G. Balducci, Taming multiple valency with density functionals: A case study of defective ceria, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 041102.
98. G. W. Watson, E. T. Kelsey, N. H. de Leeuw, D. J. Harris and S. C. Parker, Atomistic simulation of dislocations, surfaces and interfaces in MgO, J. Chem. Soc., Faraday Trans., 1996, 92, 433–438.
99. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
100. N. D. Mermin, Thermal Properties of the Inhomogeneous Electron Gas, Phys. Rev., 1965, 137, A1441–A1443.
101. P. E. Blöchl, O. Jepsen and O. K. Andersen, Improved tetrahedron method for Brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 16223–16233.
102. G. Makov and M. C. Payne, Periodic boundary conditions in ab initio calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 4014–4022.
103. J. Neugebauer and M. Scheffler, Adsorbate-substrate and adsorbate-adsorbate interactions of Na and K adlayers on Al(111), Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 16067–16080.
104. G. Henkelman, A. Arnaldsson and H. Jónsson, A fast and robust algorithm for Bader decomposition of charge density, Comput. Mater. Sci., 2006, 36, 354–360.
105. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, Improved grid-based algorithm for Bader charge allocation, J. Comput. Chem., 2007, 28, 899–908.
106. W. Tang, E. Sanville and G. Henkelman, A grid-based Bader analysis algorithm without lattice bias, J. Phys.: Condens. Matter, 2009, 21, 084204.
107. J. Tersoff and D. R. Hamann, Theory of the scanning tunneling microscope, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 805–813.
108. D. E. P. Vanpoucke and G. Brocks, Formation of Pt-induced Ge atomic nanowires on Pt/Ge(001): a density functional theory study, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 241308.
109. G. Henkelman, B. P. Uberuaga and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys., 2000, 113, 9901–9904.
110. G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys., 2000, 113, 9978–9985.
111. D. Santos-Carballal, A. Roldan, N. Y. Dzade and N. H. de Leeuw, Reactivity of CO₂ on the surfaces of magnetite (Fe₃O₄), greigite (Fe₃S₄) and mackinawite (FeS), Philos. Trans. R. Soc., A, 2018, 376, 20170065.
112. J. Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Math. Comput., 1980, 35, 773.
113. M. W. J. Chase, NIST JANAF Thermochemical Tables, American Chemical Society and American Institute of Physics for the National Institute of Standards and Technology, Washington DC, 1998.
114. X.-G. Wang, W. Weiss, S. K. Shaikhutdinov, M. Ritter, M. Petersen, F. Wagner, R. Schlögl and M. Scheffler, The Hematite (α-Fe₂O₃) (0001) Surface: Evidence for Domains of Distinct Chemistry, Phys. Rev. Lett., 1998, 81, 1038–1041.
115. B. Farkaš, D. Santos-Carballal, A. Cadi-Essadek and N. H. de Leeuw, A DFT+U study of the oxidation of cobalt nanoparticles: implications for biomedical applications, Materialia, 2019, 7, 100381.
116. D. Santos-Carballal, A. Roldan, R. Grau-Crespo and N. H. de Leeuw, A DFT study of the structures, stabilities and redox behaviour of the major surfaces of magnetite Fe₃O₄, Phys. Chem. Chem. Phys., 2014, 16, 21082–21097.
117. R. D. Johnson III, NIST Computational Chemistry Comparison and Benchmark Database, https://cccbdb.nist.gov/.
118. A. J. Devey, R. Grau-Crespo and N. H. de Leeuw, Electronic and magnetic structure of Fe₃S₄: GGA+U investigation, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 195126.
119. L. Chang, B. D. Rainford, J. R. Stewart, C. Ritter, A. P. Roberts, Y. Tang and Q. Chen, Magnetic structure of greigite (Fe₃S₄) probed by neutron powder diffraction and polarized neutron diffraction, J. Geophys. Res., 2009, 114, B07101.
120. A. Roldan, D. Santos-Carballal and N. H. de Leeuw, A comparative DFT study of the mechanical and electronic properties of greigite Fe₃S₄ and magnetite Fe₃O₄, J. Chem. Phys., 2013, 138, 204712.
121. J. M. D. Coey, M. R. Spender and A. H. Morrish, The magnetic structure of the spinel Fe₃S₄, Solid State Commun., 1970, 8, 1605–1608.
122. M. R. Spender, J. M. D. Coey and A. H. Morrish, The Magnetic Properties and Mössbauer Spectra of Synthetic Samples of Fe₃S₄, Can. J. Phys., 1972, 50, 2313–2326.
123. V. Hoffmann, Greigite (Fe₃S₄): magnetic properties and first domain observations, Phys. Earth Planet. Inter., 1992, 70, 288–301.
124. L. Chang, A. P. Roberts, A. R. Muxworthy, Y. Tang, Q. Chen, C. J. Rowan, Q. Liu and P. Pruner, Magnetic characteristics of synthetic pseudo-single-domain and multi-domain greigite (Fe₃S₄), Geophys. Res. Lett., 2007, 34, L24304.
125. M. Uda, On the Synthesis of Greigite, Am. Mineral., 1965, 50, 1487–1489.
126. M. J. Dekkers and M. A. A. Schoonen, Magnetic properties of hydrothermally synthesized greigite (Fe₃S₄)-I. Rock magnetic parameters at room temperature, Geophys. J. Int., 1996, 126, 360–368.
127. Z. He, S.-H. Yu, X. Zhou, X. Li and J. Qu, Magnetic-Field-Induced Phase-Selective Synthesis of Ferrosulfide Microrods by a Hydrothermal Process: Microstructure Control and Magnetic Properties, Adv. Funct. Mater., 2006, 16, 1105–1111.
128. P. W. Tasker, The stability of ionic crystal surfaces, J. Phys. C: Solid State Phys., 1979, 12, 4977–4984.
129. E. A. Wood, Vocabulary of Surface Crystallography, J. Appl. Phys., 1964, 35, 1306–1312.
130. A. Kiejna, T. Ossowski and T. Pabisiak, Surface properties of the clean and Au/Pd covered Fe₃O₄(111): DFT and DFT+U study, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 125414.
131. U. Terranova, C. Mitchell, M. Sankar, D. Morgan and N. H. de Leeuw, Initial Oxygen Incorporation in the Prismatic Surfaces of Troilite FeS, J. Phys. Chem. C, 2018, 122, 12810–12818.
132. C. E. Mitchell, U. Terranova, A. M. Beale, W. Jones, D. J. Morgan, M. Sankar and N. H. de Leeuw, A surface oxidised Fe–S catalyst for the liquid phase hydrogenation of CO₂, Catal. Sci. Technol., 2021, 11, 779–784.
133. G. Herzberg, Electronic spectra and electronic structure of polyatomic molecules, Van Nostrand, New York, 1966.
134. A. R. Hoy and P. R. Bunker, A precise solution of the rotation bending Schrödinger equation for a triatomic molecule with application to the water molecule, J. Mol. Spectrosc., 1979, 74, 1–8.
135. K. P. Huber and G. Herzberg, Molecular Spectra and Molecular Structure, Springer US, Boston, MA, 1979, pp. 8–689.
136. T. Shimanouchi, Tables of molecular vibrational frequencies, consolidated volume I, 1972, vol. 39.
137. K. K. Irikura, R. D. Johnson and R. N. Kacker, Uncertainties in Scaling Factors for ab Initio Vibrational Frequencies, J. Phys. Chem. A, 2005, 109, 8430–8437.
138. J. A. Pople, A. P. Scott, M. W. Wong and L. Radom, Scaling Factors for Obtaining Fundamental Vibrational Frequencies and Zero-Point Energies from HF/6-31G* and MP2/6-31G* Harmonic Frequencies, Isr. J. Chem., 1993, 33, 345–350.
139. M. W. Wong, Vibrational frequency prediction using density functional theory, Chem. Phys. Lett., 1996, 256, 391–399.
140. A. P. Scott and L. Radom, Harmonic Vibrational Frequencies:  An Evaluation of Hartree–Fock, Møller−Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors, J. Phys. Chem., 1996, 100, 16502–16513.
141. M. D. Halls, J. Velkovski and H. B. Schlegel, Harmonic frequency scaling factors for Hartree-Fock, S-VWN, B-LYP, B3-LYP, B3-PW91 and MP2 with the Sadlej pVTZ electric property basis set, Theor. Chem. Acc., 2001, 105, 413–421.
142. A. K. Mishra, A. Roldan and N. H. de Leeuw, CuO Surfaces and CO₂ Activation: A Dispersion-Corrected DFT+U Study, J. Phys. Chem. C, 2016, 120, 2198–2214.
143. A. Cadi-Essadek, A. Roldan and N. H. de Leeuw, Density functional theory study of the interaction of H₂O, CO₂ and CO with the ZrO₂ (111), Ni/ZrO₂ (111), YSZ (111) and Ni/YSZ (111) surfaces, Surf. Sci., 2016, 653, 153–162.
144. A. K. Mishra, A. Roldan and N. H. de Leeuw, A density functional theory study of the adsorption behaviour of CO₂ on Cu₂O surfaces, J. Chem. Phys., 2016, 145, 044709.
145. M. G. Quesne, A. Roldan, N. H. de Leeuw and C. R. A. Catlow, Carbon dioxide and water co-adsorption on the low-index surfaces of TiC, VC, ZrC and NbC: a DFT study, Phys. Chem. Chem. Phys., 2019, 21, 10750–10760.
146. N. Y. Dzade, A. Roldan and N. H. de Leeuw, Activation and dissociation of CO₂ on the (001), (011), and (111) surfaces of mackinawite (FeS): a dispersion-corrected DFT study, J. Chem. Phys., 2015, 143, 094703.
147. D. H. Gibson, Carbon dioxide coordination chemistry: metal complexes and surface-bound species. What relationships?, Coord. Chem. Rev., 1999, 185–186, 335–355.
148. D. H. Gibson, The Organometallic Chemistry of Carbon Dioxide, Chem. Rev., 1996, 96, 2063–2096.
149. M. Aresta, A. Dibenedetto and E. Quaranta, Reaction Mechanisms in Carbon Dioxide Conversion, Springer Berlin Heidelberg, Berlin, Heidelberg, 2016, pp. 35–69.
150. M. G. Quesne, A. Roldan, N. H. de Leeuw and C. R. A. Catlow, Carbon dioxide and water co-adsorption on the low-index surfaces of TiC, VC, ZrC and NbC: a DFT study, Phys. Chem. Chem. Phys., 2019, 21, 10750–10760.
151. S. C. Parker, N. H. de Leeuw and S. E. Redfern, Atomistic simulation of oxide surfaces and their reactivity with water, Faraday Discuss., 1999, 114, 381–393.
152. N. H. de Leeuw and T. G. Cooper, Surface simulation studies of the hydration of white rust Fe(OH)₂, goethite α-FeO(OH) and hematite α-Fe₂O₃, Geochim. Cosmochim. Acta, 2007, 71, 1655–1673.
153. N. H. de Leeuw and S. C. Parker, Computer simulation of dissociative adsorption of water on CaO and MgO surfaces and the relation to dissolution, Res. Chem. Intermed., 1999, 25, 195–211.
154. N. H. de Leeuw, G. W. Watson and S. C. Parker, Atomistic Simulation of the Effect of Dissociative Adsorption of Water on the Surface Structure and Stability of Calcium and Magnesium Oxide, J. Phys. Chem., 1995, 99, 17219–17225.
155. S. C. Parker, P. M. Oliver, N. H. De Leeuw, J. O. Titiloye and G. W. Watson, Atomistic simulation of mineral surfaces: Studies of surface stability and growth, Phase Transitions, 1997, 61, 83–107.
156. N. H. de Leeuw, J. A. Purton, S. C. Parker, G. W. Watson and G. Kresse, Density functional theory calculations of adsorption of water at calcium oxide and calcium fluoride surfaces, Surf. Sci., 2000, 452, 9–19.
157. Z. Du and N. H. de Leeuw, Molecular dynamics simulations of hydration, dissolution and nucleation processes at the α-quartz (0001) surface in liquid water, Dalton Trans., 2006, 2623–2634.
158. N. H. de Leeuw, F. M. Higgins and S. C. Parker, Modeling the surface structure and stability of α-quartz, J. Phys. Chem. B, 1999, 103, 1270–1277.
159. N. H. de Leeuw, Z. Du, J. Li, S. Yip and T. Zhu, Computer Modeling Study of the Effect of Hydration on the Stability of a Silica Nanotube, Nano Lett., 2003, 3, 1347–1352.
160. N. Y. Dzade, A. Roldan and N. H. de Leeuw, DFT-D2 Study of the Adsorption and Dissociation of Water on Clean and Oxygen-Covered {001} and {011} Surfaces of Mackinawite (FeS), J. Phys. Chem. C, 2016, 120, 21441–21450.
161. A. Roldan and N. H. de Leeuw, Catalytic water dissociation by greigite Fe₃S₄ surfaces: density functional theory study, Philos. Trans. R. Soc., A, 2016, 472, 20160080.
162. N. Y. Dzade, A. Roldan and N. H. de Leeuw, DFT-D2 simulations of water adsorption and dissociation on the low-index surfaces of mackinawite (FeS), J. Chem. Phys., 2016, 144, 174704.
163. S. Haider, A. Roldan and N. H. De Leeuw, Catalytic dissociation of water on the (001), (011), and (111) surfaces of violarite, FeNi₂S₄: A DFT-D2 study, J. Phys. Chem. C, 2014, 118, 1958–1967.
164. A. Roldan and N. H. de Leeuw, A kinetic model of water adsorption, clustering and dissociation on the Fe₃S₄ {001} surface, Phys. Chem. Chem. Phys., 2017, 19, 12045–12055.
165. N. H. de Leeuw and T. G. Cooper, A computational study of the surface structure and reactivity of calcium fluoride, J. Mater. Chem., 2003, 13, 93–101.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp00333c

---