Ab initio calculation of the magnetic Gibbs free energy of materials using magnetically constrained supercells

We present a first-principles approach for the computation of the magnetic Gibbs free energy of materials using magnetically constrained supercell calculations. Our approach is based on an adiabatic approximation of slowly varying local moment orientations, the so-called finite-temperature disordered local moment picture. It describes magnetic phase transitions and how electronic and/or magnetostructural mechanisms generate a discontinuous (first-order) character. We demonstrate that the statistical mechanics of the local moment orientations can be described by an affordable number of supercell calculations containing noncollinear magnetic configurations. The applicability of our approach is illustrated by firstly studying the ferromagnetic state in bcc Fe. We then investigate the temperature-dependent properties of a triangular antiferromagnetic state stabilizing in two antiperovskite systems Mn$_3$AN (A = Ga, Ni). Our calculations provide the negative volume expansion of these materials as well as the ab initio origin of the discontinuous character of the phase transitions, electronic and/or magnetostructural, in good agreement with experiment.

We present a first-principles approach for the computation of the magnetic Gibbs free energy of materials using magnetically constrained supercell calculations. Our approach is based on an adiabatic approximation of slowly varying local moment orientations, the so-called finite-temperature disordered local moment picture. It describes magnetic phase transitions and how electronic and/or magnetostructural mechanisms generate a discontinuous (first-order) character. We demonstrate that the statistical mechanics of the local moment orientations can be described by an affordable number of supercell calculations containing noncollinear magnetic configurations. The applicability of our approach is illustrated by firstly studying the ferromagnetic state in bcc Fe. We then investigate the temperature-dependent properties of a triangular antiferromagnetic state stabilizing in two antiperovskite systems Mn3AN (A = Ga, Ni). Our calculations provide the negative volume expansion of these materials as well as the ab initio origin of the discontinuous character of the phase transitions, electronic and/or magnetostructural, in good agreement with experiment.

I. INTRODUCTION
Solid-state magnetic materials form the basis of many current and upcoming technologies, and are thus subject to a field of research that is in continuous development. Some examples are caloric refrigeration [1][2][3] and gas liquefaction [4], permanent magnets in engines and other electronic equipment, and spintronic technology including magnetic field sensors, neuromorphic computing, and different devices for storage and processing of information [5][6][7][8][9][10]. The central component for their functionality is the formation and exploitation of complex magnetic structures, such as topologically protected skyrmions [11,12], antiferromagnetic [13] and helimagnetic [14,15] states, and high coercivity in ferromagnets [16,17] and ferrimagnets [18,19]. The fundamental understanding and prediction of this broad range of magnetic states is, therefore, crucial to advance and create new functional mechanisms.
The properties of these magnetic states can strongly depend on temperature. In fact, sometimes their stabilization is linked to the presence of thermal fluctuations [11,20] and their functionality usually involves magnetic phase transitions driven by temperature changes. The development of a first-principles tool describing temperature-dependent magnetic states and phase transitions among them is a challenging task but also a necessity to reliably predict and explain magnetic phenomena in materials. Major work has been done for the description of thermal fluctuations of local magnetic moments, including both their transverse [21][22][23] and longitudinal degrees of freedom [24,25]. Many of these theoretical developments also account for a finitetemperature magnetoelastic or magneto-phonon coupling [26][27][28][29]. However, most of them focus on the paramagnetic state, i.e. a high-temperature limit in which the local magnetic moment orientations are fully disordered away from a description of intermediate temperatures.
In this work we present a computational approach to calculate the ab initio Gibbs free energy of a magnetic material, (1) whose minimization provides the dependence on temperature of the most stable magnetic states for given values of an external magnetic field H and an applied mechanical stress σ αβ , the latter causing a material deformation described by a strain tensor ε αβ . {ê n } are unit vectors prescribing the orientations of the local magnetic moments, located at lattice sites {n} and with magnitudes {µ n }. The central tenet of our approach is to assume that {ê n } vary very slowly in comparison with {µ n } and with the underlying interacting electron system. This means that an ab initio magnetic energy can be specified for different configurations {ê n } and that G tot can be obtained by carrying out ensemble averages over {ê n }, denoted by . . . . For example, the internal magnetic energy is where P ({ê n }) ∝ exp[−βE int ] is a Boltzmann probability distribution for the local moment orientations. An associated magnetic entropy forms the second term in the right hand side of Eq. (1) and is given by S mag = −k B ln [P ({ê n })] , k B being the Boltzmann constant.
The novel aspect of our work is that the thermodynamic averages to compute G tot are performed as superpositions of supercell ab initio calculations that are magnetically constrained to different sets of {ê n }. Figure 1 shows magnetically constrained supercells for a thermodynamically excited ferromagnetic state of bcc Fe using this method, where higher and lower temperatures correspond to a set of largely disordered and ordered local moment orientations {ê n }, respectively. Most importantly, Eq. (1) is a general expression that we can compute for an arbitrary probability distribution. In our approach, therefore, different theories for P ({ê n }) can be constructed by restricting the thermal configurations to a meaningful magnetic phase space described by a trial magnetic Hamiltonian that can include the effect of nonlocal correlations, which will be explained in section II [see Eq. (9)]. Another important outcome of this work is the demonstration that the statistics of noncollinear magnetic configurations can be described even with a low number of constrained calculations for supercells of relatively small size, thus being computationally affordable.
A key aspect of our approach is to provide the dependence of G tot on an expansion in powers of magnetic order parameters akin to a Ginzburg-Landau theory. This describes a hierarchy of local moment correlation functions [30], accounting for pairwise and multisite magnetic interactions, and their magnetostructural coupling (see section II C). Higher than pairwise free energy terms can generate discontinuous and ordered-to-ordered magnetic phase transitions, fully quantified and provided by our approach. As a demonstrator of this, we study the noncollinear, geometrically frustrated, magnetism of Mn-based antiperovskite materials [31], which are driving strong interest owing to a discontinuous phase transition to a triangular antiferromagnetic state with a giant barocaloric effect [32][33][34]. We provide the ab initio origin of this transition, how it depends on the chemical composition, and compare our results with experiment. This paper is organized as follows. In section II we set the theoretical framework and present our approach. Section II B explains how to compute Eq. (2) using magnetically constrained supercells. In sections II C and II D we describe the effect of multisite magnetic interactions and magnetovolume coupling on the character of magnetic phase transitions. Computational details are given in section II E. In section II F we apply our approach to the well studied ferromagnetic state of bcc Fe as an instructive example, while section III focuses on the magnetism of antiperovskite materials. Conclusions and outlook are given in section IV. The paper finalizes with appendix A, which shows the performance of the Gibbs free energy calculation using smaller supercells.

II. AB INITIO THEORY OF THE MAGNETIC GIBBS FREE ENERGY
A. Theoretical framework The basis of our ab initio modelling is density functional theory (DFT), formally extendable to finite temperatures [35]. DFT describes the local moment orien- Magnetically constrained supercells for different local moment configurations {ên}, represented by the arrows, for the case study of ferromagnetic bcc Fe. Different temperatures correspond to different orientational probability distributions P ({ên}). The method is illustrated here for P ({ên}) fully prescribed by single-site internal magnetic fields {βhn = βh}, or equivalently by single-site local order parameters {mn = m}, where βh is set alonĝ z (see section II B 1). Note that these quantities are siteindependent for a ferromagnetic state. The smaller the value of m, the larger the magnetic thermal fluctuations. Results are shown for {m = 0, 0.50, 0.90, 0.99}, which correspond to {βh = 0, 1. 80,9.99,99.94} via Eq. (20). Red and blue arrows are used for orientations with a polar angle being θn < π/2 and θn ≥ π/2, respectively. tations {ê n } as emerging quantities from the local spinpolarization of the magnetic moment density. This establishes a framework to model transverse magnetic fluctuations at different temperatures, at and below the paramagnetic state, known as the disordered local moment (DLM) theory [21]. Our approach is based on this DLM theory and allows its application generally with any ab initio machinery able to employ magnetic constrains, as for example developed in Refs. [27,[36][37][38]. Hence, we open the opportunity to use DFT codes based on a plane wave basis instead of the Korringa-Kohn-Rostoker (KKR) formulation of DFT, and its bond to the singlesite coherent potential approximation [39,40] to carry out the ensemble averages.
The central component of DLM theory is an adiabatic approximation for the orientations of the local magnetic moments, {ê n }. These orientations are thus considered to evolve very slowly in comparison with the fast electronic motions forming the underlying many-electron in-teracting subsystem. Such an assumption means that thermally-induced transverse magnetic excitations can be obtained from DFT calculations constrained to comply with different magnetic configurations {ê n }, i.e., where µ(r; {ê n }) is the magnetic moment density, and Ω n is the volume defining the region at site n in which a local moment of size µ n emerges. Formally, a self-consistent calculation of the DFTbased total energy for given magnetic configurations {ê n } can be used to obtain the magnetic energy, E int ({ê n }) [21], introduced in Eq. (2). Such a calculation should be available by, for example, applying constraining single-site magnetic fields that establish {ê n } as minima of the system. We consider that the magnitudes {µ n } evolve in a much shorter timescale than the local moment orientations so that {µ n } instantaneously adapt to a given configuration {ê n }. We thus obtain {µ n } selfconsistently with respect to a given orientational configuration {ê n }. However, their thermal fluctuations could be described by performing additional averages [24].
E int ({ê n }) provides the corresponding partition function where β = 1 kBT , k B and T being the Boltzmann constant and the temperature, and is the magnetic energy also including a coupling of the local magnetic moments with an external magnetic field H. Note that {ê n } are treated as classical quantities and so Z and ensemble averages are obtained by performing integrals over all solid angle space, dê n = dθ n dφ n sin θ n . The probability of the system being in a magnetic state with local moment orientations {ê n } is, therefore, which gives the exact magnetic Gibbs free energy of the system by carrying out ensemble averages over all appropriately weighted magnetic configurations, The complex many-electron origins of E int ({ê n }) in metallic systems makes its dependence on the orientational magnetic configurations very complicated. Consequently, E int ({ê n }) can contain pairwise as well as higher order magnetic interaction terms, where E int,0 is a reference energy. In principle, the full complexity and myriad of magnetic interactions in E int ({ê n }) could be computed using DFT calculations if given enough computer time. However, the necessary computational costs are very far from available highperformance computational resources. Alternatively, we use a trial Hamiltonian generally given as n,n ,n (ê n ,ê n ,ê n ) + · · · , where {h n } are single-site internal magnetic fields, also called Weiss fields, forming the simplest mean-field theory for the interactions. In this work we only consider the single-site terms [see Eq. (14)], but by including higher order terms, {H (2) n,n (ê n ,ê n ), . . . }, our treatment can be systematically improved to account for different types of nonlocal magnetic correlations [21]. In section II B 1 we invoke the Peierls-Feynman inequality to find an upperbound of G [41,42], in analogy to the Gibbs-Bogoliubov inequality typically applied for the free energy of an ionic Hamiltonian [43,44]. This establishes a theory to calculate the Gibbs free energy and U mag = E int tr , which contains the effect of all magnetic interactions, with respect to the associated trial probability distribution, where Z tr is the corresponding partition function. The consideration of Eq. (10) reduces the thermally induced magnetic configurations to be studied to an affordable and meaningful magnetic phase space. A central task of our approach is to provide the internal magnetic energy as a function of the local magnetic order parameters of the system for this probability distribution, i.e. U mag ({m n }). · · · tr indicates an average over {ê n } with respect to P tr ({ê n }). As will be explained in section II B 1 and as illustrated in figure 1, in this work the trial orientational probability distribution is prescribed by the values of the order parameters, P tr ({ê n }; {m n }).
In Table I we advance and summarize the central steps, and their relation to the different sections and equations in the paper, for the computation of U mag ({m n }). These steps can be repeated for different lattice attributes to also obtain the dependence on the crystal structure and the consequent magnetostructural coupling. Once U mag and its dependencies have been computed, the Gibbs free energy is directly obtainable from Eq. (1). parameters are identical at all sites. Different values spanning from mn = 0 (paramagnetic) to mn = 1 (fully ordered ferromagnetic) can then be chosen for a given spin-polarization axis. Table I. Central steps to obtain the dependence of the internal magnetic energy on the local order parameters.
B. Computation of Umag = E int tr using magnetically constrained supercell calculations U mag = E int tr is obtained by carrying out an average of the magnetic energy over {ê n } with respect to a trial probability distribution P tr ({ê n }) through Eq. (2). Formally, a Monte Carlo integration can be employed, where V Ω = dê n = 4π is the solid angle volume, {ê n } rand i is a set of fully random local moment orientations, and N MC is the number of sets generated. E int ({ê n } rand i ) can be calculated using constraining magnetic fields in DFT calculations for supercells with fully random local moments. The number of atoms within the supercell, N sc , must be large enough to reduce boxsize effects satisfactorily. Once N MC magnetically constrained supercell calculations have been performed, Eq. (12) provides an approximated value of E int tr for an arbitrary shape of P tr ({ê n }). Note that different P tr ({ê n }) correspond to different states of magnetic thermal disorder, see Fig. 1. However, the minimum value of N MC to obtain accurately enough results increases for probability distributions P tr ({ê n }) that describe magnetic configurations away from fully random local moment orientations. It is to be presumed, therefore, that N MC is prohibitively large.
To circumvent this computational constrain, we approximate Eq. (2) instead by where now N MC magnetically constrained DFT calculations are performed for sets of local moment orientations that, instead of being fully random, directly obey the probability distribution We have found (see sections II F and III) that computationally affordable values of N MC can be achieved using Eq. (13) because the constrained magnetic configurations here already follow P tr ({ê n }).
The downside, however, is that a whole set of N MC configurations needs to be computed for each probability distribution that one aims to study, for example for each panel in Fig. 1. In the following section we show the case of the simplest single-site mean-field theory.

Peierls-Feynman inequality
Albeit our supercell approach to compute Eq. (13) is in principle implementable for a general expression of Eq. (9), we choose a trial Hamiltonian that contains only single-site magnetic fields, or Weiss fields, as originally used by Győrffy et al. [21], This H tr describes local magnetic moments whose orientations are sustained by mean fields {h n }, which can contain the effect of an external magnetic field if present. The partition function associated with H tr is and so the corresponding probability for a magnetic configuration {ê n } becomes where h n is the magnitude of h n . Owing to the singlesite nature of H tr , Eq. (16) defines single-site probability distributions for each local moment orientation, The Peierls-Feynman inequality [41,42] provides an upper bound of G given in Eq. (7), which we refer to as G u , where G tr = −k B T ln Z tr is the associated Gibbs free energy of the trial Hamiltonian, and we recall that · · · tr indicates an average over {ê n } with respect to P tr ({ê n }). From Eq. (5) and Eq. (18) one can finally write The previous result becomes Eq. (1) after applying a Legendre transform accounting for the effect of mechanical stress. It also introduces the magnetic local order parameters as the averages of the local moment orientations, as previously given in Eq. (11), as well as an expression for the magnetic entropy associated with the orientational configurations {ê n }, where Note that {P n (ê n )} is prescribed by {m n } since both of them are univocally given by {βh n }. This mean-field theory produces a typical classical behavior of m n following the Langevin function, as shown in Eq. (20) and plotted in the inset of Fig. 2. At zero temperature the local order parameter, which is always parallel toĥ n , has a magnitude equal to one, m n =ĥ n . This means thatê n does not thermally fluctuate. Note that in this limit, T → 0, {βh n } → ∞ and {h n } has a finite value. On the other hand, at high enough temperatures the thermal fluctuations can become sufficiently large to establish local moments fluctuating randomly over all possible orientations. This corresponds to a fully disordered, paramagnetic, state described by m n = 0 ({h n = 0}) and being present above Curie or Néel transition temperatures. We highlight that these transition temperatures as well as the dependence of the order parameters on temperature, {m n (T )}, are obtained by minimizing the free energy at different values of T , as explained in section II D. The snapshots shown in Fig. 1 of different magnetically constrained configurations correspond to different probabilities {P n (ê n )} (or {m n }) for a thermally fluctuating ferromagnetic state of bcc Fe.

Magnetically constrained supercell configurations
As indicated above, the probability distribution at each magnetic site n follows Eq. (17) independently owing to the single-site nature of Eq. (14), where 0 ≤ θ n ≤ π is the angle betweenê n and h n (or the magnetic local order parameter m n ). In Fig.  2 we show the behavior of P n (ê n ) against θ n for different values of βh n . While βh n = 0 describes a magnetic site that is fully disordered (m n = 0), βh n → ∞ corresponds to a non-fluctuating local moment orientation e n =ĥ n (m n = 1). We use a standard generator of uniformly distributed random numbers to obtain different magnetic configurations, which are the inputs in our magnetically constrained supercell calculations. Rotating theẑ-axis such that it is aligned with h n allows one to writeê n = (sin θ n cos ϕ n , sin θ n sin ϕ n , cos θ n ), where 0 ≤ ϕ n ≤ 2π is an in-plane angle. Since the probability integral must be normalized, we can obtain θ n by considering where x θ n is a random number ({0 ≤ x θ n ≤ 1}) and the 2π factor comes from the in-plane integration. Solving Eq. (24) gives (25) On the other hand, Eq. (23) does not depend on ϕ n , which means that ϕ n is fully random and can be obtained directly using another random number The local moment orientations for each magnetic configuration {ê n } {Pn(ên)} i are thus obtained from N pairs of random numbers, {x θ n , x ϕ n }, using Eq. (25) and Eq. (26) for given values of {βh n } (or {m n }) defining the probability distribution. Here N is the number of magnetic sites in the supercell. This is how the orientations shown in Fig. 1 have been obtained.

C. Hierarchy of local moment correlation functions and their magnetoelastic coupling
Since in our mean-field theory P tr ({ê n }) is uniquely prescribed by {βh n }, which in turn unequivocally map to {m n }, the average of E int has a direct dependence on {m n }. From Eq. (2) and Eq. (8) we, therefore, have where {U (2) ij , U (4) ijkl , . . . } form a hierarchy of second and higher order local moment correlations that are characteristic of the magnetic material investigated. Here h.o. stands for higher order terms. These local moment correlations describe expansion coefficients of the internal magnetic energy in terms of the magnetic local order parameters, in the spirit of a Ginzburg-Landau theory. As mentioned above, the dependence of E int on {m n } is performed through the computation of Eq. (13) using magnetically constrained supercells. We realize this calculation for a reduced phase space of {m n } containing the magnetic states of interest only, e.g. the ferromagnetic state in bcc Fe and a triangular antiferromagnetic state in section III. See Table I.
E int tr can be calculated for different lattice parameters and crystal structures and so its dependence on deformation can be also obtained. An expression of the Gibbs free energy in Eq. (19) that accounts for the effect of an applied hydrostatic pressure p is given by the following Legendre transform Here ij , U ijkl , · · · } inside E int tr depend on the volume V . Considering a parabolic dependence on V for the internal magnetic energy in the paramagnetic limit we can write where V PM is the volume that minimizes the free energy in this paramagnetic limit, ω = V −VPM VPM is the relative volume change, and γ is the inverse of the compressibility (bulk modulus) describing the energy cost caused by the application of p in the paramagnetic state. γ is computed, therefore, by calculating the change of the internal magnetic energy U mag = E int tr at different volumes in supercells containing fully random local moment orientations. Note that Eq. (28) can be formally generalized to account for more complicated deformations caused by the application of different types of stresses, as given by Eq. (1).
D. Minimization of the Gibbs free energy: First-order magnetic phase transitions generated by electronic and/or magnetoelastic mechanisms Once G tot has been obtained for a given magnetic material, it can be minimized with respect to {m n } for different values of T , H, and the crystal structure. This enables the ab initio study of temperature-dependent magnetic anisotropy [19,23], magnetic phase transitions and diagrams of different sort [22,[45][46][47][48], and (multi-)caloric effects for solid-state refrigeration [33,49,50]. Note that the second term in the right hand side of Eq. (19), i.e. the magnetic entropy, follows analytical expressions that are directly provided by the chosen trial Hamiltonian [see Eq. (21)]. U mag = E int tr is in consequence the only part of G tot that is material-dependent, obtained from first-principles through the computation of Eq. (13).
Higher than quadratic free energy coefficients in Eq. (27) fundamentally generate discontinuous (first-order) magnetic phase transitions from the paramagnetic state. This is a well known result formerly studied in MnAs, for which the emergence of fourth-order coefficients is driven by a magnetovolume coupling [51] [we also show this here through Eq. (34)]. However, our Gibbs free energy can also depend on quartic and higher order local moment correlation functions ({U (4) ijkl (ω), . . . }). These correlations are higher order free energy terms with a purely electronic origin, evidenced by the fact that they can emerge for calculations where the crystal structure is fixed. They describe the feedback to the magnetic interactions from the response of the underlying electronic structure to different states of thermally fluctuating magnetic order {m n } [30]. We have already shown how these multisite interactions generate first-order magnetic phase transitions experimentally observed in several materials [33,47,50]. An itinerant electron metamagnetism producing first-order transitions was already described by Wolfarth and Rhodes [52], and examples have been reported and explained by Fujita et al. for La(Fe x Si 1−x ) 13 compounds [53,54].
To exemplify how a first-order transition can be produced by both purely electronic and/or magnetoelastic effects, we consider materials that are composed of equivalent magnetic lattice sites. This implies that {m n = mm n }, i.e. the size of the local order parameters, m, is site-independent. In this case and in the absence of ex-ternal stimuli (H = 0, p = 0), the Gibbs free energy in Eq. (28) becomes where . . .
play the role of effective coefficients compactly containing the overall effect of all the local moment correlations, θ ij being the relative angle between m i and m j . We now assume that the magnetovolume coupling is described by a linear dependence for the second and fourth order terms, but that higher order coefficients are constant, , where n c is the order of the interaction. This is adopted here for illustrative purposes, but our ab initio results presented in sections II F and III are obtained without making any assumption in this regard. Minimizing Eq. (30) with respect to ω then gives which introduced back into Eq. (30) finally provides where N is the number of magnetic atoms. To obtain the previous equation we have used the following expansion of the magnetic entropy  by lowering the temperature, occurs when ∂ 2 Gtot ∂m 2 | m=0 = 0, which gives and that such a transition becomes discontinuous when Eq. (37) means that a first-order character arises when the overall fourth order coefficient at T = T tr,sec is negative, giving rise to the following condition which we plot in Fig. 3. Eq. (38) shows that a magnetovolume coupling α (2) always contributes to enhance the first-order character of the magnetic phase transition regardless of its sign. On the other hand, U 0 has to be positive to do so, and at least 30% of the magnitude of U (2) 0 if the discontinuous character has a purely electronic origin. Both interactions are possible mechanisms driving magnetic phase transitions of different types and between different magnetic states. In section III we study the Mn-based antiperovskite materials class, which exhibits a first-order phase transition to a triangular antiferromagnetic state whose origin is mainly electronic or magnetoelastic depending on the chemical composition of the material.

E. Computational details
All our fully non-collinear DFT calculations have been performed using the Vienna Ab Initio Simulation package (VASP) [55][56][57][58]. VASP already includes a method to magnetically constrain supercells to different local moment orientations {ê n }, developed by Ma and Dudarev [36]. Albeit we use such an implementation in this work, our approach can be applied together with any other DFT code suitable to perform magnetic constraints. Spin-orbit effects have not been considered. The parameter λ (referred to as a Lagrange multiplier in Ref. [36]) describing the energy penalty term E p to the total energy functional has been increased from low to high values as λ = 1, 3, 5, 7, 10, 15, 20, 25, 30 eV to preserve numerical stability within the DFT calculation. We have found that the final value of λ = 30 eV provides a satisfactory energy convergence with respect to the energy scale of the disordered and partially disordered magnetic states studied. We have constrained the orientation of the local moments only so their sizes have been allowed to relax within the self-consistent calculation for each magnetic configuration {ê n }, which is a fundamental feature of the time-scale separation of our theory (see section II A). The Wigner-Seitz radius needed as input for magnetic constraints was 1.323Å, 0.741Å, 1.217Å, and 1.286Å for Mn, N, Ga, and Ni, respectively.
We have employed the projector augmented-wave (PAW) method [59] implemented in VASP within the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation [60] to describe exchange and correlation effects. The Brillouin zone has been sampled using a 2 × 2 × 2 Monkhorst-Pack k-grid for supercells containing 432 atoms for bcc Fe (6 × 6 × 6 supercell of a cell of two atoms) and 320 atoms for Mn 3 AN (A = Ga, Ni) (4 × 4 × 4 supercell of a cell of five atoms). The cutoff energy has been 400 eV and 500 eV, respectively. Results for these materials are shown in sections II F and III. In appendix A we show results obtained for smaller supercells. The width of the smearing determining the partial orbital occupancies has been 0.1 eV. Finally, all the performed magnetically constrained DFT calculations have been supported, accelerated, and managed by Pyiron, an integrated development environment [61].

Dependence on ferromagnetic disorder
To illustrate how the magnetic Gibbs free energy of a material is obtained and minimized using our supercellbased ab initio approach, first we present its application to the referent bcc Fe. Focusing on the study of the ferromagnetic state, it follows that all the lattice sites of bcc Fe are equivalent. This means that in our mean field theory all the local magnetic moment orientations {ê n } fluctuate with the same single-site probability dis-tribution. The thermal fluctuations of a ferromagnetic state pointing along theẑ-direction are, therefore, described by {βh n = βhẑ}, or by {m n = mẑ}. Note that m = −1 βh + coth(βh) from Eq. (20). The values of the single-site magnetic fields and local order parameters are consequently the same for all sites.
We first compute the dependence of the internal magnetic energy on m as summarized in Table I. Panels (a), (b), and (c) in Fig. 4 show the total energies of 10 different orientationally constrained snapshots in 6 × 6 × 6 supercells of 2-atom cells of bcc Fe obtained for m = 0, m = 0.50, and m = 0.90, respectively (red dotted lines). Pictures of magnetic configurations for these values of m are represented in Fig. 1. The lattice parameter used is a = 2.83Å, which is the value that minimizes the energy in the ferromagnetic state within our DFT approximations [62]. We have found that the average of the total energy over N MC = 10 different magnetic configurations for each value of m suffices to compute U mag = E int tr using Eq. (13) within an error of only a small number of meV per atom (continuous blue line). The energy scale of the magnetic thermal fluctuations for the ferromagnetic state of bcc Fe is such that E int tr (m = 0) − E int tr (m = 1) ≈ 200meV/atom. We, therefore, can conclude that the application of Eq. (13) provides a satisfactory accuracy to obtain the dependence of the internal magnetic energy on sizable changes of m. However, the description of quantities with a smaller energy scale, such as the magnetocrystalline anisotropy, would require one to increase N MC and/or the size of the supercell.
The magnetic disorder characterizing the highly fluctuating paramagnetic state (m = 0) covers a wide range of magnetic excitations, which directly impacts the underlying electronic structure. Disorder typically broadens the spectral density function of electronic states in perfect periodic metals [63][64][65]. Consequently, a broadening of the total density of states (DOS) is also expected, as shown in Fig. 4(d). Indeed, decreasing the amount of magnetic disorder by increasing the value of m directly reduces the number of magnetic excitations thus sharpening the peaks of the DOS in Fig. (4)(e,f).
The time-scale separation between slowly evolving local moment orientations {ê n } and a rapidly adapting electronic structure is a central tenet of our DLM theory. Within this framework we allow the magnitudes of the local moments {µ n } to relax for a given magnetic configuration of the orientations within the self-consistent DFT cycle. In other words, we assume that these longitudinal fluctuations are much faster than the transverse ones and that their energy scale is much larger. In Fig. 4(g,h,i)   ments surrounding a given site n. The mean value of the magnetic moment magnitude, µ n , decreases slightly by lowering the order parameter from m = 0.50 to m = 0. Interestingly, we have found a non-monotonous behavior around m = 0.80 at which a maximum of µ n occurs. This is shown in Fig. 5(a), where we plot µ n against m for different volumes. The vertical bars indicate the standard deviation, which becomes larger by decreasing m in accordance to Fig. 4(g,h,i). To obtain the temperature dependence, which is encoded in m(T ), we need to minimize the free energy with respect to m at different values of T , which we address in the following.

Minimization of the free energy at different temperatures
The construction of the free energy, G tot , requires one to firstly obtain U mag = E int tr for different values of m, which must follow a polynomial dependence of even terms [see Eq. (27)]. In Fig. 5(b) we show E int tr against m computed for different volumes V , which we can express in terms of ω (see section II C). We have found that the ab initio data obtained can be described very well by  Fig. 5(c,d,e) with data points. Fig. 5(c) demonstrates that the internal magnetic energy in the paramagnetic state, U (0) , follows a pronounced parabolic behavior with a minimum at V PM , as stated in Eq. (29). On the other hand, U (2) and U (4) show a linear dependence, following the behavior in Eq. (32).
We can perform another regression, now of {U (0) , U (2) , U (4) } with respect to ω, to obtain {V PM , γ, U (2) 0 , α (2) , U (4) 0 , α (4) } using Eqs. (29) and (32). We show the results obtained in Table II. The bulk modulus in the paramagnetic state that we have com-   Table II. Coefficients describing the internal magnetic energy, the bulk modulus in the paramagnetic limit (γ), and the magnetovolume coupling. These are obtained by performing a least-squares regression of U (0) (V ) = U (0) (VPM) + 1 2 VPMγω 2 and of Eq. (32). The value of VPM given corresponds to a unit cell containing a single atom. Error estimations associated to the regression are also given. Figure 6. Dependence on temperature of (a) the total magnetization, M , the mean value of the magnitude of the local moments, µn , and of (b) the relative volume change, ω, obtained for bcc Fe.  Table II is already given in units of energy per atom). The computed value of T c is somewhat larger than the experimental one [67], but in satisfactory agreement considering the mean-field nature of our theory. A similar value has also been obtained using a mean-field and random phase approximation [68]. We also note that the leading magnetovolume coefficient, α (2) , is fairly small, which also agrees with experiment [67].
Introducing the values obtained for {V PM , γ, U (2) 0 , α (2) , U (4) 0 , α (4) } into Eq. (39) allows one to construct G tot [Eq. (30)], whose analytical minimization with respect to ω is given by Eq. (33). We can, therefore, compute the dependence on temperature of the magnetization and other properties by numerically finding the value of βh that minimizes Eq. (30) at different values of T , which is equivalent to perform a minimization with respect to m. Fig. 6 shows the results obtained for µ n (T ) and ω(T ), as well as for the total magnetization per atom, M (T ) = m(T ) · µ n (T ). Indeed, the paramagnetic-ferromagnetic phase transition is continuous and occurs at T c = 1590K. Albeit µ n (T ) presents a non-monotonic dependence on T , it is negligible and so the total magnetization behaves monotonously because its major contribution comes from m(T ). A prediction of our calculations is that α (2) is positive, which directly implies that the volume increases by lowering the temperature through T c , the so-called negative volume expansion (NVE). However, the computed effect is very small in comparison with the intrinsic lattice thermal expansion measured experimentally [67]. On the other hand, α (4) is negative and sufficiently large to cause a decrease of ω at lower temperatures when the fourth order free energy terms become important.

III. AB INITIO ORIGIN OF FIRST-ORDER MAGNETIC PHASE TRANSITIONS IN ANTIPEROVSKITE NITRIDE MATERIALS
Mn-based antiperovskite nitride systems Mn 3 AN, where A can be one or a solution among different transition metals and semiconductor elements, form a magnetic materials class that has been largely studied over the last decades owing to their interesting magnetic properties [69][70][71][72][73][74][75][76][77][78][79]. The corresponding cubic crystal structure belongs to the perovskite space group 221 in which A, N, and Mn atoms sit on corners, centers, and face centers of the unit cell, respectively. In these materials a magnetic phase transition from the paramagnetic state to a triangular antiferromagnetic state occurs at different temperatures and emerges from the geometrical frustration of antiferromagnetic interactions between nearest neighbor Mn atoms. The transition is accompanied by a substantial NVE [31,71,72] and has a first-order (discontinuous) character, which provides a giant barocaloric effect with great potential in the field of caloric refrigeration [32,34,80].
We have recently made an analysis of Mn 3 AN materials based on the exploration of free energy coefficients in Eq. (30) that reproduce experimental data describing the change of the magnetic phase transition under applied hydrostatic pressure. Such an experimentallyguided study has shown that the first-order character of Mn 3 AN arises from the combined effect of the magnetovolume coupling and multisite magnetic interactions, and that the choice of atom A selects which one of these sources is the predominant mechanism [33]. Here we use our supercell DLM approach to provide a full fristprinciples demonstration of this phenomenon.
We focus on the study of Mn 3 GaN and Mn 3 NiN, whose first-order character should be mainly driven by a magnetovolume coupling and multisite interactions, respectively [33]. The triangular antiferromagnetic state that stabilizes in these materials is a non-modulated (q = 0) magnetic structure. Here the three Mn magnetic local moments lie within the same plane and form relative angles of 120 degrees at T = 0K when magnetic anisotropy is not considered. If the cubic crystal structure is not distorted, such a triangular state is fully compensated, i.e. there is no overall magnetization. This symmetry also means that the thermal fluctuations of the local moment orientations at the three sites can be described by triangular local order parameters that have the same mag- nitude m tri . Since we do not consider spin-orbit effects we simply choose an arbitrary plane and apply Eq. (25) and Eq. (26) to find the polar and azimuthal angles for each moment in a supercell. The local axis frames are rotated accordingly to the triangular antiferromagnetic state such that {h n } (or {m n }) correctly form angles of 120 degrees between the three non-equivalent magnetic sub-lattices.
In Fig. 7 we show the dependence of µ n and the internal magnetic energy on m tri for different volumes. One can see in panels (a) and (b) that the local moment sizes in Mn atoms are statistically smaller and have a larger standard deviation in Mn 3 GaN than those computed in Mn 3 NiN. However, robust local moments emerge in both cases for larger volumes. We note that some of our calculations presented difficulties to converge for Mn 3 GaN at a = 3.78Å and m < 0.2, for which the local moment magnitudes are relatively small. Contrary to the ferromagnetic state of bcc Fe in section II F, here µ n follows a full monotonic behavior with m at all studied volumes. Very different outcomes have been obtained regarding the calculation of U mag = E int tr : Its curvature found in Mn 3 NiN is fairly insensitive to the volume, which is in sharp contrast to Mn 3 GaN. This means that Mn 3 GaN and Mn 3 NiN contain comparatively large and small magnetovolume coupling, respectively, as observed experimentally [33].
The dependence of E int tr on the order parameter can be described rather well by an expansion up to a fourth order, We can carry out again two consecutive least-squares regressions. The first for Eq. (40) to obtain the corresponding internal energy coefficients, and the second to fit these against the different volumes studied.
We show the results obtained from the first regression in Fig. 8 as data points. An inspection of the figure reveals that in both materials U (0) and U (4) follow parabolic behaviors. However, we have found a distinctive feature regarding U (2) , which is fairly constant in Mn 3 NiN while in Mn 3 GaN it shows a substantial linear dependence. We thus write which contains an additional term for a quadratic magnetovolume coupling described by β (4) , in comparison 0 , α (4) , β (4) }), obtained by performing a regression of Eq. (41). Errors associated to the regression are provided.
to Eq. (32). In Table III we show the coefficients of Eq. (41) obtained by performing the second regression with respect to ω. We use continuous lines in Fig. 8 to plot the corresponding results. The bulk moduli that we have computed in the paramagnetic limit are 84 GPa and 117 GPa for Mn 3 GaN and Mn 3 NiN, respectively. These values are somewhat smaller compared with DFT calculations made at zero temperature [73], i.e. magnetic thermal fluctuations reduce the bulk modulus. Most importantly, the leading second order magnetovolume term, α (2) , is very large for Mn 3 GaN but negligible for Mn 3 NiN, as already seen in Fig. 7(c,d) and reported in experiment [33].
The free energy per formula unit of both materials is where N sc is the number of Mn 3 AN cells with volume V PM (containing five atoms) forming the supercell, and so here {U (0) , U (2) , U (4) } are also given per formula unit. The factor of 3 accompanying the entropy term comes from the fact that one crystallographic unit cell contains three Mn atoms. Minimizing Eq. (42) provides the total magnetization, µ n , and ω as functions of T . An analytical expression for the latter is directly given by In Fig. 9 we show the results obtained after performing this minimization. The figure demonstrates that firstorder magnetic phase transitions from the paramagnetic state to the triangular antiferromagnetic state are directly described by our theory and that a NVE at the transition is correctly captured. The application of Eq. (38), which requires the calculation of α (2) 2 /2V PM γ (given in Table III), shows that the origin of the firstorder character in Mn 3 NiN is the purely electronic term U (4) , as found experimentally [33]. On the other hand, in Mn 3 GaN both U (4) and the magnetovolume coupling α (2) are important factors behind the discontinuous character of the transition (see also Fig. 3).
Values computed for both the transition temperature, T tr , and the spontaneous volume change at the transition, ∆ω, provided in Table IV, agree well with experiment for Mn 3 NiN. However, we underestimate the former and overestimate the latter for Mn 3 GaN [33,78]. We point out that T tr strongly increases with U (2) 0 . Hence, according to the results given in Fig. 8(b) a lattice thermal expansion, not considered so far but always present in materials, should substantially enhance T tr in Mn 3 GaN while having a little impact in Mn 3 NiN. This is indeed what we have found when we model this effect by modifying the third term in the right hand side of Eq. (42) as where C TE > 0 drives an increment of ω when T increases. Note that now Eq. (43) becomes ω = ω| CTE=0 + C TE T . Dashed curves in Fig. 9 show how the magnetic properties change when C TE = 40 × 10 −4 meV/K and C TE = 35 × 10 −4 meV/K are used for Mn 3 GaN and Mn 3 NiN, respectively. These estimations are chosen to yield a lattice parameter that approximately matches the experimental values in the paramagnetic state [31]. We have found that such an effect enhances the transition temperature of Mn 3 GaN up to 315K. This value is now higher than the transition temperature obtained for Mn 3 NiN, which correctly captures the experimental trend (T tr,Mn3GaN > T tr,Mn3NiN , see Table IV). We conclude that accounting for the lattice thermal expansion is crucial for Mn 3 GaN owing to its very large magnetovolume coupling. A major part of our results above agree qualitatively, and sometimes quantitatively, with experiment. However, our theory seems to substantially overestimate U (4) in Mn 3 GaN [33]. Several factors can be behind this significant discrepancy. It is well known that Mn 3 AN materials usually present a deficiency of nitrogen occupation [31,72] that potentially impacts the magnetic properties. Furthermore, minor chemical or positional disorder at the Mn sites can greatly modify the geometrical frustration of the magnetic interactions. We also expect, therefore, that a magneto-phonon coupling is another important component to accurately describe the magnetism of Mn 3 AN. All these are aspects that we have not included in our calculations and that could explain the disagreement. Nevertheless, our approach distinguishes between electronic and magnetovolume mechanisms, describes how these depend on chemical composition, and correctly captures a major part of complicated features, such as the first-order character and the NVE.
As a final remark, we highlight that we have described the dependence of the internal magnetic energy in Eq. (40) by considering the lowest possible orders of the expansion coefficients. However, we observed that such a dependence can be also described numerically well by taking a sixth order coefficient, U (6) , instead of the fourth order, U (4) . Considering U (6) provides the same qualitative and quantitative behavior for the paramagnetic-triangular antiferromagnetic phase transition in Mn 3 GaN, with only tiny numerical differences. On the other hand, for Mn 3 NiN the transition becomes continuous and falls right below a tricritical point separating first-and second-order behaviors, i.e. it becomes nearly discontinuous.

IV. CONCLUSIONS AND OUTLOOK
The development of new ab initio theories describing the temperature evolution of magnetic materials is a major challenge greatly demanded by and strongly impacting the research of solid-state magnetism. It broadly includes the investigation of functional magnetic phase transitions between complex magnetic structures that are exploited in a wide range of technological applications. In this work we address this challenge by developing an approach that provides the ab initio magnetic Gibbs free energy of a material from magnetically constrained supercell calculations. Its basis is the description of the statistical mechanics of local magnetic moments, assumed to evolve very slowly following the disordered local moment (DLM) picture [21]. We compute the internal magnetic energy of the material by performing averages of a density functional theory-based first-principles magnetic energy over a large but affordable number of noncollinear local moment configurations.
We have applied our supercell approach to study the ferromagnetic state of bcc iron and the triangular antiferromagnetic state present in the geometrically frustrated antiperovskite systems Mn 3 AN (A = Ga, Ni). Our results are in good qualitative, and sometimes quantitative, agreement with experiment. Most importantly, we describe correctly the character of the magnetic phase transitions from the paramagnetic state, either continuous or discontinuous, and quantify its origin in terms of purely electronic and/or magnetostructural sources. We have found that the mechanism giving rise to the firstorder character of Mn 3 NiN arises purely from multisite magnetic interactions while this effect as well as a magnetovolume coupling play a major role in Mn 3 GaN. Potential explanations for disagreements with experiment have been discussed.
Magnetically constrained supercell calculations are the principal computational component of our approach. This enables the application of our DLM theory using density functional theory codes based on a plane-wave basis (VASP in this work), i.e. beyond the Korringa-Kohn-Rostoker formalism and the coherent potential approximation. Our approach is computationally expensive but it is already affordable by existing supercomputers. Furthermore, satisfactory qualitative results are obtained using relatively small supercells.
An important advantage of this DLM theory is that the trial, mean-field, Hamiltonian prescribing the local moment averages [see Eq. (9)] can be naturally extended beyond the simplest single-site Weiss field parameters [21]. In other words, magnetically constrained supercell calculations can be directly used to account for nonlocal magnetic correlations akin to the nonlocal coherent po-tential approximation [64,65,81]. Along these lines, spin-cluster expansions, which also follow an adiabatic approximation for the local moments, can be used to efficiently construct complex magnetic Hamiltonians [82]. Furthermore, magnetically constrained calculations directly output the magnetic torques added to sustain the transient orientational magnetic configurations. Our approach, therefore, also offers opportunities to be combined with spin-dynamics and related finite-temperature methods [83,84]. Finally, additional averages of the internal DFT energy can be performed over other degrees of freedom, such as the local moment magnitudes and the atom vibrations [85][86][87]. Pertinent time-scale separations for the lattice dynamics could then be considered [28,88] to account for a magneto-phonon coupling [89].  Fig. 10 shows the second and fourth order internal energy coefficients computed for ferromagnetic bcc Fe and for the triangular antiferromagnetic state of Mn 3 NiN with smaller supercells. These calculations have been carried out for lattice parameters equal to a = 2.825Å and a = 3.86Å, respectively. In Table V we provide the number of supercell snapshots used to carry out the average of the magnetic energy for every value of the magnetic order parameter. This number becomes larger for smaller supercells if similar statistical accuracy for the energy average is required. The Table also shows the number of k-points used within the Monkhorst-Pack grid sampling.
The same qualitative results are obtained even after reducing the supercell size to contain only a few tenths of atoms. For example, the ratio U (4) /U (2) remains approximately zero and close to 1 for bcc Fe and Mn 3 NiN, respectively. These values are below and above the critical condition U (4) /U (2) = 3 10 , which directly implies that the corresponding character of the magnetic phase transition from the paramagnetic state is second-order and first-order, respectively [see Eq. (38) and Fig. 3]. On the other hand, we observe a non-negligible quantitative change of U (2) for bcc Fe, i.e. there is a dependence of the computed transition temperature, T tr , on the size of the supercell [see Eq. (36)]. In this case, T tr decreases when the supercell becomes larger. We conclude that qualitatively correct results can be already obtained for relatively small supercells, but that very accurate quanti- Figure 10. Dependence on the supercell size of the second and fourth order internal energy coefficients, {U (2) , U (4) }, calculated for the ferromagnetic state of bcc Fe and the triangular antiferromagnetic state of Mn3NiN. The lattice parameters are a = 2.825Å and a = 3.86Å, respectively. The critical line (dashed) separating regions between second-( U (4) U (2) ≤ 3 10 ) and first-( U (4) U (2) > 3 10 ) order magnetic phase transitions from the paramagnetic state is indicated in the lower panels.
tative calculations require a large number of atoms. However, careful reassessment in this regard should be made for the evaluation of other magnetic materials.  Fig. 10. The higher the number of atoms in the supercell, the higher the number of necessary snapshots to achieve satisfactory accuracy for the average.