Creating new layered structures at high pressures: SiS2

Old and novel layered structures are attracting increasing attention for their physical, electronic, and frictional properties. SiS2, isoelectronic to SiO2, CO2 and CS2, is a material whose phases known experimentally up to 6 GPa exhibit 1D chain-like, 2D layered and 3D tetrahedral structures. We present highly predictive ab initio calculations combined with evolutionary structure search and molecular dynamics simulations of the structural and electronic evolution of SiS2 up to 100 GPa. A highly stable CdI2-type layered structure, which is octahedrally coordinated with space group surprisingly appears between 4 and up to at least 100 GPa. The tetrahedral-octahedral switch is naturally expected upon compression, unlike the layered character realized here by edge-sharing SiS6 octahedral units connecting within but not among sheets. The predicted phase is semiconducting with an indirect band gap of about 2 eV at 10 GPa, decreasing under pressure until metallization around 40 GPa. The robustness of the layered phase suggests possible recovery at ambient pressure, where calculated phonon spectra indicate dynamical stability. Even a single monolayer is found to be dynamically stable in isolation, suggesting that it could possibly be sheared or exfoliated from bulk -SiS2.


INTRODUCTION
SiS 2 is a member of an important family of group IV-VI AB 2 compounds made of light elements including well-known materials such as CO 2 1-14 , SiO 2 [15][16][17][18][19][20][21][22] , GeO 2 [23][24][25][26][27] and CS 2 [28][29][30][31] .Limited to relatively low pressures, the structural evolution of SiS 2 has also been accurately described. 32The ambient-pressure stable phase known as NP-SiS 2 has orthorhombic Ibam structure and consists of distorted edge-sharing tetrahedra forming 1D chains which interact via weak van der Waals forces [32][33][34][35][36][37][38][39][40] .At 2.8 GPa, a first high-pressure phase HP1-SiS 2 32,36 appears, with monoclinic space group P 2 1 /c, by interconnection of the NP-SiS 2 chains to form 2D layers 32 .The connectivity pattern changes from each tetrahedron sharing two edges in the NP phase to sharing one edge and two corners in HP1 32 .The very same structure was recently predicted to be stable at pressures beyond 30 GPa in CS 2 31 .Further increase of pressure leads to the HP2-SiS 2 phase at 3.5 GPa 32,36 , again with P 2 1 /c space group, and again with edge-and corner-sharing tetrahedra, in this case, however, with a 3D connectivity network with large cavities 32 .Finally, at 4 GPa, a tetragonal I 42d structure denoted as HP3-SiS 2 takes over, a phase formed by strictly corner-sharing tetrahedra that are slightly distorted and span the whole threedimensional space [32][33][34][36][37][38] . The sae structure as HP3-SiS 2 is adopted by CO 2 phase V at high pressures 5 and it can be viewed as a partially collapsed version of SiO 2 β-cristobalite 6 .The phase diagram of silicon disulfide at pressures higher than 6 GPa remains unknown experimentally.Evers et al. 32 reasonably hypothesized that at higher pressures such as 10-20 GPa, SiS 2 might adopt a six-fold, octahedral Si coordination, accompanying a density increase.An alternative scenario might be the creation of a denser tetrahedral phase such as coesite in SiO 2 .Work is clearly needed to explore the high pressure phases of SiS 2 .While of course future experimental work is called for, experience of the last decade has shown that state-ofthe-art density-functional based crystal structure search can be extremely predictive.
By means of a well-tested, highly reliable protocol consisting of ab initio electronic structure calculations combined with evolutionary search for crystal structure prediction 41,42 and with constant-pressure molecular dy-namics (MD) simulations, we undertook a fresh theoretical exploration of high-pressure phases of SiS 2 up to 100 GPa.While confirming first of all the presence and stability of the known low-pressure tetrahedral phases, our study predicts three new structures with octahedral, six-fold Si coordination at higher pressures.While that confirms previous expectations, the surprise is that the new structures are layered.We analyze their structural and electronic properties in detail at all pressures, including properties of a single octahedral monolayer which is predicted to survive in a metastable state, should one succeed to shear off and exfoliate it away from the bulk layered structure.More generally, these findings suggest the idea to create, shear and possibly stabilize one or more SiS 2 semiconducting monolayers with over 2 eV band gap on a substrate in some future realizations.

A. Phases of SiS2 above the tetrahedral regime
From the extensive evolutionary search up to 100 GPa, MD simulations and structural optimizations (details given in Methods) we found three most promising candidates for high-pressure phases of SiS 2 beyond 6 GPa.These three lowest-enthalpy phases have space groups P 3m1, P 6 3 mc and R 3m -Fig. 1 and are very similar to each other.They all consistently emerged at every investigated pressure (10, 30, 60 and 100 GPa) in the evolutionary search.All three are layered, formed by separate sheets of edge-sharing SiS 6 octahedra with different stacking -Fig. 1.The layered character, which is of potential interest in view of a possibly facile frictional sliding under high-pressure shear, is remarkable, since in most materials one tends to associate higher density with "more 3D" networks -precisely the trend observed in the SiS 2 tetrahedral region (from NP to HP1, HP2 and HP3 between 0 and 4 GPa).
In these novel layered phases each of the octahedra shares six out of its twelve edges with six surrounding octahedra.The most stable phase is P 3m1 -shown in Fig. 2 and in Fig. 1 (left pictures), which contains only one SiS 2 unit in the primitive cell.It is isostructural to CdI 2 (C6 structure) 44 where sheets are stacked directly above each other.The same structure is found at ambient pressure in chemically similar systems with heavier atoms such as SnS 2 45 or SiTe 2 46 as well as in several other chalcogenides, iodides, chlorides and bromides 44,[47][48][49][50] , or even in BeH 2 51 .The second lowest-enthalpy phase P 6 3 mc contains two layers per unit cell which are mutually staggered and reflected with respect to each other -Fig. 1 (middle pictures).Finally, the least stable R 3m phase shown in Fig. 1 (right pictures) contains three sheets per unit cell and its structure is identical to ambient-conditions CdCl 2 .The stacking pattern of the structures is illustrated in Fig. 1 by lateral (upper pictures) and axial (lower pictures) views with respect to sheets.The arrangement of sheets in all three phases may be represented by specifying the stacking pattern of individual single-atomic layers of Si and S atoms, as far as these form planes of closedpacked atoms that are placed regularly one on top of another.The stacking pattern of the P 3m1 structure in this single-atom layer notation is /AbC/ 44 , while for P 6 3 mc it is /AbC/AcB/ and for R 3m it is /AbC/BcA/CaB/, where small letters denote layers of Si atoms and capital letters are for planes of S atoms (see legend in Fig. 1 lateral view on P 3m1).Structural data of the three phases at 10 GPa, their densities and band gaps calculated with generalized-gradient approximation (GGA) are summarized in Table A.

Stability of high-pressure octahedral layered structures
Calculated T =0 enthalpies show that the sequence of transitions for tetrahedral phases occurring at 3.5 GPa (HP1→HP2) and 4 GPa (HP2→HP3) experimentally 32 at room temperature is reasonably well reproduced within our calculations, which predict these two transitions to take place (at zero temperature) at 2.8 and 5.4 GPa, respectively -Fig. 3 (upper panel).However, all three six-coordinated phases become more stable than HP2-SiS 2 above cca.4 GPa.Hence, as far as HP3 crosses the HP2 curve only at 5.4 GPa, the calculations predict that the HP3 phase should be metastable at T = 0 at all pressures.Its experimental observation between 4-6 GPa 32 might be due to temperature, or to kinetics.The NP-HP1 transition occurred at 2.8 GPa experimentally 32 while theoretically these two phases have practically the same enthalpy at p=0 while that of NP rises considerably faster with pressure.These discrepancies between the experimental and theoretical phase boundaries, including the predicted metastability of the HP3 phase at zero temperature, might have several causes.For example, one of the most probable reason for why octahedral FIG. 1. 2 × 2 × 1 supercell of P 3m1 structure with single sheet per period (left pictures), P 63mc with two sheets (middle pictures) and R 3m with three sheets per unit cell (right pictures).Upper pictures show different stackings along the lateral direction parallel to planes, while lower pictures are views from direction axial to the sheets, which corresponds to c-axis vector view in all three phases.Black boxes represent the respective primitive cells with one, two and three SiS2 units for P 3m1, P 63mc and R 3m phases, respectively.Along the axial view, sheets in P 3m1 overlap.geometry was not yet found in laboratory could be the effect of hysteresis that is often exhibited by first-order phase transitions.The transition HP2(HP3)→ P 3m1 is indeed strongly first-order (see later), which implies possible large energy barriers hindering the transformation between the two.In addition to kinetic effects, also ther-mal effect of vibrational entropy might play a role here and perhaps it may explain the fact that normally only the NP phase is observed at ambient conditions.Moreover, inaccuracies of density functional-based methods as well as pressure measurements might also be the source of disagreements.Nevertheless, calculations showed that in the pressure range of 2.8 -5.4 GPa there exist a number of structures with very similar enthalpy and so accurate determination of equilibrium phase boundaries in this region might be complicated from both experimental and theoretical perspective.
More generally, even if this level of agreement between calculations and experiments below 6 GPa may seem somewhat imperfect, it is important to stress that it is on the contrary quite good, because at low pressures all structural enthalpy differences are generally small, and more dependent upon approximations.Both the predictive quality of calculations and the enthalpy differences increase with pressure, and that gives us confidence in our high-pressure study.
Strong structural similarity between octahedral phases implies similar thermodynamical properties.The enthalpy difference between the most stable P 3m1 and the least stable R 3m layered octahedral phase is cca.15 meV/atom at 10 GPa -Fig. 3 (upper panel), which indicates that in experiments one could possibly find SiS 2 in a state of irregularly stacked polytypes.This enthalpy difference is composed of an internal energy difference of less than 5 meV/atom and p∆V term of about 10 meV/atom.The energy differences remain similar at all pressures, probably because both inter-sheet repulsion and strong intraplanar covalent bonding change with pressure in the same way in all phases.The overall enthalpy difference, however, grows with pressure and reaches 73 meV/atom at 100 GPa -Fig. 3 (lower panel).That can be attributed to the pV term which favors phases with more efficient packing, the best being the simplest P 3m1 structure with one sheet per cell (see densities in Table A).
During the evolutionary search, we identified numerous different non-octahedral phases of SiS 2 , however, all of these were found to be grossly metastable at pressures below 100 GPa.Most of these structures are tetrahedral and among them, the HP3 phase was found to possess the lowest enthalpy at lower pressures, yet at 10 GPa it is still about 90 meV/atom higher in enthalpy compared to the three layered octahedral forms and this difference rises quickly with pressure -Fig.3 (lower panel).Some of the identified metastable structures (the quasicubic and monoclinic versions of the HP3-SiS 2 phase and orthorhombic-like version of the rutile phase) were included into the 0-100 GPa enthalpy graph to show strong preference of SiS 2 to form six-coordinated forms in the investigated pressure range.
At 4 GPa, the calculated density of octahedral P 3m1 is 3.077 g.cm −3 , while that calculated for tetrahedral HP2 is 2.530 g.cm −3 and that of HP3 is 2.577 g.cm −3 .The large density jump between P 3m1 and HP2 (HP3) -19.4% (21.6%) shows that despite the layered character, that was not predicted, other expectation of an octahedral phase by Evers et al. 32 are well borne out.The strong first-order character of this structural transition, the large enthalpy gain and density jump suggest some more qualitative but interesting points.First, this transition (or a very similar) is absolutely inevitable in SiS 2 under pressure.Second, it will necessarily involve a large hysteresis, with many possible metastability phenomena en route.Third, once created the new phase will be protected by large free energy barriers.With some qualitative analogy with graphite-diamond (where of course the density jump is a much larger 55% ) these large barriers and the associated nucleation costs might permit to the layered phases to survive in a metastable state, once recovered at ambient pressure and low temperatures.

Electronic structure
At their appearance above 4GPa all new octahedral layered phases are indirect gap semiconductors.The values of GGA band gaps at 10 GPa for the investigated phases, of order 2 eV, are given in Table A. At 30 GPa, the gaps close and all layered forms are predicted to metallize.However, it is a well-known fact that GGA approximation tends to underestimate the value of the band gap and the pressure of metallization and therefore we repeated the electronic band structure calculations for P 3m1 also using the HSE06 hybrid functional that includes exact exchange.Within HSE06, probably more reliable in this respect, all three layered octahedral phases metallize around 40 GPa.In Fig. 4, the HSE06 electronic band structure of P 3m1 is shown for the semiconducting state at 10 GPa (upper panel) and at the pressure of metallization at 40 GPa (middle panel), along the cor-  FIG. 6. Phonon dispersion curves for the P 3m1 structure at zero pressure along the selected high-symmetry lines calculated with the phonopy code 52 .The structure shows no imaginary frequencies and might be therefore quenchable to ambient conditions.

Dynamical and elastic stability
By calculating phonon dispersions at 0, 10, 30, 60 GPa and at 100 GPa, we found that the P 3m1 phase shows no imaginary frequencies from ambient pressure up to 1 Mbar, indicating its dynamical stability, at least at zero temperature.In Fig. 6, phonon dispersion relations of the P 3m1 phase are shown for zero pressure.Modes for the two other phases look, as one would expect, very similar.
Elastic constants C αβ (in Voigt notation) listed in Table A were calculated at fixed volumes corresponding to zero pressure and at T = 0 K.The Table also shows point groups of our investigated phases and their corresponding Laue groups as well as the respective symmetry relations between elastic constants and elastic stability conditions.For all structures there is about a factor of eight difference in the value of C 33 (corresponding to the direction perpendicular to planes) with respect to C 11 = C 22 .This indicates highly anisotropic mechanical properties as expected for layered systems.All three investigated structures were found to be mechanically stable at zero pressure by fulfilling the necessary and sufficient Born stability conditions for their corresponding Laue groups 53 , which come from the requirement of positive definiteness of the C αβ matrix.Therefore, recovery of these phases in a metastable state at ambient pressure is at least conceptually possible.

Monolayer SiS2
The dynamical and elastic stability of the P 3m1 layered phase at all pressures down to zero has several consequences.First, it should be possible to shear the layers relatively easily, as one does in graphite and in MoS 2 .
Again by similarity, one can hypothesize the possibility to exfoliate a single octahedral SiS 2 monolayer from bulk P 3m1 and recover it at ambient conditions.To pursue this idea, we carried out a separate study of an SiS 2 monolayer.After geometrical optimization in a cell with 15 Å vacuum between monolayers, we calculated its electronic structure along k z =0 path -Fig.7. The band structure shows that the monolayer is again an indirect semiconductor with gap of about 2.4 eV.In addition, we also calculated phonon dispersions of the monolayer and, as in the case of 3D layered octahedral phases, we found no imaginary modes indicating its dynamical and mechanical stability in isolation.taken as one layer from HP1-SiS 2 has an indirect band gap larger than 3 eV.

Molecular Dynamics of non-layered-to-layered structural transformation
As in all first-order transitions, the transformation from non-layered to layered structures will proceed by nucleation.While nucleation is in itself extremely slow and hard to simulate, much can still be learned about the local mechanism of transformation by means of variable cell constant-pressure ab initio molecular dynamics 54 which artificially permits the transformation to take place uniformly and in extremely short times.
We therefore conducted N P T ab initio MD simulations starting from the initial structure of HP3-SiS 2 for 48 and for 96 atoms at low pressure and room temperature and continuing by increasing pressure and temperature.In both sizes we observed a chemical reorganization into fully octahedral forms -at the pressure of 60 GPa and temperature 600 K for the 96-atomic system and at 50 GPa and 500 K for the 48-atomic system.The initial transformations led to irregular connectivity pattern of the SiS 6 units, which shared edges and corners and also contained unpaired S atoms.Upon further compression of the 48-atomic system to 60 GPa at 500 K, the octahedra shifted and rotated, transforming into an incompletely layered state with numerous defects.Nevertheless, a strong resemblance to the P 6 3 mc phase that contains two sheets per unit cell is visible from the polyhedral view of P 6 3 mc structure and MD result at 60 GPa -Fig.8.
MD calculations therefore confirm the existence of octahedral SiS 2 with layered character at high pressures, even though due to the limited time scales of ab initio MD simulations, proper defectless crystalline structures were not fully recovered.An additional reason for incomplete transition observed in MD in this specific case might be the fact that the HP3 phase, which should transform into P 3m1, is not layered, but fully three-dimensional.The transformation into octahedral geometry starts at random places in the (simulation) sample and the emerging octahedral arrangement of silicon and sulfur atoms formed within the initial chemical reaction then remains mostly unchanged in further evolution at the picosecond time scale.Only a substantial diffusion of certain Si atoms could create proper two-dimensional sheets from irregularly connected octahedra, and therefore the original disordered polyhedral state remains stable on the accessible ab initio time scale.
From the observed results, it can be predicted that the transformation from HP3 (or HP2) to P 3m1 phase possibly proceeds in two stages -chemical and topological.The first, chemical stage is the transformation of distorted tetrahedra of HP3 into octahedra, which lasted only about 2 ps in our simulations.In the second, topological stage, octahedra must shift and rotate in order to properly organize into individual layers, but during this process some atoms must diffuse over distances that are larger than distances of second-nearest neighbors, which in a solid material takes considerably longer time than the initial chemical reaction.
These results and conclusions are important in two respects.First, they probably explain why HP3 may have been created and survived experimentally even though metastable with respect to layered P 3m1.In turn, they support the possibility that layered P 3m1, once created, might survive as a metastable phase once brought down to lower or even zero pressure.

Structural similarities with isoelectronic compounds
Finally, we should comment on structural similarities among the isovalent compounds CO 2 , CS 2 , SiO 2 and SiS 2 .With hindsight, before the advent of efficient crystal structure prediction techniques it was customary to argue based on assumption of similarity of the tetrahedral structures, in particular to those of well-known SiO 2 [15][16][17][18][19][20][21][22] .
With the current knowledge of phase diagrams many such arguments no longer stand.In Fig. 9 we show the diagrams of the T =0 stable structures of the four compounds in a wide range of pressures up to 200 GPa.Restricting to these stable structures there is in reality only a limited amount of similarity.The only tetrahedral structures that are stable in more than one compound are the monoclinic P 2 1 /c (whose existence in CS 2 is predicted, even not yet experimentally demonstrated) and tetragonal I 42d (if it is indeed stable and not metastable in SiS 2 and stable in CS 2 , which also remains to be experimen-tally demonstrated).Considering also metastable structures, the only tetrahedral structure which seems to be universal at least topologically is the β-cristobalite structure which is found in its ideal cubic form as a metastable phase of SiO 2 and in its distorted form I 42d stable in CO 2 and stable or metastable in SiS 2 (possibly stable in CS 2 too).Surprisingly, the lowest amount of structural similarity to other compounds is found in SiO 2 while the highest is seen in SiS 2 and CS 2 .In particular, the wellknown quartz structure, which in its α-form represents the stable form of SiO 2 at ambient pressure, does not seem to be realized in stable form anywhere else; and the same is true for the coesite structure.This comparison clearly shows that simple analogies are of limited value if used as heuristic guide in the family of binary compounds that we considered, and probably also elsewhere.

DISCUSSION
In summary, we predicted by a highly reliable theoretical protocol three new low-enthalpy high-pressure phases of SiS 2 , namely P 3m1, P 6 3 mc and R 3m by means of ab initio structure searching method based on evolutionary algorithm.These three novel phases are all formed by SiS 6 octahedral units forming separate sheets and differ only by the stacking pattern of the layers.All phases are predicted to become more stable than the tetrahedral forms beyond 4 GPa at zero temperature and remain to be the three most stable structures up to at least 100 GPa, out of all structures found.Upon increase of pressure SiS 2 therefore undergoes an interesting evolution as far as the dimensionality is concerned -tetrahedral structures transform from 1D to 2D to 3D and after transition to octahedral structure it settles on a 2D layered character.The associated density increase for tetrahedral-tooctahedral transformation is predicted to be around 20% at 4 GPa showing a strong first-order nature of the transition.The associated high barriers, and the large reorganization demanded by the transformation suggests the possibility to recover the new phases in metastable form even at ambient conditions.Within the employed calculation scheme, the HP3 phase appears to be metastable at zero temperature.Enthalpies of the layered octahedral phases are very similar at lower pressures, which indicates a possible presence of a layered structure with irregular stacking of sheets at experimental conditions.In constant-pressure molecular dynamics simulations, the spontaneous appearance of octahedral geometry start-FIG.9. Comparison of zero-temperature structures of CO2 1-14 , CS2 [28][29][30][31] , SiO2 [15][16][17][18][19][20][21][22] and SiS2 [32][33][34][35][36] showing a rather weak level of similarity.The two structures P 21/c and I 42d that are common for some systems are highlighted in the diagram and shown for the case of SiS2.For CO2, phase V is usually created during heating only at pressures over 50 GPa though it is predicted to be thermodynamically stable from 20 GPa 10-12 .
ing from an initial tetrahedral one was also observed at 50 GPa, even though the resulting octahedral structure remained mostly disorganized due to short simulation times accessible in first-principles studies.All three lowenthalpy octahedral polymorphs are semiconductors with indirect 2 eV order band gaps at low pressures, values potentially interesting for layered compound applications.They metallize around 40 GPa and remain relatively poor metals up to 100 GPa.The P 3m1 phase is dynamically stable from zero pressure up to at least 1 Mbar from phonon calculations.Elastic stability of all phases was also confirmed at zero pressure and the calculated elastic constants indicate highly anisotropic mechanical properties typical for layered structures.Calculations also suggest that single octahedral layers could possibly be exfoliated from bulk P 3m1-SiS 2 once quenched to ambient conditions, as far as the monolayer is predicted to be dynamically stable once formed.Results presented in this study may be useful in order to compare the future experimental compression products with the ideal phases found in our structural search.On account of their layered nature, the application of high-pressure shear may suffice to cause interlayer sliding, with transformations and the appearance of new polytype phases.Comparing the T =0 phase diagrams of selected group IV-VI AB 2 systems shows that structural similarities between these isovalent compounds are in fact rather limited.

METHODS
We employed the open-source evolutionary algorithm package for crystal structure prediction XtalOpt 42 and carried out searches at 10, 30, 60 and at 100 GPa with 6 and 12 atom cells (2 and 4 structural units) for each pressure, enabling generation of at least 1000 structures in each case.All underlying ab initio electronic calculations were performed with the density functional theory (DFT) VASP 5.3 and 5.4 codes 55,56 employing projector augmented-wave pseudopotentials 57,58 and Perdew-Burke-Ernzerhof (PBE) parametrization of GGA exchange-correlation functional 59 .Pressures above 100 GPa could not be, and were not, addressed with these pseudopotentials.The randomly generated initial structures of the evolutionary search and those created after various mutations were initially optimized in static cells (where only ions were relaxed) and within coarse kpoint grids.Thereafter, structures were relaxed in fixedvolume geometry (ions and cell shapes were relaxed) with a finer k-points mesh.Finally, full ionic and cell degrees of freedom relaxations were carried out with finest k-point sampling in order to refine the final structures.During this progressive structural optimization, the en-ergy cutoff was also gradually increased from 340 to 400 eV.
In order to describe the dispersion forces that must be included for the obtained layered phases, especially in the low-pressure region, we employed the parameterfree Tkatchenko-Scheffler approach DFT-TS 60 , in which dispersion coefficients and damping functions of the original Grimme's method 61,62 are charge-density dependent.The dispersion forces are therefore described adaptively accounting for changing electronic structure upon compression.To investigate electronic structure of octahedral SiS 6 , we employed the HSE06 hybrid functional 63 , which tends to provide more reliable results for the band gap estimate and electronic structure of solids compared to PBE.
The second-order elastic stiffness constants C αβ were calculated using VASP by a two-steps procedure, in which the stress tensor is determined for both distorted supercells without relaxing the ions and for undistorted supercells with separately shifted ions 64 .For each phase, we used 21 different values of strain -from 0.005 to 0.025 in 0.001 intervals and from these, only results that were quite similar (formed plateaus on C αβ vs. strain dependences) were used to obtain final average -this was from strain values 0.015 to 0.02 for P 3m1 phase, between 0.011-0.021for the case of P 6 3 mc and for 0.013-0.016in-terval for R 3m.From the final elastic tensor of each phase, individual C αβ constants were extracted by requirements of symmetry according to each of the phases' Laue groups 65 .For example, if, say, C 11 should by symmetry be equal to C 22 , then C 11 was taken to be (C 11 + C 22 )/2, and so on.Elastic constants for all phases were calculated at zero temperature at volumes corresponding to zero pressure and with the use of Grimme's van der Waals correction (DFT-D2) 61,62 .
Phonon dispersion curves were calculated in the harmonic approximation by the supercell method using the phonopy code 52 .A 6 × 6 × 4 supercell (432 atoms) for the P 3m1 phase, a 6 × 6 × 2 supercell (432 atoms) for the P 6 3 mc phase and a 6 × 6 × 2 supercell (648 atoms) for R 3m phase were used for calculations of force constant matrices with phonopy.Correctness of the dispersion curves was also tested in some important cases by independent calculations using density functional perturbation theory as employed in the Quantum Espresso package 66 .
MD simulations were performed in the constant pressure-temperature NPT ensemble with the use of Parrinello-Rahman barostat 54 and Langevin stochastic thermostat 67 .The simulation samples contained 48 and 96 atoms and the supercells were generated as 4 × 4 × 4 and 4 × 4 × 8 unit-cells of the HP3-SiS 2 structure.The time step for MD was set to 2 fs.

FIG. 2 .
FIG. 2. The most stable and the most dense P 3m1 phase (C6 structure) of SiS2 with one sheet per unit cell.Polyhedral view visualized by VESTA 43 shows interbonded edge-sharing octahedra forming individual sheets.
FIG. 3. (Upper panel)Enthalpies of octahedral structures P 3m1 (black solid curve), P 63mc (red dashed curve) and R 3m (green dotted curve) and of tetrahedral NP (brown dashed-double dotted), HP1 (violet dashed-dotted) and HP2 (orange dashed) forms relative to HP3 (blue solid horizontal line) calculated up to 10 GPa.Layered octahedral phases become more stable than tetrahedral forms over 4 GPa.(Lower panel) Enthalpies relative to P 3m1 of the three octahedral phases, HP3 and of three newly identified metastable forms of SiS2 (orthorhombic version of rutile structure, quasi-cubic and monoclinic versions of HP3) in 1 Mbar range showing a strong tendency towards forming octahedral forms.

FIG. 4 .
FIG. 4. HSE06 electronic band structure of P 3m1 in the semiconducting state at 10 GPa (upper panel) and at 40 GPa (middle panel), where the system begins to metallize.In the lower panel, the Γ → K → M → Γ → A → H → L → A path in the first Brillouin zone is visualized.

FIG. 7 .
FIG. 7. Electronic band structure of a single octahedral monolayer calculated with HSE06 functional.

= 190. 9 ;FIG. 8 .
FIG. 8. Comparison of the ideal P 63mc structure and result of MD at 60 GPa and 500 K. Planes arising from irregularly connected octahedra are well-visible from the dynamical simulations.