Chemical and Lattice Stability of the Tin Sulfides.

The tin sulfides represent a materials platform for earth-abundant semiconductor technologies. We present a first-principles study of the five known and proposed phases of SnS together with SnS2 and Sn2S3. Lattice-dynamics techniques are used to evaluate the dynamical stability and temperature-dependent thermodynamic free energy, and we also consider the effect of dispersion forces on the energetics. The recently identified π-cubic phase of SnS is found to be metastable with respect to the well-known orthorhombic Pnma/Cmcm equilibrium. The Cmcm phase is a low-lying saddle point between Pnma local minima on the potential-energy surface and is observed as an average structure at high temperatures. Bulk rocksalt and zincblende phases are found to be dynamically unstable, and we show that whereas rocksalt SnS can potentially be stabilized under a reduction of the lattice constant the hypothetical zincblende phase proposed in several previous studies is extremely unlikely to form. We also investigate the stability of Sn2S3 with respect to SnS and SnS2 and find that both dispersion forces and vibrational contributions to the free energy are required to explain its experimentally observed resistance to decomposition.

Introduction

The tin sulfides are a technologically important family of earth-abundant optoelectronic materials comprising tin monosulfide (SnS), tin disulfide (SnS2), and tin sesquisulfide (Sn2S3). SnS has long been explored as a sustainable candidate photovoltaic (PV) absorber material.[1−3] SnS2 is a 2D semiconductor that can be exfoliated to give individual nanosheets,[4] and Sn2S3 is an example of a comparatively rare mixed-valence binary compound, possesses an unusual 1D bonding structure,[5] and is predicted to show ambipolar dopability.[6]

The SnS system has a rich phase space.[7] SnS is known to form several (meta)stable phases, which has led to a degree of uncertainty over which phase(s) are obtained under different growth conditions. There are currently five reported or proposed phases, viz. the ground-state orthorhombic Pnma phase,[8] a high-temperature Cmcm phase,[9] and three cubic phases: rocksalt,[10] zincblende,[11−13] and the recently reported P213 (“π-cubic”) phase with a 64-atom primitive cell.[14−17] The sesquisulfide Sn2S3 has also been a source of confusion, as it is relatively easy to prepare and clearly distinguishable from the other tin sulfides[18] yet is frequently predicted to be unstable with respect to decomposition into SnS and SnS2 in theoretical studies (e.g., as in the current Materials Project[19] entry, mp-1509[20]). Both issues are important for contemporary PV research because not only are phase impurities highly likely to play a role in the underwhelming performance of current SnS-based devices[18,21] but also tin sulfides may also form as impurities during the growth of more complex multicomponent semiconductors such as Cu2ZnSnS4 (CZTS).[22]

First-principles materials modeling, for example, within the ubiquitous Kohn–Sham density functional theory formalism,[23,24] affords a versatile means of exploring the subtleties of the equilibria between competing phases. However, to enumerate such a complex phase space, two common approximations made in contemporary modeling studies may need to be revisited, namely: (1) the poor description of weak dispersion interactions afforded by most generalized-gradient approximation (GGA) density functionals, which may be an issue in describing the sulfide phases with prominent nonbonding interactions,[5,25] and (2) the omission of vibrational contributions to the temperature-dependent free energy in thermodynamic models.[26,27] Consideration of the lattice dynamics would also enable the assessment of the bulk dynamical stability of the different phases, providing important insight into which phases are likely or unlikely to form under typical growth conditions.

In this work, we address both issues through a consistent set of first-principles lattice-dynamics calculations on the five known and proposed phases of SnS together with SnS2 and Sn2S3, performed using both the PBEsol GGA[28] and the dispersion-corrected[29] PBEsol + D3 functionals. We establish the dynamical and thermodynamic stability ordering of the competing SnS phases and provide conclusive evidence that the various experimental reports of cubic SnS are likely to be the rocksalt or π-cubic phases rather than a zincblende polymorph. We also consider the thermodynamics of the decomposition of Sn2S3 and show that to predict its stability requires an accurate description of the weak interactions in its structure as well as accounting for the vibrational contributions to its free energy.

Computational Methods

First-principles calculations were carried out using the pseudopotential plane-wave density functional theory (DFT) formalism, as implemented in the Vienna Ab initio Simulation Package (VASP) code.[30]

We performed calculations using the PBEsol functional[28] and PBEsol with the DFT-D3 dispersion correction[29] applied (i.e., PBEsol + D3). PBEsol was selected due to its being shown in a number of studies to provide a good description of the structural and vibrational properties of solids at a moderate computational cost.[31−33]

Initial models built from published crystal structures[8,9,15,34−37] were fully optimized to a tolerance of 10–2 eV Å–1 on the forces, using a plane-wave cutoff of 550 eV and carefully converged k-point sampling (Table ). A tight energy-convergence criterion of 10–8 eV was applied during the electronic minimization. Projector augmented wave (PAW) pseudopotentials[38,39] were used to model the core electrons, with the Sn 5s, 5p, and 4d and the S 3s and 3p electrons being included in the valence shell.

Table 1

Phonon Supercells and k-Point Sampling Used for the Geometry Optimizations and Supercell Finite-Displacement Phonon Calculations Performed in This Worka

		k-points
system	phonon supercell (no. atoms)	optimization	phonon
α-Sn		12 × 12 × 12 Γ-MP
β-Sn		16 × 16 × 32 Γ-MP
S		3 × 1 × 1 Γ-MP
SnS (Pnma)	6 × 1 × 6 (384)	8 × 4 × 8 MP	2 × 4 × 2 MP
SnS (cubic)	2 × 2 × 2 (512)	2 × 2 × 2 Γ-MP	1 × 1 × 1 Γ-MP
SnS (Cmcm)	6 × 1 × 6 (384)	8 × 4 × 8 MP	2 × 4 × 2 MP
SnS (rocksalt)	4 × 4 × 4 (128)	12 × 12 × 12 Γ-MP	3 × 3 × 3 Γ-MP
SnS (zincblende)	4 × 4 × 4 (128)	8 × 8 × 8 Γ-MP	2 × 2 × 2 Γ-MP
SnS2	6 × 6 × 2 (144)	8 × 8 × 6 Γ-MP	2 × 2 × 3 Γ-MP
Sn2S3	2 × 4 × 2 (320)	4 × 8 × 3 Γ-MP	2 × 2 × 2 Γ-MP

MP: Monkhorst–Pack mesh;[40] Γ-MP: Γ-centered MP mesh.

The PAW projection was performed in reciprocal space. Nonspherical contributions to the gradient corrections inside the PAW spheres were accounted for, and the precision of the charge-density grids was set automatically to avoid aliasing errors.

Lattice-dynamics calculations were performed on the optimized structures using the Phonopy package,[41,42] which was used to set up and post-process supercell finite-displacement phonon calculations.[43] VASP was used as the force calculator, and an additional charge-density grid containing eight times the number of points as the standard one was used to evaluate the augmentation charges. The supercell expansions used to determine the second-order force-constant matrices are listed in Table . For all models apart from Cmcm SnS, the force constants were calculated in expansions of the primitive unit cells. The force constants for Cmcm SnS were calculated in an expansion of the conventional cell, and a transformation matrix to the primitive cell was applied when evaluating its phonon dispersion curves.

Phonon density of states (DoS) curves were obtained by interpolating the phonon frequencies onto a uniform 48 × 48 × 48 Γ-centered q-point mesh and using the linear tetrahedron method for the Brillouin-zone integration. The phonon dispersions were obtained by interpolating the phonon frequencies along lines of q-points passing between the high-symmetry points in the Brillouin zones of the primitive unit cells. Thermodynamic functions were computed from a set of phonon frequencies evaluated on the same 48 × 48 × 48 q-point grids as were used to obtain the phonon DoS curves.

Results and Discussion

Figure shows the structures of the seven tin sulfide compounds examined in this study after structural optimization with PBEsol.

Figure 1

PBEsol-optimized structures of the Pnma, Cmcm, π-cubic, rocksalt, and zincblende phases of SnS, SnS2, and Sn2S3. These images were prepared using the VESTA software.[44]

The structural diversity in the SnS phase space is primarily driven by the ability of Sn to adopt two different oxidation states, viz. Sn(II) and Sn(IV).

The Pnma, Cmcm, and π-cubic phases of SnS are built up from distorted Sn(II) tetrahedra with a stereochemically active lone pair occupying one of the four coordination sites. The two orthorhombic phases are composed of layered 2D SnS sheets, with the lone pairs projecting into the interlayer space and facilitating a weak interaction along the crystallographic b axes. The low-temperature Pnma phase can be thought of as a small distortion of the more symmetric high-temperature Cmcm phase, and a phase transition between them occurs above 875 K.[45] The local coordination in the π-cubic phase is similar to that in the orthorhombic phases, but the Sn–S bonding in this structure instead forms a 3D network. In the higher-symmetry rocksalt and zincblende structures, the Sn lone pair must be accommodated within the close-packed structure, and the stability of this arrangement depends sensitively on the local chemical environment of the cation.[46]

In SnS2, symmetric edge-sharing Sn(IV) octahedra form 2D layers separated by a van der Waals’ gap. Sn2S3 is a mixed oxidation-state phase containing equal proportions of Sn(II) and Sn(IV) and as such displays structural motifs from both SnS and SnS2. The Sn2S3 structure consists of 1D chains of Sn(IV) octahedra capped by tetrahedral Sn(II), with weak interchain interactions facilitated by the Sn(II) lone pairs.

The optimized lattice parameters obtained with both exchange-correlation functionals (Table ) were found to agree well with experimental measurements, and the PBEsol results were a good match to previous all-electron DFT calculations using the same functional.[5] With the exception of zincblende SnS, the dispersion correction consistently yielded a smaller lattice volume, with a particularly notable contraction in the nonbonded c direction of SnS2. The reduction in volume implies that the attractive part of the dispersion correction is dominant for most of the structures, whereas for the zincblende phase the dispersion correction produces a net repulsive effect. This can be explained by the interaction of the Sn(II) lone pair with the symmetric coordination environment: The Sn–S bond length in the rocksalt structure (PBEsol: 2.874 Å, PBEsol + D3: 2.856 Å) is considerably longer than in the zincblende structure (PBEsol: 2.778 Å, PBEsol + D3: 2.780 Å), which would lead to a stronger (repulsive) interaction between the lone pair and the bonding electrons.

Table 2

Optimized Lattice Constants of the Models Considered in This Work Obtained with the PBEsol and PBEsol + D3 Exchange-Correlation Functionalsa

	PBEsol			PBEsol + D3			expt.
	a (Å)	b (Å)	c (Å)	a (Å)	b (Å)	c (Å)	a (Å)	b (Å)	c (Å)
α-Sn	6.538			6.506			6.491[35]
β-Sn	5.829		3.163	5.810		3.143	5.820[35]		3.175[35]
S	10.556	12.942	24.408	10.075	12.493	23.966	10.437[34]	12.845[34]	24.369[34]
SnS (Pnma)	4.250	11.082	3.978	4.220	10.976	3.958	4.33[8]	11.18[8]	3.98[8]
SnS (cubic)	11.506			11.405			11.603[15]
SnS (Cmcm)b	4.037	11.282	4.039	4.018	11.186	4.009	4.148[9]	11.480[9]	4.177[9]
SnS (rocksalt)	5.747			5.712
SnS (zincblende)	6.416			6.420
SnS2	3.651		6.015	3.639		5.721	3.638[36]		5.880[36]
Sn2S3	8.811	3.766	13.813	8.633	3.760	13.663	8.878[37]	3.751[37]	14.020[37]

Experimental values are given where available for comparison.

Lattice constants in ref (9) are reported at 905 K and are therefore expected to include significant thermal expansion.

The formation energies obtained with both functionals, referenced to solid S and β-Sn, are compared in Table . These values were used to draw a convex hull including the newly discovered π-cubic phase (Figure ).

Table 3

Energetic and Dynamical Stability of the Pnma, π-Cubic, Cmcm, Rocksalt, and Zincblende Phases of SnS, SnS2, and Sn2S3a

	EF (kJ mol–1 per F.U.)
system	PBEsol	PBEsol + D3	expt.	bulk dynamically stable?
SnS (Pnma)	–90.59	–95.00	–100 to −108[47,48]	yes
SnS (π-cubic)	–88.40	–92.46		yes
SnS (Cmcm)	–88.83	–93.36		no
SnS (rocksalt)	–82.05	–84.24		no
SnS (zincblende)	–18.36	–14.45		no
SnS2	–120.00	–127.47	–148 to −182[48−50]	yes
Sn2S3	–208.90	–222.32	–249 to −297[48−50]	yes

Calculated formation energies (EF) of each compound are compared to experimental data where available, while the last column lists the bulk dynamical stabilities assessed from the harmonic phonon dispersions (see Figure ).

Figure 2

Thermodynamic convex hulls showing the calculated phase stability of the seven compounds in the SnS system based on formation energies obtained with (a) PBEsol and (b) PBEsol + D3. ZB and RS refer to the zincblende and rocksalt phases of SnS. The optimized crystal-structure parameters are reported in Table .

Figure 3

Phonon dispersion and density of states curves for the Pnma (a), Cmcm (b), rocksalt (c), zincblende (d), and π-cubic (e) phases of SnS, SnS2 (f), and Sn2S3 (g) calculated within density functional theory with electron exchange and correlation described by PBEsol (blue) and PBEsol + D3 (red).

Both functionals predict the ground state of SnS to be the Pnma phase, with the high-temperature Cmcm phase calculated to be 1.76 kJ mol–1 per SnS formula unit (F.U.; PBEsol) and 1.63 kJ mol–1 per F.U. (PBEsol + D3) higher in energy. The π-cubic phase is metastable with respect to this equilibrium, being 2.19 kJ mol–1 per F.U. higher in energy with PBEsol and 2.54 kJ mol–1 per F.U. higher in energy with PBEsol + D3. The rocksalt phase is around four times higher in energy than the π-cubic phase, at 8.54 kJ mol–1 per F.U. with PBEsol and 10.75 kJ mol–1 per F.U. with PBEsol + D3. The zincblende phase has a much smaller formation energy and is thus very high in energy relative to all the other phases (PBEsol: 72.23 kJ mol–1 per F.U., PBEsol + D3: 80.55 kJ mol–1 per F.U.).

As in other calculations,[6,20] PBEsol predicts Sn2S3 to be above the convex hull, making it unstable with respect to decomposition into Pnma SnS and SnS2.

Whereas PBEsol predicts it to be 1.68 kJ mol–1 above the hull, PBEsol + D3 predicts it to be much closer at 0.14 kJ mol–1; this suggests that the dispersion correction selectively stabilizes the sesquisulfide, most likely due to its improved description of the weak interactions between the bonded chains. Both values, but in particular the PBEsol + D3 one, are on the order of differences in the vibrational zero-point energy, a point that we return to below.

Finally, we note that the similar calculations carried out previously,[5] although not including the π-cubic phase, were referenced to α-Sn; in the present calculations, both functionals predicted β-Sn to be lower in energy, but a similar comparison of the formation energies and a convex hull referenced to α-Sn yielded the same qualitative results (Table S1 and Figure S1 in the Supporting Information).

We next consider the dynamical stability of the seven sulfides. In the harmonic phonon model, displacement of atoms from their equilibrium positions along a normal-mode coordinate Q leads to a change in potential energy , where ω is the vibrational frequency. If a structure is a minimum on the potential-energy surface, then U and ω are ≥0. (The three acoustic modes have ω = 0 at the center of the Brillouin zone (Γ).) If, on the contrary, a structure is a potential-energy maximum (e.g., a saddle point), atomic motion along one or more vibrations will lead to a lowering of the energy and hence to a complex (imaginary) frequency.

In the present calculations, an analysis of the phonon dispersions (Table , Figure ), that is, the frequencies of the 3 N phonon bands as a function of the phonon wavevector, indicates the Pnma and π-cubic phases of SnS to be dynamically stable, whereas the Cmcm, rocksalt, and zincblende phases all have imaginary modes.

The instability of the high-temperature Cmcm phase is anticipated, as the isostructural Cmcm phase of SnSe exhibits a similar phonon dispersion with imaginary modes at the same reciprocal-space wavevectors.[51] The Cmcm structure is effectively a low-lying saddle point on the potential-energy surface, connecting equivalent distorted Pnma minima. Above the phase-transition temperature, sufficient thermal energy is available for the system to rapidly “hop” between the minima, and the higher-symmetry phase is observed as a crystallographic average structure with significant local distortions. The implication of this is that the Cmcm phase is only observed above the phase-transition temperature and cannot be isolated under ambient conditions, for example, by rapid quenching.

The phonon instabilities in the rocksalt phase correspond to a four-fold enlargement of the two-atom primitive cell, suggesting a transition toward the Pnma phase. Whereas the imaginary modes in this structure (and also in the Cmcm phase) are restricted to well-defined phonon wavevectors and represent a small fraction of the overall phonon density of states (DoS), the zincblende instabilities account for a significant proportion (∼1/3) of the integrated phonon DoS. Taken together with its predicted high formation energy, this phase is extremely unlikely to form, even under nonambient conditions such as high pressure.

To confirm this, we calculated phonon dispersions of the rocksalt and zincblende phases for volume compressions of up to 20%, corresponding to a 0.41 to 0.46 Å reduction in the lattice constants (Figure , see also Figure S2 in the Supporting Information). Under moderate compression, the imaginary modes in the rocksalt dispersion harden, indicating that the bulk phase can be dynamically stabilized under pressure. In contrast, compression of the zincblende phase was found to induce further phonon softening.

Figure 4

Phonon dispersion and density of states curves for the rocksalt (a) and zincblende (b) phases of SnS calculated at cell volumes ranging from the athermal equilibrium (cyan) to a 20% compression (blue). Both sets of curves were computed with PBEsol + D3, and the results from identical calculations with PBEsol are presented in Figure S2 (Supporting Information).

These results strongly suggest that cubic SnS has been misassigned as zincblende in various past studies. For example, the cubic lattice constant reported in ref (13), 5.846 Å, is much closer to those of our rocksalt models (PBEsol: 5.747 Å, PBEsol + D3: 5.712 Å) than to those of our zincblende ones (PBEsol: 6.416 Å, PBEsol + D3: 6.420 Å). The rocksalt lattice constants are also a better match for the lattice constant of 6.00 Å measured for SnS films grown heteroepitaxially on NaCl. While it is possible that a rocksalt phase could be grown on certain substrates (or under pressure), we would argue that bulk cubic SnS is more likely to be the π-cubic phase, as (1) our calculations show that this is dynamically stable in bulk; (2) the tetrahedral coordination environment in this structure is more likely to produce the tetrahedral morphology commonly observed in SnS nanoparticles;[11,16] and (3) it has been shown that the powder X-ray diffraction pattern of this phase may have been misinterpreted as that of a zincblende phase in the past.[15,17]

Both SnS2 and Sn3S3 were found to be dynamically stable, as is expected given that these phases can be prepared as bulk single crystals.[18]

Comparing the phonon dispersion and DoS curves calculated using PBEsol and PBEsol + D3 reveals a general shift of the latter to higher frequencies due to the smaller predicted unit-cell volumes.[32] The exception is zincblende SnS, where the slightly larger lattice constant predicted with PBEsol + D3 leads to a corresponding softening of the phonon frequencies. For all seven compounds, however, the main features, in particular, the presence or absence of imaginary modes, are qualitatively similar.

As well as assessing dynamical stability, lattice-dynamics calculations also allow vibrational contributions to the free energy to be taken into account when assessing thermodynamic stability. Within the harmonic model, the constant-volume (Helmholtz) free energy, A, is given bywhere Ulatt is the (athermal) lattice energy (i.e., the total energy calculated using DFT, in the present study), Uvib is the vibrational internal energy, and Svib is the vibrational entropy. A(T) is typically evaluated from the lattice energy and phonon frequencies using the bridge relation from statistical mechanicswhere Z(T) is the thermodynamic partition function, kB is Boltzmann’s constant, the phonon frequencies ω are indexed by a reciprocal-space wavevector q and a band index v, and the sum over the phonon Brillouin zone is normalized by the number of unit cells in the crystal, N, equivalent to the number of wavevectors included in the summation.

The Helmholtz energies of the Cmcm, π-cubic, and rocksalt phases of SnS, relative to the Pnma phase, are shown as a function of temperature in Figure , and free-energy differences calculated at 0, 300, and 900 K are listed in Table .

Figure 5

Thermodynamic stability of the Cmcm (red), π-cubic (blue), and rocksalt (yellow) phases of SnS as a function of temperature relative to the low-temperature ground-state Pnma phase, based on the constant-volume (Helmholtz) free energies, A, calculated with PBEsol (solid lines) and PBEsol + D3 (dashed lines).

Table 4

Relative Energies of the Five Phases of SnS Compared with the Ground-State Pnma Phasea

	(kJ mol–1 per F.U.)
	PBEsol				PBEsol + D3
	ΔUlatt	ΔA0K	ΔA300K	ΔA900K	ΔUlatt	ΔA0K	ΔA300K	ΔA900K
Pnma	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
Cubic	2.19	2.23	2.01	1.61	2.54	2.56	2.30	1.79
Cmcm	1.76	1.39	0.45	–1.93	1.63	1.29	0.49	–1.55
rocksalt	8.54	7.63	6.24	2.11	10.75	9.85	8.53	4.51
zincblende	72.23				80.55

Lattice energies, ΔUlatt, and constant-volume (Helmholtz) free energies, ΔA, at 0, 300, and 900 K calculated with PBEsol and PBEsol + D3 are compared. Energies are given in kJ mol–1 per SnS formula unit. Because of the large proportion of imaginary phonon modes in its phonon density of states, we did not calculate free energies for the zincblende phase.

Although, as noted above, the high-temperature Cmcm phase is an average structure, because the imaginary modes form a small part of the overall phonon DoS we make the approximation that the static configuration provides a reasonable representation of its (average) internal energy and lattice dynamics. Doing so predicts the free energies of the Pnma and Cmcm phases to cross at around 400–450 K, which, although qualitatively correct, is considerably lower than the 875 K at which the phase transition is observed experimentally. This can be ascribed to the approximations made in the present calculations, including the approximate treatment of the physics of the phase transition, and possibly also the neglect of volume expansion and other forms of anharmonicity.[51]

The π-cubic phase is predicted to remain metastable up to 1000 K, with an approximately constant free-energy difference with respect to the Pnma phase. In contrast, the free energy of the rocksalt phase, computed under the same assumptions as for the Cmcm phase, falls with temperature but remains above the Pnma phase in energy across the range of temperatures considered.

Finally, we also consider the effect of vibrational contributions to the free energy on the stability of Sn2S3. Figure compares the temperature dependence of the free energy, ΔA, of the disproportionation reactionIn light of the free-energy crossing evident in Figure , we consider decomposition to both Pnma and Cmcm SnS. The calculated reaction energies at representative temperatures are given in Table .

Figure 6

Constant-volume (Helmholtz) free energy, A, for the decomposition of Sn2S3 into SnS2 and SnS as a function of temperature, calculated with PBEsol (blue) and PBEsol + D3 (red). On the basis of the stability ordering of the SnS phases (Figures /5, Tables /4), we consider decomposition to both the Pnma (solid lines) and Cmcm (dashed lines) phases of SnS.

Table 5

Reaction Energies for the Decomposition Reaction Sn2S3 → SnS2 + SnS, Calculated Based on the Differences in the Athermal Lattice Energies (ΔUlatt) and the Constant-Volume (Helmholtz) Free Energies (ΔA) at 0, 300, and 900 Ka

	(kJ mol–1 per F.U.)
	PBEsol	PBEsol + D3
ΔUlatt	–1.68	–0.14
ΔA0K	–1.51	0.12
ΔA300 K	–0.98	0.69
ΔA900 K	0.37	2.26
Cmcm SnS: ΔA900 K	–1.57	0.71

Negative energy implies that Sn2S3 is unstable, while a positive energy indicates the reverse. The energies are given with respect to Pnma SnS, as this was determined to be the energetic ground state, while the 900 K free energy calculated with respect to the Cmcm phase, which has the lowest free energy of the SnS phases above ∼400 K, is also given.

PBEsol and PBEsol + D3 both predict the decomposition to the Pnma phase to become less favorable with temperature (i.e., the free energy increases), whereas decomposition to the Cmcm phase becomes more favorable (i.e., the energy decreases). The PBEsol free energies predict decomposition to the Pnma phase to become unfavorable above ∼650 K; however, decomposition to the Cmcm phase is a lower energy pathway above ∼400 K, and the reaction remains favorable.

With PBEsol + D3, on the contrary, the inclusion of zero-point energy renders decomposition unfavorable at 0 K, and it remains so up to 1000 K.

Given the theoretical reports predicting Sn2S3 to be unstable with respect to Pnma SnS and SnS2, it is of interest to ascertain whether this result can also be obtained with other dispersion correction methods. We therefore performed calculations on Pnma and Cmcm SnS, SnS2, and Sn2S3 with the newer variant of PBEsol + D3 with the Becke–Johnson damping scheme[52] (see Table S2 and Figures S3 and S4 in the Supporting Information). These calculations yielded the same qualitative result, that is, that the dispersion correction brings Sn2S3 closer to the convex hull than uncorrected PBEsol, and the subsequent inclusion of the vibrational contributions to the free energy predicts it to be thermodynamically stable up to 1000 K. Although the Becke–Johnson-damped DFT-D3 is similar to the “bare” D3, we note that the majority of the other well-known dispersion corrections, for example, DFT-D2,[53] the Tkatchenko–Scheffler method (DFT-TS),[54] and the more sophisticated many-body dispersion technique,[55,56] do not have scaling constants optimized for PBEsol. To investigate other common dispersion corrections would thus also require the exploration of other GGA functionals, which we defer to a future study.

Although, as noted above regarding the predicted Pnma-to-Cmcm phase-transition temperature, these calculations may still be missing some important effects, the predicted stability of Sn2S3 with PBEsol + D3 is in keeping with experimental findings. These results therefore highlight the potential importance of both a sufficient description of dispersion interactions and of vibrational contributions to the free energy in assessing the phase stability of the tin sulfides.

Conclusions

This work has provided detailed insight into the thermodynamic and dynamical stability of the seven currently known and proposed compounds in the SnS phase space, addressing several key outstanding questions.

The recently discovered π-cubic phase is metastable with respect to the orthorhombic Pnma/Cmcm equilibrium. The bulk rocksalt phase is higher in energy and is dynamically unstable but could potentially be stabilized under pressure or epitaxial strain. Our calculations show conclusively that the hypothetical zincblende phase is both energetically and dynamically unstable, and we suggest that reports of this phase be reassessed as either of the other two cubic phases.

Finally, we also show that accurately modeling the tin sulfide phase diagram, in particular, reproducing the experimentally observed stability of the sesquisulfide, may require theoretical techniques that afford a good treatment of dispersion interactions and possibly also the inclusion of contributions to the energetics from lattice dynamics. This result provides a baseline for further theoretical characterization of this less-well-studied tin sulfide phase.