Structural evolution and molecular dissociation of H2S under high pressures

Wen-Ji Shen(沈文吉)， Tian-Xiao Liang(梁天笑)， Zhao Liu(劉召)， Xin Wang(王鑫)，De-Fang Duan(段德芳)， Hong-Yu Yu(于洪雨)，?， and Tian Cui(崔田)2，，?

1College of Physics，Jilin University，Changchun 130012，China

2Institute of High Pressure Physics，School of Physical Science and Technology，Ningbo University，Ningbo 315211，China

Keywords: high pressure，crystal structure，phase transition，superconductivity

1. Introduction

In our work， we explore the phase transition of H2S viaab initiostructure searching andab initiomolecular dynamics (AIMD) to complement its important phase transition sequence. Three new structures of H–S compounds are discovered for the first time， and the relevant thermodynamic properties and the superconductivity properties are revealed. Remarkably，we also investigate the dissociation process of H2S at high pressure， and our calculation results can help to further clarify the phase diagrams of H2S at higher pressures and temperatures.

2. Computational details

Structure of H2S were investigated withab initiocalculations as implemented in the Universal Structure Predictor:Evolutionary Xtallography(USPEX)[37–39]in a pressure range of 0 GPa–30 GPa in pressure steps of 5 GPa. The Viennaab initiosimulation packages (VASP) code[40]was used to optimize crystal structures and calculate the electronic properties， where the Perdew–Burke–Ernzerhof (PBE)[41]of generalized gradient approximation(GGA)[42]was chosen and used as an exchange–correlation functional. And the projectoraugmented wave method(PAW)[43]potentials with the 1s1and 3s23p4configurations treated as valence electrons for H and S were used. Kinetic cutoff energy of 400 eV was selected and Monkhorst–Packk-mesh with grid spacing of 2π×0.03 ?A-1was then adopted to ensure that the enthalpy converges to less than 1 meV/atom. The calculations on phonons were carried out by the PHONOPY code，[44]which showed good agreement with those computed using the Quantum-ESPRESSO(QE) code.[45]Electron phonon coupling (EPC) matrix elements were calculated using density functional perturbation theory as implemented in the Quantum ESPRESSO package. Ultra-soft pseudopotential[46]was used with an 80 Ry(1 Ry=13.6056923(12)eV)kinetic energy cutoff and an 800-Ry charge density cutoff. A 12×12×20k-mesh grid was used for Brillouin zone(BZ)integration in electronic calculations， and EPC matrix elements were calculated in 3×3×5q-points grid. Electronic smearing in self-consistent calculation adopted Methfessel–Paxton first-order spreading[47]with a smearing width of 0.025 Ry. Furthermore， we also performed extensiveab initiomolecular dynamics(AIMD)simulations based on the Born–Oppenheimer approximation implemented in VASP，and set the temperature range to be from 100 K to 250 K. The AIMD simulation for H2S phases were performed in 2×2×2 supercells (384 atoms) forP42/nandI41/a，2×2×3 supercells(144 atoms)forP2/c，and 2×3×3 supercells (108 atoms) forP-1 with 2×2×2k-mesh. We adopted theNPT(N-constant is the number of particles，Pconstant is the pressure， andT-constant is the temperature)ensemble， lasting for 12 ps in time steps of 1 fs， and we allowed 2-ps thermalization and then extracted data from the last 10 ps.

3. Results and discussion

Here， the global structural evolution algorithm variablecomposition predictions are realized in over 2×104structures，and 2， 4， 6， 8， and 16 H2S formula units are employed in a pressure range from 0 GPa to 30 GPa. Two new phases of H2S with space groups ofP42/nandI41/aare uncovered as shown in Fig. 1. The H atoms occupy 8f sites， and S atoms 8f sites in theP42/nlattice. And H atoms account for 16f sites， and S atoms 16f sites in theI41/alattice， as shown in Table S1 in the Supplementary materials (SM). As shown in Figs.1(a)and 1(b)，P42/nandI41/acrystal lattice contain isolated H2S units inside，and we can observe from the structure intuitively thatP42/nandI41/aare molecular crystals. The enthalpies of two new structures together with the ambientpressurePbcm(phase III)and another four low-enthalpy structures mentioned by Liet al.[22]are plotted as a function of pressure in Fig.2. With pressure increasing，Pbcm(phase III)of H2S crystallize into aP42/nphase at about 3 GPa. Above 10 GPa，P42/nphase turns intoP2/cphase，[22]and then transforms intoI41/aat about 20 GPa.

Fig.1. Predicted energetically stable crystal structures of(a)P42/n and(b)I41/a in structural evolution.

Fig. 2. Curves of calculated enthalpy versus pressure for phases P21m，P42/n， I41/a， Pbcm (phase III)，[21] P-1，[22] Pmc21，[22] and Pc[22] relative to P2/c.[22]

Phonon dispersion calculations ofP42/nandI41/astructures in Figs.3(a)and 3(b)do not give any imaginary frequencies. Therefore，P42/nandI41/aare dynamically stable at corresponding pressure points. The electronic band structures of H2S are calculated at 10 GPa forP42/nphase and at 30 GPa forI41/aphase as shown in Figs. 3(c) and 3(d). The results of electrical properties show that the new phases are insulators because band gaps between valence band and conduction band are quite large，specifically，3.3 eV forP42/nand 2.5 eV forI41/a. The p-electrons of S dominate the density of states nearby the fermi surface for bothP42/nandI41/a，and dominates the density of states of valence bands significantly.

Fig.3. Calculated phonon dispersion curves for(a)P42/n at 10 GPa and(b)I41/a at 30 GPa，and corresponding electronic band structure and projected density of states for(c)P42/n at 10 GPa and(d)I41/a at 30 GPa.

To further study the dynamical properties of phasesP42/nandI41/a，ab initiomolecular dynamics(AIMD)simulations within finite temperature are performed. By calculating the radial distribution function(RDF)gSHof the averaged S-H positions， we find that crystals ofP42/nandI41/ado not change obviously with temperature rising(see Figs.S1(d)–S1(f)and Figs.S2(d)–S2(f)in the SM).To sum up the results discussed above， the final phase transition sequence at zero temperature isPbcm(phase III)→P42/n →P2/c →I41/a →Pcbelow 30 GPa.

Fig.4. Snapshot of MD trajectory at 250 K and(a)15 GPa，(b)75 GPa，and(c)155 GPa. (d)New P21/m phase of H2S.

Considering the energetically favored candidate for phase IV is theP2/c，whose x-ray diffraction patterns can reproduce the main peaks of phase IV below 2θof 20°， as well as the weak peaks around 25°. We perform the AIMD simulations in 2×2×3 supercells forP2/c-H2S at 100， 200， and 250 K in a pressure range of 15 GPa–175 GPa to study the molecular dissociation and superconducting phase of H2S.By observing the snapshots of MD trajectory at 250 K(See Figs.4(a)–4(c))，we can clearly observe the dissociation process ofP2/c-H2S and the reconstruction process of a new polymeric structure.By statistically averaging the trajectories of S and H atoms at 155 GPa and 250 K，we obtain the symmetry and atomic coordinates of this new structure(see Table 1).

To further examine the dissociation and reconstruction process of H2S，the values of RDFgSHare obtained from the AIMD simulations at 250 K and different pressures， and the results are shown in Fig.5. When the pressure increases up to 45 GPa，the first peak ofgSHremains almost unchanged but the second peak shifts left forward obviously for each curve. With the snapshots of MD trajectory，we can see that the H2S keepsP2/csymmetry but the distances between S and H atoms in different H2S molecules keep closer to each other as the pressure increases. From 55 GPa to 95 GPa， the first peak shifts slightly from 1.36 ?A to 1.37 ?A，and the first peak and the second peak tend to merge together (see Fig. 5 and Fig. S3 in the SM). It is due to the nearest S and H atoms in different H2S molecules get closer enough and the interactions between these S and H atoms get stronger. As the pressure increases，some protons occupy two symmetric positions along the S–S bond with equal probability gradually， just like the process of symmetric hydrogen bond of ice phase VII，[48]meanwhile some H2S molecules start to decompose into HS+H3S. At even higher pressure，we can see a shoulder clearly in thegSHat 115 GPa， indicating that H2S molecules start to be reconstructed. Especially when pressure increases from 115 GPa to 135 GPa， the bond length of H–S gradually increases to 1.43 ?A，forming a new structure.

Table 1. Predicted structural parameters of H2S P21/m.

Fig.5. RDFs of averaged S–H positions under different pressures of H2S–P2c at 250 K.

In the same pressure interval， combining thegSHand AIMD snapshots ofP2/cphase at 100 K and 200 K， we can see that the symmetric hydrogen bonds appear at 75 GPa and 85 GPa，respectively. More seriously，at 100 K and 200 K(see Figs. S4 and S5 in the SM)， although partial H2S molecules are dissociated at higher pressure，the reconstruction is not obvious until 155 GPa. In conclusion， the AIMD results indicate that the newP21/mstructure of H2S can be derived from dissociation ofP2/cphase， indicating the superconductivity at 160 GPa. In our further investigation， the superconducting transition temperatureTCforP21/mat 160 GPa was estimated through the McMillan equation:

whereλis the EPC constant，ωlogis an appropriately defined logarithmic average frequency，μ*is the Coulomb pseudopotential，which is a parameter accounting for the effective Coulomb repulsionn generally takesμ*=0.1 for hydride，λis the EPC parameter，and logarithmic average frequencyωlogis the logarithmic average frequency:λandωlogare calculated from

The superconducting transition temperatureTCforP21/mis estimatedat 82.5 K under 160 GPa， obtainingλ=1.93 andωlog=588.1 K correspondingly.

As shown in Fig. 2， at 80 GPa–160 GPa，P-1 phase is more energetically stable than others. Hence， we also perform the AIMD simulations in 2×3×3 supercells forP-1-H2S at 100 K， 200 K， and 250 K in a pressure range of 90 GPa–200 GPa. TheP-1 phase will finally turn into theCmcaphase[22]with pressure increasing. From the snapshots of MD trajectory at 250 K(see Figs.6(a)–6(c))，we can see that theP-1 phase remains stable in a pressure range of 90 GPa–120 GPa. TheP-1 phase contains the H3S–SH3units which are connected mutually through planar symmetric hydrogen bonds of two H atoms. Under further compression， theP-1 phase eventually transforms into theCmcaphase，and each S atom in theCmcastructure bonds with two H atoms in adjacent two layers，and all H atoms bond with two S atoms，which fits well with our simulations.

Furthermore， the RDFgSH， are also obtained from the AIMD simulations at 100 K， 200 K， and 250 K， and the results are shown in Fig.7. With pressure increasing at 250 K，the first peak at 1.4 ?A and the second peak at 1.9 ?A simultaneously shift rightward and leftward，for each curve，respectively.When the pressure increases to 140 GPa，the second peak becomes a shoulder at about 1.7 ?A and completely merge together with the first peak above 160 GPa. This is similar to the dissociation and reconstruction process ofP2/c，and the same results are obtained at 200 K and 100 K (see Figs. 6(d) and 6(e)). And as the temperature decreases， the pressure of the phase transition fromP-1 toCmcaincreases. With the RDFgSH，as shown in Fig.7，we can see that the transition pressure at 250 K，200 K，and 100 K are about 140 GPa，160 GPa，and 180 GPa，respectively.

Fig.6. Snapshot of MD trajectory at(a)90 GPa and 250 K，(b)120 GPa and 250 K，(c)160 GPa and 250 K，(d)180 GPa and 100 K，and(e)180 GPa and 200 K.

Fig.7. Dynamical behaviors of H and S atoms in S–H compounds，obtained from AIMD simulations for P-1 phase in solid states，for RDFs of averaged S–H positions at(a)100 K，(b)200 K，and(c)250 K.

4. Conclusions

In this work， we predict two new structures of H2S withP42/nandI41/ain a lower pressure range (0 GPa–30 GPa).TheP42/nphase and theI41/aphase are dynamically stable and insulators. The final phase transition sequence at zero temperature and high pressure is fromPbcm(phase III) toP42/nat 3 GPa，toP2/c(phase IV)at 10 GPa，then toI41/aat 20 GPa， and finally toPcat 40 GPa. Our calculations for thermodynamic，electronic，and dynamic properties reveal this newly discovered phase transition.We also preform the AIMD simulation forP2/cof H2S in different pressure and temperature ranges. At 250 K，the radial distribution function calculations obviously show the symmetric hydrogen bond and the reconstruction of H2S as the pressure increases. Besides， we find a newP21/mphase at 160 GPa and 250 K. TheP21/mphase may be a potential superconductor withTCof 82.5 K at 160 GPa. This study may be helpful in studying the lowtemperature superconducting phase of H2S. And this result may be conducive to further clarifying the phase diagram of H2S and other H–S compounds under higher pressure and temperature.

Acknowledgements

The computation of this work was carried out at the High-Performance Computing Center (HPCC) of Jilin University，Beijing Super Cloud Computing Center， and TianHe-1(A) at the National Supercomputer Center in Tianjin.

Project supported by the National Natural Science Foundation of China(Grant Nos.11704143，11804113，11604023，and 12122405).