Evaluating structure selection in the hydrothermal growth of FeS2 pyrite and marcasite.

While the ab initio prediction of the properties of solids and their optimization towards new proposed materials is becoming established, little predictive theory exists as to which metastable materials can be made and how, impeding their experimental realization. Here we propose a quasi-thermodynamic framework for predicting the hydrothermal synthetic accessibility of metastable materials and apply this model to understanding the phase selection between the pyrite and marcasite polymorphs of FeS2. We demonstrate that phase selection in this system can be explained by the surface stability of the two phases as a function of ambient pH within nano-size regimes relevant to nucleation. This result suggests that a first-principles understanding of nano-size phase stability in realistic synthesis environments can serve to explain or predict the synthetic accessibility of structural polymorphs, providing a guideline to experimental synthesis via efficient computational materials design.

Nucleation and growth from solution remains one of the most experimentally and geologically important synthesis methods for crystalline solids. Hydrothermal growth, which involves precipitation from a superheated aqueous solution of precursor salts, is a particularly common route for natural mineral formation and synthetic single-crystal growth1. Despite the importance of this method, recipes for the hydrothermal growth of target solid phases remain largely empirical. The absence of predictive approaches to synthesis impedes the realization of novel compounds predicted by highly developed and high-throughput computational materials design approaches2345. In addition, the computational screening of novel materials requires an efficient filter for synthesizability, which has in the past been restricted to the thermodynamic stability of the material of interest and empirical ‘structural similarity' arguments36.

We propose a general framework within which to view hydrothermal growth in a computationally accessible and systematic manner. We hypothesize that with a judicious choice of conditions, a synthesis can be designed such that it leads to the formation of phases that are instantaneous ground states of the relevant thermodynamic potential, when that thermodynamic potential is extended to include less common variables, such as particle size and environment-dependent surface energy. Beyond a certain stage of growth, the material can be kinetically frozen in an intermediate phase, leading to a bulk-metastable synthesis product7. Consequently, metastable phases that are thermodynamic ground states at some intermediate stage of growth are deemed synthetically accessible. This quasi-thermodynamic vision of hydrothermal synthesis provides a baseline of synthesizability that can be readily evaluated computationally to isolate driving forces that favour the formation of a target phase. Knowledge of these driving forces can then guide experimental synthetic efforts, help to computationally select metastable phases that are likely to be synthetically accessible and identify deviations from classical thermodynamic behaviour8, such as in non-classical nucleation9.

We base our study of hydrothermal phase selection on the FeS2 mineral system due to its engineering relevance1011, geologic importance121314 and unresolved structure selection mechanism13. The FeS2 system contains two common phases—pyrite and marcasite—and while hydrothermal recipes for the synthesis of both phases are established12151617, the underlying forces governing phase selection during the growth are not understood13. It is known that marcasite can be grown as the dominant phase below pH=5 (refs 12, 15, 16), despite pyrite being the thermodynamic ground state of bulk FeS2 (ref. 18). However, the mechanism by which pH influences phase selection in FeS2 is unclear as it does not affect the relative stability of bulk pyrite and marcasite13151920.

Here we quantify phase selection during the hydrothermal growth of FeS2 by evaluating the full thermodynamic potential governing the evolution of the system throughout the growth process. The thermodynamics describing particle growth in solution are given by the sum of bulk and surface energy, which scales as , where gb is the volumetric bulk Gibbs free energy of formation, is the particle-averaged surface energy and r is the particle size8. Following the formalism outlined above, we propose that the effect of pH can be understood in terms of nucleation and growth from solution that incorporates this competition between bulk and surface stability21222324. Contrary to the bulk, the surface energies of the two phases vary with pH due to the adsorption of H+and OH− ions252627. By accounting for the adsorption in the evaluation of surface energy, we are able to fully account for the effects of solution chemistry in a theoretical treatment of synthesis, accounting for ‘spectator ions' that influence the growth through the surface of the material, but are not represented in the chemical formula of the bulk product. Thus, we evaluate , the driving force for the formation of the marcasite phase with respect to pyrite, at all stages of growth and as a function of the growth environment. The bulk energy of the growing crystal is determined by the energy of the pure crystal, along with contributions from defect formation, off-stoichiometry and strain. In this work, however, we focus on the growth of pure FeS2 pyrite and marcasite, assuming the bulk energy of both phases to be that of their stoichiometric configuration28. The energy of the solid–liquid interface is governed by a combination of bulk-like bond breaking and off-stoichiometry due to adsorption and segregation. It is convenient to approximate the interface energy by the sum of a solid–solvent interface energy and the free energies of adsorption for solute species within the electrostatic double layer, neglecting in this case segregation from the bulk solid. Thus, we separately compute the free energy of a pristine interface between the stoichiometric solid and solvent, , and the free energy of adsorption of solutes from the solution, , giving us all the information necessary to obtain the free energy of the solid–liquid interface, . While in principle, adsorption-induced segregation can be included29, we do not include this coupling for FeS2 as no significant segregation is to be expected in this compound. By evaluating the thermodynamics relevant to FeS particle growth from first principles, we find that the transition from pyrite to marcasite growth under acidic conditions may be explained by the pH-dependent stability of the surfaces of the two phases, suggesting that the quasi-thermodynamic vision of synthesis described here may serve as a valid and computationally accessible metric of the synthesizability of metastable materials.

Results

Surface thermodynamics

We first evaluate the relative energies of the various crystallographic facets in pyrite and marcasite and their tendency to adsorb OH− and H3O+ions, as given in Table 1 and illustrated in Fig. 1. In pyrite, the (100) and (210) facets are dominant in the vacuum and solvated cases, in line with the common occurrence of these facets in natural cubic and pyritohedral habits of pyrite3031. In marcasite, we find that the spread of surface energies between the different facets is smaller than in pyrite, with the (010), (101), (110) and (111) facets all having low energies, in agreement with their prevalence in natural marcasite crystals31. We then calculate the particle-averaged surface energy of each phase, which gives the total energetic contribution of the surface to the free energy of the solid (Fig. 1b). One way to verify the accuracy of the surface energy curve is through the experimentally measured isoelectric point, which corresponds to a transition from a positively charged (clean or HO+ adsorbed) to a negatively charged (OH− adsorbed) surface. The onset of OH− adsorption onto surface Fe sites in pyrite around pH=1, seen as the point at which the surface energy of pyrite begins to decrease, agrees well with the experimentally measured isoelectric point (IEP) in pyrite at pH=1.4 (refs 25, 26). Similarly, the absence of a maximum in surface energy in marcasite down to pH=0 suggests that some of the marcasite surfaces are always hydroxylated, even at very low pH, which agrees with the lack of an experimentally observed IEP in marcasite within an experimentally accessible pH range32. In both cases, the agreement of the equilibrium particle morphology and adsorption character of negative ions with experimental observations of the IEP indicates that the interface free energies of the low-energy facets in the OH−-adsorption regime are captured reasonably accurately.

Table 1

Surface and adsorption energetics of FeS2 in water.

Phase	Facet
		(J m−2)	(J m−2)	(eV/pH=0, 473 K)	(eV/pH=0, 473 K)
Pyrite	(100)	1.38	1.11	0.300	0.789
	(110)	2.14	1.26	0.757	0.206
	(111)	1.80	1.71	0.072	0.353
	(210)	1.82	1.07	0.899	0.299
Marcasite	(001)	1.70	1.45	0.215	0.210
	(010)	1.54	1.10	0.485	0.884
	(100)	2.12	*	*	*
	(011)	1.75	1.74	0.097	−0.246
	(101)	1.07	0.94	0.684	0.215
	(110)	1.68	1.19	0.471	−0.254
	(111)	1.67	1.21	0.443	−0.097

Calculated surface energies of various facets of pyrite and marcasite in vacuum () and in contact with pure non-dissociated water (), as well as the calculated adsorption energy of H3O+ and OH− at infinite dilution, with respect to the chemical potential of each ion in solution at pH=0, 473 K and a hydrated adsorption site. More precisely, captures the strength of the adsorbate–solid interactions with respect to the free energy of the ion in solution and a hydrated solid surface, but does not include the contribution of adsorbate–adsorbate interactions or configurational entropy on the surface. *, Omitted due to convergence issues on the hydration reference state.

Figure 1

Surface energies of FeS2 pyrite and marcasite.

(a) Equilibrium particle shapes (Wulff shapes) for pyrite and marcasite in vacuum and in solution at pH=0 and pH=7. (b) Surface energies of pyrite and marcasite averaged over the equilibrium Wulff shape across a range of pH levels. The solid lines give the surface energy of pyrite and marcasite per unit area, while the dashed line gives , the effective molar surface energy of pyrite scaled to account for the higher density of pyrite relative to that of marcasite. The density adjustment is necessary for a direct comparison of the effect of surface energy on stability, as it accounts for the fact that the relevant free energy for determining the relative stability of pyrite and marcasite is the molar free energy. Thus, the surface energy of a pyrite particle with an equal mole number to that of a marcasite particle must be scaled down to account for the smaller size of the denser pyrite structure.

From the particle-averaged surface energies given in Fig. 1b, it is clear that while marcasite surfaces are more stable than that of pyrite under highly acidic conditions, a rapid onset of OH− adsorption onto pyrite under more basic conditions stabilizes the pyrite surface. The origin of this transition lies in the variation in OH− adsorption strength among the various facets of pyrite and marcasite. In pyrite, OH− adsorbs onto the (210) facets covering the majority of the Wulff shape, while in marcasite, stabilization due to OH− adsorption is initially limited to the otherwise unstable (110) and (111) facets.

Phase selection during synthesis

The influence of surface energy on phase selection in FeS2 synthesis can be viewed from both a thermodynamic and a kinetic standpoint. Combining the bulk and surface energy of pyrite and marcasite across all sizes and pH levels, we construct the size–pH phase diagram of FeS2, given in Fig. 2. We can immediately see that marcasite is the lowest-energy phase in acid at small particle sizes, giving rise to a driving force for the formation of marcasite under these conditions. At this stage of growth, the system is significantly influenced by nucleation kinetics, which we can analyse within the scope of classical nucleation theory. As shown schematically in Fig. 3, nucleation from solution proceeds over a nucleation barrier, which arises from the energetic penalty of forming a high-surface-area critical nucleus and scales as , where is the average surface energy of the nucleating phase and gb is the volumetric bulk-driving force for precipitation8. The relative rates of pyrite and marcasite nucleation are exponential in the difference between their nucleation barriers, such that even a small decrease in surface energy from pyrite to marcasite can lead to a large excess of marcasite nucleation. Taking the experimentally reported supersaturation for FeS2 hydrothermal growth15, we can immediately see that at the critical nucleus size, marcasite has a lower free energy than pyrite below pH=4, and thus nucleates exponentially faster. This transition to marcasite-dominant nucleation agrees with the experimentally observed onset of marcasite formation between pH=4 and pH=6 (refs 12, 15), and the absence of marcasite when FeS2 is grown under more basic conditions. Thus, we can conclude that marcasite growth in acidic media may be explained by its finite-size thermodynamic stability and preferential nucleation under these conditions.

Figure 2

Thermodynamics of nanoscale FeS2.

The finite-size phase diagram of FeS2 across a range of pH values, illustrating the low-particle-size, low-pH region of thermodynamic stability for marcasite. Note that we report a single critical nucleus size based on the experimentally reported supersaturation15 for both pyrite and marcasite because the difference between the two is negligible.

Figure 3

Relationship between finite-size energetics and nucleation kinetics.

A schematic illustration of nucleation kinetics leading to the formation of metastable marcasite due to a lower kinetic barrier to nucleation and a correspondingly exponentially higher nucleation rate Jm, relative to that of pyrite ( and Jp, respectively)8.

Discussion

The agreement between the experimentally observed stabilization of marcasite in acid and our computational results lends significant credibility to the model of synthesis derived here. Considering that the energy scale of pH is small compared with typical energy scales involved in solid-state chemistry, the agreement of the transition point between pyrite and marcasite within 1–2 pH units is quite remarkable. We speculate that the primary reason underlying this result is the systematic cancellation of error between the chosen adsorbed and reference states. For example, we neglect dispersion interactions in our model due to computational constraints, despite the fact that they certainly play a significant role in determining the behaviour of the real FeS2-water interface. However, the error due to this necessary simplification likely cancels between the adsorbed and reference states of the ions, giving a sufficiently accurate estimate of the behaviour of the system.

One other potential issue in the analysis presented here is the validity of the classical energy decomposition into bulk and surface terms to obtain the energy of the nucleus8. Accounting for the exact dynamics and free energies of the growing nucleus is extraordinarily difficult given any current computational or experimental method, and impractical given the goal of obtaining a scalable ‘synthesizability filter' for computational materials discovery. Instead, the semi-continuum analysis presented here aims to provide a first-order extrapolation of finite-size free energies from the bulk to the size scales relevant to nucleation. Of course, given the small energy scales involved, more detailed studies of the small-scale thermodynamics and nucleation kinetics in this system would help clarify the validity of the approximations made here, as well as identify any non-classical nucleation and growth behaviour that may occur.

We expect that the approach used here to study phase selection between pyrite and marcasite FeS2 can be widely applied to studying polymorphism in other chemical systems synthesized by the hydrothermal method. Indeed, finite-size stability of CaCO3 mediated by Mg uptake from solution has been used to explain the preferential nucleation of metastable aragonite over calcite in present-day oceans24. The general lack of empirical parameters in the derivation of the finite-size phase diagram allows this approach to be extended to chemical spaces with scarce to no experimental data, as is often necessary for computational materials discovery. By constructing the finite-size phase diagram with the inclusion of any arbitrary set of ‘spectator ions' adsorbing on the solid surface, it is possible to identify approximate solution conditions under which there may be a driving force for the formation of a target metastable phase. With this knowledge, one can design syntheses that would allow the system to express the identified driving force, nucleating within the desired region of the phase diagram. Thus, we believe that both the general model proposed here, and the analysis of FeS2 can serve as a useful thermodynamic baseline for predicting phase selection during synthesis and assist the realization of novel materials.

Methods

Thermodynamic model of an aqueous interface

The defining feature of an aqueous interface is the existence of an electrostatic double layer due to the adsorption of charged species to the solid33. However, much of the complexity of the double layer can be neglected when calculating the total energy of the interface. To derive an approximate treatment of this structure, it is helpful to break down the electrostatic double layer into three components: the electronic ‘space charge' region, the chemisorbed region within the Helmholtz plane and the physisorbed ‘diffuse' region outside the Helmholtz plane, as shown schematically in Fig. 4a. The adsorption energy of chemisorbed species is likely large, and thus must be represented accurately, accounting among other features for the charge state of the adsorbate. The space charge region in the solid arises due to the changes in the electronic structure of the solid associated with bond breaking and adsorption from the liquid, and thus will be captured intrinsically by any accurate first-principles model of the surface and the chemisorbed layer. In contrast, the contribution of the diffuse layer to the total energy of the system is not always significant. The electrostatic potential at the Helmholtz plane, experimentally measured as the ξ potential, is typically on the order of 40 mV (refs 25, 26), while the capacitance of the double layer region is on the order of 0.5 F m−2 (refs 25, 26), meaning that the total energy stored in the diffuse layer is on the order of 10−4 J m−2. This quantity is negligible on the scale of total interfacial energies in ceramic–aqueous systems, which are in the range of 1 J m−2 (ref. 34). Thus, a model of the chemisorbed species within the Helmholtz plane, constructed to ensure the correct charge state of the adsorbed species, but neglecting the details of the diffuse layer, yields an accurate estimate of the true interfacial energy, even in the presence of an electrostatic double layer.

Figure 4

Scheme for computing the free energy of an aqueous interface.

(a) The structure of an electrostatic double layer, consisting of a space charge in the solid, a tightly bound chemisorbed layer within the Helmholtz plane and weakly bound physisorbed ions in the diffuse layer. (b) The continuum solvation model provided by VASPsol, which accurately captures the interactions between solute molecules and the solvent, as well as the solvent with itself, but does not accurately represent the interactions between the solvent and an extended solid surface. We introduce a solvation scheme to correct the unreliable solid–solvent interaction with a more physical model calculated from the explicit adsorption of water molecules. (c) To calculate the adsorption energy of charged ions at infinite dilution, we choose a calculation scheme that allows electron transfer between the cationic and anionic species to occur self-consistently. In both the adsorbed and reference states, the ions are sufficiently separated to allow the continuum solvent to partially screen the interactions between them, such that the effect of the electrostatics can be subtracted out analytically as a post-calculation correction. Here we show the adsorption of H3O+ onto a sulfur site on the marcasite (001) surface as an example of this calculation, within the VASPsol continuum solvent model.

To fully model a solid–liquid interface, it is generally necessary to consider the adsorption of all ions present in solution. A more tractable simplification is to consider only the effect of known potential-determining ions, as these ions are by definition those which adsorb strongest and thus determine the structure of the double layer. In the case of FeS2, the potential-determining ions are H3O+ and OH− (refs 25, 26), suggesting that in an aqueous medium, the tightly bound chemisorbed layer consists primarily of these species, in addition to the H2O solvent molecules. Thus, it is only necessary to consider the adsorption of OH−, H2O and H3O+, accounting for all other ions only to the extent that they set the ionic strength and pH of the solution. The energy of adsorbing H2O is equivalent to the energy of solvating the solid surface, and will be addressed in a later section. Further, we assume that at a given pH, OH− and H3O+ will never be adsorbed simultaneously. On the basis of this set of assumptions, the adsorption energy of the two ions can be written based on the minimization of the free energy of adsorption, given by the enthalpy and entropy of adsorption with respect to the number of adsorbates :

We approximate the enthalpy of adsorption by accounting for adsorbate–solid, adsorbate–solvent and adsorbate–adsorbate interactions. The adsorbate–solid and adsorbate–solvent interactions are captured by the enthalpy of adsorption at infinite dilution , where Eads and Eref are the density functional theory (DFT) energies of the ion adsorbed onto the solid and in solution respectively, shown schematically in Fig. 4c. Note that this energy of adsorption includes the energy of desorbing a water molecule, as the adsorption process is competitive with the pure solvent. We approximate adsorbate–adsorbate interactions with the Debye–Huckel model of screened electrostatics in an electrolytic medium, given here by Vel. Thus, we can write the enthalpy of adsorption for OH− and H3O+ as:

where r are the positions of the adsorbed ions on the solid surface, λ is the Debye screening length of the solution and ɛ is the dielectric constant of the solution near the interface. In our model, we use a Debye screening length of λ=1.0 nm, based on an average over screening lengths for reported synthesis recipes for pyrite and marcasite15. As we are considering the two-dimensional electrostatic interactions between adsorbates, screened by adsorbed water molecules, we set the dielectric constant ɛ=12, estimated from reported experimental and computed values of the dielectric constant of interfacial water in similar systems3536. Finally, we set our temperature to 473 K in accordance with the experimental conditions commonly reported for FeS2 hydrothermal growth1516.

To obtain the entropy of adsorption, we consider the entropy of the adsorbed and solution states of the ion, sads and ssoln, respectively. The entropy of the adsorbed ion is well approximated by the configurational entropy over adsorption sites. The entropy of the ion in solution is given by the configurational entropy over the translational degrees of freedom of the ion in solution, which, assuming that OH−, H2O and H3O+ all have approximately the same volume, is given by kBlog[x], where x is the mole fraction of the ion of interest. We then relate the entropy of H3O+ to pH by treating pH as an activity with respect to a 1 M solution of H3O+ at standard state and assuming that the solution behaves ideally, which yields:

where the Mw is the molarity of water and T0=300 K is temperature in the reference state. Following the same assumptions, as well as the fact that OH−, H2O and H3O+ are in equilibrium, we derive the entropy of OH− in solution in terms of the calculated formation enthalpy of H3O+ and OH− from 2H2O, which we denote :

A detailed derivation of these results is given in Supplementary Note 1. Combining the solution references with the configurational entropy of the adsorbed state, we have the entropy of adsorption for both H3O+ and OH− in terms of pH:

We have thus obtained a thermodynamic picture of ion adsorption that is efficiently computable from first principles and captures the primary trends we could expect to see at the solid–liquid interface as a function of pH. A similar analysis can be readily performed for other dissolved ions, based on their computed solubility product Ksp, generalizing this approach to a solid–aqueous interface with any ideal or near-ideal aqueous solution.

Computational implementation

On the basis of the thermodynamic framework derived above, it is clear that to obtain a full quasi-thermodynamic picture of the hydrothermal growth of FeS2 pyrite and marcasite, only a few density functional theory calculations are necessary. First, we must calculate the bulk energy and structure of pyrite and marcasite. Then, for each low-energy crystallographic facet of each phase we must obtain the interfacial energy between the FeS2 solid and water, or equivalently, the solvation energy of each crystal facet, . Finally, for each solvated facet, we must calculate the enthalpy of adsorbing dilute H3O+ and OH− ions onto all likely adsorption sites, . An example calculation, illustrating the thermodynamic formalism can be found in Supplementary Note 2.

Computational details

All calculations were done using the Vienna Ab-Initio Simulation Package (VASP)3738 implementation of DFT, using projector augmented wave (PAW) pseudopotentials3940 with a plane-wave basis set using an energy cutoff of 520 eV. Consistently with previously reported results, we find that the PBEsol exchange-correlation functional41 provides an accurate and computationally efficient model of FeS2 (refs 42, 43), correctly stabilizing pyrite over marcasite as the ground state of the system in agreement with experiment1244 and higher-order functionals, although not quite reaching the experimentally measured transition enthalpy between the two phases18 (Supplementary Table 1). To ensure consistency between the high-symmetry bulk calculations and low-symmetry adsorption calculations, we remove all symmetry restrictions from the calculation, giving the system identical relaxation degrees of freedom across all calculations. Finally, for bulk calculations, we choose a Γ-centred k-point mesh (6 × 6 × 6 for pyrite and 6 × 6 × 8 for marcasite) based on previously optimized calculation parameters in similar systems45.

In our surface calculations, we consider the symmetrically distinct low-index facets of pyrite and marcasite previously reported to be significant. Specifically, in pyrite, we consider the (100), (110), (111) and (210) facets4647, while in marcasite, we consider the (100), (010), (001), (110), (101), (011) and (111) facets31, defined with respect to the unit cells given in Supplementary Table 1. To generate the surface structures, we choose surface terminations that minimize the number and strength of bonds broken, evaluated based on the integral of charge density associated with each bond in question, and are maximally non-polar, following the Tasker surface stability criterion48. In the case where several surface terminations satisfy these criteria, we consider all such terminations. The resulting most stable (under solvated conditions) surface structures are shown in Supplementary Fig. 1. While there is limited experimental data available to verify the accuracy of this approach, in the case of the well-characterized pyrite (100) surface, our approach leads to a surface structure consistent with that derived from low-energy electron diffraction (LEED) characterization49. Finally, we neglect the contribution of the solid to the solid–liquid interface entropy as it is known to be negligibly small in similar ceramic systems50.

Solvation model

To account for solvation, we rely on the VASPsol continuum solvation model51 to avoid the computationally prohibitive sampling of explicit solvent configurations. The VASPsol model serves two important purposes–it reproduces the mean-field interactions between ions and bulk solvent, and provides a dielectric medium that screens electrostatic interactions between charged adsorbates, their counterions and their periodic images. However, while the VASPsol model is known to accurately reproduce the energy of solvating isolated molecules51, its performance with respect to the solvation of solid surfaces is uncertain.

To correct any unphysical interactions between the VASPsol continuum solvent and the solid slab, we introduce a solvation correction scheme. We assume that the VASPsol model accurately reproduces all solvent–solvent and solvent–ion interactions, but fails to capture solvent–solid interactions, as shown schematically in Fig. 4b. To correct this error, we first remove the energy associated with the interaction of the continuum solvent and the solid by subtracting out the difference between the energy of the clean surface (solid ‘slab') in contact with vacuum and in contact with the continuum solvent , which we will refer to as . We then add back the interactions between the solvent and the slab by explicitly calculating the energy of adsorbing isolated water molecules within the continuum solvent, , for each adsorption site on the solid, where is the experimentally measured entropy of bulk water. Note that we neglect the entropy of the adsorbed water as we assume that the interfacial water layer is relatively constrained and ice-like, significantly reducing its entropy relative to that of the bulk solution52. Having obtained this shift for each adsorption site, we can correct any calculation done with only the continuum solvent to capture the solid–solvent interactions potentially misrepresented by VASPsol.

For example, to calculate the energy Einterfac of a surface with Nsites identical adsorption sites, of which Nads are occupied by some adsorbing ions and Nsites−Nads are filled by water, we calculate the energy of a periodic slab with explicit adsorbed ions (but not water molecules) within VASPsol to get . We then apply the solvation correction to get the true interface energy:

In the case where there are distinct adsorption sites, the energy of solvation becomes site specific. For FeS2, we assume that the site-specific solvation energy is determined by the local chemistry, giving separate solvation energies for Fe and S sites on the surface. As it is impossible to adsorb H2O simultaneously to adjacent Fe and S sites due to steric constraints, we take the lower energy of the Fe and S adsorption sites for each facet as the facet-specific solvation energy , and the density of these sites as the number of adsorption sites Nsites. Combining these terms, we obtain the solid–solvent interfacial energy term, with Nads=0:

where simplifies to in the case of zero adsorbed ions.

Charged adsorption

To obtain a full picture of interfacial stability across various solution conditions (here, pH levels), it is necessary to calculate the enthalpy of adsorption of all relevant ions (here, OH− and H3O+) at infinite dilution, as discussed in the thermodynamic derivation earlier. To do so, it is necessary to ensure that in both the adsorbed and reference state, the charge state of the adsorbing ion is physical. Under periodic boundary conditions imposed by plane-wave DFT, the charge state can be set explicitly by removing a number of electrons from the system, and compensating the resulting charge with a homogeneous background. A more physical model is to include a counter charge in the system and allow charge transfer from the cationic to the anionic species to occur self-consistently. However, in this case, it is important to take care to ensure that no unphysical electrostatic interactions between the anion, cation and their periodic images contribute to the adsorption energy.

One approach to ensure that unphysical electrostatic interactions do not contribute to the adsorption energy is to construct a supercell large enough, such that the electrostatic interactions between ions decay to a negligible level. A more computationally tractable approach is to choose a reference state, such that the electrostatic interactions either cancel out between the adsorbed and reference states, or can be subtracted out analytically, giving an accurate enthalpy of adsorption at infinite dilution . One such choice of reference state is given in Fig. 4c. In this set-up, the nearest neighbour cation-cation and anion-anion image interactions cancel out between the adsorbed and reference state, leaving only the cation-anion interactions and ion–solid interactions. In both cases, the continuum solvent medium provided by VASPsol ensures rapid decay of the electrostatic interactions as a function of distance, such that even at a 5 Å minimum separation (with a VASPsol dielectric constant ɛ=80 for an aqueous system), electrostatic interactions between the ions are small enough that they can be subtracted out from the total energy as a post-calculation correction. If we assume that the ion–solid interaction in the reference state is small, which is reasonable considering that the solid is not charged and the separation between the ion and solid is over 10 Å, the only remaining interaction is the one we are interested in—the ion–solid interaction in the adsorbed state. Note that in this case, the ion charge state is self-consistently set to the correct value in both the adsorbed and reference configurations, as can be confirmed by a Bader charge analysis53. Relaxed adsorption geometries for H3O+ and OH− derived using this approach for all facets of pyrite and marcasite can be found in Supplementary Figs 2 and 3, respectively.

Data availability

Complete computational data to support the findings of this study is available from the authors on reasonable request.

Additional information

How to cite this article: Kitchaev, D. A. & Ceder, G. Evaluating structure selection in the hydrothermal growth of FeS2 pyrite and marcasite. Nat. Commun. 7, 13799 doi: 10.1038/ncomms13799 (2016).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.