A universal chemical potential for sulfur vapours

Abstract

The unusual chemistry of sulfur is illustrated by the tendency for catenation. Sulfur forms a range of open and closed S_n species in the gas phase, which has led to speculation on the composition of sulfur vapours as a function of temperature and pressure for over a century. Unlike elemental gases such as O₂ and N₂, there is no widely accepted thermodynamic potential for sulfur. Here we combine a first-principles global structure search for the low energy clusters from S₂ to S₈ with a thermodynamic model for the mixed-allotrope system, including the Gibbs free energy for all gas-phase sulfur on an atomic basis. A strongly pressure-dependent transition from a mixture dominant in S₂ to S₈ is identified. A universal chemical potential function, μ_S(T,P), is proposed with wide utility in modelling sulfurisation processes including the formation and annealing of metal chalcogenide semiconductors.

---

1 Introduction

Sulfur is an abundant resource exploited by industry on a scale of tens of millions of tonnes per year.¹ 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.² (Carbon now has a higher number of known allotropes but the majority of these are not naturally-occurring.) 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.³ Of these, only S₂ and S₈ 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₆.⁴ That 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.⁵ Copper 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₂ZnSn(S,Se)₄ (CZTS) and Cu₂SnS₃ (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 terawatt-scale photovoltaics.⁶ As such, thin-film processing in sulfur atmospheres is of considerable interest, as the inherent safety of industrial processing may be improved by eliminating the use of toxic H₂S. In addition to chalcogen annealing, which is used to increase grain size, substitute other elements or directly form chalcogenides from elements, high-quality single-crystal samples may be produced using chemical vapour transport of elemental chalogens.^7–9 Previous work on the thermodynamics of such processing has tended to assume that sulfur adopts one particular gaseous allotrope (either S₂ or S₈), but the validity of this assumption has not been explored in depth.^10–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–1100 K.⁴

Mass spectrometry at a relatively mild 105 °C has observed a series of charged clusters with the form (S_8n)⁺.¹³ In 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.¹⁴ An ab initio study was carried out for S₂ through to S₁₃ in an early application of the Car–Parrinello simulated annealing method.¹⁵ Energies 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₄ 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,17 The C_2v structure was ruled out in the simulated annealing study with LDA, although the authors noted the experimental evidence for its existence.¹⁵ A 2003 review by Steudel et al.¹⁸ 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.¹⁸ The work compares several sets of enthalpies relative to S₈ 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₂–S₈, 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.

2 Methods

2.1 Density functional theory

Energies and forces of arbitrary clusters of sulfur atoms were computed within Kohn–Sham density-functional theory (DFT).¹⁹ A range of exchange–correlation functionals were used in this work: PBE is a popular and elegant implementation 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;²³ PBE0 is applies similar principles to the parameter-free PBE functional.²⁴ (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,26 As 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,28 All calculations were open-shell with S₂ 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.

2.2 Global structure search

Global structure optimisation was carried out with the USPEX package, which was originally developed for crystalline systems and has been adapted for use with clusters.^29–31 At this stage, molecules with n > 8 were disregarded, as experimental results anticipate high- and low-temperature limits dominated by S₂ and S₈, 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⁻⁴ eV Å⁻¹ and converge energy to within 10⁻⁶ eV. Point groups were assigned to the structures using Materials Studio version 6.0, a proprietary package developed by Accelrys.

2.3 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,33 The 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₈ and S₂.

2.4 Thermochemistry

2.4.1 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₀ given a set of vibrational energies ε, the rotational constant σ, moment of inertia Iwhere

C_v = C_v,trans + C_v,vib + C_v,rot (2)

S = S_vib + S_trans + S_rot (4)

These were applied as implemented in the Atomic Simulation Environment (ASE) Python package.³⁴ (Note that the expressions for monatomic and linear molecules are slightly different.) The rotational constants σ were assigned from the point groups.

2.4.2 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₈:

ΔH_Sx = H_Sx − xH_Sα (6)

The preferred experimental value for ΔH_sub is 100.416/8 = 12.552 kJ mol⁻¹, from experiments at 298 K.³ 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.

2.4.3 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 well-established and based on key work in ref. 36–38.

We attempt to minimise the Gibbs free energy

subject to the mass balance constraintwhere 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 formedand differentiated to form a set of equations defining 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 where P is the total pressure and the total quantity . 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°, μ^°_i(T) = μ_i(T,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 obtaincombining (21) and (22) we eliminate n_tand clean up the notation by denoting as Φ_i

Finally, to obtain the chemical potential of the mixture we note from (12) that for all i. Therefore

λ = μ_S, (25)

the normalised chemical potential of sulfur vapour on an atom basis. (A mathematical derivation is given in the Appendix.)

3 Results

3.1 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.

3.1.1 S₂. 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.

Table 1 Calculated and experimental bond length r in S₂. Experimental value is NIST/JANAF-recommended distance³

DFT functional	r/Å
PBE	1.911
PBEsol	1.903
LDA	1.895
PBE0	1.884
Experiment	1.889

---

3.1.2 S₃. The evolutionary algorithm process eliminated all but a C_2v non-linear chain for S₃. 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°).⁴¹ We have also considered the simple triangular allotrope, which is ∼0.5 eV higher in ground-state energy.

3.1.3 S₄. 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 and 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,43 Various cyclic and tetrahedral candidate structures yielded a relatively flat puckered ring with D_2d symmetry.

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

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

3.1.5 S₆. 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₃ cycles, or alternatively as a cluster of S₂ diatoms. This appears to be the D_3h “prism” structure identified by Wong et al.;⁴⁴ 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.

3.1.6 S₇. 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.

3.1.7 S₈. No evolutionary algorithm study was applied for S₈, 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,⁴⁵ and relaxed to form an isolated D_4d ring.

3.1.8 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⁻¹ atom, which is comparable to the energy difference between neighbouring cluster sizes.

Fig. 2 Ground-state energies from DFT of clusters included in study. Energies are relative to the energy for S₈ with each functional, and normalised to the number of atoms. A point is also included from reference data;³ 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₂ energy from PBE0 calculations agrees closely with this reference data.

3.2 Vibrational properties

Vibrational frequencies were calculated for all of the allotropes listed in Section 3.1; frequencies for S₂ and S₈ are listed in Table 2.

Table 2 Calculated and experimental vibrational frequencies for S₂ and S₈.³ All frequencies in cm⁻¹

	LDA	PBEsol	PBE0	PBE0 (scaled)	B3LYP	Expt
S2	716	713	751	721	714	724
S8	73	73	74	71	74	56
	73	73	75	72	74	56
	136	136	150	144	145	152
	136	136	150	144	145	152
	188	187	197	189	191	191
	188	187	197	189	191	191
	217	215	223	214	214	218
	228	228	248	238	242	243
	248	247	256	246	249	248
	248	247	256	246	249	248
	391	382	434	417	381	411
	418	411	454	436	407	437
	418	411	454	436	407	437
	473	467	492	472	455	471
	473	467	492	472	455	471
	479	474	493	473	461	475
	479	474	493	473	461	475
	486	482	497	477	470	475

---

3.2.1 Empirical corrections. Empirical scale factors were determined by fitting the frequencies to the experimental spectrum for S₈. Note that frequencies are linearly proportional to their corresponding zero-point energies and hence this may also be seen as fitting to zero-point 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⁻¹ error) with the stretching frequency of S₂, which was not used in the fit (Table 2).

Table 3 Optimal scale factors for exchange–correlation functionals, fitting to ground-state frequencies of S₈.³ 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

Functional	Scale factor	s/cm^−1	s/eV (ZPE)
LDA	1.0085	11.57	0.00072
PBEsol	1.0201	12.39	0.00077
PBE0	0.9596	6.41	0.00040
B3LYP	1.0332	11.05	0.00068

---

Fig. 3 Vibrational frequencies of S₈ calculated with various DFT functionals, compared with recommended experimental values.³

Least-squares fitting was carried out with the Levenberg–Marquardt algorithm as implemented in Scipy.^39,40

3.3 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 “lower-level” calculations predict a diverse mixture of phases, hybrid DFT strongly supports the dominance of S₈ and S₂, 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.

Fig. 4 Compositions of modelled S_x mixtures over range of equilibrium temperatures and pressures. Results are given for 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.

The corresponding free energies are also plotted in Fig. 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₈ 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₂ 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).

Fig. 5 Chemical potential of S vapours per mole of atoms, given at several pressures according to range of calculation methods. Data for S₂ and S₈ are also provided from the thermochemical literature.³ At low pressures, the free energy diverges by more than 50 kJ mol⁻¹ S atoms between the S₂ and S₈ allotropes at high temperatures, while at high pressures there is less variation. Results from hybrid DFT calculations with scaled frequencies closely track the minimal value from the literature, while the local and semi-local exchange–correlation functionals diverge from this data due to over-estimation of the formation energy of S₂.

Fig. 6 Chemical potential of S vapours over range of T,P, compared with individual allotropes. The equilibrium mixture is lower in energy than any single allotrope, but in most T/P regimes lies close to the chemical potential of S₂ or S₈. Data from vibrational calculations with PBE0 and empirically-corrected frequencies.

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₂ and S₈. The resulting plot (Fig. 8) shows that this has an impact ranging from around 1–4 kJ mol⁻¹, 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⁻¹ S atoms.

Fig. 7 Temperature-pressure map of approximations to free energy of mixture. At dashed line ; in shaded region the error in chemical potential μ associated with assuming a single phase S₂ or S₈ 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 Depression in chemical potential of sulfur vapour μ_S due to mixing and presence of minor allotropes. .

3.4 Parameterisation

For convenience, a parameterised fit has been generated for the chemical potential of S over the T, P range 400–1500 K, 10⁰ to 10⁷ Pa, incorporating an error function “switch” between S₂ and S₈ 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 iswhere μ_S8(T,P) = 7.620 × 10⁻¹ − 2.457 × 10⁻³T − 4.012 × 10⁻⁶T² + 1.808 × 10⁻⁹T³ − 3.810 × 10⁻¹³T⁴ + , μ_S2(T,P) = 1.207 − 1.848 × 10⁻³T − 8.566 × 10⁻⁷T² + 4.001 × 10⁻¹⁰T³ − 8.654 × 10⁻¹⁴T⁴ + . T_tr, the transition temperature obtained by solving is approximated by the polynomial T_tr = 5.077 × 10² + 7.272 × 10¹ log₁₀ P − 8.295(log₁₀ P)³ + 1.828(log₁₀ P)². The height of the Gaussian correction a(P) = 1.414 × 10³ − 2.041 × 10²log₁₀ P + 6.663 × 10¹(log₁₀ P)², and the more arbitrarily assigned width and offset parameters b = 10, c = 80, w = 100.

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 μ_S2, μ_S8 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⁻¹.

Fig. 9 Error of parameterisation in kJ mol⁻¹. Error is reduced to less than 1 kJ mol⁻¹, but is highly non-uniform. Parameterisation is recommended for convenient application over wide T–P ranges; the full equilibrium solution is required to correctly capture fine detail.

4 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₂ and S₈. 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⁻¹ atoms, while the minor phases contribute free energy of the order 1 kJ mol⁻¹ atoms. The resulting chemical potential data is made available through tabulated data, a parameterised model with error of the order 0.5 kJ mol⁻¹ atoms and through open-source code; the reference energy is compatible with the NIST-Janaf thermochemical tables for the solid α-sulfur phase.³ This phase is frequently used as a reference state for thermodynamic studies of defects and stability 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.

5 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-model and 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₈ 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 FHI-aims 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.

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

T/K	log10(p/Pa)
	1.00	1.67	2.33	3.00	3.67	4.33	5.00	5.67	6.33	7.00
100	4.73	4.88	5.04	5.20	5.36	5.52	5.68	5.84	6.00	6.16
150	2.29	2.53	2.77	3.01	3.25	3.49	3.72	3.96	4.20	4.44
200	−0.39	−0.07	0.25	0.57	0.89	1.21	1.53	1.85	2.17	2.49
250	−3.27	−2.87	−2.47	−2.08	−1.68	−1.28	−0.88	−0.48	−0.08	0.32
300	−6.34	−5.86	−5.39	−4.91	−4.43	−3.95	−3.47	−2.99	−2.51	−2.03
350	−9.58	−9.02	−8.46	−7.90	−7.34	−6.78	−6.23	−5.67	−5.11	−4.55
400	−12.97	−12.33	−11.69	−11.05	−10.41	−9.77	−9.13	−8.49	−7.85	−7.21
450	−16.50	−15.77	−15.05	−14.33	−13.61	−12.89	−12.17	−11.45	−10.73	−10.01
500	−20.20	−19.37	−18.56	−17.75	−16.94	−16.14	−15.33	−14.53	−13.73	−12.93
550	−24.24	−23.17	−22.22	−21.31	−20.40	−19.51	−18.62	−17.73	−16.85	−15.96
600	−29.74	−27.46	−26.12	−25.03	−24.01	−23.01	−22.03	−21.05	−20.08	−19.11
650	−37.54	−33.52	−30.62	−29.01	−27.78	−26.65	−25.56	−24.49	−23.42	−22.36
700	−45.63	−41.17	−36.83	−33.61	−31.81	−30.45	−29.22	−28.04	−26.87	−25.72
750	−53.78	−49.00	−44.23	−39.63	−36.36	−34.48	−33.03	−31.71	−30.43	−29.18
800	−61.99	−56.89	−51.79	−46.72	−41.99	−38.90	−37.03	−35.51	−34.10	−32.74
850	−70.27	−64.84	−59.43	−54.02	−48.67	−44.06	−41.31	−39.46	−37.88	−36.39
900	−78.59	−72.85	−67.11	−61.38	−55.67	−50.16	−46.04	−43.61	−41.79	−40.15
950	−86.97	−80.91	−74.85	−68.80	−62.75	−56.78	−51.43	−48.04	−45.84	−44.01
1000	−95.39	−89.01	−82.64	−76.26	−69.90	−63.57	−57.48	−52.84	−50.06	−47.98
1050	−103.86	−97.17	−90.47	−83.77	−77.09	−70.43	−63.88	−58.14	−54.50	−52.07
1100	−112.38	−105.36	−98.34	−91.33	−84.32	−77.34	−70.42	−63.91	−59.21	−56.29
1150	−120.94	−113.60	−106.26	−98.93	−91.60	−84.29	−77.03	−70.00	−64.26	−60.68
1200	−129.53	−121.88	−114.22	−106.57	−98.92	−91.29	−83.70	−76.25	−69.65	−65.25
1250	−138.17	−130.19	−122.22	−114.24	−106.28	−98.33	−90.41	−82.60	−75.33	−70.03
1300	−146.84	−138.54	−130.25	−121.96	−113.67	−105.40	−97.16	−89.01	−81.23	−75.04
1350	−155.55	−146.93	−138.32	−129.71	−121.10	−112.51	−103.95	−95.46	−87.25	−80.27
1400	−164.29	−155.36	−146.42	−137.49	−128.57	−119.66	−110.77	−101.95	−93.36	−85.72
1450	−173.06	−163.81	−154.56	−145.31	−136.07	−126.84	−117.63	−108.49	−99.53	−91.33

---

6 Appendix

6.1 Proof that G = λ

We define the molar Gibbs free energy of sulfur atoms in a molecular gas mixture asand substitute in (24)

(Notation the same as Section 2.4.3). From (12), μ_i = a_iλ and hence

Acknowledgements

The authors thank J. M. Skelton and J. M. Frost for useful discussions. USPEX/VASP calculations with PBE were carried out using the University of Bath's High Performance Computing facilities. Hybrid DFT calculations were carried out using ARCHER, the UK's national high-performance computing service, via our membership of the UK's HPC Materials Chemistry Consortium, which is funded by EPSRC (grant no. EP/L000202L000202). A. J. J. is part of the EPSRC Doctoral Training Center in Sustainable Chemical Technologies (grant no. EP/G03768XG03768X/11). The contribution of D. T. was supported by ERC Starting Grant 277757.

References

1. W. Nehb and K. Vydra, Sulfur, Ullmann’s Encyclopedia of Industrial Chemistry, 2006 DOI:10.1002/14356007.a25_507.pub2.
2. G. N. Lewis and M. Randall, J. Am. Chem. Soc., 1914, 36, 2468–2475.
3. M. J. Chase, J. Phys. Chem. Ref. Data, Monogr., 1998, 9, 1–1951.
4. H. Rau, T. Kutty and J. Guedes de Carvalho, J. Chem. Thermodyn., 1973, 5, 833–844.
5. P. Y. Yu and M. Cardona, Fundamentals of Semiconductors, Springer-Verlag, Berlin, Heidelberg, New York, 3rd edn, 2003.
6. D. M. Berg, R. Djemour, L. Gütay, G. Zoppi, S. Siebentritt and P. J. Dale, Thin Solid Films, 2012, 520, 6291–6294.
7. M. Lichtenstriger, H. U. Bölsterli and R. Nitsche, J. Phys. Chem. Solids, 1961, 21, 199–205.
8. D. Colombara, S. Delsante, G. Borzone, J. Mitchels, K. Molloy, L. Thomas, B. Mendis, C. Cummings, F. Marken and L. Peter, J. Cryst. Growth, 2013, 364, 101–110.
9. L. A. Burton, D. Colombara, R. D. Abellon, F. C. Grozema, L. M. Peter, T. J. Savenije, G. Dennler and A. Walsh, Chem. Mater., 2013, 25, 4908–4916.
10. A. J. Jackson and A. Walsh, J. Mater. Chem. A, 2014, 2, 7829.
11. V. Kosyak, N. B. Mortazavi Amiri, A. V. Postnikov and M. A. Scarpulla, J. Appl. Phys., 2013, 114, 124501.
12. J. J. Scragg, T. Ericson, T. Kubart, M. Edoff and C. Platzer-Björkman, Chem. Mater., 2011, 23, 4625–4633.
13. T. P. Martin, J. Chem. Phys., 1984, 81, 4426.
14. R. Steudel, Studies in Inorganic Chemistry, Elsevier, Amsterdam, 1984, vol. 5, pp. 3–37.
15. D. Hohl, R. O. Jones, R. Car and M. Parrinello, J. Chem. Phys., 1988, 89, 6823.
16. G. E. Quelch, H. F. Schaefer and C. J. Marsden, J. Am. Chem. Soc., 1990, 112, 8719–8733.
17. M. W. Wong and R. Steudel, Chem. Phys. Lett., 2003, 379, 162–169.
18. R. Steudel, Y. Stuedel and M. W. Wong, Top. Curr. Chem., 2003, 117–134.
19. W. Kohn and L. Sham, Phys. Rev. A, 1965, 140, 1133–1138.
20. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
21. J. Perdew, A. Ruzsinszky, G. Csonka, O. Vydrov, G. Scuseria, L. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
22. R. H. Hertwig and W. Koch, Chem. Phys. Lett., 1997, 268, 345–351.
23. A. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
24. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
25. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
26. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
27. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175–2196.
28. V. Havu, V. Blum, P. Havu and M. Scheffler, J. Comput. Phys., 2009, 228, 8367–8379.
29. A. R. Oganov and C. W. Glass, J. Chem. Phys., 2006, 124, 244704.
30. A. R. Oganov, A. O. Lyakhov and M. Valle, Acc. Chem. Res., 2011, 44, 227–237.
31. A. O. Lyakhov, A. R. Oganov, H. T. Stokes and Q. Zhu, Comput. Phys. Commun., 2013, 184, 1172–1182.
32. J. P. Merrick, D. Moran and L. Radom, J. Phys. Chem. A, 2007, 111, 11683–11700.
33. I. M. Alecu, J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2010, 6, 2872–2887.
34. S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng., 2002, 4, 56–66.
35. W. R. Smith and R. W. Nissen, Chemical Reaction Equilibrium Analysis: Theory and Algorithms, John Wiley & Sons, New York, 1982.
36. S. R. Brinkley, J. Chem. Phys., 1946, 14, 563.
37. S. R. Brinkley, J. Chem. Phys., 1947, 15, 107–110.
38. W. B. White, J. Chem. Phys., 1967, 46, 4171.
39. D. W. Marquardt, J. Soc. Ind. Appl. Math., 1963, 11, 431–441.
40. SciPy, Open source scientific tools for Python, 2001, http://www.scipy.org/, online; accessed 2015-03-06.
41. M. C. McCarthy, S. Thorwirth, C. A. Gottlieb and P. Thaddeus, J. Am. Chem. Soc., 2004, 126, 4096–4097.
42. M. S. Boumedien, J. Corset and E. Picquenard, J. Raman Spectrosc., 1999, 472, 463–472.
43. P. Hassanzadeh and L. Andrews, J. Phys. Chem., 1992, 96, 6579–6585.
44. M. W. Wong, Y. Steudel and R. Steudel, J. Chem. Phys., 2004, 121, 5899–5907.
45. L. A. Burton and A. Walsh, J. Phys. Chem. C, 2012, 116, 24262–24267.

---

Footnotes

† Electronic supplementary information (ESI) available: Tabulated free energy and enthalpy data. Additional data and code available in external repositories with DOI: 10.5281/zenodo.28536; DOI: 10.6084/m9.figshare.151373; DOI: 10.6084/m9.figshare.1513833. See data access statement for more information. See DOI: 10.1039/c5sc03088a

‡ Current address: EPFL Valais Wallis, EPFL LSMO, Rue de l'Industrie 17, Case postale 440, CH-1951 Sion, Switzerland.

§ Note that the implementation of B3LYP in FHI-aims uses a parameterisation of the local density contribution based on the random phase approximation in order to match values obtained with Gaussian, another quantum chemistry code.²²

---