Novel high pressure structures and superconductivity of niobium disulfide
NNovel high pressure structures and superconductivity of niobium disulfide
Zhong-Li Liu,
1, 2, ∗ Ling-Cang Cai, and Xiu-Lu Zhang College of Physics and Electric Information, Luoyang Normal University, Luoyang 471022, China Laboratory for Shock Wave and Detonation Physics Research,Institute of Fluid Physics, P.O. Box 919-102, 621900 Mianyang, Sichuan, China Laboratory for Extreme Conditions Matter Properties,Southwest University of Science and Technology, 621010 Mianyang, Sichuan, China (Dated: January 16, 2014)We have investigated the pressure-induced phase transition and superconducting properties ofniobium disulfide (NbS ) based on the density functional theory. The structures of NbS at pres-sures from 0 to 200 GPa were predicted using the multi-algorithm collaborative (MAC) structureprediction technique. The previously known 1 T -, 2 H -, and 3 R -NbS were successfully reproduced.In addition, many metastable structures which are potential to be synthesized were also discovered.Based on the enthalpy calculations, we found that at 26 GPa NbS transits from the double-hexagonal (2 H ) structure to the tetragonal I /mmm structure with a 10.6% volume reduction.The calculated elastic constants and phonon dispersion curves of I /mmm -NbS confirm its me-chanical and dynamical stability at high pressure. More interestingly, the coordination number ofNb in I /mmm structure is eight which is larger than that in the traditional metal dichalcogenides,indicating a new type of bondings of Nb and S atoms. In the new Nb-S bondings, one Nb atomand neighboring eight S atoms form a [NbS ] hexahedron unit. Furthermore, I /mmm -NbS ex-hibits a higher superconducting critical temperature than 2 H -NbS , as is resulted from the strongerelectron-phonon coupling coefficients. I. INTRODUCTION
Transition metal dichalcogenides (TMDs) MX (M =Nb, Ta, Mo, W, X = S, Se, Te) have intriguing proper-ties, ranging from insulator to metal and superconduc-tor, and thus always attract extensive interests of ex-perimentalists and theorists. Thanks to their in-planecovalent bondings and weak interlayer van der Waalsinteractions, they could be easily exfoliated down to amonolayer which shows very exotic properties. For exam-ple, bulk MoS is an indirect-band-gap semiconductor, while the monolayer MoS is a direct-band-gap semicon-ductor. Consequently, the TMDs have shown excitingprospects for a variety of applications, such as catalystsand lubricants in the petroleum industry , promising ap-plications in nanoelectronics and optoelectronics, andenergy storage applications. The typical representative of TMDs is 2 H -NbSe ,showing a large charge density wave (CDW) (at 33 K)that coexists with superconductivity ( T c = 7.2 K). Nio-bium disulfide (NbS ) also belongs to the family of TMDcompounds. But the CDW order appeared in 2 H -NbSe is absent in 2 H -NbS , and its occurrence is sup-pressed by the large anharmonic effects. However, italso shows superconductivity at a similar transition tem-perature of T c = 6 K. More interestingly, the T c of2 H -NbS increases smoothly from 6 K at zero pressure to ∼ also similar to the behavior of T c in2 H -NbSe which increases to ∼ Theupper critical field of 2 H -NbS has an initial decreaseas pressure increases, contrary to the increase of T c , butabove 8.7 GPa it increases again with pressure. NbS is a two-gap superconductor, similar to NbSe .The heat capacity of a 2 H -NbS has been measured downto 0.6 K and in magnetic fields up to 14 T by Kaˇcmarˇcik et al . The temperature dependence of the electronic spe-cific heat can be attributed to either the existence of astrongly anisotropic single-energy gap or a two-gap sce-nario with the large gap about twice bigger than the smallone. The field dependence of the Sommerfeld coefficientinduces a magnetic field dependence of the superconduct-ing anisotropy. The two-gap scenario conclusions aresupported by the absence of in-plane gap anisotropy inrecent STM imaging of the vortex lattice in NbS . H -NbS has a layered structure and therefore haslarge anisotropic electrical, optical, and magnetic prop-erties. It has been applied as catalyst for the purifi-cation of petroleum, cathode materials in secondarybatteries, humidity sensors, and so on. In exper-iment, presently the one-layer trigonal 1 T -NbS , two-layer hexagonal 2 H -NbS , and three-layer rhombohe-dral 3 R -NbS heve been synthesized. Large-scalesynthesis of 3 R -NbS nanosheets has also been recentlyrealized. Different low dimensional structures of NbS have different physical and chemical properties. Lowdimensional materials depend on and can be exfoliatedfrom bulk materials. It is necessary to uncover as manycrystal structures of NbS as possible. From some newcrystals, it is expected to exfoliate some new-functionallow dimensional materials.It is known that pressure is able to modulate theproperties of materials through changing their crystalstructures. Furthermore, the structures of NbS un-der high pressure are fundamental to understand its su- a r X i v : . [ c ond - m a t . s up r- c on ] J a n perconductive properties. The mechanism of pressure-induced superconductivity and the superconducting tem-perature in NbS above 20 GPa still remain unknown tous. This motivates us to investigate the superconductiv-ity of NbS at higher pressures. In this work, we firstpredicted the high-pressure structures of NbS and de-termined its phase transition sequence using the multi-algorithm collaborative (MAC) crystal structure predic-tion technique combined with the density functional the-ory (DFT). Then we calculated the superconducting crit-ical temperature through electron-phonon coupling cal-culations.The paper is organized as follows. Section II containsthe computational details. The results and discussion arepresented in Sec. III. Conclusions follow in Sec. IV. II. COMPUTATIONAL DETAILS
In order to determine the high-pressure structures ofNbS , we searched its low-energy structures from 0 to200 GPa using our developed MAC crystal structure pre-diction technique. The multi algorithms including theevolutionary, the simulated annealing, and the basin hop-ping algorithms are combined to collaboratively searchthe global energy minima of materials with the fixed sto-ichiometry. The MAC algorithm and all the relevanttechniques are incorporated in the
Muse code. The re-sults were also carefully cross checked and confirmed bythe
Calypso code, which is based on the particleswarm optimization algorithm.The ab initio optimizations for every structure gen-erated by the
Muse code were performed with vasp package.
We tested the local density approxima-tion (LDA) and the generalized gradient approximation(GGA) parametrized by Perdew, Burke and Ernzerhof(PBE) for exchange correlation energy. The two ap-proximations give the similar structures order in struc-ture prediction. While the LDA calculated lattice con-stants are better than the GGA for NbS . So in the staticcalculations, we adopted the LDA exchange-correlationfunctional. The electron-ion interactions are described bythe projector augmented wave (PAW) scheme. Thepseudopotentials for Nb and S have the valence electrons’configurations of 4p and 3s , respectively. Toachieve good convergences the kinetic energy cutoff andthe k -point grids spacing were chosen to be 500 eV and0.02 ˚A − , respectively. The accuracies of the target pres-sure and the energy convergence for all optimizations arebetter than 0.1 GPa and 10 − eV, respectively. Thesearched systems contained 6, 9, 12, 15 and 18 atomsin the unit cell. TABLE I: The comparison of the calculated lattice constants a , c , and c/a of different NbS structures with the correspond-ing experimental values.Structure Method a (˚A) c (˚A) c/a P (GPa) Reference1 T -NbS VASP-LDA 3.253 5.341 1.642 0 This workExperiment 3.420 5.938 1.736 0 182 H -NbS VASP-LDA 3.287 11.421 3.475 0 This workExperiment 3.418 11.860 3.470 0 18Experiment 3.310 11.890 3.592 0 19Experiment 3.330 11.950 3.589 0 103 R -NbS VASP-LDA 3.286 17.577 5.349 0 This workExperiment 3.335 17.834 5.336 0 17
III. RESULTS AND DISCUSSIONA. Structure prediction for NbS at high pressures In our MAC structure searches, the structures weregenerated with symmetry constraints in the first gener-ation and optimized with vasp package at fixed pres-sures. The pressures applied to crystal structures inoptimizations go from 0 to 200 GPa with the intervalof 20 GPa. At each fixed pressure, the enthalpies ofthese structures were calculated and compared to findthe proper path towards the lowest-enthalpy structure.Results show that the previously known 2 H -NbS hasthe lowest enthalpy at lower pressures (below 20 GPa)and the I /mmm structure has the lowest enthalpy athigher pressures (above 20 GPa). The previously known1 T - and 3 R -NbS structures were also easily reproduced.The large-scale 3 R -NbS nanosheets are synthesized veryrecently. More interestingly, we found a new two-layerhexagonal structure whose energy is very close to 3 R -NbS at ambient pressure. We refer to it as 2 H (cid:48) -NbS .According to energy criterion, it is potential to be syn-thesized in experiment. Meanwhile, many other struc-tures were found to be energetically competitive, includ-ing P m P P
22 structures, and so on. Amongthese structures, the trigonal P m R -NbS in thewhole pressure range of interest. So, it is also expectedto be synthesized in experiment.The calculated lattice constants of 1 T -, 2 H -, and 3 R -NbS are listed in Table. I, in comparison with experi-mental values. We note that the lattice constants a and c of the three structures are all slightly under-estimated in our LDA calculations. But the calculated c/a values are all in good agreement with experiments.To further examine the three structures, we also simu-lated their X-ray diffraction (XRD) patterns and com-pare them with experimental data. The calculated XRDpatterns of the three structures are all in excellent agree-ment with corresponding experiments (Fig. 1), indicatingthat each structure is identical to the known one. (b) P / mmc Expt.: 2H-NbS I n t e n s it y ( a r b . un it ) (c) R m Expt.: 3R-NbS
2θ ( ○ )10 20 30 40 50 60 70 (a) Expt.: 1T-NbS P m FIG. 1: (color online). Simulated XRD patterns of 2 H -NbS ,3 R -NbS , and 1 T -NbS , in comparison with the correspond-ing experimental results (1 T : 18, 2 H : 18, 3 R : 17). The new 2 H (cid:48) -NbS crystal has the 2 H -MoS structureand can be formed by shifting one layer of atoms in 2 H -NbS . The shifting distance is 0.577 lattice constant a along typical direction. 2 H (cid:48) -NbS has six atoms in prim-itive cell with the lattice constants of 3.28, 3.28 and 11.65˚A at ambient pressure. The Nb and S atoms are at Wyck-off’s 2c positions (1/3, 2/3, 1/4) and 4f positions (1/3,2/3, 0.62), respectively. We show the 2 H (cid:48) -NbS struc-ture and the shifting direction in Fig. 2. The shiftingdirection is parallel to the layer plane (Fig. 2 f). Thatis to say that the structures of the two layers are thesame. The unique difference between 2 H (cid:48) - and 2 H -NbS is the relative positions of the two layers. The coordi-nation numbers of Nb atoms in both 2 H (cid:48) - and 2 H -NbS are six. One Nb atom and the neighboring six S atomsform a [NbS ] trigonal prismoid. Accordingly, the coor-dination numbers of S atoms in both 2 H (cid:48) - and 2 H -NbS are three.The new I /mmm structure is plotted in Fig. 3. It hassix atoms in conventional unit cell (three in primitive cell)with the lattice constants of 3.15, 3.15 and 7.91 ˚A at 26GPa. The Nb and S atoms are at Wyckoff’s 2a positions(0.0, 0.0, 0.0) and 4e positions (0.0, 0.0, 0.34), respec-tively. More interestingly, the coordination number of Nbin I /mmm is eight. In this new type of bondings, oneNb atom and neighboring eight S atoms form a [NbS ]hexahedron. The coordination number of S is four. Toour knowledge, this type of covalent bondings has notbeen reported in TMD crystals. In general, in TMDs themetal atom has traditional six nearest neighbors. Thisnew type of eight nearest neighbors in I /mmm -NbS implies new potential chemical and physical properties,especially in two-dimensional crystals. (b) top view (c) front view(a) 2H'-NbS (d) left view(e) 2H-NbS (f) shift of 2H-NbS FIG. 2: (color online). The 2 H -NbS and 2 H (cid:48) -NbS crystalstructures. (a) The structure of 2 H (cid:48) -NbS . (b) Top view of2 H (cid:48) -NbS . (c) Front view of 2 H (cid:48) -NbS . (d) Left view of 2 H (cid:48) -NbS . (e) The structure of 2 H -NbS . (f) The shift directionof one layer of 2 H -NbS to form 2 H (cid:48) -NbS (top view). B. Phase transition and structural stability ofNbS In order to obtain the phase-transition sequence ofNbS under compression, we calculated the energies forits different phases at 0 K and pressures from 0 to 200GPa. The enthalpies vs pressure data of different struc-tures with respective to 2 H -NbS are plotted in Fig. 4,from which we note at 0 K the previously known hexag-onal 2 H -NbS is stable up to 26 GPa. Above 26 GPa,NbS transits to the tetragonal I /mmm structure whichremains stable up to a very high pressure, 200 GPa, theupper limit of our interest. Upon compression, NbS ex-hibits a volume reduction of 10.6% at 26 GPa (Fig. 5).This volume reduction directly results in the decrease ofinterlayer distance and the aggregation of S atoms aroundNb atoms. Although at ambient conditions, 1 T - and 3 R -NbS have relatively higher energies than 2 H -NbS , theyhave been synthesized successfully. The energies of 2 H (cid:48) -and P m are close to that of 3 R -NbS , so we be-lieve they are both potential to be synthesized in exper-iment. After all, the trigonal P m R -NbS in the whole pressurerange. The energies of P
22 and P
22 structures are (a) I4/mmm(b) top view (c) front view(d) left view
FIG. 3: (color online). Predicted I /mmm crystal structure.(a) The structure of I /mmm -NbS . (b) Top view. (c) Frontview. (d) Left view. both much higher than that of I /mmm structure. Sothe Gibbs free energy barrier is to high for NbS to over-come.
26 GPa 2 HP m P T RP I mmm H' E n t h a l py d i ff e r e n ce ( e v / f . u . ) −1.5−1−0.500.51 Pressure (GPa)0 50 100 150 200 FIG. 4: (color online). Enthalpy differences of predictedstructures relative to 2 H -NbS structure under high pressure. The mechanical stability of 2 H - and I /mmm -NbS -10.6% 2 HI mmm V o l u m e ( Å / f . u . ) FIG. 5: (color online). The equation of states of NbS . Thevertical dash curve indicates the volume reduction of NbS atthe phase transition point, 26 GPa. are confirmed by their elastic constants (shown in Ta-ble II and Fig.6) according to the elastic criteria of thehexagonal systems, C > | C | , ( C + 2 C ) C > C , C > , (1)and tetragonal systems, C > , C > , C > , C > ,C > C , C + C > C , C + C ) + C + 4 C > , (2)respectively. The increasing of the elastic constants of I /mmm -NbS with pressure reflect its enhanced stabil-ity as pressure increases (Fig.6). The new 2 H (cid:48) structureis also mechanically stable at ambient and high pressureaccording to the elastic criteria of the hexagonal crys-tals. While, the P m C . It is also worthy to note that the shearmodulus C of 2 H -NbS increases with pressure, butthe C values of 2 H (cid:48) -NbS remain small as pressure in-creases. This implies 2 H -NbS is more stable than 2 H (cid:48) -NbS . So, it is easier to synthesize 2 H -NbS in experi-ment other than 2 H (cid:48) -NbS .To further check the dynamical stability of the newstructures, 2 H (cid:48) - and I /mmm -NbS , we determinedtheir vibrational frequencies using density functional per-turbation theory (DFPT), as implemented in theQUANTUM-ESPRESSO package. For the exchange-correlation functional we used the Perdew Zunger local-
TABLE II: The elastic constants of different NbS structuresunder high pressure. The pressure ( P ) and elastic constantsare all in GPa.Structure P C C C C C C H H (cid:48) P m C C C C C C E l a s ti c c on s t a n t s ( G P a ) FIG. 6: The high-pressure elastic constants of I /mmm -NbS . density approximation (LDA) and ultrasoft pseudopo-tential. We applied the Vanderbilt ultrasoft pseudopo-tentials for Nb and S with the valence electrons con-figurations 4s and 3 s , respectively. Theultrasoft pseudopotentials were generated with a scalar-relativistic calculation.We careful tested on k and q grids, the kinetic energycutoff, and other technical parameters to ensure goodconvergence of phonon frequencies. The kinetic energycutoff, the energy cutoff for the electron density, andthe k grids were chosen to be 40 Ryd., 400 Ryd., and16 × ×
16 Monkhorst-Pack (MP) meshes in both to-tal energy and phonon dispersion calculations, respec-tively. We applied the Gaussian smearing method withthe smearing width of 0.05 Ryd. For the dynamical ma-trices of the I /mmm structure, we used a 2 × × q grid,giving 8 wave vectors q in the irreducible wedge of thefirst BZ. For the 2 H - and 2 H (cid:48) -NbS , the q grid meshes (b)(a) ZY PNΣ ZΓΣYXΓ F r e qu e n c y ( c m - ) PNΣ ZΓΣYXΓ E n e r gy ( e V ) −10−7.5−502.55 FIG. 7: (color online). The band structure and phonon dis-persion curve of I /mmm -NbS at 60 GPa. (a) the bandstructure, (b) the phonon dispersion curve (left) and phonondensity of states (right). are 2 × ×
4, also giving 8 wave vectors.Phonon dispersion curves (Figs. 7(b) and 8(a)) do notshow any imaginary frequencies, indicating dynamicalstability of I /mmm - and 2 H (cid:48) -NbS . So we believe I /mmm - and 2 H (cid:48) -NbS are both mechanically and dy-namical stable. The phonon dispersion curve of 2 H -NbS are also presented in Fig. 8(b), compared to the experi-mental data. The agreement of our calculated phononfrequencies and the 300 K experimental data is quitegood. By comparing the phonon dispersion curves of2 H - and 2 H (cid:48) -NbS , we note the phonon frequencies of2 H (cid:48) -NbS exhibit softening near to A point (along theΓ-A, A-L and H-A directions), implying its metastabilitycompared to 2 H structure. This is consistent with theconclusions from the elastic constants calculations. C. Electron-phonon coupling andsuperconductivity
We calculated the superconducting transition temper-ature T c of NbS using the Allen-Dynes form of theMcMillan equation, T c = ω ln . (cid:20) − . λ ) λ − µ ∗ (1 + 0 . λ ) (cid:21) , (3) (b)(a) AHLAΓKMΓ F r e qu e n c y ( c m - ) F r e qu e n c y ( c m - ) FIG. 8: (color online). The phonon dispersion curve of 2 H (cid:48) -(a) and 2 H -NbS (b) at 0 GPa (lines). The experimentaldata [Ref. 10] of 2 H -NbS at 2 K (open diamonds) and 300K (filled circles) are also plotted for comparison. λ T C ( K ) P (GPa)0 20 40 60 80 100 FIG. 9: (color online). The comparison of calculated super-conducting critical temperatures with experimental results37. The calculated electron-phonon coupling coefficients λ are also plotted. where λ (= 2 (cid:82) ∞ α F ( ω ) /ω d ω ) is the electron-phononcoupling constant, ω ln the logarithmic average frequency,and µ ∗ the Coulomb pseudopotential. The logarithmicaverage frequency is calculated by ω ln = exp { λ (cid:90) ∞ dωα F ( ω ) lnω/ω } . (4) The Eliashberg spectral function, α F ( ω ), which mea-sures the contribution of the phonons with frequency ω to the scattering of electrons, can be written as, α F ( ω ) = 12 πN ( (cid:15) F ) (cid:88) qν γ qν ω qν δ ( ω − ω qν ) , (5)where N ( (cid:15) F ) is the EDOS at the Fermi level. Thelinewidth of the phonon mode was calculated from, γ qν = 2 πω qν (cid:88) kjj (cid:48) | g qνk + qj (cid:48) ,kj | δ ( (cid:15) kj − (cid:15) F ) δ ( (cid:15) k + qj (cid:48) − (cid:15) F ) , (6)where g qνk + qj (cid:48) ,kj is the electron-phonon coupling matrixelement. The Coulomb pseudopotential µ ∗ was takenthe typical value 0.10 in all the superconducting criticaltemperatures ( T c ) calculations.The calculated T c of 2 H - and I /mmm -NbS are plot-ted in Fig. 9, compared with recent experimental data. The resulting T c s of 2 H -NbS are in very good agreementwith experiment and increase with pressure. It is inter-esting that the T c of I /mmm structure is higher thanthat of 2 H structure and decreases with pressure. Thisis resulted from the stronger electron-phonon couplingcoefficients λ in I /mmm -NbS (Fig. 9). The phononcalculations indicate that I /mmm is unstable below 10GPa. The highest T c of I /mmm (at 10 GPa) is 17.83 K.From the electronic energy band structure of I /mmm -NbS (Figs. 7(a)), we note it is metallic. It is previouslyknown that 2 H -NbS is also metallic, so pressure doesnot change the metallic properties of NbS , but enhancesthe electron-phonon coupling effects and thus increasesthe the superconducting critical temperature. IV. CONCLUSIONS
In conclusion, we predicted three new 2 H (cid:48) -, P m I /mmm -NbS structures using the MAC crystalstructure prediction technique. The new 2 H (cid:48) -NbS canbe formed by shifting the layer of atoms along typicaldirection parallel to the layer plane. Based on enthalpycalculations, we found 2 H -NbS transits to the tetrago-nal I /mmm structure at 26 GPa. The new bondings in I /mmm form a [NbS ] hexahedron, which has not beenreported in TMD crystals. More interestingly, the super-conducting temperature of I /mmm -NbS is higher thanthat of 2 H -NbS and decreases as pressure increases, re-sulted from the stronger electron-phonon coupling coef-ficients λ in I /mmm -NbS . In the stability region of I /mmm structure, the highest T c is 17.83 K. V. ACKNOWLEDGMENTS
The research was supported by the National NaturalScience Foundation of China (11104127, 11104227), theNSAF of China under grant No. U1230201/A06, theProject 2010A0101001 funded by CAEP, and the ScienceResearch Scheme of Henan Education Department under Grand No. 2011A140019. ∗ Electronic address: [email protected] L. Wei, C. Jun-fang, H. Qinyu, and W. Teng, Physica B , 2498 (2010). Y. Ding, Y. Wang, J. Ni, L. Shi, S. Shi, and W. Tang,Physica B , 2254 (2012). S. Leb´egue and O. Eriksson, Phys. Rev. B , 115409(2009). P. Raybaud, J. Hafner, G. Kresse, and H. Toulhoat, J.Phys.: Condens. Matter , 11107 (1997). Q. Wang, K. Kalantar-Zadeh, A. Kis, J. Coleman, andM. Strano, Nature Nanotech. , 699 (2012). M. Chhowalla, H. S. Shin, G. Eda, L.-J. Li, K. P. Loh,and H. Zhang, Nature Chem. , 263 (2012). D. Moncton, D. JAxe, and F. DiSalvo, Phys. Rev. B ,801 (1977). I. Guillam´on, H. Suderow, S. Vieira, L. Cario, P. Diener,and P. Rodi`ere, Phys. Rev. Lett. , 166407 (2008). J. Ka˘cmar˘cik, Z. Pribulov´a, C. Marcenat, T. Klein,P. Rodi`ere, L. Cario, and P. Samuely, Phys. Rev. B ,014518 (2010). M. Leroux, M. Le Tacon, M. Calandra, L. Cario, M.-A. M´easson, P. Diener, E. Borrissenko, A. Bosak, andP. Rodi`ere, Phys. Rev. B , 155125 (2012). J. Wilson, F. DiSalvo, and S. Mahajan, Adv. Phys. ,117 (1975). H. Y. and R. Aoki, J. Phys. Soc. Jpn. , 1327 (1986). V. G. Tissen, M. R. Osorio, J. P. Brison, N. M. Nemes,M. Garc´ıa-Hern´andez, L. Cario, P. Rodi`ere, S. Vieira, andH. Suderow, Phys. Rev. B , 134502 (2013). C. Geantet, J. Afonso, M. Breysse, N. Allali, andM. Danot, Catal. Today , 23 (1996). N. Kumagai and K. Tanno, Electrochim. Acta , 935(1991). W. M. R. Divigalpitiya, R. F. Frindt, and S. R. Morrison,J. Phys. D: Appl. Phys. , 966 (1990). W. Ge, K. Kawahara, M. Tsuji, and H. Ago, Nanoscale , 5773 (2013). C. J. Carmalt, T. D. Manning, I. P. Parkin, E. S. Peters,and A. L. Hector, J. Mater. Chem. , 290 (2004). F. Jellinek, G. Brauer, and H. M¨uller, Nature , 376(1960). W. B. Clark, J. Phys. C: Solid State Phys. , L693 (1976). S. Onari, T. Arai, R. Aoki, and S. Nakamura, Solid State Commun. , 577 (1979). Z. L. Liu, arXiv:1303.2802 . Y. Wang, J. Lv, L. Zhu, and Y. Ma, Phys. Rev. B ,094116 (2010). Y. Wang, J. Lv, L. Zhu, and Y. Ma, Comput. Phys. Com-mun. , 2063 (2012). J. Lv, Y. Wang, L. Zhu, and M. Y., J. Chem. Phys. ,084104 (2012). G. Kresse and J. Hafner, Phys. Rev. B , 14251 (1994). G. Kresse and J. Furthm¨uller, Comput. Mater. Sci. , 15(1996). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994). G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999). Z.-j. Wu, E.-j. Zhao, H.-p. Xiang, X.-f. Hao, X.-j. Liu, andJ. Meng, Phys. Rev. B , 054115 (2007). S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. , 1861 (1987). S. Baroni, S. D. Gironcoli, A. D. Corso, and P. Giannozzi,Rev. Mod. Phys. , 515 (2001). P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni,I. Dabo, A. D. Corso, S. D. Gironcoli, S. Fabris, G. Fratesi,R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj,M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri,R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto,C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen,A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.Condens. Matter , 395502 (2009). J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981). D. Vanderbilt, Phys. Rev. B , 7892 (1990). V. G. Tissen, M. R. Osorio, J. P. Brison, N. M. Nemes,M. Garc´ıa-Hern´andez, L. Cario, P. Rodi`ere, S. Vieira, andH. Suderow, Phys. Rev. B , 134502 (2013). H. J. Monkhorst and J. D. Pack, Phys. Rev. B , 5188(1976). B. Allen and R. C. Dynes, Phys. Rev. B , 905 (1975). W. L. McMillan, Phys. Rev. , 331 (1968). G. M. Eliashberg, Sov. Phys.-JETP16