A universal chemical potential for sulfur vapours

The equilibrium thermochemistry of mixed vapour-phase sulfur allotropes is computed and solved to obtain a single chemical potential function.


Introduction
Sulfur is an abundant resource exploited by industry on a scale of tens of millions of tonnes per year. 1 While it may be found in its elemental form, the primary industrial source is hydrogen sulfide, a byproduct of the oil and gas industry.The vast majority of industrial sulfur is converted to sulfuric acid or sulfur dioxide before further use; this may explain the surprising shortage of data in the thermochemical literature regarding the vapour phase of elemental sulfur.
Historically, the thermochemistry of sulfur has been studied experimentally and has been understood to be associated with a variable composition for over a century; Lewis and Randall remarked in 1914 that "no other element is known to occur in as many different forms as sulfur" while studying the free energy of a number of these forms. 2(Carbon now has a higher number of known allotropes but the majority of these are not naturallyoccuring.)However, contemporary reference data for sulfur still does not present a complete picture; the NIST-JANAF Thermochemical Tables (1998) give thermochemical data for two solid phases, one liquid phase, the ions S + and S -and eight gas allotropes S 1-8 . 3Of these, only S 2 and S 8 are from spectroscopic data.The allotropes S 3-7 are assumed to exist and are assigned energies following an interpolation scheme suggested by Rau et  al. (1966), which also makes use of experimental data for S 6 . 4hat paper rules out the significant presence of tautomers, finding little evidence of a tautomer contribution and assuming that they have relatively high energy.The authors generally reserve speculation on the actual structures of the components of their equilibrium model.
In recent years considerable attention has turned to metal chalcogenides; II-VI semiconductors such as ZnS, CdS, PbS are widely studied in many contexts. 5Copper indium gallium selenides (CIGS) and cadmium telluride (CdTe) are used as the basis for "second-generation" thin-film photovoltaic devices, and have seen a dramatic rise in production.Cu 2 ZnSn(S,Se) 4 (CZTS) and Cu 2 SnS 3 (CTS) devices have so far struggled to match these materials in terms of energy conversion efficiencies, but hold significant long-term promise due to their use of highly abundant elements; such availability is a prerequisite for terawattscale photovoltaics. 61][12] It is undermined however by the model derived by Rau et al., which predicts that no one component makes more than 50% of the gas mixture at temperatures between 800-1100K. 4ss spectrometry at a relatively mild 105 • C has observed a series of charged clusters with the form (S 8n ) + . 13In the mid 1980s, a number of cyclic allotropes had been identified by crystallisation and X-ray diffraction, but this only covered the range n = 6 -20. 14An ab initio study was carried out for S 2 through to S 13 in an early application of the Car-Parrinello simulated annealing method. 15Energies were calculated using density-functional theory with the local density approximation (LDA).While limited by the inherent difficulties in exploring the entire potential energy surface of the atomic positions, this thorough study generated 21 allotropes, finding a local maximum in the atomisation energy at n = 8.A later (1990) paper used coupled-cluster electronic structure calculations to study the proposed tautomers of S 4 in depth, concluding that the planar structure with C 2v symmetry is lowest in energy, with a trans (C 2h ) structure also visible in experimental spectra; a more recent ab initio study reached similar conclusions regarding stability while challenging the spectroscopic assignment of the phases. 16,17The C 2v structure was ruled out in the simulated annealing study with LDA, although the authors noted the experimental evidence for its existence. 15A 2003 review by Steudel et al. 18 collects more recent data, including both experimental and theoretical studies of vapour-phase allotropes; this review notes the weakness of the widespread assumption that each size is represented by a single species. 18The work compares several sets of enthalpies relative to S 8 that have been obtained experimentally; variability is high for the smaller allotropes while there is fairly good agreement for the larger allotropes.Studies are generally carried out at a single temperature, such that the temperature and pressure dependence of the thermochemistry must be derived from statistical mechanics and analysis of vibrational information.
In this study, we develop a set of structures for S 2 -S 8 , compute their Gibbs free energy from first-principles and with empirical corrections, and solve the temperature-dependent chemical potential to describe the gaseous mixture.The potential function will be important for quantitative investigations of defect formation and phase stability in metal sulfide materials.

Density functional theory
Energies and forces of arbitrary clusters of sulfur atoms were computed within Kohn-Sham density-functional theory (DFT). 19 range of exchange-correlation functionals were used in this work: PBE is a popular and elegant implentation of the Generalised Gradient Approximation (GGA) and PBEsol restores a periodic exchange contribution leading to improved performance for solids; 20,21 B3LYP * is a widely-used "hybrid" functional which combines pre-existing gradient corrections with "exact" Hartree-Fock exchange; 23 PBE0 is applies similar principles to the parameter-free PBE functional.24 (While PBE is generally preferred to PBEsol for molecular calculations, PBEsol was included in this study for its compatibility with other all-electron work using this functional.) Calculations for the evolutionary algorithm search used the Vienna Ab Initio Simulations Package (VASP) with the PBE exchange-correlation functional and a plane-wave basis set with a 500 eV energy cutoff. 25,26As calculations in VASP employ a periodic boundary condition, orthorhombic bounding boxes were employed with 10Âȃ of vacuum between each molecule and its periodic images.Electronic structure iteration used only the Γpoint of this large cell.
Further calculations used the Fritz Haber Institute ab initio molecular simulations package (FHI-aims) to carry out all-electron DFT calculations with numerically-tabulated basis sets. 27,28All calculations were open-shell with S 2 adopting its low-energy triplet spin configuration.The recommended "tight" basis set was employed for initial relaxation and study with PBEsol, which extends the minimal set of occupied orbitals with 6 additional functions.This was extended further to the full "tier 2" set of 9 additional functions for calculations with the LDA, PBE0, and B3LYP functionals.

Global structure search
0][31] At this stage, molecules with n > 8 were disregarded, as experimental results anticipate high-and low-temperature limits dominated by S 2 and S 8 , respectively.Clusters were generated for S 3-7 , and refined with an evolutionary algorithm to minimise the ground-state energy until a number of seemingly distinct clusters were identified by inspection.The atomic positions of these clusters were then optimised in FHI-aims calculations with PBEsol, using the BFGS algorithm to minimise the atomic forces to less than 10 -4 eV Å -1 and converge energy to within 10 -6 eV.Point groups were assigned to the structures using Materials Studio version 6.0, a proprietary package developed by Accelrys.

Vibrational frequencies
Vibrational frequencies were calculated within the harmonic approximation by making finite displacements to each atomic position to obtain the local potential wells, and diagonalising the resulting dynamical matrix to obtain the normal modes and their frequencies.This is implemented as a script and diagonalisation routine provided with FHI-aims.
Improved vibrational frequencies may be obtained by applying an empirically-derived scale factor to the vibrational eigenvalues computed using DFT; collections of such scale factors have been published for large test-sets of molecules. 32,33The use of these factors is somewhat problematic when creating a systematic, transferable set of data but offers an opportunity to create the most realistic thermochemical model possible.Given that the calculations in this work involve a more limited subset of atomic interactions, we choose to fit a scaling factor to the experimentally-reported frequencies of S 8 and S 2 .

Thermochemistry of individual gas species.
Thermochemical properties were calculated within the ideal gas, rigid-rotor and harmonic vibration approximations.A set of textbook equations forms the chemical potential µ for a nonlinear molecule from the ground-state electronic energy E 0 given a set of vibrational energies ε, the rotational constant σ , moment of inertia where T 0 These were applied as implemented in the Atomic Simulation Environment (ASE) Python package. 34(Note that the expressions for monatomic and linear molecules are slightly different.)The rotational constants σ were assigned from the point groups.

Reference energies.
A number of ab initio methods have been applied.In order to compare the energies, a reference point is needed.Conventionally the enthalpy of the ground state is zero; however, in this case the ground state phase α-sulfur is relatively expensive to compute.We therefore use the experimental sublimation enthalpy to obtain a reference from the calculated enthalpy of S 8 : The preferred experimental value for ∆H sub is 100.416/8= 12.552 kJ mol −1 , from experiments at 298K. 3 Note that the physical system does not in fact sublime at high temperatures, but passes through a molten phase.Nonetheless, it is more practical (and perfectly valid) to retain α-S as the reference state over the whole temperature range studied.

Equilibrium modelling.
The following derivation closely follows the approach and notation of Ref. 35, which describes a generalised "non-stoichiometric method" for solving chemical equilibria.This approach is wellestablished and based on key work in Refs.36-38.
We attempt to minimise the Gibbs free energy subject to the mass balance constraint where N is the number of unique species i with stoichiometric coefficient a i ; n is the quantity of species i and b is the total number of sulfur atoms.The classic approach for a constrained optimisation is the method of Lagrange multipliers.The Lagrangian is formed and differentiated to form a set of equations definining the equilibrium state. and The species chemical potential µ i calculated as in Section 2.4.1 is a function of both temperature and the partial pressure p i = P n i n t where P is the total pressure and the total quantity n t = ∑ N i n i .The temperature dependence is complex and we are willing to solve the equilibrium at each temperature of interest, so we form a temperature-dependent standard free energy at a reference pressure P From here we drop the parenthetical indication that µ • i is a function of temperature, and define the unit of pressure as the reference pressure, such that P • = 1.Substituting ( 16) into (12), we obtain and summing over i The only unknown variable in this expression is λ ; rearranging slightly we form a polynomial which is suitable for solving by standard numerical methods.The method employed in this work is the Levenberg-Marquardt least-squares algorithm, as implemented in Scipy. 39,40 To recover the composition n, we rearrange (18): and substitute into the second equilibrium condition (10) to obtain combining ( 21) and ( 22) we eliminate n t and clean up the notation by denoting exp Finally, to obtain the chemical potential of the mixture we note from (12) that µ i a i = λ for all i.Therefore the normalised chemical potential of sulfur vapour on an atom basis.(A mathematical derivation is given in Appendix A.)

Sulfur allotropes
A variety of candidate structures were generated in the evolutionary algorithm study with the PBE functional.The low-energy candidates following geometry optimisation are discussed in this section.

S 2 .
Diatomic sulfur has the point group D ∞h , in common with other homonuclear diatomics.The atoms were initially set 2 Å apart, and relaxed to a bond length of 1.91 Å.Studies with other functionals were relaxed either from this distance or from 2 Å.The resulting bond lengths are given in Table 1. ). 41 We have also considered the simple triangular allotrope, which is ∼ 0.5 eV higher in ground-state energy.

S 4 .
A range of branched and cyclic structures were generated in the evolutionary algorithm.The structures included in the equilibrium modelling are shown in Fig. 1.The lowest-energy structure identified was the 'eclipsed' C 2v chain; this is in agreement with the high-level theoretical studies in Ref. 16,17.These studies identified a 'trans' C 2h structure as being likely to exist; there is some spectroscopic evidence for the viability of this isomer as well as a branched chain, but we were not able to reproduce stable structures corresponding to these allotropes through geometry optimisation. 42,43Various cyclic and tetrahedral candidate structures yielded a relatively flat puckered ring with D 2d symmetry.

S 5 .
Although a wide range of branched and chain structures were generated, the main candidate is the 5-membered ring with C s symmetry.

S 6 .
In addition to a cyclic C 2v allotrope, relatively low-energy branched and chain variations were identified.Of considerable interest is also a structure which may be viewed as a stack of two S 3 cycles, or alternatively as a cluster of S 2 diatoms.This appears to be the D 3h "prism" structure identified by by Wong et al. 44 ; the characteristic S-S bond lengths from that study were 190.1 and 276.2 pm, while the corresponding average distances from optimisation with the same hybrid XC functional (B3LYP) in this work were 189.0 and 275.7 pm.It is worth stressing that no explicit dispersion terms were included in any of the electronic structure calculations.The evolutionary algorithm results rapidly provided the same C s cyclic structure as that obtained by energy minimisation from a regular polygon.A branched structure, generated early in the progress of the algorithm, was also selected as an interesting alternative to include.This was about 1 eV lower in energy than the other candidates at that stage.Geometry optimisation by force relaxation yielded a compact structure, also with C s (mirror-plane) symmetry.

S 8 .
No evolutionary algorithm study was applied for S 8 , as its ring structure is quite well-known.The initial geometry was extracted from the crystal structure for the condensed α-S phase used in a previous study, 45 and relaxed to form an isolated D 4d ring.

Ground-state energies.
An inspection of the ground-state energies from DFT reveals a trend of smoothly decreasing energy per atom with cluster size for the minimum-energy configuration at each size (Fig. 2).The variation within the clusters included at each size is of the order 10 kJ mol −1 atom −1 , which is comparable to the energy difference between neigbouring cluster sizes.

Vibrational properties
Vibrational frequencies were calculated for all of the allotropes listed in section 3.1; frequencies for S 2 and S 8 are listed in Table 2.

Empirical corrections.
Empirical scale factors were determined by fitting the frequencies to the experimental spectrum for S 8 .Note that frequencies are linearly proportional to their corresponding zero-point energies E ZPE = 1 2 hν and hence this may also be seen as fitting to zeropoint energy on a per-mode basis.The factors were calculated for each functional (Table 3); scaling the frequencies from PBE0 by 96% was found to give the best overall fit, and is employed here  as the reference "empirically-corrected" method.The resulting set of frequencies is illustrated in Fig. Âȃ3 alongside the uncorrected and experimental values.Using this scale factor also gives good agreement (< 4 cm -1 error) with the stretching frequency of S 2 , which was not used in the fit.(Table 2) Least-squares fitting was carried out with the Levenberg-Marquardt algorithm as implemented in Scipy. 39,40

Equilibrium model
Equilibrium compositions and free energies were computed as a function of temperature and pressure for all the data sets computed (Fig. 4).There is significant disagreement between the predictions of the local exchange-correlation functionals LDA and PBEsol and the predicted composition from the hybrid functional PBE0, both before and after frequency scaling.While the "lowerlevel" calculations predict a diverse mixture of phases, hybrid DFT strongly supports the dominance of S 8 and S 2 , at low and high temperatures respectively.In all cases, this simplicity is strongest at low total pressure.The other phases which are present in any  significant quantity are the cyclic allotropes where N = 4-7, in the range 600-1000 K.
The corresponding free energies are also plotted in Figure 5; we note that agreement between the methods is much stronger at low temperatures where the mixture is dominated by larger molecules.This may be an artefact of aligning the free energies of the S 8 atoms; divergence in the energies of the smaller molecules leads to the disagreement at high temperatures.The other trend of note is the presence of a sharp bend in the µ − T curve, particularly at low pressure, corresponding to the presence of S 2 molecules.The point of onset depends on the data source, but the curve for PBE0 with empirical corrections closely tracks the minimum of the two curves from reference data.This represents a challenge to the formation of a simple parameterised model function, as it suggests the presence of a spike in the second derivative.Popular parameterisations of thermochemical properties, such as those in the NIST "WebBook", employ multiple temperature regions.This is usually viewed as a limitation, as it introduces non-physical discontinuities; with care, they could be aligned to an apparently physical discontinuity in the function.Taking the PBE0 results with empirical corrections as our preferred model, the free energy of the mixture is plotted with the chemical potentials of its component species on an atomic basis (Fig. 6).
The depression in free energy due to mixing of allotropes and  presence of minor components can be quantified by subtracting the chemical potential of the mixture from the minimum of the chemical potentials of the majority components S 2 and S 8 .The resulting plot (Fig. 8) shows that this has an impact ranging from around 1-4 kJ mol −1 , depending on the pressure.This is illustrated as a contour plot in Fig. Âȃ7; within each unshaded region a single-phase model is adequate to within 1 kJ mol −1 S atoms.

Parameterisation
For convenience, a parameterised fit has been generated for the chemical potential of S over the T, P range 400-1500K, 10 0 -10 7 Pa, incorporating an error function "switch" between S 2 and S 8 dominated regions and a Gaussian correction for the free energy depression where there is substantial mixing of phases.In eV per S atom, for T in K, the form of the parameterisation is It is noted that this parameterisation contains many fitting parameters; however, given its physically-motivated form the resulting function is smooth and well-behaved over the region studied, while the fits to µ S 2 , µ S 8 and T tr have some value in their own right.The fitting error is plotted in Fig. 9, and while somewhat irregular remains below 1 kJ mol −1 .

Conclusions
The chemical potential of sulfur vapours has been studied by solving the thermodynamic equilibrium of 13 gas-phase allotropes, including the dominant components S 2 and S 8 .Thermochemical data was obtained from first-principles calculations and corrected with an empirical scaling factor for the vibrational frequencies.The transition between these dominating phases is highly pressure-dependent, and the free energy is further depressed at the transition temperature by the presence of additional phases, especially at elevated pressures.Selection of an inappropriate gas phase can lead to errors of the order 50 kJ mol −1 atoms, while the minor phases contribute free energy of the order 1 kJ mol −1 atoms.The resulting chemical potential data is made available through tabulated data, a parameterised model with error of the order 0.5 kJ mol −1 atoms and through open-source code; the reference energy is compatible with the NIST-Janaf thermochemical tables for the solid α-sulfur phase. 3This phase is frequently used as a reference state for thermodynamic studies of defects and sta- bility in metal chalcogenides; the application of this gas-phase potential may allow such studies to examine a wide range of reactions involving sulfur vapours, taking into account the equilibrium within the vapour phase.The selection of appropriate chemical potentials is also critical for the development and interpretation of phase diagrams.

Data Access Statement
The reference implementation of this model, complete with Python 2.7 code to generate all the plots in this paper as well as tabulated data in the form of Table 4, is available online at https://github.com/WMD-Bath/sulfur-modeland a snapshot of the code at the point of submission of this article is hosted by Zenodo and available with the DOI: 10.5281/zenodo.28536.In addition, full tables are provided with this paper in the ESI † for the composition, enthalpy and chemical potential from the calculations with PBE0 and empirical corrections; one set of enthalpy and chemical potential data follows Table 4 and uses the enthalpy of α-S as a reference energy (for use with other tabulated data) while the other employs the ground state of S 8 as a reference energy (for use with first-principles calculations.)The code and its dependencies are Free Software, using a range of licenses.Input and output files from DFT calculations with FHIaims have been deposited with Figshare and are available with the DOI: 10.6084/m9.figshare.1513736.A set of data generated during the evolutionary search, consisting of candidate structures and the DFT energies used to rank them, has been deposited with Figshare and is available with the DOI: 10.6084/m9.figshare.1513833.(Notation the same as Sec.2.4.3).From (12), µ i = a i λ and hence ĜS (T, P) = ∑ N i=1 a i λ Φ i ∑ N j=1 a j Φ j = λ (29)

Fig. 1
Fig. 1 Predicted low-energy sulfur clusters with symmetry assignment

Frequency / cm − 1 Fig. 3
Fig.3Vibrational frequencies of S 8 calculated with various DFT functionals, compared with recommended experimental values.3

Fig. 4 Fig. 5 Fig. 6
Fig.4 Compositions of modelled S x mixtures over range of equilibrium temperatures and pressures.Results are presented for density functional theory with one local (LDA), one semi-local (PBEsol) and one non-local exchange-correlation functional with empirical corrections.Composition is given in units of atom fraction.It is expected that the most accurate results are obtained using PBE0 with scaled frequencies.

Fig. 7
Fig.7Temperature-pressure map of approximations to free energy of mixture.At dashed line1 2 µ S 2 = 1 8 µ S 8 ; in shaded region the error in chemical potential µ associated with assuming a single phase S 2 or S 8 exceeds 1 kJ/mol S atoms; in unshaded regions the corresponding single-phase free energy is close to the energy of the mixture.

Fig. 8 Fig. 9
Fig. 8 Depression in chemical potential of sulfur vapour µ S due to mixing and presence of minor allotropes.∆µ mixture = µ S − min µ S 2 2 , µ S 8 8 Centre for Sustainable Chemical Technologies and Dept. of Chemistry, University of Bath, Claverton Down, Bath BA2 7AY, UK b Global E 3 Institute and Department of Materials Science and Engineering, Yonsei Uni-Additional data and code available in external repositories with DOIs: 10.5281/zenodo.28536;10.6084/m9.figshare.151373;10.6084/m9.figshare.1513833.See Data Access Statement for more information.

Table 1
3alculated and experimental bond length r in S 2 .Experimental value is NIST/JANAF-recommended distance.3Theevolutionaryalgorithmprocess eliminated all but a C 2v nonlinear chain for S 3 .This corresponds to "thiozone", which has a well-characterised structure by rotational spectroscopy (bond length 1.917(1) Å and angle 117.36(6)• ; the values from optimisation with PBE0 in this study are 1.901 Å and 118.2 •

Table 2
Calculated and experimental vibrational frequencies for S 2 and S 8 .
3round-state energies from DFT of clusters included in study.Energies are relative to the energy for S 8 with each functional, and normalised to the number of atoms.A point is also included from reference data3; this is derived from the enthalpies of formation at zero temperature, based on spectroscopic observations and equilibrium studies.While the energies from different exchange-correlation functionals diverge across the series, the S 2 energy from PBE0 calculations agrees closely with this reference data.

Table 3
Optimal scale factors for exchange-correlation functionals, fitting to ground-state frequencies of S 83 .Standard deviations s for the least-squares fit are given over the set of frequencies in units of frequency and their corresponding zero-point energies per sulfur atom.

Table 4
Gibbs free energy of S vapours, tabulated from calculations with PBE0 and empirical corrections, with reference state (H=0) α-sulfur at 298.15K.Energies in kJ mol -1 , column headers in log 10 (pressure/Pa).Tables are provided with more values and greater decimal precision in the supplementary information.