Metastability relationship between two- and three-dimensional crystal structures: a case study of the Cu-based compounds

Some of the three-dimensional (3D) crystal structures are constructed by stacking two-dimensional (2D) layers. To study whether this geometric concept, i.e., using 2D layers as building blocks for 3D structures, can be applied to computational materials design, we theoretically investigate the dynamical stability of copper-based compounds CuX (a metallic element X) in the B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_h$$\end{document}h and L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_1$$\end{document}1 structures constructed from the buckled honeycomb (BHC) structure and in the B2 and L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_0$$\end{document}0 structures constructed from the buckled square (BSQ) structure. We demonstrate that (i) if CuX in the BHC structure is dynamically stable, those in the B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_h$$\end{document}h and L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_1$$\end{document}1 structures are also stable. Using molecular dynamics simulations, we particularly show that CuAu in the B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_h$$\end{document}h and L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_1$$\end{document}1 structures withstand temperatures as high as 1000 K. Although the interrelationship of the metastability between the BSQ and the 3D structures (B2 and L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_0$$\end{document}0) is not clear, we find that (ii) if CuX in the B2 (L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_0$$\end{document}0) structure is dynamically stable, that in the L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_0$$\end{document}0 (B2) is unstable. This is rationalized by the tetragonal Bain path calculations.

the (111) axis the B1 is transformed into the B3 structure. In this way, the existence of CuX in the B h structure has not yet been explored. Note also that among binary metallic phases, L1 1 as well as B h structures are quite rare structure, as has been pointed out in Ref. 25 and observed in AlSn solid solutions 26 . Therefore, it will be interesting to study whether the design of CuX is possible based on the geometric concept (i.e., from 2D to 3D) combined with the dynamical stability analysis beyond the formation energy analysis 25 .
In this paper, we investigate the interrelationship of the dynamical stability between the 2D and the 3D structures by focusing on CuX with a metallic element X. By calculating the phonon density-of-states (DOS) from first principles, we demonstrate that if CuX in the BHC structure is dynamically stable, CuX in the B h and L1 1 structures are also stable. This indicates that the BHC structure serves as a building block for the B h and L1 1 structures, while in some compounds the dynamical stability of the B h and L1 1 structures relies on the coupling between the BHC structures along the c axis. By performing the molecular dynamics simulations, we particularly show that CuAu in the B h and L1 1 structures are stable up to 1000 K.
We also study the dynamical stability of CuX in the the buckled square (BSQ), B2, and L1 0 structures, where the latter two structures can be constructed by stacking the BSQ structures (see Fig. 1). However, we find it difficult to discuss the metastability relationship between the 2D and 3D tetragonal structures because most of CuX in the BSQ structure are unstable. Analogous to the instability of the bcc (fcc) elemental metals in the fcc (bcc) structure 27 , we can find that if CuX in the B2 (L1 0 ) structure is dynamically stable, that in the L1 0 (B2) structure is unstable. This is visualized by means of the tetragonal Bain path calculations.
The metastability relationships between different crystal structures are summarized in Table 1: The first two relationships between the fcc, hcp and bcc structures can be applicable to the stability of elemental metals 27 ; the third is the relationship between the 2D elemental metals; the fourth and fifth are the relationships between the 2D and 3D elemental metals; and the others are the relationships between the 2D and 3D compounds. It has been shown that the (m, n) Lennard-Jones crystals also satisfy the first and the fifth relationships when m + n ≥ 18 , i.e., the short-ranged potentials 28 . We hope that the present study paves the way to computational materials design based on the interrelationship between the 2D and 3D structures. Making the metastability database for other compounds XY and exploring other 2D structures as building blocks for the 3D structures are left for future work.
The present study will serve as one of the geometry-based approaches that have been proposed in computational materials science. Hart et al. have generated superstructures as a derivative of a parent crystal structure 29 . In contrast, Kolli et al. have mapped complex crystal structures onto several parent structures 30 . They have pointed out that 15 of the top 20 parent structures can be generated by vertically stacking simple 2D structures including the close-packed, honeycomb, kagome, square, and more complex structures. As an extension of the present investigation, these will be candidate structures for creating the metastability database of the ordered compounds XY. Figure 1. Illustrations of the 2D (BHC and BSQ) and 3D structures (B h , L1 1 , B2, and L1 0 ) and its interrelationships: Stacking the 2D structures vertically produces the 3D structures and the tetragonal and trigonal distortions of the B2 structure yield the L1 0 and L1 1 structures, respectively. BSQ, B2, and L1 0 . We next study the dynamical stability of CuX in the BSQ, B2, and L1 0 structures. The phonon DOS of these structures are also provided in the Supplementary Information. As shown in Fig. 3, only CuBe, CuRh, and CuCu in the BSQ structure are dynamically stable, while most of CuX in the B2 or L1 0 structures are dynamically stable. Therefore, it seems to be difficult to find a useful relationship of the dynamical stability between the BSQ and the 3D structures. Interestingly, one can find that if CuX in the B2 structure is dynamically stable, that in the L1 0 structure is unstable, and vice versa: if CuX in the L1 0 structure is dynamically stable, that in the B2 structure is unstable. This is analogous to the metastability relationship between the bcc and the fcc structures in the elemental metals listed in Table 1. This will be rationalized by the tetragonal Bain path calculations below. www.nature.com/scientificreports/ In the present study, CuX in the B2 structure for X = group 1 (K and Rb), group 2 (Be, Mg, Ca, Sr, and Ba), group 3 (Sc, Y, and Lu), and Pd have been identified to be dynamically stable. In contrast, CuBe 16 , CuSc 18 , CuY 19 , CuZr 21 , CuPd 17 , and CuZn 20 in the B2 structure have been synthesized experimentally at ambient condition. For the case of the L1 0 structure, 20 CuX in the L1 0 structure have been predicted to be dynamically stable, whereas only CuAu has been synthesized experimentally 23 . It has been shown that the lowest energy phonon at the R point, (0, π/a, π/c) , of the L1 0 CuAu is stabilized by the long-range interatomic interactions from the first to more than the fifth nearest neighbors 31 . If the same scenario holds, the stabilization of the phonons at the N point, (0, π/a, π/a) , and the R point will be a key to synthesize the ordered compounds in the B2 and L1 0 structures, respectively. Formation energy. As a measure of relative stability compared to elemental metals Cu and X, Fig. 4 plots the formation energies E j (CuX) defined as Eq. (1) for j = BHC , B h , and L1 1 (triangles) and for j = BSQ , B2, and L1 0 (circles). The CuX in the BHC and BSQ structures have positive E j and those in the 3D structures are energetically more stable. The B h and/or L1 1 structures have negative E j for X = Li , Be, Sc, Y, Lu, Pd, Pt, Zn, Al, and Ga, while the B2 and/or L1 0 structures have negative E j for X = Li , Be, Mg, Ca, Sc, Y, Lu, Ti, Zr, Pd, Pt, Au, Zn, Al, and Ga. The group-dependence of the energetic stability is observed, that is, the CuXs with the group 2, 3, 12, and 13 elements of X tend to have low energy, while those with the group 6 elements of X have high energy. This tendency is similar to that predicted for the Pb-and Sn-based compounds 32 .  www.nature.com/scientificreports/ Tetragonal and trigonal paths. To understand the metastability relationship between the B2 and the L1 0 structures proposed above, we show the total energy curves along the tetragonal Bain path in Fig. 5, where the volume of CuX is fixed to that of the B2 structure. If the total energy takes the minimum value at c/a = 1 (B2), no energy minimum is observed at c/a > 1 (L1 0 ), as for X = K , Rb, group 2, and group 3. In contrast, if the total energy takes the minimum value at c/a > 1 , then that takes the maximum or higher values at c/a = 1 . It should be noted that the anomalous stability of CuPd in the B2 structure (see the group 10 in Fig. 3) can be understood as the shift of the energy minimum to c/a = 1.
The B2 structure can be transformed into the L1 1 structure through the trigonal path, elongating or shortening the cubic structure along the (111) direction with the volume fixed to that of the B2 structure. Figure 6 shows the total energy as a function of cos γ , where γ is the angle between the primitive lattice vectors (see Methods for the details). For most of X, one can observe a double well-like energy curve with a potential energy barrier height of about 1 eV, indicating that the B2 ( cos γ = 0 ) and L1 1 ( cos γ ≃ 0.8 ) structures are energetically stable. The energy maximum at cos γ = 0.5 corresponds to the B1 structure. For X = Ca , Sr, Ba, Sc, Y, and Lu, the  www.nature.com/scientificreports/ energy curve around cos γ = 0.5 is almost flat. In these compounds, the L1 1 structure was transformed into the B1 structure after the geometry optimization (see the Supplementary Information). The stability properties of CuX along the tetragonal and trigonal paths in Figs. 5 and 6 are almost consistent with those in Figs. 2 and 3. However, the energetic stability along the tetragonal and trigonal paths may not guarantee the dynamical stability of the structure 27,33 . This is because the phonon calculations investigate whether the Hessian of the potential energy is positive definite around the crystal structure considered, whereas the transformation path calculations restrict atomic displacements to only one distortion direction assumed. For X = Ti , Zr, and Hf (group 4), although the trigonal path calculations predict that the B2 and L1 1 structures are stable, CuX in those structures are unstable (see Figs. 2 and 3). This implies that other transformation paths may exist, i.e., the compound will show a phase transformation due to a negative curvature in the potential energy surface around the B2 and L1 1 structures. It will be desirable to perform more systematic calculations that investigate the interrelationship between metastable structures, as demonstrated in elemental metals 34 . CuAu in the B h and L1 1 structures. The experimental synthesis of the BHC CuAu has been reported recently 6 . Based on the phonon calculations, we confirmed that such a structure is dynamically stable 15 and predicted that the B h and L1 1 structures are also dynamically stable, as they are constructed from the 2D structures. We discuss the temperature effect on the stability in the 3D structures below.
For the B h structure, the lattice parameters are a = 2.858 Å and c/a = 1.526 . For the L1 1 structure, the lattice parameter is a = 4.726 Å and the primitive vectors are (u, v, v), (v, u, v), and (v, v, u) in units of a/ √ 3 with u = 0.239 and v = 0.669 . The total energy of the B h and L1 1 structures are energetically higher than that of the L1 0 structure, the ground state structure of CuAu, by 0.176 eV and 0.138 eV per a cell. These values are much higher than the thermal energy at the ambient condition. As shown in Fig. 6, the potential barrier height ( ∼ 1 eV) between the L1 1 and the B2 ( ≃ L1 0 ) structures is also high enough not to yield the phase transition. Therefore, once synthesized, these metastable structures will be stable. Figure 7(a) and (b) shows the temperature (T)-dependence of the atomic distribution in the B h supercell at the end of the molecular dynamics (MD) simulation (5 ps). The B h structure is stable up to T = 1000 K because the structure of each layer (i.e., the 2D hexagonal layer of Au and Cu) is still preserved. When T = 1500 K, the hexagonal symmetry in each layer is broken. When T is increased to 2000 K, the atomic displacements are www.nature.com/scientificreports/ significant, so that the layered structure is not observed. Similar tendency is observed in the L1 1 structure after the MD simulation of 4 ps, as shown in Fig. 7(c) and (d).

Methods
The total energy and the optimized structures of CuX were calculated within the density-functional theory implemented in Quantum ESPRESSO (QE) 35 . We used the exchange-correlation functional of Perdew, Burke, and Ernzerhof (PBE) 36 and the ultrasoft pseudopotentials in pslibrary.1.0.0 37 . We used the cutoff energies of 80 and 800 Ry for the wavefunction and the charge density, respectively, used 20 × 20 × 1 k grid and 20 × 20 × 20 k grid for 2D and 3D structures, respectively 38 , and set the smearing parameter of σ = 0.02 Ry 39 . For the BHC and BSQ structures, the size of the unit cell along the c axis was set to be 14 Å to avoid the spurious interaction between different unit cells. For the geometry optimization of the BHC, B h , BSQ, B2, and L1 0 structures, the initial guess of the lattice parameters was set to the same as the optimized value obtained by our previous calculations 15 . The optimization calculations of the L1 1 structure was started by assuming cos γ = 0.82 ( γ = 35 • between the primitive lattice vectors) and the unit cell volume equal to that of the B h structure. For the tetragonal Bain path calculations, we first calculated the unit cell volume of the B2 structure, 0 = a 3 0 with the optimized lattice parameter a 0 , and next deformed the shape of the unit cell but by keeping the volume 0 = a 2 c fixed, where the primitive lattice vectors are defined as (a, 0, 0), (0, a, 0), and (0, 0, c). For the trigonal Bain path calculations, we also fixed the unit cell volume to be 0 . By setting the input parameter ibrav to −5 in QE 35 , we defined the primitive lattice vectors as (u, v, v), (v, u, v), and (v, v, u) in units of a/ √ 3 . Given a parameter c = cos γ = (v 2 + 2uv)/(u 2 + 2v 2 ) that is equal to the cosine of the angle between any pairs of primitive lattice vectors, in turn, we determined the values of a, u, and v by using the . We performed the spin-polarized calculations, where the total energy and forces are converged within 10 −5 Ry and 10 −4 a.u., respectively, which are smaller than the default values by a factor of ten. In our previous study, we performed the spin-unpolarized calculations by using the default values of the convergence criteria 15 . In the spin-polarized calculations, the finite magnetic moment was observed for CuX with X = Cr , Mn, Fe, Co, and Ni when a ferromagnetic state was assumed as an initial guess. As listed in Table 2, the size of the magnetic moment in the BHC and BSQ structures seems to be larger than that in their 3D counterparts. We also found that a few compounds show a finite magnetic moment: For the BHC structure, CuTi (0.87); for the B2 structure, CuHg (0.47), CuIr (0.47), CuRe (2.28), CuRh (1.02), and CuZr (0.01); and for the L1 0 structure, CuHf (0.01), CuRh (0.72), and CuZr (0.05), where the figure in parenthesis is in units of µ B per a cell. It has been shown that the 2D elemental metals of Ba, group 2 (Sc and Y), group 3 (Ti, Zr, and Hf), V, group 8 (Ru and Os), group 9 (Rh and Ir), and group 10 (Pd and Pt) can have a finite magnetic moment 12 . This has been attributed to a decrease in the coordination number as well as a change in the electron DOS shape. Therefore, it will be helpful to study the effect of the electronic band structure on the dynamical stability of CuX, in addition to the geometric concept investigated in the present work. It will also be interesting to study the competition of the energetic stability between the ferromagnetic and antiferromagnetic phases in a systematic manner. However, these are beyond the scope of the present work.
The formation energy of CuX in the structure j per a unit cell including two atoms was calculated by where ε j (CuX) is the total energy of CuX in the structure j. ε min (X) is the minimum total energy of X, where this is determined by calculating the energy of the bcc, fcc, and hcp structures. Note that the functional dependence of the E j of ordered compounds has been studied in detail [41][42][43][44] . The use of the PBE functional can predict the value of E j accurately for strongly bonded systems, while it may underestimate E j for weakly bonded systems such as CuAu. The calculated values of E j as well as those in the Materials Project (MP) database 40 were provided in the Supplementary Information. To show the numerical accuracy within the PBE calculations, we compared the lattice parameters with those in the MP database 40 , and the experimental values for the B2, L1 0 , and L1 1 structures, as listed in Table 3. Except for CuPt, the agreement between the present calculations, the MP database 40 , and the experiment is good (the deviation is less ± 0.03 Å). The lattice parameters of CuX in the structure j were also provided in the Supplementary Information. www.nature.com/scientificreports/ The dynamical stability and the instability of CuX in the structure j were determined by the phonon DOS calculations within the density-functional perturbation theory 45 implemented in the QE code 35 . We used more than a 6 × 6 × 1 q grid for the 2D structures (seven q points for the BHC and ten q points for the BSQ structures), 3 × 3 × 3 q grid (six q points) for the B h and L1 0 structures, 3 × 3 × 3 q grid (four q points) for the B2 structure, and 2 × 2 × 2 q grid (four q points) for the L1 1 structure.
For the case of CuAu in the B h and L1 1 structures, the MD simulation was performed by using the QE 35 . The 3 × 3 × 3 supercell (54 atoms) was used by considering only the Ŵ point in the Brillouin zone. The ionic temperature T is controlled by the velocity scaling method, where T = 500 , 1000, 1500, and 2000 K were assumed. The time step of the MD simulation was set to 20 a.u. (0.96 fs) and the number of the steps was set to 5250 (5.08 ps) and 4200 (4.06 ps) for the B h and L1 1 structures, respectively.
Received: 21 May 2021; Accepted: 6 July 2021 Table 3. The lattice parameters of CuX in the B2 structure, CuAu in the L1 0 structure, and CuPt in the L1 1 structure. For the L1 1 structure, the parameters of the conventional unit cell (including six atoms) are listed. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.