Crystal Field Splitting is Limiting the Stability and Strength of Ultra-incompressible Orthorhombic Transition Metal Tetraborides

The lattice stability and mechanical strengths of the supposedly superhard transition metal tetraborides (TmB4, Tm = Cr, Mn and Fe) evoked recently much attention from the scientific community due to the potential applications of these materials, as well as because of general scientific interests. In the present study, we show that the surprising stabilization of these compounds from a high symmetry to a low symmetry structure is accomplished by an in-plane rotation of the boron network, which maximizes the in-plane hybridization by crystal field splitting between d orbitals of Tm and p orbitals of B. Studies of mechanical and electronic properties of TmB4 suggest that these tetraborides cannot be intrinsically superhard. The mechanical instability is facilitated by a unique in-plane or out-of-plane weakening of the three-dimensional covalent bond network of boron along different shear deformation paths. These results shed a novel view on the origin of the stability and strength of orthorhombic TmB4, highlighting the importance of combinational analysis of a variety of parameters related to plastic deformation of the crystalline materials when attempting to design new ultra-incompressible, and potentially strong and hard solids.


Stabilization from oI10[71]- to oP10[58]-TmB
. We first searched the most stable structures of stoichiometric TmB 4 by geometry optimizations of 26 commonly observed Tm-B, Tm-C, Tm-N, Tm-Al, Tm-Si ICSD structure types 34 . The geometries were optimized using DFT as implemented in the VASP code 35,36 . The details of the adopted initial structures of TmB 4 are given in the supplementary materials. The formation energy/enthalpy 1 5 4 was calculated from the chemical reaction Tm + 4B = TmB 4 . The energies of α − B and Tm in magnetic state were adopted as the reference ground states. Figure 1 shows the calculated formation enthalpy of CrB 4 and FeB 4 . The most negative formation enthalpies of oP10 [58] indicate that it is the most energetically favored structure. As listed in Table S1 of the supplementary materials, eight typical structures, such as oP10 [58], oI10 [71], oP10 [59], hP20 [194], hP10 [194], mS10 [12], mS20 [14] and mS30 [12], have been reported for transition metal tetraborides. In the present study, the oP10 [58] is found to be the most stable one for both CrB 4 and FeB 4 . The formation enthalpies obtained in the present study for the oI10[71] structure of CrB 4 (ΔH = − 0.301 eV/atom) and of FeB 4 (ΔH = − 0.143 eV/atom) are slightly higher than those calculated in previous publications of ΔH = − 0.3098 eV/atom for CrB 4 37 , and of ΔH = − 0.1698 eV/atom for FeB 4 37 . Note that the initial structure with the mS20 [14] symmetry transforms into oP10 [58] for both CrB 4 and FeB 4 . Our evolutionary structure search scheme 38,39 confirms further that oP10 [58] is the most stable one. A systematic thermodynamic investigation of possible ground state structures with consideration of vibrational entropy has been performed on FeB 4 by Kolmogorov et al. 33 Their results indicate that the oP10 [58]-FeB 4 are thermodynamically stable in the considered temperature range, but it lies 3 meV/atom above the α− B <−> oP12-FeB 2 tieline at T = 0 K while 10 meV/atom below the tieline at T = 900 K 33 . The decomposition energy from FeB 4 to oP12-FeB 2 and α-B are calculated to be positive (~0.01 eV/atom), while the energy from FeB 4 to oP8-FeB and α− B are negative (~ − 0.02 eV/atom), in agreement with previous calculations 33 .
A similar trend of thermodynamic stability is also observed for MnB 4 in oI10 [71] and oP10 [58] structures. The formation enthalpy of oI10[71]-MnB 4 of − 0.2704 eV/atom decreases slightly by 0.01 eV/atom (0.0097 eV/atom 37 ) due to the transformation to the oP10 [58] structure whose formation enthalpy is − 0.2804 eV/atom. All these calculations suggest that oI10[71]-TmB 4 is thermodynamically unstable with respect to oP10 [58]-TmB 4 . Note that a Peierls distortion mechanism has been proposed as the origin of stabilization of mS20 [14]-MnB 4 . However, such distortion will be hindered by temperature and the phonon assisted crossover observed from nonmagnetic Peierls insulator (mS20) to a magnetic Stoner metal (oP10) 30 . Taking CrB 4 as a prototype, Fig. 2a [58] transformation, the molar volume changes a bit and the simulated XRD figures show a similarity between oI10[71] and oP10 [58], but some minor XRD peaks show up for oP10 [58] structure which may be identified by experiments 40 .
As the next step we calculated the phonon dispersion curves as shown in Fig. 2c 33 It should be noted that the oI10[71]-TmB 4 structure comprised of tetragonal boron network has been previously examined within the extended Huckel method. It has been concluded that maximum binding in the 3d series is achieved for Cr, and that the electron-rich Mn, Fe, Co and Ni tetraborides may be unstable in this configuration 41  With focus on the phonon-assisted transformation from oI10[71] to oP10 [58], we attempted to find the electronic origin of the stabilization from higher symmetry structures to those with lower symmetry. After comparing the electronic structure of oI10[71]-and oP10 [58]-CrB 4 , we found that the transformation from oI10[71] and oP10 [58] induces a highly directional electronic partition at two inequivalent sites of the boron atoms in the oP10 [58] structure: one with minor charge transfer of ~0.12 electrons, and the other one with a higher charge transfers of ~0.32 electrons. Eight equivalent boron sites are found in oI10[71] structure with an average Bader charge of ~0.22 electrons. The stronger anisotropy of charge transfer in oP10 [58] structure indicates its higher electronic polarization and localization as compared to oI10[71] one. The Bader charge density analysis 42 shown in Fig. 2 confirms the different charge transfer at crystallographic sites of Cr for oI10[71] structure (+ 0.90 electrons) and oP10 [58] structure (+ 0.88 electrons), suggesting a slightly higher ionic contribution for the former. The anisotropy of Bader charges can be attributed to the changes from regular to distorted boron network where metal atoms are donating electrons to stabilize the oP10 [58] structure. The analysis of crystal orbital overlap population (COOP) curves of oI10[71]-and oP10 [58]-CrB 4 shown in Fig. 2e,f, provides a further confirmation of slightly higher COOP(Cr-B) of 0.28 in oP10 [58] structure as compared with oI10[71] structure (0.26). Note that an average COOP value of bonds is used to do the comparison for convenience. Similar higher values of COOP(Fe-B) (0.23) and COOP(Mn-B) (0.20) are also shown for oP10 [58] structure as compared to oI10[71] structure (0.21 and 0.18 for Fe-B and Mn-B bonds respectively) (refer to Figs S3 and S4 in supplementary materials). This suggests a slightly stronger bond strength for oP10 [58] structure. When observing the structures in Fig. 2a,b one notices much "dense" bond network in oP10 [58] structure as compared with oI10[71], which is a qualitative but illustrative confirmation of the analysis of the electronic structure. The calculated valence charge density difference (VCDD) of oI10[71]-and oP10 [58]-CrB 4 shown in Fig. 2 provide additional support of their stabilization by crystal field splitting of the d orbitals. By comparing the isosurface of VCDD of the oI10[71]-CrB 4 ( Fig. 2a) with that of oP10 [58]-CrB 4 (Fig. 2b) we see a big difference in the shape of VCDD isosurface at the metal sites: the VCDD isosurface show anisotropic shape in oP10 [58]-CrB 4 , whereas they are relatively isotropic in oI10[71]-CrB 4 . This suggests that although oP10 [58] has a lower lattice symmetry, it possess a stronger electronic directionality (see Fig. 2a,b), in agreement with the analysis of charge transfer and COOP analysis in the preceding paragraph. In order to get such electronic partition of d orbitals, the boron network in oI10[71] structure needs to reorganize its positions to transform into oP10 [58] structure by rotation operations. Furthermore, the boron-boron bonds indicated by the black arrows in oP10 [58] 4 , the DOS around E F (0.88 states/eV cell) is slightly higher than that in oP10 [58]-CrB 4 (0.80 states/eV cell), indicating a stronger "splitting" of bonding and antibonding states for oP10 [58]-CrB 4 , in agreement with the preceding COOP analysis. A similar feature is observed for MnB 4 , but does not apply for oI10 [71] and oP10 [58]-FeB 4 in which the Fermi level is located at the left shoulder of a peak in the DOS. Such differences will be shown to be responsible for their different shear moduli in the following section. Interestingly, we find in Fig. 3 a distinct change of the majority of − d y x 2 2 and d xy orbitals (which mainly contribute to the in-plane rotation) below and above E F when the structure changes from oI10[71] to oP10 [58] for both CrB 4 and FeB 4 , suggesting that crystal field splitting of d orbitals is responsible for the phonon-assisted stabilization from oI10[71] to oP10 [58] as discussed above. In case of CrB 4 , the first peak of d x 2 −y 2 below the E F move upward by 0.1 eV during the transformation from oP10 [58] to oI10[71], while the first peak of d xy below E F (at − 1.167 eV) nearly disappears. For FeB 4 however, the unexpected peaks show up at the Fermi level for oP10 [58] structure, but not for oI10[71] structure. 4 . We next explore the mechanical properties of the oP10 [58]-TmB 4 with Tm changing from Cr to Mn and to Fe. The calculated elastic constants of oP10 [58]-TmB 4 are listed in Table S2 of the supplemental materials and compared with the results from other publications. It is shown that our present calculations are in good agreement with the others, and all the three orthorhombic tetraborides meet the requirements of elastic stability 43 The Voigt average bulk and shear moduli are derived from elastic constants and compared with other hard materials in Table 1 19,23,27,28,[31][32][44][45][46][47][48] . It is seen that all three oP10 [58] tetraborides are ultra-incompressible materials owing to their high bulk moduli of more than 250 GPa. A relatively high value of shear moduli is observed for oP10 [58]-CrB 4 (267 GPa) and oP10 [58]-MnB 4 (247 GPa), suggesting they are stiff. However, a much lower value found for oP10 [58]-FeB 4 (192 GPa), indicates that it may be less stiff than the other two tetraborides, in agreement with the aforementioned DOS analysis that suggests a stronger metallic bonding for FeB 4 .

Mechanical properties of oP10[58]-TmB
Well established criterion for whether a solid is ductile or brittle is if dislocation embryos can be nucleated from an atomically sharp crack prior to the crack propagation 49,50 . If this is the case at a given temperature, the solid is considered to be intrinsically ductile 49,50 . This distinction of the mechanical behavior is characterized by the ratio of the shear to bulk modulus G/B, based on the consideration that B is a measure of the resistance to fracture and G the resistance to plastic deformation 49,51,52 . The critical G/B ratio which separates ductile and brittle materials is around 0.57, i.e. if G/B < 0.57 the material is ductile, and it is brittle when G/B > 0.57 51,52 . In view of  [58] corresponds to a change of DOS shape and relative position of the decomposed d orbitals, which is more pronounced for FeB 4  the relatively high ratio of G/B ≈ 0.957 for oP10 [58]-CrB 4 , 0.879 for oP10 [58]-MnB 4 , 0.674 for oP10 [58]-FeB 4 , we conclude that all these tetraborides are brittle 49,51 . Even in brittle solids, a brittle to ductile transition is possible at an elevated temperature where the nucleation of dislocation embryos becomes possible prior to cleavage propagation 50 .
It is known that high elastic moduli do not guarantee high resistance against plastic deformation because upon finite shear electronic instabilities and transformations to softer phases with lower shear resistance often occur 9,26 . The anisotropic ideal shear strength of oP10 [58]-TmB 4 may be obtained from the calculated stress-strain relationships along different crystallographic planes and directions. The slip systems of (110) [1][2][3][4][5][6][7][8][9][10] and (010)[101] were reported as the weakest ones for oP10 [58]-CrB 4 by Li et al. 23 and by Niu et al. 19 , respectively. A similar disagreement exists also for the weakest slip system in FeB 4 , where (110)[001] has been reported by Li et al. 28 but (111)  by Zhang et al. 48 . We studied all the four slip systems in more detail. Our results suggest that the (110) [1][2][3][4][5][6][7][8][9][10] slip systems is the weakest one for CrB 4 , whereas (110)[001] is the weakest for FeB 4 and MnB 4 . Figure  4a-c show the calculated dependence of the stress, magnetic moment and total energy on strain for oP10 [58]-CrB 4 , oP10 [58]-MnB 4 and oP10 [58]-FeB 4 along their weakest slip systems. The minimum strength for the three TmB 4 is shown in Table 1 and compared with previous results for these tetraborides 19,23,27,28,[30][31][32] , and with those of WB 3 45,47 . It is seen that the minimum shear strengths of oP10 [58]-CrB 4 of 36.2 GPa in the (110) [1][2][3][4][5][6][7][8][9][10] slip system are comparable to those of WB 3 (37.7 GPa) and ReB 2 (34.4 GPa), but they are much lower than those of c-BN (58.3 GPa 46 ). This suggests that CrB 4 cannot be intrinsically stronger than WB 3 in view of the fact that its strength is between that of WB 3 and ReB 2. These results are in agreement with the recent experiments in which the load-invariant ("asymptotic") Vickers hardness of CrB 4 (23.3 GPa 21 and 28.6 GPa 22 ) was found to be comparable to that of ReB 2 . The high hardness of 48 GPa for CrB 4 "predicted" by Niu et al. 19 on the basis of the "hardness model" is not justified. Similar conclusions apply also for MnB 4 and FeB 4 because a much lower ideal shear strength is found for oP10 [58]-MnB 4 (35.1 GPa) and oP10 [58]-FeB 4 (25.7 GPa), although Niu et al. 19 predicted a high value of hardness of 41.5 GPa for MnB 4 . The reason is simply the fact that the "theoretical models of hardness" compares in fact elastic stiffness, and do not account for electronic instabilities upon a finite shear 18,26 . Bond deformation mechanism of oP10[58]-TmB 4 . In view of the mechanical weakness of oP10 [58]-TmB 4 , we investigated the bond deformation paths and electronic instability mode of oP10 [58]-TmB 4 along the weakest shear path. Figure 4a,c show that both stress-strain and total energy-strain curves are smooth, i.e. no distinct electronic instabilities occur as in the cases of ReB 2 18 and WB 3 6 , but not for the change of magnetic moment for FeB 4 . Figure 4d-f show the variation of the VCDD from the equilibrium (γ = 0.0000), to a shear strain γ = 0.1717 corresponding to the stress maximum and at shear strain γ = 0.4002 after the shear instability for CrB 4 in the weakest (110)[1-10] slip system. One can see the change of in-plane lengths of the B1-B2 and B3-B4 bonds with respect to (001) planes from ~1.85 Å at equilibrium to ~2.11 Å at maximum stress. With the elongation of B1-B2 and B3-B4 bond lengths, the magnitude of charge depletion between them increases gradually before the stress reaches the peak value of about 25.7 GPa (see Fig. 4e). After the peak stress, the B2 atom and B4 atom come closer together and form B2-B4 bonds with length of ~1.97 Å at a strain of 0.4002, while the B1 and B3 atoms further separate. The shear deformation after the peak stress is accompanied by a stress decrease to about 5 GPa (Fig. 4a).
In FeB 4 and MnB 4 , the weakest link does not appear along the (110) [1][2][3][4][5][6][7][8][9][10] shear. Therefore we concentrate on the deformation within the (110)[001] slip system. In view of the similarity between FeB 4 and MnB 4 , we shall take FeB 4 as an illustrative example. Figure 4g-i show the variation of the VCDD isosurfaces in equilibrium (γ = 0.0000), at a shear strain of γ = 0.2434 corresponding to the maximum stress, and at γ = 0.4002 after the shear instability of FeB 4 , respectively. No significant bond weakening or breaking appear before the shear stress reaches the maximum value at strain of about 0.25. After the shear instability, a multiple out-of-plane bond weakening and breaking appears between the {001} planes (see the regions marked by black arrows in Fig. 4h,i). In the following orbital analysis we shall show that such bond weakening is attributed to the electronic instability between the d orbitals of Fe and p orbitals of B, which are responsible for the out-of-plane bonding between {001} planes.

Electronic origin of lattice instability under shearing.
For that purpose, we analyze the change of orbital-decomposed DOS of d orbitals of Tm and p orbitals of boron in TmB 4 before and after lattice instability. Figure 5a-f shows the orbital-decomposed DOS of CrB 4 at equilibrium, at a strain of 0.1717, and at a strain of 0.4002. Upon shear but before the lattice instability, the occupied d xy states located below E F move towards E F (see the red curve between − 2 eV and 0 eV in Fig. 5a,b), while the peak of d xy above E F move downwards and split into two parts between 0 and 1 eV, and between 1 and 2 eV. The splitting feature is also observed for the d z 2 states above E F (see blue curve between 0 eV and 2 eV in Fig. 5a,b). During this process, the shape of the other three d-states around E F does not change significantly. After the lattice instability (see Fig. 5c), most of the d xy and d z 2 states between − 2 eV and 0 eV have vanished, and both contribute substantially to peaks in the DOS between 0 and 2 eV above E F . At the same time, the pseudogap in the p z states of boron became narrower and finally vanished, while p x and p y increase their contribution at E F after the lattice instability. At the peak stress, all three p orbitals contribute to the states at E F . Therefore, the change of the DOS arising from the d xy orbitals of Cr and from p x and p y orbitals of B are mainly responsible for the in-plane bond weakening of B1-B2 and B3-B4 bonds within the x-y plane, as shown above. Figure 6 shows the orbital-decomposed EDOS of the (110)[001] slip system of FeB 4 in equilibrium, at a strain of 0.2434 corresponding to the stress maximum, and at strain 0.4002 after the shear instability. In equilibrium, dz 2 states show a pseudogap at E F whereas the other four d orbitals contribute mostly to the finite value at E F . Before the instability, the peaks of all five orbitals between − 6 eV and 0 eV move towards E F , and at the maximum stress one profound peak of d z 2 character is seen at the Fermi level giving large contribution to the states at E F , while nearly all peaks of other spin-up orbitals above the Fermi level vanish. After the lattice instability, the states of all spin-down orbitals move further upwards, whereas the d xz spin-down states appear above Fermi level. The appearance of boron p z states at the Fermi level and the decrease of peaks of p x and p y states can be clearly seen in Fig. 6d-f. Therefore we can conclude that the electronic instability of FeB 4 during (110)[001] shear is mainly due to the changes of the energies of states of Fe d z 2 and d xz orbitals, and of p z orbitals of B. MnB 4 shows a similar instability as FeB 4 , therefore we do not discuss it here in any further detail (see Fig. S8 in supplemental materials).
The electronic instability of TmB 4 (Tm = Cr, Mn and Fe) upon shear due to d orbitals of Tm and p orbitals of B does not resemble the cases of ReB 2 18 and WB 3 6 . The important difference is that in TmB 4 the covalent bond networks are mainly responsible for the shear instability, while metal-boron or metal-metal bonds are the carrier of the shear instability in ReB 2 and WB 2 . Although the three dimensional covalent network has been proposed to be responsible for high hardness in transition metal borides for a long time, our results show uniquely that the 3D boron network in the tetraborides may not be strong enough to provide these materials with high plastic resistance. Because the strength of a real crystalline material is limited by a variety of defects, such as dislocations, cracks, grain boundaries and others, it is usually orders of magnitude smaller than that of an ideal crystal. The Peierls-Nabarro (PN) stress of dislocations provides a realistic explanation on plastic resistance of a real crystal and can be correlated to shear moduli and ideal shear strength [53][54][55] where α = 2π and K = G/(1− ν ) for edge dislocation and K = G for screw dislocation, b is the Burgers vector, ν is the Poisson ratio, and ζ is the dislocation half-width which can be expressed as ζ πτ = Kb/4 max , where τ max is the ideal shear strength. Thus, the plastic resistance of a crystal depends mainly on the shear modulus, ideal shear strengths, lattice topology, bonding nature and deformation mode. The present work emphasizes the necessity to conduct a combinational analysis of the lattice stability, of the values of shear moduli, of the anisotropic shear strength, of the complexity of lattice topology and its changes during shear, of the nature of chemical bonding, and others which are critical parameters in Peierls-Nabarro dislocation model of plastic lattice resistance [53][54][55] . Despite a century of research on dislocation mediated plasticity, a systematic connection to those quantities derived from first principle calculations has not been rationalized, inducing massive applications of a single parameter e.g. elastic moduli or ideal shear strength, to directly quantify the mechanical strength and hardness of a real material 3,4,[9][10][11][12]19,23,24,28 .

Discussion
In summary, we carried out comprehensive density functional theory calculations to determine the thermodynamic, mechanical and phonon stabilities, the deformation paths and the electronic instability modes of orthorhombic TmB 4 . The electronic structure calculations reveal that the transformation from oI10[71]-TmB 4 to oP10 [58]-TmB 4 can be explained by the variation of the in-plane electronic structure within (001) planes and the formation of new boron-boron bonds at hollow sites by distortion of boron networks. These processes significantly enhance the electronic hybridization of d orbitals of Tm with p orbitals of B by crystal field splitting. Depending on the different deformation paths (slip systems), different orbitals are responsible for the electronic instability upon finite shear. The relatively low shear moduli and ideal strengths of TmB 4 suggest that these materials cannot be intrinsically superhard. The instability of weak 3D boron covalent networks is found to be responsible for the weakness of oP10 [58]-TmB 4 . These results illustrate the importance and necessity of a combinational analysis of a variety of parameters related to plastic deformation of the crystalline materials, and their combination in the attempt to design new intrinsically hard and superhard materials.

Methods
First principles calculations. The present first-principles density functional theory (DFT) calculations of the formation energy of TmB 4 , Tm = Cr, Mn and Fe were done by means of the VASP code 35 with the projector augmented wave method 36 employed to describe the electron-ion interaction and the generalized-gradient approximation of Perdew-Burke-Ernzerhof (PBE) 56 for the exchange correlation term. The integration in the Brillouin zone has been done on special k grids for the phases that were under consideration, which were determined according to the Monkhorst-Pack scheme, the energy cutoff of 600 eV, and the tetrahedron method with Blöchl corrections for the density of states and smearing methods for the stress calculations. The verification of the reliability of our calculations has been done by calculating the total energies, equilibrium lattice parameter, and bulk modulus of TmB 4 and compared with previous theoretical and experimental values in the present work, and of a number of other materials in our earlier work. In the present studies, all DFT calculations are performed with spin-polarized scheme because of the appearance of magnetic momentum during the deformation. Nevertheless, in cases where no magnetic momentum has been found, additional non-spin-polarized calculation is used to check its possible influence.
The phonon dispersion. The phonon calculations were performed within the harmonic approximation using the direct method 57 based on the calculated non-vanishing Hellman-Feynman forces employing the Phonopy code 58 . To confirm our results, we also used a linear response method based on the perturbation theory as implemented in VASP code.
2 , (r8)ε = (δ,0,δ,0,0,0) with relation τ τ δ δ (r9)ε = (0,δ,δ,0,0,0) with relation τ τ δ δ The stress-strain relationships. First, the atomic basis vectors of a given crystal cell were projected onto the Cartesian coordinates R with one axis vector being parallel to the imposed strain direction for the tension. For the shear, one axis vector was perpendicular to the slip plane and another one was parallel to the slip direction in that plane. Afterwards, the crystal has been incrementally deformed by transforming the unstrained atomic basis vector matrices R to the strained ones R′ using the deformation matrices. where e 1 = e xx , e 2 = e yy , e 3 = e zz , e 4 = e zy + e yz , e 5 = e zx + e xz , and e 6 = e yx + e xy are the strain components in the Voigt notation. The diagonal strain components e xx , e yy , and e zz represent the tension while the off-diagonal components represent the shear. In order to keep the crystal under a stress state of uniaxial tension or shear, the strained cell has been relaxed for both the atomic basis vectors and for the atom coordinates inside the unit cell by keeping the applied strain component fixed and relaxing the other five strain components until their conjugate stress components i.e., Hellmann-Feynman stresses reached negligible values. Such a relaxation scheme is accomplished by using a slightly modified VASP code with specific constraints of strain components. To ensure that the strain path is continuous, the starting position at each strain step has been taken from the relaxed coordinates of the previous strain step. In the instance of having a large strain, the crystal symmetry may be changed and the Brillouin zone significantly deformed. Therefore we adopt a high energy cutoff of 600 eV and verified the convergence of the calculations of the stress-strain curves by using different meshes of k points. Although the spin-polarized calculation does not have big impact on the shear strength on FeB 4 and CrB 4 , we performed the spin-polarized calculations for a comparison because of the significant magnetic effect on the results of MnB 4 and possible appearance of magnetic moment during severe deformation. In each run the same low enthalpy structure was found at least 3 times (often more) before the run was terminated. A six-step structural-optimization scheme was used for all runs, and each step employed the geometry of the structure from the previous step for the initial geometry. Only the ions were allowed to relax in the first two steps, while the last four steps also allowed the lattice parameters to relax. The precision of the calculation was increased at each step. The calculations were performed using the PBE exchange-correlation functional, the PAW method and an energy cut-off of 450 eV in the final step. The lowest energy structures from each search were optimized using the aforementioned computational settings, to obtain a more accurate energy ranking.