High-throughput compatible approach for entropy estimation in magnetocaloric materials: FeRh as a test case Journal of Alloys and Compounds

Aiming to predict new materials for magnetic refrigeration from high-throughput calculations asks for an accurate, transferable, and resource-wise balanced approach. Here, we analyze the in ﬂ uence of various approximations on the calculation of key properties of magnetocaloric materials, while revisiting the well-known FeRh system for benchmarking our approach. We focus on the entropy change and its contributions from the electronic, lattice, and magnetic degrees of freedom. All approximations considered are based on ﬁ rst-principles methods and have been tested, and compared for FeRh. In particular, we ﬁ nd that in this context, the Debye approximation for the lattice entropy fails, due to the presence of soft phonon modes in the AFM phase. This approximation is frequently used in the literature as a simple alternative to full phonon calculations. Since soft modes are likely to occur also among promising magnetocaloric materials where structural transformations are common, the use of the Debye approximation should be discarded for these systems treatment. This leaves the calculations of the lattice contribution the most demanding task from the computational point of view, while the remaining contributions can be approximated using more ef ﬁ cient approaches. The entropy change D S shows a peak around 370 K, for which the total entropy change is given by 24.8 JK (cid:1) 1 kg (cid:1) 1 ( D S ele ¼ 7.38, D S lat ¼ 7.05, D S mag ¼ 10.36 JK (cid:1) 1 kg (cid:1) 1 ) in good agreement with magnetocaloric materials show magnetostructural or magnetoelastic transitions indicating strong coupling between lattice and magnetic degrees of freedom. A consequence of neglecting these


Introduction
The idea of replacing convectional room temperature cooling devices by solid-state magnetic devices, which have the potential for better energy efficiency without producing harmful greenhouse gases, has promoted the interest in magnetocaloric materials. The search for new materials with a more attractive performance/cost ratio or tuning of known compounds is crucial in order to use such devices in mass production and everyday applications [1e5].
First principles high-throughput calculations can be a powerful approach to identify suitable candidates with desired properties by screening a large body of data. To be able to do that, large databases and screening parameters, which are carefully selected to achieve a balance between the accuracy and the cost of the computation are required [6,7].
The performance of materials used in a magnetocaloric cycle can be characterized by the refrigerant capacity, RC ¼ DS iso DT adi , where DS iso is the isothermal entropy variation and DT adi is the adiabatic temperature change. None of these parameters can be easily estimated from first principles electronic structure calculations, instead, a more tailored approach is necessary to take into account the finite temperature effects. Analogous to Ref. [8], we propose the use of the entropy variation between the involved magnetic phases (DS) as an approximation of DS iso . In this way, the transition itself and the field contribution are not included on the description, simplifying considerably the calculation. The magnetocaloric effect is in general small unless it is operated at temperatures in the vicinity of a phase transition, whereas it is strongly enhanced by the entropy variation of a transition, which justifies our approach [2,9]. In a simplified model, entropy can be described by the sum of three independent contributions: the electronic entropy, the magnetic entropy and the lattice entropy: S ¼ S ele þ S mag þ S lat . This is a simplification of the real processes, since most of the magnetocaloric materials show magnetostructural or magnetoelastic transitions indicating strong coupling between lattice and magnetic degrees of freedom. A consequence of neglecting these coupling terms (or taking them as constants) is that, their contribution are "double-counted" when summing the three contributions. However, as shown here, the simplified approach still provides a reliable estimation for DS without overburdening the calculations [1,10].
By using DS as a screening parameter, we are likely limiting our search to materials with first order transitions, since they have enhanced entropy variation [2,4,10]. These materials show better magnetocaloric performance but can also be more challenging to operate in practice, due to hysteresis losses. As pointed out in Ref. [4], first order transitions have hysteresis that can reduce drastically the performance in multi-cycle processes and thus make the materials less attractive for real cooling devices applications.
However, even with the above-mentioned limitations in mind, DS is a natural choice for screening potential magnetocaloric materials, when attempting high-throughput approaches.
In order to be used in high-throughput calculations we need to explore the degree of complexity needed to get reliable estimations for the different contributions for DS. Therefore, FeRh, a wellknown magnetocaloric system, is used as a test case keeping in mind that the approach should be as general as possible in order to be transferable to other systems. Starting from simple models, the different conventional approaches are compared relatively to their performance and applicability for high-throughput calculations. We would like to stress that the focus of our study is on the methodology used for first-principles entropy estimations and not on the test material, FeRh, itself that was chosen by thorough studies available in the literature [8,11e22].
Over the years, the unusual metamagnetic first-order transition of ordered FeRh alloys with CsCl structure has caught huge attention, which is reflected, in a larger number of experimental and theoretical studies [11e13, 18,23e31]. An isostructural transition from a lowtemperature antiferromagnetic (AFM) phase to a high-temperature ferromagnetic (FM) phase occurs near room temperature (around 340K), accompanied by a volume increase of about 1%. The transition is also characterized by the considerable gain of the magnetic moment of the Rh atoms (z 1 m B ) from a nonmagnetic state in the AFM phase (type G), which stabilizes the FM phase [13].
Early attempts to determine the origin of the transition e.g. using the exchange-inversion model of Kittel [32,33] were incompatible with the large entropy variation observed in FeRh. Based on the measured electronic contribution to the entropy variation Tu et al. proposed that the transition might be driven by changes in the electronic structure [23], however this explanation did not compare to previous results for Ir doped FeRh [34]. Later, it was proposed by Gruner et al. that the transition is driven by magnetic fluctuations [13], and the same conclusion was obtained by Gu et al. [18] and Staunton et al. [16] using different approaches.
Nowadays, there is a renewed interest in these compounds due to their magneto-and barocaloric properties. Examples of such studies are e.g. the performance of the magnetocaloric effect (MCE) under cyclic conditions [29] and the variation of the magnetocaloric response between FeRh based ternary compounds [35]. Very recently, the existence of an orthorhombic low-temperature phase of FeRh has been predicted from first-principles calculations [12,14,36] as well as a martensitic transformation under strain [14,15,37].
The existence of such broad knowledge and detailed information in the literature together with the complex metamagnetic behaviour that demands a careful treatment makes FeRh an ideal test system for our purpose to identify a method that can be applied in a high-throughput study for finding new magnetocaloric materials.
The discussion in the present work is divided in two parts. In the first part, we discuss thoroughly the single entropy contribution in terms of electronic, lattice and magnetic components. This is done for FeRh using different approximations, albeit without considering thermal effects on the structure. In the second part, we include volume expansion/contraction from thermal effects and compare with the previous results, using the approximations we found to be adequate to describe the system. From this we are able to conclude which is the most viable approach to be applied in high-throughput calculations.

Computational details
The structural properties as well as structure relaxations were performed using the VASP (PAW) code [38e40] with PAW potentials [41,42] while the PHONOPY [43] code was used to obtain the vibrational density of states and the phonon spectra. Magnetic and electronic properties needed for entropy calculations were derived from a full-potential linear muffin-tin orbital method (FP-LMTO) using the RSPt code [44], and respective temperature dependent quantities such as the adiabatic magnon density of states or the Curie temperature were computed using the UppASD code [45]. In all the DFT calculations, the functional GGA-PBE [46] was used, since it shows in general a good performance in transition metals and compounds, which represents the substantial part of the future database to screen.
Both FM and AFM phases were relaxed on cubic cells of 16 atoms (8 f.u.) taking 4s, 4p, and 3d for Fe as well as 5s, 5p and 4d orbitals for Rh as valence states. A kinetic energy cutoff of 500 eV, roughly 2 times bigger than the default value, was used. For sampling the Brillouin zone we used a k-mesh 12 Â 12 Â 12 generated with the Monkhorst-pack scheme in combination with a smearing of 0.05 eV according to the Methfessel-Paxton scheme (2 nd order). Tests with the inclusion of the Fe 3p and Rh 4p semi-core states in the valence, as well the usage of a higher cutoff energy (750 eV) revealed that the calculations are converged with respect to these parameters. The relaxed lattice parameters of 2.99 Å (AFM) and 3.01 Å (FM) are in good agreement with previous calculations, e.g. 2.99 (3.01) Å [12], 3.00 (3.01) Å [8] for the AFM (FM) phase. They are also in good agreement with experimental measurements, 3.00 Å [47], and 2.98 (3.00) [48] [11] in the AFM phase. At T ¼ 0K, the AFM phase is lower with 26.9 meV/ atom (VASP calculation 1 ) compared to the metastable FM phase.
The phonon calculations were performed within the harmonic approximation employing the finite displacement method in a similar setup as used on the structural relaxation. We used displacements of 0.01 Å for these calculations. A 2 Â 2 Â 2 supercell from the relaxed structures with 128 atoms (64 f.u.) was employed. For this supercell we used a coarser k-mesh 6 Â 6 Â 6. No improvement was observed by increasing the cutoff energy to 750 eV and neither a significant change of the phonon spectrum with the inclusion of 3p (Fe) and 4p (Rh) orbitals in convergence tests performed for cells of 16 atoms. To make it easier to compare results with previous calculations [8,12], we employed also a similar setup, with a cutoff of 500 eV and the inclusion of the semi-core states in all phonon calculations.
For calculations performed with the RSPt code, we used fcc-like structures of 4 atoms (2 f.u.), with the previously relaxed lattice parameter on a 36 Â 36 Â 36 k-mesh with related integrated quantities broadened by Fermi smearing of 1 mRy. The exchange parameters J ij were calculated using the Liechtenstein method [51,52], implemented in the RSPt code, as described in Ref. [53]. The Curie temperature (T C ) calculated via mean field theory according to the obtained values of J ij for the FM phase is of 804 K, which is comparable to the experimentally measured value of 675 K [33,47]. This agreement is good, given the fact that mean-field theory tends to overestimate T C about $20% (as discussed e.g. in Ref. [45]). For analysis of the long range behaviour with distance, the exchange parameters were calculated on a denser k-mesh of 64Â 64Â 64 to assure convergence of the results.

Results
As initial approach, anharmonic effects raised by thermal expansion were neglected and we consider only the DFT groundstate volumes for both magnetic phases. We extend the use of this terminology for elastic/structural properties for this approach to distinguish clearly that the volumes were fixed. The assumption of purely harmonic forces between atoms is insufficient to describe the thermal expansion or contraction of a material, and it may be important to consider anharmonic effects, for accurate calculations of phase stability and entropy estimates. To compare improvements obtained by this description, relatively to the previous "Harmonic" approach, we used the quasiharmonic approximation (QHA) to include the effects of thermal expansion on the entropy estimates (see more details further) [54].

Electronic structure, and its contribution to the entropy
The density of states (DOS) of FeRh is shown in Fig. 1a, for the FM and AFM configuration. Note that for the AFM configuration we show the spin-polarized DOS of only one Fe atom.
The figure shows the atom with more spin-down electrons occupied, representing a Fe atom with a negative atomic moment. The Fe atom with positive moment has exactly the same DOS, although with opposite spin-projection to that shown in Fig. 1. In agreement with previous findings in the literature, a strong hybridization between iron and rhodium orbitals is observed [17]. In particular, a strong hybridization between Fe t 2g and Rh e g orbitals near the Fermi energy (E F À E ¼ 2 meV) occurs in the AFM phase, where it also can be assumed some hybridization between Fe e g and Rh t 2g orbitals in the peak around E F À E ¼ 2 meV, see Fig. 1a. For the FM state the hybridization seems to weaken, and be confined on the minority spin channel, mainly observed between t 2g orbitals of Fe and Rh near the Fermi energy. This observation may emphasize the picture of quenched Rh magnetic moments due to the competing influence of neighbouring iron atoms on the AFM phase. The hybridization that is diminished between FeeRh on the FM phase can be directly ascribed to the lifting of the anti-parallel alignment of the surrounding iron atoms. On the other hand, it can also be related to the increased volume, which can reduce orbital superposition or a combination of both effects.
In the FM phase (Fig. 1b), it is possible to distinguish a significant difference between t 2g and e g orbitals of Fe at Fermi level which can be an indication of different magnetic behaviour of these orbitals similar as observed for bcc-Fe in Ref. [55]. There it was found that t 2g orbitals, with likewise bigger contribution for the electronic density of states DOS(ε F ), were related to the long-range Ruder-maneKitteleKasuyaeYosida (RKKY) interactions while the e g were associated with direct exchange with nearest neighbours. The similarities, observed in the projected Fe DOS, might hint for the existence of some similarities between magnetic behaviour of Fe atoms of both compounds.
The contribution of electronic excitations to the entropy is given by the mixing entropy of occupied and unoccupied states: In Fig. 2, the results obtained from both models are compared, showing a good agreement for temperatures till 300K. Outside this range the results differ, resulting in a small deviation observed for temperatures close to the transition temperature for the transition from antiferromagnetism to ferromagnetism. The deviations arise from the AFM phase and can be explained by the absence of peak structures in DOS at the Fermi energy ( Fig. 1) as the Sommerfeld approximation assumes. Nevertheless, the deviation between the discussed models for S ele estimation, in the range of the transition, is small being the obtained difference of 1.84 J K À1 kg À1 for DS ele , and 26 K for the transition temperature (discussed in further sections). Although there is not significant loss of accuracy estimating S ele using Eq. (2), using the definition of mixing entropy does not imply extra computational effort. Thus to avoid eventual inaccuracies that may arise by using Sommerfeld approximation, we use the definition in Eq. (1) as the standard method for calculating the electronic entropy.

Magnetic contribution to the entropy
For materials with order-disorder magnetic transitions, the maximum magnetic entropy variation between phases can be roughly estimated from DS mag ¼ Nk B ln½2S þ1 (in a quantum description) with N being the number of magnetic atoms [24]. This comes about since very few microstates are available for highly ordered states, and the entropy of this configuration can be neglected in the limiting case T / 0. In contrast, for the disordered configuration we have ð2S þ 1Þ N arrangements for the spins for T /∞, which results in the entropy change across the order-disorder transition as described above.
According to the analysis above, it is expected that order-order transitions at finite temperature have a considerably smaller entropy change from the magnetic subsystem.
Based on this and the argument that considering the itinerant nature of magnetism of FeRh, the magnetic contribution to DS is already included in the electronic entropy computed from the DOS, some reports argue that the magnetic contribution of the entropy does not need to be considered separately [8]. To some extent, it is possible that some magnetic contribution is captured by DS ele , since some coupling between the degrees of freedom is expected. However, taking into consideration the increase of the Rh magnetic moment from 0 (AFM) to z1m B (FM), it seems that the magnetic entropy for this transition must be considered specifically. A good reason for that are the new two-site interactions between Fe and Rh atoms of the ferromagnetic phase, described e.g. by the Heisenberg Hamiltonian, which should be considered for a proper system description (see Fig. 3a). Interestingly, Fig. 3a shows that the close range of the FeeFe exchange is quite similar for the AFM and FM configurations. The nearest neighbour interaction is antiferromagnetic in both phases, although the strength is larger for the AFM configuration. In addition, the general trend of the FeeFe interaction is quite similar for both configurations. The interaction that stabilizes the FM phase is hence not found in the FeeFe Heisenberg exchange. Instead, as Fig. 3a shows, the strong ferromagnetic FeeRh interaction is what makes the FM configuration stable at all. This represents an interesting boot-strapping effect, when the FM configuration is what allows for a sizeable Rh moment, and the sizeable Rh moment is what ensures a large FeeRh exchange interaction that makes the FM (meta-)stable [13].
To estimate the magnetic entropy variation we started by using a simple approximation analogous to the one used in Ref. [57] for LaFe 13Àx Si x alloys. From the fundamental thermodynamic relation dU ¼ TdS À PdV one can, for isochoric processes, approximate the entropy as DS ¼ DU=T. Although crude, this approximation should give an acceptable estimate for the entropy variation in first-order transitions where the entropy varies discontinuously at the transition temperature. Using this and describing the magnetic energy of the system by the classic Heisenberg Hamiltonian:    [13,25,58,59] that an extension of the Heisenberg exchange model can be made, using e.g. higher-order interactions, to obtain a proper magnetic description of the system with very satisfactory results on Monte Carlo simulations. Here we took a different route to avoid the use of a tailored model and evaluated the magnetic entropy from spin-wave fluctuations, similar to the work in Refs. [18]. This approach is possible since the AFM-FM transition happens at considerable lower temperatures (z340 K) than the Curie temperature of the FM phase and the spin fluctuations can still be considered to be relatively small [12]. It is also necessary to guarantee, in order to use this approach that Stoner excitations are not dominating, as was shown in Refs. [17,18]. For this reason, we calculated the magnon density of states (MDOS) from the adiabatic magnon spectrum. This calculation relied on Heisenberg exchange parameters, J ij , estimated from DFT calculations. Due to the bosonic nature, the entropy of the magnons is given by: where gðεÞ is the MDOS and n ¼ ½expðε=½k B TÞ À 1 À1 is the Bose-Einstein distribution. For these calculations a perfectly aligned configuration was assumed for the spin moments when calculating the magnon dispersion (m ! k m ! z ). Analogous calculations performed for a thermally relaxed (at 300K) magnetic configuration do not deviate significantly from these results.
In contrast to the observation of the electronic entropy contribution, a magnetic entropy maximum (10.92 JK À1 kg À1 ) is obtained at around 315 K, as it can be seen in Fig. 4. This peak is of major importance since it hints to the existence of the phase transition. At least, it shows that the magnetic entropy will favour the ferromagnetic phase. Also, the lack of similar peaked behaviour around the transition temperature in the other entropy contributions (see discussion of lattice contributions, below) suggests that the transition is triggered by the magnetic features of the system. Thus, based on the applied magnetic model (spin-wave fluctuations) and the obtained MDOS for the AFM phase, we suggest that at low temperatures the Rh atoms are magnetically suppressed by the anti-parallel alignment of the surrounding Fe atoms configuration. This generates a vanishing local Weiss field on the Rh atom, that at the transition temperature has its symmetry broken by the spin fluctuations, which allow Rh to become magnetically polarized and thus stabilizes the FM configuration. A similar picture of a transition driven by small magnetic fluctuations, as the one described above, is concluded in other works, both with similar methods [18] as well as from different approaches [13,16,18,58].
It was already pointed out in the previous section, that iron atoms in FeRh might possess some features that are similar to the features they have in elemental form (bcc Fe). For instance, the existence of oscillating long-range interactions, must be taken into consideration when calculating the MDOS and thermodynamic properties from it. In Fig. 5, it is shown how the entropy peak varies with the range of magnetic couplings J ij included in the calculation of the MDOS 2 . DS mag and the peak temperature show a significant dependence on the cutoff radius for the J ij such that a considerable long range of interactions must be included to a fairly converged estimation. This is a consequence of long-range magnetic interactions of J FeÀFe , visible on Fig. 3b) that oscillate significantly till a range of 10 primitive cells approximately. Based on these observations we included interactions up to 12 lattice parameters in our calculations.  The fact that the sensitivity with respect to the range of the interactions, is mainly due to interactions in the FM phase, reflect that long-range oscillating interactions are stronger on this phase (see Fig. 3b), and more significant on t 2g orbitals, underpins the similarity between bcc Fe and the FM phase of FeRh (see results for bcc Fe in Ref. [55]), when it comes to understanding the Heisenberg exchange. We point out that the sensitivity on the cutoff of the Heisenberg interactions is important for many prospective magnetocaloric materials, since many of them are metallic and have Fe as a key element, and the long-range magnetic interaction between Fe atoms seems to be of particular importance.
Our results of the magnetic entropy change across the AFM -FM transition are in agreement with previous calculations, see Table 1. The quite large difference between our entropy calculation and the results obtained by Gu et al., who used a similar computational approach [18], are most likely caused by the shorter range of exchange interactions considered in their work. This might also partly explain the small deviation between our results and the ones from the models used in Refs. [13,58]. Relative to the transition temperature, the higher result obtained by Staunton et al. stands out from the remaining values [16]. Such deviation might be related with the method itself -finite temperature spin density functional theory is implemented in the disordered local moment approachwhich differs significantly from the other approaches. The calculations of Ref. [16] were done from an electronic structure theory that allows a random distribution of spin-orientations, and therefore neglects short-range correlations. This approach is well established and is argued [16] to describe better the electronic structure at finite temperatures.

Lattice contribution
The calculation of properties related to the crystal lattice can become very demanding regarding computational resources. In order to calculate such properties in an efficient way, it is imperative to minimize the numeric effort by using expedient models, without compromising significantly the accuracy. To verify which approximation is appropriate to estimate the lattice entropy, we compare the results of models of various accuracy and complexity. As contribution for the lattice entropy, only the vibrational entropy was considered.

Debye model
In the Debye model the phonon dispersion relation is treated as where v s is the speed of sound in the material.
Therefore, the vibrational density of states (VDOS) is given by: up to the cutoff Debye's frequency. The entropy then becomes [54]: where Q is the Debye temperature. An important consequence of this model is that at a fixed temperature the variation of the Debye temperatures between phases (DQ) has opposite sign to the respective variation of lattice entropy, DQ=jDQj ¼ À DS lat =jDS lat j, this can be used to understand the nature of the lattice entropy, i.e., if it is collaborative (has same sign) or detrimental (opposite sign) relative to all other entropy contributions. The Debye temperature can be computed as Q ¼ Zð6p 2 nÞ 1=3 v s =k B with the atomic density n and v s being the average velocity of sound in the crystal. For isotropic crystals, the later is approximated as the average value of the shear and longitudinal sound velocity. It is generally expressed in terms of the bulk modulus (B), the density (r), and a correction The correction parameter x depends on the elastic properties of the system. However, in Ref. [60] it was proposed that for a given class of materials x might be universal and can be derived from elastic constants. To verify if this approximation could be used in calculations of magnetocaloric materials we extracted the shear and bulk moduli, from data found in literature [61e65], for magnetic materials with structural transitions (expected to be present in an important class of interesting candidate materials for magneto caloric applications). Fig. 6 shows that any possible linear trend as obtained in Refs. [60] is not reasonable if compounds with different structure types are compared. This makes the approach of Ref. [60] less appropriate for high-throughput calculations and data-mining algorithms. For FeRh, in particular, this approach is specially unappealing considering the only materials with similar structure and properties are alloys very close to the stoichiometric  compound.
We conclude that however inconvenient it is, the elastic properties have to be calculated for each material that one includes in any data set for high-throughput calculations, when searching for new magnetocaloric materials. For FeRh this exercise leads to a very interesting result when comparing the two magnetic phases (AFM and FM). Since the two phases have the same crystal structure, one might naively assume very similar elastic properties for both phases. However, this assumption leads to a lattice entropy contribution of 7.9 JK À1 kg À1 at 328 K, which deviates from a more accurate calculation that takes into account the difference in elasticity of the two systems (discussed more in detail below) that yield a value of À30.1 JK À1 kg À1 at the same temperature. This later approach gets closer to the extracted from calorimetric measurements z À 33 JK À1 kg À1 (328K) using the same model [66].
In order to describe accurately the difference in elasticity of the two phases of FeRh, we evaluated the elastic constants using the RSPt software, for both phases. We used the stress-energy response as described in Ref. [67,68]. The values of C 11 ¼ 194.9 (257.1) C 12 ¼ 194.9 (165.2) and C 44 ¼ 135.3 (115.6) GPa were estimated for the AFM (FM) phase, with qualitative agreement with previous calculations [15].
Comparing the Debye temperatures derived by using the same Poisson ratio (n) for both phases and from the calculation with the different ratios for the AFM and FM phase, demonstrates the sensitivity of this model for lattice entropy to small deviations (Dn ¼ 0:05) of the elastic properties. Taking into account the differences in the elastic behaviour of the AFM and FM phase, the change in the Debye temperature DQ is in good agreement with the experimental results, see Table 2.

The Debye-Grüneisen model
Taking a more sophisticated approach to estimate the lattice entropy, by use of full phonon calculations [8,12], leads different values of DS lat , compared to the findings from the Debye model. This is discussed in detail in the following subsection. To investigate whether a simplified approach can be improved, we first extended the Debye model to the Debye-Grüneisen model, where effects of volume variation are taken into consideration for the lattice properties. The Grüneisen parameter, needed for this model, is calculated from where B 0 is the volume derivative of the bulk modulus. The parameter g is an additive factor, usually taken as g ¼ 1 for low temperatures and g ¼ 2=3 for high temperatures [60,69]. Considering the volume expansion, V AFM /V FM , an increase is obtained for jDQj which implies an increase in the magnitude of the lattice entropy variation comparatively to the previous estimate, and does not lead to theoretical values closer to the observed data.
The Debye model is known to be accurate in the limits T ≪ Q and T [ Q 4 . Outside this temperature interval it is less reliable. It is for such temperatures that the magnetic transition for FeRh happens, which partly explains the difference obtained for S lat using full phonon calculations. As the discussion in the next section shows, the existence of soft vibration modes has a major role in explaining these contradictory results between the simple Debye model and the results from full phonon calculations.

Entropy from full phonon calculations
The presence of soft phonon modes, reaching imaginary frequency, leads to a structural transition, which leads to enhancement of DS lat . Even if these soft phonon modes do not reach imaginary frequency (indicating structural instability) they may provide a hint for possible transition. It is reasonable to expect that a fair amount of magnetocaloric candidates will show this behaviour. Soft modes of the acoustic branch can give raise to energy peak structures in the vibrational density of states at low energies that are not captured by the Debye model and lead to inaccuracies even at low temperatures. Since the Debye model fails to describe the thermodynamic properties of such materials, full phonon calculations must at least be tried in order to compare with more approximate methods, and to assess if more efficient avenues exist for the calculation of the lattice entropy.
The calculated phonon dispersion, displayed in Fig. 7, shows that for most of the reciprocal space, the acoustic modes of FeRh behave quite similarly in both the FM and AFM phase. However, the AFM phase has conspicuous soft modes that even become imaginary as in previous calculations [8,12], which also showed imaginary frequencies around the M point. Such behaviour points to a dynamical instability [8,12].
This result was thoroughly discussed for FeRh in Ref. [12] and a competing low temperature monoclinic structure was proposed. Nevertheless, near the experimental transition temperature between the AFM and FM phase, the structure is known to be cubic,  possibly caused by an entropy driven stabilization of the cubic phase, e.g. as discussed in Ref. [70]. Since the part of reciprocal space that contains imaginary frequencies is very small, as observed by their minor contributions to the VDOS in Fig. 7, its influence on thermodynamic properties is expected to be negligible [12]. We therefore neglected this contribution, to the estimation of thermal properties to avoid numerical complications. The entropy variation derived from the VDOS (using the same expression as Eq. (4) but with gðεÞ as the VDOS) is shown in Fig. 8. It has the same sign and order of magnitude as the electronic contribution. To be precise the difference in calculated lattice entropy is 7.05 J K À1 kg À1 at T ¼ 373 K. 5 Comparing the estimates in Fig. 8 it is interesting to note that the trends for DS mag between Debye model and the full-phonon calculations start to differ around 40 K, when the result for the later approach displays a small entropy peak. Tracing this to Fig. 7 we can relate it to the flattening of the phonon spectra around 0.8 THz for the AFM phase, which explains the small entropy peak obtained for the full-phonons approach as the excitation of the soft phonon branches [12]. The indicated observation also underlines the role of the presence of the soft mode to the failure of the Debye model application in this material. The difference between entropy results for phonon calculations and for Debye model in this material is as remarkable as surprising, especially when considering that it is an isostructural transition we consider. A priori there are no indications pointing to the need of a more complex approach, and it is clear from the calculations discussed here for the lattice entropy, that the applicability of any simplified method, such as the Debye model, should be carefully verified for lattice contributions of the entropy variation. This shortcoming of simplified models, should be taken into consideration when estimating entropy variations of any material.

Total entropy variation
The sum of all hitherto discussed contributions to the entropy, defines the total entropy variation between FM and AFM phases, according to our model. In Fig. 9 (solid line) it is clear that the total entropy difference between the FM and AFM phase has a pronounced temperature dependence. In addition, the figure shows a major, broader peak around $ 370 K with maximum entropy difference, DS max , of 24.8 J K À1 kg À1 .
In Table 3 we list all calculated contributions to the entropy change between the FM and AF phase, at T ¼ 350K. It may be seen from the table that all contributions are collaborative and comparable in magnitude [12].
In the taken approach, both the type of phase transition and stoichiometry of the compound are described without the existence of losses, associated to the transition (e.g. coexistence of phases) and without defects of the material, and with this in mind, it is not unexpected that theory overestimates somewhat the entropy contributions and therefore, in this case, the total entropy change. Besides, it is important to note that experimentally is very difficult to achieve the equiatomic concentration and very close alloys are measured instead, for which the entropy variation varies slightly [35]. Nevertheless, the comparison between the here calculated value of DS and experimental results is quite satisfactory [26,29,34,48,66,71e73] being of the same sign and order of magnitude. Thus, if one has the ambition to make theoretical screening approaches in combination with first principles highthroughput calculations, the level of approximation employed here seems to be the simplest way that is capable of a fairly accurate prediction while maintaining computational efficiency.
Experimentally it is not straight forward (or accurate [8]) to disentangle the entropy contributions. In Table 3 the entropy contributions of the present calculations are listed, together with the ones extracted from calorimetric measurements of Ref. [66]. Cooke et al. [66] extracted the lattice entropy by naively fitting the lowtemperature data to the Debye model. This approach fails for FeRh, as discussed above, and consequently, the estimated huge magnetic entropy contribution, calculated from S mag ¼ S tot À S ele À S lat , is not seen as a realistic contribution. Both the extracted DS lat and DS mag are unusually high in magnitude, compared to the usual total entropy values [4]. These values are also unexpectedly high, given the isostructural nature of the volume expansion and the order-order nature of the magnetic transition. Taking into account these considerations, it is more plausible that the high magnetic entropy variation listed as an experimental values in Table 3, is really due a collaborative sum of all entropy contributions [10].
To associate with certain the peak entropy of Fig. 9 to the AFM / FM phase transition, we compared the free energies of both phases. However, we did not obtain an intersection of the free energies, at least not in the considered range of temperatures (0e500K). This is in agreement with results of Ref. [12], but in disagreement with the data of Ref. [18]. Theoretically, our results imply that no phase transition can be associated to the discussed entropy peak, making it as pertinent and interesting as the minor entropy dip around 40K. To our knowledge, the latter does not indicate any known transition and most likely reflects the soft phonons of the AFM phase. A comparison between our results and calculations in Ref. [8,12] as well as [18] reveal that the later reference achieves a significantly smaller energy difference between FM and AFM states, around 2.80 meV/atom in comparison to our value; 27 meV/atom. This energy is 35.4 meV/atom in Ref. [12] and 29.1 meV/atom in Ref. [8] for similar calculations. When compared to experiments, the value of Ref. [18] is clearly closer to experimental estimates, which lie around 2.7 meV/atom [26,66]. This improvement of the energy difference estimation between magnetic phases seems to be due to the unique exchange and correlation functional used in Ref. [18]. The authors of this work employed the Langreth-Mehl-Hu functional [74,75], which appears to have as a feature the reduction of energy between phases [76], and suppression of the magnetic moment [77]. Although this functional provides reasonable results for FeRh, it is a less tested functional for general investigations that involve a large group of compounds. In absence of a firm test, this functional is difficult to apply in a predictive study. Another possibility for the too large energy difference between the AFM and FM phase could be due to dynamical correlations of the electronic structure.
If the estimated DE 0 from DFT is used to estimate the entropy variation as DS ¼ DU=T (as attempted for DS mag ) we obtain a value of 87.99 JK À1 kg À1 . This strong disagreement with experimental measurements also underlines that DE 0 is not properly estimated by DFT.
An important point from this discussion is the difficulty to predict with certainty the temperature for the AFM / FM phase transition. Instead of comparing the free energies of the phases, it is of interest to take a simpler approach and consider the transition to be caused by thermal energy from T ¼ DE 0 =k B . Using this approach, an estimate of 346 K was obtained in Ref. [8] and 350 K in Refs. [58]. Although this value is within the experimental value, applying the same approach using data from the total energy calculations presented here, or from other calculations [11e13], reveals that this simplified method is very sensitive to the details of the calculations, meaning that its use introduces a non-negligible degree of uncertainty while not describing necessarily the physical picture.

Quasiharmonic approximation
To account for anharmonic interactions we use the QHA, which minimizes at each temperature the volume-dependent Gibbs free energy: In our approach we consider magnetic, electronic and lattice contributions to the free energy FðT; VÞ ¼ E 0 ðVÞ þ F mag ðV; TÞ þ F lat ðV; TÞ þ F ele ðV; TÞ. Contributions to the lattice and magnetic energy were obtained at constant volume via a calculation of the respective density of states (gðεÞ mag=lat ) and Bose-Einstein distribution function (nðε; TÞ): ε gðεÞ mag=lat nðε; TÞdε À TS mag=lat : Note that we added the entropy contribution also to this term. These contributions were evaluated, following the same procedure as before, for a series of volumes, and then fitted by cubic splines. For each phase, 9 volumes were considered, including the relaxed volume, ranging the lattice parameter between 2.98 Å(3.00 Å) and 3.00 Å(3.03 Å) for the AFM (FM) phase. F ele ðV; TÞ was calculated similarly, using instead the Fermi-Dirac distribution function, which includes the temperature dependence of the energy. F ele ðV; 0Þ was used as reference level for the electronic free energy since the DFT ground state energies are already included in E 0 (V). The later energies were fitted by the Murnaghan equation of state [78] for the internal energy.
The linear thermal expansion (LTE) obtained from the QHA is shown in Fig. 10. It may be noticed that there is decent agreement with experimental measurements [47] and theory. Similar to the experimental data, there is in our calculations a jump in the LTE at the magnetic phase transition.
We found that contrary to the electronic and lattice contributions, the magnetic contribution to the free energy opposes the volume expansion. Nevertheless, the lattice is the dominant contribution in the considered temperature range and dictates the thermal expansion, and the observed behaviour arises dominantly from the vibrational properties.
The LTE coefficient a l , can be estimated from a linear fit (Dl=l ¼ aT) in the same temperature range as for the experimental data. From the theory, we obtain a slope of the Dl=l curve that for both the FM and AFM phase is similar to the experimental data. The main difference between theory and experiment is the size of the volume expansion at the magnetic phase transition that is smaller in theory compared to the experimental values [47]. This disagreement is not surprising given the simplicity of the model used and how similar are the values between the magnetic phases.
The volume and temperature dependent free energies of the AFM and FM phases allow for the most accurate estimate of the phase stability and entropy, among the calculations presented in the paper. We compare the QHA and Harmonic results in Table 3, Table 3 Comparison of estimated entropy contributions at Tz 350 K for the harmonic and quasi-harmonic approaches with previous calculations in literature and experimental measurements. It is also indicated for the "Harmonic" (QHA) approach, estimated values at the temperature for which the entropy variation has a peak -T ¼ 373 (316)  together with previously reported data. It may be seen that the total entropy change of the AFM / FM phase transition is almost insensitive to the level of approximation, while for the individual contributions there is a more significant difference between the QHA and Harmonic approximations.
We find that there is compensation of DS lat and DS mag , which vary similarly but in opposite direction, as can be seen by comparing the entropy contributions on Table 3. This is caused by variation of J ij parameters with volume, which decreases for bigger volumes. As reported also in Ref. [19], we also observe that couplings between iron moments are significantly more sensitive to this variation than couplings between iron and rhodium magnetic moments (data not shown).
The predicted volume variation in the AFM is responsible for the loss of the monotonous behaviour of DS lat in the transition range.
Since the linear thermal expansion behaves as experimentally measured, and the entropy peak reassembles more the discontinuity expected for first order transitions, we consider that there is a qualitative improvement of the physical description. Besides becoming sharper, the broader entropy peak shifts from 373K to 316K when using the QHA, as seen in Fig. 9.
Our results are close to previous, similar first-principles calculations, combining the QHA results of Ref. [12] (or the ones from Ref. [8]) for DS ele þ DS lat with the results of Ref. [18]. Although there is a small deviation, it is accurate enough to be used in highthroughput calculations, keeping in mind that a more accurate result would need a more tailored and computationally expensive method. Also, we treated entropy contributions as independent, which is a simplification of the problem.

Conclusions
The aim of this paper is to derive a reliable approach based on first principles calculations to determine the entropy change in materials with first order phase transition that can be used in high throughput studies. Thus we have to balance between computational effort and accuracy and a detailed study, concerning estimates for electronic, lattice and magnetic entropy contributions according to different models, was performed using the wellknown MCE system FeRh as test material. Based on the assumption that we can treat the entropy as a sum of three independent parts, i.e. handle the electronic, lattice, and magnetic contributions to the entropy separately, we tested different approximations for each entropy contribution. It turned out that the entropy, or the entropy change, in our test case, FeRh, is very sensitive to approximations made, e.g., even small alterations on elastic properties between the two magnetic phases need to be taken into account for a reliable estimation of the entropy. This means that the Debye model is not adequate and it should not be considered for the highthroughput applications. Although the simplicity of the Debye model is appealing in terms of computational efficiency, it fails to estimate the vibrational entropy in the presence of soft phonon modes. We believe that this sensitivity regarding the vibrational properties is not exclusive to FeRh. Rather, we expect that many of the magnetocaloric candidates will show similar behaviour, which means that accurate phonon calculations are necessary for a reliable description of the entropy.
From the results of the magnetic properties, we conclude that it is necessary to consider an appropriate cutoff for the exchange interactions when using a spin-wave description due to the possibility of long-range interactions, which can have a considerable influence. However, this aspect should be less relevant in case of order-disorder magnetic transitions where spin-flip like excitations are dominant and the Heisenberg model can be used in combination with Monte Carlo simulations to estimate the magnetic entropy.
We observe an entropy peak around the expected transition temperature raised solely from the magnetic contribution, which allows us to support the picture of a magnetism driven transition as discussed in previous works using different approaches [13,16,18,58].
Although a DS peak is regularly observed in phase transitions, it is necessary to compare the free energies of the phases to associate a DS peak to a phase transition. For FeRh, it was not possible to establish this association due to the overestimation of the energy difference between magnetic phases by traditional DFT. Such difficulty raised our awareness of this limitation in our first-principles approach. Nevertheless, beyond DFT methods offer tools for circumventing this problem, and can be used to improve DE 0 estimation to verify the transition occurrence, if needed.
Adopting the QHA approach allows for a more complete description of the systems and a qualitative improvement of the entropy variation is obtained, by the sharpening of the DS peak as expected for a first-order transition. Despite this, no quantitative improvement of the entropy variation is obtained that justifies the significant increase of computational effort required for this treatment.
Therefore, it can be stated that the "Harmonic" approach balances in a very satisfactory way the accuracy and the computational effort. The obtained results DS ele ¼ 7:38 JK À1 kg À1 , DS lat ¼ 7:05 JK À1 kg À1 , and DS mag ¼ 10:36 JK À1 kg À1 are in good agreement with previous calculations and the total entropy variation DS ¼ 24:78 JK À1 kg À1 is close to the experimental range. This establishes the cornerstones for a reliable entropy estimation at high-throughput scale computations, while allowing for reasonable computational effort that allows to avoid possible pitfalls of the calculations.   10. Variation of the linear thermal expansion with temperature for AFM phase (blue line) and FM phase (orange line) obtained within the quasi-harmonic approximation along with experimental measurements (purple circles) [47]. Ttr and T0tr indicate respectively the estimated and the measured transition temperature. For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.