Introduction

Predicting the ground-state structure of crystalline materials, initially thought to be an unsolvable problem, became an active area of research with the advent of efficient numerical implementation of computational total energy methods. Various approaches to exploring potential energy surface of solids (PES) from first principles (ab initio thermodynamics) have been developed,1,2,3 leading to exciting discoveries such as superconducting dense hydrogen sulfide,4 new and intriguing forms of matter at elevated pressures,5 new functional materials,6,7,8 and other.9,10 Beyond the ground-state structures, efforts in exploiting polymorphism and extending structure prediction and materials by design strategies to metastable systems have also been pursued.11,12,13,14 Describing accurately the PES was proven essential also in the efforts to predict from first principles new thermodynamically stable topological insulators15,16 and Weyl-Kondo semimetals.17 The space of strongly correlated electron systems, on the other hand, represents a virtually untapped territory for finding new materials exhibiting potentially ground-breaking physical properties. However, exploring the PES of these materials (usually d- or f-electron systems) poses significant challenges.18,19 This includes both the “choice” of the ground-state structure and the energy ordering of different PES local minima.

The fact that strong electron correlations can influence the ground-state crystal structure has been demonstrated previously.20,21,22,23,24,25,26 For example, including correlations at an appropriate level of theory (Random Phase Approximation and Quantum Monte Carlo) was shown to be critical in reproducing the known rocksalt ground-state structure of MnO.21,22,23 Contrary to experiments, correlation-deficient approaches such as the classic approximations to density functional theory (DFT)27,28 would suggest zincblende or wurtzite structure, both featuring tetrahedrally coordinated atoms, to be lower in energy than the octahedrally coordinated rocksalt. However, the general physical mechanism through which electron correlations influence relative energies of different crystal structures remains elusive, and a systematic investigation of the role of correlations in determining the thermodynamically stable crystal structure of transition metal compounds is presently missing.

In order to investigate these fundamental questions, in this paper we study the influence of strong electronic correlations present in six transition metal binary oxides and chalcogenides (CrO, MnO, FeO, CoO, CoS, and CoSe) in four common crystal structure types shown in Fig. 1 (rocksalt, NiAs-type, zincblende, and wurtzite). This particular selection of crystal structures covers both the change in local coordination of the atoms and in the long-range order, because it combines octahedral coordination with cubic symmetry (rocksalt), octahedral coordination with hexagonal symmetry (NiAs-type), tetrahedral coordination with cubic symmetry (zincblende), and tetrahedral coordination with hexagonal symmetry (wurtzite). Furthermore, the selected set includes the experimental ground-state structure for all studied systems. To perform the necessary quantum simulations we combine the local density approximation (LDA)29 with the rotationally-invariant slave-boson mean-field theory (RISB)30,31,32,33 (Supplemental Material), which is equivalent to the multi-orbital Gutzwiller approximation24,34,35,36,37 and includes electron correlations beyond the single-particle picture, providing us with predictions typically almost as accurate as dynamical mean-field theory (DMFT)38,39,40—while being orders of magnitude less computationally demanding.

Fig. 1
figure 1

Representation of generic energy profile as a function of the crystal configuration. Crystal structures considered in this work: NiAs-type, rocksalt, wurtzite, and zincblende

The main findings of our work are the following: (1) The strong electron correlations influence dramatically many important features of the PES in all d-electron materials considered, such as the energy ordering of different polymorphs and the thermodynamically stable crystal structure. (2) Available theories describing the strong electron correlations beyond a mean-field single-particle picture, see, e.g., refs. 24,32,34,35,37,38,39,40, provide us with effective tools for simulation-based structure prediction studies of d-electron materials. (3) The influence of Mott localization on charge transfer is the key physical mechanism determining structural stability in all of the Mott systems considered, while d-electron covalency effects are also essential for predicting structural stability in all metallic systems.

By employing LDA + RISB, we were able to reproduce the experimentally known ground-state structure of all compounds, as well as to uncover the dominant physical mechanisms by which strong correlations influence the energy ordering of different crystal structures in these systems.

Results and discussion

We have performed DFT, DFT + U41 and DFT + RISB30,31,32 ground-state calculations of CrO, MnO, FeO, CoO, CoS, and CoSe in four different crystal structures (rocksalt, NiAs-type, zincblende and wurtzite, see Fig. 1). The LDA + RISB calculations have been performed assuming a Hund’s coupling constant strength J = 0.9 eV and scanning different values of Hubbard interaction strength U. Since all of these transition metal compounds are paramagnetic at room temperature, the LDA + RISB simulations were performed assuming from the onset paramagnetic solutions.

In Fig. 2 is provided a bird’s eye view of the main properties of the materials considered, inherent in their zero-temperature thermodynamically stable phase, i.e., (1) the crystal configuration, (2) the equilibrium density, and (3) whether the system is a metal or an insulator. The theoretical results (triangles) are shown in comparison with the experiments (circles). The reported experimental data were obtained from refs. 26,42,43,44,45,46. The crystal structures are color coded by blue, green, red, yellow, for NiAs-type, rocksalt, wurtzite, and zincblende, respectively. The insulating phases are indicated by half-filled symbols, while the metallic phases are indicated by fully filled symbols.

Fig. 2
figure 2

Theoretical (triangles) unit cell volumes and crystal structures in comparison with the experiments (circles). The colors, blue, green, red, and yellow, correspond to NiAs-type, rocksalt, wurtzite, and zincblende structures, respectively. The metallic and insulating solution are labeled by filled and half-filled symbols, respectively

From the experimental data we observe that all of the oxides considered favor the rocksalt structure, while the thermodynamically stable lattice configuration of CoS and CoSe is NiAs-type. All of the materials have a metallic ground-state except MnO, FeO and CoO, which are Mott insulators.

The LDA fails to reproduce the experimental crystal structure (rocksalt) for all of the oxides considered.47,48 The method is also unable to capture the fact that MnO, FeO and CoO are insulators, and the predicted equilibrium volumes are generally inaccurate. Inclusion of spin polarization (ordering) at the level of the generalized gradient approximation (GGA) to DFT49 plus U (GGA + U) or LDA + U straightens out the limitations of unpolarized LDA only in part. In particular, as shown in the Supplemental Material, the spin polarized GGA + U suggests different ground-state structures of both CoS and CoSe depending on the value of U, never reproducing the experimentally known NiAs structure as the lowest energy one. The lowest energy structure changes from the rocksalt derived (U = 0 eV) to zincblende (U = 6 eV) and wurtzite (U = 12 eV). Furthermore, calculated equilibrium volumes are very different from the experimental ones. Hence, the experimentally observed change in the long-range order from the rocksalt to NiAs-type structure, that preserves octahedral coordination of atoms and occurs as the anion is replaced from CoO to CoS and CoSe, is overall poorly described by spin polarized GGA + U (or LDA + U). These results follow from an exhaustive enumeration of different spin configurations constructed on all symmetry inequivalent supercells with up to four formula units (~700 calculations for both CoS and CoSe). The spin polarized GGA + U calculations were performed following the approach described in ref. 50, which employs PAW treatment of valence electrons51 and the Perdew-Burke-Ernzerhof (PBE) functional form for the exchange-correlation functional within the GGA approximation to DFT,49 as implemented in the VASP code.52 These results constitute unequivocal evidence of the fact that the strong electron correlations influence substantially the behavior of both the electronic structure and the total energy of these materials.

Remarkably, the LDA + RISB theory provides results in very good quantitative agreement with the experiments for all six transition metal binary oxides and chalcogenides considered, simultaneously. The simulations are particularly accurate for U = 13 eV. In fact, for this value of the Hubbard interaction strength, the method captures, at the same time, all of the physical properties examined, including the crystal structure and the insulating nature of MnO, FeO and CoO. Furthermore, the experimental equilibrium volumes of all materials are reproduced within 4% error. We note that the method captures the correct crystal structure of all materials also for U = 8 eV, although the overall accuracy of the results is not as satisfactory as for U = 13 eV. In particular, the equilibrium volume of FeO and CoO is underestimated for this value of U. The reason why varying the value of U influences considerably the equilibrium volume for these two materials is that the MIT of these two systems, which occurs at U 13 eV, is first order. Consequently, the equilibrium volume evolves discontinuously at the critical point.53,54 Note that the values of the MIT critical U of FeO and CoO reported above are likely overestimated (this is a well-known systematic limitation of the RISB approximation).

These results indicate clearly that taking into account the strong electron correlations beyond the single-particle picture, e.g., utilizing the LDA + RISB approach, results in a remarkably effective tool for structure prediction work of d-electron materials.

In order to investigate in further detail the influence on the crystal structure of the strong d-electron correlations and assess its relevance for the prediction of new polymorphs, it is also interesting to inspect the energy profiles of the thermodynamically unstable crystal configurations. For illustration purposes, the behavior of the theoretical total energy as a function of the volume is shown in Fig. 3 only for CoS and MnO, both in LDA and in LDA + RISB. The analysis of the other materials is reported in the Supplemental Material.

Fig. 3
figure 3

LDA and LDA + RISB behavior of the energy as a function of the volume for CoS a, b and MnO c, d. The blue circle, green triangle, red square, and yellow diamond curves correspond to NiAs-type, rocksalt, wurtzite, and zincblende structures, respectively. The experimental equilibrium volumes for each compounds are marked by vertical solid lines, while the theoretical equilibrium volumes are marked by in dashed lines

Interestingly, the RISB correction modifies dramatically the LDA energy order of the crystal structures, both for CoS (which is a metal) and MnO (which is a Mott insulator). We note also that the LDA + RISB relative energies between the different crystal structures of both MnO and CoS, evaluated at their respective experimental equilibrium volumes, are considerably larger with respect to the Nèel temperature of the materials examined, that are all paramagnetic at room temperature. In particular, the energy difference between the octahedrally coordinated configurations (Rocksalt and NiAs-type) and the tetrahedrally coordinated configurations (Wurzite and Zincblende) is of the order of ~1 eV for both of these materials.

The critical role of the strong electron correlations in determining the crystal structure is well exemplified by the calculations of MnO displayed in panel (d) of Fig. 3. In fact, MnO is a Mott insulator in the octahedrally coordinated structures, while it is a metal in the tetrahedrally coordinated structures. As shown in the Supplemental Material, the same behavior is observed in FeO and CoO, which (at their respective experimental equilibrium volumes) are also Mott insulators in their thermodynamically stable crystal configuration, while they are metals in the tetrahedrally coordinated structures. These results indicate clearly that computing the behavior of the PES requires to take into account the subtle competing mechanisms underlying the strongly correlated regime around the Mott transition—which can only be accomplished by many-body techniques able to take into account the strong electron correlations beyond a mean-field single-particle picture, such as DMFT38,39,40 and the RISB.30,31,32 In this respect, we note also that the interplay between crystal structure and d-electron correlations in the materials examined here is considerably more complex with respect to f-electron systems such as elemental Pr and Pu, where the RISB correction to the total energy was shown to be very similar for all phases.24

It is also interesting that turning on the d-electron Hubbard interaction has a qualitatively distinct effect on the relative energies between different crystal structures with respect to increasing the volume. In fact, for all of the materials considered, the LDA Wurtzite and Zincblende energies tend to become smaller than the Rocksalt and NiAs-type energies at larger volumes. Instead, turning on the interaction parameters (U, J) has the opposite effect, i.e., it favors energetically the Rocksalt and NiAs-type structures. This implies that the pronounced physical effects of the strong electron correlations observed on the behavior of the total energy, which are accounted for by the LDA + RISB theory, are inherently “orbitally-selective”; in the sense that their influence on the crystal structure cannot be simply traced back into a uniform renormalization of the bandwidth, as one might naively expect. As we are going to show, in fact, one of the main physical effects at play is charge transfer, which is a multi-orbital phenomenon decisively influenced by the strong electron correlations.

In order to investigate the connection between Mott physics and the energy order of the crystal structures, it is insightful to take a step back and analyze the energy-ordering problem from the point of view of a simple classic point-ion electrostatic (PIE) model. The PIE model and the Madelung energy were successful in explaining the structure and the order-disorder phenomena in spinels55,56 as well as in non-isovalent perovskite alloys.57 Here we employ the PIE model and Ewald summation described in ref. 55, and evaluate the Madelung energies at the experimental equilibrium volumes using as inputs the theoretical local d-electron occupations nd reported in Table 1, which we calculated as outlined in the Supplemental Material. The point charges used to compute the Madelung energy are (a) the transition metal sites are assumed to carry a positive charge equal to the total number of valence electrons minus the value of nd, and (b) the anion sites are assumed to have a negative charge such as to make the system charge-neutral.

Table 1 LDA and LDA + RISB d-electron occupations nd computed at the experimental equilibrium volumes

Interestingly, the values of nd are considerably influenced by the strong d-electron correlations, especially for the Mott insulating phases of MnO, FeO and CoO, which are realized in the NiAs-type and Rocksalt crystal structures. For these materials, we find that the simple PIE electrostatic model supplemented by the LDA + RISB occupations is sufficient to capture the correct LDA + RISB energy order between the different crystal structures at their respective experimental volumes, while this is not the case if the LDA occupations are used, see Fig. 4. This analysis enables us to prove that the tendency of electron correlations to favor energetically the Rocksalt and NiAs-type structures is a consequence of Mott physics and electrostatics. In fact, the key mechanism at play is that the metallic phases (which are tetrahedrally coordinated) display smaller charge transfer from the cation to the anion with respect to the Mott insulators (which are octahedrally coordinated). Note that, not surprisingly, the PIE model is unable to capture the energy order of the metallic systems CrO, CoS and CoSe (not shown), as the d electrons retain a significant mixed-valence character in these systems (see Supplemental Material), and the covalent effects are very important.

Fig. 4
figure 4

Energy differences between different crystal structures, computed within LDA + RISB and the PIE model supplemented by the LDA + RISB d occupations (LDA + RISB + PIE) and the LDA d occupations (LDA + PIE). The results are shown for the Mott insulators MnO, FeO and CoO, and are computed at their respective experimental equilibrium volumes. The Rocksalt LDA + RISB + PIE and LDA + RISB energies are both conventionally assumed to be 0 in all materials considered

Our results indicate that the main reason why the strong d-electron correlations affect dramatically the crystal structure of all Mott systems considered (MnO, FeO and CoO) is that they influence substantially the charge transfer mechanism in these systems, which, in turn, controls the electrostatic interactions. As we have shown, the energy order is considerably influenced by the d-electron correlations also in the metallic systems examined. On the other hand, the qualitative features of the chemical bonds in these metals cannot be captured by electrostatics alone, as the covalent effects are also very important. In this respect, it is also interesting to note that for metallic systems (CrO, CoS and CoSe) the energy differences between crystal structures are much smaller than in the case of Mott insulators, implying a more subtle balance between electrostatic and covalent effects—which is accurately captured by a proper treatment of strong correlations.

In summary, we have analyzed theoretically the “choice” of the ground-state structure of six transition metal binary oxides and chalcogenides among four different crystal structures, utilizing LDA in combination with the rotationally-invariant slave-boson (RISB) theory, finding very good quantitative agreement with the experiments. These simulations demonstrated that the subtle competing mechanisms underlying the strongly correlated regime influence dramatically the behavior of the PES—and, in particular, the energy order between different crystal configurations—in all of the transition metal compound considered. Remarkably, for all Mott systems considered we found that the main qualitative features inherent in the influence of the strong electron correlations on the crystal structure can be understood in terms of a simple electrostatic model based on the sole knowledge of the d-electron occupations. In particular, these results indicate that one of the key physical effects to be simulated accurately for predicting the energy ordering of these materials is the variation of charge transfer induced by the strong d-electron correlations around the Mott point, which is dramatically influenced by strong electron correlations. On the other hand, we have shown that electrostatics alone cannot capture the energy order of the metals considered, where the d electrons retain a significant mixed-valence character—which indicates that covalent effects influence dramatically the chemical bonds in these systems.

From the computational standpoint, these findings suggest that simulation-based structure prediction work of transition metal compounds requires theoretical tools able to describe simultaneously: (i) the details of the bands structure, and (ii) the d-electron correlations around the Mott transition—which can influence substantially the electronic structure and the behavior of the PES in these materials. For this purpose, taking into account the strong electron correlations beyond single-particle mean-field schemes, such as LDA and GGA + U, is an absolute necessity. In this respect, the LDA + RISB technique is particularly appealing, as it is both sufficiently accurate and computationally convenient to be applied to high-throughput computational materials design.

Another important conclusion arising from our study is that the typical energy differences between the crystal structures considered, which have been all computed assuming paramagnetic solutions from the onset, are generally much larger with respect energy scales characterizing magnetism in all of the materials considered. This observation indicates that magnetic order has a negligible effect on the total energy with respect to the d-electron atomic scales originating the orbitally-selective correlations. This physical insight results also in an appealing simplification from the computational standpoint, as it suggests that it may be generally possible to perform accurate structure prediction work assuming from the onset paramagnetic solutions (as we did in the LDA + RISB calculations of the present work), i.e., without breaking translational invariance or time reversal symmetry.

Methods

The LDA and LDA + RISB calculations were performed utilizing the DFT code WIEN2k.58 The LDA + RISB solver was implemented following ref. 32. The LAPW interface between WIEN2k and the RISB was implemented as described in ref. 59, utilizing the fully-localized limit (FFL) double-counting functional.41 All calculation were performed setting 50,000 k-points and RKmax = 8.