Protective layer enhanced the stability and superconductivity of tailored antimonene bilayer

TL;DR

This study predicts that electron doping and Ca-intercalation can induce superconductivity in antimonene bilayers, and demonstrates that protective layers like graphene and h-BN improve stability and raise the superconducting transition temperature to 9.6 K.

Contribution

It introduces protective layers to enhance stability and superconductivity in antimonene bilayers, a novel approach for two-dimensional superconductors.

Findings

Antimonene bilayers can become superconducting via doping and intercalation.

Protective layers significantly improve stability in air.

Superconducting transition temperature reaches 9.6 K.

Abstract

For two-dimensional superconductors, the high stability in ambient conditions is critical for experiments and applications. Few-layer antimonene can be non-degradative over a couple of months, which is superior to the akin black phosphorus. Based on the anisotropic Migdal-Eliashberg theory and maximally-localised Wannier functions, this work predicts that electron-doping and Ca-intercalation can transform $\beta$-Sb bilayer from a semimetal to a superconductor. However, the stability of antimonene bilayer in air trends to be decreased due to the electron doping. To overcome this drawback, two kinds of protective layers (graphene and $h$-BN) are proposed to enhance the stability. Interestingly, the superconducting transition temperature will also be enhanced to $9.6$ K, making it a promising candidate as nanoscale superconductor.

Equations

$$\lambda =\sum _{\mathbf{q}\nu }{\lambda }_{\mathbf{q}\nu }=2{\int }_0^{\infty }\frac{{\alpha }^{2}F(\omega )}{\omega }d\omega .$$

$$T_{\mathrm{C}}=\frac{{\omega }_{ln}}{1.2}\exp [-\frac{1.04(1+\lambda )}{\lambda -{\mu }^{\ast }(1+0.62\lambda )}],$$

$$\Delta (\text{k},i{\omega }_n)=\frac{\varphi (\text{k},i{\omega }_n)}{Z(\text{k},i{\omega }_n)},$$

$$\frac{N_s(\omega )}{N({\epsilon }_F)}=Re[\frac{\omega }{\sqrt{{\omega }^2-{\Delta }^2(\omega )}}].$$

Full text

I Introduction

Recently, few-layer black phosphorus (BP) as an emerging member of two-dimensional (2D) materials has drawn lots of research attentions due to its unique properties, such as high carrier mobility (up to $1000$ cm2V*-1s-1* at room temperature) Qiao et al. (2014) and anisotropic in-plane properties Tran et al. (2014). However, the stability of few-layer BP remains a big challenge in experiments, on account of decomposition in air within several hours Castellanos-Gomez et al. (2014); Island et al. (2015). Antimony, as another member of nitrogen group, has several 2D phases (referred to as antimonene), such as $\alpha$-, $\beta$-, $\gamma$-, and $\delta$-Sb Zhang et al. (2016a). Among these allotropes, the $\beta$-Sb owns the lowest energy Zhang et al. (2016a), which has been successfully prepared in recent experiments Ares et al. (2016); Ji et al. (2016). Comparing with black phosphorus, few-layer antimonene is more stable in ambient conditions, which can be non-degraded over a couple of months Ares et al. (2016); Ji et al. (2016). Interestingly, $\beta$-Sb transforms from a semimetal to a semiconductor with decreasing thickness to few layers, and the monolayer shows an indirect band gap of $\sim 2$ eV Zhang et al. (2016a); Kamal and Ezawa (2015). It is also predicted that antimonene owns high carrier mobility Zhang et al. (2016a), wide range optical absorption Ares et al. (2016), and high refractive index Singh et al. (2016), which may lead to potential applications in optoelectronics. Thus, lots of experimental and theoretical works have been done on antimonene recently.

Superconductivity in 2D materials has also attracted enormous interests Ge et al. (2015). Recently, Penev et. al. predicted that 2D boron monolayer would exhibit intrinsic phonon-mediated superconductivity Penev et al. (2016). Our recent study also predicted Mo2C monolayer to be a quasi-2D superconductor, whose superconducting transition temperature ($T_{\rm C}$) could be adjusted by functional groups Zhang and Dong (2017). Moreover, superconductivity could also be induced and enhanced in other semiconducting or semimetal 2D materials via various methods, e.g. by improving carrier density and electron-phonon coupling (EPC) Zhang et al. (2016b); Zhang and Dong (2016); Chi et al. (2015). Even though, the high stability in ambient conditions of 2D superconductors is critical for experiments and applications, which has been rarely touched in these theoretical studies. In this sense, tailored antimonene may be a good candidate for 2D superconductor. Thus, it is highly desired to carefully investigate the stability and superconductivity of tailored antimonene.

In this work, using density functional theory (DFT), anisotropic Migdal-Eliashberg theory, and maximally-localised Wannier functions (MLWFs), the stability and superconductivity of electron-doped and Ca-intercalated $\beta$-Sb bilayer have been investigated. Our calculation finds that it will transform $\beta$-Sb bilayer from a semimetal to superconductor by increasing the carrier density. However, when the doping density is $>0.8$ electrons/cell, the strong EPC would lead to dynamic instability of the $\beta$-Sb bilayer structure. To pursuit stable and superconducting $\beta$-Sb bilayer, two kinds of protective layers, graphene and $h$-BN, have been considered, which can improve the superconducting $T_{\rm C}$ to $9.6$ K.

II Model & methods

The electronic structure, molecular dynamics (MD), and climbing-image nudged elastic band (CI-NEB) calculations have been performed using the Vienna ab initio simulation package (VASP) with projector-augmented wave (PAW) potentials Kresse and Furthmüller (1996); Kresse and Joubert (1999). The generalized gradient approximation of Perdew-Burke-Ernzerhof (GGA-PBE) formulation are used with a cutoff energy $500$ eV. The phonon dispersion and electron-phonon coupling calculations are carried out using the Quantum-ESPRESSO distribution, which are calculated within the framework of density function perturbation theory (DFPT) Giannozzi et al. (2009). The full relativistic norm-conserving potential are generated by $ld1.x$, which have been well tested for further calculations. The cutoff energy for the expanding wave function is $65$ Ry. Structure optimization and electronic structure are repeated by using Quantum-ESPRESSO and the results are consistent with those obtained using VASP. In order to deeply understand the electronic and superconducting properies, the recently developed Wannier interpolation technique is also used in our calculation Giustino et al. (2007); Noffsinger et al. (2010). More computational details are given in the Supplementary Materials Sup .

III Results & discussion

III.1 Pristine antimonene bilayer

Among all known allotropes of antimonene, $\beta$-Sb is the ground state of monolayers Zhang et al. (2016a). Different from planar graphene/$h$-BN, the monolayer $\beta$-Sb owns the buckled structure due to $sp^{3}$ hybridization. The calculated electronic structure and phonon dispersion of $\beta$-Sb monolayer are shown in Fig. S1 Sup , in agreement with previous studies Zhang et al. (2015); Aktürk et al. (2015); Wang et al. (2015).

The $\beta$-Sb bilayer adopts the A-B type configuration (CG) [Fig. 1(a)], i.e. the stacking sequence of atomic planes is $AB~{}CA$ Zhang et al. (2016a), while the A-A CG with a mirror symmetry between two layers [Fig. 1(b)] owns a little higher energy (about $12.5$ meV per atom) than the A-B CG. These two CGs own the same point group, i.e. $D_{3d}$ with $12$ symmetry operations, and both are dynamically stable according to their phonon spectrums [Fig. S2(a-b)]. The optimized lattice constants are $4.142$ Å and $4.075$ Å for A-B and A-A CG, respectively.

From their corresponding Raman intensity [Fig. S2(c)], the interlayer breathing mode ($A_{\rm 1g}^{1}$) and in-plane shear mode ($E_{\rm g}^{1}$) of $\beta$-Sb bilayer locate in low frequency range (less than $50$ cm*-1*), due to the weak interlayer interaction. Since the calculated interlayer distance ($d$) is shorter in the A-B CG ($2.39$ Å for A-B and $3.07$ Å for A-A), the stronger interlayer interaction leads to harder breathing mode and in-plane shear mode. The results imply that considerable proportion of chemical interaction between layers in the A-B CG, beyond the pure vdW interaction Aktürk et al. (2015). Such chemical interaction leads to significant changes in the electronic structure of $\beta$-Sb bilayer, as shown in Fig. 2. Specifically, the band structure of A-B CG [Fig. 2(a)] shows the semimetallic characters (band gap $\sim 0$ eV), while the A-A CG remains an indirect semiconductor with a gap of $\sim 0.4$ eV. The valence band maximum (VBM) and conduction band minimum (CBM) in the A-A and A-B CG’s are mainly contributed by Sb2’s $s$- and $p_{z}$-orbitals. Furthermore, using the MLWFs, the nearest-neighbor interlayer hopping of $p_{z}$ orbitals reaches $0.41$ eV, a considerable magnitude. To further verify the crucial role of interlayer interactions in the electronic structure, the interlayer distance is artificially lifted to $3.0$ Å and $3.5$ Å, which gives rise to finite band gaps (Fig. S3).

III.2 Superconductivity of tailored antimonene bilayer

Doping is a powerful method to improve carrier density in materials. To verify the maximum doping concentrations in $\beta$-Sb bilayer, pure electron doping will be tested first. Five doping concentrations ($0.2$, $0.4$, $0.6$, $0.8$, and $1.0$ electrons/cell) are considered here, and the corresponding maximal carrier density is about $6.6\times 10^{14}$ cm*-1*. After doping, the in-plane lattice constants and atomic positions are further optimized, and the A-B CG still owns the lowest energy. As expected, electron doping makes $\beta$-Sb bilayer transform from a semimetal to a metal [Fig. S4 and Fig. 3(a-b)]. It should be noted that the spin-orbit coupling (SOC) has very limited influence on the electron structure of electron-doped A-B CG [Fig. S4 and Fig. 3(a-b)], although Sb is a heavy element.

The calculated phonon dispersions and EPC constant $\lambda_{\bf{q}\nu}$ for electron-doped antimonene bilayer are shown in Fig. S5. In addition, the detailed definition and calculation descriptions of EPC are given in the Supplementary Materials Sup . When the doping concentration is up to $1.0$ electrons/cell, there is a quadratic phonon of the $ZA$ mode (out-of-plane acoustical mode) becoming imaginary round the $\Gamma$ point due to the strong EPC [Fig. 3(c-d)]. According to the density of states (DOS) [Fig. 3(a-b)], the Fermi level at $1.0$ electrons/cell doping is moved to near the Van Hove singular point, resulting in the structural instability. Therefore, the maximum doping concentration is around $7.0\times 10^{14}$ cm*-1* for free-standing $\beta$-Sb bilayer. The biaxial tension can remove this negative frequency around the long-wavelength limit, where $ZA$ mode becomes linear dispersion ($\omega\sim|{\bf q}|$) in the vicinity of $\Gamma$ point Penev et al. (2016); Zhou et al. (2015).

To study the superconductivity of electron-doped $\beta$-Sb bilayer, the EPC constant $\lambda$ is also calculated by summing over the first Brillouin zone (BZ), or integrating the $\alpha^{2}F(\omega)$ in the $\bf{q}$ space Allen and Dynes (1975); Allen (1972):

[TABLE]

The obtained $\alpha^{2}F(\omega)$, phonon density of states (PDOS, $F(\omega)$) and $\lambda(\omega)$ are plotted in Fig. 4(a) and Fig. S6. Comparing with $F(\omega)$’s, all $\alpha^{2}F(\omega)$’s have similar shape, implying that all vibration modes contribute to the EPC and corresponding frequency region is very remarkable (red area of Fig. 4(a) and Fig. S6). Due to the factor $1/\omega$ in the definition of $\lambda$ (see Eq. 1), the contributions from high-frequency is negligible, which are apparently shown in Fig. 4(a) and Fig. S6. Exactly, the frequency regions of $A_{\rm 1g}^{1}$ and $E_{\rm g}^{2}$ correspond to area-A and area-B, respectively. The obvious peaks of $\alpha^{2}F(\omega)$ exist in these two areas, indicating that these two vibrational modes have great effects on the superconductivity of electron-doped antimonene bilayer. The $T_{\rm C}$ can be estimated using the Allen-Dynes modified McMillan equation Allen and Dynes (1975):

[TABLE]

where $\mu^{*}$ is the Coulomb repulsion parameter and $\omega_{ln}$ is the logarithmically averaged frequency. When taking a typical value $\mu^{*}=0.1$, the estimated $T_{\rm C}$ and corresponding superconductive parameter are shown in Fig. 4(b). When the doping concentration is less than $0.6$ electrons/cell, all $T_{\rm C}$’s are under $1$ K, and the effect of SOC is not obvious. The EPC is rapidly enhanced when more electrons are added to $\beta$-Sb bilayer and the SOC may also promote EPC. The highest $T_{\rm C}$ can reach $2.6$ K.

As aforementioned, the strong EPC in electron-doped $\beta$-Sb bilayer restricts its doping concentration ($<7.0\times 10^{14}$ cm*-1*). In the following, the superconductivity of $\beta$-Sb bilayer under the real Ca-doping conditions will be investigated. Considering the symmetry, there are four mostly possible configurations for Ca intercalation: AB-I CG [Fig. 5(a)], AB-II CG [Fig. S8(a)], AB-III CG [[Fig. 5(b)])] and AA CG [Fig. S8(b)]. After full optimization of all configurations (details are given in Supporting Materials Sup ), the AB-I CG owns the lowest energy, while the AB-III CG is a little higher in energy ($\sim 3.4$ meV/cell) and others are even higher.

The calculated band structures without/with SOC for AB-I CG (AB-III CG) are shown in Fig. 6(c) and (e) (Fig. S9). As expected, the $\beta$-Sb bilayer transforms from a semimetal to a metal by Ca intercalation. Differing from above pure electron-doped case, there are two bands (blue and orange curves in Fig. 6(a)) crossing the Fermi level, and their corresponding Fermi surfaces are shown in Fig. 6(b). According to the projection, the hole pocket Fermi surfaces around the $K$ points are mainly contributed by the Ca layer, while the sun-shape Fermi surface around the $\Gamma$ point is constructed by the combination of $\beta$-Sb bilayer and Ca layer. When SOC is considered in calculations for the AB-I CG, the spin degeneracy of the valence and conduction bands along the line $\Gamma-K$ is broken, as shown in Fig. 6(e-f). Generally, the Kramer theorem leads to the degeneracy at the $\Gamma$ and $M$ points. For other $k$ points, the spin degeneracy is determined by the time-reversal and $D_{3h}$ point-group symmetry. Due to lack of inversion-reversal symmetry in $D_{3h}$, the spin splitting ($\Delta(k,\varphi)=\beta(k)\left|\sin 3\varphi\right|$) is induced at the $K$ point Henk et al. (2003), as found in the MoS2 monolayer Zhu et al. (2011). However, for the AB-III CG (Fig. S9) with both the inversion-reversal and time-reversal symmetry, SOC can not induce spin splitting in the whole BZ, which is similar to aforementioned electron-doped antimonene bilayer.

Like aforementioned electron-doped $\beta$-Sb bilayer, all Raman modes are soften after Ca intercalation (Table S2) in comparison to pristine $\beta$-Sb bilayer. The charge transfer from Ca layer to Sb layer [Fig. 7(a)] results in ionic interaction between Ca and Sb. In addition, the breathing mode $A_{\rm 1g}^{1}$ of $\beta$-Sb bilayer is not sensitive to the Ca intercalation, differing from other 2D materials (e.g. MoS2 Zhang et al. (2016b)) with full interlayer vdW interaction. Interestingly, through electron localization function (ELF) analysis [Fig. 7(b)], the generated Coulomb repulsion between charged Sb-Sb pairs weakens its covalent bond, causing remarkable red shift of resting Raman modes.

The phonon dispersion and corresponding EPC of AB-I and AB-III are also calculated, as shown in Fig. 8 and S10. Due to the Fermi contour nesting, the strong EPC is happened around $K$ point, where $E_{\rm g}^{1}$ and $E_{\rm g}^{2}$ make great contributions to the EPC. Comparing AB-I with AB-III, the $\lambda_{\bf{q}\nu}$ is obviously different for frequency below $70$ cm*-1* due to their different stacking and group points. The vibrational modes at low frequency in AB-III CG have stronger EPC than those in AB-I CG.

Our results for $\alpha^{2}F(\omega)$, PDOS, and $\lambda(\omega)$ of AB-I and AB-III CG’s are shown in Fig. 7(c) and Fig. S11, respectively. Similar to the pure electron-doping case, $\alpha^{2}F(\omega)$ and $F$($\omega$) of these two have similar peaks, indicating all vibration modes contribute to EPC. Specially, there are three areas (blue shading in Fig. 7(c)) where obvious peaks of $\alpha^{2}F(\omega)$ exist in the AB-I CG. These frequency ranges of area-A, area-B, and area-C correspond to the ZA, $E_{\rm g}^{1}$, and $E_{\rm g}^{2}$ modes, respectively. Therefore, these three modes promote the superconductivity in Ca-intercalated $\beta$-Sb bilayer. The superconducting $T_{\rm C}$ is also estimated based on Eq. 2, which are $4.9$ and $5.3$ K ($6.8$ and $6.9$ K) for the AB-I CG (AB-III CG) without and with SOC, respectively. Thus, the SOC can slightly improve the superconductivity of Ca-intercalated $\beta$-Sb bilayer, particularly for the case without inversion symmetry, which in agree with the previous study Sklyadneva et al. (2013).

From above electronic and phonon calculations, one can notice that the Fermi surface and phonon dispersion are anisotropic due to the remarkable difference between $\Gamma-M$ and $\Gamma-K$, which may raise the question whether the EPC is also anisotropic. Hence, the proper treatment of anisotropic EPC based on Wannier function is adopted in following study. Moreover, the superconducting gap is also calculated, which is given by Margine and Giustino (2013):

[TABLE]

where $Z(\textbf{k},i\omega_{n})$ and $\phi(\textbf{k},i\omega_{n})$ are the mass renormalization function and the order parameter, respectively. From this quantity, the leading edge of the superconducting gap for electron momentum k on the Fermi surface is obtained by solving for $\omega$ in $Re\Delta(\textbf{k},i\omega_{n})=\omega$. The $T_{\rm C}$ of superconductor is obtained from $\Delta_{0}\neq 0$ at the highest temperature, while $\Delta_{0}=0$ in the whole temperature range for usual metal. To calculate $\Delta$, the DFT band structure is fitted by using Wannier functions. The obtained superconducting gap of AB-I is shown in Fig. 7(d) as a function of temperature. When taking a typical value $\mu^{*}=0.1$, there is nearly only one superconducting gap, which is closed around $7.2$ K. These results indicate the superconducting gap function of Ca-intercalated $\beta$-Sb bilayer is very weakly anisotropic and the corresponding $T_{\rm C}$ is about $7.2$ K, only slightly higher than that estimated using McMillan’s approach. To accurately verify the $T_{\rm C}$, real part of the normalized quasiparticle density of states is calculated below and above $T_{\rm C}$, which is based on isotropic approximation:

[TABLE]

From Fig. 7(e), one can see that Van Hove singularity exists in the normalized quasiparticle density of states at $6.5$ K, which disappears at $7.5$ K. Hence, $T_{\rm C}$ should be between $6.5$ K and $7.5$ K. The ratio between the gap and $T_{\rm C}$ gives: $2\Delta_{0}(k_{\rm B}T_{\rm C})^{-1}=3.42-3.94$ ($k_{\rm B}$: Boltzmann constant), very close to the ideal BCS value of $3.53$, which is also an evidence for phonon-medium BCS superconductivity in Ca-intercalated $\beta$-Sb bilayer.

In the end of this subsection, the way how to insert Ca into the interlayer of $\beta$-Sb bilayer is proposed. First, it is difficult for Ca to directly cross through the atomic gaps of $\beta$-Sb due to the large energy barrier (about $2.5$ eV according to the NEB result). Even though, there remain two possible ways to insert Ca, i.e. using boundary and vacancy. By inserting Ca from boundary, our MD simulation indicates that Ca atoms could spontaneously go into the interlayer of $\beta$-Sb bilayer. However, due to the stronger interlayer interaction than those in graphene and $h$-BN, there is a large energy barrier (about $1.2$ eV) preventing Ca’s freely moving in interior at room temperature. By inserting Ca from vacancy of $\beta$-Sb, there are mainly three processes, including generating vacancy, inserting Ca via vacancy, and restoring vacancy. According to our MD simulation, Ca atoms could cross through the vacancy with the help of thermal fluctuation, but they could not spontaneously move in the interlayer, similar to above boundary-inserting case. Hence, additional driving forces (e.g. electric current) are required for internal moving within the interlayer of $\beta$-Sb bilayer, which deserve further studies.

III.3 Protective layers

The superconductivity of $\beta$-Sb bilayer has been predicted upon Ca intercalation. The stability in ambient conditions is the first condition for experiments, which needs careful theoretical investigations. To simulate the Ca-intercalated $\beta$-Sb bilayer in air, here the interaction between O2 and AB-I CG is studied using a $4\times 4$ supercell. There are three most possibly adsorption sites for O2 on its surface, as shown in Fig. S12, where the C-site owns the lowest energy. The corresponding binding energy ($E_{b}=E_{\rm total}-E_{\rm AB-I}-E_{\rm O_{2}}$) is calculated to be $-0.9$ eV, indicating strong absorption. The absorption process can be intuitively simulated as the atomic relaxation process in DFT calculation, as shown in Fig. 9, suggesting no energy barrier from the physical absorption to chemical absorption.

In order to improve the stability and protect the superconductivity of Ca-intercalated $\beta$-Sb bilayer, two kinds of protective layers (graphene and $h$-BN) are considered to avoid the effect of oxygen or water (Fig. 10(a)). Considering the physical similarity between AB-I and AB-III CG, here only the AB-I CG is selected as a represent, while the conclusions can be naturally extend to cover AB-III CG. In addition, to reduce the lattice mismatch between interlayer and protective layers, $(\sqrt{3}\times\sqrt{3})R30^{\circ}$ supercell of graphene and $h$-BN are used, as shown in Fig. 10(b), where lattice mismatch is less than $2\%$. After full optimization of atomic positions, the distances between the $\beta$-Sb bilayer and protective layers are around $3$ Å, suggesting the vdW-type coupling. The corresponding electronic structure for graphene- and $h$-BN covered $\beta$-Sb layers are shown in Fig. 11(a) and  11(b), respectively. Due to the insulating nature of $h$-BN, it has negligible effect on the Fermi surface of Ca-intercalated $\beta$-Sb bilayer. However, the graphene’s Dirac cone will affect the Fermi surface around the $\Gamma$ point.

The dynamic stability of these two systems are also confirmed by phonon dispersion [Fig. S13 (a-b)], and the occupied phonon states of protective layers are mainly in high frequency region due to their small masses and strong $sp^{2}$ covalent bonds. To verify the effect of protective layers on the superconductivity, the corresponding $\alpha^{2}F(\omega)$, PDOS, and $\lambda(\omega)$ are also calculated, as shown in Fig. 11(c-d). As expected, the superconductivity is mainly concentrated at the interlayer, while the impact of covered layers are tiny according to the shape of $\alpha^{2}F(\omega)$ and PDOS. The obtained $T_{\rm C}$ based on McMillan’s approach are $7.3$ and $9.6$ K for the graphene- and $h$-BN covered layers, respectively. The improvement of $T_{\rm C}$ may be part from lattice mismatch (about $1.6\%$ tensile stress). To verify this guess, the strain effect on pure Ca-intercalated $\beta$-Sb bilayer is also calculated, which gives $T_{\rm C}$ up to $6.7$ K. In this sense, the covered layers may also partially participate in EPC, besides the pure strain effect. In particular, the superconductive properties of $h$-BN covered bilayer is superior than the graphene covered one, due to less disturbance on Fermi surface, making it a promising candidate as a nanoscale superconductor.

IV Conclusion

In summary, using the first-principles calculations, it is found that different stacking leads to varying interlayer interaction and interlayer hopping parameter, which have great influence on electronic structure of $\beta$-Sb bilayer. In addition, electron-doping and Ca-intercalation can transform the $\beta$-Sb bilayer from a semimetal to a superconductor. However, the strong EPC would lead to dynamic instability of $\beta$-Sb bilayer, when the doping density is $>0.8$ electrons/cell. Moreover, the stability of antimonene bilayer in air is decreased after the electron transfer. For further experiments and applications, two kinds of protective layers (graphene and $h$-BN) have been designed, which have little influence on the Fermi surface of Ca-intercalated $\beta$-Sb bilayer. More importantly, the superconducting $T_{\rm C}$ is enhanced to $9.6$ K with the help of these protective layers.

Acknowledgements.

Work was supported by the National Natural Science Foundation of China (Grant Nos. 11674055, 11834002), the Scientific Research Foundation of Graduate School of Southeast University (Grant No. YBJJ1769), and Fundamental Research Funds for the Central Universities. Most calculations were done on Tianhe II of National Supercomputer Center in Guangzhou (NSCC-GZ)