Methanol formation from CO₂ catalyzed by Fe₃S₄{111}: formate versus hydrocarboxyl pathways

Abstract

Carbon capture and utilisation is one of the most promising techniques to minimize the impact of the increasing amount of carbon dioxide in the atmosphere. Recently, the mineral greigite was shown to be capable of catalysing CO₂ conversion, leading to useful small organic molecules. Here, we have carried out a systematic study of the adsorption and selective reduction of CO₂ on the Fe₃S₄{111} surface. We have considered both formate and hydrocarboxyl key intermediates, leading to different reaction pathways via Eley–Rideal and Langmuir–Hinshelwood mechanisms, and we have built a kinetic model considering the wide range of intermediates in the reaction network. Our results show that the mechanism to produce formic acid takes place via formate intermediate mostly on Fe_A sites, while methanol is formed via hydrocarboxyl intermediates on Fe_B sites. From the kinetic model, we have derived a reaction constant comparison and determined the limiting step rates. The overall process takes place under very mild conditions, requiring only a small energy input that might come from a chemiosmotic potential.

---

Introduction

Our present dependence on fossil fuels means that, as our demand for energy inevitably increases, so do emissions of greenhouse gases, most notably carbon dioxide (CO₂). To avoid the obvious consequence of severe climate change,¹ the concentration of such greenhouse gases in the atmosphere must be stabilized.^2–5 Integrated carbon capture and subsequent sequestration is generally recognized as the most promising option to tackle CO₂ gases in the short to medium term. A long-term solution is its reduction and conversion to suitable compounds: CO₂ is a source of carbon and may be used as a feedstock for conversion to useful chemicals.⁶ Reductive reactions involving CO₂ are a grand challenge because they require energy to generate reduced forms of CO₂, such as formate or CO.⁷ In particular, production of CO is considered an important objective in the context of production of renewable carbon feedstock chemicals.⁸ The main challenge in CO₂ reduction is increasing the energy efficiency of the activation process of the molecule, which is hindered primarily by high CO₂ reduction over-potentials,⁹ which is also suggested to be the rate limiting step in the conversion process. In most cases, mild temperatures but high CO₂ pressures (>10 bar) and long reaction times are required.^10,11 Nevertheless, experiments on iron sulfide have shown that electroreduction under room conditions and small overpotentials is achievable.¹²

Reactions of interest for CO₂ conversion are classified into two types: reductive and non-reductive reactions, both requiring the activation of the C=O double bond, which is accompanied by charge transfer to the lowest unoccupied molecular orbital.¹³ Solid state catalysts may act as a redox center, transferring electrons from its valence band to the molecular orbital while decreasing the O–C–O bond angle.^9,14 In this process, two problems are encountered: (i) the large thermodynamic stability of the molecule¹⁵ and (ii) the appreciable kinetic energy barrier to activate CO₂.¹⁶ Several possible reaction pathways have been proposed to transform CO₂ into one-carbon organic molecules, such as CH₄ or CH₃OH. These involve the reverse water gas-shift (rev-WGS) reaction, yielding CO followed by the subsequent hydrogenation via formyl intermediate (HCO) to methanol¹⁷ or to a variety of long chain hydrocarbons via the Fischer–Tropsch process.¹⁸ This reaction mechanism on Cu showed that, under ambient conditions, the reaction takes place through carboxyl intermediate.¹⁹ Other possible pathways have shown formate intermediate (HCOO) with subsequent O–CO bond scission, yielding similar products.^11,12 Recent studies have claimed that CH₃OH is produced by the subsequent hydrogenation of HCOOH to H₃CO via C–O bond scission,²⁰ although formic acid might dissociate back to formate with a low energy barrier.¹⁷ It has also been proposed that the hydrogenation of CO₂ towards CH₃OH takes place through dihydroxycarbene (COHOH), hydroxymethylidyne (COH) and hydroxymethylene (HCOH) intermediates.²⁰

Based on such mechanistic information, we have analyzed, in detail, the CO₂ interaction and conversion on greigite (Fe₃S₄) surfaces. This magnetite (Fe₃O₄) isomorphic mineral is formed as an intermediate in the solid-state transformation of mackinawite into pyrite, playing a crucial role in the pyrite formation pathway.^21–26 Fe₃S₄ has been widely identified in marine soils and sedimentary rocks of up to a few million years old.^27,28 Greigite has also been associated as a catalyst in a number of key biochemical reactions associated with the “iron–sulfur world” hypothesis for the origin of life.^12,29–37 Inspired by the CO₂ conversion in this hypothesis and recent experimental results, we have explored the entire reaction network on the Fe₃S₄{111} surface, mapping out the complete energy profile, and implemented it into a comprehensive kinetic model under particular initial conditions. This model is based on unrestricted elementary reaction steps derived rigorously from density functional theory (DFT) calculations, which have been shown to be instrumental in unravelling multiple aspects of heterogeneous catalysis.^17,20,38,39

Computational details

Electronic structure calculations

Periodic plane-wave DFT calculations were carried out to study the CO₂ adsorption and its reactivity with adsorbed hydrogen on the greigite surface Fe₃S₄{111}. All calculations were performed using the VASP software.^40,41 Ion–electron interactions were represented by the projector-augmented wave (PAW) method⁴² and the generalized gradient approximation (GGA) with the Perdew–Wang 91 functional⁴³ within the spin interpolation formula of Vosko et al.⁴⁴ We also considered non-spherical contributions from the gradient corrections to the PAW spheres. All the calculations include the long-range dispersion correction approach by Grimme,⁴⁵ which is an improvement on pure DFT when considering large polarizable atoms.^46–51 We have used the global scaling factor parameter optimized for PBE (s₆ = 0.75). The Kohn–Sham valence states were expanded in a plane-waves basis set with a cut-off of 600 eV for the kinetic energy.⁵² This high value for the cut-off energy ensured that no Pulay stresses occurred within the cell during relaxations.

The initial magnetic moment was described by high spin distributions in both types of Fe, octahedral and tetrahedral, with antiparallel orientation as reported previously.^53,54 For an accurate treatment of the electron correlation in the localized d-Fe orbital, we have used the U approximation^55,56 (U_eff = 1 eV),^57–59 which improves the description of localized states in this type of system where standard LDA and GGA functionals fail.⁶⁰ The choice of the U parameter is rather empirical, a feature that also appears when using computationally expensive hybrid functionals since the amount of Fock exchange is system-dependent.^60–63 Calculations were described by a Monkhorst–Pack grid of 5 × 5 × 1 K-points, which ensures electronic and ionic convergence.⁶⁴ The geometry of all stationary points was found with the conjugate-gradient algorithm and was considered converged when the force on each ion dropped below 0.03 eV Å⁻¹. The energy threshold defining self-consistency of the electron density was set to 10⁻⁵ eV. In order to improve the convergence of the Brillouin-zone integrations, the partial occupancies were determined using the tetrahedron method with Blöch corrections smearing, with a set width for all calculations of 0.02 eV. All total energies have been extrapolated to K_BT = 0 eV.

Computational model

The greigite unit cell consists of eight Fe₃S₄ subunits with a cubic lattice parameter of ∼9.8 Å,^65,66 which is close to the calculated parameter resulting from the cell optimization (9.671 Å). The inverse thio-spinel arrangement is reflected by the formula AB₂S₄, where there are two possible locations for the Fe ions: the tetrahedral sites (A), filled by Fe³⁺ ions, and the octahedral sites (B), where both Fe³⁺ and Fe²⁺ ions co-exist.^{24,25,58,67–70} In agreement with experiments, the electronic structure shows a high-spin configuration for both Fe sites in antiparallel alignments, resulting in a ferromagnetic material.^{24,25,59,67,68,70–72} The orbital spin-splitting in the valence region results in localized outermost 3d-electrons and in ordered magnetism.^69,73,74 Good agreement with experimental evidence has been obtained by using the same computational details as listed above.^53,54

We prepared the Fe₃S₄{111} surface as a slab model by cutting the bulk structure with the METADISE code,⁷⁵ which considers periodicity in the plane direction and provides atomic layer stacking, resulting in non-dipolar reconstructions.⁷⁶ The slab contains 56 atoms (24 Fe and 32 S) per unit cell, where the atoms are positioned in four layers of two Fe₃S₄ units each, exposing an area of 81.0 Ų and a thickness of sufficient size to relax the two uppermost layers (four Fe₃S₄ units) until energy convergence, keeping the bulk structure frozen at the bottom. We added a vacuum width of 12 Å between periodic slabs, big enough to avoid the interaction between images. Isolated molecules were placed in the center of a 15 × 16 × 17 ų cell and optimized with the same criteria.

Every intermediate is linked to the next minimum along the reaction mechanism by a saddle point in the pathway of minimum energy across the potential surface. These saddle points are the reaction transition states (TS), which determine the kinetics of the process. We have determined the TS by means of either the dimer method^77,78 or Climbing Image Nudge Elastic Band (CI-NEB),^79,80 depending on how far the starting structure is from the products. For instance, CI-NEB was used when in the initial structure hydrogen is further than 2.5 Å from the hydrogenated centre. The identified TS were further confirmed by a vibrational frequency calculation, in which only one imaginary frequency is obtained, corresponding to the reaction coordinate. The dimer image is relaxed to both local minima, leading to the initial state and the final state.

All binding energies (E_B) are given with respect to the naked slab and the reactants in the gas phase, see eqn (1). In the discussion of the DFT results, we have referred to energy values where the zero point energy (ZPE) correction is neglected unless noted otherwise.

E_{Mol+Fe3S4{111}} corresponds to the total energy of the adsorbate(s) on the Fe₃S₄{111} slab, E_Fe3S4{111} is the energy of the Fe₃S₄{111} slab and E_Mol is the energy of the isolated molecule(s) in vacuum. The activation energy (E_A) of a certain step is the energy required to surmount the potential barrier characteristic of the transition state. We have defined this energy barrier as the difference between the initial state and transition state energy for the forward process. The reaction energy (E_R) of each step is calculated as the total energy difference between the final state (product(s)) and the initial state (reactant(s)). We have discussed the feasibility of particular steps and the overall reaction trends (ΔE) as a function of the reference energy, i.e. clean surface and gas phase molecules, H₂ and CO₂.

Kinetic model

We have presented an extensive kinetic model^81–83 for the reduction of CO₂ on Fe₃S₄{111} surfaces on the basis of 140 elementary steps, with the migration of the CO intermediate between the two different active sites, Fe_A and Fe_B. No assumptions regarding the mechanism or the rate-limiting step are made, all the processes are bi-directional, and all model parameters are rigorously derived from DFT calculations, which we have compared with experimental data collected under realistic conditions, see Fig. S1 in ESI.† In the frame of transition state theory,^84,85 this kinetic model may poorly describe the surface kinetics when local ordering is important and diffusion of adsorbents as H ad-atoms is slow. In these cases, computationally expensive approximations can be used, for instance kinetic Monte Carlo (kMC) is capable of accounting for the correct local structure.⁸⁶ Unfortunately, such an analysis remains completely impractical for reaction networks as complicated as the one described here. Lateral adsorbate–adsorbate interactions can also influence the stability of intermediates.⁸⁷ Although we have included in the reaction network the effect of H ad-atoms near the catalytic center and the diffusion of CO₂ and CO molecules among Fe-sites, we have nevertheless considered the total surface H coverage (θ_H) to be low under typical conditions, and we have also included the H₂ evolution reaction as a function of θ_H. The maximum surface coverage was restricted to 1 ML, and multilayer adsorption was not considered. We have considered three different adsorption sites, a surface sulfur for the H ad-atom adsorption (eight per unit cell) and the metal centers, Fe_A or Fe_B (one site per cell).

Results and discussion

In order to develop a comprehensive model, firstly we have systematically studied the adsorption–desorption processes for the different reactants, followed by their co-adsorption and reaction leading to an extensive reaction network. Secondly, we have implemented these results in a micro-kinetics model to obtain a description of the process dynamics. Given the large amount of data, we have not described all the results in detail; the interested reader can find more relevant data in the ESI.† Finally, we have discussed the overall reaction trends and their implications in the CO₂ conversion process.

Density functional theory results

Greigite (Fe₃S₄) is an inverse spinel-structured material whose particles in hydrothermal synthesis expose {001} and {111} surfaces.⁸⁸ We have focused our attention on the Fe₃S₄{111} surface as this surface has shown conversion of carbonate species.³⁷ The Fe₃S₄{111} surface contains eight S atoms and three Fe atoms at the top layer of a unit cell, two from a tetrahedral bulk position (Fe_A) and one from an octahedral bulk position (Fe_B).

CO₂ adsorption and activation. Carbon dioxide interacts with the Fe centers as well as with a surface depleted in S through the oxygen lone pair electrons. The most stable configurations led to practically undisturbed molecules, see Table 1. The weak interactions with the surface, and the lack of charge transfer and molecule distortion, indicate physical adsorption on the naked surface, see Fig. 1. The activation of CO₂ from the physisorbed state is an endothermic process, but it is more likely to take place on Fe_B than on Fe_A, by 0.2 eV, perhaps due to the different ability of Fe_A and Fe_B to exchange electron density. The atomic charge analysis indicates that the transferred charge by Fe_A is twice that of Fe_B. We have noticed that, although most of the charge is transferred directly through the Fe site, its spin density remains unchanged. Fig. S3† represents the electron density depletion on the sulfur and appearance on both S–C and O lone pairs; rearrangement around the Fe atoms displays changes in the atomic orbital hybridization. Thus, the transferred charge is originally from the sulfur atoms near the metal centers. We have also considered the presence of atomic hydrogens on the CO₂ activation process. Accordingly to the adsorption energies, the presence of H ad-atoms on the surface stabilizes the activated species (CO^#₂) on both metal centers owing to the hydrogen providing labile electrons on the sulfur, which transfers them to the CO₂ antibonding orbital. CO₂ activation produces notable changes in the molecular structure: the O–C distance elongates by 0.07 Å, the O–C–O angle bends (∼130–160°), and the ν(C–O) shifts by 300 cm⁻¹ on Fe_B and up to 660 cm⁻¹ on Fe_A.

Table 1 Summary of binding energies (E_B) and optimized geometry parameters for the different adsorbents as well as for the activated CO₂ on Fe₃S₄{111} (# denotes the activated state). The sorption sites are described as Fe_B, Fe_A and S vacancy (v_S). The distances (d), angles (∠), and asymmetric stretching frequencies (ν_C–O) are given. The transferred charge (q) refers to the number of electrons transferred from the surface. Adsorbed structures are represented in the ESI, Fig. S3.1†

System	Site	E B (eV)	d O–C (Å)	∠O–C–O (°)	ν C–O (cm^−1)	q (e^−)
CO2	FeA	−0.02	1.176	178.3	2359	0.0
CO2	FeB	−0.19	1.175	179.4	2363	0.0
CO2	vS	−0.20	1.175	179.4	2366	0.0
CO^#2	FeA	0.40	1.233	138.1	1879	0.4
CO^#2	FeB	0.19	1.207	157.4	2103	0.2
CO^#2 + 2H	FeA	−0.33	1.250	131.9	1700	0.8
CO^#2 + 2H	FeB	−0.31	1.213	153.9	2050	0.8

---

Fig. 1 Schematic representation of non-activated CO₂ on Fe_A (left) and on Fe_B (right). The light grey and yellow balls represent Fe and S, respectively, and the red and dark grey balls represent O and C, respectively.

The positive E_B values of the activated CO₂ with respect to the energy reference indicate the preference of the molecule to desorb without modification before being activated. However, the presence of co-adsorbed hydrogens stabilizes CO^#₂ and its transition state, as shown by the reaction energies in Table 2. The lack of H ad-atoms on the surface makes the activation of CO₂ thermodynamically unfavourable and high pressures are required to increase the CO₂ activity. On the other hand, the presence of electron-donating species, i.e. H ad-atoms, enhances the CO₂ activation by decreasing the energy barrier, especially on electron deficient centers such as Fe_B.

Table 2 Elementary reaction energies for CO^*₂ activation on Fe₃S₄{111} (* denotes the adsorbed state and # denotes the activated state). E_A is the activation energy, E_R is the product reaction energy and ΔE indicates the product energy with respect to isolated molecules and the Fe₃S₄{111} slab. The structures are represented in the ESI, Fig. S3.1†

Label	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
2a	CO^*2 → CO^#2	FeA	0.43	0.39	0.36
2b	CO^*2 → CO^#2	FeB	0.41	0.38	0.20
2c	CO^*2 + 2H* → CO^#2 + 2H*	FeA	0.62	0.06	−0.33
2d	CO^*2 + 2H* → CO^#2 + 2H*	FeB	0.36	0.32	−0.30

---

Formate and hydrocarboxyl intermediates. From CO₂ adsorption and activation in the presence of H ad-atoms, we have systematically studied its hydrogenation on both Fe_A and Fe_B, leading to formate and hydrocarboxyl intermediates via Langmuir–Hinshelwood (LH) reaction mechanism, where both reactants lie on the surface. The hydrogenation with the lowest activation energy is on the oxygen (3c from Table 3), leading to COOH, which binds to S and Fe_A by C and O, respectively, see Fig. 2. The product is practically in thermodynamic equilibrium with the reactants by an energy difference of only 0.02 eV. The most favourable process on Fe_B is reaction 3d, which proceeds through the hydrogenation of the oxygen coordinated to Fe, at −0.16 eV below co-adsorbed CO^#₂ and two H ad-atoms.

Table 3 Elementary reaction energies for the CO^#₂ hydrogenation via Langmuir–Hinshelwood (LH) and Eley–Rideal (ER) mechanisms on the Fe₃S₄{111} surface. The structures are represented in the ESI, Fig. S3.2†

Label	Mechanism	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
3a	LH	CO^#2 + 2H* → HCOO* + H*	FeA	1.21	−0.28	−0.60
3b		CO^#2 + 2H* → HCOO* + H*	FeB	1.15	−1.25	−1.56
3c		CO^#2 + 2H* → COOH* + H*	FeA	0.10	0.02	−0.31
3d		CO^#2 + 2H* → COOH* + H*	FeB	0.71	−0.16	−0.46
3e	ER	CO2 + 2H* → HCOO* + H*	FeA	0.94	0.03	−0.60
3f		CO2 + 2H* → HCOO* + H*	FeB	1.86	−0.92	−1.56
3g		CO2 + 2H* → COOH* + H*	FeA	0.84	−0.17	−0.81
3h		CO2 + 2H* → COOH* + H*	FeB	1.28	0.19	0.65

---

Fig. 2 Schematic representation of hydrocarboxyl (from 3c) on Fe_A (left) and formate (from 3b) on Fe_B (right). The light grey and yellow balls represent Fe and S, respectively, and the red and dark grey balls represent O and C, respectively.

We have also investigated the Eley–Rideal (ER) mechanism, describing the interaction between gas phase CO₂ molecules and pre-adsorbed H on Fe₃S₄{111}. We have brought the molecule from a vacuum to the Fe center, exposing independently carbon or oxygen to the H ad-atom, which produces formate or hydrocarboxyl intermediates on the surface. The reaction energies are summarized in Table 3 (3e–3h). HCOO is derived from reaction 3e on Fe_A, which is the most suitable pathway for the formation of formic acid after surmounting a barrier of 0.3 eV above the energy reference. The generation of HCOO on Fe_B has an energy barrier of over 1 eV. Optimized HCOO binds to both Fe_A and Fe_B with an Fe–O distance of ∼1.8 Å, see Fig. 2. The formate intermediate was also identified experimentally on metallic Cu upon exposure to CO₂ and H₂ gases (ref. 89 and references therein).

The barrier energies for COOH formation via the ER mechanism are 0.1 (3g) and 0.6 eV (3h) lower than those for 3e and 3f on Fe_A and Fe_B, respectively. The formation of COOH as an intermediate is not surprising as it has been identified as the key reaction intermediate in the water gas shift reaction.¹⁹ Hydrocarboxyl has two distinguishable configurations depending on the direction of the hydroxyl group: trans-COOH and cis-COOH. The energy difference makes the trans-configuration the most stable by 0.4 eV.

Considering both mechanisms, i.e. LH and ER, the most favorable paths leading to HCOO intermediates are 3evia ER and 3bvia LH on Fe_A and Fe_B, respectively. In the transition state structure for the 3e process, almost linear CO₂ hovers close to the H ad-atom at a distance (C–H) of 1.526 Å and the molecule starts bending with an angle of 145°. The activation energy for the process 3b is 1.15 eV above that of the reactants and the C–H distance (1.692 Å) indicates an early transition state. The lowest activation energy for the formation of COOH is reaction 3c, although the energy required for the adsorption of the reactants and products on Fe_A is higher than that for the ER reaction 3g. The early saddle point for the hydrogenation of CO^#₂via3c places the H at 1.300 Å from the oxygen on Fe_A, while on Fe_B, the CO₂ in the 3d transition state lies at a distance of 1.374 Å from the H. Upon COOH formation, the HO–C bond becomes weaker and elongates by 0.1 Å compared with CO₂. We have further investigated the COOH dissociation to CO and OH, and found that it is thermodynamically and kinetically likely (E_A = 0.60 eV; E_R = −0.11 eV), proceeding to the formation of H₂O + CO.

The small difference between the energy barriers (2c and 3e < 0.05 eV) to HCOO and COOH intermediates indicates the low selectivity of the process on Fe_A, while a difference of 0.44 eV between the activation barriers of 3b and 3d indicates the selectivity for COOH intermediates on Fe_B. However, external conditions and the presence of dopants and promoters may influence the stability of intermediates, increasing the energy difference and therefore the selectivity.

Formic acid, formaldehyde and carbon monoxide. Formate intermediates remain adsorbed on the surface, allowing further changes such as C–O bond scission or additional hydrogenation. We have studied the breaking of the C–O bond but the process is thermodynamically and kinetically unfavourable (E_R > 1.5 eV; E_A > 1.8 eV). We have also studied the H migration process from formate to hydrocarboxyl but the process is kinetically unaffordable (E_A ∼ 3 eV). A consecutive hydrogenation of the formate intermediate leads to compounds such as formic acid (HCOOH) and formaldehyde (H₂CO), see Table 4.

Table 4 Elementary reaction energies for the HCOO hydrogenation on the Fe₃S₄{111} surface. Notice that we have differentiated O¹ and O² from the same CO₂ molecule. The structures are represented in the ESI, Fig. S3.3†

Label	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
4a	HCOO* + H* → HCOO^1H*	FeA	0.08	−0.45	−1.05
4b	HCOO* + H* → HCOO^1H*	FeB	0.56	0.28	−1.28
4c	HCOO* + H* → HCOO^2H*	FeA	1.28	0.19	−0.42
4d	HCOO* + H* → HCOO^2H*	FeB	2.12	0.71	−0.84
4e	HCOO* + H* → H2CO + O*	FeA	1.83	1.51	1.74
4f	HCOO* + H* → H2COO*	FeB	2.12	1.30	−0.26

---

We have explored the energy profile for HCOOH formation, see Fig. 3, and found that, from HCOO, the transition states leading to hydrogenation of the uncoordinated oxygen (4a,b) are more favourable than those leading to hydrogenation of the oxygen coordinated to Fe (4c,d). We have also noticed that the process is thermodynamically and kinetically more likely on Fe_B than on Fe_A. Between the two sites, pathway 4a is more feasible, leaving the HCOOH molecule adsorbed on Fe_A with a binding energy of −0.86 eV. The formation of HCOOH is considerably more likely than the formation of formaldehyde (4e and 4f), although this species has been identified on other catalysts as a potential intermediate in methanol production from CO₂.^17,90 Once HCOO is hydrogenated to form formaldehyde, its precursor (H₂COO) remains adsorbed on Fe_B, but on Fe_A, the C–O bond is completely broken and H₂CO is released from the surface, leaving an O ad-atom strongly adsorbed on Fe_A, thereby poisoning the active site. Formaldehyde is also a product of CO hydrogenation as described below.

Fig. 3 Schematic representation of HCOO¹H (from 4a) on Fe_A (left) and CO + H₂O (from 5f) on Fe_B (right). The light grey and yellow balls represent Fe and S, respectively, and the red and dark grey balls represent O and C, respectively.

Similarly to the formate intermediate, further hydrogenation of COOH leads to thermodynamically favorable HCOOH, but the activation energy (E_A > 1.2 eV) makes the reaction kinetically implausible (reactions 5a and 5b). As summarized in Table 5, dihydroxycarbene (HOCOH) is an alternative hydrogenation route as experimentally proposed in the methanol synthesis on oxide catalysts,⁹¹ but its formation is highly inaccessible from a kinetic point of view (5c and 5d). Nevertheless, we have considered HOCOH formation and dissociation to COH* + OH*, as suggested on Cu{111}.¹⁷ However, on Fe₃S₄{111}, the reaction pathway is inhibited not only by the activation energy forming HOCOH but also by the unlikely scission of HO–C, which has E_R > 0.8 eV and E_A > 1.4 eV on both Fe sites. Under certain conditions, HOCOH may be observable experimentally but it is not part of the main reaction pathway on greigite. The hydrogenation of the OH group on the COOH intermediate takes place mostly on Fe_B and leads to the formation of co-adsorbed H₂O and CO (reactions 5e and 5f). While the process is slightly endothermic on Fe_A with a large energy barrier above the reference of 1.6 eV, on Fe_B it is energetically favourable, see Fig. 3. The co-adsorbed H₂O stabilizes the system by 0.58 eV and, compared with HCOOH formation via reaction 3e, its energy barrier is higher by less than 0.05 eV.

Table 5 Elementary reaction energies for the COOH hydrogenation on the Fe₃S₄{111} surface. The structures are represented in the ESI, Fig. S3.4†

Label	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
5a	COOH* + H* → HCOOH*	FeA	1.21	−0.74	−1.05
5b	COOH* + H* → HCOOH*	FeB	1.38	−0.24	−0.70
5c	COOH* + H* → HOCOH*	FeA	1.39	−0.13	−0.44
5d	COOH* + H* → HOCOH*	FeB	1.67	0.07	−0.40
5e	COOH* + H* → CO* + H2O	FeA	1.84	0.07	−0.24
5f	COOH* + H* → CO* + H2O*	FeB	0.73	−0.49	−0.95

---

We have evaluated the movement of CO between both Fe_B and Fe_A from 5f and 5e, respectively, and also considered the presence of co-adsorbed H₂O (on Fe_B) or hydrogen ad-atoms, see Table 6. Before and after H₂O desorption from Fe_B, the migration of CO to Fe_A is endothermic and kinetically impractical (E_A > 1.5 eV), see Fig. 4. We have found that once water desorbs, the CO molecule shrinks the Fe_B–C bond by 0.08 Å and lies perpendicular to the surface at 1.754 Å. The intramolecular bond is 1.161 Å, which, together with a shift in its stretching frequency of 103 cm⁻¹ with respect to the isolated molecule (d_(C–O) = 1.142 Å, ν_(C–O) = 2129 cm⁻¹), indicates a weak C–O bond. The presence of co-adsorbed hydrogens decreases the energy barrier by up to 0.82 eV, making the process more likely, although _(A)CO is 0.16 eV less stable than _(B)CO. The effect of the hydrogens is similar to the CO₂ activation process, which is related to the availability of charge density. The stabilization of the CO adsorption on Fe_A by the H ad-atoms is observable in the Fe–C distance differences – the Fe_A–C bond is 0.18 Å larger than the Fe_B–C bond compared to 0.12 Å in the presence of hydrogens – and the asymmetric stretching frequencies vary by 23 and 11 cm⁻¹, respectively, without and with H ad-atoms. We have carried out the hydrogenation of CO on both sites, Fe_A and Fe_B, independently.

Table 6 Elementary reaction energies for the CO movements between Fe_B and Fe_A on the Fe₃S₄{111} surface. We have labelled each species regarding its position: (A) or (B). The structures are represented in the ESI, Fig. S3.5†

Label	Reaction	E A (eV)	E R (eV)	ΔE (eV)
6a	(B)CO* + (B)H2O* → (A)CO* + (B)H2O*	1.57	0.46	−0.49
6b	(B)CO* → (A)CO*	1.84	0.39	0.02
6c	(B)CO* + 2H* → (A)CO* + 2H*	0.82	0.16	−0.37

---

Fig. 4 Schematic representation of CO migration from Fe_B (left) to Fe_A (right) according to the process 6a. The light grey and yellow balls represent Fe and S, respectively, and the red and dark grey balls represent O and C, respectively.

Methanol. We have studied the hydrogenation of CO yielding CH₃OH on the Fe₃S₄{111} surface via two possible intermediates, CH_yOH and CH_xO, where the last intermediate may also produce formaldehyde. We have explored both energy profiles, summarized in Table 7 and 8, respectively. We have also connected both pathways with the hydrogenation of formaldehyde: CH₂O + H → CH₂OH.

Table 7 Elementary reaction energies for the CO hydrogenation on the Fe₃S₄{111} surface via the CH_yOH intermediate. The structures are represented in the ESI, Fig. S3.6†

Label	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
7a	CO* + 2H* → HCO* + H*	FeA	0.52	0.31	−0.06
7b	CO* + 2H* → HCO* + H*	FeB	0.37	0.17	−0.36
7c	HCO* + H* → HCOH*	FeA	0.55	0.26	0.20
7d	HCO* + H* → HCOH*	FeB	1.09	0.26	−0.11
7e	HCOH* + 2H* → H2COH* + H*	FeA	1.11	−0.81	−1.04
7f	HCOH* + 2H* → H2COH* + H*	FeB	0.86	−0.58	−1.08
7g	H2COH* + H* → H3COH*	FeA	1.17	−0.75	−1.79
7h	H2COH* + H* → H3COH*	FeB	0.85	−0.99	−2.07

---

Table 8 Elementary reaction energies for the HCO hydrogenation on the Fe₃S₄{111} surface via the CH_xO intermediate. The structures are represented in the ESI, Fig. S3.7†

Label	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
8a	HCO* + H* → H2CO*	FeA	0.94	0.75	0.69
8b	HCO* + H* → H2CO*	FeB	1.11	−0.57	−0.93
8c	H2CO* + 2H* → H3CO* + H*	FeA	0.34	−1.68	−1.23
8d	H2CO* + 2H* → H3CO* + H*	FeB	1.04	−0.58	−1.79
8e	H3CO* + H* → H3COH*	FeA	0.42	−0.98	−1.79
8f	H3CO* + H* → H3COH*	FeB	1.14	−0.28	−2.07

---

The carbon monoxide hydrogenation step is slightly unfavorable thermodynamically, E_R(Fe_B) < 0.2 eV, although the overall energy is lower than the reference energy (ΔE), indicating that it is a viable process. On Fe_A, the reaction 7a is endothermic by 0.3 eV whereas it is also kinetically controlled as the transition state lies above the energy reference. The transition states for the formation of HCOH (7c,d) and H₂CO (8a,b) also lie above the reference energy on both Fe_A and Fe_B, see Fig. 5. These processes, thus, require energy from external sources, such as a natural chemiosmotic or applied potential to surmount the energy barrier. The over-potential needed for the formation of HCOH on Fe_A (7c) is the smallest one (0.5 eV), leading to a downhill conversion towards CH₃OH. The HCOH intermediate has also been detected during CO reduction on transition metal surfaces.⁹² The presence of co-adsorbed hydrogens stabilises HCOH, whereas the production becomes thermodynamically unfavourable in pathway 8. The formation of HCOH (7c) has an energy barrier 0.4 eV lower than that for the production of H₂CO (8a), indicating selectivity of this process on Fe_A. On the Fe_B adsorption site, the E_A values are bigger than those for Fe_A because the reactants start from a more stable position (by 0.30 eV), although the transition states for 7d and 8b are 0.1 and 0.2 eV more favourable than their counterparts on Fe_A.

Fig. 5 Schematic representation of HCOH (from 7c) on Fe_A (left) and H₂CO (from 8a) on Fe_B (right). The light grey and yellow balls represent Fe and S, respectively, and the red and dark grey balls represent O and C, respectively.

Although CO reduction is likely to take place following process 7c, the hydrogenation of HCOH (7e) has a slightly higher energy barrier than the 8c process, by 0.1 eV, leading to competing processes. The reduction of H₂CO via the 9a process is 0.2 eV lower than that via7e, allowing the hydrogenation of the oxy group and lying −1.04 eV below the energy reference (Table 9). However, the 7f, 8d and 9b processes on Fe_B require less energy to surmount their transition states, leading to intermediates stabilised by up to 1.8 eV below the reference. Once methanol is formed through either the 7 or 8 pathway, it binds to Fe_A or Fe_B with an adsorption energy of −0.54 and −0.82 eV, respectively. Overall, the CH₃OH formation is more feasible on Fe_B where only 7d and 8b are limiting steps.

Table 9 Elementary reaction energies for the H₂CO hydrogenation on the Fe₃S₄{111} surface. The structures are represented in Fig. S3.8†

Label	Reaction	Site	E A (eV)	E R (eV)	ΔE (eV)
9a	H2CO* + 2H* → H2COH* + H*	FeA	0.18	−1.49	−1.04
9b	H2CO* + 2H* → H2COH* + H*	FeB	1.14	0.13	−1.08

---

Reaction energy profile. We have expressed, in energy profiles, all the steps described above for the CO₂ reduction at both Fe_A and Fe_B sites on the Fe₃S₄{111} surface, see Fig. 6 and 7, respectively. The activation of CO₂ is known to be an endothermic process, although the presence of atomic H on the surface stabilizes the activated molecule, making this a favourable process on Fe_B (2d). Upon CO₂ activation (2), the molecule may undergo hydrogenation on either C or O. Although reactions such as 3b may be driven by favourable thermodynamics, the activation energy is higher than the reference energy, indicating that the molecule would desorb before being hydrogenated. Reaction 3c is an exception as it has an activation energy of only 0.1 eV with a transition state 0.2 eV below the reference. Processes that do not require CO₂ pre-activation and proceed via the ER mechanism (3e–h) still require an energy input to surmount the activation energy and not all of them are thermodynamically driven. The difference in activation energies for the formation of HCOO and COOH is only 0.1 eV (3evs.3d) and therefore the selectivity for specific products, i.e. HCOOH and CH₃OH, is very low. Many of the following processes studied are energetically inhibited under mild conditions and the pathways are limited to reactions 4a on Fe_A, and 4b, 4d and 5f on Fe_B.

Fig. 6 Reaction energy profiles for the CO₂ hydrogenation process taking place on Fe_A. The blue line highlights the formate pathway. The labels corresponding to elementary steps are placed at the transition state. Processes such as water molecule release (−H₂O) and surface hydrogenation (+H₂) are also included. The horizontal brown line indicates the energy reference of the naked surface and isolated molecules.

Fig. 7 Reaction energy profiles for the CO₂ hydrogenation process taking place on Fe_B. The blue line highlights the formate pathway. The labels corresponding to elementary steps are placed at the transition state. Processes such as water molecule release (−H₂O) and surface hydrogenation (+H₂) are also included. The horizontal brown line indicates the energy reference of the naked surface and isolated molecules.

The formate adsorbed on Fe_A, the product of reaction 3e, is followed by a fast process producing HCOOH, which remains adsorbed, suggesting the possibility for secondary reactions, e.g. C–C formation or combining with NH₂ to give formamide.¹² This agrees with previous experiments, which observed large amounts of formate on metal surfaces upon exposure to CO₂ and H₂ gases.^89,93 Thus, the rate-limiting step for HCOOH formation is the first CO₂ hydrogenation with a TS energy 0.3 eV above the reference.

Recently identified as a key intermediate in the water gas shift reaction,^19,94 we have identified hydrocarboxyl as the first intermediate in the CO formation prior to methanol yield. The reverse water gas shift reaction takes place when oxygen from the CO₂ molecule reacts with H to produce CO and H₂O, i.e. processes 3d and 5f. This pathway on Fe_B led to carbon monoxide bonded to Fe_B with an adsorption energy of −1.1 eV, which is stabilized by the presence of H ad-atoms (by −0.17 eV). Although the formation of HCO is thermodynamically feasible, its hydrogenation has an energy barrier of at least 0.75 eV, processes (7 and 8), which are driven by very stable intermediates. The reduction of HCO is, hence, the rate-limiting step for CH₃OH formation from CO₂ on Fe₃S₄{111}. Once the limiting step is surmounted, the methoxy route (8) is the most favourable despite the fact that pathway 7 is also downhill energetically. Methoxy intermediates are very stable and remain on the surface for longer, which may lead to interactions with other molecules or intermediates, leading to e.g. longer chain hydrocarbons.¹²

Kinetic model results

We have applied a kinetic model, including the adsorption–desorption of the molecules and all the steps described above, in a consistent reaction network. We have characterized the desorption of CO₂, H₂O, H₂, HCOOH and CH₃OH individually and modelled the temperature programmed desorption (TPD) starting with full coverage and raising the temperature by 1 K s⁻¹. Molecules like CO₂, HCOOH and CH₃OH that have different adsorption sites, i.e. Fe_A and Fe_B, were equally distributed among them. For example, there are two configurations of HCOOH on Fe_A and two on Fe_B, and we set each configuration at an initial coverage of 0.25 ML. We have also considered secondary reactions like the exchange between CO₂ adsorption sites and the dissociation–association of the water molecule. The formation of H₂ is considered as a function of the H-coverage, for which we have investigated the associative desorption with different coverages, i.e. 0.25, 0.5, 0.75 and 1 ML, over the eight sulfur sites of an Fe₃S₄{111} cell unit. The overlap of these individual TPDs led to Fig. 8.

Fig. 8 Temperature programmed desorption of the molecules involved in the reduction mechanism of CO₂ on Fe₃S₄{111}. The temperature range considered is from 50 to 700 K, evaluated every 2 K each 2 s, leading to a heating rate of 1 K s⁻¹.

Physisorbed CO₂ desorbs from the surface with a minimal thermal energy, in agreement with the weak interaction with either Fe_B or Fe_A. Separately, H ad-atoms remain scattered on the surface until the temperature reaches 500 K, when they recombine and evolve as H₂ giving a sharp peak, which is linked to a fast process. Water molecules desorb from Fe_B relatively easily at 135 K with a narrow signal; note that it is a single molecule per unit cell. Formic acid is desorbed in three signals instead of the four expected from the different adsorbed conformations due to the overlapping of two signals around 378 K, from reactions 4a and 4d. Methanol TPD also presents multiple signals, at 119 K and at 185 K, related to the molecule adsorbed on Fe_A and Fe_B, respectively.

We have derived the reaction constants for the forward and the reverse reactions, tabulated in Table S1.† The relationship between these constants (K = K₊₁/K₋₁) provides important information with respect to a particular step. For instance, although reaction 3f has a high energy barrier compared with other CO₂ hydrogenation processes, the driving force of the HCOO formation makes K₊₁ bigger than K₋₁ and therefore the equilibrium constant (K) shifts towards formate. We have plotted the logarithm of K for every hydrogenation process in Fig. 9.

Fig. 9 Relationship between the equilibrium constants (log[K₊₁/K₋₁]) at a temperature of 300 K for the CO₂ reduction towards HCOOH (blue) and CH₃OH (red) on Fe_A (circles) and Fe_B (squares) sites on the Fe₃S₄{111} surface. The labels indicate the processes according to Fig. S3.† Formic acid and methanol desorption/adsorption constants are included on Fe_A and Fe_B, respectively. The solid symbols show the limiting step for each mechanism.

Although the activation of CO₂ on the surface is unlikely (reactions 2a–d), the presence of H ad-atoms enhances its molecular adsorption on either Fe_A or Fe_B as the top left values of Fig. 9 indicate. However, the Eley–Rideal mechanism for the formate formation on Fe_A is preferred and this is followed by a second hydrogenation (4a), leading to formic acid, whose desorption, however, is the slowest step, thereby reducing its production. On Fe_B, the formation of formic acid is capped by the second hydrogenation (4b) with the reverse process being more favourable, which means that formic acid released from Fe_A would adsorb on Fe_B and remain as formate, limiting the site availability. The processes 5c, 7a and 8a on Fe_A have similar consequences with respect to methanol formation. The methanol formation on Fe_B may take place via the highly favourable Eley–Rideal mechanism (3h) or through CO₂ activation and hydrocarboxyl formation (2b and 3d, respectively). The only slightly reverse step is the hydrogenation of CO (7b), which limits the process via the H₂CO intermediate, mechanism 8.

The mechanism showed HCOOH and CH₃OH as products of the CO₂ conversion with the reaction constants indicating that intermediate species like formate, carbon monoxide and CH_xO have a long enough residence time on the surface to allow them to react with other C-species, leading to e.g. acetic and pyruvic acid. Indeed these C₂ and C₃ products were detected experimentally after exposing Fe₃S₄ particles to a saturated aqueous solution of CO₂ under a mild reduction potential.¹²

Conclusions

We have carried out an extensive density functional theory study to map the CO₂ reduction mechanism network on the mineral surface Fe₃S₄{111}. Among the multiple reaction pathways explored, we have discerned the most favourable by using arguments based on energy profiles and reaction constants. Our results show distinctive behaviour for the different adsorption sites: Fe_A and Fe_B. On Fe_A, the CO₂ conversion occurs through a formate intermediate via Eley–Rideal mechanism, yielding HCOOH. On Fe_B, the formation of the hydrocarboxyl intermediate is favoured, which dissociates and yields CO on the surface plus a water molecule. Further hydrogenations of CO are limited by the formation of the HCO intermediate; despite being energetically favourable compared with the energy reference, the reaction constant towards HCO is smaller than that of the reverse processes. Following reductions are driven by the stability of their products. Methanol, in contrast to formic acid, desorbs from the surface easily, releasing an active site for further CO₂ conversion. Although the production of formic acid or methanol is determined by the adsorption site, CO₂ and intermediates are able to move from one site to another, as shown in step 6.

Thus, Fe₃S₄{111} catalyzes the CO₂ transformation by hydrogen, yielding mostly formic acid and methanol, where the rate-limiting step is the CO₂ activation via Langmuir–Hinshelwood mechanism, and indicating that the production of methanol is through the formation of the HCO intermediate. The selectivity towards particular products increases, e.g. by blocking selectively an Fe site, where capping Fe_A, for instance, leads to the reduction of CO₂ into methanol.

Acknowledgements

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

References

1. M. C. Urban, Science, 2015, 348, 571–573.
2. J. Rogelj, J. Nabel, C. Chen, W. Hare, K. Markmann, M. Meinshausen, M. Schaeffer, K. Macey and N. Hohne, Nature, 2010, 464, 1126–1128.
3. H. Arakawa, M. Aresta, J. N. Armor, M. A. Barteau, E. J. Beckman, A. T. Bell, J. E. Bercaw, C. Creutz, E. Dinjus, D. A. Dixon, K. Domen, D. L. DuBois, J. Eckert, E. Fujita, D. H. Gibson, W. A. Goddard, D. W. Goodman, J. Keller, G. J. Kubas, H. H. Kung, J. E. Lyons, L. E. Manzer, T. J. Marks, K. Morokuma, K. M. Nicholas, R. Periana, L. Que, J. Rostrup-Nielson, W. M. Sachtler, L. D. Schmidt, A. Sen, G. A. Somorjai, P. C. Stair, B. R. Stults and W. Tumas, Chem. Rev., 2001, 101, 953–996.
4. L. Plasseraud, ChemSusChem, 2010, 3, 631–632.
5. M. Meinshausen, N. Meinshausen, W. Hare, S. C. B. Raper, K. Frieler, R. Knutti, D. J. Frame and M. R. Allen, Nature, 2009, 458, 1158–1162.
6. M. Aresta and A. Dibenedetto, Dalton Trans., 2007, 2975–2992.
7. S. C. Roy, O. K. Varghese, M. Paulose and C. A. Grimes, ACS Nano, 2010, 4, 1259–1278.
8. G. Centi and S. Perathoner, Catal. Today, 2009, 148, 191–205.
9. H. J. Freund and M. W. Roberts, Surf. Sci. Rep., 1996, 25, 225–273.
10. M. R. Kember, A. Buchard and C. K. Williams, Chem. Commun., 2011, 47, 141–163.
11. N. S. Spinner, J. A. Vega and W. E. Mustain, Catal. Sci. Technol., 2012, 2, 19–28.
12. A. Roldan, N. Hollingsworth, A. Roffey, H. U. Islam, J. B. Goodall, C. R. Catlow, J. A. Darr, W. Bras, G. Sankar, K. B. Holt, G. Hogarth and N. H. de Leeuw, Chem. Commun., 2015, 51, 7501–7504.
13. H.-J. Freund and M. W. Roberts, Surf. Sci. Rep., 1996, 25, 225–273.
14. J. P. Pradiser and C. M. Pradier, Carbon Dioxide Chemistry: Environmental Issues, Woodhead Publishing, 1994.
15. M. Aresta and A. Dibenedetto, Dalton Trans., 2007, 2975–2992.
16. U. Jayarathne, P. Chandrasekaran, H. Jacobsen, J. T. Mague and J. P. Donahue, Dalton Trans., 2010, 39, 9662–9671.
17. Y. F. Zhao, Y. Yang, C. Mims, C. H. F. Peden, J. Li and D. H. Mei, J. Catal., 2011, 281, 199–211.
18. T. Riedel, M. Claeys, H. Schulz, G. Schaub, S. S. Nam, K. W. Jun, M. J. Choi, G. Kishan and K. W. Lee, Appl. Catal., A, 1999, 186, 201–213.
19. A. A. Gokhale, J. A. Dumesic and M. Mavrikakis, J. Am. Chem. Soc., 2008, 130, 1402–1414.
20. L. C. Grabow and M. Mavrikakis, ACS Catal., 2011, 1, 365–384.
21. S. Hunger and L. G. Benning, Geochem. Trans., 2007, 8, 1–20.
22. L. G. Benning, R. T. Wilkin and H. L. Barnes, Chem. Geol., 2000, 167, 25–51.
23. R. T. Wilkin and H. L. Barnes, Geochim. Cosmochim. Acta, 1996, 60, 4167–4179.
24. M. J. Dekkers, H. F. Passier and M. A. A. Schoonen, Geophys. J. Int., 2000, 141, 809–819.
25. M. J. Dekkers and M. A. A. Schoonen, Geophys. J. Int., 1996, 126, 360–368.
26. A. R. Lennie, S. A. T. Redfern, P. E. Champness, C. P. Stoddart, P. F. Schofield and D. J. Vaughan, Am. Mineral., 1997, 82, 302–309.
27. U. Frank, N. R. Nowaczyk and J. F. W. Negendank, Geophys. J. Int., 2007, 168, 904–920.
28. U. Frank, N. R. Nowaczyk and J. F. W. Negendank, Geophys. J. Int., 2007, 168, 921–934.
29. S. Mann, N. H. C. Sparks, R. B. Frankel, D. A. Bazylinski and H. W. Jannasch, Nature, 1990, 343, 258–261.
30. M. Posfai, P. R. Buseck, D. A. Bazylinski and R. B. Frankel, Am. Mineral., 1998, 83, 1469–1481.
31. M. Posfai, P. R. Buseck, D. A. Bazylinski and R. B. Frankel, Science, 1998, 280, 880–883.
32. W. Martin, J. Baross, D. Kelley and M. J. Russell, Nat. Rev. Microbiol., 2008, 6, 805–814.
33. M. J. Russell and A. J. Hall, Evolution of Early Earth's Atmosphere, Hydrosphere, and Biosphere: Constraints from Ore Deposits, 2006, vol. 198, pp. 1–32.
34. W. Martin and M. J. Russell, Philosophical Transactions of the Royal Society of London Series B-Biological Sciences, 2003, vol. 358, pp. 59–83.
35. C. Huber and G. Wachtershauser, Science, 1997, 276, 245–247.
36. G. Wachtershauser, Prog. Biophys. Mol. Biol., 1992, 58, 85–201.
37. G. Wachtershauser, Microbiol. Rev., 1988, 52, 452–484.
38. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757–10816.
39. J. K. Norskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37–46.
40. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
41. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15.
42. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
43. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671.
44. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
45. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
46. S. Irrera, A. Roldan, G. Portalone and N. H. de Leeuw, J. Phys. Chem. C, 2013, 117, 3949–3957.
47. N. Y. Dzade, A. Roldan and N. H. de Leeuw, J. Chem. Phys., 2013, 139, 124708.
48. S. S. Tafreshi, A. Roldan, N. Y. Dzade and N. H. de Leeuw, Surf. Sci., 2014, 622, 1–8.
49. S. Haider, A. Roldan and N. H. de Leeuw, J. Phys. Chem. C, 2014, 118, 1958–1967.
50. N. Dzade, A. Roldan and N. de Leeuw, Minerals, 2014, 4, 89–115.
51. F. Zhang, J. D. Gale, B. P. Uberuaga, C. R. Stanek and N. A. Marks, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 054112–054117.
52. N. Mermin, Phys. Rev., 1965, 137, 1441–1443.
53. A. Roldan, D. Santos-Carballal and N. H. de Leeuw, J. Chem. Phys., 2013, 138, 204712–204716.
54. D. Santos-Carballal, A. Roldan, R. Grau-Crespo and N. H. de Leeuw, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 195106.
55. V. I. Anisimov, M. A. Korotin, J. Zaanen and O. K. Andersen, Phys. Rev. Lett., 1992, 68, 345–348.
56. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
57. A. J. Devey, R. Grau-Crespo and N. H. de Leeuw, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 195126–195133.
58. L. Chang, B. D. Rainford, J. R. Stewart, C. Ritter, A. P. Roberts, Y. Tang and Q. W. Chen, J. Geophys. Res.: Solid Earth, 2009, 114, 1–10.
59. L. Chang, A. P. Roberts, Y. Tang, B. D. Rainford, A. R. Muxworthy and Q. Chen, J. Geophys. Res., 2008, 113, 1–16.
60. I. D. R. Moreira, F. Illas and R. L. Martin, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 155102–155116.
61. I. Ciofini, F. Illas and C. Adamo, J. Chem. Phys., 2004, 120, 3811–3816.
62. F. Illas and R. L. Martin, J. Chem. Phys., 1998, 108, 2519–2527.
63. D. Munoz, N. M. Harrison and F. Illas, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 085115–085124.
64. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
65. J. M. D. Coey, M. R. Spender and A. H. Morrish, Solid State Commun., 1970, 8, 1605.
66. M. R. Spender, J. M. D. Coey and A. H. Morrish, Can. J. Phys., 1972, 50, 2313–2326.
67. D. J. Vaughan and J. A. Tossell, Phys. Chem. Miner., 1983, 9, 253–262.
68. D. J. Vaughan and J. A. Tossell, Am. Mineral., 1981, 66, 1250–1253.
69. D. J. Vaughan and J. R. Craig, Am. Mineral., 1985, 70, 1036–1043.
70. K. K. Surerus, M. C. Kennedy, H. Beinert and E. Munck, Proc. Natl. Acad. Sci. U. S. A., 1989, 86, 9846–9850.
71. D. J. Vaughan and M. S. Ridout, J. Inorg. Nucl. Chem., 1971, 33, 741–746.
72. M. Braga, S. K. Lie, C. A. Taft and W. A. Lester, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 10837–10851.
73. D. J. Vaughan and J. A. Tossell, Am. Mineral., 1981, 66, 1250–1253.
74. D. Rickard and G. W. Luther, Chem. Rev., 2007, 107, 514–562.
75. G. W. Watson, E. T. Kelsey, N. H. deLeeuw, D. J. Harris and S. C. Parker, J. Chem. Soc., Faraday Trans., 1996, 92, 433–438.
76. P. W. Tasker, J. Phys. C: Solid State Phys., 1979, 12, 4977–4984.
77. G. Henkelman and H. Jonsson, J. Chem. Phys., 1999, 111, 7010–7022.
78. A. Heyden, A. T. Bell and F. J. Keil, J. Chem. Phys., 2005, 123, 224101–224115.
79. G. Henkelman and H. Jonsson, J. Chem. Phys., 2000, 113, 9978–9985.
80. G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000, 113, 9901–9904.
81. I. Chorkendorff and J. W. Niemantsverdriet, Concepts of Modern Catalysis and Kinetics, Wiley-VCH, Weinheim, 2005.
82. A. Roldan, G. Novell, J. M. Ricart and F. Illas, J. Phys. Chem. C, 2010, 114, 5101–5106.
83. P. Stoltze, Prog. Surf. Sci., 2000, 65, 65–150.
84. H. Eyring, J. Chem. Phys., 1935, 3, 107–115.
85. M. G. Evans and M. Polanyi, Trans. Faraday Soc., 1935, 31, 0875–0893.
86. H. S. Bengaard, J. K. Nørskov, J. Sehested, B. S. Clausen, L. P. Nielsen, A. M. Molenbroek and J. R. Rostrup-Nielsen, J. Catal., 2002, 209, 365–384.
87. L. C. Grabow, B. Hvolbaek and J. K. Norskov, Top. Catal., 2010, 53, 298–310.
88. N. Hollingsworth, A. Roffey, H.-U. Islam, M. Mercy, A. Roldan, W. Bras, M. Wolthers, C. R. A. Catlow, G. Sankar, G. Hogarth and N. H. de Leeuw, Chem. Mater., 2014, 26, 6281–6292.
89. Y. Yang, C. A. Mims, R. S. Disselkamp, J. H. Kwak, C. H. F. Peden and C. T. Campbell, J. Phys. Chem. C, 2010, 114, 17205–17211.
90. A. A. Peterson, F. Abild-Pedersen, F. Studt, J. Rossmeisl and J. K. Norskov, Energy Environ. Sci., 2010, 3, 1311–1315.
91. L. Z. Gao and C. T. Au, J. Catal., 2000, 189, 1–15.
92. P. Ferrin and M. Mavrikakis, J. Am. Chem. Soc., 2009, 131, 14381–14389.
93. E. Vesselli, M. Rizzi, L. de Rogatis, X. L. Ding, A. Baraldi, G. Comelli, L. Savio, L. Vattuone, M. Rocca, P. Fornasiero, A. Baldereschi and M. Peressi, J. Phys. Chem. Lett., 2010, 1, 402–406.
94. J. Yoshihara and C. T. Campbell, J. Catal., 1996, 161, 776–782.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5fd00186b

---