This is the accepted version of the journal article: Pantaleone, Stefano; Corno, Marta; Rimola, Albert; [et al.]. «Ab-initio computa- tional study on Fe2NiP schreibersite : bulk and surface characterization». ACS Earth and Space Chemistry, Vol. 5, Issue 7 (July 2021), p. 1741-1751. DOI 10.1021/acsearthspacechem.1c00083 This version is available at https://ddd.uab.cat/record/256635 under the terms of the license Ab-initio computational study on Fe2NiP schreibersite: bulk and surface characterization Stefano Pantaleone,*1,2 Marta Corno,1 Albert Rimola,3 Nadia Balucani,2,4,5 and Piero Ugliengo*1 1 Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, I-10125, Torino, Italy 2 Dipartimento di Chimica, Biologia e Biotecnologie, Università degli Studi di Perugia, Via Elce di Sotto 8, I-06123 Perugia, Italy 3 Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalonia, Spain 4 Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy 5 Université Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), F-38000 Grenoble, France Keywords: Meteorites, phosphorus problem, DFT, prebiotic chemistry, surface modelling 1 Abstract Phosphorus is ubiquitous in planet Earth and plays a fundamental role in all living systems. Finding a reasonable prebiotic source of phosphorus is not trivial, as common sources where it is present nowadays are in the form of phosphate minerals, which are rather insoluble and non-reactive materials, and, accordingly, unavailable for being readily incorporated in living organisms. A possible source of phosphorus is from the exogenous meteoritic bombardment and, in particular, in iron/nickel phosphides. These materials, by simple interaction with water, produce oxygenated phosphorus compounds, which can easily react with organic molecules, thus forming C-O-P bonds. In the present work, periodic ab-initio simulations at PBE level (inclusive of dispersive interactions) have been carried out on metallic Fe2NiP-schreibersite, as a relative abundant component of metallic meteorites, in order to characterize structural, energetics and vibrational properties of both bulk and surfaces of this material. The aim is to study the relative stability among different surfaces, to also characterizinge both the nanocrystal morphology of the mineral and the reactivity towards water molecules. 2 Introduction Among the main biogenic elements (SPONCH), phosphorus represents only ⁓1% in mass of living organisms. Nevertheless, it is ubiquitous in the form of phosphate (PO43-): ATP (as well as ADP and AMP), nucleic acids, and phospholipids are paradigmatic compounds containing phosphates essential for life. Phosphorus abundance and its availability in prebiotic scenarios are still under debate, baffling the scientific community for decades. The major phosphorus reservoir on Earth is the mineral apatite (Ca5(PO4)3(F,OH,Cl)). However, it suffers from poor solubility and reactivity, and, accordingly, it can hardly be considered as a source of phosphorus ready to be incorporated in “living” phosphorus.1 In 1955 Gulick proposed, for the first time, as possible exogenous source of phosphorus the mineral schreibersite ((Fe,Ni)3P),2 undergoing corrosion by coming in contact with water. The proposal remained speculative up to the discovery of amino acids trapped in meteoritic materials from the first analysis on the Murchison meteorite.3–6 Because of this important finding, i.e. that biological matter could be delivered on Earth from exogenous sources, the role of phosphorus in this context was also reconsidered. Indeed, schreibersite is particularly interesting because Fe and Ni are among the first heavy elements to condense from the solar nebula7,8 (where schreibersite is one of the minor phases9).10 The commonly accepted hypothesis is that, during the heavy meteor bombardment in the Archean era (4.0-3.8 billion years ago), significant amounts of reduced phosphorus have been incorporated in the Earth’s crust.11–13 According to recent models, it was calculated that ca. 1-10% of the early Earth crust was composed of phosphide minerals.14 Very recently, Hess et al.15 proposed a new route to provide available phosphorus on the early Earth through lighting strikes as a major facilitator of prebiotic phosphorus reduction. Indeed, lighting strike hitting clay-rich soil of the early crust would form highly reduced glasses called fulgurites, 3 in which schreibersite is synthesized in situ by the high energy provided during the electric discharge. Clearly, this mechanism is independent of meteorite flux and provides prebiotic reactive phosphorus on Earth-like planets for a much longer time than that of the late meteoritic bombardment of the early Earth. The amount of schreibersite produced during the Hadean and early Archean era with that mechanism has been estimated to be 10-1000 kg/year of phosphide and 100-10000 kg/year of phosphite and hypophosphite. Therefore, the mineral schreibersite may really play a central role as a source of phosphorus due to its relatively high abundance and wide distribution in the continental crust. Schreibersite, irrespective of its origin, would then suffer corrosion by the weathering processes. Indeed, in 2005 Pasek et al. applied Gulick’s hypothesis by performing a corrosion experiment on iron phosphide (Fe3P) with water solutions of different compositions.16 The impressive results they obtained are not limited to the fact that many oxygenated phosphorus compounds were found after Fe3P corrosion, but, also, these species reacted with organic molecules. Moreover, the authors dealt with the abundance of the generated reactive P: they estimated that from a 60 ton meteorite, at least 1 ton of reduced phosphorus would be produced. Since then, the same research group carried out several experiments under different conditions of the substrate (by introducing a proper amount of nickel), temperature,17 pH,18 and reactants19 (also adding organic species).20 All these experiments clearly show the formation of oxygenated phosphorus compounds, in particular phosphites (HPO32-/H2PO3-), more reactive and soluble than phosphates, and also the phosphorylation of organic compounds, including nucleosides.18 An accurate study based on RAIRS (Reflection−Absorption InfraRed Spectroscopy) demonstrates that H2O molecules bind preferentially to P atoms on the schreibersite surface at low temperatures (120-140K),21 this mechanism being also confirmed by isotopically enriched (18O) water.22 4 Recently, it has been proven that schreibersite is a good catalyst to also activate the reaction among carbohydrates to give more complex sugars (formose reaction network), accompanied by well-known oxygenated phosphorus species (phosphites and, in minor quantity, phosphates).23 Additionally, schreibersite can also be of interest for industrial purposes. It is a transition metal phosphide (TMP) and these materials have been reported to be highly active electrocatalysts in water splitting reactions24–26 and to present high potential for electrochemical energy storage and conversion.27–29 To the best of our knowledge, very few computational studies have been performed on 30–32 schreibersite, and all of them regard the phase stability of different Fe3P polymorphs at different conditions of high temperature and pressure, in order to study the stability and possible phase changes to understand the presence of these materials in planetary cores. No computational work has been found on Fe2NiP-schreibersite, and, in particular, no theoretical studies on the surfaces of Fe3P and Fe2NiP have been carried out so far, which would be the first and foremost aspect to focus on if one is interested in corrosion reactions operated by water. Therefore, in the present work we have carried out first principles calculations by adopting density functional theory (DFT) on the bulk structure of Fe2NiP-schreibersite. Moreover, a detailed analysis of the low-Miller-index surfaces has been performed to predict elucidate the stability of the studied surfaces and predict the nanocrystal equilibrium morphology, which would be useful to assess and the potential reactivity of the exposed faces towards water, (which will be addressed specifically in a separate work). 5 Computational Details Bulk Periodic DFT calculations were carried out with the Vienna Ab-initio Simulation Package, (VASP) code,33–36 which uses projector-augmented wave (PAW) pseudopotentials37 to describe the ionic cores and a plane wave basis set for the valence electrons. This approach is particularly suitable, compared with localized basis functions, for describing delocalized electrons schreibersite due to its the intrinsic metallic character of schreibersite. Geometry optimizations and frequency calculations were performed with the gradient corrected PBE functional,38 which has already proven to provide good results on metallic systems.39,40 Moreover, several different schemes were adopted to add the contribution from dispersive (London) interactions (not accounted for by plain PBE) on the bulk structure to assess their performance in predicting the cell parameters with respect to the experimental ones. Specifically, the following dispersion corrections have been applied: • No dispersion correction (i.e. pure PBE functional) • D2 correction with Grimme’s original parameters;41 • D* correction, i.e. the Grimme’s D2 correction modified for solid systems, where van der Walls radii (R0) are rescaled by 1.05 (1.3 for H atoms);42 • D*0 correction, i.e. the D* correction where the C6 coefficients for the transition metals Fe and Ni are set to 0;43 • D*n correction, i.e. the D* correction where the C6 coefficients for the transition metals Fe and Ni are set to the values of the preceding noble gas (Argon);43 • D3 Grimme’s correction44 with Becke-Johnson (BJ) damping.45 6 The cutoff energy of plane waves (which control the accuracy of the calculations) was set to 500 and 1000 eV, in order to check the stability of the basis set when the cell volume of the bulk is free to relax. The self-consistent field (SCF) iterative procedure was converged to a tolerance in total energy of ΔE = 10-6 eV for geometry optimizations, while for frequency calculations the tolerance was increased to ΔE = 10-8 eV. As regards the optimization procedure, both atomic positions and cell parameters were free to relax. The tolerance on gradients was set to 0.01 eV/Å for each atom in each direction. The Monkhorst-Pack sampling of the Brillouin zone was used for the k-points mesh. Shrinking factors have been set on (8 8 8) and (16 16 16) for a total number of 76 and 552 irreducible k-points, respectively, in the Brillouin zone, ensuring a good accuracy in computing the electronic structure of a metallic system. The difference between the two k-grids in the total energy is 0.3 meV, and, accordingly, the (8 8 8) sampling was used to save computer resources when dealing with calculations on surfaces. Vibrational frequencies were computed on the bulk structure at Γ point, by numerical differentiation of the analytical first derivatives, using the central difference formula (i.e. two displacements of 0.02 Å for each atom in each (x,y,z) direction), in order to confirm that the optimized structure is a minimum (all real frequencies). To calculate Raman intensities, the Phonopy46 and Phono3py47,48 codes were used both for generating atomic displacements (only for symmetry irreducible atoms), and for processing VASP outputs. For a deeper analysis of the computed normal modes and a better assignment of low-frequency modes, whose atomic contributions are difficult to discriminate, the CRYSTAL17 code was used.49,50 In the specific, the CRYSTAL17 code has powerful tools to make isotopic substitutions, thus calculating a posteriori (i.e. at zero computational cost) the frequency shift due to the change of the atomic masses. This allows to separate the contribution of each atom to a certain normal mode on the basis of the 7 frequency shift when the isotopic substitution is applied. As this is a post-processing tool, (i.e. it does not require any additional ab-initio calculation as it relies on the stored force constant matrix), it does not depend on the chosen methodology, and therefore the isotopic shift calculated by CRYSTAL17 can be directly applied also to the VASP frequencies. Here we want to highlight that, despite the differences in the two codes (CRYSTAL17 uses localized gaussian functions as basis set, which are not effective in treating metals due to their intrinsic delocalized electrons), vibrational frequencies are in good agreement with that computed by VASP. We also checked, by visual inspection, the correspondence between VASP and CRYSTAL17 for each normal mode (see ESI). Surfaces After the careful check on the accuracy of the chosen computational parameters on the bulk structure, for the calculations on surfaces only two dispersion schemes were adopted: • D*0 correction, i.e. the D* correction where the C6 coefficients for transition metals Fe and Ni are set to 0; • D3 Grimme’s correction with Becke-Johnson (BJ) damping. The cutoff energy for plane waves was set to 500 eV (differences between 1000 and 500 eV on the bulk structure are negligible, vide supra). The self-consistent field (SCF) iterative procedure was converged to a tolerance in total energy of ΔE = 10-6 eV, as for the bulk. During the optimization procedure, only atomic positions were free to relax, keeping the cell parameters fixed to their bulk optimized values in order to simulate the rigidity imposed by the semi-infinite crystal underneath the surfaces. Shrinking factor for surfaces was set to (8 8 1), the last one being the number of k-points in the non-periodic direction (i.e. no sampling of the reciprocal space). As VASP relies on plane wave basis sets and, accordingly, surfaces are replicated also along the non- 8 periodic direction, the vacuum space among fictitious replicas was set to at least 20 Å in order to have no interactions among them. Therefore, the final c cell axis was set from 50 to 70 Å, depending on the surface thickness. Visualization and manipulation of the structures and figures rendering have been done with the MOLDRAW,51 VESTA,52 VMD,53 and POVRAY54 programs. Results Bulk structure The bulk structure of schreibersite is reported in Figure 1. The empirical formula is (Fe2,Ni)P. It belongs to the 𝐼4̅ (n° 82) space group (tetragonal family, with a=b≠c and α=β=γ=90°). It exhibits 16 irreducible atoms and 4 symmetry operators bringing a total of 32 atoms in the crystallographic unit cell.55 Table 1 shows the comparison between experimental and calculated cell parameters. As a general comment, all the computational parameters lead to a more compact structure with respect to the experimental one. The small changes between different cutoffs for the plane waves proves that 500 eV leads to converged results. In contrast, the system is sensitive to different dispersion schemes. It has already been demonstrated that for bulk transition metals pure PBE gives a good description of structure data, cohesive energies and bulk moduli.39,40 This is confirmed by our results: while pure PBE gives the best results, adoption of original Grimme’s dispersion parameters (D3 and D2) and the corrected one for solids (D*) leads to the largest deviations from experimental data. As D2 correction consists of few parameters, it is easy to change the dispersive contribution on selected atoms, and, in our specific case, the C6 coefficients on Fe and Ni have been reduced to those of the previous noble gas (Argon, D*n) or even switched off (D*0), this later case 9 presenting a very small deviation of the bulk parameters with respect to pure PBE. Details on the geometrical neighborhood of each irreducible atom are reported in the SI (see Table S1 and Figure S1). As our future purposes include the adsorption of water on the surfaces of this material, where in particular phosphorus plays a central role, instead of pure PBE for the surface analysis we used the D*0 correction where dispersion is left only on phosphorus (i.e., C6 coefficients were set to 0 on transition Fe and Ni metal atoms). From the point of view of the electronic structure, schreibersite is an open-shell conducting system (metallic behavior). Figure 2 clearly shows that there is no band gap among valence and conducting bands, i.e. there is an electron transfer as a result of the frontier band crossing. The metallic behavior of this system is confirmed by Bader charges56–58 reported in Table 2 which reveal an almost zero charge on Ni showing an atomic character, while on the two kind of Fe atoms Bader’s charges are more or less the same showing slightly positivized atoms (0.321e). Phosphorus charge keeps the cell neutrality by assuming a significant anionic character (-0.640e). These results are in line with prior studies of X-ray photoelectron spectroscopy, where the oxidation state of phosphorus is close to -1,17,59–62 as well as the higher “nobility” of the Ni atoms (which present a zero charge) compared with the Fe atoms.60,61 This confirmation is particularly important to disentangle the debate about the oxidation state of phosphorus, which in many cases is assumed to be -3, similar to that in compounds like AlP and Ca3P2.20 Fe and Ni in their atomic state have d6 and d8 configuration, respectively. From Table 2, it turns out that Ni atoms are basically in their atomic low spin electron configuration, i.e. all the electrons are coupled. On Fe atoms the high spin state would result in a quintet multiplicity (i.e. two paired and four unpaired electrons): however, the total spin moment is close to 2 for Fe1, which presents four paired and two unpaired electrons. On Fe2 instead, one of the two unpaired electrons is shared 10 with the other metals: the sum of the spin moment fractional part over all the metal atoms is close to 1 (0.277e (Fe1) + 0.386e (Fe2) + 0.245e (Ni), for the 500/NoD case), which represents the “missing” spin fraction on Fe2. As expected, phosphorus does not carry significant spin density. The same considerations on the accuracy of different dispersion schemes done for geometrical parameters are also valid for electronic data on Table 2, i.e. almost insensitivity towards the plane wave cutoff while being affected by different dispersion schemes. This is not surprising, as different dispersion schemes do not alter the electron density directly but shrink the structure ultimately affecting the final spin coupling. Finally, we computed the Raman spectrum only with the D*0 dispersion scheme, due to the high cost of the simulation. Experimental (adapted from Fig. 15e of ref. 62) and simulated (FWHM 50 cm-1) spectra are shown in Figure 3a. To account for the anharmonicity and the adopted functional systematic errors, we blue-shifted the computed spectrum by 40 cm-1 to align it with the experimental first peak at 212 cm-1 while conserving the internal band shifts. In ref.62 different samples of meteoritic schreibersite were extensively characterized by means of Raman spectroscopy, among other techniques. In the following, we will provide a full assignment of the computed peaks in comparison with those of synthetic schreibersite Fe2NiP.62 In Figure 3b the same computed Raman spectrum of Figure 3a is shown with a better resolution to highlight the fine structure of the spectrum. Our predicted spectrum does not show bands above 6500 cm-1 (the highest computed band is at 466 cm-1, last weak peak of Figure 3b), at variance with experimental spectra showing specific signals above 600 cm-1.62 The 466 cm-1 band is unambiguously attributed to the antisymmetric stretching of the lightest atoms present in the structure (Fe--P): therefore, the only way to justify higher frequency peaks in the experiments is to invoke the role of impurities, probably oxygen 11 atoms bringing new stretching Fe-O bands into play. Indeed, from the comparison with other works on schreibersite, but also on iron oxides, some weak signals between 500 and 700 cm- 1 can be attributed to inclusions of hematite, magnetite and goethite.61,63 Figure 4 shows the most relevant normal modes of the spectrum, highlighting the nuclear displacements to allow classifying each mode (GIF figures of each vibrational normal modes are available in the SI). The low-frequency modes (panels from Figure 4a to 4f, and first two peaks of Figure 3a) correspond to phonon modes involving heavy atoms (Fe and Ni). These modes are difficult to disentangle in terms of simple geometrical changes due to the coupled nature of these vibrations in the metallic system. Remarkably, these signals can be interpreted as the fingerprint of schreibersite. In Table 3, the various isotopic substitutions made on the structure are shown in traffic-light colors, from red (large shifts) to green (small shifts), which help to assign the contributions of each species to these highly coupled phonon modes. As one can see, P has basically no contributions to low-frequency modes, as its isotopic substitution (33P) does not affect those bands. In contrast, substitutions on Fe and Ni (i.e., 58Fe and 60Ni) produce frequency shifts up to 3 cm-1. Panels from Figure 4g to 4j (the shoulder at 300 cm-1 in the computed spectrum of Figure 3a) correspond to combining Fe1—Fe1 stretching and P-Fe2-P bending, while in panels of Figure 4k and 4l the P-Fe2-P bending modes are coupled with the Fe2-Ni stretching (310 cm-1 and 314 cm- 1 ). The last row of Figure 4 (from 4m to 4p) corresponds to Fe2-P stretching with different symmetry, with very small contributions from Fe1 and Ni atoms. 12 Surface models Starting from the crystal bulk structure of schreibersite (see Figure 1), periodic slab models with low Miller indices were built up. Specifically, the following five surface models were generated (see Figure 5): (100) = (010), (001), (101) = (011), (110) and (111). As it is well known, to be physically acceptable, a surface must be stoichiometric (i.e. electro neutral) and non-polar (i.e. no dipole moment across the non-periodic direction), which is achieved by adopting non-polar repeat units (RUs). The later point is achieved by modelling surfaces exhibiting symmetric top/bottom terminations. Therefore, fFrom all the possible facets cut out from the bulk structure, only the (100), (001), and (110) surfaces satisfy the abovementioned requirements (see Figure 5). The other surfaces, namely (011) and (111), in spite of being stoichiometric, exhibit a dipole moment across the non-periodic direction, and, accordingly, they will not be considered in the surface analysis due to their electronic instability as a function of the slab thickness (see Figure S4). For each physically acceptable surface, different terminations can be possible, and all of them have been investigated in order to find the most stable one. Different surface terminations are indicated using the following notation: (XXX)Y, where XXX are the Miller indices, and Y is the termination we are referring to. Specifically, the (001) facet presents 3 possible surface terminations ((001)1, (001)2, (001)3), while the (100) and (110) facets only present 2 terminations ((100)1, (100)2, (110)1, (110)2). From the different RUs represented in Figure 5, slab models of different thickness, ranging from 1 to 5 RUs, were created, in order to identify the minimum thickness which leads to converged properties. As a reference property to check for the convergence, we focused on the surface formation energy Es as a function of the surface thickness, adopting the following formula: (𝐸𝑠𝑙𝑎𝑏,𝑛 − 𝑛[𝐸𝑠𝑙𝑎𝑏,𝑛 − 𝐸𝑠𝑙𝑎𝑏,𝑛−1 ]) 𝐸𝑠 = 2𝐴 13 where Es is the surface formation energy, Eslab,n is the energy of the n-layer slab, Eslab,n-1 is the energy of the (n-1)-layer slab, and A is the surface area, which is accounted for twice (2) because of the top and bottom surfaces. In this equation the term [Eslab,n - Eslab,n-1] replaces the more commonly used Ebulk (the total energy of the bulk unit cell). This choice avoids numerical instabilities due to comparing energies of systems with different periodicity (see details in ESI).64– 66 Indeed, the classical formula referring to a bulk energy, is numerically unreliable when dealing with metals, for which the k-point sampling for system of different dimensionality affects the Es with large errors or instabilities. As a consequence, we only discuss results derived from the above equation. In Figure 6a, the converged surface formation energies (i.e. the values of the 5-layer slabs) of all the facets studied in this work is shown. The convergence of the surface energy versus the slab thickness for all the surfaces is shown in the ESI. It turns out that the most stable surface is the (110)2 with Es of 1.51 (D*0) and 2.28 (D3) J m-2. The majority of the other surfaces, namely (001)2, (001)3, (100)1, (100)2, exhibit Es values in the range around 2.01-2.06 (D*0) and 2.74-2.80 J m-2 (D3). The most unstable surfaces are the (001)1 and (110)1, with Es values of 2.40 (D*0) and 2.56 (D*0) J m-2, and 3.15 (D3) and 3.12 (D3) J m-2. Interestingly, there is a qualitative relationship between the Es values and the roughness of the surface, i.e. the rougher the surface the higher the Es. This is shown by the lateral views of the slabs (Figure 5, second and third columns) from which the (110)2 exhibits the flattest surface, while the (001)1 and (110)1 are more rumpled. In all cases, the Es rapidly converges to stable values with slab thickness of 3/4- layers. The D3 calculated Es values are all systematically larger by 0.7 J m-2 than those with the D*0 correction. This is due to the over-stabilization of the Grimme D3 dispersive energy of the bulk structure, which increases the cost of extracting the surface form the bulk compared to plain PBE and PBE-D*0. 14 The convergence of the Es is paired with that of other surface properties. For instance, we computed the Bader charges of the innermost atoms of the most stable (110)2 surface, which indeed converge to the bulk values (see Figure S5), i.e. the Fe atoms exhibit positive Bader charges, Ni atoms are almost neutral, while P atoms carry a large compensating negative charge. These values are expected to play a significant role in the reactivity of the different surfaces towards water. From the converged surface formation energies of all the studied facets, we built the corresponding crystallographic shape, as produced by the Wulff construction (Figures 6b).67 Interestingly, only the (110) and (001) faces appear represented in the crystal morphology, due to the high surface formation energy, and accordingly instability, of the (100) one. The large difference in the surface formation energies Es due to the adoption of different dispersion schemes (D*0 and D3) does not alter the shape of the Wulff polyhedron, as the relative surface formation energies Es for the different facets are method independent. 68 Conclusions This work deals with the quantum mechanical treatment of the mineral schreibersite Fe2NiP, a relevant material for the exogenous delivery of phosphorus to the early Earth through iron rich meteorite bombardment. Indeed, phosphorus may be released from schreibersite by contact with water in the form of phosphide/phosphates. A variety of biochemically relevant molecules have been shown to be formed when schreibersite is contacted with simple organic molecules, forming the fundamental the C-O-P bond. 15 Here, DFT simulations using the gradient corrected PBE functional, supplemented by the Grimme’s dispersion D2 and D3 corrections has been adopted to theoretically characterize some physicochemical properties of the bulk and different surfaces of this material. The blind application of the dispersion (London) contribution has been found to be detrimental for the prediction of the bulk schreibersite, as the cell parameters become underestimated compared to the plain PBE ones, which in turn are in good agreement with experiment. Following previous literature advices, we restore the good agreement with the experiment by setting to zero the C6 parameter of Fe and Ni elements (no contribution to dispersion), while keeping the original value for phosphorus. This ensures that dispersion will be properly included when adsorption of molecules is taken into account as an essential component for any non-covalent interactions. Band structure and project density of states revealed a metallic nature of schreibersite, while Bader charges reveal an almost zero charged Ni atom, a positively charged Fe ions (+0.321 e) compensated by a strong anionic character of the P atom (-0.640 e). The open shell character of schreibersite brings a spin moment value of almost nil on P atom, and strongly localized on Fe1 (2.271 e) and Fe2 (1.394 e) with a minor role of Ni atom (0.245 e). The simulated and experimental Raman spectra, which are in good agreement upon applying a proper scaling factor, were compared, and the principal bands were assigned with the purpose of helping the assignment of the experimental bands, as the samples are often plagued by impurities. Due to the metallic character of the material, the normal modes are fully coupled and extend to almost all atoms of the unit cell. Only modes involving P are well separated, falling in the highest frequency region of the spectrum. As regards the surfaces of schreibersite, low-Miller-indiced facets were studied with the aim of identifying those which are likely exposed to the external environment. Calculations show that 16 only the (110) and the (001) surfaces appear as extended faces when the Wulff theorem is applied to construct a representative structure of a likely nanoparticle. This study constitutes the first step towards further investigations in studying the interaction and reactivity of this material with molecules. Particularly appealing is to study reaction with water, as water reactions with the schreibersite surfaces are the first steps in getting P atoms in the soluble form of phosphites/phosphates that may be incorporated in biologically relevant molecules in prebiotic environments. Dedicated simulations for this purpose are ongoing. This is particularly useful for future purposes in studying the interaction of this material with whatever molecule of interest and in particular with water. Indeed, water reactions with the schreibersite surfaces are the first steps in getting P atoms in the soluble form of phsophites/phosphates able to be incorporated in biological relevant molecules in prebiotic environments. 17 a) Bulk structure – ab view b) Fe1 neighbors – bc view c) Fe2 neighbors – bc view d) Ni neighbors – bc view Figure 1. PBE-D*0 optimized geometry of the schreibersite bulk structure (crystallographic cell) in coordination polyhedral representation. Atom legend colors: P is in yellow, Fe in light grey, and Ni in dark grey. Unit cell is represented as thin black lines. 18 Spin down Spin up Figure 2. Calculated projected density of states (PDOS) and electronic band structure. Spin up and down are represented with the right and left portion of the image. 19 a) Experimental vs computed spectra b) High-resolution computed spectrum Figure 3. Top: Raman spectra of schreibersite: in red the experimental one62, adapted with permission from [C. Pirim, M.A. Pasek, D.A. Sokolov, A.N. Sidorov, R.D. Gann, T.M. Orlando, Investigation of schreibersite and intrinsic oxidation products from Sikhote-Alin, Seymchan, and Odessa meteorites and Fe3P and Fe2NiP synthetic surrogates, 2014, 140, 259-274]. Copyright 20 [2014] [Elsevier], in blue the computed one, with FWHM of 50 cm-1. Bottom: high-resolution computed Raman spectrum (FWHM 5 cm-1). 21 a) 211 cm-1 b) 216 cm-1 c) 238 cm-1 d) 252 cm-1 e) 260 cm-1 f) 273 cm-1 g) 275 cm-1 h) 279 cm-1 i) 297 cm-1 j) 307 cm-1 k) 310 cm-1 l) 314 cm-1 m) 354 cm-1 n) 421 cm-1 o) 426 cm-1 p) 441 cm-1 Figure 4. Atomic displacements (green arrows) related to the normal modes associated with the Raman spectrum of Figure 3b. Atomic legend colors: P is in yellow, Fe in light grey, and Ni in dark grey. 22 a) (001)1 – ab top b) (001)1 – ac side c) (001)1 – bc side d) (001)2 – ab top e) (001)2 – ac side f) (001)2 – bc side g) (001)3 – ab side h) (001)3 – ac side i) (001)3 – bc side j) (110)1 – ab top k) (110)1 – ac side l) (110)1 – bc side 23 m) (110)2 – ab top n) (110)2 – ac side o) (110)2 – bc side p) (100)1 – ab top q) (100)1 – ac side r) (100)1 – bc side s) (100)2 – ab top t) (100)2 – ac side u) (100)2 – bc side Figure 5. Physically sounded Stoichiometric and non-polar ((001), (110), (100) and unsounded ((011), (111), i.e. stochiometric but polar) facets cut out from the bulk structure. Atom legend colours: P is in yellow, Fe in light grey, and Ni in dark grey. 24 Figure 6. Left: Converged (5-layer slab) surface formation energy Es (J m-2) of all the studied surfaces. The surfaces which do not contribute to the Wulff shape are represented in grey-scale colors. Right: Wulff shape superimposed to a representative nanoparticle structure. Colors in the chart correspond to those of the Wulff shape. 25 Table 1. Experimental (EXP) vs calculated (with different basis sets and dispersion schemes) cell parameters (units of Å) of the schreibersite bulk structure, with the corresponding percentage deviation. cutoff/dispersion a c diff % a diff % c 500/NoD 8.984 4.377 -0.327 -1.062 1000/NoD 8.996 4.382 -0.185 -0.946 500/D2 8.878 4.332 -1.496 -2.073 500/D* 8.887 4.329 -1.402 -2.144 500/D*n 8.925 4.350 -0.977 -1.675 500/D*0 8.984 4.374 -0.320 -1.121 500/D3 8.882 4.331 -1.453 -2.109 1000/D3 8.892 4.335 -1.341 -2.022 EXP55 9.013 4.424 26 Table 2. Bader charge analysis and spin moment (in elementary charge units, e) of the irreducible atoms in the bulk primitive unit cell calculated with different basis sets and dispersion schemes. Net charge Atom 500/NoD 1000/NoD 500/D2 500/D* 500/D*n 500/D*0 500/D3 1000/D3 P -0.640 -0.635 -0.607 -0.609 -0.626 -0.641 -0.606 -0.603 Fe 0.321 0.320 0.308 0.309 0 .315 0.321 0.307 0.308 Ni -0.001 -0.006 -0.008 -0.010 -0.005 -0.002 -0.009 -0.012 Spin moment P -0.054 -0.054 -0.051 -0.051 -0.053 -0.054 -0.052 -0.052 Fe1 2.277 2.287 2.204 2.208 2.231 2.271 2.199 2.207 Fe2 1.386 1.400 1.287 1.285 1.333 1.394 1.284 1.293 Ni 0.245 0.245 0.238 0.238 0.242 0.245 0.244 0.244 27 Table 3. Harmonic vibrational frequencies (cm-1) for the bulk Fe2NiP mineral and their shifts (cm-1) due to isotopic substitutions. Color coding highlights large (red, from -10.6 to -6.9 cm-1), medium (yellow, from -3.2 to -1.1 cm-1) and small (green, from -1.0 to -0.8 cm-1) shifts. Blanked cell for shift smaller than -0.8 cm-1. No subst. 33P 58Fe1 58Fe2 60Ni 138 -0.9 160 -0.9 183 -0.9 -0.8 190 -1.9 193 -0.9 -1.7 211 -0.8 -1.2 -1.8 212 -1.2 216 -2.1 217 -1.8 -1.1 219 -1.0 -0.9 -1.3 221 -1.0 -1.3 238 -1.0 -1.0 -2.0 240 -1.6 -1.1 252 -1.1 -1.9 260 -1.4 -1.6 261 -2.2 -1.5 273 -0.8 -2.7 275 -1.0 -2.1 -1.8 279 -2.8 283 -0.8 -1.4 -1.4 -1.3 285 -1.0 -1.7 -1.6 297 -3.2 307 -3.0 -3.0 310 -1.7 -0.9 314 -1.1 -1.4 -2.2 354 -7.0 366 -6.9 -1.3 387 -8.0 -1.0 421 -9.3 -0.8 426 -10.6 434 -9.5 -1.2 441 -10.2 -0.9 453 -10.4 -0.8 467 -9.2 -1.5 -0.9 28 ASSOCIATED CONTENT Supporting Information. The following files are available free of charge. PDF file containing: computational details on computing surface formation energies, bulk neighbor analysis, VASP vs CRYSTAL17 frequency comparison, surface energy charts, Bader charge analysis of bulk and surfaces, optimized bulk and surface structures (in POSCAR format). ZIP file containing: GIF animated figures of normal modes atomic displacements. AUTHOR INFORMATION Corresponding Author *To whom correspondence should be addressed: stefano.pantaleone@unito.it piero.ugliengo@unito.it Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Funding Sources The authors acknowledge the Italian Space Agency for co-funding the Life in Space Project (ASI N. 2019-3-U.O). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 865657) for the project “Quantum Chemistry on Interstellar Grains” (QUANTUMGRAIN). 29 ACKNOWLEDGMENT AR is indebted to the “Ramón y Cajal” program. MINECO (project CTQ2017-89132-P) and DIUE (project 2017SGR1323) are acknowledged. This paper has been benefit of discussion with the members of Marie Sklodowska-Curie” Astro-Chemical Origins” (ACO) project. 30 References (1) Schwartz, A. W. Phosphorus in Prebiotic Chemistry. Philos. Trans. R. Soc. B Biol. Sci. 2006, 361 (1474), 1743–1749. https://doi.org/10.1098/rstb.2006.1901. (2) Gulick, A. Phosphorus as a Factor in the Origin of Life. Am. Sci. 1955, 43 (3), 479–489. (3) Kvenvolden, K.; Lawless, J.; Pering, K.; Peterson, E.; Flores, J.; Ponnamperuma, C.; Kaplan, I. R.; Moore, C. Evidence for Extraterrestrial Amino-Acids and Hydrocarbons in the Murchison Meteorite. Nature 1970, 228, 923–926. (4) Zhao, M.; Bada, J. L. Extraterrestrial Amino Acids in Cretaceous/Tertiary Boundary Sediments at Stevns Klint, Denmark. Nature 1989, 339, 463–465. (5) Engel, M. H.; Macko, S. A.; Silfer, J. A. Carbon Isotope Composition of Individual Amino Acids in the Murchison Meteorite. Nature 1990, 348, 47–49. https://doi.org/10.1038/348047a0. (6) Engel, M. H.; Macko, S. A. Isotopic Evidence for Extraterrestrial Non-Racemic Amino Acids in the Murchison Meteorite. Nature 1997, 389, 265–268. https://doi.org/10.1038/38460. (7) Grossman, L. Condensation in the Primitive Solar Nebula. Geochim. Cosmochim. Acta 1972, 36 (5), 597–619. https://doi.org/10.1016/0016-7037(72)90078-6. (8) Kelly, W. R.; Larimer, J. W. Chemical Fractionations in Meteorites-VIII. Iron Meteorites and the Cosmochemical History of the Metal Phase. Geochim. Cosmochim. Acta 1977, 41 (1), 93–111. https://doi.org/10.1016/0016-7037(77)90190-9. 31 (9) Lehner, S. W.; Buseck, P. R.; McDonough, W. F. Origin of Kamacite, Schreibersite, and Perryite in Metal-Sulfide Nodules of the Enstatite Chondrite Sahara 97072 (EH3). Meteorit. Planet. Sci. 2010, 45 (2), 289–303. https://doi.org/10.1111/j.1945-5100.2010.01027.x. (10) Pasek, M. A. Phosphorus as a Lunar Volatile. Icarus 2015, 255, 18–23. https://doi.org/10.1016/j.icarus.2014.07.031. (11) Pasek, M.; Lauretta, D. Extraterrestrial Flux of Potentially Prebiotic C, N, and P to the Early Earth. Orig. Life Evol. Biosph. 2008, 38, 5–21. https://doi.org/10.1007/s11084-007-9110-5. (12) Pasek, M. A. Rethinking Early Earth Phosphorus Geochemistry. Proc. Natl. Acad. Sci. U. S. A. 2008, 105 (3), 853–858. https://doi.org/10.1073/pnas.0708205105. (13) Pasek, M. A.; Harnmeijerbc, J. P.; Buickb, R.; Gulla, M.; Atlasa, Z. Evidence for Reactive Reduced Phosphorus Species in the Early Archean Ocean. Proc. Natl. Acad. Sci. U. S. A. 2013, 110 (25), 10089–10094. https://doi.org/10.1073/pnas.1303904110. (14) Pasek, M. A. Schreibersite on the Early Earth: Scenarios for Prebiotic Phosphorylation. Geosci. Front. 2017, 8 (2), 329–335. https://doi.org/10.1016/j.gsf.2016.06.008. (15) Hess, B. L.; Piazolo, S.; Harvey, J. Lightning Strikes as a Major Facilitator of Prebiotic Phosphorus Reduction on Early Earth. Nat. Commun. 2021, 12, 1–8. https://doi.org/10.1038/s41467-021-21849-2. (16) Pasek, M. A.; Lauretta, D. S. Aqueous Corrosion of Phosphide Minerals from Iron Meteorites: A Highly Reactive Source of Prebiotic Phosphorus on the Surface of the Early Earth. Astrobiology 2005, 5 (4), 515–535. https://doi.org/10.1089/ast.2005.5.515. (17) La Cruz, N. L.; Qasim, D.; Abbott-Lyon, H.; Pirim, C.; McKee, A. D.; Orlando, T.; Gull, 32 M.; Lindsay, D.; Pasek, M. A. The Evolution of the Surface of the Mineral Schreibersite in Prebiotic Chemistry. Phys. Chem. Chem. Phys. 2016, 18 (30), 20160–20167. https://doi.org/10.1039/c6cp00836d. (18) Gull, M.; Mojica, M. A.; Fernández, F. M.; Gaul, D. A.; Orlando, T. M.; Liotta, C. L.; Pasek, M. A. Nucleoside Phosphorylation by the Mineral Schreibersite. Sci. Rep. 2015, 5, 2–7. https://doi.org/10.1038/srep17198. (19) Pasek, M. A.; Dworkin, J. P.; Lauretta, D. S. A Radical Pathway for Organic Phosphorylation during Schreibersite Corrosion with Implications for the Origin of Life. Geochim. Cosmochim. Acta 2007, 71 (7), 1721–1736. https://doi.org/10.1016/j.gca.2006.12.018. (20) Pasek, M.; Herschy, B.; Kee, T. P. Phosphorus: A Case for Mineral-Organic Reactions in Prebiotic Chemistry. Orig. Life Evol. Biosph. 2015, 45 (1–2), 207–218. https://doi.org/10.1007/s11084-015-9420-y. (21) Qasim, D.; Vlasak, L.; Pital, A.; Beckman, T.; Mutanda, N.; Abbott-Lyon, H. Adsorption of Water, Methanol, and Formic Acid on Fe2NiP, a Meteoritic Mineral Analogue. J. Phys. Chem. C 2017, 121 (25), 13645–13654. https://doi.org/10.1021/acs.jpcc.7b01312. (22) Bryant, D. E.; Kee, T. P. Direct Evidence for the Availability of Reactive, Water Soluble Phosphorus on the Early Earth. H-Phosphinic Acid from the Nantan Meteorite. Chem. Commun. 2006, 4 (22), 2344–2346. https://doi.org/10.1039/b602651f. (23) Pallmann, S.; Šteflová, J.; Haas, M.; Lamour, S.; Henß, A.; Trapp, O. Schreibersite: An Effective Catalyst in the Formose Reaction Network. New J. Phys. 2018, 20, 055003. 33 https://doi.org/10.1088/1367-2630/aabb99. (24) Wang, Y.; Kong, B.; Zhao, D.; Wang, H.; Selomulya, C. Strategies for Developing Transition Metal Phosphides as Heterogeneous Electrocatalysts for Water Splitting. Nano Today 2017, 15, 26–55. https://doi.org/10.1016/j.nantod.2017.06.006. (25) Kim, B.; Kim, T.; Lee, K.; Li, J. Recent Advances in Transition Metal Phosphide Electrocatalysts for Water Splitting under Neutral PH Conditions. ChemElectroChem 2020, 7 (17), 3578–3589. https://doi.org/10.1002/celc.202000734. (26) Weng, C. C.; Ren, J. T.; Yuan, Z. Y. Transition Metal Phosphide-Based Materials for Efficient Electrochemical Hydrogen Evolution: A Critical Review. ChemSusChem 2020, 13 (13), 3357–3375. https://doi.org/10.1002/cssc.202000416. (27) Shi, Y.; Li, M.; Yu, Y.; Zhang, B. Recent Advances in Nanostructured Transition Metal Phosphides: Synthesis and Energy-Related Applications. Energy Environ. Sci. 2020, 13 (12), 4564–4582. https://doi.org/10.1039/d0ee02577a. (28) Feng, L.; Xue, H. Advances of Transition Metal Phosphide Application in Electrochemical Energy Storage and Catalysis. ChemElectroChem 2017, 4, 20–34. (29) Dinh, K. N.; Liang, Q.; Du, C. F.; Zhao, J.; Tok, A. I. Y.; Mao, H.; Yan, Q. Nanostructured Metallic Transition Metal Carbides, Nitrides, Phosphides, and Borides for Energy Storage and Conversion. Nano Today 2019, 25, 99–121. https://doi.org/10.1016/j.nantod.2019.02.008. (30) Gu, T.; Fei, Y.; Wu, X.; Qin, S. High-Pressure Behavior of Fe3P and the Role of Phosphorus in Planetary Cores. Earth Planet. Sci. Lett. 2014, 390, 296–303. 34 https://doi.org/10.1016/j.epsl.2014.01.019. (31) Sagatov, N. E.; Gavryushkin, P. N.; Banayev, M. V.; Inerbaev, T. M.; Litasov, K. D. Phase Relations in the Fe-P System at High Pressures and Temperatures from Ab Initio Computations. High Press. Res. 2020, 40 (2), 235–244. https://doi.org/10.1080/08957959.2020.1740699. (32) Zhao, Z.; Liu, L.; Zhang, S.; Yu, T.; Li, F.; Yang, G. Phase Diagram, Stability and Electronic Properties of an Fe-P System under High Pressure: A First Principles Study. RSC Adv. 2017, 7 (26), 15986–15991. https://doi.org/10.1039/C7RA01567D. (33) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47 (1), 558–561. (34) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal- Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49 (20), 14251– 14269. (35) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B - Condens. Matter Mater. Phys. 1996, 54 (16), 11169–11186. (36) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15–50. (37) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented - Wave Method. Phys. Rev. B 1999, 59 (3), 1758–1775. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. 35 Phys. Rev. Lett. 1996, 77 (18), 3865–3868. (39) Janthon, P.; Kozlov, S. M.; Viñes, F.; Limtrakul, J.; Illas, F. Establishing the Accuracy of Broadly Used Density Functionals in Describing Bulk Properties of Transition Metals. J. Chem. Theory Comput. 2013, 9 (3), 1631–1640. https://doi.org/10.1021/ct3010326. (40) Janthon, P.; Luo, S.; Kozlov, S. M.; Viñes, F.; Limtrakul, J.; Truhlar, D. G.; Illas, F. Bulk Properties of Transition Metals: A Challenge for the Design of Universal Density Functionals. J. Chem. Theory Comput. 2014, 10 (9), 3832–3839. https://doi.org/10.1021/ct500532v. (41) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27 (15), 1787–1799. (42) Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. B3LYP Augmented with an Empirical Dispersion Term (B3LYP-D*) as Applied to Molecular Crystals. Cryst. Eng. Comm. 2008, 10 (4), 368–376. (43) Cutini, M.; Maschio, L.; Ugliengo, P. Exfoliation Energy of Layered Materials by DFT-D: Beware of Dispersion! J. Chem. Theory Comput. 2020, 16 (8), 5244–5252. https://doi.org/10.1021/acs.jctc.0c00149. (44) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (45) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465. 36 (46) Togo, A.; Tanaka, I. First Principles Phonon Calculations in Materials Science. Scr. Mater. 2015, 108, 1–5. https://doi.org/10.1016/j.scriptamat.2015.07.021. (47) Togo, A.; Chaput, L.; Tanaka, I. Distributions of Phonon Lifetimes in Brillouin Zones. Phys. Rev. B - Condens. Matter Mater. Phys. 2015, 91 (9), 094306. https://doi.org/10.1103/PhysRevB.91.094306. (48) Skelton, J. M.; Burton, L. A.; Jackson, A. J.; Oba, F.; Parker, S. C.; Walsh, A. Lattice Dynamics of the Tin Sulphides SnS2, SnS and Sn2S3: Vibrational Spectra and Thermal Transport. Phys. Chem. Chem. Phys. 2017, 19, 12452–12465. https://doi.org/10.1039/c7cp01680h. (49) Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Rérat, M.; Casassa, S.; Baima, J.; Salustro, S.; et al. Quantum-Mechanical Condensed Matter Simulations with CRYSTAL. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, 8 (4), 1–36. https://doi.org/10.1002/wcms.1360. (50) Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Rerat, M.; Casassa, S.; Baima, J.; Salustro, B.; et al. CRYSTAL17 User Manual. Int. J. Quantum Chem. 2014, 114 (19), 1287–1317. (51) Ugliengo, P.; Viterbo, D.; Chiari, G. MOLDRAW: Molecular Graphics on a Personal Computer. Zeitschrift fur Krist. - New Cryst. Struct. 1993, 207 (1–2), 9–23. (52) Momma, K.; Izumi, F. VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44 (6), 1272–1276. https://doi.org/10.1107/S0021889811038970. 37 (53) Humphrey, W.; Dalke, A.; Schulten, K. VMD : Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. (54) Povray http://www.povray.org/. (55) Skála, R.; Císařová, I. Crystal Structure of Meteoritic Schreibersites: Determination of Absolute Structure. Phys. Chem. Miner. 2005, 31 (10), 721–732. https://doi.org/10.1007/s00269-004-0435-6. (56) Henkelman, G.; Arnaldsson, A.; Jónsson, H. A Fast and Robust Algorithm for Bader Decomposition of Charge Density. Comput. Mater. Sci. 2006, 36 (3), 354–360. https://doi.org/10.1016/j.commatsci.2005.04.010. (57) Sanville, E.; Kenny, S. D.; Smith, R.; Henkelman, G. Improved Grid-Based Algorithm for Bader Charge Allocation. J. Comput. Chem. 2007, 28 (5), 899–908. https://doi.org/10.1002/jcc. (58) Tang, W.; Sanville, E.; Henkelman, G. A Grid-Based Bader Analysis Algorithm without Lattice Bias. J. Phys. Condens. Matter 2009, 21 (8), 084204. https://doi.org/10.1088/0953- 8984/21/8/084204. (59) Grosvenor, A. P.; Wik, S. D.; Cavell, R. G.; Mar, A. Examination of the Bonding in Binary Transition-Metal Monophosphides MP (M = Cr, Mn, Fe, Co) by X-Ray Photoelectron Spectroscopy. Inorg. Chem. 2005, 44 (24), 8988–8998. https://doi.org/10.1021/ic051004d. (60) Bryant, D. E.; Greenfield, D.; Walshaw, R. D.; Evans, S. M.; Nimmo, A. E.; Smith, C. L.; Wang, L.; Pasek, M. A.; Kee, T. P. Electrochemical Studies of Iron Meteorites: Phosphorus Redox Chemistry on the Early Earth. Int. J. Astrobiol. 2009, 8 (1), 27–36. 38 https://doi.org/10.1017/S1473550408004345. (61) Bryant, D. E.; Greenfield, D.; Walshaw, R. D.; Johnson, B. R. G.; Herschy, B.; Smith, C.; Pasek, M. A.; Telford, R.; Scowen, I.; Munshi, T.; et al. Hydrothermal Modification of the Sikhote-Alin Iron Meteorite under Low PH Geothermal Environments. A Plausibly Prebiotic Route to Activated Phosphorus on the Early Earth. Geochim. Cosmochim. Acta 2013, 109, 90–112. https://doi.org/10.1016/j.gca.2012.12.043. (62) Pirim, C.; Pasek, M. A.; Sokolov, D. A.; Sidorov, A. N.; Gann, R. D.; Orlando, T. M. Investigation of Schreibersite and Intrinsic Oxidation Products from Sikhote-Alin, Seymchan, and Odessa Meteorites and Fe3P and Fe2NiP Synthetic Surrogates. Geochim. Cosmochim. Acta 2014, 140, 259–274. https://doi.org/10.1016/j.gca.2014.05.027. (63) de Faria, D. L. A.; Lopes, F. N. Heated Goethite and Natural Hematite: Can Raman Spectroscopy Be Used to Differentiate Them? Vib. Spectrosc. 2007, 45 (2), 117–121. https://doi.org/10.1016/j.vibspec.2007.07.003. (64) Boettger, J. C. Nonconvergence of Surface Energies Obtained from Thin-Film Calculations. Phys. Rev. B 1994, 49 (23), 16798–16800. https://doi.org/10.1103/PhysRevB.49.16798. (65) Boettger, J. C.; Smith, J. R.; Birkenheuer, U.; Rösch, N.; Trickey, S. B.; Sabin, J. R.; Apell, S. P. Extracting Convergent Surface Energies from Slab Calculations. J. Phys. Condens. Matter 1996, 10, 893–894. https://doi.org/10.1088/0953-8984/8/36/005. (66) Fiorentini, V.; Methfessel, M. Extracting Convergent Surface Energies from Slab Calculations. J. Phys. Condens. Matter 1996, 8, 6525–6529. https://doi.org/10.1088/0953- 8984/8/36/005. 39 (67) Wulff, G. Zur Frage Der Geschwindigkeit Des Wachstums Und Der Auflosung Der Krystallflachen. Z. Krist. 1901, 34, 449–530. (68) Zamirri, L.; Corno, M.; Rimola, A.; Ugliengo, P. Forsterite Surfaces as Models of Interstellar Core Dust Grains: Computational Study of Carbon Monoxide Adsorption. ACS Earth Sp. Chem. 2017, 1 (7), 384–398. https://doi.org/10.1021/acsearthspacechem.7b00041. 40 Table of Content Graphics (TOC) 41