Does Fe 2 + in olivine-based interstellar grains play any role in the formation of H 2 ? Atomistic insights from DFT periodic simulations †

Understanding how H2 is formed in the interstellar medium (ISM) is of fundamental relevance for several reasons. Among the different molecular species detected, H2 is the most abundant molecule and is observed in every ISM environment. From an astrophysical perspective, H2 is an effective coolant of interstellar clouds during their gravitational collapse, thus facilitating the star formation, whereas from an astrochemical point of view, it plays a crucial role in most reaction networks to form other molecules, through neutral (H2) and ionized (H2 + and H3 ) forms. It is well established that the conversion of atomic H into molecular H2 requires the presence of interstellar dust grains (IDGs), sub-micron sized particles consisting of silicaceous and carbonaceous materials. This is because the gas-phase radiative association of two ground state H atoms (i.e., H + H H2 + hn) has a very low rate coefficient (E10 29 cm s ) due to spin forbidden transitions from the dissociative (Su) to the ground electronic state (Sg). Several experimental reports 4–14 and astrochemical modelling studies have shown that IDGs allow efficient H recombination to justify the large observed H2 abundances. An excellent review on this subject was recently published by Vidali. Silicates, and in particular olivines and pyroxenes (Mg2xFe(2x 2)SiO4 and MgxFe(x 1)SiO3 (x = 0–1), respectively), are important constituents of IDGs. Several theoretical studies have addressed H adsorption and H2 formation on silicate surfaces, using either cluster or periodic approaches to mimic the IDG surfaces, as well as using different DFT methods. Nevertheless, practically all of them focused on Mg-pure silicates, with only one addressing the H adsorption on Fe2SiO4. 22 Fe has a 3d electronic configuration that may lead to the formation of different low-lying electronic states and thus, it exhibits a more complex electronic structure than Mg, leading to the inference of a more interesting chemistry towards H adsorption and, accordingly, towards its reactivity. In the following, we report the results of H adsorption and subsequent recombination to form H2 on a Fe-containing olivine model surface, treated by ab initio periodic calculations performed using the CRYSTAL09 code, with the main goal of assessing whether Fe plays any significant and/or different role than Mg-pure olivines. Crystalline silicates in IDGs are Fe-poor in composition (10–15% of the total metal content), whereas amorphous silicates are richer in Fe (up to ratios of Fe/Mg E 1) but most of them either belong to the silicate matrix or are in the form of admixtures of metallic Fe. Accordingly, the amount of Fe placed at the IDG surfaces is relatively small. In this work, we adopted a periodic surface model of (010) Mg2SiO4 by replacing one surface Mg 2+ by Fe per unit cell (see Fig. 1A), giving rise a surface with the formula Mg1.875Fe0.125SiO4, hereafter referred to as Ol(010). The presence of only one surface Fe cation will avoid us working with complex electronic configurations and, as Ol(010) has one Fe and one Mg at the outermost positions, it will help us to gain a better understanding of the role of Fe and a proper comparison with Mg. In a previous work, we have described that Fe in its quintet state placed at the outermost positions of the (010) slab gives the most stable Fe-containing surfaces. Nonetheless, in the present work, the adsorption of one H atom on Fe leads to change in the oxidation state, resulting in either a quartet or a sextet spin configuration. Therefore, we performed a previous calibration study comparing the energy difference between the sextet and quartet states calculated with different DFT methods and at the CCSD(T) level on a cluster model consisting of the first Fe coordination sphere (the ESI† shows the cluster model adopted). This calibration study was performed using the Gaussian09 code. a Departament de Quı́mica, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. E-mail: albert.rimola@uab.cat b Dipartimento di Chimica and NIS – Nanostructured Interfaces and Surfaces – Interdepartment Centre, Università degli Studi di Torino, Via P. Giuria 7, 10125 Torino, Italy † Electronic supplementary information (ESI) available: Energetic values and structures, and extended computational details. See DOI: 10.1039/c6cc02313d Received 16th March 2016, Accepted 14th April 2016


Extended Computational Details
All periodic calculations have been performed with the ab initio CRYSTAL09 code. 1,2is code implements the Hartree−Fock and Kohn−Sham self-consistent field method based on localized Gaussian Type Orbitals (GTO) for periodic systems. 3The selfconsistent field (SCF) calculations and geometry optimizations were performed with the B3LYP-D2* functional, which includes an empirical a posteriori correction term proposed by Grimme 4 to account for dispersion forces (missed in the pure B3LYP 5,6 method), but whose initial parametrization (D2) was modified for extended systems (D2*), 7 to provide accurate results for the calculations of cohesive energies of molecular crystals and of adsorption processes within a periodic treatment; 8-10 and also with the BHLYP 11 functional, because it better describes the electronic structure of Fecontaining systems, as it is shown with the calibration study (see Table S1).Transition state (TS) search has been performed using the distinguished reaction coordinate (DRC) technique as implemented in CRYSTAL09, which has been proven to be robust and efficient enough for the proton jump of non-hydrated and hydrated acidic zeolites. 12The activated complex structures corresponding to the TS have been checked by ensuring that only one imaginary frequency resulted by the Hessian matrix diagonalization.All calculations involving one H atom have been run as open-shell systems based on the unrestricted formalism, whereas for dihydrogen adsorption the starting guess was open shell broken symmetry but collapsed to the closed shell system.Geometry optimizations have been performed in the P1 group symmetry (no symmetry), in order to ensure the maximum degrees of freedom during the optimization.Net charges and electron spin densities on the atoms were derived from the Mulliken population analysis.
The multi-electron wave function is described by linear combination of crystalline orbitals, which in turn are expanded in terms of GTO basis sets.Two different Gaussian basis sets have been adopted: (i) a B1 basis set described by the following all-electron contractions: (8s)-( 831sp (62111111s)-( 331111p)-(311d) for Fe; these basis functions were already used in previous works focused on the forsterite [13][14][15] and fayalite 16 bulk properties.For all calculations, a TZP basis set from Ahlrichs and coworkers 17 has been used for the H atoms.All the geometry optimizations have been carried out using the B1 basis sets and the energy is refined with single-point energy calculations at B2 onto the optimized B1 geometries (hereafter referred as B2//B1).In a previous work, 18 we showed that B2//B1 energy values are almost indistinguishable from those at B2//B2 level.
We set the shrinking factor of the reciprocal space net, defining the mesh of k points in the irreducible Brillouin zone, 19 to 5 and 20 for B1 and B2 calculations, respectively, requiring the diagonalization of the Hamiltonian matrix in 3 and 6 k points, respectively.
The accuracy of both Coulomb and exchange series was set to values of overlap integrals of 10 -6 and 10 -16 for both B1 and B2.A pruned (75, 974) grid has been used for the Gauss-Legendre and Lebedev quadrature schemes in the evaluation of functionals. 2,20The condition to achieve SCF convergence between two subsequent cycles was set to 10 -7 Hartree.Relaxations of both the internal atomic coordinates and the unit cell parameters for the bare surface, on the one hand, and relaxations of only the internal atomic coordinates keeping the lattice parameters fixed at the bare surface for the rest of the calculations, on the other hand, were carried out within the same run by means of analytical energy gradients 21 using a quasi-Newton algorithm, in which the quadratic step (Broyden−Fletcher−Goldfarb−Shanno Hessian updating scheme, BFGS) [22][23][24][25] is combined with a linear one as proposed by Schlegel. 26atom adsorption and recombination was only considered on the top surface of the slabs, since this greatly simplifies the localization of TS structures.We are conscious that this approach breaks the symmetry of the system with, nevertheless, a negligible effect on the energy profiles.The adsorption energies (ΔE) per mole of an H atom and per unit cell were computed as: where is the energy of the relaxed unitary cell containing the forsterite surface S in interaction with the H atom, is the energy of the relaxed unitary cell of the free forsterite surface, and is the energy of the free H atom.
CRYSTAL09 computes the zero-point energy (ZPE) corrections and thermodynamic quantities using the standard statistical thermodynamics formulae based on partition functions derived from the harmonic oscillator approximations, which are used to correct the adsorption energy values by temperature effects.The corresponding vibrational frequencies are calculated by obtaining the eigenvalues from diagonalization of the mass-weighted Hessian matrix at Γ point (point 0 in the first Brillouin zone, called the central zone).The mass-weighted Hessian matrix was obtained by numerical differentiation (central-difference formula) of the analytical first energy derivatives, calculated at geometries obtained by displacing, in turn, each of the 3N equilibrium nuclear coordinates by a small amount, 0.003 Å (N is the number of atoms). 27For the considered systems in this work, building up the full mass-weighted Hessian matrix would have been very expensive, so that only a portion of the dynamical matrix was computed by considering the displacements of a subset of atoms; i.e., the H atoms and the first and second-layer atoms of the surface.
Quantum tunnelling, a purely quantum mechanical effect, can play a significant role in the studied processes allowing chemical reactions to occur at significant rates at low temperatures, which classically would have negligible rates, 28 and the fact that they involve H atoms.The probability for a given system to tunnel through a reaction barrier depends primarily on the curvature of the barrier, which is controlled by the transition vibrational frequency and, to a lower degree, on the height of the barrier.The importance of tunnelling for a specific reaction can be estimated by employing the tunnelling crossover temperature , which can be calculated using the formula by (3) where ∆ ‡ is the free-energy barrier calculated at the temperature T.      Table S5.Adsorption energies (in kJ mol −1 ) calculated at the B3LYP-D2*/B2//B3LYP-D2*/B1 and BHLYP-D2*/B2//BHLYP/B1 level, respectively, for the global H adsorption processes on Ol(010) to form the 010-Fe1-Mg1, 010-Fe1-O1, 010-Fe1-O2, 010-O1-Mg1, 010-O2-Mg1, 010-O1-O1, 010-O1-O2 and 010-O2-O2 complexes at quintet spin state.

Reaction
Method -192.8 -5.1  -197.9   ) for H 2 formation from the most significant doubly-H adsorption complexes.All the structures are calculated in the quintet state, which are more stable than the heptet state (see Table S7).Adsorption energies (values in italics above the reactants) are referenced with respect to the Ol(010) + 2H zero-energy asymptote, whereas values of the energy profiles are referenced with respect to the corresponding reactants.Bond distances in Å.

Figure S1 .
Figure S1.(A) Cluster model adopted for the calibration study, in which one H atom is adsorbed on the Fe ion.(B) The cluster derives from the first coordination sphere around the Fe ion of the 010-Fe1 complex (inset structure of the first image) optimized at BHLYP level.The dangling bonds were saturated by H atoms at distances of 1.000 Å with the corresponding O atoms.The formula of the cluster is [Fe(OH) 2 (H 2 O)], in which two OH -groups are present to ensure the electroneutrality of the cluster as in the periodic surface.

Figure S3 .
Figure S3.B3LYP-D2*-optimized geometries of the different complexes for the H adsorption on the Fe-containing surface at the sextet state.Bond distances (in Å) and adsorption energies including zero-point energy corrections (in kJ mol -1 ) are also included.

Figure S5 .
Figure S5.B3LYP-D2*-energy profiles including zero-point energy corrections (in kJ mol -1 ) for H 2 formation from the most significant doubly-H adsorption complexes.All the structures are calculated in the quintet state, which are more stable than the heptet state (see TableS7).Adsorption energies (values in italics above the reactants) are referenced with respect to the Ol(010) + 2H zero-energy asymptote, whereas values of the energy profiles are referenced with respect to the corresponding reactants.Bond distances in Å.

Table S1 .
Energy differences between the sextet and quartet electronic states (E sextet -E quartet ) obtained at different DFT methods and at CCSD(T) (in kJ mol −1 ).Basis set employed in the DFT methods: 6-31G(d,p) Pople basis set for H and O, and a Wachters basis set with f functions for polarization for Fe.Basis set employed at CCSD(T): augcc-pVTZ for H and O atoms, and Roos Augmented Triple Zeta ANO for Fe.Calculations are based on single-point energy calculations onto the cluster model at the geometry of the 010-Fe1 complex optimized with the BHLYP method.

Table S2 .
Second ionization energies (in kcal mol −1 ) for Mg and Fe calculated at B3LYP, BHLYP, and CCSD(T) by employing Roos Augmented Triple Zeta ANO as a basis set.Experimental values from NIST are also included.

Table S4 .
Net charges and electronic spin densities on the H, Fe and Mg atoms, and the sum of the spin density values of the O atoms closest to H, computed at BHLYP-D2*/B2//BHLYP/B1 level for the different singly-H adsorptions on Ol(010) at sextet spin state.

Table S6 .
Net charges and electronic spin densities on the H, Fe and Mg atoms, and the sum of the spin density values of the O atoms closest to H, computed at BHLYP-D2*/B2//BHLYP/B1 level for the different doubly-H adsorptions on Ol(010) at quintet spin state (H 1st and H 2nd denotes the first and second H position according to the nomenclature of the complex, respectively).