Ab initio thermodynamic properties of certain compounds in Nd-Fe-B system

In this work, we report the results of ab initio calculations of thermochemical properties of several compounds in the Fe-Nd, B-Nd and B-Fe-Nd systems. We have performed DFT+U calculations to compute the enthalpy of formation of the compounds NdB6, NdB4, Nd2B5, Nd2Fe17 and Nd5Fe2B6. It was found that the values obtained with an effective Hubbard U correction have better agreement with the experimental data. We have also computed the vibrational contribution to the heat capacity (Cp) of the compounds as a function of temperature was computed using the quasharmonic approximation. For most of the compounds these properties have not been experimentally determined until now. Hence, the computed ab initio thermodynamic properties will serve as useful input for the Gibbs energy model parameter assessment using the Calphad method.


I. INTRODUCTION
Ab initio calculations based on density functional theory (DFT) has gained much popularity mainly because of the very reliable results that can be obtained through these computations. Amongst many properties that can be computed by this method, the thermochemical properties such as enthalpy of formation, enthalpy of mixing, heat capacity, etc. are of particular interest. They serve as important input in the Gibbs energy modelling in the context of the CALPHAD (acronym for CALculation of PHAse Diagrams) approach, especially when the corresponding experimental data are not available. This greatly improves the quality of the generated Gibbs energy functions. With the advancements in high-performance computing, it is now possible to carry out these calculations for systems with a large number of atoms. In this work, we report the results of ab initio calculations of thermochemical properties of the compounds NdB 6 , NdB 4 , Nd 2 B 5 , Nd 2 Fe 17 and Nd 5 Fe 2 B 6 .
The rare-earth and transition metal compounds with localized 3d/4 f orbitals constitute what is known as the "strongly correlated" systems. Compounds of this class contain localized electronic orbitals that introduce significant same-site Coulomb repulsion [1]. Mott insulators are also a well-known group of materials belonging to this category.
A shortcoming of DFT is the failure of the conventional exchange-correlation functionals (LDA/GGA) in describing the quantum mechanical nature of the strongly correlated systems. For instance, the DFT calculations using LDA/GGA functionals identifies the Mott insulators as metals. On the other hand, Mott insulators are well described by the Hubbard or Anderson models based on the tight-binding approach [1]. The U parameter in the Hubbard model takes care of the effect of same-site repulsion. This parameter is adopted within DFT and gives rise to the so-called DFT+U method, which introduces a correction to the total energy [2]. This approach has thus increased the accuracy of the results of DFT calculations concerning strongly correlated materials. These orbitals are specially treated in the DFT+U scheme with the correction using a Hubbard U parameter.
To the best of our knowledge, ab initio calculations to compute the finite temperature thermodynamic properties of these compounds have not been reported till now. Li et al. [3] has reported ab initio thermodynamic properties of several borides relevant to the Fe-B system. Recently, another ab initio investigation [4] reported enthalpies of formation of Nd-B compounds using GGA functionals.
The thermodynamic modelling of Nd-Fe-B system by van Ende et al. [5] made use of the available experimental C p data for Gibbs energy modelling of NdB 6 and NdB 4 , whereas C p of Nd 2 B 5 , Nd 2 Fe 17 and Nd 5 Fe 2 B 6 were approximated using Neumann-Kopp rule due to lack of experimental data. It is envisaged that the results of the present work could be used in order to further improve the thermodynamic modelling of Nd-Fe-B system.

II. MODEL AND METHODOLOGY
In this work, DFT calculations were carried using the QUANTUM ESPRESSO [6] code, which uses a plane-wave arXiv:1811.06744v2 [cond-mat.mtrl-sci] 22 Mar 2019 basis set to expand the wavefunctions. Core electrons were represented with projector augmented waves (PAW) [7] pseudopotentials. These pseudopotentials use 3 valence electrons (2s 2 2p 1 ) for B, 8 valence electrons (4s 2 3d 6 4p 0 ) for Fe, and 14 valence electrons (5s 2 5p 6 6s 2 4f 4 5d 0 ) for Nd. The GGA exchange-correlation functional parametrized by Perdew-Burke-Ernzerhof [8] (PBE) was used. Marzari-Vanderbilt smearing [9] was used to smooth the discontinuities in integrals due to the Fermi surface. To aid convergence in calculations, the largest smearing width σ was used, which does not introduce ground state energy deviation larger than 0.002 Ry/atom. The width of Marzari-Vanderbilt smearing used in this study ranged from 0.01 Ry up to 0.05 Ry.
The calculations were done on k-points selected using the Monkhorst-Pack scheme [10]. Convergence tests were carried out for every structure up to an energy convergence of 0.001 Ry/atom. The same convergence test was applied to determine the appropriate cut-off energy for the plane-wave expansion of the wavefunction. The Hubbard correction (GGA+U) was employed to account for localized orbitals due to the presence of 3d/4 f orbitals. The effective U corrections are obtained using an effective U parameter U eff = U − J as outlined by Dudarev et al. [11].
The Hubbard U parameter allows for the self-consistent determination as outlined and implemented by Cococcioni et al. [12]. Since DFT treats orbitals as delocalized, partial occupation is favoured. Hubbard U correction can be considered as the energy cost to enforce full occupation of the Hubbard sites. Cococcioni et al. formulated a linear response approach by perturbation of the occupation matrices for the Hubbard sites. By introducing a finite perturbation in the occupation matrix, the gradient of energy with respect to orbital occupation may be established. Hence, the energy cost of full occupation may be obtained, which is the effective Hubbard U parameter.
Calculations were performed using GGA+U with selfconsistently determined values for U eff . These were obtained for every inequivalent atomic sites occupied by the Nd and Fe. The Nd or Fe atoms occupying each inequivalent site were perturbed separately and U eff obtained for the corresponding Hubbard site. The magnitude of U eff is related to the atom-like nature of electrons occupying the Hubbard sites.
The enthalpy of formation for different compounds were calculated with respect to (αNd), (αFe) and (αB) using the following equation.
In order to compute the vibrational contribution to C p of the compounds, phonon frequencies were computed using the quasiharmonic approximation (QHA), as implemented in the phonopy [13] code. The approach uses finite displacements in order to compute force constants, which are subsequently used to construct the dynamical matrices and phonon frequencies. Vienna Ab Initio Simulation Package (VASP) [14][15][16][17] was used as the force calculator, since it provides more reliable fixed-volume optimization for QHA purposes. Phonon calculations are performed for these compounds using varying supercell sizes. Choices regarding computed supercell sizes are made with respect to computational cost and the number of calculations needed over the period of study.
There are more than one formulations of specific heat (C p ) as implemented in phonopy [18]. The phonon contribution to Helmholtz free energy is computed as in Equation (2): where ν denotes band index and q denotes the q-points. The Gibbs energy is used in order to derive properties at constant pressure. It may be calculated as the minimum on a volume V, as seen in Equation (3): which requires the computation of the Helmholtz free energy at multiple volumes. This is applied over a volume range of −5% to +5% of the equilibrium volume, with 1% increments (a total of 11 points). C p can then be calculated from the second derivative of Gibss free energy with respect to temperature: The result of Equation (4) is collated in the file Cp-temperature.dat within phonopy code. A different derivation may produce a smoother result in the Cptemperature polyfit.dat, as described in Equation (5): where S denotes the entropy, again taken at minimum with respect to volume V at temperature T . It is common practice among CALPHAD researchers to obtain high temperature data of up to T = 3000 K for a complete database. For visualization of crystal structures, the XCRYSDEN code [19] was used.

III. RESULTS AND DISCUSSION
A. Elements

Nd
The ground state crystal structure of Nd (αNd) belongs to P6 3 /mmc space group. Nd atoms are arranged in a double hexagonal close-packed (dhcp) arrangement with two distinct symmetry-irreducible positions with an antiferromagnetic ordering. In order to account for the antiferromagnetic ordering in the structure, the initial magnetic moments for the two positions are specified as antiparallel to one another. The crystal structure is shown in Figure 1. Initial data on the lattice parameters and atomic positions are taken from the work of Wyckoff [20] (Table I), along with the calculated U eff .

Fe
The ground state crystal structure of Fe (αFe) belongs to the Im3m space group. The crystal structure is shown in Figure 1, with a ferromagnetic ordering. The initial lattice constants were taken from Kohlhaas et al. [21]. The details of the crystal structure is listed in Table II. A single Hubbard site was treated with U eff of 3.9182 eV.
From geometry optimization using GGA, we obtained an underestimated lattice constant of around 2.8253 Å. This result is comparable to a host of past ab initio calculations of ground state Fe, which shows a slight underestimation [22,23]. However, from geometry optimization using GGA+U an overestimated lattice constant of around 2.9127 Å was obtained, leading to a large difference in ground state energy between the two optimized structures. This would play a big part in the enthalpy of formation calculations for Nd 5 Fe 2 B 6 , discussed in section III F.

B
There is uncertainty regarding the ground state structure of boron, whether it is αB or βB. Both αB and βB belong to the space group R3m. While αB was found to be energetically more stable in an earlier work by Shang et al. [24], another study found βB to be more stable by an energy difference of  [25]. The αB is regarded as the ground state in this work due to this marginal difference in energy and its simpler structure compared to βB. The initial lattice parameters and atomic positions required for the calculations are taken from the experimental work of Will and Kiefer [26]. These are listed in Table III. A rhombohedral primitive cell was constructed based on the crystal structure data found in their work. It contains 12 atoms (Figure 1), compared to the 105 atoms in the case of βB.
B. NdB 6 NdB 6 crystallizes with a cubic crystal structure. Its space group is is Pm3m (Figure 2). The initial lattice parameters are taken from McCarthy and Tompson [27], while the atomic positions are from the prototype structure CaB 6 [28]. The crystal structure data is listed in Table IV. No significant change in lattice parameters or atomic positions were observed after the geometry optimization.
As a part of our initial calculations, the enthalpy of formation of NdB 6 was calculated using several methods. One of the calculations was done with GGA-PBE exchange-correlation functional without any correction. The other calculations involve Hubbard U corrections and with inclusion of the J term (GGA+U +J). The last calculation was performed with the effective U parameter (GGA+U eff ). This approach was done in order to establish the reliability of the effective Hubbard correction, in comparison to the non-corrected GGA and stricter Hubbard corrections.
The GGA+U + J calculation was performed with uniform correction parameters. The U value was taken as 4.8 eV and the J value as 0.6 eV, both taken from a cRPA (con- strained random phase approximation) calculation by Nilsson et al. [29]. The effective Hubbard parameter U eff was estimated to be 5.337 eV using the method proposed by Cococcioni [12]. The calculated enthalpy of formation using different methods are listed in Table V, along with the experimental data from Storms [30]. The enthalpy of formation calculated using GGA+U + J and GGA+U eff methods are in good agreement with the ex-perimental data. GGA-PBE calculations show slightly more negative value. These results show that while conventional GGA is sufficient to predict the enthalpy of formation of NdB 6 to a certain extent, both the simplified and strict Hubbard corrections shows improvement. The Hubbard correction allows for a better representation of the localized orbitals on the Hubbard site (atomic-like 4 f orbitals), leading to a more reliable prediction of energetics of the system.
The smaller primitive cell of NdB 6 allows for a larger, 2x2x2 supercell to be used in the phonon calculations. We would expect the best result for this supercell relative to the other compounds within the scope of this work run with 1x1x1 supercells. In addition to the required C p data, we have also obtained the 1x1x1 phonon dispersion in order to make a rough comparison with other lanthanide series hexaborides. This is shown in Figure 3. Comparison is made between the phonon dispersion patterns of NdB 6 and LaB 6 , (Figure 4), as well as NdB 6 [31] Eryigit [31]. It is evident that the phonon dispersions are similar, which is expected due to the similarity in their crystal structures. The dispersion curve includes a recognizable flat acoustic modes across the high q-points, which is indicative of vibrations of large rare-earth atoms between octahedral B atoms. It is clearly seen in all the rare-earth hexaborides RB 6 .
Comparison is made to theoretical prediction based on heat capacity expression suggested by Bolgar et al. [35], as shown in Figure 6. The calculated C p of NdB 6 in this work has good FIG. 5. Comparison between calculated phonon dispersion of 1x1x1 supercell NdB 6 (red) in this work and CeB 6 . Two data sets for CeB 6 are shown, DFPT calculation results by Gürel and Eryigit [31] (black) and the experimental work of Kunii et al. [34] (circles). Modified from Figure 3 of Gürel and Eryigit's work [31] agreement with the measurements of Reiffers et al., with the exception in the very low temperature region T < 20 K (see Fig.1 in reference [36]). Magnetic transitions that are not captured by QHA, largely contributes to C p in this temperature region.
We see in Figure 6 that the formulation of C p from Equation (5) is quantitatively similar to those from Equation (4), and is qualitatively less prone to numerical fluctuation. This leads us to believe that Equation (5) is more suited to model C p . Therefore, the C p derived from Equation (5) is favored and is used from this point on in this work.
A major shortcoming of our calculation is the instability of the numerical results from T > 1500 K, the initial stages of which may also be seen in Figure 6 (more visible for the results derived from Equation (4), although for both formulations instability arises in even higher temperatures). This high temperature unreliability is due to the approximations made when computing the phonon frequencies, notably the periodic boundary condition.
The frozen phonon approach uses small atomic displacements to calculate the force constant matrices. These displacements are repeated over the periodic boundary condition however, which negatively impacts the reliability of calculation. Using larger supercells diminish this effect, and high temperature regions will amplify the resulting unreliability as temperature is a factor in both specific heat formulations. It is clear that the 2x2x2 supercell used in this calculation is insufficient to produce reliable values for over T > 1500 K. For other compounds in this work which utilized 1x1x1 supercells, only values for T < 1000 K are considered reliable for CALPHAD. This falls short of the initial expectation to provide C p values up to T = 3000K.  The crystal structure of NdB 4 belongs to the P4/mbm space group. It can regarded as a structure obtained from a distorted cubic NdB 6 lattice. The structure has B atom pairs between distorted cubic sublattices of Nd atoms and B octahedras. It exhibits interesting magnetic properties [37]. Details concerning its crystal structure are given in Table VI. There is only one Hubbard site in the NdB 4 unit cell. The initial lattice parameters and atomic positions are taken from Salamakha et al. [38].

Equation (4) Equation (5) Estimation
The geometry optimization did not result in significant changes in the lattice parameters and atomic positions from the initial values. The calculated enthalpy of formation for NdB 4 are listed in Table VII, along with the experimental value from Meschel and Kleppa [39]. In the GGA+U + J calculations, the values of U and J are from [29] as in the case of NdB 6 . Similar to NdB 6 , GGA+U + J and the simplified GGA+U eff results are in closer agreement with the experimental data. From these calculations it is evident that simplified  GGA+U eff method is quite reliable.
Novikov et al. [40] studied the tetraborides of rare-earth elements lanthanum, dysprosium, holmium, and lutetium for 2 K < T < 300 K. It was shown that the C p value for rare-earth tetraborides is up to 80 J/mol/K at 300 K. The calculated C p for NdB 4 using QHA is shown in Figure 7. It is evident that it corresponds to the expected range for a rare-earth tetraboride as stated. These values are also compared with the results of Bolgar et al. [35].
Calculated values are also comparable to the low temperature measurement performed by Watanuki et al. [41], and there is good agreement of both for T > 20 K. Small peaks in T < 20 K seen in the experimental results, was not reproduced in our calculations. As is with NdB 6 , these peaks correspond to the complex magnetic ordering in NdB 4 reported in other works, most recently by Yamauchi et al. [37]. Figure 7 shows that the QHA result underestimates C p from the experimental result for T > 298.15 K. There is also a slight underestimation for T > 20 K from comparison with low temperature measurement listed above. This likely indicates that the lattice contribution alone cannot accurately represent the total value for C p . The electronic contribution is also required in order to obtain reliable values for the C p of this compound. Crystal structure of Nd 2 B 5 belongs to C2/c space group. It has a monoclinic unit cell, similar to Gd 2 B 5 . The initial lattice parameters and atomic positions are from Roger et al. [42] ( Table VIII). There are two Hubbard sites in the Nd 2 B 5 . The calculated enthalpy of formation of Nd 2 B 5 calculated using different schemes are listed in Table IX. The enthalpies from GGA and GGA+U eff calculations are in good agreement with the experimental data from Meschel and Kleppa [43].
Despite having only slightly distorted structure than NdB 4 , experimental enthalpies are significantly different. However, the computed values are not able to reflect this difference in enthalpy. The same result has been obtained in the recent DFT calculations for the Nd-B binary system by Colinet and Tedenac [4]. Similar difference were also seen for isostructural Gd 2 B 5 , investigated in their work. This discrepancy between experimental and DFT results suggests a characteristic of the R 2 B 5 structure which is not being captured by either GGA or GGA+U calculations.
The calculated C p for Nd 2 B 5 is shown in Figure 8. It is worth noting that the stoichiometry R 2 B 5 is not common among rare-earth borides.  of Crystallographic Data for Intermetallic Phases [45]. Both Nd and Fe sites are treated as Hubbard sites, with five distinct U eff parameters for each site. These parameters, as well as the crystal structure data, are listed in Table X. This structure contains three formula units (57 atoms). A rhombohedral primitive cell containing one unit formula (19 atoms) was used for the calculations in this work.
The calculated enthalpy of formation for Nd 2 Fe 17 and comparison with experimental data is given in Table XI. It is evident that Hubbard correction results in significant deviations from regular GGA calculations. We apply the same treatment for the enthalpy of formation as applied to the compound Nd 5 Fe 2 B 6 (see section III F).
Unlike Nd 5 Fe 2 B 6 , only treating Nd as Hubbard sites (GGA+U eff , Nd) still results in an overestimated value of enthalpy of formation. This implies that the bonds that occur within the structure is sufficiently delocalized such that the introduction of Hubbard U correction does not properly reflect the electronic state. This suggests that the 4 f electrons in Nd participate in the bonding within Nd 2 Fe 17 sufficiently to lose its atomic-like nature. We postulate that the abundance of coordinated Fe atoms around the Nd atoms as the source of this behavior. In this compound, Nd atoms are present in small amounts and always surrounded by 6 Fe atoms, a unique situation compared to other compounds investigated in this study. This might have lead to a hybridization between the 4 f electron in Nd with one of the 3d levels in a coordinated Fe atom. Such a mechanism would explain how the loss of atomic-like nature of the 4 f orbitals and why the Hubbard U correction does not reflect an accurate electronic structure. This is a likely explanation, especially given that Nd 4 f electrons only slightly contribute to bonding in metallic Nd [46], in contrast with most other lanthanides.
With this in mind, we calculated the enthalpy of formation of Nd 2 Fe 17 by regular GGA, and adding the effective Hubbard U correction only for ground state Nd. Instead of the uniform treatment of Hubbard correction for all phases involved in the calculation (GGA+U eff ) or uniformly using regular GGA (GGA) results, only ground state Nd is calculated with the Hubbard correction.
This resulted in a slightly negative enthalpy of formation, as shown in Table XI. It agrees well with the experimental data Meschel et al. [47]. The calculated C p using QHA for this compound is shown in Figure 9. As such there are no experimental measurements available for comparison.
F. Nd 5 Fe 2 B 6 Nd 5 Fe 2 B 6 phase is one of three stable ternary compounds in Nd-Fe-B system. It is also referred to as T3 phase [48]. It has a rhombohedral structure with R3m space group. A primitive cell with 13 atoms was used in the present work as input structure. The details about the unit cell, containing 3 formula units (39 atoms), is summarized in Table XII. The initial lattice parameters and atom positions were taken from Yartys et al. [49].  For this compound, implementing the Hubbard correction to Fe results in a much larger shift of ground state energy for the elements and far less so for the Nd 5 Fe 2 B 6 structure itself. This imbalance resulted in a large shift in the calculated enthalpy of formation. Another GGA+U calculation on Nd 5 Fe 2 B 6 is performed by treating the Nd sites only as Hubbard sites. We take the value of U eff for Fe to be zero. The calculated enthalpy of formation is listed in Table XIII. It can be seen that the modified GGA+U eff gives better agreement with the enthalpy obtained from GGA-PBE value. The inclusion of Fe as Hubbard sites overestimates the enthalpy value. The calculated C p for Nd 5 Fe 2 B 6 using QHA is shown in Figure 10. However, without any available experimental results for comparison, caution must be exercised in its eventual use in the Calphad modelling of the system.

IV. SUMMARY
We have performed DFT calculations on several stable binary and ternary compounds in Nd-Fe-B and its constituent binary systems in order to obtain values of enthalpy of formation and C p . The simplified Hubbard U correction U eff was adopted in order to account for the localized 4 f orbitals present. The results of our calculations are compared with theoretical predictions and experimental measurements from literature. It is evident that Hubbard U correction, treating Nd atoms as Hubbard sites, successfully corrects for the regular GGA method on the prediction of enthalpy of formation for several compounds. Phonon calculations were performed in order to obtain lattice vibrations contribution to C p . It was seen that the results achieve good agreement with low temperature measurements except for magnetic ordering in T < 20 K. These computed values will serve as useful input for CAL-PHAD type optimization.  The ground state of Fe converges to ferromagnetic ordering. As previously stated, GGA+U eff geometry optimization of ground state Fe results in an overestimation of the lattice constant to 2.9127 Å. Meanwhile, GGA geometry optimization results in a slight underestimation of the lattice constant of 2.8253 Å, in good agreement with previous ab initio works which utilized the GGA-PBE exchange correlation functional [22,23]. B atoms are not counted as Hubbard sites and are nonmagnetic. As such, geometry optimization are performed as a non-spin polarized calculation with regular GGA, with the results shown in Table XVI. Optimized calculation parameters are as follows: Cutoff energy (wavefunction) : 40 Ry Cutoff energy (electron density) : 320 Ry k-mesh : 9 × 9 × 9 Maximum smearing width : 0.01 Ry Unlike ground state Nd, we found for NdB 6 ferromagnetic ordering to be energetically more favorable compared to antiferromagnetic ordering. GGA+U eff geometry optimization produces results shown in Table XVII  Optimized calculation parameters are as follows: