Stability of B2 compounds: Role of the $M$ point phonons

Although many binary compounds have the B2 (CsCl-type) structure in the thermodynamic phase diagram, an origin of the structural stability is not understood well. Here, we focus on 416 compounds in the B2 structure extracted from the Materials Project, and study the dynamical stability of those compounds from first principles. We demonstrate that the B2 phase stability lies in whether the lowest frequency phonon at the $M$ point in the Brillouin zone is endowed with a positive frequency. We show that the interatomic interactions up to the fourth nearest neighbor atoms are necessary for stabilizing such phonon modes, which should determine the minimum cutoff radius for constructing the interatomic potentials of binary compounds with guaranteed accuracy.

Introduction. Since the characterization of the crystal structure of CsCl in 1921 [1], many of the B2 (called CsCl-type) compounds have been synthesized experimentally, resulting in that the B2 structure is the most common phase in the thermodynamic phase diagram of binary compounds. Analyses of materials database have shown how frequent they appear [2][3][4]. More recently, Kolli et al. identified 267 parent crystal structures that can generate their derivative ordered phases, and showed that the body-centered cubic (bcc) structure is the most common parent crystal structure [4] by using the Materials Project (MP) database [5]. Among them, the B2 structure is the most common ordering on bcc structure. Apart from the ambient condition, the B2 structure may also appear: the B2 structure can be transformed from the B1 (NaCl-type) structure under high pressure (such as alkali halides and alkali-earth oxides [6,7]) and from the L1 0 (CuAu-type) structure in warm dense matter regime [8].
Although the elastic stability of the B2 structure has been studied in a wide variety of compounds [9][10][11], such a stability does not always yield the dynamical stability against the zone boundary phonon excitations [12]. For the parent bcc structure, it is well known that the transverse acoustic phonon at the N point in the Brillouin zone (BZ), propagating along the diagonal direction of any two axes in the cubic cell, has relatively low frequencies [13,14]. It has been shown that such a phonon is stabilized by long-range interatomic interactions, allowing alkali metals to form the bcc structure at the ambient condition [15]. By considering the fact that the B2 structure is equivalent to the bcc structure when two species are assumed to be the same element, we expect that a similar scenario holds: the lowest frequency phonon at the M point that is stabilized by the long-range interactions determines the stability of the B2 compounds. We have recently confirmed that the long-range interatomic interactions up to the fifth nearest neighbor (5NN) atoms are needed to understand the dynamical stability of the CuAu in the L1 0 structure [8].
The range of the interatomic interactions is of prime importance in the field of atomistic modeling of con-densed matters. For the bcc elemental metals, the cutoff radii should be larger than the 3NN or 4NN distances [16][17][18]. However, it has not been understood why more than 3NN distances are required to describe the potential energy surface accurately.
In this paper, we study the dynamical stability of 416 compounds in the B2 structure from first principles, by assuming zero temperature and pressure. We show that 266 out of 416 compounds are dynamically stable, and demonstrate that the stability of the B2 compounds is mainly determined by whether the lowest frequency phonon at the M point in the BZ is endowed with a positive frequency. In addition, we develop a force constant model taking into account the interatomic interactions up to the 6NN atoms, and demonstrate that the interatomic interactions up to the 4NN atoms are enough to stabilize the lowest energy phonons at the M point. This should determine a minimum cutoff radius of the interatomic potentials for binary compounds. The present work unveils the microscopic mechanism of the B2 compounds stability, which stimulates other studies involving different crystal structures.
Methods. By using the MP [5] and the pymatgen [19], we first extracted the list of the B2 compounds. By setting the space group to P m3m with the number of atoms in a unit cell being two and excluding the atoms in the actinide series, 416 B2 compounds having the inorganic crystal structure database (ICSD) IDs were found. We next optimized the lattice parameter a for the 416 compounds. All the density-functional theory (DFT) calculations were performed with the Quantum ESPRESSO (QE) [20] using the Perdew-Burke-Ernzerhof (PBE) [21] functionals of the generalized gradient approximation for the exchange-correlation energy and using the ultrasoft pseudopotentials provided in the pslibrary.1.0.0 [22]. Spin-polarized approximation was used for all calculations. The cutoff energies for the wavefunction and the charge density are 60 Ry and 600 Ry, respectively, and 20×20×20 k grid was used in the self-consistent field (SCF) calculations [23]. The SCF convergence threshold was set to be 10 −8 Ry and the smearing parameter of σ = 0.02 Ry [24] was used for all calculations. The total energy and forces were converged within 10 −5 Ry and 10 −4 a.u., respectively. The accuracy of the present calculations was checked by comparing the optimized a's with those in the MP [5]. We have confirmed that the optimized a's agree with the reference values within an error of 1 % except for the 11 Ce-based compounds, ClO, and NCl. The optimized a's in the Ce-based compounds are larger than the reference values by a few percent. This may be due to the absence of the f -electrons in the present calculations, resulting in no magnetic moments, whereas the Ce-based compounds show ferromagnetic phase in the MP [5]. The error of the a's for ClO and NCl were 4.6 and 11.1 %, respectively.
Although the formation energies E form of the 416 compounds can be obtained from the MP [5], the dynamical stability properties are not always obtained. We thus performed phonon dispersion calculations based on density-functional perturbation theory (DFPT) [25] implemented in QE. The threshold parameter for the selfconsistency (tr2 ph) was set to be 10 −14 , and 4×4×4 q grid (10 q points) was used. We calculated the phonon dispersions along the symmetry lines Γ-X-M -Γ-R-X. When the phonon frequency ω is imaginary, the phonon energy is represented as a negative value, − |ω|, with the Planck constant . In the present work that adapts a finite size q grid, we identify the B2 compound as dynamically stable if the lowest phonon energy is larger than ε min = −1 meV.
Phonon dispersions. We have found that 266 out of 416 compounds are identified to be dynamically stable. of the transverse acoustic branch at the M point tends to be small, compared to that at the other symmetry points (except for Γ) in the BZ. In a similar manner, the instabilities of the ZnCu and ZrCu are due to the phonon softening around the M point, leading to negative phonon energies. The CuN shows strong phonon softening over the entire BZ. This may be ascribed to the positive formation energy (1.22 eV/atom [5]), implying that the B2 CuN is decomposed into the Cu crystal in the face-centered cubic structure and the N 2 molecules.
For the B2 compounds including the platinum-group metals (Os, Ru, Ir, Rh, Pt, and Pd), Hart et al. have identified that 16 compounds have the B2 structure as the ground state by using the DFT calculations within the PBE [26]. We have also identified that the 16 compounds are dynamically stable. In addition, we find that the B2 FeRh is dynamically stable, which is consistent with the experimental synthesis [27], whereas the FeRh in the B2 structure has been predicted to be unstable [26].
Overall, the dynamical stability of the B2 compounds is determined by the M point phonon energy. The phonon dispersions for the 416 compounds are shown in the Supplemental Material [28].
The M point phonons. In order to understand to what extent the M point phonon determines the dynamical stability of the B2 compounds, we plot the minimum phonon energy within the entire BZ ( ω min ) as a function of the lowest phonon energy at the M point ( ω M ) in Fig. 2. When ω M ≤ 0, the relationship of ω M = ω min holds, whereas when ω M > 0, the equality of ω min = 0 holds. These show that the stability of the M point phonon is a good descriptor for identifying the dynamical stability of the B2 compounds.
In Fig. 2, we can find some exceptions that do not follow the B2 stability criteria above. There are 16 compounds (AlPt, CsN, KS, LuMg, MgHg, MgTl, MnAu, MnHg, NaS, PrMg, SiRh, TcB, TiRh, VRu, YMg, and ZrRh) that satisfy both ω M > 0 and ω min < ε min . Although the M point phonon does not have the minimum phonon energy in these compounds, negative phonon energies appear along the Γ-M line in the BZ (see the Supplemental Material [28]). The instability of VRu is anomalous because negative phonon energies are observed only along the Γ-R line in the BZ, while the experimental synthesis of the VRu in the B2 and tetragonally distorted B2 structures [29,30] has been reported.
The value of E form may be another descriptor for understanding the stability of the B2 compounds: the dynamically stable compounds have negative E form and, in turn, the unstable compounds have positive E form . However, many exceptions have been found: A similar issue has been found in a wide variety of materials, such as ordered alloys [31,32] and many two-dimensional materials [33][34][35]. The metastability of materials has to be studied in detail.
It is noteworthy that the B2 compounds including a semiconducting element (group 14-17) tend to be unstable (see Table S1 and Figure S1 [28]). Such compounds include the strongly bonded systems mentioned above. The instability of these may be due to the presence of different ground state structure such as the B1 structure. The understanding for the group dependence of the dynamical stability is left for future work. Effect of long-range interatomic interactions. We next study an origin of the positive value for the phonon energies at the M point and identify the role of the longrange interatomic interactions. Based on standard lattice dynamics [36], we derived analytical expressions for the phonon frequencies for the wavevector q = (π/a, π/a, 0), at the M point in the BZ, by assuming the B2 compounds that consist of atoms κ = A and B with the masses M κ . By using the translational symmetry of the crystal, the dynamical matrixD is given bỹ where D κκ ′ αβ (R j , 0) is the force constant matrix, i.e., the force along the direction of α acting on the atom κ in a unit cell characterized by the lattice vector R j when the atom κ ′ in the cell of R = 0 is displaced along the direction of β. With the q grid used in the present work, the R j can take the vectors of 3 i=1 m i a i with m i = −1, 0, 1, 2, where a i 's are the primitive lattice vectors in the cubic cell. We thus considered the force constants up to the 6NN sites of B2 AB. By assuming that the atom A is located at the origin, the position of the first, second, third, fourth, fifth, and sixth NNs are B(a/2, a/2, a/2) , A(a, 0, 0), A(a, a, 0), B(3a/2, a/2, a/2), A(a, a, a), and A(2a, 0, 0), respectively, and these equivalent sites, where the numbers of the equivalent sites are 8, 6, 12, 24, 8, and 6, respectively. Due to the equivalence between the x and y directions, four different phonon modes are present: the z-polarized modes having ω 2 zκ =D κκ zz (q) (κ = A or B) and the x-and y-polarized modes having with The value of ω 2 ± s in Eq. (2) does not change when the indexes x and y are replaced with y and x, respectively, leading to the doubly degenerated modes. The displacement patterns of the normal modes with ω ± are described by a combination of the x (y)-polarized vibration of the atom A and the y (x)-polarized vibration of the atom B, which are shown in Fig. 3. From the expressions of Eqs. (2) and (3) Fig. 4. For clarity, the cases of the 15 Cu-based compounds are shown. The phonon frequencies are denoted as ω k with k = 1, · · · , 6 in an ascending order. The lowest energy phonons (ω 1 and ω 2 ) correspond to the normal modes having the frequency ω − except for CuN. For the low energy phonon modes, the deviation from the DFPT is large when p ≤ 3. The p = 4 is a critical value for the convergence of the phonon energy, and such a p corresponds to the interatomic interactions between different atoms A and B, which is consistent with our expectation above. On the other hand, p = 2 is found to suffice for the high energy phonons. These results indicate that the low and high energy phonons at the M point are stabilized by the long-range and short-range interatomic interactions, respectively. The instability of the M point phonon modes in the CuN was due to the negative value of ω 2 zκ =D κκ zz (q) with κ = Cu. The interatomic distance with the p = 4 atoms is √ 11a/2 ≃ 1.66a, and the total number of the NN atoms up to p = 4 is 50. This might be a minimum criteria for determining the cutoff radius (R c ) of the interatomic potentials in the AB-type compounds with guaranteed accuracy. Interestingly, Seko et al. have constructed an interatomic potential for the bcc K with a = 5.284Å, and shown that the energy, force, and stress calculated by using the interatomic potential agree well with those obtained by the DFT calculations when the value of R c is set to be more than 9Å [17]. Such a R c is quite similar to √ 11a/2 = 8.762Å.
Some compounds have a large discrepancy of the lowest phonon energy between the 6NN model and the DFPT. This implies that more long-range interactions are required to achieve the convergence to the DFPT results. It should be noted that the stability of the d electron compounds might not be described accurately within the PBE approximation. The value of E form tends to be overestimated within the PBE when one studies the B2 compound that consists of the atoms having the completely filled d orbitals such as Cu, Ag, Au, Zn, and Cd [37]. More analysis using other functionals is beyond the scope of this work. The comparisons between the pNN and the DFPT for the 416 compounds are provided in the Supplemental Material [28].
Conclusions. In conclusion, we studied an origin of the dynamical stability for 416 B2 compounds extracted from the MP [5], and identified that the 266 out of 416 compounds are dynamically stable. The lowest phonon energy at the M point determines whether the B2 compound is dynamically stable, and the interatomic interactions up tp the 4NN atoms are necessary for an accurate description of the B2 phase stability, which should determine the minimum cutoff radius for constructing the interatomic potentials for the AB-type compounds. The instability of the other compounds may be attributed to (i) the positive formation energy of the B2 structure (such as CuN) and/or (ii) the presence of different structures (such as CsF).
The present work assumes the B2 ordering on bcc structure. A competition between the B2 ordering and the finite temperature effects play an important role to understand the stability of binary [38] and high-entropy alloys [39]. Exploring a phonon stabilization under such situations is left for future work. In addition, the role of the spin-orbit coupling on the metastability [40] needs to be studied when heavy elements are included. We hope that this work stimulates other studies involving different crystal structures and gives a useful insight for developing the interatomic potentials of binary systems. This work was supported by JSPS KAKENHI (Grant No. JP21K04628). A part of numerical calculations has been done using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo.