Development of molecular cluster models to probe pyrite surface reactivity

The recent discovery that anaerobic methanogens can reductively dissolve pyrite and utilize dissolution products as a source of iron and sulfur to meet their biosynthetic demands for these elements prompted the development of atomic‐scale nanoparticle models, as maquettes of reactive surface sites, for describing the fundamental redox steps that take place at the mineral surface during reduction. The given report describes our computational approach for modeling n(FeS2) nanoparticles originated from mineral bulk structure. These maquettes contain a comprehensive set of coordinatively unsaturated Fe(II) sites that are connected via a range of persulfide (S22−) ligation. In addition to the specific maquettes with n = 8, 18, and 32 FeS2 units, we established guidelines for obtaining low‐energy structures by considering the pattern of ionic, covalent, and magnetic interactions among the metal and ligand sites. The developed models serve as computational nano‐reactors that can be used to describe the reductive dissolution mechanism of pyrite to better understand the reactive sites on the mineral, where microbial extracellular electron‐transfer reactions can occur.


| INTRODUCTION
Pyrite (FeS 2(p) ) or "fool's gold", is the most abundant iron-sulfur (Fe-S) mineral in Earth's crust. 1,2It can host a number of transition metals of economic importance, including nickel, cobalt, and molybdenum. 3In the absence of oxidative weathering, FeS 2(p) has long been thought to be a sink for Fe and S even in strong acidic conditions. 4In turn, this was thought to render FeS 2(p) unavailable as an energy source or source of elements for biology.However, several recent studies have shown that FeS 2(p) can be reductively dissolved through both abiotic and biotic processes at anoxic, physiological conditions. 5For the former, dissolved hydrogen gas in micromolar concentrations and at biologically relevant temperatures (approximately 40 C) can abiotically reduce FeS 2(p) , resulting in what has been hypothesized to be a pyrrhotite (Fe 1-x S (p) ) surface-associated phase and release of HS À (aq) into solution. 6][10][11][12][13] The insoluble nature of FeS 2(p) indicates 14 that reduction by methanogens takes place extracellularly either via direct contact with the mineral surface to catalyze reduction 15 or with the asistance of soluble redox shuttles for long range reduction, such as quinone analogs (i.e., anthraquinone-2,6-disulfonate). 16Mineral reduction by methanogens does not appear to be energy-yielding and rather is carried out to generate bioavailable Fe and S in order to meet the cells' biosynthetic demands for these elements.Collectively, these observations require a shift in how we consider the stability of FeS 2(p) and the mobilization of elements from this mineral in anoxic environments, including the vast sedimentary marine basins, the source of most FeS 2(p) produced today. 1,179][20][21][22][23][24] Focusing on the Fe-S 2 unit of the mineral bulk (Figure 1; Figure S1), the coordination environment of Fe 2+ ion can be described as slightly distorted octahedral, while each S center is surrounded by three Fe 2+ ions and covalently bound to the other S of the persulfide ligand in a tetrahedral arrangement.The similarity in energy and shape of the Fe 3d and S 3p orbitals creates a rich ensemble of electronic structure features that are present in the form of S-ligand to Fe donation (S !Fe), Fe to S-ligand backdonation (Fe !S) as covalent interactions, and Fe… S/Fe…Fe/S…S ionic interactions.In order to capture these structural features, a common approach for describing mineral surfaces at the atomic level is using slab models that are periodic with respect to the basal plane.Representative surface models are shown in Figure 1 that highlight a range of coordination environments of the surface-exposed Fe and S ions.The four dominant faces for FeS 2(p) are the cubic (100), octahedral (111), and pyritohedral (102 and 210) crystal morphologies.
As shown in the surface models in Figure 1, we can identify 5-, 4-, and 3-coordinate Fe sites with square pyramidal, seesaw, and trigonal pyramidal coordination environments, respectively.From coordination chemistry principles, we anticipate the crucial importance of the unsaturated Fe and S sites in reactivity toward substrate binding (physisorption), activation (chemisorption), and concomitant bond breaking/formation (reactivity).
F I G U R E 1 Overview of the most dominant FeS 2(p) surface planes (generated using CSD Mercury 25 ) as macroscopically displayed by the specimens (top), schematic representations of the three crystal morphologies (yellow polygons), and atomic-scale chemical environments of Fe (dark red atoms) and S sites (yellow golden atoms) at (100), (111), (102), and (210) surfaces with varying degrees of displacement to expose 5-, 4-, and 3-coordinate Fe sites with square pyramidal, seesaw, and trigonal pyramidal coordination geometries, respectively, as well as persulfide-rich termination of the mineral surface.
During the reductive dissolution of FeS 2(p) , formation of Fe-bound hydride, protonated persulfide, bisulfide/hydrogen sulfide intermediates are expected.The examination of these intermediates and related transient species with periodic boundary condition (PBC) models would pose prohibitive computational costs and severe limitations for the applied level of theory.Instead, we advocate for the development of molecular cluster models that feature the reactive chemical environments; thus, forming maquettes for the Fe sites and persulfide-rich regions.Additional advantage of using FeS 2(p) nanoparticle maquettes is the utilization of a broad range of potential energy surface mapping methodologies for describing localized reactivity and molecular orbital-based electronic structure at affordable computational costs, while employing a high level of theory (range-separated hybrid functionals with dispersion corrections) and accessing a large variety of chemical environments (including polarizable continuums) for the same model.We already demonstrated the benefits of computational maquettes for aluminosilicates 26 and carbon materials. 27e novelty of this work is a maquette construction strategy, where instead of identifying inner and outer spheres of interactions around a central unit, we construct nanoparticle models using a repeat pattern (a coarse grain) of structural units.For example, in the mackinawite mineral (FeS (m) ), each mineral layer can be built up from [2Fe-2S] rhombs.9][30][31] The FeS 2 (p) crystal structure offers a structural unit with 2(FeS 2 ) pentagons that contain 1,2-μ 2 (end-on) and 1,1-μ 2 (side-on, dangling) bridging persulfide ligands.These pentagonal units can be expanded into nanoparticles leading to n(FeS 2 ) nanoparticles with n = 8, 18, and 32 defined as generation 1, 2, and 3 maquettes herein, respectively.The specific compositions (n) were selected to maintain consistency with layered, mackinawite-like nanoparticles of n(FeS) (aq) stoichiometry, as plausible final products of the overall reductive dissolution process.
The n(FeS) (m) mackinawite nanoparticles are composed of single phase.

| METHODS
All calculations were carried out using the Gaussian16 suite of programs. 32The employed density functional theory was first selected by comparing the experimental crystal structure (proximal Fe-S p , distal Fe…S d , S-S, and Fe…Fe distances, and lattice constants) with their calculated counterparts (see Supporting Information Analysis 1).The survey of the Inorganic Crystal Structure Database (ICSD) 33 resulted in several entries for pyrite from different measurements under various conditions (Table S1).In our work we used the structure with the lowest R factor as reference crystal structure (ICSD#109377). 19The GGA and meta-GGA functionals selected for validation are shown in Table S2.
Overall, the most reasonable model for FeS 2(p) was obtained using the meta-GGA M15L functional 34 with negligible structural differences between double and triple-ζ quality basis sets.However, upon further evaluation for molecular n(FeS 2 ) nanoparticles, we identified chronic convergency issues that could be mitigated by switching to the MN15 hybrid meta-GGA functional. 35A detailed description of the basis set dependence and the functional evaluations in comparison to ab initio correlated molecular orbital calculations are given in Supporting Information Analysis 2 (Figures S4-S8, Table S3).We alternated between triple-ζ (def2TZVP) and double-ζ (def2SVP) basis sets 36 for smaller and larger models, respectively, as has been evaluated for Fe-S systems earlier. 37Given the presence of coordinatively unsaturated sites, all models were embedded into a solvent density-based (SMD) polarizable continuum model 38 parameterized for ambient aqueous conditions.The reference level herein is MN15/def2SVPjSMD unless specified otherwise.
The non-standard, ultra-fine integration grid was used.Spin densities are from Mulliken population analysis. 39A consistent color scheme was adopted with red, blue, and gray colors representing excess spin up/α, spin down/β densities on the Fe centers and negligible spin density on the S centers, respectively.Open-shell, diamagnetic spin states (S t = 0, where S t refers to broken-symmetry calculations 40,41 with contributions from all possible m s = 0 sub-levels) were generated with antiferromagnetic-(AFC) and ferromagnetic-coupled (FC) high-spin Fe 2+ ions in diamagnetic ground states in comparison to models containing only closed-shell diamagnetic ferrous centers (S = 0).The AFC and FC states were constructed by merging Fe 2+ and S 2 2À ionic fragments. 37The spin densities at contour level of 0.0035 a.u.Fe 2+ ions at the surface-exposed sites.We carried out a set of single point energy calculations for mononuclear models with the [Fe (S 2 ) x ] 2-2x composition in crystallographic positions (Figure S2) in order to evaluate the energetic ordering of LS, IS, and HS states.In addition to the model composition, we also evaluated the impact of using hybrid GGA functionals by considering a gradually varied amount of the HF exchange term using the PW91 functionals, [42][43][44] which showed decent bulk structure description (Table S2).
This level of theory assessment is relevant for bridging the pure MN15L functional employed for the highest chemical accuracy for periodic FeS 2(p) models versus the meta-hybrid functional (MN15) for achieving convergency for the molecular cluster models.
The lowest energy spin states (entries marked with bold 0.0 values) in Table 1 for the coordinatively unsaturated Fe environments overwhelmingly correspond to the HS(S = 2) state independently of the amount of HF exchange applied.The diamagnetic LS(S = 0) state becomes stable as the number of ligands and the mixing of the ionic character (HF exchange) increases at the expense of the covalent delocalization (DF exchange).The 5-coordinate model with pure GGA and hybrid functionals with less than 25% HF exchange shows energetic preference for the IS(S = 1) state.The same trends can be obtained as well, when using a smaller def2SVP basis set.It is then important to underscore that the surface-exposed, coordinatively unsaturated Fe sites manifest a significantly different electronic structure than that of the bulk.Thus, we must consider the presence of paramagnetic centers at the FeS 2(p) surface and build up shell-like models interfacing the paramagnetic surface and diamagnetic bulk.sphere, since all levels of theory including the pure HF method with a compact def2SVP basis set show LS to be the ground state (Figure S3).
Furthermore, Ca 2+ and Na + as neutralizing counterions in the position of the outer sphere Fe 2+ ions from the bulk FeS 2(p) structure show preference for the low-spin state only at low HF exchange mixing or in pure GGA functionals.Protons, unless their positions are optimized as in [Fe(S 2 H 2 ] 6 ] 2+ , and point charges, even in the presence of polarizable continuum, do not favor the LS state for either hybrid or GGA functionals.Therefore, it is important to highlight that, if the composition of the computational model captures the chemical bonding description, the predicted energetic order of the spin states will be accurate regardless of whether hybrid or pure GGA functional is employed. A brief summary of additional modeling efforts to create coordination models for a central Fe-S 2 unit with inner, outer, and peripheral spheres of ions is provided in Supporting Information Analysis 3, Figures S9-S11 and Tables S4 and S5.Given the lack of chemical similarity of the spherical nanoparticles and the surface-exposed Fe sites in Figure 1, we focused on models with 2(FeS 2 ) pentagon units, as coarse grain building blocks.

| Magnetic interactions in 2(FeS 2 ) models
The pentagonal 2(FeS 2 ) building block of the FeS  Note: The PW91 functional set [42][43][44] was utilized for both exchange and correlation interactions with the def2TZVP basis set 36 due to the good agreement in describing the bulk periodic model (Table S2).All complexes were embedded into the SMD implicit solvation model 38 with water parameters.F I G U R E 2 Initial (top row) and optimized (bottom row) structures of 2(FeS 2 ) models at the MN15/def2TZVPjSMD level of theory with selected geometric parameters, atomic spin densities (red: spin up/α; blue: spin down/β densities), and spin density contour plots for the lowest energy spin-coupling schemes.The coupling constants (J) were determined using Yamaguchi's spin projection method (see Table S6 for details). 40he relative energies for the optimized structures were calculated from electronic energies (ΔE SCF ).
of the optimized structures reveal two energetically close arrangements and a chair/boat conformational difference of 18 kJ mol À1 .
Side-on persulfides at a single Fe site or bridged between multiple Fe sites are expected for the FeS 2(p) surface.It is also important to note that the initial and optimized structures differ considerably.This highlights the considerable differences between the bulk and surface Fe-S, Fe…S, Fe…Fe, S-S, and S…S distances and related bond angles.
Furthermore, the structural reorganization has a strong impact on the energetic preference of the spin-coupling schemes.These pose caveats when we worked with larger maquettes (vida infra), including the relative energy differences of spin-coupling schemes among initial, bulk FeS 2 -based crystal structures, and their optimized counterparts.
Larger structural displacements for smaller maquettes correspond to greater deviations in the relative order of various spin states as illustrated in Figure 2. The initial structures of the four-and five-membered ring models show the AFC state to be the most stable, while the sixmembered ring ferromagnetic (FC) state becomes the lowest in energy due to lack of an efficient M-L-M super-exchange path in the 1,2-bridging persulfide with M-L-L'-M pattern.The structural optimization reorders this energetic preference, since the short Fe…Fe distance of 2.66 Å (bottom left corner, Figure 2) enables the direct exchange interaction between the paramagnetic Fe centers.This renders the four-membered ring structure to have FC state, while the other two favor the AFC state.The spin-coupling constants vary from À245 cm À1 to +191 cm À1 with the AFC interactions outcompeting the FC interaction as can be seen from the initial five-membered ring model.
We generated 2(FeS 2 ) models with explicit water solvent molecules at each Fe site for addressing the influence of unsaturated coordination environment.Figure S12 summarizes the initial fivemembered ring models, their optimized counterparts, and the relative energies in the AFC state.Neither the tetra-nor hexa-coordinate environments prevented the rearrangements of the initial bulk FeS 2 crystal structure as the position of persulfide ligands significantly shifted.This suggests that the role of explicit solvent water molecules is expected to be negligible in comparison to the Fe 2+ /Fe 2+ , Fe 2+ /S 2 2À , and S 2 2À / S 2 2À interactions.Thus, we focused on expanding our study to tetranuclear clusters with an increased number of Fe-S 2 interactions.
Figure 3 shows the large-scale changes among the initial and geometry optimized structures of the flattened 4(FeS 2 ) model with one dangling persulfide and three conjoined pentagons.
The quasi-planar pentagons folds up to an egg-shaped cluster with features that are unprecedented in the FeS 2 bulk, such as the short Fe…Fe distance of 2.80 Å.The double end-on coordination of the persulfide (S11 and S12) in Figure 3B allows for the direct exchange to dominate the magnetic bonding between the FC pair of Fe4 and Fe10, as also seen for 2(FeS 2 ) models in Figure 2. Due to the different arrangements of the persulfide ligands at the opposite sides of the cluster, no similar interaction can be identified between Fe1 and Fe7 centers.The optimized structures in Figure 3  The pure spin state from all FC Fe centers with S = 8 was found to be 18 kJ mol À1 higher in energy than the S t = 0 ground state in Figure 3D.The expectation value of the spin operator of <S 2 > = 72.05for the S(S + 1) value indicated close to pure high-spin state.Thus, the pentagonal prism model for the 4(FeS 2 ) can already be considered as a minimalist building block for larger nanoparticles.

| Maquettes with 8(FeS 2 ) composition
The first generation (G1) of molecular cluster model as a computational nano-reactor is the 8(FeS 2 ) maquette (Figure 4).It features a diverse set of Fe and S chemical environments and exhibits structural stability of its framework.The initial structures were strictly based on the FeS 2(p) bulk (Figure 4A).The 8(FeS 2 ) model can be described as two stacked 4(FeS 2 ) planar sheets (Figure 3A) or three fused 4(FeS 2 ) prisms (Figure 3C) that impacts the topology of spin-coupling schemes.The coordination numbers of Fe centers vary from 3 to 5 with trigonal, seesaw, or tetragonal pyramidal environments as seen in the slab models of Figure 1.Similarly, the persulfide ligands vary from being dangling (S23-S24) or bridging that connect 3 + 1, 2 + 1, or 2 + 2 Fe centers, respectively.
The diamagnetic, closed-shell electronic structure for the "surface-only" maquette without a 6-coordinate, bulk Fe site is more than 1 MJ mol À1 higher energy than the FC structure in Figure 4B.On the contrary, an arbitrarily selected AFC structure confirms the essential need for considering models containing Fe sites with alternating signs of spin densities via broken-symmetry density functional approach.
The importance of a balanced distribution of FC/AFC pairs for achieving the lowest energy structure has already been shown for the 4 (FeS 2 ) models (vida supra).In Figure 4C, the central/bottom and left prism have more FC pairs than AFC pairs, which suggests that lower energy structures may exist when considering alternative spincoupling schemes.In the following, we present a set of rules for approximating the lowest energy structure and validating these by considering all possible spin-coupling combinations.

| Atom-by-atom spin-flip search
The approach of evaluating the energy impact of altering the α/β-spin densities of each Fe center takes advantage of the well-defined nature of the FC structure (Figure 4B) with the pure spin state of S = 16.The orange triangles in Figure 5 show the energy consequence of lowering the spin state to S t = 12 (<S 2 > ≈ 160.0) containing a single high-spin Fe site with negative spin density (m s = À2).
A single spin flip causes the energy to drop by at least 15 kJ mol À1 with a significant preference for the Fe of the highest coordination number (tetragonal pyramidal, Fe4) to be AFC (À37 kJ mol À1 ) in comparison to other Fe centers (À20 kJ mol À1 on average).The S = 1 intermediate-spin state was found to be the ground state for the square pyramidal mononuclear complex using hybrid and pure GGA functionals in Table 1.Therefore, we considered the energetic preference of Fe4 being in the IS state.The calculation  of the S t = 14 (<S 2 > =210.8)state was found to be even higher (+143 kJ mol À1 ) than the initial pure FC state and thus Fe sites with m s = ±1 IS state will not be considered.Continuing the search (light blue circles, Figure 5), we altered the spin for four-and threecoordinate Fe centers and found that, for the S t = 8 state (<S 2 > = 79.9), the lowest energy structure was obtained when the tetracoordinate Fe16 (seesaw ligand environment) formed multiple AFC interactions with the rest of the Fe ions at À52 kJ mol À1 relative to the FC state.Notably, this state is already lower energy than the arbitrarily constructed structure in Figure 4C.The other tetracoordinate site (Fe22) is the second most favorable Fe center for being involved in AFC interactions.As seen in Figure 5 for the S t = 4 state (purple square, <S 2 > = 31.8),the spin flip is energetically favorable for the trigonal Fe19 site (À57 kJ mol À1 ).
In contrast, the lowest energy structure (À59 kJ mol À1 ) was localized by inverting the spin on Fe4, Fe13, Fe16, Fe19 from the FC structure to S t = 0 open-shell, diamagnetic state (<S 2 > = 15.9) with maximized number of 14 AFC and minimized number of 2 FC interactions (Figure 5, green diamonds).The lowest energy 8(FeS 2 ) nanoparticle from the atom-by-atom spin-flip approach provides us the most significant guideline with clear preferences for the highest coordinated Fe site to be involved in AFC interactions and for a greater number of AFC than FC pairs.When we considered the prism-based topology with pentagonal sides of the 8(FeS 2 ) model, the equal division of the AFC/FC pairs improved the atom-by-atom spin-flip approach in comparison to the arbitrary AFC structure shown earlier in Figure 4C.However, only the prism with Fe1-Fe4-Fe19-Fe22 and the Fe4-Fe7-Fe16-Fe22 have the equal number of αand β-spin centers, while the third prism with Fe4-Fe10-Fe13-Fe16 has more β-spin than α-spin centers.The evaluation of all 35 S t = 0 spincoupling schemes (Supporting Information Analysis 4, Table S7) confirms that it is not possible to have all prisms containing the same number of αand β-spin centers due to the shape of the 8(FeS 2 ) cluster.
Upon optimization of the three lowest energy structures (within 10 kJ mol À1 , Figure S13), we found that the relative energetic order of the crystal structure-based, initial, and optimized structures was structure.This is the result of decreased backdonation and thus decreased Fe-S bond order and vice versa increased S-S bond order.The large variability in the electronic structure between the surface-exposed and bulk Fe sites provides another rationale for using our maquettes in comparison to periodic surface models (Figure 1) with boundary conditions that limit the presence of spin polarization and structural elasticity.

| Maquettes with 18(FeS 2 ) composition
Our proposal for a second generation (G2) maquette is described by three stacked 6(FeS 2 ) layers or 12 fused 4(FeS 2 ) prisms (Figure 7).The novelty of the structure in comparison to G1 is the presence of a single bulk Fe site (Fe34, Figure 7) that has the full hexakis(persulfido) coordination environment.Similarly, the S25-S26 persulfide ligand has a complete double tetrahedral environment with three Fe sites at each S. As shown above for the mononuclear [Fe(S 2 ) x ] 2-2x models, the ground state for the bulk Fe site is closed-shell diamagnetic (S = 0).
The odd number of the surface Fe sites will render the total spin state to be S t = 2 due to the presence of an uncompensated high-spin Fe 2+ site with four unpaired electrons (m s = +2).We created various spincoupling schemes with mixed FC electronic structure, where the lowspin Fe34 is surrounded by seven high-spin Fe sites in FC-only arrangement (Figure 7A), a layered AFC arrangement where the bottom 6 and 2 Fe sites from the right-hand side were spin flipped (Figure 7B), in addition to the guided AFC structure (Figure 7C), as the first application of the preference of spin-coupling schemes developed for the G1, 8(FeS 2 ) maquette.
In the lowest energy G2 model (Figure 7C), the Fe sites with the highest coordination number (Fe31) are involved in the greatest number of AFC interactions, the total number of AFC pairs are maximized as can be seen from 6 AFC and 1 FC in Figure 7C, in contrast to 3 AFC and 4 FC at Fe31 center in Figure 7B.Furthermore, the distribution of αand β-spin sites are adjusted as we walk through square pyramidal, sea-saw, and pyramidal coordination environment from the mixed FC state with a diamagnetic bulk Fe site.The importance of maximizing the number of AFC pairs is further substantiated by Figure 7D, which plots the four different numbers of AFC/FC pairs that we can define for the 18(FeS 2 ) model.Since multiple spin-coupling schemes exist for the maximized number of AFC pairs, the guided approach cannot result in a unique solution.Without considering all the permutations for selecting the highest coordination environments, we could not localize the global minimum in the electronic structure space for the initial structure.
Upon geometry optimization of the 18(FeS 2 ) maquettes with AFC/FC pairs based on the guidelines developed, we obtained a stable nanoparticle without large-scale structural arrangements as shown in Figure 8.The RMSD of the optimized structure and the initial structure is 0.039 Å, which is remarkably small.This can be attributed to the presence of a complete coordination environment of the central Fe34-S35-S26 unit.This model is proposed to be a realistic computational nano-reactor for investigating the interconversion of bulkto-surface Fe sites and the formation of a plausible surface-associated Fe 1-x S (p) phase.It is important to highlight that the Fe-S distances increased by 0.11-0.24Å around the bulk center relative to the FeS 2 (p) crystal structure, while the central S-S distances decreased as an already described pattern for the 8(FeS 2 ) model (vide infra).This deviation further reinforces the concept that surface-exposed Fe sites behave significantly differently from the bulk Fe sites.
It is also noteworthy that out of the 12 prisms composing the 18(FeS 2 ) model, 8 prisms involve the bulk Fe site (right-hand side, Figure 8), but each of them has 2AFC/1FC or 1AFC/2FC pairs; while 2 of the remaining prisms (left-hand side in Figure 8) have balanced number of AFC/FC pairs.
F I G U R E 6 Equilibrium structures for the three most stable 8(FeS 2 ) nanoparticles at the MN15/def2SVPjSMD level of theory.Energies (ΔE SCF ) are given relative to the initial structure of the FC state.
Therefore, based on the 8(FeS 2 ) G1 and 18(FeS 2 ) G2 models, we further refined the guidelines for obtaining a low energy spincoupling scheme for the minimal S t state, without systematically going through a full combinatorial spin-flip search among the unique spin-coupling schemes.The construction steps are as follows: 1. Hexakis(persulfido) iron sites with coordinatively saturated persulfide ligands have low-spin state (bulk Fe-S 2 sites); 2. Fe sites with the highest coordination number must be involved in antiferromagnetic interactions; 3. The total number of antiferromagnetic interactions needs to be maximized; 4. The spin up (α, m s = +2) and spin down (β, m s = À2) centers need to be distributed to maximize the number of balanced AFC/FC pairs in prisms with pentagonal sides; and 5.The presence of side-on coordinated persulfide ligand may increase cluster stability.

| Maquettes with 32(FeS 2 ) composition
The third generation (G3) model is our most realistic FeS ) units with a coordinative saturated environment for both Fe and S sites in low-spin environment in all models (hence the term "mixed").Among the numerous S t = 0 broken-symmetry solutions, we highlight the layered, blockish, and guided mixed AFC state as the lowest energy structure obtained for the G3 model of egg-shaped, spherical morphology.The energies (ΔE SCF ) are relative to the mixed FC state at the top.
which is notably small for close to 100 atoms in comparison to the energy changes observed during the relaxation of smaller models.
When considering the prism-based topology of the model, the shell-like structure for the high-spin outer surface and the core bulk structure reduces the number of complete 4(FeS 2 ) prisms to 8, among which only 3 do not have a balanced number of AFC and FC pairs.When completely removing the 6(FeS 2 ) core (Figure 10B), surface-exposed Fe-S 2 units form a continuous monolayer around the bulk core (Figure 10C).Overall, the third generation model

| DISCUSSION
The main deliverable of our systematic study is the set of three generations of computational nano-reactors future investigations of FeS 2(p) surface reactivity.While our motivation is to probe the reductive dissolution of FeS 2(p) , it is important to highlight that these models are also ideal for studying the atomic-scale mechanism(s) of oxidative of FeS 2(p) [48] The first generation model (G1) contains 8 surface-exposed Fe sites with high-spin, antiferromagnetic spin/spin coupling (Figure S16).The second generation model (G2) contains a single bulk Fe-S 2 moiety with 17 surface-exposed Fe sites, where the Fe 2+ and S 2 2À ions are in three layers with maximized number of AFC pairs.Fe sites of the highest coordination number are preferentially involved in AFC interactions (Figure S17).The third generation and most realistic F I G U R E 1 0 Optimized third generation, 32(FeS 2 ) model (S t = 0 with 6(FeS 2 ) bulk core in low spin, m s = 0 state) with selected bond lengths obtained by following the spin-coupling guidelines for a low energy structure at the MN15/def2SVPjSMD level of theory.The energy stabilization due to optimization is 240 kJ mol À1 (ΔE SCF ) and is given relative to the mixed FC state in Figure 9.
Using the computational models (see summary of n(FeS) (aq) nanoparticles in Figure S16 and Table S9), we can estimate the overall reaction thermodynamics.These thermochemical calculations define the overall free energy and enthalpy (in parenthesis) for n

(n = 8 )
, double (n = 18), or triple (n = 32) rings of [2Fe-2S] rhombs that surround a central rhomb.Note that the overall thermodynamics of the reduction 6 in Equation (1) excludes an intermediate step of the formation of a surface-associated Fe 1-x S (p) phase.The main outcome of employing the maquette construction study is the set of computational nano-reactors that feature Fe and S sites closely mimicking the surface formations of FeS 2(p) from Figure 1.The molecular cluster building strategy can be generalized for creating Fe-S nanoparticles of other Fe-S mineral phases.The computational nano-reactors open the possibility for systematically mapping the energetically favorable sites for interactions with reducing agents, formation of intermediates during the proton-and electron-transfer processes, and the release of the sulfide/bisulfide ions in the reductive dissolution process of FeS 2(p) to FeS (m) though a postulated Fe 1-x S (p) Considering coordination chemistry principles around a central Fe-S 2 unit, Figure S3 details various termination strategies for the [Fe (S 2 ) 6 ] 10À , a highly charged model to address the impact of uncompensated charge.The cooperativity between the S !Fe donation and Fe !S backdonation interactions is well demonstrated by the models containing additional Fe 2+ or [Fe(OH 2 )] 2+ ions in the outer 2(p) crystal structure shown in the middle of Figure 2 inspired the creation of a rhomb (left-hand side, as in biological Fe-S clusters) and a hexagon (right-hand side) shaped binuclear complex for the evaluation of the FC and AFC magnetic bonding among the m s = ±2 Fe centers.The three different ring systems are plausible for FeS 2(p) surface formations when considering the bis(1,1-μ 2 ), mixed (1,1-μ 2 ) and (1,2-μ 2 ), and bis(1,2-μ 2 ) bridge of the persulfide ligands.Starting from the FeS 2 crystal structure (hence the identical bond lengths, interatomic distances, and bond angles in the top row of T A B L E 1 Relative energies (ΔE SCF , eV) of high-spin (HS, S = 2), intermediate-spin (IS, S = 1), and low-spin (LS, S = 0) states for [Fe(S 2 ) x ] 2-2x mononuclear complexes (Figure S2) as a function of the number of persulfide ligands and the amount of HF exchange (HFX) [Fe(S 2 ) x ]

Figure 2 )
Figure 2), the optimized structures (bottom row) clearly indicate how persulfide manifests a remarkable flexibility when attached to coordinatively unsaturated Fe sites.The two distinctly different equilibrium structures for the 2(FeS 2 ) models are the testament to the high electrophilicity and coordination flexibility of the HS Fe sites and S 2 2À

4 (
Figure 3C, the energetically most stable form for the non-planar 4 (FeS 2 ) nanoparticle model is shown in Figure 3D.After repeated optimization steps for eliminating the imaginary normal modes, the dangling S12 folds in and binds to Fe10, otherwise the overall pentagonbased prism structure remains intact.The energy difference between the flattened and the prism models indicates the preference of threedimensional network structures when developing computational models for the FeS 2(p) surface.Unexpectedly, the energetic variations were only within 4 kJ mol À1 among the three different coupling schemes involving front/rear (Fe1-Fe7/Fe4-Fe10), top/bottom (Fe1-Fe4/Fe7-Fe10), and left/right (Fe4-Fe7/Fe1-Fe10) pentagons (FC Fe pairs/AFC to each other), as shown in Figure 3C and 3D.

F I G U R E 4
Single point, FeS 2(p) crystal structure-based 8(FeS 2 ) nanoparticles in (A) closed shell, (B) ferromagnetically coupled, and (C) antiferromagnetically coupled electronic structures at the MN15/def2SVPjSMD level of theory.Relative energies were calculated from electronic energies (ΔE SCF ).The numbering scheme of atoms is maintained throughout for all 8(FeS 2 ) models discussed in the manuscript.
reasonably well maintained.The lowest energy initial structures approximate well the lowest energy optimized structures.The modest reordering in stability is due to the structural relaxation of the persulfide ligands, which can result in changes in the Fe coordination numbers, the number of the AFC and FC interactions, and, most importantly, the emergence of the strongly covalent side-on persulfide coordination.The latter provides another criteria for obtaining the lowest energy structure, since increasing the number of Fe-S covalent/ionic interactions further increases the cluster's stability.Three examples for the lowest energy stationary structures of the 8(FeS 2 ) models are shown in Figure6, where the calculated F I G U R E 5 Energy stabilization (ΔE SCF ) relative to the full FC state (S = 16) as a result of changing individual Fe spin states from m s = +2 to m s = À2 (atom-by-atom spin-flip steps) for the 8(FeS 2 ) maquette using the MN15/def2SVPjSMD level of theory.Orange triangles, blue circles, purple squares, and green diamonds correspond to S t = 12, 8, 4, and 0 states with 1, 2, 3, and 4 Fe ions having negative spin densities (m s = À2), respectively.The number of AFC and FC interactions are shown next to the lowest energy spin-coupling schemes.root-mean-square deviations (RMSD) of the structures in Figures 6B and 6C relative to the lowest energy Figure 6A are 0.636 and 0.578 Å, respectively.The RMSD values for only the Fe sites are 0.280 and 0.192 Å, which show that the Fe framework is quite similar; however, the position of the S23-S24 persulfide deviates significantly.The RMSD between the two lowest energy structures (within 1 kJ mol À1 ) is only 0.217 Å, when the spin polarization of the Fe1 and Fe19 sites are switched.Thus, any of these structures are reasonable candidates for a computational nano-reactor for screening chemical reactivity.It is noteworthy that the S-S distances were decreased by 0.02-0.05Å, but the Fe-S distances increased by 0.13-0.18Å relative to the FeS 2(p) crystal

Figure 10 .
Figure 10.As seen for the G2 model, the RMSD deviation between the initial and the optimized structures is 0.443 Å.The RMSD values are mainly due to the rearrangement of the dangling persulfides S95-S96 and S77-S78, forming side-on coordinated ligands.Interestingly, persulfide S41-S42 remains in the 1,2-bridging position between Fe40 and Fe43, surrounded by a S-rich environment, as for the (111) surface in Figure1.In contrast to the surface-exposed Fe-S distances, which deviate more from the FeS 2(p) crystal structure, the Fe-S distances varied by À0.02 -+0.12 Å around the Fe13 bulk center in the 32(FeS 2 ) nanoparticle.An overlay of the initial and optimized structures can be seen in FigureS15and the detailed metric description is provided in TableS8.The increased size of the bulk core in G3 versus G2 models allows for better reproduction of the crystal structure for the Fe-S distances (<0.1 Å) around the diamagnetic hexacoordinated Fe centers.The optimization stabilizes the cluster by 240 kJ mol À1 ,

(
G3) with 28 exposed Fe and 56 S sites displays the structural changes that are expected in going from bulk to surface.Thus, it forms a sufficient model that can realistically capture the interactions for a nanoparticle composed of the Fe 2+ and S 2 2À ions in thermodynamically stable arrangements as in the various surface planes shown in Figure1.In addition, this model will also be ideal for monitoring the structural changes along the proposed FeS 2(p) !Fe 1-x S (p) !FeS (m) interconversion.

FeS 2 ( 6 (
Figure 11 provides a direct comparison of three surface cuts (potentially 3 periodic models) of the most studied surfaces (100) of FeS 2(p) and the analogous chemical environments from a single 32(FeS 2 ) model with bond lengths and interatomic distances.Table2summarizes the diversity of Fe and S environments for each of the models that serve as metal-or ligand-centered reactive sites for redox chemistry, small molecule interaction, and heterometal substitution sites for the surface models in Figure1.Although 4(FeS 2 ) model can be used for exploring physisorption, it is apparent from Table2that the 4FeS 2 model contains only a limited number of environments.On the other hand, the 18(FeS 2 ) model has a comprehensive set of chemical environments and can be used for mapping a broader range of reactivity.Further, the 32(FeS 2 ) model offers the most comprehensive set of Fe and S environments that can be used to probe the reactivity of all surface sites shown in Figure1, in addition to the bulk-to-surface interconversion.In the Introduction, Equation (1) defined the overall thermodynamics for the FeS 2(p) reductive dissolution process.We selected the compositions (n = 4 and 8) to connect the n(FeS 2 ) nanoparticles with diamond-shaped n(FeS) nanoparticles that are organized around a central [2Fe-2S] rhomb.Experimentally, the standard Gibbs free energies of formation was established for FeS 2(p) , Fe 1-x S (p) and FeS (m) minerals as À160 kJ mol À1 , À136 kJ mol À1 , and À 89 kJ mol À1 , respectively.6Supplementing these with the standard free energy of H 2 S (aq) of À28 kJ mol À1 ,49 the thermodynamics for the first step of reductive dissolution process (FeS 2(p) !Fe 1-x S (p) )50 is spontaneous by À4 kJ mol À1 .However, the second step (Fe 1-x S (p) !FeS (m) ) is non-spontaneous (+43 kJ mol À1 ) at standard state, which can be made thermodynamically feasible by considering coupling with acid/base chemistry and elevated temperatures.

= 4 and 8
to be -219(À162) and À468(À356) kJ mol À1 , respectively Extrapolation of these values to experimental conditions of approximately 40 C and 1 μM H 2(aq) concentration for G2 model (n = 8) gives an overall spontaneous reaction for each FeS 2 unit of À11 kJ mol À1 .The apparent disagreement between the mineral phase-based estimate for interconversion (+39 kJ mol À1 ) and calculated values (À11 kJ mol À1 ) highlight the significant energetic difference in thermodynamics between the bulk phase and surface-only reactivity.The large negative numbers clearly show how much more reactive the FeS 2(p) surface is toward reduction than the bulk.

5 |
CONCLUSIONSAtomic-scale assessment of the plausible Fe and S sites at the most dominant FeS 2(p) surface planes revealed the presence of 3-, 4-, and 5-coordinate Fe sites with trigonal pyramidal, seesaw, and tetragonal pyramidal coordination geometries and a wide range of bridging persulfide environments.We developed molecular cluster models or maquettes for investigating a wide variety of localized reactive sites using a single model, as an alternative to employing multiple periodic models for each environment.Specifically, we composed three generations of models with none, one, and six low-spin, diamagnetic Fe sites, as models for pyrite bulk that are surrounded with 8, 17, and 28 surface-exposed, paramagnetic, ferro-and antiferromagnetically coupled, coordinatively unsaturated Fe sites, respectively.Our study F I G U R E 1 1 Highlights of redox reactive sites as a function of Fe coordination environments in the computational nano-reactors (bottom row) and (100) surface planes (upper row) with interatomic distances optimized at the MN15/def2SVPjSMD level of theory.The FeS 2(p) bulk structural parameters are Fe-S(proximal) = 2.26 Å, Fe … S(distal) = 3.45 Å, S-S = 2.16 Å, and Fe … Fe = 3.85 Å. T A B L E 2 Summary of chemical speciation of the surface-exposed Fe (coordination number in parenthesis) and S (number of Fe centers attached in parenthesis) sites for the computational nano-reactors n(FeS 2 ) developed in the given work.
ligand interactions documented the remarkable versatility of the persulfide ligand as a stabilizing force by forming various end-on and side-on bridging coordination environments.In addition, the highspin Fe sites form a complex spin/spin-coupling network that can be energetically weighed by the set of guidelines for obtaining low energy equilibrium structures.The 8(FeS 2 ) G1 model is ideal for screening of redox chemistry and small molecule interactions with the 3 different chemical environments for the Fe and S sites.Extension to 18(FeS 2 ) and 32(FeS 2 ) models allows for probing more diverse Fe and S environments, in addition to allowing for the evaluation of the role of the bulk FeS 2(p) core interconversion among a postulated Fe 1-x S (p) and FeS (m) mineral phases.The largest model with sixteen 2(FeS 2 ) pentagons is our most realistic atomic-scale model for FeS 2(p) framboids.Preliminary thermodynamic calculations employing our models suggest spontaneous reductive dissolution of the FeS 2(p) surface and interconversion into n(FeS) (m) , mackinawite-like nanoparticles using dissolved H 2 .Notwithstanding, the computational nano-reactors are readily employable for studying pyrite oxidation reactions, radical coupling, and the influence of heterometal substitution in modulating reactivity and reaction thermodynamics.