Alkali metals, with the simplest single-valence-electron configuration and remarkably high chemical reactivity, have long been considered as a paradigm in fundamental physics, especially in the development of the electronic theory of solids.1 Under ambient conditions, all alkali metals form a body-centered cubic lattice and exhibit intrinsic metallicity.2–4 However, under high pressure, they undergo complex structural phase transitions, accompanied by exotic electron behaviors, such as s → p/d charge transfer5,6 and valence electrons located in lattice interstitials and acting as anionic electrons (i.e., electrides7), as well as interesting electronic properties, including metal–semiconductor/insulator transitions,7–9 superconductivity induced by interstitial electrons,10 and so on. Thus, alkali metals remain an important and active topic in the field of high-pressure physics.11–16

On the other hand, alkali metals have the ability to form unconventional compounds, especially at high pressures, and this has played an increasingly important role in fundamental research.17–22 For instance, at high pressure, Na can react with He atoms, which are chemically very inert, resulting in the first thermodynamically stable compound Na₂He with an electride character, which can be attributed to the eight-center two-electron bonds of the empty Na₈ cubes.19 Recently, pressure-induced Li-rich Li–Al compounds have also been found to be electrides, in which Al atoms exhibit diverse oxidation states owing to differences in Li → Al charge transfer. Even more interestingly, the Li₆Al electride, with two-dimensional sheets of Li composed of five- and eight-membered rings, exhibits an unexpected coexistence of superconductivity and superionicity.20 Consequently, alkali metals have unique advantages in triggering novel physical and chemical properties.23

Alkali-metal borides with B-rich stoichiometry, stabilized at high pressures, have also attracted great interest, because the weak ionic interaction between alkali metals and B units facilitates the stabilization of new B allotropes, and the electrons donated by alkali metals can compensate for the electron deficiency of B atoms.24–26 This type of study has not only promoted the discovery of materials with interesting properties (e.g., hardness27 and superconductivity24), but has also become a routine way to access new B allotropes. Our group proposed two hitherto unknown B allotropes (B₄ and B₁₇) with the highest superconducting transition temperature among bulk boron allotropes at that time, which were successfully obtained by removing the Na atoms from the predicted high-pressure phases of NaB₄ and Na₂B₁₇.24,28 In the K/Rb–B system, the metastable KB₄ and the stable RbB₈ compounds were theoretically identified and considered as ideal precursors to obtain a two-dimensional B sheet (a double layer consisting of six-membered rings) and the bulk o-B₁₆ phase in an experimental synthesis.25,26 On the other hand, Ma and co-workers29 proposed the idea of designing Li-rich borides (Li_xB, x = 1, 1.5, 2, 2.5, 3, …, 8) and identified the Li₆B compound, in which B atoms have a full-shell electronic configuration and there are no electrons located in the interstitial regions.

The superionic state, as an intermediate between the solid and liquid states in which some atoms still vibrate around their equilibrium positions while the others can diffuse in certain ways, has broad application prospects in solid-state batteries.30 Electrides, characterized by the presence of electrons residing at interstitial regions, have also shown unusual properties (e.g., ultra-low work function,31 high catalytic activity,32 and superconductivity33–35). In addition, interstitial electrons can also mediate superionic behavior.36 However, superionic electrides are still rather rare. Considering that superionicity usually appears in compounds with weak bonds or light elements, alkali-metal-rich borides have an inherent potential to exhibit coexisting superionic and electride behaviors.

Sodium, with unusual chemical attributes and relatively light mass, has demonstrated a strong ability to stabilize many allotropes and compounds (e.g., NaSi₆37 and Na₆Au38) under different conditions. Therefore, we have wondered if modulation of the Na/B ratio combined with pressure could achieve the objective mentioned above. In this work, we perform an extensive first-principles structural search of Na/B-rich Na–B compounds at 100, 200, and 300 GPa, and identify eight previously unknown Na–B compounds with diverse electronic properties, including semiconductivity, metallicity, and superconductivity. Unexpectedly, we find that the Na-richest Na₉B compound is an electride that also exhibits superionic properties at 200 GPa and 2000 K.

The intelligence-based particle-swarm optimization algorithm, implemented in the CALYPSO code,39,40 has become an effective way to predict stable or metastable elemental solids and binary and ternary compounds,18,41–44 some of which have already been synthesized,45–49 thereby accelerating the discovery of new materials. Here, we also employed this algorithm to explore the low-enthalpy Na_xB_y phases (x = 1, y = 1, 2, …, 20; x = 2, y = 1, 3, 5, …, 19; x = 3, y = 1, 2, 20; x = 4, 5, …, 12, y = 1) at 100, 200, and 300 GPa. The calculations of structural relaxations and electronic properties were performed with the Vienna Ab Initio Simulation Package (VASP).50 The projector augmented wave (PAW)51 method was adopted, with 2p⁶3s¹ and 2s²2p¹ as valence electrons for Na and B atoms, respectively. The Perdew–Burke–Ernzerhof (PBE)52 generalized gradient approximation was used as the exchange-correlation potential. The validity of the pseudopotentials has been confirmed by the perfect fit of the Birch–Murnaghan equation of states with the full-potential linearized augmented-plane-wave (LAPW) method (Fig. S0, supplementary material).53 A plane-wave cutoff energy of 850 eV and a Monkhorst–Pack54 scheme with a k-point grid of 2π × 0.03 Å⁻¹ were adopted to obtain a well-converged enthalpy. Phonon calculations were performed using the supercell finite displacement method in the PHONOPY code.55 Electron–phonon coupling (EPC) calculations were conducted within density functional perturbation theory (DFPT),56 as carried out in the QUANTUM ESPRESSO package.57Ab initio molecular dynamics (AIMD) NVT simulations using a 2 × 2 × 1 supercell (160 atoms) for Imma Na₉B at 200 GPa were performed for 10 000 steps with a step size of 1 ps at given temperatures of T = 500, 1000, 2000, and 3000 K with the Nosé–Hoover thermostat,58 and data of the last 2 ps were extracted.

Inspired by recently reported pressure-induced unconventional stoichiometric compounds (e.g., Li₆Al20 and CeH₁₈59), we have carried out an extensive structure search with 42 stoichiometries of Na_xB_y with up to four formula units per simulation cell to find new Na–B compounds at the selected pressures of 100, 200, and 300 GPa. Based on the identified structures, the energetic stability of the Na_xB_y system [Fig. 1(a)] can be verified by calculating their formation enthalpies relative to the known Na/B elemental solids according to the formulaΔH(NaxBy)=H(NaxBy)−xH(Na)−yH(B)x+y.

Generally, the phases that fall on the solid lines of the convex hull in Fig. 1(a) are considered thermodynamically stable, and they do not decompose into elemental solids or other binary compounds, while those on the dashed lines of the convex hull are thermodynamically metastable or unstable, depending on the values of the decomposition enthalpy and the dynamic stability. Besides reproducing the phases (namely, C2/m Na₃B₂, I4/mmm Na₂B₃, I4/mmm NaB₄, and Pm Na₂B₁₇) already reported by our group at 100 GPa,24 we have identified eight novel compounds as being stable, including Na-rich phases (Imma Na₉B, P-1 Na₅B, Cmcm Na₂B, and P2₁/m Na₃B₂) and B-rich ones (Imm2 Na₂B₃, Imma Na₂B₃, Cmmm NaB₂, and Pmm2 NaB₇), whose stable pressure ranges are shown in Fig. 1(b). It should be noted that Na₉B becomes the alkali-metal-richest stoichiometry among the known high-pressure alkali-metal borides. On the other hand, the absence of imaginary frequencies on the phonon dispersion curves in the whole Brillouin zone demonstrates their dynamical stability (Fig. S1, supplementary material).

To evaluate the synthetic feasibility of these newly identified compounds, we have computed their enthalpies of formation in their stable pressure regions relative to Na/B elemental solids and Na–B compounds (NaB₃, NaB₁₅, and Na₃B₂₀) stabilized at ambient pressure. As depicted in Fig. S2 (supplementary material), each predicted compound exhibits the potential for synthesis as indicated by a negative enthalpy of formation. Additionally, the optimal synthetic pathway is clearly identified with a pentagram symbol. Furthermore, it is worth noting that the precursor materials required for synthesizing each compound vary. For instance, Na and NaB₁₅ serve as the starting materials for Na₉B, while Na₃B₂₀ and NaB₁₅ are the precursors for NaB₇.

Despite differences in chemical composition and structural symmetry [Figs. 2(a)–2(g)], the crystal structures of stable Na–B compounds exhibit a common feature of a host–guest structure, in which host structures prefer the Na/B-rich stoichiometry [Fig. 2(h)]. This is distinguished by comparing the coordination numbers of the Na/B atoms in NaB_n/BNa_n polyhedra, in which the central atom with the lower coordination number is the host. In detail, the degree of aggregation of the building block associated with Na atoms increases with Na content, which changes from a single atom in Pmm2 NaB₇ to a zigzag chain in Cmmm NaB₂ and a wrinkled ladder chain in Imma Na₂B₃, to a Na₈ ring in P2₁/m Na₃B₂, to a quadrangular chain and a Na₆ ring in Cmcm Na₂B, and eventually to a Na₁₅/Na₁₆ polyhedron in P-1 Na₅B/Imma Na₉B. The same trend of variation can be found with increasing B content. On the other hand, with increasing Na content in these compounds, the amount of charge obtained by B atoms gradually increases and finally converges [i.e., the B atom has a closed-shell electronic configuration, Fig. 3(a)], as observed in Li–B compounds.29

Na-rich Na₉B and Na₅B are stabilized in orthorhombic [space group Imma, Fig. 2(a)] and triclinic [space group P-1, Fig. 2(b)] structures, respectively. For the Na-richest Na₉B, a three-dimensional Na lattice is formed through Na₁₆ polyhedra sharing vertices and faces, in which irregular voids between these polyhedra can accommodate electrons donated by Na atoms, as will be discussed below. In Na₅B, the Na atomic motif is connected by Na₁₅ polyhedra through shared vertices and edges, and the size and number of voids are also decreased compared with those of Na₉B owing to the smaller polyhedron. Furthermore, the B atoms are located at the center of the polyhedra of these two structures and are hypercoordinated (16- and 15-fold coordination for Na₉B and Na₅B, respectively). The average lengths of the Na–Na/B bonds in Na₉B at 200 GPa and Na₅B at 300 GPa are 2.21/2.31 and 2.20/2.20 Å, respectively, which are comparable to the 2.17 Å of the Na–Na bond in Na₆Au38 and the 2.18 Å of the Li–B bond in Li₆B at 100 GPa.29

Cmcm Na₂B, P2₁/m Na₃B₂, Imm2 Na₂B₃, Imma Na₂B₃, and Cmmm NaB₂ [Figs. 2(c)–2(f) and S3, supplementary material], consist of a metallic network connected by Na atoms and a covalent structure formed by B atoms (Figs. S3 and S4, supplementary material). Interestingly, Cmcm Na₂B [Fig. 2(c)], can be considered as a variant of the host–guest structure tI19-Na.2 The host-type metallic Na network with an Na–Na distance of 1.80–2.63 Å at 200 GPa is the same as that of the tI19-Na phase. For the guest structure, the B atoms form a zigzag chain with a bond length of 1.55 Å, which not only places the guest position of the Na linear chain in the tI19-Na phase, but also inserts into the unoccupied octagonal channels, which can be attributed to the smaller radius of the B atom.

Na₃B₂ adopts a monoclinic structure with the P2₁/m space group at pressures ranging from 255 to 300 GPa [Fig. 2(d)]. Specifically, Na atoms arrange themselves into a flat plane formed by edge-sharing Na₈ rings in the ac plane. This plane is then staggered along the b axis to create an interconnected network of Na atoms with helical channels. On the other hand, B atoms form a stripe-like arrangement along the b axis, composed of B₃ triangles that share edges and are positioned at the centers of the helical Na channels.

Considering that the Imm2 and Imma phases of Na₂B₃ belong to the orthorhombic crystalline system [Figs. 2(e) and S3, supplementary material], here we have analyzed their phase-transition process, which can be understood from the variations of Na and B units. The low-pressure Imm2 phase is composed of alternating Na₄ tetrahedra and Na₂ dimers via vertex-sharing. With pressure, the Na₂ dimers along the a axis rotate toward the b axis, accompanied by ruptures of the Na₄ tetrahedra, leading to the formation of a folded ladder shape. With regard to the B units, Imm2 Na₂B₃ consists of alternating B₄ tetrahedra and B₄ quadrilaterals (two B₃ triangles sharing edges) through shared vertices. The increased pressure results in a greater degree of folding of the B₄ quadrilaterals, and eventually they transform into B₄ tetrahedra, which can be viewed as a linear chain made up of vertex-sharing B₄ tetrahedra. Meanwhile, the decreasing distances between Na/B atoms under high pressure play an important role in making the structure more compact.

NaB₂ is predicted to adopt a Cmmm symmetry [Fig. 2(f)], where the B atoms form a three-dimensional framework featuring octagonal channels along the b axis. This framework can be understood as being composed of two building blocks: a B₆ octahedron and a B₆ hexagon. Moreover, a Na network is formed within the octagonal channels, with an Na–Na distance ranging from 1.98 to 2.17 Å at 200 GPa. The Na network is embedded within the structure, contributing to the overall stability and arrangement of the compound.

The B-richest NaB₇ stabilizes in an orthorhombic structure with Pmm2 symmetry [Fig. 2(g)]. There, the B atoms form a three-dimensional covalent framework composed of B₁₈ hexadecahedra and B₁₄ icositetrahedra through face-sharing, and each Na atom is located at the center of the B₁₈ hexadecahedron. The B–B bond length is in the range of 1.54–1.74 Å at 300 GPa, which is shorter than the 1.80 Å in γ-B,60 confirming its covalent bonding characteristics, as evidenced by the electron localization function (ELF) analysis [see Fig. 5(a)]. Furthermore, the interaction between Na and B atoms, with a bonding length of 1.96–2.48 Å, is ionic, which is in agreement with the Bader analysis [Fig. 3(a)].

Diverse elemental compositions and atomic arrangements in compounds often result in varied and distinct properties.17 In the case of our predicted compounds, their properties exhibit a dependence on B content. Specifically, as the B content increases, these compounds undergo a transformation from being semiconductors (Na₉B and Na₅B) to metals (Na₂B, Na₂B₃, and NaB₂) or superconductors (Na₃B₂ and NaB₇), as observed in the calculated electronic band structures, the projected density of states (PDOS), and the superconductivity analysis (Fig. S5 and Table S3, supplementary material). Additionally, the B-poor compound Na₉B exhibits unexpected electride and superionic features. Therefore, we discuss their properties in order of increasing B content.

For Na₉B, the electron localization function (ELF) reveals that some of the valence electrons are localized around the B atoms, and the others, known as interstitial anion electrons (IAEs), are located in the BNa₁₆ polyhedral cavities in an irregular shape, showing an electride feature [Fig. 3(b)]. The Na–B bond has an ionicity, which is mainly attributed to the interaction between Na 3s and B 2s/2p_x/2p_z [Fig. 3(c)]. On the other hand, considering the local IAEs in the polyhedra, we can conclude that the nonmetallicity of Na₉B is attributable not only to the superimposed Na–B interaction caused by hypercoordination, but also to the interaction between Na cations and IAEs located in the polyhedral cavities. With pressure, Na atoms donate more electrons to interstitial regions and B atoms, and the localization of these electrons is enhanced, resulting in an increased bandgap in Na₉B [Figs. 3(d) and S6 and Table S2, supplementary material].

Pressure-induced superionicity in materials has recently attracted research interest,35,61 especially because of the observed coexistence of electride and superionic states in Li₆Al.20 Here, we have explored the potential superionicity in the Na₉B electride via AIMD simulations. According to the MSD and the atomic trajectories, we can calculate the diffusion coefficients D and analyze the motion behaviors of Na/B atoms at 500–3000 K to determine if the superionic state appears. At 500–1250 K, D_Na ≈ D_B ≈ 0, indicating that Na₉B is a solid state where both Na and B atoms vibrate around their equilibrium positions (Figs. S7 and S8, supplementary material). On heating to 2000 K, the heavier Na atoms exhibit obvious diffusion (D_Na = 2.6 × 10⁻⁵ cm²/s > 0, D_B ≈ 0) [Fig. 3(e)], and the compound enters a superionic state, with a diffusion coefficient comparable to those of other Na-based superionic compounds.62,63 When the temperature rises to 2850 K [Fig. S7(d), supplementary material], it begins to transform into a liquid state, while complete melting requires a higher temperature.

Among the superionic compounds studied so far, diffusion of lighter atoms is more commonly observed.20,35,61 However, there are also some compounds that exhibit superionic behavior in heavier atoms.64,65 It is important to highlight that the diffusion behavior of the heavier Na atoms in this work can be explained by an unusual mechanism, in which Na atomic motions are influenced by IAEs. Owing to the irregularity of the interstitial electronic morphology, the Na1 atom, located at the position where two polyhedra are connected [black dashed oval in Fig. 4(a)], is closer to the IAEs, and so it is easily attracted and deviates from its equilibrium position to position A. Meanwhile, this behavior causes the appearance of a vacancy and induces other Na atoms (e.g., Na2) to move toward it, and the Na atom at position A will also move to another equilibrium position (e.g., Na3) owing to repulsion by Na cations. This phenomenon explains the initial motion of Na atoms within a representative polyhedron. Given that the structure consists of multiple such polyhedra, numerous Na atoms within these polyhedra begin to move owing to the attraction of IAEs, resulting in the formation of multiple vacancies. Subsequently, Na atoms diffuse along various pathways within and between the polyhedra, ultimately leading to diffusion throughout the lattice [Figs. 4(b) and S9, supplementary material]. It is noteworthy that the superionic behavior in Na₉B is also influenced by the weak Na–B/Na–Na interactions. This suggests that achieving superionic behavior in the hP4 Na electride with strong Na–Na interactions is challenging (Table S1, supplementary material). In addition, at any given moment, only a small portion of the Na atoms in the BNa₁₆ polyhedra diffuse into adjacent polyhedra, and so the B atoms are still bound.

In condensed matter physics, in addition to studying macroscopic atomic movements, it is crucial to gain a deeper understanding of temperature-dependent variations of electronic properties.66,67 Interestingly, the electride characteristics of Na₉B persist from the solid to the liquid region, accompanied by morphological changes of the IAEs as the temperature increases (up to 3000 K) (Fig. S10, supplementary material). This behavior is similar to that observed in the hP4 Na electride at high temperatures.66 In the solid region, Na₉B is a semiconductor, and its bandgap gradually decreases with increasing temperature (Figs. S8 and S11, supplementary material). However, in the superionic region, it becomes a metal, and the electronic states at the Fermi level E_F increase with temperature (Figs. S8 and S11, supplementary material). The metallic behavior is primarily attributed to the contribution of Na 3p/3s electrons, which is related to the charge transfer from 3s to 3p orbitals of Na atoms [Fig. 3(i)]. It is worth noting that the absence of an Na d orbital in Na₉B distinguishes it from IAEs with d-orbital attributes in the hP4 Na electride.66

Compared with Na₉B, Na₅B contains fewer valence electrons contributed by Na atoms and smaller polyhedral cavities, and so the valence electrons are located mainly around the boron atoms, as shown by the ELF [Fig. 3(f)]. On the other hand, the calculated crystal orbital Hamilton population (COHP) of the Na–B bond is comparable to that of Na₉B [Fig. 3(g)], which is associated with the hybridized state of Na 3s and B 2s/2p_x/2p_y. Therefore, its nonmetallicity is a consequence of the superimposed ionic interaction between Na and B atoms. Furthermore, the bandgap of Na₅B gradually increases (Table S2, supplementary material), and many more electrons, especially from the B 2p and Na 3p orbitals, occupy the upper part of the valence band with pressure [Fig. 3(h)]. The former is associated with a band flattening induced by the strong electronic localization (Figs. S12 and S13, supplementary material), and the latter is attributed to the increasing amount of Na → B charge transfer (Table S2, supplementary material). The pressure-dependent variation of atomic orbital energy levels also supports the analysis above [Fig. 3(i)].

With increasing B content, compounds such as Na₂B, Na₃B₂, Na₂B₃, NaB₂, and NaB₇ exhibit metallic properties, primarily due to the contribution of the B 2p electrons, as indicated by the PDOS (Fig. S5). Among these compounds, the Na-rich Na₃B₂ shows superconductivity, albeit with a very low T_c of 0.01 K at 300 GPa, which can be attributed to the weak electron–phonon coupling (Table S3, supplementary material).

As far as the B-rich NaB₇ is concerned, this has a unique B cage and a high DOS at the E_F, which is promising with regard to enhanced superconductivity, and thus we focus on analyzing its electronic and superconducting properties. In the projected electronic band structure [Fig. 5(b)], there is a flat band along the Γ–Z direction and around the X point near the E_F, resulting in a Van Hove singularity with a large DOS at E_F [Fig. S5(g), supplementary material], and two steep bands cross the E_F along the X–U direction. It is worth noting that there are distinct nesting behaviors observed in the Fermi surfaces along the Γ–Z direction and around the X point, as evidenced by the pronounced peaks in the Fermi surface nesting function ξ(Q) (Fig. S14, supplementary material). According to Bardeen–Cooper–Schrieffer (BCS) theory,68 this obvious “flat band–steep band” characteristic, the large DOS at the E_F, and the Fermi surface nesting behaviors imply that Pmm2 NaB₇ is a potential superconductor.

The calculated electron–phonon coupling (EPC) parameter λ of NaB₇ is 0.60 at 300 GPa, which is comparable to the 0.70 of MgB₂ at 0 GPa,69 and higher than the 0.36 of C2/m ScB₆ at 500 GPa.70 According to the Eliashberg spectral function and the projected phonon density of states (PHDOS) [Fig. 5(c)], the vibrations of Na and B atoms are obviously coupled in the whole phonon frequency range. Low-frequency vibrations below 11.3 THz contribute 73.9% of the total λ, while high-frequency vibrations contribute 26.1%. The strong EPC dominated by low frequencies can be attributed to the phonon softening along the Γ–Z direction and around the X point [Fig. 5(d)], which is associated with the flat-band feature near the E_F and the behavior of the Fermi surface nesting (Fig. S14, supplementary material). According to the Allen–Dynes-modified McMillan equation,71 the estimated superconducting transition temperature T_c is 10.2 K at 300 GPa with the typical Coulomb pseudopotential parameter of μ* = 0.1, which is comparable to the 9.0 K of α-BeB₆ at 1 atm,72 and higher than the 1.0 K of Amm2 Mg₃B₁₀ at 40 GPa73 and the 0.2 K of Immm SrB₈ at 150 GPa.74

For pressure-induced superconductivity, dynamically stabilized pressure of materials is crucial for further experimental verifications and applications. Interestingly, the dynamic stability of NaB₇ can be maintained as low as 55 GPa (Fig. S15, supplementary material), which prompts us to investigate its pressure-dependent superconductivity [Fig. 5(e) and Table S3, supplementary material]. As pressure decreases, T_c gradually increases, reaching a maximum value of 16.8 K at 55 GPa (Fig. S16, supplementary material). Since the λ and ω_log are the two key parameters determining the value of T_c, their variations with pressure are also analyzed. In detail, as pressure decreases, λ increases and ω_log decreases. Thus, λ plays a leading role in the evolution of the value of T_c. The higher value of λ at low pressure can be attributed to the increased DOS at the E_F, especially from the 2s and 2p electrons of B atoms [Fig. 5(f)].

In summary, we have performed first-principles swarm-intelligent structure search calculations to construct the high-pressure (100–300 GPa) phase diagram of the Na–B system with various chemical compositions. Eight novel Na–B compounds (two semiconductors and six metals) have been identified as being stable. All compounds can be seen as host–guest structures, and they exhibit novel Na/B atomic building blocks regulated by the Na/B content. In particular, Imma Na₉B is the alkali-metal-richest boride and exhibits coexisting superionic and electride characteristics at 200 GPa and 2000 K. With increasing temperature, it also experiences a transition from a semiconductor to a metal related to Na 3p/3s states. Additionally, B-rich NaB₇, with unique B cages, has an estimated T_c of 16.8 K at 55 GPa, which is mainly due to the coupling of B 2p electrons and low-frequency phonons. Our work represents an important step in furthering understanding the structures and properties of high-pressure alkali-metal borides.

Acknowledgment. This work was supported by the Natural Science Foundation of China under Grant No. 21573037, the Postdoctoral Science Foundation of China under Grant No. 2013M541283, the Natural Science Foundation of Hebei Province (Grant No. B2021203030), and the Science and Technology Project of Hebei Education Department (Grant Nos. JZX2023020 and QN2023246). A.B. acknowledges financial support from the Spanish Ministry of Science and Innovation (Grant No. PID2019-105488GB-I00) and the Department of Education, Universities and Research of the Basque Government and the University of the Basque Country (Grant No. IT1707-22). This work was carried out at the National Supercomputer Center in Tianjin, and the calculations were performed on TianHe-1(A).