Calphad Modeling of LRO and SRO Using ab initio Data

: Results from DFT calculations are in many cases equivalent to experimental data. They describe a set of properties of a phase at a well-deﬁned composition and temperature, T , most often at 0 K. In order to be practically useful in materials design, such data must be ﬁtted to a thermodynamic model for the phase to allow interpolations and extrapolations. The intention of this paper is to give a summary of the state of the art by using the Calphad technique to model thermodynamic properties and calculate phase diagrams, including some models that should be avoided. Calphad models can decribe long range ordering (LRO) using sublattices and there are model parameters that can approximate short range ordering (SRO) within the experimental uncertainty. In addition to the DFT data, there is a need for experimental data, in particular, for the phase diagram, to determine the model parameters. Very small differences in Gibbs energy of the phases, far smaller than the uncertainties in the DFT calculations, determine the set of stable phases at varying composition and T . Thus, adjustment of the DFT results is often needed in order to obtain the correct set of stable phases.


Introduction
In materials physics the scientists performing density functional theory (DFT) calculations have a preference for the cluster variation method (CVM) or Monte Carlo (MC) methods to calculate the equilibria in their systems because these can describe both long range ordering (LRO) and short range ordering (SRO) in the crystalline phases. However, both CVM, developed by Kikuchi [1], and MC, are computationally heavy for multi-component systems, and this limits their use for real alloys. Thus there is interest in using Calphad type models for equilibrium and phase diagram calculations and this paper gives an introduction to how these models handle LRO and SRO.
The Calphad models based on the compound energy formalism (CEF) [2], as described in this paper, are computationally highly efficient and make it possible to calculate multi-component equilibria and phase diagrams for systems with many components and phases. Such phases can feature four or more sublattices to describe complex states of LRO and SRO. CEF also includes an approximate way to handle SRO [3]. The use of data from DFT calculations in Calphad type models has been discussed in many papers [4][5][6][7][8] and an example of how to use DFT data in CEF models will be given below together with some advice for selecting the appropriate models for different phases. Recently, a further modification of CEF has been proposed by Dupin et al. [9] to improve both extrapolations and calculation speed.
where a s is the site ratio, y si is a constituent fraction and b iA is a stoichiometric factor. Phase has a defined lattice or is a gas, liquid or amorphous phase and has a thermodynamic model. Reciprocal model within CEF has two sublattices with two constituents in each.
Site ratio, denoted a s for sublattice s, is the number of sites in the sublattice in the unit cell.
Species is a molecule-like aggregate with one or more elements with fixed ratios-for example, H 2 O. The vacancy, denoted Va, is a special species used to represent a vacant site in a phase. A species may also have an electric charge and be called ion. SRO: Short range ordering is present when the local environment of a given species is not randomly surrounded by other species. However, the deviation from randomness vanishes as the distance increases. Stoichiometric factor, denoted b iA is the ratio of element A in species i. Sublattice: A specific set of sites in the unit cell of the lattice. Only sites which are, at least partially, occupied by atoms in a crystal structures are recognized as sublattices. In some cases a sublattice may combine sites which are crystallographically different.
A substitutional regular solution is a model where all constituents are assumed to mix randomly with energetic interaction parameters between sets of two or more constituents. The thermodynamic model describes the Gibbs energy of a phase as a function of T, P and its constitution or cluster probabilities.

Modeling Long and Short range Ordering
In the CVM, configurational entropy associated with SRO and LRO is treated approximately within a certain spatial range defined by the so-called maximal clusters. A well-known set of maximal clusters with a spatial range including both nearest and next-nearest neighbors are the octahedron-tetrahedron clusters for the fcc structure and its superstructures. The maximal clusters contain smaller subclusters. The configurational entropy is then evaluated as a sum of entropy contributions over all clusters from the maximal clusters down to the smallest subclusters. The complicating factor that occurs in this sum, is that smaller subcluster contributions are already included to various extents in the larger clusters. Therefore, a sophisticated weighting scheme through the so-called Kikuchi-Barker coefficients is required to assure that each entropy contributions enters in the sum with a proper weight. The selection of the maximal clusters is not trivial. There exist maximal clusters that do not yield properly weighted entropy sums; however, well-chosen maximal clusters can give remarkably accurate configurational entropy values. The CVM is not well suited for treating cases with many species, or when long-ranged interactions necessitate large maximal clusters. In both these cases the number of (sub)cluster configurations becomes overwhelmingly large and the numerical solution of the CVM equations impractical. For that reason the CEF is the method of choice within the Calphad approach.
In the CEF, a phase with several sublattices (sites in the unit cell) is thought of as a mixture of several endmembers that at least hypothetically can exist and which specifies a constituent on each sublattice. The configurational entropy of the phase is calculated by assuming that the constituents on each sublattice mix randomly weighted by the site ratio for each sublattice. This can be considered as the LRO contribution to the configurational entropy of the phase.
For a phase which always has LRO and would require complex clusters in the CVM, it is possible to use a CEF model and adjust some of the endmember energies to compensate for the lack of SRO. In some cases it may be necessary to choose a cell that is a multiple of the unit cell, requiring more endmember energies to be adjusted or to use a method discussed by Liu [24]. For very small unit cells, for example, in phases with order-disorder transitions such as A1/L1 2 or A2/B2, which can exist with considerable SRO contributions without any LRO (and thus without additional endmember parameters in a CEF model), an approximate SRO contribution can be added, as explained in Section 2.4. Clearly, the CEF model provides an easily solvable model that requires as input a few fraction variables pertaining to sublattices, and a limited number of endmember energies. The endmember energies can be computed with DFT methods, and in Section 5 the methodology to perform an assessment combining experimental and DFT data is discussed.

Using the Reciprocal Model in CEF to Model LRO and SRO
CEF and CVM are mean field models which both aim to describe the thermodynamic properties of crystalline phases. Nevertheless they are based on completely different concepts: CEF originates from a model describing the entropy of molten salts without any crystalline structure. The reciprocal model, which is the basic modeling tool within CEF, is based on studies of reactions in salt systems such as: NaCl + KBr <=> NaBr + KCl and Temkin [25] proposed in 1945 a model mixing cations and anions on two separate sets of sites, a two-sublattice model without any SRO. Hillert and Staffanson [26] adapted this model 1970 to describe interstitial solutions of C in steels and it was generalized in 1981 to multiple sublattices and constituents by Sundman and Ågren [27] as part of a long-term project to develop databases for steels, oxides and superalloys. In 1998 an approximate SRO contribution was added to CEF by comparing the configurational entropies for a four-sublattice CEF model and a tetrahedron CVM model for FCC in an assessment of Au-Cu by Sundman et al. [3]. In 2001, Hillert [2] summarized this development which has also resulted in the Thermo-Calc [28] software and databases.
All CEF models with sublattices are extensions of the reciprocal system with two sublattices and two constituents in each, denoted: where A and B are constituents on a sublattice with the site ratio a and C and D are constituents on the other sublattice with the site ratio c. Note that this model does not assume any geometric arrangement of the sites. The same species can be constituents in both sublattices and the site ratios are the only references to the crystal structure.
The reciprocal model is important to understand because phases with several sublattices and with three or more constituents in each will have many model parameters, particularly endmembers. To understand the relations between the endmembers it is useful to divide the phase into a sequence of reciprocal subsystems. See, for example, how to find the endmembers related to energy for a Frenkel defect in Section 4.3. We thus start by exploring the properties of the reciprocal model, and for further details of the CEF model, see the book by Hillert [29] or Lukas et al. [11].
The constitutional square of the reciprocal model is shown in Figure 1a with the four endmembers (A:C), (A:D), (B:C) and (B:D) in the corners and four interaction parameters, L AB:C , etc., along the edges. In the center there is a reciprocal interaction parameter, L A,B:C.D , representing interactions on both sublattices. In Figure 1b an example of the Gibbs energy surface for a reciprocal system is shown.  The Gibbs energy for a reciprocal model for the phase α is: where G M is the Gibbs energy per mole formula unit and T the absolute temperature. The terms srf G M (the surface of reference), cfg S M (configurational entropy) and E G M (excess Gibbs energy) are given in Equations (4)- (6). The subscript M in Equation (3) is used to indicate moles per formula units of the phase because a constituent may contain more than one atom or represent a vacancy, whereas m in G m means per mole of components. The phase superscript is omitted in the following unless needed.
srf G M = y 1,A y 2,C • G A:C + y 1,A y 2,D • G A:D + y 1,B y 2,C • G B:C + y 1,B y 2,D • G B:D (4) cfg S M = −R(a(y 1,A ln(y 1,A ) + y 1,B ln(y 1,B )) + c(y 2,C ln(y 2,C ) + y 2,D ln(y 2,D ))) (5) E G M = y 1,A y 1,B (y 2,C L A,B:C + y 2,D L A,B:D ) + y 2,C y 2,D (y 1,A L A:C,D + y 1,B L B:C,D ) + y 1,A y 1,B y 2,C y 2,D L A,B:C,D (6) where y si in the fraction of constituent i on sublattice s; • G A:C is the Gibbs energy of the endmember (A:C), where the "pre-superscript" • means the parameter is independent of the constitution. R is the gas constant and cfg S M the configurational entropy assuming ideal mixing on each sublattice. The contribution depends on the number of sites a and c on each sublattice. The interaction parameters L A,B:C , etc., give the interaction between two constituents in one sublattice with a specific constituent in the second sublattice. Finally L A,B:C,D is an interaction parameter for all four constituents, and all interaction parameters may depend on T.
The interaction parameters L A,B:C may depend on the fractions of the interacting constituents as a Redlich-Kister series: where ν is a power and n should normally not exceed three and each ν L A,B:C may depend linearly on T. The reciprocal parameter, L A,B:C,D , represents a simultaneous exchange energy between all four constituents and is frequently used to approximate the SRO contribution, as discussed in Section 2.4. Understanding the properties of the reciprocal system is essential for modeling any multi-sublattice phase and it occurs in many different phases; for example:

•
An interstitial solution of C in the FCC phase in a C-Fe-Ti system is modeled as (Fe,Ti) 1

The Endmember Concept and Metastability
An endmember parameter • G A:C represents the Gibbs energy of a compound with a given structure. These are sometimes stable compounds, but in a solution phase they most frequently are metastable at all T. For example, in a Laves_C14 phase modeled as (Fe, Mo) 2 (Fe, Mo) 1 only the endmember, Fe 2 Mo, represents a stable phase. Two of the other endmembers represent pure Fe and pure Mo as metastable in the Laves_C14 structure and the fourth endmember is a metastable antistructure Laves_C14 phase with the wrong stoichiometry.
The Laves_C14 phase is an example of a phase with very small solubility in the Fe-Mo system and it is frequently modeled using defects. The defects can be of different types, such as vacancy, anti-site atom or interstitial. In Calphad models the Laves phases are normally modeled with anti-site atoms as Laves phases occur in many other binary systems, including either Fe or Mo. In order to develop multi-component databases it is necessary to treat all phases with the same structure with the same model because the binary Laves phase often dissolves other elements. It may even occur inside a ternary system without being stable in any of the binaries. It is thus simpler to use the same type of defect in all cases, especially to reduce the number of endmembers. In order to model the Laves phase in a particular binary with no intention to extend to multi-component system, one should of course use the physically correct defects.
The enthalpy of formation of the metastable endmembers for many pure elements in several intermetallic phases has been calculated by Sluiter [30], and for compatibility between assessments it is important to use the same values for the same structures and elements in different assessments. These have almost the same relevance for multi-component modeling as the lattice stabilities of the pure elements discussed in Section 1.1.

The Reciprocal Energy and Miscibility Gap
The Gibbs energy difference between the endmembers on the diagonals in a reciprocal system is very important: This reciprocal energy, ∆G recip AB:CD , corresponds to the energy difference of the diagonals in the center of the reciprocal square, as shown in Figure 2 for two different cases of ∆G recip AB:CD . In Figure 2a the diagonals have the same energy in the middle but in Figure 2b there is a difference, and if this is sufficiently large there is a miscibility gap, even if all excess parameters are zero.  In many reciprocal systems the miscibility gap is real-for example, in the C-Fe-Ti system (see Section 4.1) there is one between austenite where C dissolves interstitially in FCC-Fe and the cubic TiC carbide, with considerable fraction of Va, both having the B1 structure. However, in other systems the miscibility gap can appear because a metastable endmember parameter is badly estimated. This can be problematic because it is not possible to suppress this miscibility gap using interaction parameters.
When modeling a single phase with the reciprocal model it is always possible to set three of the endmember energies to zero to plot a Gibbs energy surface, as in Figure 2b. If the Gibbs energies of the endmembers vary independently on T there is a risk of forming a miscibility gap when T varies. However, frequently one or more of the endmembers are metastable and can be estimated in a way that prevents a reciprocal miscibility gap. For example if three are known, the fourth, • G B:D , can be estimated by forcing ∆G recip AB:CD = 0: A reciprocal miscibility gap can usually be avoided in an assessment but may appear when combining two independent assessments of a phase modeled with sublattices to a ternary system because they share several endmembers. If experimentally there should be no miscibility gap, in this phase it can be quite difficult to suppress that without modifying the original assessments.
The reason to use the reciprocal model should be that A and B mix on one sublattice and C and D on another. We can compare this with a substitutional model mixing four species, AC, AD, BC and BD. In such a model there is no tendency to form a miscibility gap, but the crystallographic information that A and B mix on one sublattice and C and D on the other is lost.

The Reciprocal Model for Order-Disorder Transitions
A reciprocal model (A,B) a (A,B) a can be used for a phase which can be totally disordered, i.e., with y 1,A ≡ y 2,A in the disordered state. These endmember energies must be equal: where u AB as the bond energy and z is the number of bonds between the sublattices. Ansara et al. [31] used a reciprocal model with different numbers of sites on the sublattices in their first assessments of the L1 2 ordering in FCC. Such a model means some atoms have nearest neighbors in the same sublattice and require additional interaction parameters related to the bond energy in that sublattice. This model for FCC ordering is now discouraged, as discussed in Section 4.5.2, and a four-sublattice model where all nearest neighbors are on different sublattices is recommended. As already discussed, such a four-sublattice model can be considered as a superposition of several reciprocal models.
(a) y 1,B and y 2,B vs T at x B =0. 5 ( The reciprocal model by itself can be used for an order-disorder transition when all nearest neighbor atoms are on the opposite sublattice-for example, an A2/B2 ordering modeled without vacancies. Figure 3 has several diagrams for a bond energy parameter u AB = −300R. The diagrams show how the constituent fractions of the sublattices vary with T or composition, the resulting configurational heat capacity and the 2nd order transition line in the phase diagram. In the model u AB R = 600 K at the equi-atomic composition and in Figure 3d the configurational heat capacity at the equiatomic composition is plotted as a function of T. For T > T o/d it becomes zero as there is no SRO in this version of the reciprocal model; see Section 2.4. In Figure 3e the configurational heat capacity is plotted as a function of composition at T = 400 K in the ordered state and there is a peak in this heat capacity at the ideal ordering configuration which has been explained in [32,33].

An Approximation of the SRO Contribution
In an assessment of the Au-Cu system, Sundman et al. [3] derived a first order approximation of the SRO contribution to CEF by comparing the configurational entropy in a quasichemical model for BCC to that of a recipocal model. They could show that a reciprocal parameter, with a value related to the bond energy, gives a topologically correct phase diagram with separate maxima for ordered L1 2 and L1 0 also for a four-sublattice CEF model of Au-Cu.
The approximate SRO contribution for a reciprocal system (A,B)(A,B) was derived as: where ∆G recip AB:AB is the reciprocal Gibbs energy defined in Equation (8) when the constituents on the two sublattices are identical and z is the number of nearest neighbors. ∆G SRO M depend on the constitution in the same way as the reciprocal parameter, L A,B:A,B , and the approximate SRO contribution, u SRO AB , can be included in the reciprocal model parameter: In a BCC phase ∆G recip AB:AB = z 2 u AB for a model with one atom/(unit cell) and at the temperature for the equiatomic order-disorder transition, u AB R , and these values can be inserted into Equation (13) to obtain the SRO parameter at the equiatomic order-disorder T: Note that u AB must be negative to promote ordering and u SRO AB is always negative in Equation (13). The bond energy u AB will normally depend on T but the constant value evaluated at T o/d is a practical approximation at all T. There is a factor 1 T in Equation (13), which could give a reasonable decrease of the approximate SRO entropy at higher T, but such a factor would make the disordered state stable at low T. A constant u SRO AB does not contribute to the configurational heat capacity in the disordered state and will slightly overestimate the SRO at high T, but this can be compensated for by other model parameters.
The value of a reciprocal interaction parameter in an assessment will depend on many different kinds of experimental and theoretical data, and using a constant bond energy for SRO is a reasonable first approximation. In a CEF model, few parameters have a direct physical meanings but contain contributions from many different phenomena.

The Disordered Fraction Set
A phase α modeled with sublattices may have some properties which are independent of its constitution, i.e., y α si , but others which may depend on the composition, i.e., the mole fractions x α i , as discussed in [11]. The mole fractions can be calculated from the constituent fractions, y, using Equation (2). If the constituents are the same as the elements, and we have vacant sites, denoted Va, in some sublattices this becomes: Most implementations of the CEF model have introduced a "disordered fraction set," which grants a possibility of including a separate set of Gibbs energy parameters, G α/dis M (x), which are added to the Gibbs energy for the ordered part, G α/ord M (y) using the equation: where the set of mole fractions, x, are calculated from the set of constituent fractions y, using Equation (15). The configurational entropy of the disordered fraction set, cfg S α/dis M , is subtracted because the configurational entropy is included in G α/ord M using the constituent fractions. Further details can be found in [11].
The physical reason to add a disordered fraction set is that many phases with LRO have properties that are independent of the constitution and parameters, for these can be described using this set. However, the disordered fraction set also simplifies the construction of multi-component databases for phases with order-disorder transitions, such as BCC, FCC and HCP, which can exist both as completely disordered (with equal fractions on all sublattices) and as ordered. Parameters from binary and ternary assessments with only disordered BCC, FCC and HCP are simply added to the disordered fraction set when combined with systems assessed using the ordered model for these phases.
For historical reasons there is also a slightly different way to add the ordered and disordered parts: With these equations, introduced by Ansara et al. [34], it is possible to assess the disordered part independently from the ordered. It can be useful when assessing a system with a phase such as FCC in Au-Cu that exists both as ordered and disordered. The parameters in G α/dis M in Equation (17) are the complete description of the phase when it is disordered and the parameters in G α/ord M can be assessed independently because they have no influence when the phase is disordered; i.e., y = x.
In addition to the sublattices for order-disorder, there may also be a sublattice for interstitial constituents in both the ordered and disordered fraction set.

CEF Calculations for Prototype FCC Ordering
In 1938, Shockley [35] presented a phase diagram (as in Figure 4a) for fcc ordering in an FCC system using a four-sublattice model with an ideal configurational entropy. This diagram has the wrong topology, as the maxima for the ordered phases, L1 2 and L1 0 , are not separated and the reason for this is that the SRO contribution is ignored. In 1951 Kikuchi [1] derived CVM, a mean field model for LRO and SRO and when computers were becoming available some 20 year later, several calculations of topologically correct phase diagram for FCC with separated maxima as in Figure 4b using CVM models were published [36][37][38]. For binary systems with simple lattices these calculations are trivial today, but for more complex structures, CVM models which may require a large number of clusters to describe SRO even, in systems with few components, there is still an interest in simpler models to describe SRO. An example of a simpler model is Figure 4c, a prototype FCC ordering with correct topology calculated using a four-sublattice CEF model with same bond energy as in (a) and (b) and with an SRO contribution according to Equation (14).  (a) The CEF model without SRO is used (as by Shockley), in (b) a tetrahedron CVM model and in (c) a CEF model with SRO according to Equation (14). Note also that the T axis in (a) is double that in (b)-(d).
In (d) and (e) a configuration independent parameter has been added.
The phase diagrams in Figures 4c,d use the four-sublattice CEF model with the parameters in Equation (19) and in Figures 4d,e also the interaction parameter in the disordered fraction set, L A1 A,B , in Equation (20).
where u AB is the bond energy, R the gas constant and the • G ijkl are the energies for the endmember tetrahedra occupied with different sets of atoms. The · · · indicate all parameters with the same set of atoms permutated on the sublattices are equal. In a reciprocal parameter, L A,B:A,B: * : * , the "*" indicates that the parameter is independent of the constituents in the sublattice. The configuration independent parameter in Equation (20) can be used to fit experimental data with no or minor effect on the topology. The reciprocal parameters are an approximate contribution to the Gibbs energy due to SRO, and when set equal to the bond energy u AB , in accordance with Equation (14) they give the correct topology in FCC as well, with separate L1 2 and L1 0 ordering regions. However, Equation (13) was derived by comparing the quasichemical model with a two-sublattice model for BCC and it means that in a four-sublattice CEF model for ordering in FCC the relation for SRO between the reciprocal parameter and reciprocal energy is: when ∆G recip AB:AB: * : * = 2u AB , z = 12 and Several examples of ternary FCC ordered phase diagrams using this model are discussed in a paper by Kusoffsky et al. [39] and a simplified version is used in a commercial multi-component database for Ni-based superalloys [40].

Fe-Ni as an Example for Modeling with DFT
In the CE-CVM the SRO and LRO are described as Gibbs energy based on the atomic interactions. First the methodology of CE-CVM is explained and then the FCC phase in the Fe-Ni system is used as an example. This is then compared with the four-sublattice CEF model and with an assessment using also experimental data. The Fe-Ni system has a disordered FCC phase stable across the whole system at high T. At lower T the Fe-rich side transforms to BCC, but on the Ni-rich side there is an order-disorder transition to an L1 2 ordered phase. This is also enhanced by a ferro-magnetic transition. Recent data [41,42] indicate also that there can be an L1 0 ordered phase at low T. This system is thus interesting to use as a comparison between CVM and CEF modeling using only DFT data and a full assessment using also experimental data.

Effective Cluster Interaction
The CVM requires the effective cluster interactions (ECIs), and these values are evaluated by using the cluster expansion method (CEM). In the CEM, each ordered structure, R is represented as an arrangement of spin operators (R = (σ 1 , σ 2 , σ 3 , ..., σ i , ...)). Here, the spin operators (σ i ) correspond to each elemental species and index i represents the atomic site position in the ordered structure. By using a combination of spin operators extracted from different sites, it is possible to represent clusters in the ordered structure. For example, (σ i , σ j ) is a pair cluster consisting of sites i and j, and (σ i , σ j , σ k ) is a three-body cluster consisting of sites i, j, k.
Furthermore, the correlation functions (ξ α ) represented by Equation (22) are introduced by taking the average values of spin products for clusters such as points, pairs, three-bodies, and so on.
In the ξ α , α represents the type of cluster, and N point , N pair , N tri are the total numbers of points, pairs, and three-body sites in the structure, respectively. The energy of the regular structure can be expressed by the sum of the products of the correlation functions and ECIs (J null , J point , J pair , J tri , ...).
Ideally, it is possible to reproduce the energy of any ordered structure strictly by using the infinite number of clusters in Equation (23); however, in the actual calculation, it is necessary to truncate it into a finite number of clusters. The cluster interactions tend to be stronger at shorter distances and weaker at longer distances. Therefore, a method to set a threshold for cluster size is generally adopted. The cluster of the maximum size is given as α max , and the product of the correlation function and ECIs are summed up for the sub-clusters contained in the α max .
Since the correlation function is obtained from the ordered structure, and the total energy of the left side is obtained from DFT calculation, the unknown parameter is only ECIs. The ECIs can be determined by the method of least squares from the many number of relation expression of Equation (24).

CVM Calculation
Once ECIs are determined, the energy of any atomic configuration can be determined within the accuracy of the cluster expansion without DFT calculations. In addition, the free energy of mixing, including the entropy of the configuration, can be expressed as follows.
This equation is the second term of entropy added to Equation (24), and γ α S α is the contribution of entropy from the cluster α using the Kikuchi-Barker coefficient [1,43]. The free energy is calculated by applying the variational method to Equation (25) so that it is minimized.

CEM and CVM Calculations for Fe-Ni
The DFT calculations were performed with the Vienna ab initio simulation package (VASP) [44,45] in which ion-electron correlation is described by projector augmented wave method (PAW) [46,47]. Correlation and exchange functions were given by the generalized gradient approximation proposed by Perdew et al. (PBE-GGA) [48]. We used pseudo-potentials for Fe and Ni with 8 and 10 of valence electrons, respectively. The k-point sampling meshes of the Brillouin zone for a FCC primitive cell were 1000 points and altered the k-points meshes so that the spanning of k-points (dk) are comparable. The Brillouin-zone integration was performed using the Methfessel-Paxton technique [49] and smearing of the electron was set to 0.1 eV. The plane-wave cutoff for the wave function was 350 eV. In all calculations, we assume an initial ferromagnetic order and the unit cells of the Fe-Ni ordered structures are fully relaxed with respect to the volume, shape of the unit cell and atomic positions. The formation enthalpy (∆H) of a structure is defined as the difference between its total energy and the concentration weighted total energy of pure Fe and Ni of FCC structure.
For CEM and CVM calculation, we used the code developed in [50]. We used the ordinary tetrahedron-octahedron (TO) cluster as the basic cluster. An occupation variable of each possible site i (σ i ) takes the value +1 if Fe occupies that site and 0 if Ni occupies that site. Figure 5 shows the FCC structure and Table 1 lists the sub-clusters of the TO cluster. The specific values of the correlation functions of the L1 0 and L1 2 ordered phases are also shown in Table 1. The formation energy of each ordered phase is as follows from the values of these correlation functions and Equation (24). The optimal ECIs values for a set of 11 sub-clusters using 47 enthalpies of formation are the following (meV/atom): The ECIs with finite values consist of points, pairs and triangles. The ECIs of the other sub-clusters, including 4-6 points clusters, were given to 0.
The formation energies of L1 2 -Fe 3 Ni, L1 0 -FeNi and L1 2 -FeNi 3 from the ECIs are: The predictive error [51], averaged over structures, becomes 17 meV/atom. The information of the crystal structures used in CEM is shown in Appendix A. The calculated phase diagram is shown in Figure 6a. Table 1. Sub-clusters of the tetrahedron-octahedron (TO) clusters in FCC structure. Letters correspond to those in Figure 5. Specific values of correlation functions (ξ α ) of L1 2 and L1 0 are also listed.

CEF Calculation for Fe-Ni
Modeling based on CEF is described in detail in Lukas et al. [11]. For the FCC ordering with L1 0 and L1 2 on the Fe and Ni-rich sides we use a four-sublattice model with the following parameters taken from Equation (29) The ev2j value is the Faraday's constant to convert eV/atom to Joule/mol. A configuration independent contribution L Fe,Ni = 12, 000 J/mol was used to have the same maximum order-disorder T for L1 2 as for the CVM calculation. Equation (16) was used to calculate the phase diagram in Figure 6b and the agreement with the CVM diagram is satisfactory. The small perturbation on the L1 0 phase is due to the formation of an F' phase, with the ordering (Fe:X:Ni:Ni) where X is close to Fe, formed by a second order transition from L1 0 .

The Experimentally Assessed Fe-Ni System
The Fe-Ni system has been assessed with the Calphad method using the four-sublattice CEF model for the L1 2 and L1 0 ordering in the FCC phase by Cacciamani et al. [41], and more recently by Ohnuma et al. [42], who included new experimental data. In Figure 6c the latter assessment was used to calculate the phase diagram. The dashed line is the ferromagnetic transition line in FCC and there is a tri-critical point due to a magnetically induced miscibility gap in the FCC phase. The magnetic contribution to the Gibbs energy is based on an empirical model proposed by Inden [52] and explained in Lukas et al. [11], which depends on the experimental Curie temperature, T C and Bohr magneton number, β. With a CEF model it is possible to combine DFT data and experimental information to adjust endmember energies.

Examples Using DFT Data for CEF Models with LRO and SRO
A large number of the crystalline phases have some form of LRO, and some of them are explained in more detail below. In general it is simple to use DFT data in a CEF model when the crystalline structure is known. Phases with ionic constituents are described in Section 4.3, and how to model defects using CEF is described in Section 4.4.
For phases which always have LRO, the contribution to the Gibbs energy due to SRO is usually small and can be integrated in the LRO description. Thus SRO is only discussed in connection with phases which can also be stable while being totally disordered, such as FCC, BCC and HCP.

Interstitial Solutions, Carbides and Nitrides
A small atom like C may dissolve in the interstitial sites in metals, for example, the ferrite and austenite of Fe. This is a very important case as it gives unique properties to steel. The C atoms mix with vacancies in the octahedral interstitial sublattice in the austenite and this is a B1 structure. This means austenite is in fact the same phase as the cubic carbide in, for example, Ti-C. In the C-Fe-Ti system, austenite and the cubic TiC phase are thus modeled as the same phase. In the assessment by Ohtani et al. [53] it was described as a reciprocal system (Fe,Ti) 1 (Va,C) 1 where the endmember (Fe:Va) is the pure iron and (Ti:C) the stable stoichiometric cubic carbide. The other two endmembers (Fe:C) are a metastable cubic carbide and (Ti:Va) a metastable FCC-Ti. In the C-Fe-Ti system there is a miscibility gap between Fe (with some dissolved C) and the TiC carbide (with a significant amount of Va).
Most carbides and nitrides can also be described using one sublattice for the metallic elements and one for C and N, sometimes also including Va and DFT calculations can be used to estimate the metastable endmembers. Sometimes more complex models are used, for example the κ carbide, with a perovskite structure, was assessed in an in the Al-C-Fe system by Connetable et al. [54] using 4 sublattices for FCC on the Al-Fe side (which represent an L1 2 structure) and an interstitial sublattice site with C and Va. A full assessment would require 8 sublattices to describe also ordering for the C and Va. In a recent reassessment of this system by Zheng et al. [55] the same models for ordering were used with an improved fit to the experimental data.

Intermetallic Phases
There are a large number of intermetallic phases, Laves, A15, σ, µ, χ, etc., with different types of structures and solubility ranges. In some cases the structure is not even known but in most cases it is possible to define a CEF model which describes the stable composition range. But frequently the model has a larger composition range than the range of stability of the phase and there are often many endmembers that cannot be determined experimentally. Before DFT was available the number of crystallographic sites were often reduced to make it possible to estimate the endmember energies by extrapolations from the stable range. Even today it can be a cumbersome task to determine all endmembers in a multi-component system for a phase with more than three sublattices when there are three or more components. As already mentioned, many of these endmembers may be mechanically unstable which means the value calculated by DFT is meaningless, see Section 5.2.1.
In order to combine separate assessments of systems including intermetallic phases in a database the endmember energies for the pure elements must be the same in all assessments in such a phase. As already mentioned several such energies for TCP phases have been published by Sluiter [30].

The σ phase and other TCP phases
The σ phase has five sublattices but in many commercial databases the σ phase is modeled using three sublattices, (A,B) 10 (B) 4 (A,B) 16 with a restricted composition range. Here the elements A are those with preferred FCC lattice, B those with preferred BCC lattice. The reason for this simplification was to reduce the number of endmembers and thus difficulties assessing endmembers values and to reduce calculation time.
Today this restriction is less valid, as the thermodynamic software are better and computer hardware is also much faster. Thus a five-sublattice model is the natural model to use to take into account the crystallographic structure for the σ phase. In a multi-component system this may mean several 1000 endmembers but only a fraction of these are needed in the description of the stable range of the σ phase because one can ignore those that have small chances to be a stable endmember as discussed in Section 5.2.
Many TCP phases such as µ, χ, A15, Laves, etc., are modeled differently in different assessments. This makes it difficult to combine the assessments to databases and the general recommendation is to use the crystallographically correct number of sublattices even if one can simplify the model in a particular system.

The Effective Bond Energy Formalism
An interesting new development is the Effective Bond Energy Formalism (EBEF) recently proposed by Dupin et al. [9]. EBEF assumes ideal mixing of the constituents on each sublattice but the endmember parameters are replaced by bond energy pairs between constituents in different sublattices. This can reduce the time for equilibrium calculations and the number of DFT calculations needed for endmember parameters multi-component systems. It also has the potential of improving the predictions from lower order to higher order systems but more work is needed to test this model.

CEF Models with Ionic Constituents
Oxygen behaves differently from C or N because it will normally capture electrons from the metallic elements. S is also prone to capture electrons and compounds with F, Cl etc are also strongly ionic but there are fewer solution phases with these elements. There has been several assessments for different oxide systems using CEF-for example, [56][57][58][59][60][61][62][63] using both experimental and DFT data.
In this context we will only look at some aspects of modeling the UO 2 phase with C1 structure which is important as nuclear fuel. The phase diagram for UO 2±x shown in Figure 7a is calculated from an assessment by Guéneau et al. [64] as part of the development of the Thermodynamics of Advanced Fuels -International Database (TAF-ID) [65,66]. The model for the UO 2±x phase with a C1 structure, in Figure 8a, takes into account the wide composition range and that U can have several valences and that a small amount of O can also occupy octahedral interstitial sites: where the third sublattice represent the interstitial sites. There is an enormous amount of experimental data on this system and there are also DFT calculations [60,67,68] of defects but so far there is no CVM based model which can consistently describe all data with the same accuracy as the CEF model. For simulations of the fuel behaviour it is also important that the model can describe the solubility of 10 or more fission products and equilibria with new phases formed during the burnup of the fuel. The CEF model for UO 2±x has 12 endmembers but only one of these is electrically neutral, all other has a net charge. To understand the relations between endmembers the prisms in Figures 8b,d can be helpful. They both represent the C1 phase but the endmembers are joined in different order. The endmembers are denoted "4OV," "4OO," "5VV," etc., to specify the constituents in all three sublattices. The shaded surfaces inside the prisms are the regions where the phase is electrically neutral. With some efforts one can figure out that the two separate shaded areas in Figure 8d are the same area as in Figure 8b.   Figure 8c illustrate a defect when one oxygen ion moves from the tetrahedral to the octahedral sublattice which can be related to the difference between endmember energies in Equation (32). The defect in Figure 8e represent a change of the valency of two U ions: 2 U +4 = U +3 +U +5 . The use of these defects is explained in more details in [69].

CEF Models for Defects
There is no real difference between modeling defects and solubilities in a phase. Defects can be substitutional, interstitials or antisite and the latter two can be modeled as a reciprocal system. Traditionally very simple models like the Wagner-Schottky [70] have been used to handle small solubilities but it cannot be extended to multi-component systems, to higher fractions of defects or when several defects may interact, as is the case for the defects in UO 2±x in the previous section. The reciprocal model for the Laves phase was discussed in Section 2.1 and it can be extended to higher fractions of defects by additional parameters. Defect models and in particular thermal vacancies are discussed by Rogal et al. [71].

Phases with Order-Disorder Transitions
The FCC, BCC and HCP phases have a large number of ordered superstructures. Most of them, such as D0 22 , D0 23 , etc., can be treated as separate phases but in a few cases the equilibrium between the ordered and disordered structures are technically important and an accurate description requires that both the ordered and disordered phases are described with the same Gibbs energy model.

BCC ordring using Two sublattices
The B2 phase is important in many alloys and the A2/B2 transition can be modeled with two sublattices. But if the Heusler phase, L2 1 , is important one must have four-sublattices [8]. A particular problem with the B2 phase is that one frequently has vacancies on one side of the ideal ordering. For example in the Al-Ni system the B2 is stable in the middle as shown in Figure 9a and on the Al rich side, the Ni atoms are replaced by vacancies, not by Al atoms. This requires that one includes thermal vacancies already in the model for the disordered BCC phase as explained in the assessment of Al-Ni by Ansara et al. [34].

Never model FCC ordering with Two sublattices
The two-sublattice model for FCC ordering using (A,B) 0.75 (A,B) 0.25 is obsolete today. It cannot describe the L1 0 and in multi-component systems it requires a very complicated set of additional ternary and higher order parameters, related to the different binary bond energies. If these are not correct the "disordered state" with exactly the same constituent fractions on both sublattices will not exist. It is very easy to forget some of these parameters and end up with an FCC phase that is never completely disordered and this is not evident in a calculated phase diagram. A careful check must be made that the constituent fractions in the disordered state are identical in both sublattices.

The Four-Sublattice model for order-disorder
The four-sublattice model can represent a tetrahedron in the FCC lattice, positions A-B-F-G in Figure 5. It is symmetrical with all nearest neighbors on another sublattice and it can describe the ordered L1 2 and L1 0 phases and calculate multi-component equilibria almost as fast as the unsymmetrical two-sublattice model if the software takes the symmetries into account. It is possible to add an interstital sublattice to this model representing the octahedral interstitial sites for C and N together with Va. New assessments of FCC ordering should always use the four-sublattice model to provide a better description of FCC ordering for future applications. The topology for the order-disorder transitions should be assessed also in the metastable composition ranges as for the FCC phase in Figures 6b and 9b in order to have reasonable extrapolations to higher order systems. If there is a need to be compatible with an old database using the two-sublattice model this can be assessed at the same time.
The ordered HCP phase can use the same four-sublattice model for B19 and D0 19 ordering. In addition to the Au-Cu system there are a few more assessments using the four-sublattice model, the Al-Ni system shown in Figure 9a from [72], Al-Ni-Pt by Lu et al. [73] and Fe-Ni by [41,42]. The SRO contribution using four sublattices is obtained from the reciprocal energy, as described in Section 2.4, but as there are several such subsystems, the "constituent cube" in Figure 9d is an attempt to explain how this can be done. All the ordered structures A 3 B, A 2 B 2 and AB 3 are included in the Figure 9d (34) using the endmember parameters in Equation (19). This reciprocal energy can be used in the reciprocal parameters to approximately describe SRO in this system and because all sites are equivalent there are two more unique reciprocal systems which give three reciprocal parameters to describe SRO in the four-sublattice CEF model: L X:X:A:A , L X:X:A:B = L X:X:B:A and L X:X:B:B , where X indicate a sublattice with interaction. These can be different because the bond energy, u AB , may be different for x B = 0.25, 0.5 and 0.75. In Equations (19) and (31) u AB was assumed to be a constant.
For the BCC phase the four-sublattice model is needed for the D0 3 and Heusler phase, L2 1 . It is more complicated than FCC because the BCC tetrahedron is unsymmetrical and the four-sublattice BCC model can also describe the B32 structure and there are four independent reciprocal parameters. The first assessment of the Al-Fe phase using four sublattices for BCC ordering and the ferro-magnetic ordering is shown in Figure 10a, assessed by Sundman et al. [23]. The calculated heat capacity including magnetism and configuration at x A = 0.3 as a function of T is shown in Figure 10b and in Figure 10c the variation of the constituent fractions. Several crystal structures and related endmembers based on the four-sublattice BCC model are explained in [8,74].

Modeling SRO in Liquids
Liquids are outside the scope of this paper, but the liquid phase is present in most phase diagrams and it is a kind of baseline for the crystalline phases. Metallic liquids with no or small SRO can usually be modeled as a substitutional regular solution, and there are three basic modeling techniques for SRO in liquids: 1.
The associated model which add constituents with a stoichiometry close to the SRO composition and include interactions between these and the constituents representing the elements, all mixing randomly. This overestimates the configurational entropy and requires many interaction parameters between the constituents. It includes substitutional regular solutions as a subset.

2.
The ionic two-sublattice liquid model [75] is a reciprocal model with cations on one sublattice and anions, vacancies and neutrals on the other. It is similar to a CEF model but the site ratios are variable. This model includes the substitutional regular solution and some associated models as subsets but can have problems with reciprocal miscibility gaps. 3.
The quasichemical model which is the simplest CVM based model has been modified by Pelton et al. [76] to exclude LRO and use only two bonds per constituent to avoid negative configurational entropy. The stoichiometries of the constituents are selected according to the SRO composition.
It has some problems incorporating substitutional regular solutions.
Each of these models has their advantages and disadvantages. More research based on both experiments and Molecular Dynamics is needed to find a liquid model which is really satisfactory.

Assessments Using CEF Models Including ab initio Data
In this section the methodology using CEF models is discussed using previous studies. There are several criteria for selecting a model for a crystalline phase when assessing a system: Physical soundness.

2.
Compatibility with existing assessments that should be combined with the new assessment.

3.
Amount of experimental and theoretical data. 4.
Extrapolations to higher order systems.
It is straightforward to combine results from DFT calculations with experimental information in CEF models. The DFT values are treated as any other experimental data on the thermodynamic properties such as enthalpies of formation or mixing, chemical potentials or heat capacities as well as phase diagram data when determining the model parameters. The most frequently used DFT results are energies of formation of endmembers at 0 K, but it is possible to use heat capacities and other properties calculated by DFT.

The Number of Sublattices to Use
There are still problems, including all crystallographic features in a sublattice model and in some cases simplifications being reasonable. For an assessment of a system which will not be used in any database, one should of course use the best physical model possible but if it is part of an ongoing database development there is little reason to use better models than already used in the database. However, new assessments of order-disorder transitions in FCC, BCC or HCP phases should always use four sublattices for the metallic sublattices, as discussed in Section 4.5.2. This will give improved extrapolations to multi-component systems, as shown in [8,73].
For TCP phases it is now possible to use the crystallographically correct number of sublattices combined with a disordered fraction set. As explained below, one can avoid calculating a large number of endmembers by DFT because only a few are needed to have a reasonable multi-component Gibbs energy surface for the stable phase.
In addition to the endmemer parameters, one can use interaction parameters in the disordered fraction set. Interaction parameters in the sublattices should be used only in phases with order-disorder transitions or CEF models describing phases with interstitials or defects, such as the Laves phase, where the two sublattices have very different coordination numbers, as discussed by Crivello et al. [77].

Reducing the Number of Endmember Parameters
A phase with several sublattices and many constituents will have several thousand endmembers in a multicomponent system. However, an endmember with a Gibbs energy far above the disordered Gibbs energy surface (sometimes called the convex hull), at the composition of the endmember, will never have any influence on this surface and can be ignored. Only those endmembers with a Gibbs energy below or close to the disordered Gibbs energy surface of the phases are needed. In Figure 11a all 47 DFT calculations in Appendix A are plotted against the composition in mole fraction. For the CEF calculation in Figure 6b only those representing the lowest endmember energies are needed. In Figure 11b the energies of all 32 endmembers in a five-sublattice model for Cr-Fe are plotted from a calculation by Korzhavyi et al. [5], and in Figure 11c there are a large number of DFT calculations for several phases in the Re-Ta system from a paper by Palumbo [78]. Very few of them are needed for the lowest Gibbs energy curve for each phase using a CEF model and it should be possible to adapt the software for DFT calculations to use crystallographic information and already calculated bond energies to avoid spending effort on calculating properties of endmembers that are likely to be significantly above the "convex hull" of a phase.

Mechanically unstable endmembers
When calculating the formation energy for an endmember by DFT it is not unusual that the result represents a mechanically unstable structure. This is evident if the structure changes by relaxing it or if the calculated phonon frequencies of the endmember are negative. As an example, Figure 12 shows the calculation results for the Laves_C14 of (Fe,Mo) 2 (Fe,Mo). The calculation results of Fe 2 Mo and Mo 2 Fe, which are the end members of (Fe,Mo) 2 (Fe,Mo), are shown in Figure 12a,b, respectively. Fe 2 Mo, which is a stable structure, has positive frequencies for all phonon dispersions. On the contrary, due to the mechanically unstable, negative frequency, modes appear in the phonon dispersions for Mo 2 Fe. Any thermodynamic value calculated for such a mechanically unstable endmember is meaningless. If the endmember is necessary for the model its values must be estimated in some other way, for example, using Equation (10), as discussed in Section 5.2.2. There are cases when a phase is mechanically unstable at 0 K but can become stable at finite T, such as BCC for pure Ti. This is discussed in the paper by Antolin et al. [79] and also by [17].

Estimating Endmember parameters
Before DFT became an easily available tool, it was a complicated task to estimate metastable endmembers, in particular, endmembers such as • G Laves_C14 Fe:Fe which is shared among several different binaries and must have the same value in all. A long time ago it was agreed [80] to assume that the pure elements in the Laves phases should have an enthalpy difference between the stable pure element, GHSER, and the Laves phases equal to 5000 J/mol atoms, but this has now been replaced by the DFT calculations by Sluiter [30].
In a paper by Hillert et al. [81] it is explained how one can divide a complex CEF model into separate reciprocal systems and by setting ∆G recip AB:CD = 0 and using Equation (10) to express the metastable or unstable endmembers as functions of three other endmembers. See also Section 4.3 for how to estimate some endmembers in the C1 phase.

Conclusions
Based on the explanations and discussions above we make the following conclusions: • DFT is a powerful method to obtain many different kinds of information about materials from basic physical theories. • The Calphad method uses a mean field theory to model the Gibbs energies of different phases, mainly based on experimental data, in order to calculate phase diagrams and extract the thermodynamic properties of a system.