Chemical bonding in dense A^IVX^VI and A^V₂X^VI₃ chalcogenides: electron-deficient multicenter bonds in electron-rich elements

Abstract

Chemical bonding in dense chalcogenides of the A^IVX^VI and A^V₂X^VI₃ families has been a subject of significant debate, with competing models trying to explain their unique and exceptional properties, in particular those corresponding to phase change materials. We present a paradigm-shifting advancement by demonstrating that electron-deficient multicenter bonds (EDMBs) occur in the dense crystalline chalcogenides of these two families with electron-rich elements, trying to resolve the long-standing controversy between metavalent and hypervalent bonding models—a fundamental question that has made difficult to get a rational materials design for decades. Through a comprehensive theoretical investigation of monochalcogenides (GeS, SnSe, PbTe) and sesquichalcogenides (As₂S₃, Sb₂Se₃, Bi₂Se₃, Bi₂Te₃) under compression, we discuss and explain why heavier chalcogenides naturally exhibit phase change, thermoelectric, and topological properties at room pressure while lighter compounds require compression to exhibit these exceptional properties—a phenomenon that has puzzled researchers for years. Our systematic analysis uncovers a universal pattern in EDMB formation following two or three well-defined stages, depending on the specific elements involved, that agrees with the recently published unified theory of multicenter bonding. Our unprecedented application of three local energy densities (kinetic, potential, and total) at bond critical points provides quantitative validation of the EDMB model and provides a complementary and unified explanation for the anomalous sixfold coordination in these compounds that defies the conventional 8 − N rule. Most significantly, our work establishes clear structure–property relationships that provide a direct pathway to rational design of new materials with tailored properties.

---

Introduction

Heavy chalcogenides of A^IVX^VI and A^V₂X^VI₃ families (A = group 14 and 15 elements; X = chalcogens) have attracted significant attention due to their exceptional properties and diverse applications. These materials serve as the foundation for various technological applications, including high-performance thermoelectric devices based on tetrel-chalcogen compounds (such as PbTe,¹ PbSe,² GeTe,³ and Bi₂Te₃⁴) and topological insulators utilizing the trigonal tetradymite-like R3̄m phase found in several sesquichalcogenides (including Sb₂Te₃, Bi₂Se₃, and Bi₂Te₃^5,6).

The remarkable potential of these crystalline materials stems from their unique property portfolio, which includes moderate electrical conductivity and unusually high values of Born effective charges (Z*), optical dielectric constants (ε_∞), and Grüneisen parameters for transverse optical modes (γ_TO). Additionally, these materials exhibit Effective Coordination Numbers (ECoN) that deviate from predictions based on the 8 − N rule, also known as Pearson's rule.^7,8 These distinctive properties have long been attributed to unconventional chemical bonding mechanisms within these materials.^9,10

The nature of chemical bonding in binary chalcogenides of the A^IVX^VI and A^V₂X^VI₃ families has become a focal point of recent research. At room pressure (RP), the crystalline structures of most compounds in these families (excluding phase change materials (PCMs)) typically show geometries that can be described as Jones-Peierls distorted versions of the sixfold coordinated cubic rocksalt (rs) phase.^11–13 This distorted structural arrangement is consistent with the simultaneous presence of primary covalent bonds and secondary non-covalent bonds. These secondary bonds are typically associated with stereochemically-active lone electron pairs (LEPs) that result in low-symmetry crystalline structures.¹¹ PCMs, however, present more symmetric crystalline structures due to the lack of stereochemical activity of the LEPs, however, the nature of the chemical bonds is more complex. Their unique and unconventional properties have recently led to the development of two competing bonding models, which aim to explain the exceptional properties of PCMs: the metavalent bonding model and the hypervalent bonding model.

The metavalent bonding model has proposed that the crystalline phases of heavy A^IVX^VI and A^V₂X^VI₃ chalcogenides exhibit a novel “brand” of chemical bonding: the metavalent bonding.^8,14 The authors proponent of the metavalent bonding model previously considered that the PCMs featured “resonant covalent bonds”.¹⁵ Metavalent bonds are proposed to be electron-deficient two-center-one-electron (2c-1e) bonds, in contrast to conventional single covalent bonds that are two-center-two-electron (2c-2e) bonds.^16,17 Proponents of this model have also suggested that the unique properties of PCMs are directly linked to the effective masses of charge carriers, which in turn relate to the mechanism of metavalent bond formation. On the other hand, the hypervalent bonding model suggests that the crystalline phases of heavy A^IVX^VI and A^V₂X^VI₃ chalcogenides are characterized by the century-old electron-rich multicenter bonding (ERMB). Electron-rich multicenter bonds (ERMBs), also referred to as hypervalent bonds or hyperbonds, encompass a variety of bonding schemes characterized by electron delocalization involving more than two centers. While the 3-center-4-electron (3c-4e) configuration serves as a canonical example, proponents of the ERMB model emphasize a broader spectrum of delocalized interactions, potentially involving higher center counts.^18–20 In this context, the properties and geometries of PCMs are attributed to this electron-rich multicenter delocalization.

These two competing models are supported by different methodological approaches: the metavalent bonding model is primarily supported by density-based methods (analyzing the topological properties of electron density), while the hypervalent bonding model is backed by orbital-based methods (examining orbital wavefunctions). Proponents of the two models have used specialized computational tools for these analyses, such as CRITIC2²¹ for topological analysis of electron density and the local orbital basis suite towards electronic-structure reconstruction (LOBSTER)^22–24 for studying orbital wavefunctions. It is worth noting that independent researchers have also identified multicenter bonds in tetradymite-like A^V₂X^VI₃ chalcogenides, such as Bi₂Te₃,²⁵ and that both methodologies are based on theoretical calculations performed according to the density functional theory (DFT).

While both density-based and orbital-based methods offer valuable insights, disagreements between these perspectives have persisted^26–28 and recent analyses suggest that the metavalent and hypervalent bonding models for PCMs may be converging toward a common understanding.²⁹ This convergence has led to the proposal of a third bonding model, the electron-deficient multicenter bonding (EDMB) model, that aims to resolve the controversy between the previous two approaches.³⁰ Electron-deficient multicenter bonds (EDMBs) are also known as three-center-two-electron (3c-2e) bonds, which have been found to form also 2c-1e bonds in electron-rich systems.^9,10,30 The EDMB model proposes that the unique properties of the crystalline phases of heavy chalcogenides stem from the presence of EDMBs in their crystalline structures.^9,10,30 This unified approach acknowledges both the electron-deficient character of this bonding (as proposed by the metavalent model) and the multicenter nature of this bonding (as suggested by the hypervalent model). By integrating these perspectives, the EDMB model provides a more comprehensive framework for understanding chemical bonding in these materials. In this context, and to clarify the EDMB nature of bonding in PCMs, the proponents of the EDMB model have emphasized the importance of distinguishing between the two types of multicenter bondings: ERMBs and EDMBs. These considerations have culminated in the recent formulation of a new unified theory of multicenter bonding that provides hints to solve the issue.^9,10

High pressure (HP) serves as a powerful tool for investigating changes in chemical bonding, as it enables the gradual modulation of interatomic distances. The application of HP typically induces crystalline phase transitions that transform low-coordinated structures observed at RP into high-coordinated structures at HP. Therefore, HP leads to an increase in electronic density in materials. As emphasized in the unified theory of multicenter bonding, the increase in the electronic density of a material can be also induced by chemical reduction (contrary to chemical oxidation) and by the chemical substitution of elements by their heavier counterparts, so HP plays a similar role to chemical reduction and chemical substitution.^9,10 However, HP offers a distinct advantage over reduction and substitution: it allows for a gradual and controlled increase in electronic density due to the subtle modification of interatomic distances. Therefore, HP studies enable us to observe the fine details of chemical bond transformation as a continuous process, with each pressure point representing a photogram of a movie.

In this study, we explore in depth the chemical bonding characteristics of the A^IVX^VI and A^V₂X^VI₃ chacogenides, providing detailed explanations to clarify each aspect thoroughly. We demonstrate that the dense phases of these chalcogenides exhibit EDMBs, as predicted by the unified theory of multicenter bonding.^9,10,30 While the low-coordinated and low-symmetry crystalline structures of A^IVX^VI and A^V₂X^VI₃ compounds with light elements are mainly characterized by distinct primary iono-covalent bonds plus secondary non-covalent interactions, the high-coordinated and high-symmetry crystalline structures in these families (typical of PCMs) can feature EDMBs or mixed bonding that includes EDMBs. To reach this conclusion, we build upon the methodological approaches previously employed in our theoretical framework,⁹ which integrated concepts from both the metavalent (density-based) and hypervalent (orbital-based) bonding models to examine their respective advantages and limitations for describing PCMs. The present study extends this approach to binary chalcogenides, following our earlier analysis of elemental pnictogens and chalcogens, which allowed us to eliminate the complexities introduced by cation–anion interactions and develop a clearer understanding of the bonding mechanisms in these materials. In the present contribution, we conduct a detailed theoretical investigation of chemical bonding in selected binary chalcogenides of the A^IVX^VI and A^V₂X^VI₃ families under compression, applying the principles of the new theory of multicenter bonding. We aim to demonstrate that the EDMBs previously identified in dense pnictogens and chalcogens also occur in dense chalcogenides, including PCMs, and that the pressure-induced formation of EDMBs in these materials shows the stages observed in pnictogens and chalcogens and described in the unified theory of multicenter bonding.^9,10

Crystal structure

The A^IVX^VI semiconductors (A = Ge, Sn, Pb; X = S, Se, Te) are isoelectronic with group 15 (or V) elements (N, P, As, Sb, Bi) and display three distinct structural arrangements at RP. This structural diversity closely parallels that observed in group 15 elements, reflecting similar bonding mechanisms. The lighter compounds in this family (GeS, GeSe, SnS, and SnSe) adopt the orthorhombic GeS-type structure (space group Pnma), which serves as the binary counterpart to the black phosphorus structure (Fig. 1(a)).³¹ In the GeS-type structure, both the group 14 atoms and the chalcogen atoms typically form three primary bonds with their counterparts, creating layered arrangements. In contrast, the heavier compounds (SnTe, PbS, PbSe, and PbTe) crystallize in the rs structure (space group Fm3̄m),³² characterized by a six-fold coordination. This higher coordination reflects the increased tendency toward multicenter bonding in heavier elements.³³ Finally, GeTe adopts a rhombohedral structure at RP analogous to the A7 structure exhibited by group-15 elements such as As, Sb, and Bi.^34,35

Fig. 1 Crystal structures of GeS and As₂S₃ at ambient conditions. (a) and (b) Perspective views of the Pnma phase of GeS along the b-axis. Both Ge and S atoms are threefold coordinated with S and Ge atoms, respectively. Solid lines correspond to primary Ge–S intralayer covalent bonds of similar length (∼2.4 Å) and red dashed lines correspond to secondary Ge–S intralayer non-covalent interactions (>3.0 Å). (c) and (d) View of the monoclinic P2₁/n structure of As₂S₃ perpendicular to the b- and a-axis, respectively. The AsS₃ pyramids with short As–S bonds resulting in the threefold coordination of As and twofold coordination of S. The structure can be described as AsS spiral chains (following the As1–S1–As2–S2 sequence) interconnected by S3 atoms positioned between the chains. As atoms are threefold coordinated with S atoms (see the AsS₃ pyramids) and S atoms are twofold coordinated with As atoms. Solid and red dashed lines represent the primary short covalent and secondary long non-covalent intrachain bonds, respectively, while blue dashed lines represent the secondary long non-covalent interchain bonds.

On the other hand, the A^V₂X^VI₃ chalcogenides (A = As, Sb, Bi; X = S, Se, Te) exhibit greater structural diversity at RP. Light sesquichalcogenides show low-coordinated and low-symmetry crystalline structures, in which cations and anions exhibit threefold coordination. In particular, As₂S₃ typically crystallizes in a monoclinic phase (P2₁/n),^36,37 consisting of zigzag layers primarily stacked along the b-axis (see Fig. 1(c)). Similarly, Sb₂Se₃ adopts a zigzag layered orthorhombic structure (Pnma, see Table S1 and Fig. 1(a)). In contrast, heavy sesquichalcogenides, like Bi₂Se₃ and Bi₂Te₃, are narrow bandgap semiconductors that crystallize in a highly-symmetric layered rhombohedral (R3̄m) tetradymite-like structure composed of Te–Bi–Te–Bi–Te sublayers (stacked along the c-axis of the hexagonal unit cell) and where the central Te atoms and Bi atoms exhibit a higher (sixfold) coordination.

In all these sesquichalcogenides, layers were traditionally thought to be held together by weak van der Waals (vdW) interactions,³⁸ but recent research suggests more complex bonding mechanisms may be involved in these sesquichalcogenides and related layered materials.^39–42 In particular, these studies have revealed that the interlayer bonds in heavy sesquichalcogenides, and related layered materials are significantly stronger than simple vdW interactions. This enhanced bonding strength arises from additional interlayer charge transfer originating from unconventional intralayer bonds, consistent with the EDMB model. In particular, As₂Te₃ presents an interesting case of polymorphism, crystallizing in two distinct modifications. α-As₂Te₃ crystallizes in a quasi-layered monoclinic structure (C2/m),⁴³ while β-As₂Te₃ adopts a tetradymite-like structure, similar to Bi₂Te₃.^44,45 Both polymorphs show distinct pressure dependence of their structures that are related to their different chemical bonding.^43,46

Computational details

In this study we have performed DFT theoretical calculations as implemented in the Vienna Ab initio Simulation Package (VASP).^47–49 The electronic structure calculations employed projector augmented wave (PAW) potentials⁵⁰ with Ge (4s² 4p²), Sn (4d¹⁰ 5s² 5p²), Pb (5d¹⁰ 6s² 6p²), As (4s² 4p³), Sb (5s² 5p³), Bi (5d¹⁰ 6s² 6p³), S (3s² 3p⁴), Se (4s² 4p⁴), and Te (5s² 5p⁴) electrons treated as valence electrons. The calculations were performed with the generalized gradient approximation (GGA) of Perdew–Burke–Ernzerhof revised for solids (PBEsol).⁵¹ The wavefunction was expanded in plane waves with a kinetic energy cutoff sufficient to ensure convergence of the total energy within 10⁻⁶ eV per atom. The Brillouin-zone (BZ) sampling was performed using Γ-centered Monkhorst–Pack⁵² grids with k-point meshes optimized for each structure. All details of theoretical calculations as well as the optimized structures at RP are provided in Table S1. The computational parameters were carefully selected to balance computational efficiency with the accuracy required for reliable bonding analysis. All structural visualizations presented in this work were generated using the VESTA program,⁵³ which allows a detailed examination of atomic arrangements and bonding environments.

As discussed in the introduction, the chemical bonding in A^IVX^VI and A^V₂X^VI₃ chalcogenides has been previously analyzed through various frameworks, including conventional bonding models (covalent, ionic, and metallic)^54–56 as well as the more specialized hypervalent and metavalent bonding models.^{14–16,29,37,57} In the present study, we study the chemical bonding in these systems to provide robust evidence for the EDMB model in dense A^IVX^VI and A^V₂X^VI₃ chalcogenides. For that purpose, we combine DFT calculations with the use of the Quantum Theory of Atoms in Molecules (QTAIM) framework.⁵⁸ This approach allows us to analyze the topology of the electron density and extract quantitative descriptors of chemical bonding. Specifically, we obtain the Bader's atomic charges and the delocalization indexes (DI) between any two atoms. With the DI values we calculate the number of electrons shared (ES) between any two atoms as twice the DI, which gives an idea of the covalency of the interatomic bonds. On the other hand, we estimate the ionicity of the interatomic bonds with the normalized number of electrons transferred (ET) between two atoms.^9,10 The combination of DFT calculations with QTAIM analysis creates a powerful methodology for investigating the fundamental bonding mechanisms, particularly as they transform under pressure. This methodology was used in the two papers that describe the fundamentals of the unified theory of multicenter bonding,^9,10,30 and has allowed us to identify EDMBs in electron-rich elements, such as polyiodides,⁵⁹ AX₃ halides⁶⁰ and even AIO₃ oxide perovskites.⁶¹ In this context, it must be stressed that our unified theory of multicenter bonding is contrary to the common belief that EDMBs are not possible in electron-rich elements. In this work, we provide additional pieces of evidence of EDMB formation in electron-rich systems by analyzing the kinetic G(r), potential V(r), and total H(r) energy densities using the Thomas–Fermi approximation, incorporating the semiclassical gradient correction introduced by Kirzhnits.^62,63 These topological properties provide complementary perspectives on the nature of chemical bonding in materials, allowing us to characterize the formation and properties of EDMBs under varying conditions of pressure and composition.

To study the bonding properties of materials with a density-based method as QTAIM, the DFT charge densities, ρ(r), were calculated using the Yu–Trinkle method,⁶⁴ as implemented in CRITIC2 program.²¹ VASP output files (CHGCAR and AECCAR) were used to calculate Bader atomic charges and atomic volumes. At each pressure, the number of critical points satisfies the Morse zero-sum rule. Quantum ESPRESSO (version 6.5)⁶⁵ was utilized for this analysis, in conjunction with wannier90⁶⁶ and the CRITIC2 programs.²¹ Single-point calculations were performed at the VASP equilibrium geometries, using the same uniform k-point grids mentioned above. A plane-wave cutoff of 100 Ry and a density cutoff of 400 Ry were consistently applied. Norm-conserving pseudopotentials for the Kohn–Sham states and PAW data sets for the all-electron density were sourced from the pslibrary.⁶⁷ Delocalization index (DI) calculations, used to determine the number of electrons shared (ES) as 2 × DI, were conducted using a Wannier transformation as described in ref. 68.

The local energy densities are calculated at the bond critical points (BCPs), where ∇ρ(r) = 0, using the Kirzhnits approximation^62,63 as implemented in CRITIC2 program. The ratios of the local energy densities to the corresponding charge density yield dimensionless values that serve as reliable indicators for classifying chemical bonding. The first approach, introduced by Espinosa,⁶⁹ establishes H_b/ρ_b < 0 for shared-shell (SS) interactions, such as covalent and polar bonds, and H_b/ρ_b > 0 for closed-shell (CS) interactions, including ionic bonds, hydrogen bonds, and vdW interactions. The second one was proposed to more generically distinguish between SS interactions, characterized by G_b/ρ_b < 1, and CS interactions, where G_b/ρ_b > 1.⁷⁰ Lastly, the |V_b|/G_b ratio defines three distinct regions of chemical bonding: (a) a pure CS interaction, where |V_b|/G_b < 1, ∇²ρ_b > 0; (b) a pure SS interaction, where |V_b|/G_b > 2, ∇²ρ_b < 0; and (c) a transitional CS interaction, where 1 < |V_b|/G_b < 2, ∇²ρ_b > 0. Therefore, by analyzing these energy densities and their ratios, we can assess the pressure dependence of G_b and V_b, enabling the accurate characterization of chemical bonds and their variations under compression. Moreover, these ratios facilitate comparisons across different systems, as they are intrinsic properties, unlike extrinsic properties such as ρ(r) and its Laplacian, ∇²ρ(r).⁷¹

Results and discussions

Chemical bonding in A^IVX^VI chalcogenides

GeS. At RP, Ge and S atoms in the Pnma phase (B16) are threefold coordinated, consistent with other isoelectronic chalcogenides of group 15 (pnictogen) elements. The threefold coordination satisfies the 8 − N rule (also known as Pearson's rule) since these chalcogenides, which have an average of five valence electrons per atom, require three covalent bonds to complete their valence shell and consequently satisfy the octet rule (also known as the Lewis–Langmuir's rule). At HP, GeS—along with similar compounds GeSe, SnS, and SnSe—undergoes one or more phase transitions from the orthorhombic Pnma phase toward the orthorhombic Cmcm phase.^72,73 The Cmcm (B33 or TlI-like) phase also has a layered structure, but in this case, Ge and S atoms are primarily fivefold coordinated (sevenfold coordinated when considering secondary bonds to neighboring layers). As a result, Ge and S atoms in this phase do not conform to the 8 − N rule. The apparent violation of this rule seems to result in the violation of the octet rule because it looks as if atoms would be establishing more covalent bonds than those required to complete the octet (5 bonds × 2 valence electrons = 10 valence electrons). This apparent violation of the octet rule is not true. This can be understood if we recognize that most bonds in the Cmcm phase are not conventional simple covalent 2c-2e bonds since they contain fewer than two electrons per atomic pair, similar to what has been observed in HP phases of group 15 elements.⁹ If this interpretation holds, the octet rule itself may not be truly violated.¹⁰

We selected GeS as our primary test case because it exhibits the most distorted structure relative to the octahedral coordination present in several PCMs, such as PbS, PbSe, SnTe, and PbTe. This makes GeS an excellent reference point, as our findings can be readily extrapolated to other light A^IVX^VI compounds, like GeSe, GeTe, SnS, and SnSe showing similar structures but with less distorted octahedral arrangements at RP. Notably, beyond 9 GPa, GeS transitions to a monoclinic P2₁/m phase,⁷³ which is only a slight distortion of its orthorhombic Pnma phase at RP. Since this monoclinic phase has not been observed in GeSe, we have instead simulated the evolution of the Pnma phase of GeS up to its second-order transition towards the Cmcm phase, as reported in ref. 74. Importantly, our simulated pressure-dependent structural parameters for the Pnma and Cmcm phases of GeS up to 45 GPa (Fig. 2(a) and (b)) align well with experimentally reported data for the Pnma, P2₁/m, and Cmcm phases.⁷³ This consistency reinforces the reliability of our approach, allowing us to explore bond descriptors in greater detail and gain deeper insight into the pressure-induced changes in chemical bonding in GeS up to 45 GPa.

Fig. 2 GeS: pressure dependence of structural parameters: (a) unit-cell volume; (b) lattice parameters; and (c) Ge–S bond distances corresponding to the two nearest intralayer Ge–S bonds (axial (d₁) and equatorial (d₂)) and the nearest intralayer Ge⋯S distance (d₃). In (a), (b), and (c), lines and symbols correspond to theoretical and experimental⁷³ results, respectively. Vertical dashed lines at 20 and 36 GPa indicate the limits of three stages of EDMB formation mechanism. (d) Pressure dependence of the calculated number of electrons shared (ES) in the different Ge–S interactions. (e) Relationship between the equatorial distances: primary covalent (d₂) and secondary non-covalent (d₃) at different pressures. (f) Pressure dependence of charge density, ρ_b, at the bond critical point (BCP). Points in (d)–(f) correspond to simulated data and are plotted to show the dense number of calculated pressure points in the studied materials.

The crystalline structure of GeS features three primary short intralayer Ge–S bonds: one covalent axial bond (with distance d₁) and two covalent equatorial bonds (with distance d₂). Additionally, there are two secondary non-covalent intralayer equatorial bonds (with distance d₃), as illustrated in Fig. 1(b) and 2(c). Since A^IVX^VI chalcogenides are isoelectronic to pnictogens, the three covalent bonds in Pnma-type GeS are best described as dative or coordinative covalent bonds, where S atoms (with excess electrons relative to pnictogens) share more electrons than Ge atoms (with fewer electrons relative to pnictogens). Furthermore, these Ge–S bonds are not purely covalent but rather iono-covalent or polar covalent due to the electronegativity difference between Ge and S. At HP, we observe a progressive equalization of the equatorial short (d₂) and long (d₃) Ge–S bond lengths, while the axial short Ge–S bond (d₁) exhibits distinctly different behavior (Fig. 2(c)).

The decrease of the axial short Ge–S bond length (d₁) is the typical trend of an isolated covalent bond as in C (diamond), Si, and Ge. However, the trend toward bond length equalization between primary short iono-covalent bonds (d₂) and secondary long non-covalent equatorial bonds (d₃), and in particular the increase of the primary short iono-covalent equatorial bond (d₂), is certainly not the typical trend of a covalent bond under compression. This trend to equalization of primary and secondary bonds at HP concomitant with the anomalous increase of the short primary bond at HP has been previously observed in the pressure-induced formation of multicenter bonds (in both ERMBs and EDMBs).^9,10 Therefore, we conclude that the trend observed in GeS at HP is indicative of the pressure-induced formation of multicenter bonds and it remains to be determined whether the resulting multicenter bonds correspond to ERMBs or EDMBs.

To demonstrate that EDMBs form in dense GeS, we calculate the ES and ET values between any two atoms for the different two-center bonds present in Pnma-type GeS. Fig. 2(d) shows that the ES value of the axial Ge–S bond shows an almost flat trend with increasing pressure, while the ES value of the primary equatorial Ge–S bond decreases and that of the secondary equatorial Ge–S interaction increases concurrently with pressure until the ES values of both equatorial bonds equalize. This behavior confirms that the short axial and equatorial bonds in the HP Cmcm phase of GeS are fundamentally different, despite both short bonds are considered classical iono-covalent bonds in the Pnma phase at RP. Fig. 2(d) also confirms that all equatorial bonds in the HP phase (above 35 GPa) become equivalent. The equalization of equatorial bond distances in the HP phase of GeS is associated to a charge transfer from the primary equatorial bonds toward the secondary equatorial bonds (known as the trans influence in chemical terminology). This leads us to conclude that the equatorial Ge–S bonds in the HP phase of GeS are multicenter bonds.^9,10 Moreover, an ES value around 0.8 combined with an ET value of 0.41 at 40 GPa reflects the EDMB nature of these equatorial Ge–S bonds in the HP Cmcm phase.

In previous works, we used the 2D ES vs. ET map (Fig. 3) as a tool for classifying chemical bonds in various materials.^9,10,30,41 This map plots the number of electrons shared (ES) against the normalized number of electrons transferred (ET) between any two atoms, thus providing a comprehensive framework for understanding different bonding types. The classification of bonds in Pnma- and Cmcm-type GeS mentioned in the previous paragraph is consistent with our two-dimensional ES vs. ET map.^9,10,30 In addition, the axial Ge–S bond in the Cmcm phase can be classified as iono-covalent based on its ES (1.1) and ET (0.41) values at 40 GPa (not shown in the map), unlike the equatorial Ge–S bonds in the Cmcm phase (see green hexagon for GeS at 40 GPa in Fig. 3). This classification is similar to the different primary short Ge–S bonds in Pnma-type GeS at RP (ES ≈ 1.3, ET = 0.41, see red hexagon at 0 GPa in Fig. 3).

Fig. 3 Updated 2D map illustrating the number of electrons shared (ES) versus the normalized number of electrons transferred (ET), used for classifying chemical bonds in materials. The map highlights the orange and green regions corresponding to materials with electron-rich multicenter bonds (ERMBs) and electron-deficient multicenter bonds (EDMBs). The materials analysed in this study are represented by red and green opened hexagons corresponding to the covalent bonds and EDMBs, respectively. Closed circles represent materials analyzed in previous studies.^9,10,28

The transformation of GeS from a state with primary covalent equatorial bonds plus secondary weak interactions (Pnma phase) to a state with pure equatorial EDMBs (Cmcm phase) occurs through three distinct stages, as illustrated by the pressure dependence of the d₂ and d₃ bond distances and the bond charge densities at BCPs, ρ_b, in these equatorial bonds as shown in Fig. 2(e) and (f). This three-stage transformation mechanism parallels our recent observations of EDMB formation in dense group 15 elements.⁹ In stage 1 (0–20 GPa): both equatorial Ge–S bond distances (the strong primary covalent bond d₂ and the weak secondary non-covalent bond d₃) decrease gradually with pressure, thus leading to an increase in charge density along the two bonds. In stage 2 (20–34 GPa): d₃ continues to decrease while d₂ anomalously increases. This behavior provides clear evidence of the trans influence phenomenon occurring in stage 2,^9,10 where the strengthening of the secondary bond due to the increase in charge density affects the primary bond, which evidences a notable decrease in charge density. The trans influence process in stage 2 accounts for the multicenter character of the emerging EDMB.⁹ Finally, stage 3 (above 34 GPa) shows that the primary and secondary bonds become equivalent in the orthorhombic Cmcm phase. This stage is characterized by a quasi-linear geometry of EDMBs in the equatorial plane and by a normal decrease in bond distance and increase in bond charge density of EDMBs with increasing pressure.

In addition to the analysis of the chemical bonds in different pnictogens and chalcogens as a function of pressure using the values of ES, ET, bond lengths, and bond charge densities in ref. 9, we also characterized the formation of unconventional bonds in dense pnictogens and chalcogens with the Laplacian of the charge density, ∇²ρ_b, at the BCPs. Fig. 4(a) shows the evolution of the Laplacian of different bonds in GeS under pressure. At RP, the three Ge–S iono-covalent bonds exhibit small positive values of ∇²ρ_b, as expected for polar covalent bonds. Similarly, the secondary non-covalent bonds show small positive Laplacian values typical of weak bonds. Since all Ge–S bonds show similar pressure dependence in their ∇²ρ_b values we cannot distinguish the change in chemical bonding with ∇²ρ_b. Instead, we turn to the pressure dependence of energy density ratios for a clearer picture of the pressure-induced changes in chemical bonding. Fig. 4(b)–(d) illustrate the evolution of three energy density ratios in GeS under compression. At RP, the three initial short iono-covalent Ge–S bonds show similar but slightly different values, with a significant gap between these values and those for the two secondary equatorial non-covalent interactions. This gap decreases with increasing pressure, while the differences between the ratios of the two polar covalent bonds increase (Fig. 4(b) and (d)).

Fig. 4 GeS: pressure dependence of topological parameters along the three Ge–S distances (axial (d₁), equatorial (d₂), and secondary equatorial (d₃)) at the BCP: (a) Laplacian of charge density, ∇²ρ_b; (b) and (c) ratios of total and kinetic local energy densities, H_b and G_b, to the charge density, ρ_b, respectively; (d) |V_b|/G_b ratio. The vertical dashed line indicates the theoretical phase transition pressure between the Pnma and Cmcm phases related to the formation of EDMBs.

The different behavior of the energy densities in the different bonds in GeS with increasing pressure allows us to classify the different bonds. The Ge–S bonds at RP can be classified as polar covalent bonds since they have positive ∇²ρ_b, H_b/ρ_b < 0, G_b/ρ_b > 0, and 1 < |V_b|/G_b < 2.^69–71 Among these, the three primary short polar covalent Ge–S bonds at RP exhibit the highest absolute values of these ratios, consistent with the shared-shell (SS) nature of these interactions. In contrast, the secondary long interactions at RP show very low values of these ratios, consistent with the weak nature of closed-shell (CS) interactions.^69–71 Interestingly, this picture changes at the Pnma to Cmcm phase transition in GeS occurring around 38 GPa, as already suggested in previous simulations.⁷⁴ As already commented, a square pyramidal arrangement around Ge and S atoms is obtained in the Cmcm phase. Our calculations show a continuous increase in the absolute values of the three ratios for the secondary interactions under compression (Fig. 4(b)–(d)). In other words, secondary bonds are no longer secondary bonds in the Cmcm phase. In contrast, the covalent equatorial bonds show a flat pressure behavior in the H_b/ρ_b and |V_b|/G_b ratios up to 25 GPa, after which these values decrease due to the strong increase in G_b (Fig. 4(b)–(d)). This behavior differs markedly from that of the covalent axial bond, which shows a significant increase in the |V_b|/G_b ratio and the absolute value of H_b/ρ_b. Since Ge–S bonds can be classified as transitional CS interactions (1 < |V_b|/G_b < 2, ∇²ρ_b > 0), a higher absolute value of the H_b/ρ_b ratio indicates a more covalent and stronger bond. In this context, the axial Ge–S bond exhibits greater covalent character and strength compared to the equatorial Ge–S bonds.⁶⁹ This comparison confirms that the axial and equatorial polar covalent bonds are not equivalent when they are compressed, thus reflecting the changing nature of the equatorial bonds under compression—specifically, the formation of EDMBs at HP in the equatorial plane of the Cmcm phase.

The equatorial Ge–S EDMBs in Cmcm-type GeS can be considered to be formed by infinite linear [Ge–S] chains along two directions (a and b) in an infinite square planar lattice. These bonds are characterized by small but positive ∇²ρ_b values, high |H_b|/ρ_b and G_b/ρ_b ratios, and relatively small |V_b|/G_b values compared to the polar covalent Ge–S bonds. Based on these topological properties and their pressure dependence, we can distinguish these bonds and classify the equatorial Ge–S bonds in Cmcm-type GeS as EDMBs. Two important observations should be noted here. The first is that the |V_b|/G_b > 1 (H_b < 0) value in all bond types of Pnma-type GeS indicate that these interactions are stabilized by a local concentration of charge, despite the small positive values of ∇²ρ_b suggesting charge depletion at the interatomic surface. The same applies to the two types of bonds present in Cmcm-type GeS. The second is that EDMBs in Cmcm-type GeS exhibit small values of both ∇²ρ_b and |V_b|/G_b, contrasting with the very large negative and positive values of ∇²ρ_b and |V_b|/G_b, respectively, observed in molecules with well-known ERMBs (such as FHF⁻).⁶⁹ These results suggest that the energy density values may help distinguish between ERMBs and EDMBs in a complementary way to ES and ET values within the ES vs. ET map.

The behavior here described for GeS is expected to be similar in isostructural GeSe, SnS, and SnSe, though the pressure required to achieve EDMB formation in the Cmcm phase is much lower in these compounds because they have heavier elements than those present in GeS. Similar trends in structural parameters have been observed in all four compounds. For instance, SnSe, which is the closest compound to PCMs, undergoes the Pnma-to-Cmcm phase transition around 10 GPa, as reported in previous studies.^75–77 Our theoretical simulations of structural parameters for SnSe up to 18 GPa (Fig. S1) show good agreement with experimental results. The formation of equatorial EDMBs in the Cmcm phase from the original primary equatorial iono-covalent and secondary equatorial non-covalent bonds is demonstrated in Fig. S2 and S3 in the SI. These figures show similar features to those already discussed for GeS. The main difference is that the pressure-induced EDMB formation in SnSe proceeds in two stages rather than three stages as in GeS (Fig. S1 and S2). The reason is that SnSe at RP is already in stage 2 of the process of multicenter bond formation and not in stage 1 as GeS.

PbTe. This compound represents an important reference point in our study of EDMBs in chalcogenides. Unlike the lighter compounds in the A^IVX^VI family (GeS, GeSe, SnS, and SnSe) that require HP to increase the atomic coordination and form EDMBs, rs-PbTe naturally exhibits EDMBs at RP. In the rocksalt structure, each Pb atom is surrounded by six Te atoms in an octahedral arrangement, and similarly, each Te atom is coordinated by six Pb atoms. This sixfold coordination exceeds what would be expected from the conventional 8 − N rule, which would predict a coordination number of 3 for both Pb and Te atoms since PbTe is isoelectronic to Sb and Bi. This “hypercoordination” of atoms in rs-PbTe is a characteristic feature of materials containing multicenter bonds,^9,10 and in particular, in PbTe they correspond to EDMBs as we will show. The bonding in PbTe can be understood as a natural extension of the pressure-induced transformations we have already commented in GeS and SnSe. PbTe already possesses this unconventional bonding character at RP and not at HP as in GeS and SnSe. The existence of EDMBs in rs-PbTe at RP can be attributed to the increased electronic density in Pb and Te as compared to lighter elements. This increased electronic density promotes the formation of multicenter bonds with electron-deficient character at smaller pressures than in lighter chalcogenides.^9,10

Topological analysis of the electron density in PbTe (see Table S2 in the SI) reveals bond characteristics consistent with EDMBs at RP: relatively low ES (around 0.87) and ET (0.31) values, positive values of ∇²ρ_b at the BCPs, and energy density ratios similar to those observed for the EDMBs in the Cmcm phase of GeS and SnSe (Fig. S3). These properties place PbTe firmly in the EDMB region of our ES vs. ET classification map (Fig. 3). Interestingly, Krebs⁷⁸ offered an interpretation of the chemical bonding in group-V and group-VI elements, such as arsenic and selenium, as well as in binary IV–VI compounds like PbS, in a manner that closely aligns with the framework we present here and in our previous work⁹ (see Fig. S4). In fact, the presence of EDMBs in PbTe at RP explains many of its unique properties, including its exceptional thermoelectric performance, high Born effective charges, and large optical dielectric constant. These properties arise directly from the distinctive electronic structure associated with electron-deficient multicenter bonding.

Chemical bonding in A^V₂X^VI₃ crystals

The chemistry of A^V₂X^VI₃ chalcogenides and their behavior under compression has been extensively studied by our research group.^{37,57,79–82} In one of our previous works, we conducted an integrated experimental and theoretical study of layered monoclinic α-As₂S₃ (mineral orpiment) under applied pressure.³⁷ At that time, we reported that the coordination of As atoms increases from threefold at RP to over fivefold at pressures exceeding 20 GPa, attributing this “hypercoordination” of As atoms to the formation of metavalent bonds. We now recognize that this interpretation—which we also applied to the HP behaviour of other sesquichalcogenides, such as Sb₂S₃, Sb₂Se₃, and Bi₂S₃⁵⁷—requires revision. At the time of those studies, our understanding of multicenter bonding was less developed, and we considered a good approach the metavalent bonding model.^8,14 In the present article, we aim to reinterpret the chemical bonding in α-As₂S₃ under compression using our recently proposed unified theory of multicenter bonds in which we show that the new metavalent bond is just the century-old EDMB observed in materials with electron-rich elements.^9,10

As₂S₃. The structural and vibrational properties of α-As₂S₃ at RP and HP conditions were thoroughly discussed in a previous work.³⁷ At RP, α-As₂S₃ forms layers in which the two independent As atoms (As1 and As2) are threefold coordinated to S atoms (S1, S2, and S3) through iono-covalent bonds (of approximately 2.3 Å bond length) within the same layer (see green pyramids in Fig. 1(c)), so both As and S atoms satisfy the 8 − N and octet rules. Additionally, As1 and As2 atoms have two more distant S neighbours within the same layer (intralayer secondary interactions) at distances greater than 3.0 Å (see red dashed lines in Fig. 1(d)) and two further S neighbours in adjacent layers (interlayer secondary interactions) at distances exceeding 3.5 Å (see blue dashed lines in Fig. 1(d)). On increasing pressure, the hypercoordination of As1 and As2 atoms is reproduced in Fig. S5a and b, which are similar to those already reported in ref. 37 and it is clear how As1 atoms change from threefold to fivefold coordination (plus two) above 16 GPa, while As2 atoms change from threefold to fivefold coordination at the same pressure. In this work, we will focus on analyzing the bonding framework in α-As₂S₃ using bonding descriptors that go beyond those previously employed. In the following sections, we analyze the chemical bonding nature of the two non-equivalent As atoms (As1 and As2) at HP, using the same bonding descriptors previously discussed for GeS.

Fig. 5(a) presents the ES values along the various bonds associated with the As1 atom. For the three covalent As1–S bonds, the ES value is approximately 2 at RP, as expected for almost pure covalent bonds, while the intralayer and interlayer interactions display significantly lower ES values, as expected for non-covalent interactions. This behavior can be explained in terms of bond distances and charge densities at the BCPs (see Fig. 5(b) and (c)). As pressure increases, we observe a trans influence of the secondary intralayer As1–S3 (pink) and As1–S1 (green) bonds into the primary intralayer As1–S3 (black) and As1–S1 (blue) bonds, respectively, since these two bond pairs tend to equalize at approximately 16 GPa (Fig. S5a). Contrarily, the primary intralayer As1–S2 (red) bond and the two secondary As1–S2 interlayer interactions (open symbols), which lack neighbouring secondary and primary bonds, respectively, do not exhibit a trans influence in their distance and charge density relationships. The weak nature of the two interlayer As1–S2 interactions is confirmed by their very low ES values, which correspond to very low charge density at the BCPs (see Fig. 5(b) and (c)). In this way, we can conclude that the ES values of As1 atoms show three different bonding behaviors beyond 16 GPa: one primary covalent As1–S2 bond, two pairs (As1–S1 and As1–S3) of EDMBs, and two secondary non-covalent As1–S2 interactions. This explains the five plus two atomic coordination of As1 atoms. Similarly, the coordination of As2 increases from threefold at RP to fivefold above 16 GPa, driven by the equalization of the two As2–S1 and two As2–S2 intralayer bond distances (leading to two EDMBs) and the persistence of the covalent As2–S3 bond (see Fig. S5b in the SI). Similar bonding features regarding ES, bond distances, and charge densities already described for As1–S bonds are found for As2–S bonds (Fig. 5(d)–(f)).

Fig. 5 (a) and (d) Calculated number of electrons shared (ES) along the different bonds of As1 and As2 atoms in the monoclinic structure of α-As₂S₃ up to 34 GPa. Intralayer and interlayer interactions are shown in solid and dashed lines as well as with solid and open symbols, respectively. The relationship between the primary intralayer bonds and their corresponding secondary intralayer bonds and interlayer interactions in terms of (b), (e) distances and (c), (f) charge densities at the BCP.

As already shown for GeS, the change in chemical As–S bonding in α-As₂S₃ at HP can be traced by the changes in the topological properties (∇²ρ_b as well as H_b/ρ_b, G_b/ρ_b, and |V_b|/G_b ratios). Fig. 6 and 7 show the pressure dependence of these topological properties for the As1 and As2 bonds, respectively. At RP, the three primary, short As1–Sx (x = 1, 2, 3) and As2–Sx (x = 1, 2, 3) bonds can be classified as almost pure covalent bonds based on the following topological properties: ∇²ρ_b < 0, H_b/ρ_b < 0, 0 < G_b/ρ_b < 1, and 1 < |V_b|/G_b > 2. The two longer (secondary) intralayer bonds of As1 and As2 atoms show a non-covalent character (and larger polarity) than the primary bonds (see Fig. 6(a) and 7(a)) and evidence features, like ∇²ρ_b > 0, H_b/ρ_b ≈ 0 and 1 < |V_b|/G_b < 2, that classify them as transitional CS interactions. Finally, the two long (secondary) interlayer As1–S2 and As2–S2 bonds show characteristics of pure CS interactions, with ∇²ρ_b > 0, H_b/ρ_b > 0, G_b/ρ_b > 0, and |V_b|/G_b < 1 as illustrated in Fig. 6(b)–(d) and 7(b)–(d).

Fig. 6 As₂S₃: pressure dependence of topological parameters along the different As1–S distances at the BCP: (a) ∇²ρ_b; (b) H_b/ρ_b; (c) G_b/ρ_b; (d) |V_b|/G_b. The vertical dashed line indicates the theoretical phase transition pressure where symmetrization of bond distances occurs leading to the formation of EDMBs.

Fig. 7 As₂S₃: pressure dependence of topological parameters along the different As2–S distances at the BCP: (a) ∇²ρ_b; (b) H_b/ρ_b; (c) G_b/ρ_b; (d) |V_b|/G_b. The vertical dashed line indicates the theoretical phase transition pressure where symmetrization of bond distances occurs leading to the formation of EDMBs. Unlike in the previous figure, the values of the long interlayer non-covalent interactions are not plotted because there is no BCP along these distances. The critical points along these distances are of cage/ring type but not of bond type.

As pressure increases, significant changes occur in the bonding framework around As1 and As2 atoms: (a) each pair of intralayer covalent and non-covalent bonds converge to have similar distances and strengths, as indicated by comparable values of ∇²ρ_b, H_b/ρ_b, G_b/ρ_b, and |V_b|/G_b ratios; (b) the coordination of As1 atoms increase from threefold at RP to more than fivefold (5 + 2) above 16 GPa, while that of As2 atoms increase from three to five in the same pressure range (see Fig. S5) due to the equalization of several pairs of As–S bonds. Contrarily, the primary covalent As1–S2 and As2–S3 bonds, positioned nearly perpendicular to the square planar base formed by the other four equatorial bonds remain relatively unchanged at HP. This behaviour parallels what we observe for the axial Ge–S bond in GeS. This is consistent with the covalent nature of these bonds in the whole pressure range, as already commented with the ES values. A uniform behavior of the topological properties is also found for the two secondary interlayer non-covalent bonds that remain with increasing pressure (see dashed lines in Fig. 6 and 7). On the contrary, the EDMBs formed above 16 GPa show a distinct behaviour of the topological properties before and after 16 GPa. Noteworthy, the changes in the convergent pairs of bonds are similar to those already reported for GeS, thus confirming the EDMB nature of these newly formed bonds in α-As₂S₃ above 16 GPa.

A key observation from our analysis is the slight difference in the strength of the EDMB pairs formed around As1 and As2 atoms. The two EDMB pairs formed around As1 atoms show somewhat different ES values while those around the As2 atom have almost identical ES values (Fig. 5(a) and (d)). This distinction is also evident in the H_b/ρ_b, G_b/ρ_b, and |V_b|/G_b ratios at HP. For As1, the two pairs of EDMBs exhibit two separate trends, whereas for As2, they converge into a single region (see Fig. 6 and 7). The similar bonding characteristics of the EDMBs around the As2 atom, initially exist as two primary covalent and two intralayer interactions, are also evident in the evolution of their distances and charge density ρ_b, as shown in Fig. 5(e) and (f). The reason for the different bonding properties of EDMBs in As1 and As2 atoms comes from the different atomic environment of both atoms since As1 shows 5 + 2 atomic coordination while As2 shows pure 5 coordination.

In summary, it can be concluded that above 16 GPa, the As1 (As2) polyhedral units in mineral orpiment transition from threefold coordination to 5 + 2 (5) coordination while maintaining the same space group. When comparing the EDMB formation mechanism in α-As₂S₃ to that discussed earlier for GeS and SnSe, we observe that α-As₂S₃ exhibits only two stages of this mechanism (like SnSe and unlike GeS). Note that the covalent As–S bonds that transform into EDMBs already show an increase in bond length under compression at RP. This is consistent with our previous findings that pnictogen elements (As, Sb, and Bi) typically show a three-stage mechanism, while chalcogen elements (Se and Te) show a two-stage one under compression.⁹ Also notice that A^IVX^VI compounds are isoelectronic to pnictogens, while the average number of electrons per formula unit in A^V₂X^VI₃ compounds is much closer to chalcogens.

Sb₂Se₃. This compound crystallizes in an orthorhombic Pnma phase (Z = 4, U₂S₃-type). This structure includes two crystallographically distinct Sb sites, labeled Sb1 and Sb2, and the structure can be considered as composed by Sb1Se₇ and Sb2Se₇₊₁ coordination polyhedra.⁸³ In these distorted polyhedra, Sb1 is coordinated to 3 + 3 Se atoms, while Sb2 exhibits a 3 + 2 coordination environment if we consider bond distances below 3.4 Å (Fig. S6c and e in the SI). The pressure dependence of the volume and lattice parameters is shown in Fig. S6a and b, where a good agreement of our theoretical calculations with experimental values⁸⁴ can be appreciated. It must be noted that a second-order phase transition from Pnma to β-Pnma has been recently suggested around 12.5 GPa,⁸⁵ but we have noted as Pnma all our data up to 30 GPa. The good agreement between our theoretical calculations and experimental values⁸⁴ also extends to the pressure dependence of the Sb–Se bond distances (experimental bond distances from ref. 82 have not been plotted in Fig. S6c and e for clarity reasons in order to show the large number of pressure points calculated). Longer Sb–Se bonds are significantly compressed, while shorter ones remain relatively unchanged. Notably, the bond lengths within the Sb1Se₇ polyhedron converge between 10 and 20 GPa, indicating a rather quick pressure-induced symmetrization. This behavior reflects the diminishing stereochemical activity of the Sb³⁺ LEPs under compression, which can be considered as an indication of the EDMB formation.⁸⁴ The ES values of Sb1–Se bonds are in the range of 0.8–1.2 at HP once the symmetrization occurs. Concomitantly, the ET values of Sb1–Se bonds are in the range of 0.30–0.32. Therefore, the Sb1–Se bonds above 10 GPa can be placed in the region of EDMBs in Fig. 3. A similar pressure-induced symmetrization occurs for the Sb2–Se bonds (Fig. S6e) with the short Sb2–Se3 bond and the long Sb2–Se1 bond being EDMBs at RP and the short covalent Sb2–Se1 bond tending to become an EDMB above 10 GPa (Fig. S6f). These Sb2–Se bonds have not been plotted in Fig. 2 for the sake of clarity. Interestingly, Sb1–Se and Sb2–Se bonds in Sb₂Se₃ under compression exhibit trends in their topological parameters that closely resemble those observed for As1–S and As2–S bonds in As₂S₃. In both cases, pressure progressively drives the secondary interactions toward values comparable to those of the original covalent bonds. This result supports the emergence of EDMBs in Sb₂Se₃ above 10 GPa (Fig. S7 and S8 in the SI) as it was observed for As₂S₃ above 16 GPa.

Bi₂Se₃ and Bi₂Te₃. These compounds represent particularly significant cases of study within the A^V₂X^VI₃ chalcogenide family, as they naturally exhibit EDMBs at RP in the rhombohedral R3̄m structure (tetradymite-type). Unlike lighter compounds in this family, such as As₂S₃ and Sb₂Se₃, which require HP to form EDMBs, Bi₂Se₃ and Bi₂Te₃ provide an opportunity to study these unconventional bonds without the need for pressure application. In the R3̄m phase, Bi₂Se₃ (Bi₂Te₃) exhibits a layered structure consisting of quintuple Se (Te)–Bi–Se (Te)–Bi–Se (Te) layers stacked along the c-axis of the conventional hexagonal unit cell. Traditionally, the bonding within these quintuple layers was described as covalent, while the interactions between adjacent quintuple layers were considered weak vdW forces. However, our detailed topological analysis (see Table S2) reveals a bonding scenario that aligns with the EDMB model. In other words, the intralayer Bi–Se (Bi–Te) bonds can be understood as EDMBs. This result is in line with the claim of Wuttig and coworkers that consider that intralayer bonds in tetradymite-like compounds show metavalent (electron-deficient) bonding.⁷

More specifically, our analysis of the electron density topology in the R3̄m phase of Bi₂Se₃ and Bi₂Te₃ at RP (see Table S2) reveals several key characteristics that support the presence of intralayer EDMBs in these compounds: (1) coordination environment: each Bi atom in the R3̄m phase of Bi₂Se₃ (Bi₂Te₃) is coordinated by six Se (Te) atoms (three from each adjacent Se (Te) layer within the quintuple layer), forming an octahedral-like arrangement. The three external Se (Te) atoms of the quintuple layer are Se1(Te1) atoms and the three internal Se (Te) atoms in the quintuple layer are Se2 (Te2) atoms. This six-fold coordination exceeds what would be expected from the conventional 8 − N rule, which would predict a coordination number of 3 for Bi atoms. (2) Accordingly, there are two types of Bi–Se (Bi–Te) bonds: Bi–Se1 (Bi–Te1) and Bi–Se2 (Bi–Te2). The ES values of the Bi–Se1 (Bi–Te1) bonds are around 1.1–1.2, while the ES values of the Bi–Se2 (Bi–Te2) bonds are around 0.8–0.9. These values, taken together with the corresponding ET values 0.29 (0.19) for both Bi–Se (Bi–Te) bonds, indicate that these bonds show ES vs. ET values that are significantly lower than those expected for conventional covalent bonds, thus showing the electron deficiency of these bonds (see Fig. 3). In particular, those values allow us to locate these bonds in the region of EDMBs.^9,10 More precisely, the internal Bi–Se2 (Bi–Te2) bonds are pure EDMBs, while the external Bi–Se1 (Bi–Te1) bonds with larger ES values can be considered as a mixture between pure covalent and pure EDMBs and are in the borderline between both regions (see Fig. 3). This view is consistent with the analysis of the chemical bonds already performed on tetradymite-like SnSb₂Te₄ and PbBi₂Te₄ with septuple layers in which external bonds within the layers were found to be more covalent in nature than the internal bonds within the layers, which are clearly EDMBs.^86,87 (3) The Laplacian and energy density ratios: all Bi–Se (Bi–Te) bonds show ∇²ρ_b > 0 combined with H_b/ρ_b > 0, G_b/ρ_b > 0, and 1 < |V_b|/G_b < 2. These values place them in the transitional CS interaction category since these bonds cannot be considered either as a SS interaction as covalent bonds or a CS interaction as non-covalent bonds. This result is consistent with the EDMB model which considers that multicenter bonds are longer and softer than covalent bonds and can be considered as bonds with the shared electrons belonging partly to the covalent sphere and partly to the vdW sphere.¹⁰ (4) Interlayer bonding: the interactions between adjacent quintuple layers, traditionally described as vdW forces, show non-negligible ES values (ca. 0.3–0.4) and topological properties that indicate weak but significant bonding character. The stronger character of the interlayer interactions in tetradymite-like A^V₂X^VI₃ and related chalcogenides, and other layered materials, such as BiTeX (X = Cl, Br, I) compounds, has been ascribed to the presence of delocalized electrons in the interlayer space,³⁹ as a consequence of the presence of intralayer EDMBs.^41,87

In summary, the heavier A^IVX^VI and A^V₂X^VI₃ compounds, such as PbTe, Bi₂Se₃, and Bi₂Te₃, have crystalline structures at RP that exhibit a high atomic coordination, do not satisfy the 8 − N rule, and feature EDMBs. On the contrary, the lighter A^IVX^VI and A^V₂X^VI₃ compounds, such as GeS, SnSe, As₂S₃ and Sb₂Se₃, which show layered crystalline structures at RP that exhibit a low atomic coordination, satisfy the 8 − N rule, feature covalent bonds, and develop EDMBs at HP due to the pressure-induced equalization of intralayer bonds. Therefore, pressure is a fundamental tool that allows the formation of EDMBs in these compounds due to the increase in electronic density, thus mimicking the effect of substitution of light by heavy elements. Interestingly, the pressure for the emergence of EDMBs in SnSe and Sb₂Se₃ is similar (above 10 GPa). This is a consistent result since both compounds SnSe and Sb₂Se₃ are the A^IVX^VI and A^V₂X^VI₃ compounds with the heavier atoms showing conventional covalent bonds at RP. This means that both compounds are the materials that will develop EDMBs under compression at the smallest pressures. On the other hand, GeS and As₂S₃ tend to form EDMBs at much larger pressures due the lighter elements present in these two compounds. This trend of pressure-induced formation of EDMBs in A^IVX^VI and A^V₂X^VI₃ compounds is shown in Fig. 8. The lighter compounds of both families (GeS and As₂S₃) that are farther from the heavier compounds of both families (PbTe and Bi₂Te₃) will develop EDMBs at much larger pressures than the rest of the compounds of these two families.

Fig. 8 Pressure-induced formation of EDMBs in A^IVX^VI and A^V₂X^VI₃ compounds. Blue compounds are made of light elements and show crystalline structures at RP that satisfy the 8 − N rule and tend to form EDMBs at HP. Red compounds are made of heavy elements and crystallize in structures at RP that do not satisfy the 8 − N rule and already exhibit EDMBs at RP. These structures are the rocksalt and tetradymite structures of A^IVX^VI and A^V₂X^VI₃ compounds, respectively. Pink compounds are made of a mixture of light and heavy elements and show polymorphism at RP with a stable structure at RP satisfying the 8 − N rule and having conventional covalent bonds (rhombohedral structure in GeTe and monoclinic structure in As₂Te₃) and a metastable structure at RP that violates the 8 − N rule and feature EDMBs (rocksalt structure in GeTe and tetradymite structure in As₂Te₃). The red arrow indicates that pressure leads the compounds with light elements to behave at HP as the compounds with heavy elements at RP.

Conclusions

In this work, we have conducted a comprehensive theoretical investigation of chemical bonding in chalcogenides of the A^IVX^VI and A^V₂X^VI₃ families under compression, applying the principles of our recently proposed unified theory of multicenter bonding.^9,10 Our findings provide strong evidence that electron-deficient multicenter bonds (EDMBs) are a fundamental bonding mechanism in the binary crystalline phases of case study chalcogenides, offering a revolutionary unified framework that resolves the long-standing controversy between the previously proposed metavalent and hypervalent bonding models in the crystalline phase of chalcogenide-based phase change materials. The systematic analysis of GeS, SnSe, PbTe, As₂S₃, Sb₂Se₃, Bi₂Se₃, and Bi₂Te₃ reveals a clear pattern in the formation of EDMBs. In the lighter compounds of the two families, such as GeS and As₂S₃, EDMBs form only under high pressure through well-defined transformation mechanisms. GeS exhibits a three-stage transformation process, where primary covalent bonds and secondary interactions gradually equalize to form EDMBs at pressures above 38 GPa. In contrast, α-As₂S₃ undergoes a two-stage transformation process, achieving EDMB formation above 16 GPa. These staged processes parallel our previous observations in elemental pnictogens and chalcogens, confirming the universality of the EDMB formation mechanism. This mechanism of pressure-induced EDMB formation is also observed in other light compounds in these two families. On the contrary, the heaviest compounds of the two families, such as PbTe and Bi₂Te₃, naturally exhibit EDMBs at room pressure. The different chemical bonds found in light and heavy compounds in the two families is due to the larger electronic density in the later than in the former. This provides a coherent explanation for the exceptional properties of the heavier chalcogenides that behave as phase change materials. The sixfold coordination observed in these compounds, which exceeds predictions based on the conventional 8 − N rule, is a direct consequence of their EDMB-type bonding.

Evidence supporting the EDMB model in dense chalcogenides of the A^IVX^VI and A^V₂X^VI₃ families is provided by our topological analysis of the electron density. As in previous works, we have shown that the electron sharing (ES) values for EDMBs consistently fall between 0.8 and 1.4, significantly lower than the values of ca. 2.0 observed for conventional covalent bonds. In addition, in this work we have provided additional evidence of the formation of EDMBs using the characteristic energy density ratios (H_b/ρ_b, G_b/ρ_b, |V_b|/G_b) at the BCPs. From now on, these ratios will help us to further distinguish EDMBs from both covalent bonds and electron-rich multicenter bonds (ERMBs). Particularly revealing is the trans influence observed in the stage 2 of the EDMB formation that results in the charge transfer from primary to secondary bonds, shown by the pressure-dependence of the ES values. This charge transfer and equalization of bond distances and ES values is a clear evidence of the multicenter character of the electron-deficient bonds formed in the HP phases of A^IVX^VI and A^V₂X^VI₃ compounds, as already discussed in the unified theory of multicenter bonding.^9,10 The two-dimensional ES vs. ET map we have developed provides a powerful classification framework for chemical bonding in different molecules and solids. This map clearly distinguishes unconventional EDMBs and ERMBs from conventional covalent, ionic, and metallic bonds. Moreover, it allows us to see the evolution of chemical bonds on element substitution and on increasing pressure; i.e., on changing the electronic density.^9,10 This unified approach bridges the gap between density-based and orbital-based methods, offering a more comprehensive understanding of bonding in these materials.

A notable structural characteristic across all studied compounds whose structure shows covalent bonds at room pressure is the persistence of one directional covalent bond even at high pressure, while other bonds transform into EDMBs. This mixed bonding character—combining both directional covalent bonds and directional delocalized multicenter bonds—contributes to the unique anisotropic properties of these materials at high pressure, including their exceptional thermoelectric performance, high Born effective charges, and large optical dielectric constants. These exceptional properties are somewhat different to those of phase change materials at room pressure since these compounds already show EDMBs at room pressure and exhibit isotropic properties, e.g., cubic rs-PbTe shows isotropic properties whereas Cmcm-SnSe shows anisotropic properties.

The implications of our findings extend beyond theoretical understanding. By establishing clear structure–property relationships, our work provides a foundation for designing new materials with tailored properties. The connection between EDMBs and functional properties, such as thermoelectric efficiency and topological insulator behavior, can open avenues for materials engineering through controlled bond modification. Our study also validates pressure as an invaluable investigative tool for understanding chemical bonding. By enabling the observation of continuous bond transformation processes, pressure effectively mimics the transition from lighter to heavier elements, providing insights that would be difficult to obtain through other means.

In summary, the EDMB model for phase change materials offers a unified framework for understanding the exceptional properties of chalcogenides, resolving the controversy between competing bonding models. This work extends the study of elemental materials under compression⁹ to binary compounds related to phase change materials, providing a more comprehensive understanding of chemical bonding in chalcogenides beyond conventional categories. Future research directions include extending this framework to other material families and further investigating the relationship between EDMBs and functional properties, with the ultimate goal of designing new materials with enhanced performance for technological applications.

Conflicts of interest

The authors declare no conflict of interest.

Data availability

The data supporting this article has been included as part of the SI.

Theoretical calculation methods are available in the SI linked to this article. The authors have cited additional references within the SI. See DOI: https://doi.org/10.1039/d5tc02328a

Acknowledgements

This publication is financed by the Spanish Ministerio de Ciencia e Innovación and the Agencia Estatal de Investigación MCIN/AEI/10.13039/501100011033 as part of the project MALTA Consolider Team network (RED2022-134388-T), and I+D+ i projects PID2019-106383GB-42/43, and PGC2021-125518NB-I00, and by project PROMETEO CIPROM/2021/075 (GREENMAT) financed by Generalitat Valenciana. This study also forms part of the Advanced Materials programme supported by MCIN with funding from European Union NextGenerationEU (PRTR-C17.I1) and by Generalitat Valenciana through project MFA/2022/025 (ARCANGEL).

References

1. Y. Pei, A. LaLonde, S. Iwanaga and G. J. Snyder, Energy Environ. Sci., 2011, 4, 2085.
2. H. Wang, Y. Pei, A. D. LaLonde and G. J. Snyder, Adv. Mater., 2011, 23, 1366–1370.
3. D. Wu, L.-D. Zhao, S. Hao, Q. Jiang, F. Zheng, J. W. Doak, H. Wu, H. Chi, Y. Gelbstein, C. Uher, C. Wolverton, M. Kanatzidis and J. He, J. Am. Chem. Soc., 2014, 136, 11412–11419.
4. X. Yan, B. Poudel, Y. Ma, W. S. Liu, G. Joshi, H. Wang, Y. Lan, D. Wang, G. Chen and Z. F. Ren, Nano Lett., 2010, 10, 3373–3378.
5. Y. L. Chen, J. G. Analytis, J. H. Chu, Z. K. Liu, S. K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain and Z.-X. Shen, Science, 2009, 325, 178–181.
6. H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang and S.-C. Zhang, Nat. Phys., 2009, 5, 438–442.
7. Y. Cheng, O. Cojocaru-Mirédin, J. Keutgen, Y. Yu, M. Küpers, M. Schumacher, P. Golub, J.-Y. Raty, R. Dronskowski and M. Wuttig, Adv. Mater., 2019, 31, 1904316.
8. M. Wuttig, V. L. Deringer, X. Gonze, C. Bichara and J.-Y. Raty, Adv. Mater., 2018, 30, 1803777.
9. H. H. Osman, A. Otero-de-la-Roza, P. Rodríguez-Hernández, A. Muñoz and F. J. Manjón, J. Mater. Chem. C, 2024, 12, 10447–10474.
10. H. H. Osman, P. Rodríguez-Hernández, A. Muñoz and F. J. Manjón, J. Mater. Chem. C, 2025, 13, 3774–3803.
11. P. B. Littlewood, Crit. Rev. Solid State Mater. Sci., 1983, 11, 229–285.
12. J.-P. Gaspard, A. Pellegatti, F. Marinelli and C. Bichara, Philos. Mag. B, 1998, 77, 727–744.
13. J.-P. Gaspard, C. R. Phys., 2016, 17, 389–405.
14. Y. Cheng, S. Wahl and M. Wuttig, Phys. Status Solidi RRL, 2021, 15, 2000482.
15. M. Xu, S. Jakobs, R. Mazzarello, J.-Y. Cho, Z. Yang, H. Hollermann, D. Shang, X. Miao, Z. Yu, L. Wang and M. Wuttig, J. Phys. Chem. C, 2017, 121, 25447–25454.
16. M. Cagnoni, D. Führen and M. Wuttig, Adv. Mater., 2018, 30, 1801787.
17. Y. Yu, M. Cagnoni, O. Cojocaru-Mirédin and M. Wuttig, Adv. Funct. Mater., 2020, 30, 1904862.
18. J. Hempelmann, P. C. Müller, P. M. Konze, R. P. Stoffel, S. Steinberg and R. Dronskowski, Adv. Mater., 2021, 33, 2100163.
19. P. C. Müller, C. Ertural, J. Hempelmann and R. Dronskowski, J. Phys. Chem. C, 2021, 125, 7959–7970.
20. J. Hempelmann, P. C. Müller, C. Ertural and R. Dronskowski, Angew. Chem., Int. Ed., 2022, 61, e202115778.
21. A. Otero-de-la-Roza, E. R. Johnson and V. Luaña, Comput. Phys. Commun., 2014, 185, 1007–1018.
22. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2013, 34, 2557–2567.
23. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37, 1030–1035.
24. R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier and R. Dronskowski, J. Comput. Chem., 2020, 41, 1931–1940.
25. J. C. Golden, V. Ho and V. Lubchenko, J. Chem. Phys., 2017, 146, 174502.
26. T. H. Lee and S. Elliott, Proc. Natl. Acad. Sci. U. S. A., 2024, 121(19), e2403308121.
27. J.-Y. Raty, C. Bichara, C.-F. Schön, C. Gatti and M. Wuttig, Proc. Natl. Acad. Sci. U. S. A., 2024, 121(19), e2405294121.
28. M. Wuttig, C. Schön, D. Kim, P. Golub, C. Gatti, J. Raty, B. J. Kooi, Á. M. Pendás, R. Arora and U. Waghmare, Adv. Sci., 2023, 35, 2308578.
29. R. O. Jones, J. Phys.: Condens. Matter, 2022, 34, 343001.
30. F. J. Manjón, H. H. Osman, M. Savastano and Á. Vegas, Materials, 2024, 17, 2840.
31. H. Wiedemeier, H. Georg and G. von Schnering, Z. Kristallogr. – Cryst. Mater., 1978, 148, 295–304.
32. R. W. G. Wyckoff, Crystal Structures, Interscience Publishers, New York, 2nd edn, 1963, vol. 1.
33. W. Grochala, R. Hoffmann, J. Feng and N. W. Ashcroft, Angew. Chem., Int. Ed., 2007, 46, 3620–3642.
34. T. Chattopadhyay, J. X. Boucherle and H. G. vonSchnering, J. Phys. C: Solid State Phys., 1987, 20, 1431–1440.
35. P. Bauer Pereira, I. Sergueev, S. Gorsse, J. Dadda, E. Müller and R. P. Hermann, Phys. Status Solidi B, 2013, 250, 1300–1307.
36. N. Morimoto, Mineral. J., 1954, 1, 160–169.
37. V. P. Cuenca-Gotor, J. Á. Sans, O. Gomis, A. Mujica, S. Radescu, A. Muñoz, P. Rodríguez-Hernández, E. L. da Silva, C. Popescu, J. Ibañez, R. Vilaplana and F. J. Manjón, Phys. Chem. Chem. Phys., 2020, 22, 3352–3369.
38. G. V. Gibbs, A. F. Wallace, R. Zallen, R. T. Downs, N. L. Ross, D. F. Cox and K. M. Rosso, J. Phys. Chem. A, 2010, 114, 6550–6557.
39. R. Wang, F. R. L. Lange, S. Cecchi, M. Hanke, M. Wuttig and R. Calarco, Adv. Funct. Mater., 2018, 28, 1705901.
40. R. Vilaplana, S. Gallego-Parra, E. Lora da Silva, D. Martínez-García, G. Delaizir, A. Muñoz, P. Rodríguez-Hernández, V. P. Cuenca-Gotor, J. A. Sans, C. Popescu, A. Piarristeguy and F. J. Manjón, J. Mater. Chem. C, 2023, 11, 1037–1055.
41. E. Bandiello, H. H. Osman, J. A. Sans, P. Rodríguez-Hernández, A. Muñoz, J. González-Platas, C. Popescu, E. Greenberg, V. Prakapenka, C. Drasar, A. H. Romero and F. J. Manjón, J. Mater. Chem. C, 2024, 12, 18660–18675.
42. A. Adam, Mater. Res. Bull., 2007, 42, 1986–1994.
43. G. Carron, Acta Crystallogr., 1963, 16, 338–343.
44. V. A. Kirkinskij and V. G. Yakushev, Izv. Akad. Nauk SSSR, Neorg. Mater., 1974, 10, 1431–1435.
45. H. W. Shu, S. Jaulmes and J. Flahaut, J. Solid State Chem., 1988, 74, 277–286.
46. V. P. Cuenca-Gotor, J. A. Sans, J. Ibáñez, C. Popescu, O. Gomis, R. Vilaplana, F. J. Manjón, A. Leonardo, E. Sagasta, A. Suárez-Alcubilla, I. G. Gurtubay, M. Mollar and A. Bergara, J. Phys. Chem. C, 2016, 120, 19340–19352.
47. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
48. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
49. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
50. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
51. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
52. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
53. K. Momma and F. Izumi, J. Appl. Cryst., 2011, 44, 1272–1276.
54. P. Enders, Phys. Status Solidi B, 1985, 129, 89–99.
55. C. Bellin, A. Pawbake, L. Paulatto, K. Béneut, J. Biscaras, C. Narayana, A. Polian, D. J. Late and A. Shukla, Phys. Rev. Lett., 2020, 125, 145301.
56. S. Maier, S. Steinberg, Y. Cheng, C.-F. Schön, M. Schumacher, R. Mazzarello, P. Golub, R. Nelson, O. Cojocaru-Mirédin, J.-Y. Raty and M. Wuttig, Adv. Mater., 2020, 32, 2005533.
57. E. L. da Silva, M. C. Santos, P. Rodríguez-Hernández, A. Muñoz and F. J. Manjón, Crystals, 2023, 13, 498.
58. R. F. W. Bader, Atoms in molecules: a quantum theory, Clarendon Press, Oxford, 1990.
59. M. Savastano, H. H. Osman, Á. Vegas and F. J. Manjón, Chem. Commun., 2024, 60, 12677–12689.
60. E. Bandiello, A. Lobato, F. Izquierdo, H. H. Osman and F. J. Manjón, J. Mater. Chem. C, 2025 , in revision.
61. H. H. Osman, J. L. Rodrigo-Ramón, S. Ullah, E. Bandiello, D. Errandonea, Ó. Gomis, T. García-Sánchez, P. Botella, R. Oliva, P. Rodríguez-Hernández, A. Muñoz, C. Popescu, F. G. Alabarse and F. J. Manjón, Chem. Mater., 2025, 37(11), 4187–4202.
62. D. A. Kirzhnits, Sov. Phys. JETP, 1957, 5, 1.
63. D. A. Kirzhnits and J. L. Gammel, Am. J. Phys., 1967, 35, 1166–1167.
64. M. Yu and D. R. Trinkle, J. Chem. Phys., 2011, 134, 64111.
65. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
66. A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt and N. Marzari, Comput. Phys. Commun., 2008, 178, 685–699.
67. A. Dal Corso, Comput. Mater. Sci., 2014, 95, 337–350.
68. A. Otero-de-la-Roza, Á. Martín Pendás and E. R. Johnson, J. Chem. Theory Comput., 2018, 14, 4699–4710.
69. E. Espinosa, I. Alkorta, J. Elguero and E. Molins, J. Chem. Phys., 2002, 117, 5529–5542.
70. R. F. W. Bader and H. Essén, J. Chem. Phys., 1984, 80, 1943–1960.
71. Gatti Carlo, Z. Kristallogr., 2005, 220, 399–457.
72. T. Chattopadhyay, A. Werner, H. G. von Schnering and J. Pannetier, Rev. Phys. Appl., 1984, 19, 807–813.
73. R. P. Dias, M. Kim and C.-S. Yoo, Phys. Rev. B, 2016, 93, 104107.
74. M. Durandurdu, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 144106.
75. I. Loa, R. J. Husband, R. A. Downie, S. R. Popuri and J.-W. G. Bos, J. Phys.: Condens. Matter, 2015, 27, 072202.
76. X. Chen, P. Lu, X. Wang, Y. Zhou, C. An, Y. Zhou, C. Xian, H. Gao, Z. Guo, C. Park, B. Hou, K. Peng, X. Zhou, J. Sun, Y. Xiong, Z. Yang, D. Xing and Y. Zhang, Phys. Rev. B, 2017, 96, 165123.
77. I. Efthimiopoulos, M. Berg, A. Bande, L. Puskar, E. Ritter, W. Xu, A. Marcelli, M. Ortolani, M. Harms, J. Müller, S. Speziale, M. Koch-Müller, Y. Liu, L.-D. Zhao and U. Schade, Phys. Chem. Chem. Phys., 2019, 21, 8663–8678.
78. H. Krebs, Fundamentals of inorganic crystal chemistry, McGraw-Hill, London, 1968.
79. R. Vilaplana, D. Santamaría-Pérez, O. Gomis, F. J. Manjón, J. González, A. Segura, A. Muñoz, P. Rodríguez-Hernández, E. Pérez-González, V. Marín-Borrás, V. Muñoz-Sanjose, C. Drasar and V. Kucek, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 184110.
80. R. Vilaplana, O. Gomis, F. J. Manjón, A. Segura, E. Pérez-González, P. Rodríguez-Hernández, A. Muñoz, J. González, V. Marín-Borrás, V. Muñoz-Sanjosé, C. Drasar and V. Kucek, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 104112.
81. A. Segura, V. Panchal, J. F. Sánchez-Royo, V. Marín-Borrás, V. Muñoz-Sanjosé, P. Rodríguez-Hernández, A. Muñoz, E. Pérez-González, F. J. Manjón and J. González, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 195139.
82. F. J. Manjón, R. Vilaplana, O. Gomis, E. Pérez-González, D. Santamaría-Pérez, V. Marín-Borrás, A. Segura, J. González, P. Rodríguez-Hernández, A. Muñoz, C. Drasar, V. Kucek and V. Muñoz-Sanjosé, Phys. Status Solidi B, 2013, 250, 669–676.
83. G. P. Voutsas, A. G. Papazoglou, P. J. Rentzeperis and D. Siapkas, Z. Kristallogr. – Cryst. Mater., 1985, 171, 261–268.
84. I. Efthimiopoulos, J. Zhang, M. Kucway, C. Park, R. C. Ewing and Y. Wang, Sci. Rep., 2013, 3, 2665.
85. H. Cheng, J. Zhang, P. Yu, C. Gu, X. Ren, C. Lin, X. Li, Y. Zhao, S. Wang and Y. Li, J. Phys. Chem. C, 2020, 124, 3421–3428.
86. J. A. Sans, R. Vilaplana, E. L. da Silva, C. Popescu, V. P. Cuenca-Gotor, A. Andrada-Chacón, J. Sánchez-Benitez, O. Gomis, A. L. J. Pereira, P. Rodríguez-Hernández, A. Muñoz, D. Daisenberger, B. García-Domene, A. Segura, D. Errandonea, R. S. Kumar, O. Oeckler, P. Urban, J. Contreras-García and F. J. Manjón, Inorg. Chem., 2020, 59, 9900–9918.
87. T. Garcia-Sanchez, V. P. Cuenca-Gotor, H. H. Osman, A. Muñoz, P. Rodriguez-Hernandez, J. Á. Sans, S. Gallego-Parra, C. Popescu, O. Oeckler, D. Errandonea and F. J. Manjón, J. Alloys Compd., 2025, 1010, 177010.

---