Density functional theory study of structural and thermodynamical stabilities of ferromagnetic MnX (X = P, As, Sb, Bi) compounds

Density functional theory (DFT) calculations for deriving enthalpies of formation H for ferromagnetic MnX (X P, As, Sb, Bi) compounds were made for the two competing structures, hexagonal and orthorhombic . Standard calculations were performed by using pseudopotentials with the generalized-gradient-approximation (PBE) as exchange-correlation functional. Enhanced exchange-correlation interactions were included by making use of a so-called DFTU approach which requires as a parameter. Application of PBE potentials for all compounds and elementary phases (all-PBE) resulted in negative values of H for MnP and MnAs in both structures whereby the result for MnP agrees very well with experiment. For MnSb and MnBi the all-PBE calculation gives a positive nonbonding H disagreeing with experiment. To overcome this discrepancy for MnSb and MnBi a DFTU ansatz was employed for all compounds and elemental Mn. The values for ranging between 0.7 for MnBi and 1.4 eV for MnAs were determined by fitting the DFT results to measured data of H. As a reference for pure Mn the -Mn phase was taken with eV by which choice the experimental volume is fitted. Atomic volumes and ionicities were derived applying Bader’s concept resulting in ionicities of Mn less than .

are metals and the term 'strongly correlated' is not appropriate. Nevertheless, performing standard DFT calculations strikingly fail for MnAs predicting the wrong structure and for MnSb and MnBi resulting in a positive formation enthalpy. Therefore, corrections need to be introduced which improve the approximations to the exchange-correlation functional. Such schemes exist such as the so-called LSDA+U or DFT+U method which assumes strongly localized orbitals (such as the 3d orbitals of Mn) replacing the local approximations based on the density by orbital-based many-body descriptions. However, at present there is no first-principles technique available for deriving the parameters U and J which occur in the DFT+U approach. Therefore, these parameters have to be chosen according to some meaningful criterion, e.g. from first principles by fitting cRPA calculations [2], or by matching experimental quantities, like lattice parameters or magnetic moments. It is therefore not surprising to find quite some spread of values among the various authors ranging from 1.2 to 3.0 eV, depending both on the X-ion and the quantity matched.
Magnetic materials with large magnetization, large aniso tropy constant, and high Curie temperature are always attractive due to their potential use in permanent magnet applications, like in electric vehicle motors and wind power generators, two important sustainable energy devices. These most important property is the maximum amount of magnetic energy these materials can store, where rare-earth element based materials have dominated for the last two decades. Yet there are two major concerns with which lead to an intensive search for new materials with reduced or no rare-earth element content: low operation temperature and limited supply with high and unstable prices. MnBi might be one of the promising candidates as it uniquely shows increasing coercivity with temperature with a value larger than that of the champion Nd 2 Fe 14 B above 150 °C, and does not contain critical rare-earth elements ( [29,30] and references therein). Furthermore, MnBi also shows a large Kerr rotation angle, essential in magneto-optical recording [26]. Another promising application of the magnetic properties of Mn-pnictides rests on the theoretical suggestion that MnAs and MnSb show half-metallic properties when crystallized in particular structures ( [21] and references therein). This property is strongly impacting spin transport and consequently the possible use in spintronic devices could be promising. Not only magnetic properties are interesting in these compounds, but also superconductivity has been discovered in MnP hinting towards a new family of superconductivity emerging at the border of long-range helimagnetic orders [10].

Density functional theory calculations
Density functional theory (DFT) was applied by using the Vienna ab initio simulation package (VASP) [31,32] with the projector augmented wave potential [33,34] construction. For approximating the exchange-correlation functional we used Perdew's generalized-gradient-approximation [35,36] (PBE). The pseudopotentials were constructed in such a way that for Mn semi-core 3s 2 and 3p 6 states were also treated as valence states resulting in 15 valence states in total. For P, As and Sb the five outermost valence states s 2 p 3 were considered whereas for Bi the 15 valence states 5d 10 6s 2 6p 3 were used. An energy cutoff for the plane wave basis of 600 eV was chosen. A k-point grid defining the smallest allowed space between k-points was constructed for all of the calculations. As an example, for the PBE equilibrium volume of MnAs in the B8 1 structure a grid of 14 × 14 × 8 k-points was obtained whereas for the B3 1 structure it was a 8 × 13 × 7 grid for the whole Brillouin zone. For the derivation of local properties such as the local magnetic moments spheres with a radius 2.2 Å for Mn, 2.33 Å for P, 2.3 Å for As, 2.98 Å for Sb and 3.09 Å for Bi were chosen. For the compounds MnSb and MnBi also relativistic calculations including spin-orbit coupling were undertaken.
For enhancement of exchange-correlation interactions a DFT+U method for the 3d wavefunctions of Mn was applied. These calculations were done choosing the simplified rotationally invariant approach of Dudarev et al [37]. This approach needs only one parameter, namely the difference U eff = U − J in which we set J = 0. For finding the appropriate U eff parameters for the MnX compounds as well as for γ-Mn we had to fit the enthalpy of formation to measured values. The accuracy of the fit for U eff was made within 0.1 eV, the values ranging from U eff = 0.7 eV for MnBi to U eff = 1.4 eV for MnAs. For MnP an average value of the largely scattered measured data was taken as reference. For γ-Mn the measured volume was chosen for fitting the value U eff = 1.3 eV.
The charge transfer, atomic volumes and charges were computed by analyzing charge densities in terms of the quantum theory of atoms in molecules according to Bader et al [38][39][40][41]. This concept utilizes the gradient of the DFT derived charge density of the solid by searching for surfaces of zero flux without falling back to any assumptions based on free atoms or ad hoc chosen atomic sphere sizes. By applying Bader's method one obtains space filling atomic volumes and the sum of atomic charges is equal to the total valence charge. The shapes of the Bader volumes are not spherical and therefore no decomposition into s-, p-or d-like atomic charges can be made.
For deriving a physically meaningful charge transfer and ionicity we calculated the atomic Bader charges for the selfconsistently derived charge density and the Bader charges for the superposed atomic charge densities in the corresponding atomic Bader volumes V Bader as listed in table 4. Their difference was then taken as charge transfer for each atom [42].

Elemental phases
For the energy of formation the total energies of the reference elemental solids of Mn and X (X = P, As, Sb, Bi) are needed. The groundstate phase of Mn, α Mn has a complicated noncollinear spin structure with a large number of atoms per unit cell. Calculations for this phase are very demanding and are done rather rarely as in [43,44]. Therefore, as many authors doing DFT calculation for Mn-compounds we focus on the high-temperature γ-Mn phase with a face centered tetragonal layerwise antiferromagnetic structure. If no relaxation of the c-axis is allowed then the crystal structure is fcc with a unit cell containing two atoms, and c/a tet = √ 2 or c/a cub = 1, depending if the lattice parameter a is taken for a truly tetragonal structure or referring to the cubic case of an fcc lattice. Once the total energy for this idealized structure is done one can take the energy difference to the α phase of 67 meV as published in [45]. Assuming that this energy difference is independent of any correlation parameter U eff the formation energy for any MnX compound can be calculated with respect to the α phase and directly compared to experiment.
The remaining question is if the calculation of the total energy of γ-Mn requires any U eff and which one. For solving this problem we refer to measurements of the volume of the γ phase of 12.94 Å 3 per atom [46]. Scanning through the volume as a function of U eff the value of U eff = 1.3 eV was determined which resulted exactly in the experimental volume. The value U eff = 1.3 eV was used throughout all the calculations of formation energies, apart from the case in which all total energies of the constituents were derived from PBE calculations. Early calculations applying PBE or the local density approximation (LDA) without any additional correlation by Qiu et al [47] resulted in much smaller volumes and values of c/a cub = 0.99 (LDA) and c/a cub = 0.94 (PBE). Figure 1 shows the dependeny of the atomic volume on the choice of U eff . There appears a low to high volume trans ition for U eff ≈ 1.4 eV which is accompanied by a low-spin to high-spin transition. The local magnetic moment of γ-Mn for U eff = 1.3 eV amounts to 2.93 µ B . Allowing for the relaxation of the c-axis the fct phase undergoes very small changes of volume/atom (12.92 Å 3 ) and c/a cub = 0.99. Also the total energy difference to the unrelaxed case is very small, namely 0.5 meV per atom.
For the calculation of the reference total energy for P we considered the monoclinic crystal structure of red Hittorf's P [48,49] with 84 atoms in the unit cell. Both volume and structural parameters were relaxed. Compared to black P (as often taken as reference for DFT calculation of compounds with P) the monoclinic phase is more stable by 30 meV/atom (2.85 kJ mol −1 atom). Other structures for white P were tried but their total energies were less negative than in the mentioned case. The remaining group V atoms, As, Sb and Bi were calculated with their ground state structure (prototype As(α), space group nr. 166) [50] with relaxation of all structural parameters.

Crystal structure of MnX compounds
Two structures were considered for the ferromagnetic MnX compounds. The unit cell of the hexagonal B8 1 structure (proto type: NiAs, space gr. no 194) [51] is spanned by the three basis vectors with atomic positions in lattice coordinates ). (2) The volume of the unit cell for the orthorhombic B3 1 structure (prototype MnP, space group nr. 62) [52] is twice as large as for B8 1 with a rectangular basal plane. In the usual crystallographic notation the z-axis for B8 1 becomes the x-axis for B3 1 . The basis vectors are now The parameters a and c have the same meaning as for B8 1 . Therefore, when we will present our calculated structural results we use the same notation by which a and c for both structures can be directly compared. When the structure is fully relaxed then Because of the lower symmetry of B3 1 compared to B8 1 the first and third coordinates of the above positions are free to relax and get distorted from the ideal positions when the B3 1 structure can be stabilized. This was the case for MnP and MnAs in the PBE calculations. In all other cases switching on U eff for MnAs, MnSb and MnBi starting with the orthorhombic structure (atomic positions taken from MnP) ended always in the high symmetry B8 1 structure after stress tensor relaxation and force minimization.

Formation energy
The formation enthalpy ∆H for all structures of the compound MnX is defined in terms of the VASP total energies E MnX , E X , E Mn by The energy E MnX per formula unit is derived from the VASP total energies for the B8 1 structure by E B81 MnX = E B81 /2 and for the B3 1 structure by E B31 MnX = E B31 /4, because the B8 1 unit cell contains two Mn and two X atoms and for B3 1 it contains four Mn and four X atoms. All structural parameters for the MnX B81 , MnX B31 compounds were relaxed. For the hexagonal B8 1 structure only volume and lattice parameters a, c needed to be optim ized whereas for the B3 1 structure also b/a and atomic coordinates had to be optimized. The total energies E(X) per atom are directly derived from the corresponding total energies for the P, As, Sb and Bi unit cells. The reference energy E γ-Mn refers to the γ phase of Mn, as discussed above, which is not the ground state. For comparison with experiment we have to add the difference of the total energies between the α and γ phases, ∆H Mn = E α-Mn − E γ-Mn , defining the formation enthalpy for each structure as  Enthalpy of formation ∆H for the ground-state structure of MnX (X = P, As, Sb, Bi) compounds in kJ mol −1 . Calculated data (marked by circles) are compared to measured data, a: [53], b: [54], c: [55], d: [56]. Dashed black line: all constituents calculated with PBE, red line: results for U eff,opt ; blue line: spin-orbit coupling included.
This makes ∆H less negative than ∆H γ by 0.067 eV or 6.46 kJ mol −1 . For MnSb and MnBi by including spin-orbit coupling relativistic calculations were performed for all constituents in (6). As illustrated by figure 2 to a very good approximation the VASP total energies scale linearly with the parameter U eff for all compounds MnX and U eff 2 for each structure of MnX (see table 1). For all calculated ∆H MnX (U eff ) the parameter U eff (Mn) = 1.3 eV is fixed and the correction ∆H Mn is independent of U eff (MnX) for the 3d-states of Mn in MnX. Then we derive for the enthalpy for each structure of MnX

MnX compounds
The results of a variety of calculations and experimental data are plotted in figure 3 showing the trend of increasing enthalpy of formation of MnX compounds with increasing size of the X atom. As it is obvious from the figure doing a standard calcul ation with PBE for all constituents (MnX, X and γ-Mn as well) gives reasonable values for MnP and MnAs but it leads to rather large errors in comparison to experiment for MnSb and MnBi. Although the trend is reproduced the thus calculated ∆H leads to positive formation energies for MnSb and MnBi. For MnSb and MnBi also relativistic calcul ations including spin-orbit coupling were performed. As expected for MnBi the relativistic effect is largest increasing ∆H by about 6 kJ mol −1 . Table 2 lists the data presented in figure 3 and also PBE results for the unstable phases MnPB8 1 and MnAsB8 1 . Whereas for MnP the PBE calculation yields the right structural stability (B3 1 more stable than B8 1 ) it gives the wrong answer for MnAs (B3 1 more stable than B8 1 ). As discussed later on doing DFT+U calculations reverses the structural stability as it is correct in comparison to experiment. Table 3 shows calculated equilibrium lattice parameters for MnX compounds with B8 1 structure and in addition also the values for the B3 1 structure of MnP, which is the ground state structure. The results are obtained by DFT+U calculation with the optimized values U eff,opt . Discussing the hexagonal structures the atomic volumes increase with increasing size of X atoms with a particularly large increase for MnAs as compared to MnP and MnSb as compared to MnAs. Interestingly, the c/a ratios 1.58 and 1.51 are similar for MnP and MnAs being significantly larger than 1.37 and 1.34 for MnSb and MnBi. Table 4  One should keep in mind that the optimization of U eff for MnP was done for the B3 1 structure and the optimized value U eff,opt = 1.2 eV was also used for calculations in the B8 1 structure. For MnAs it was done for the B8 1 structure and the resulting U eff,opt = 1.4 eV was also used for calculations of the B3 1 structure. Optimization of U eff could only be done in this way, because measurements were made for the ground state structures, namely B3 1 for MnP and B8 1 for MnAs. Figure 5 shows local moments of Mn as a function of the local Mn volume as derived from Bader's charge analysis. As already noticed before most of the theoretical work on MnP is devoted to the study of the involved magnetic properties. Early on the Mn moments and its interplay with structural properties have been in the focus of research [4,5]. More complicated magnetic interactions have been studied in [7] discussing structural conditions for metamagnetism or in [8] focussing on magnetic interactions between adjacent MnP layers of MnP-based compounds. Recent work [9, 10] reveals a spiral ground state with a pressure dependent propagation vector and a ferromagnetic state in between, as observed in the experimental phase diagram. In our calculations we find that ferromagnetic MnP clearly favors the orthorhombic B3 1 structure as readily seen from table 2. The stabilization of the B3 1 structure is largest for PBE and decreases with increasing U eff switching to B8 1 for values of U eff beyond 1.6 eV. This change in stability is likely due to the increasing magnetic Mn moment for larger U eff , which also leads to larger Table 4. Ionicity (ionic.) according to Bader's charge transfer concept in units of proton charge and Bader's atomic volume (V Bader ) in Å 3 . These data were derived from PBE calculations. The local magnetic moments (m(PBE) and m(U eff,opt ) in units of Bohr magneton µ B were obtained from PBE and DFT+U(U eff,opt ) calculations with the optimized parameter U eff,opt (in eV).   atomic volumes, favored by the B8 1 structure. Comparing to other theoretical work, there is surprisingly little information available for formation energies, only [6] give a value of −85.6 kJ mol −1 for the B3 1 structure calculated with PBE for exper imental lattice parameters which is much less negative than our PBE value, mostly due to their larger cell volume of 12.27 Å 3 /atom. The values given by the Materials Project for B3 1 -MnP [57], 119 kJ mol −1 , and B8 1 -MnP [58], 89 kJ mol −1 , agree well with our PBE values for ∆H and slightly worse with the experimental values in table 2.

MnAs.
MnAs was studied mostly for magnetic ordering, e.g. by mapping on a Heisenberg Hamiltonian in order to calculate the Curie temperature [14]. An extensive and thorough study by an ab initio approach for the effective onsite Coulomb interaction parameters [2] resulted in values for MnAs of U LDA+U = 1.73 eV and J = 0.56 eV, yielding (U − J) = 1.17 eV agreeing rather well with our calculated value of 1.4 eV. It should, however, be noted that the calculations in [2] were done for the cubic zinc-blende structure with an atomic volume of 22.54 Å 3 per atom whereas in our case the volume per atom is 17.48 Å 3 for the geometrically relaxed B8 1 structure.
Because of the expected competition of the B8 1 and B3 1 structures also for MnAs both structures were studied. For this we started with the converged structural parameters as calculated for MnP and let VASP relax the structure in the MnAs calculation. For PBE both structures could be stabilized with the unexpected results, that the total energy for B3 1 was lower than for B8 1 contradicting the experimental findings of a stable hexagonal B8 1 structure below 300 K.
As shown in figure 6 for the PBE calculation the total energy for the B8 1 structure is less negative than for B3 1 with a local energy minimum around B8 1 . The structural stability changes when the VASP calculation is done with a small U eff = 0.5 eV. Then, the energy is lowest at B8 1 showing a shallow minimum at the B3 1 structure. Already at U eff = 0.2 eV the B3 1 structure changes into B8 1 when a stress tensor technique is used by VASP for structural relaxation. The closeness of both stuctures is manifested by the thermodynamic study in [54] which shows a remarkable λ-like phase transition at 315.6 K above at which temperature the B3 1 structure is stabilized. On the theoretical side, the close competition between both structures is also found in [12] where for PBE the cohesive energy was found to be the same for both structures. However no values for the formation energy were given.

MnSb.
Magnetism is also in the center of theoretical research on MnSb. Early calculations at experimental geometries [20] found a Mn magnetic moment of 3.3 µ B and a small negative moment on the Sb site in good agreement with our results and other work [4,11]. Interestingly, the stabilization of a ferromagnetic alignment was explained by covalent interactions of unoccupied states of the Sb-5p band with the Mn-3d bands [20]. DFT+U calculations [23] with U eff = 3 eV gave values for the structural parameters similar to the present ones, e.g. the atomic volume, V at = 20.96 Å 3 , is only slightly below our calculated value of 21.07 Å 3 . However, from the U eff dependence of the volume presented in [23] one would expect a reduction of V at for smaller values of U eff worsening the agreement. The results of DFT calcul ations [24] were also used in high-temperature series expansions of the magnetic susceptibility to predict quite well the Curie temperature of ferromagnetic MnSb. Calculations for cubic MnSb [22] predicted half-metallicity which makes it a promising candidate for room temperature spin injection into semiconductors. Although the cubic polymorph is predicted to be less stable by ≈80-85 kJ mol −1 and predicted to decompose [21] there is evidence for a successful epitaxial growth [22]. To our knowledge no theoretical estimates of formation energies exist besides a much too small value of −6.56 kJ mol −1 as reported by the materials project [59] yielding a non-magnetic configuration.

MnBi.
Regarding theoretical studies MnBi poses no exception from the focus on magnetic rather than thermodynamical properties. The first spin-polarised band-structure calculations [25] produced a magnetic moment of 3.6 µ B in agreement with our results and [11] and showed the Mn-3d electrons are partially delocalised and strongly hybridized with Bi-6p states. MnBi also shows a giant Kerr effect which according to the calculations in [26] is caused by the combination of the large Mn magnetic moment and the large spin-orbit coupling of Bi together with the strong hybridization between the Mn-3d bands and Bi-6p states. A reasonable description of the magnetic properties like saturation magnetization, anisotropy constant, and Curie temperature was found both in [27] and [28]. Recent DFT+U calculations with U eff = 2 eV [29] successfully explained temperature dependent anomalies of lattice constants and magnetic properties and were able to pin-point the sign change of the magnetic anisotropy energy to small increases in the basal lattice constant, but also stated that a single U eff cannot be adjusted to reproduce both lattice constants and magnetic moment perfectly. Consequently the calculated value of V at in [29] amounts to 25.23 Å 3 much larger as compared to our present value of 23.57 Å 3 for an Figure 6. VASP total energy difference E trans − E B81 for MnAs along the structural transition path going from B8 1 to B3 1 . The intermediate structures were obtained by by linearly interpolating between the PBE derived structural parameters of the B8 1 and the B3 1 structures. The same structures were used for the U eff = 0.5 eV calculation.
U eff of 0.7 eV, which, however, is only slightly larger than the PBE value of 23.06 Å 3 reported also in [29]. Yet, concerning formation energies to our knowledge only a positive DFT value of 88.5 kJ mol −1 is reported by the materials project [60] yielding a non-magnetic configuration which readily decomposes into elemental Mn and Bi.

Summary and conclusion
We tried to predict the energies of formation and structural stabilities for ferromagnetic MnX compounds. However, DFT calculations with a parameter free local exchangecorrelation functional such as the PBE aproximation, failed. The PBE results (i.e. compounds and elemental solids were all calculated with PBE potentials) in figure 3 show good agreement with experiment for MnP and MnAs, whereas for MnSb and MnBi ∆H is zero or even significantly positive. Compared to experiment and the calculations with U eff , ∆ H increases too fast when going from MnP to MnBi (i.e. larger sizes of the X atoms). When shifting the all-PBE data down by a constant by about 70 kJ mol −1 in order to get a good agreement with experiment for MnBi then the formation energy for MnP is far too negative. Only a nonconstant shift varying with X could possibly give the right correction. Assuming that the PBE calcul ations for the MnX compounds and elemental Mn are correct, then the correction could only come from the elemental phases of the X atoms. For useful corrections one would need less negative total energies for Sb and Bi and the change must be quite substantial, namely 60-70 kJ mol −1 , which we consider as an unrealistic error. Therefore, we conclude that the largest error of ∆H comes from the PBE total energy calculations for the MnX compounds.
Discussing now the DFT+U calculations one has to keep in mind, that two different U eff for Mn 3d-states were used, namely one for γ-Mn (U eff,opt = 1.3 eV) and X-dependent values of U eff,opt for fitting to the measured formation energies of the MnX compounds. Interestingly, U eff,opt = 1.2 and 1.4 eV for MnP and MnAs agree rather well with the value 1.3 eV for γ-Mn, and for these compounds the PBE results are good. One might argue, that because of the similar U eff,opt the enhanced exchange-correlation effects basically cancel each other in the energy difference by which ∆H is defined in (6) and (7). The situation is quite different for MnSb and MnBi with their smaller values of U eff,opt = 0.9 and 0.7 eV. Would we increase these values to about 1.3 eV then the total energy becomes less negative and as consequence the formation energy is wrong, as discussed for the PBE results. The values U opt,eff for MnSb and MnBi are significantly smaller than 1.3 eV for γ-Mn, and the enhanced exchange-correlation effects do not cancel. In comparison to PBE for MnSb and MnBi the total energies do not increase (i.e. getting less negative) so much as for Mn when switching on U eff . In other words, the repulsive effects introduced by DFT+U are less for MnSb and MnBi than for γ-Mn. Therefore, ∆H is lowered in comparison to the PBE results and some kind of 'bonding' effect is mimicked. Summarizing our study, it is obvious that even for not strongly correlated systems, such as the MnX compounds with U eff,opt 1.4 eV corrections beyond the standard PBE approach are needed in order to get a correct enthalpy of formation. It would be higly desirable to design such corrections free from empirical parameters, or at least being valid for some larger class of materials.