Computing solubility products using ab initio methods

The solubility product of NbC in low alloyed steel is computed from electronic density functional methods including the effects of electronic, vibrational, and magnetic excitations. Although many simplifications are made in the computations, agreement with experimental data is within the scatter of the latter. The T = 0 K terms dominate in the determination of the solubility product but vibrational and magnetic contributions play a significant role also while electronic excitations can be neglected. Supercell calculations were shown to be poorly suited for determination of embedding enthalpies of solutes in bcc Fe.


Introduction
In recent years, ab initio modeling of steel has begun in earnest. Ab initio research on steel started somewhat later than on other alloys, such as aluminum alloys, for a number of reasons: the presence of magnetism in steel presents a formidable challenge that is not present in aluminum alloys; interstitial solutes such as carbon in steel mean that computations require more memory and more computer time as compared to most other alloys; displacive structural phase transformations play an essential role in steel while such transformations are absent or of much lesser importance to other alloys where transformations might be of order-disorder type. Hence, ab initio modeling of steel provides a great challenge.
Ab initio methods have made significant inroads already in the field of computational thermodynamics, often referred to as ''CALPHAD''. In this field it is becoming increasingly common to augment the experimentally known thermochemical data with ab initio computed compound formation enthalpies. The ab initio computed formation enthalpies pertain to a temperature of 0 K, for which experimental data can be obtained only through extrapolation of finite temperature data, which has often proved less reliable than the ab initio data. The T = 0 K ab initio enthalpies strongly restrain the modeling and parameter choices at finite temperatures and so it is believed that this marriage of ab initio and experimental data gives more reliable thermodynamic representations. The CALPHAD method is a sophisticated approach in which phase equilibria and other thermodynamic properties can be computed for many-component alloys. Obviously, such versatility makes for complex models that obscure the actual role and contribution of ab initio modeling. Therefore, here a very basic, but high practical property is considered, the computation of solubility products. A solubility product for a certain compound forming from a solution is the greatest value that the product of the concentrations of the constituents of the compound can take in that solution. If the actual product of concentrations is smaller than the solubility product, the compound will not be thermodynamically stable while in contact with the solution, and vice versa if it exceeds the solubility product constituents in the solution will lead to growth of the compound phase. Of course, the solubility product is generally a function of temperature. The dependence on other variables than temperature is generally much less in solid state systems. For this reason the solubility product is T. Klymko a very basic measure that informs us whether at a certain temperature precipitates are likely to form, or are likely to dissolve.
Below, the solubility product of the compound NbC in ferrite (bcc Fe) and in austenite (fcc Fe) will be examined, because the precipitation of NbC is an important factor in achieving good properties in a wide range of low-alloyed steels. A rather simplistic thermodynamic model will be employed, that is a caricature of the more sophisticated models that are used in actual CALPHAD modeling, but with the compelling advantage that the main physical effects can be taken into account while maintaining mathematical simplicity. In addition, a simple model has relatively few variables so that the role of the ab initio input data is particularly evident.

Thermodynamic model
The free energy G of a Fe-Nb substitutional solid solution based on the bcc crystal lattice, per bcc site, can be written as, where x Fe (x Nb ) is the fraction of Fe (Nb) atoms per bcc lattice site, G a [Fe] (G a [Nb]) is the free energy of pure bcc Fe (bcc Nb), and where DG a mix is the free energy of mixing. Assuming that the number of vacancies on the bcc crystal lattice is negligible, the sum of all fractions equals unity, x Fe ? x Nb = 1 in this case. In the absence of interstitials, a fraction per bcc lattice site is equivalent to an atomic concentration.
G a [Fe] and G a [Nb] are temperature dependent, G a = H a -TS a , with implicit temperature dependence in both the enthalpy H a and the entropy S a . Recently, there has been much success in determining the temperature dependence fully ab initio by considering various thermal excitations, e.g., for fcc Al [8], bcc Fe [17,18], and cementite Fe 3 C [7]. DG a mix can be split in an ideal entropy part and an excess mixing part, where the excess part behaves as a polynomial in the composition, with S id ðxÞ ¼ Àk B ½x lnðxÞ þ ð1 À xÞ lnð1 À xÞ; ð3Þ In these equations g 0,Nb a , g 1,Nb a ,... are expansion coefficients with reference to the bcc pure elements at the same temperature, k B is the Boltzmann constant. For small x Nb the first-order term g 0,Nb a should be adequate and we shall refer to it as the embedding free energy DG a emb ½Nb of element Nb in the matrix a. The most important physical origins of the excess term are chemical interactions, elastic distortions, magnetic, and vibrational effects, as well as electronic entropy effects in the case of metallic materials, and, for more concentrated alloys also short-or long-range order.
C is an interstitial solute in bcc Fe which takes the octahedral interstice of which there are 3 per bcc site. We now define a property x C as the number of carbon atoms per bcc site, this means that x C = 3 in the hypothetical situation where all octahedral interstices are occupied by C. The free-energy G of a bcc Fe-C solid solution per bcc site, where G d [C] is the free energy of pure C diamond, and DG a mix ðx C Þ is the free energy of mixing. We have selected diamond and not graphite as our reference state for carbon because its structural and vibrational properties are very well-described by ab initio methods, and also because the free-energy difference between the two states is small compared to the formation and embedding energies we shall encounter below. In analogy with Eq. 2 DG a mix ðx C Þ can be separated into an ideal mixing entropy term and an excess term, with f C = x C /3, and In a similar fashion as for Nb we can define the embedding free energy DG a emb ½C with where we used the carbon atomic fraction Ignoring the Wagner interaction [64] of the dilute Nb and C species, the free energy of the dilute solid solution per bcc site can be written as The free energy of the dilute fcc solid solution can be derived analogously, provided that we take into account that C occupies octahedral interstices of which there is only one interstice per fcc site, The fcc free energy can be written relative to that of the bcc solid solution, where the superscripts ca are shorthand for differences between properties pertaining to fcc and bcc. The embedding free energy for a species ''Q'' can be divided in several contributions, which are described in detail below, where we recognize vibrational, magnetic, and electronic contributions. The first-term concerns an enthalpy because the configurational entropy is treated separately as its derivative is ill-behaved in the limit c i = 0. The free energy of the stoichiometric NbC phase can be represented as G [NbC]. Its formation free energy per formula unit can be found easily from the general formula for a compound P n Q m where elements P and Q are in their reference states. In actuality, NbC precipitates in fcc or bcc Fe are not stoichiometric but have a considerable concentration of vacancies on the C sublattice as can be deducted e.g. from the detailed experimental data in Tables III [24] and II [47] and as studied in detail by Balasubramanian et al. [3]. Often it is assumed that the composition is close to NbC 0.87 , i.e., about 13 % of the C sublattice is vacant. For simplicity, we will assume stoichiometry. The assumption can be overcome relatively easily, but we will leave it for a future contribution. The solubility product for stoichiometric NbC in ferrite (bcc) is now derived easily for the chemical reaction where DG vanishes at equilibrium. It is computed according to a tangent of the solid solution free energy intersecting with the stoichiometric precipitate free energy Expanding it to lowest order in x Nb and x C , and solving DG ¼ 0 then gives, where the subscript ''lmt'' indicates a solubility limit which occurs at equilibrium between solid solution and NbC.
is recognized as the formation free energy of NbC relative to the pure elements DG f ½NbC, see Eq. 14, so that we may write the solubility product as The solubility product for stoichiometric NbC in austenite (fcc) is defined by reference to ferrite as Furthermore, to facilitate comparisons with experimental data, it is useful to define the solubility products with reference to concentrations expressed in weight percent ðKÞ rather atoms per crystal lattice site (K) according to where Q may represent C, Fe, or Nb, and m Q is atomic mass of atomic species Q, so that where in the RHS we used x Nb , x C ( 1. Experimentally, solubility products appear over some temperature range to be rather well represented by straight lines in Arrhenius plots. Therefore, the parameters A and B in have been widely reported in the literature.

Magnetic effects
It should be mentioned that the magnetic disordering that occurs in ferrite in the temperature range of interest here J Mater Sci (2012) 47:7601-7614 7603 strongly affects the free energies [1,5,13]. When considering the effect of magnetic disordering on the excess mixing properties, such as DG a emb ½Nb and DG a emb ½C, the effect can be very significant also. This is illustrated by using the ''proportional conversion method'' of Nishizawa et al. [34] which has been applied to Fe-C and other systems [9,37]. The ''proportional conversion method'' describes the effect of alloying as, where c Q is the atomic fraction of the atomic species ''Q'' (Ni and Co species are not counted in the sum), T C * is the Curie temperature of the Fe-rich bcc alloy, T C (0) is the same for pure bcc Fe (1043 K), and T * = TT C (0) / T C * (this is Eq. 49 in Ref. [34]). For many sufficiently dilute alloys the Curie temperature T C * is well described by a linear function of the composition where DT C;Q is the negative slope of T C * associated with species ''Q''. Assuming that the alloying elements in their reference states are not magnetic, the magnetic free energy then is where C . To relate G mag * (T) to the mixing free energy we need to remove the linear part (1 - As r Q is the only solute-specific parameter, the magnetic effect on the embedding free energy of species Q, In Fig. 1 G mag (0) (T) is visualized using the pure bcc Fe parameters from Ref. [5]. For elements that strongly reduce T C , so that r Q is of order unity, it implies that the magnetic contribution to the embedding free energy is of the same order as the magnetic free energy of pure bcc Fe.

Electronic excitation effects
At finite temperature the Fermi-Dirac distribution broadens around the Fermi-level which signifies that some electronic states above the Fermi-level are populated and some states below it are vacated. Obviously this has consequences for the energy and entropy of the system. For insulators and semiconductors the presence of a band gap at the Fermi-level makes electronic excitation effects negligible, but for metallic materials the Sommerfeld approximation [2] gives a good representation of the free energy contribution where nð F Þ is the density of electronic states at the Fermilevel, see Table 1. The contribution to the embedding free energy is computed by approximating the composition derivative with a finite difference where Q represents either Nb or C, and m, n, and p are chosen sufficiently large. In this case we selected for Nb; m = 53, n = 54, p = 54, and for C; m = 54, n 54, p = 64.
The vibrational contribution to the formation free energy of NbC is given by a similar equation where for NbC a Nb 32 C 32 supercell was used. Vibrational effects are generally considered to be the origin of the close-packed to bcc transitions with increasing temperature in several metals because the relatively open bcc structure tends to have a greater vibrational entropy than the close-packed structures. Computationally, vibrational excitation effects are usually included in the harmonic or the quasi-harmonic approximation. In the latter approximation, the crystal structure is allowed to expand and relax as the temperature is increased but otherwise it corresponds to the harmonic approximation. In the earliest implementations the Debye model was often used [32,66].
Here, we will use the so-called direct method [39,53] where interatomic force parameters are determined in a supercell by moving each individual atom away from its equilibrium position, one atom at a time, and monitoring the resulting forces on the remaining atoms. The force parameters are then Fourier transformed to obtain a dynamical matrix which, upon diagonalization, yields the vibrational density of states g(x) as function of the vibrational frequency x. The vibrational density of states provides the free energy due to vibrations[53] where m ¼ expðÀ" hx=k B TÞ. DG vib;emb is evaluated through finite differences just as for the electronic free-energy contribution, see Eq. 29.

Ab initio model parameters
As seen in the section above various terms are needed to compute solubility products. The formation free energies of precipitating phases DG f ½NbC, embedding energies DG a , and energy parameters specific for the bcc-fcc transformation DG cÀa : Strictly speaking only ground state parameters, pertaining to T = 0 K, can be obtained in a straightforward manner for first-principles computations. Vibrational excitations, in the harmonic or quasi-harmonic approximation can be included without great difficulty, and electronic excitations also do not pose great challenges, see e.g., Ref. [8]. Including magnetic excitations ab initio is more challenging, although here too progress in being made, see e.g., [7,17,18]. Here, we will treat the vibrational and electronic excitations ab initio, but the magnetic excitations will be treated through the empirical approach of Nishizawa et al. [35]. The ground-state properties can be computed straightforwardly with any modern electronic structure code, here we used the Vienna Ab initio Software Package (VASP) 4.18 [22]. We used the projector augmented-wave (PAW) [4,23] type pseudo-potentials and a generalized gradient approximation (GGA) type exchangecorrelation functional [42]. The Nb pseudo-potential-treated 11 electrons per Nb atom as valence electrons (4p 6 , 5s 2 , and 4d 3 ), the C pseudo-potential-treated 4 electrons per C atom (2s 2 and 2p 2 ), and the Fe pseudo-potential considered eight valence electrons per Fe atom (4s 1 3d 7 ). Accuracy was set to high, the plane-wave expansion of the wave functions included waves up to 400 eV kinetic energy. Integrations in reciprocal space were carried out with a spatially homogeneous C centered grid such that the product of atoms in the cell and k-points in the first Brillouin zone was about 10000. Structural optimizations were deemed converged when the magnitude of the force on any atom core was less than 0.5 meV/nm, and each individual stress component was less than 0.3 GPa and typically much smaller than this limit. Table 1 lists the structural and energetic parameters for the pure elements and NbC in their ground states. Typically, lattice parameters are reproduced to within 1 or 2 %, and while the total energies are known to contain sizeable systematic errors, these errors largely cancel so that formation enthalpies, DG f ½NbC(T = 0 K), are generally accurate up to about a 1 kJ/mole atoms. Nevertheless, there is a very significant difference with the experimental value, even if we consider that the latter pertains to a temperature of 298 K rather than 0 K. The experimental NbC formation -137.6 Ref. [12], -133 Formula unit is abbreviated as f.u.
J Mater Sci (2012) 47:7601-7614 7605 entropy is 3.3 J/mol atom K [67]. Even if we neglect that this entropy tends toward zero as the temperature approaches 0 K, this temperature difference accounts only for a few percent of the difference between ab initio and experimental free energy change.
Embedding free energy of Nb in ferrite The embedding free energy of Nb in ferrite at T = 0 K can be computed ab initio with a variety of methods, such as via Fe bcc supercell calculations in which a single-Nb atom is placed among a great number of Fe atoms [51,55], by determining the mixing enthalpy from special quasi random structures [55], and by determining the mixing enthalpy from a cluster expansion (CE) for bcc Fe-Nb substitutional alloys [55]. As the second and third approaches are both based on the CE formalism, we will consider the first and second methods here. Supercells can be used to model dilute solutions, particularly when we consider a series of supercells of increasing size. Such a sequence of larger and larger cells then models increasingly dilute solid solutions, and by taking a limit the infinite solution can be evaluated [30,54,55] using arguments derived from linear elastic theory. Furthermore, in the case of dilute alloys one can consider both fixed-volume (energy) and -pressure (enthalpy) cells. Fixed volume cells, when viewed as being embedded in the matrix material represent the case of the embedding material being infinitely stiff, whereas fixed-pressure cells (typically with P = 0 GPa) represent embedding in an infinite-compliant matrix material. As can be expected, the most realistic representation is some intermediate in between these two extremes, although a little tilted toward the fixed pressure side [30]. Supercells can be created in a variety of ways. In Fig. 2 we show supercells with sc, bcc, and fcc, arrangements of lattice points. Supercells as a representation of dilute solid solutions has as main shortcoming artificial periodicity, and consequently highly specific directions in which the dilute species ''feel'' each other. Figure 2 illustrates that the dilute species interact primarily along h100i, h111i, and h110i directions in the sc, bcc, and fcc type supercells, respectively. In highly anisotropic materials this means that sc, bcc, and fcc type supercells will each have their own convergence toward the infinite-dilution limit. This is evident from Fig. 3 where sc, bcc, and fcc type cells clearly follow distinct trends. This is not surprising in the case of ferrite, it is elastically far from isotropic (A = 2.41 [49]) and the magnetization is strongly coupled to elastic deformation. Hence, it is to be expected that the interaction between Nb atoms with their strong elastic distortions is highly sensitive upon the direction of interaction. It appears that for constant volume cells fcc and bcc type cells are best, whereas for constant pressure cells sc type supercells are optimal. When a series of supercells is considered of various types without regard for the various supercell types a rather jumbled sequence of formation energies and enthalpies is found, as shown in a recent study of Fe-Mo alloys by Lejaeghere et al. [26]. Ozolins and Asta [38] observed oscillatory interactions in dilute solutions of Sc in Al which they attributed to a combination of geometric volume expansion and electronic Friedel oscillations. While electronic Friedel oscillations do not provide an explanation here, one could suspect oscillations in the spin density. However, this did not show up in our analysis of the spin density in the various cells considered: Nb atoms aligned antiferromagnetically with the Fe matrix and the local moment centered on the Nb atom, computed as the spin density integrated over the Nbcentered Voronoi volume, was about -0.7 l B per Nb atom. Aside from the local moment on the Nb atom the spin density in the Fe matrix appeared rather unperturbed from the case without Nb, with the spin density strongly peaking near the Fe atoms and the spin density taking moderately negative values in the interstitial regions, without showing any sign of significant anisotropy. Therefore, lattice relaxations in elastically strongly anisotropic bcc Fe are expected to be main reason for the jumbled sequence of formation energies and enthalpies. A much more troubling aspect of Fig. 3 is that while constant volume and constant pressure results do seem to converge to the same infinitedilution limit for each type of supercell, as in the study of Mishin et al. [30], the sc, bcc, and fcc type supercells each converge to their own distinct limit. It is evident that bcc, fcc, and sc type cells when extrapolated to n -1 = 0, which corresponds to DH a emb ½Nb, give values of about -2.4, -5.2, and -10.9 kJ/mol, respectively. This is the case not only for the embedding energy, but also for other properties such as the magnetization change per Fe replaced by Nb. This is in contrast to earlier defect energy supercell calculations on non-magnetic materials, where generally accurate and consistent results were obtained [30,51,54,55]. In keeping with the previously mentioned role of elastic relaxations, it should be noted that the bcc, fcc, sc sequence is seen also in the elastically stiffest directions, which are in decreasing order h111i, h110i, and h100i. When considering all constant volume results we extract a DH a emb ½Nb value of -7.01 kJ/mol with a squared correlation coefficient of 0.60, while the zero pressure results give a DH a emb value of -10.06 kJ/mol with a squared correlation coefficient of 0.17. The energies can be extrapolated more reliably than the enthalpy results-presumably because changes of volume affect the systematic error in the enthalpy calculations. As the supercell results are not conclusive, another method is used also: evaluation of the mixing enthalpy in the dilute limit using special quasi-random structures (SQS) [69]. Details of how these SQS were designed will not be reported here, the description of these structures is given in Table 2. The ab initio computed formation enthalpies are shown in Fig. 4. As described in Ref. [55], the mixing enthalpy can be obtained by fitting a low order polynomial in the composition to the SQS enthalpies. Both third-and fourth-order polynomials were considered giving the following expressions for the mixing enthalpy, SQS18bcc-A 12 B 6 (a) SQS18bcc-A 12 B 6 (c)    More important than this difference between the two fits is the fact that it differs in sign from the results obtained from the supercell calculations. In some alloys, such as Fe-Cr [29,68], it is well-known that mixing and unmixing tendencies may change with composition so that dilute and more concentrated alloys behave differently. Such tendencies have never been reported for Fe-Nb however, and the isoelectronic Fe-V system has been reported to not feature such a change with composition [28,29]. Although Fe-Nb is a compound forming system, and the bcc-based FeNb B2 and Fe 3 Nb DO 3 compounds form weakly exothermically from the pure elements, the actually occurring Fe-Nb compounds are not bcc superstructures but structures stabilized by atomic size differences. The size mismatch between the small Fe atoms and the much larger Nb atoms makes unmixing tendencies in the disordered solid solution more plausible. It should be remarked that the SQSs derived mixing enthalpy also agrees in sign with that obtained from thermodynamic modeling of the experimental phase equilibria [40] over the entire composition range. A reason for the apparent discrepancy between supercell and SQS derived mixing and embedding enthalpies may stem from the magnetic order. The magnetic order relevant for experimental observations is invariably in the neighborhood of the Curie temperature and thus includes significant magnetic disorder. Such disorder is known to strongly affect mixing and unmixing tendencies, as is apparent from studies on Fe-Cr (see Fig. 3c) [19]. One could speculate that in the SQSs the magnetic disorder which follows from the quasi-random distribution of Nb atoms, to some extent mimics the solid solution near the Curie temperature. The supercell results might then strictly speaking be correct at T = 0 K, but actually be irrelevant to the actual finite temperature Fe-Nb solution. Therefore, we prefer the SQS-based results over the supercell results. Hence, we proceed with an embedding enthalpy at T = 0 K for Nb of 20 ± 2 kJ/mol (average of the third-and fourthorder polynomials).
The effect of magnetic disordering to the embedding free energy can be evaluated with Eq. 27. The important parameter here is the rate at which the Curie temperature changes with Nb concentration. This property was measured by Stoelinga et al. [57] to be an increase of 0.8 K per atomic percent Nb. If we assume that the Curie temperature exhibits the same compositional trends as the magnetization per atom and make an estimate of the latter on the basis of the SQS results a contrasting picture emerges, see Filled, open, and gray symbols with black dot indicate sc, bcc, and fcc type supercells. Note that 1/n corresponds to the atomic fraction of Nb, so that extrapolation to 1/n = 0 yields the embedding enthalpy Fig. 5. The magnetization per atom of the Fe-Nb SQS structures in the limit of vanishing Nb concentration decays with 0.033 l B per atomic percent Nb. In amorphous Fe-Nb alloys also a significant decrease in the Curie temperature is measured [61]. Here, we will use the Stoelinga results [57]. We thus arrive at an estimate of r Nb = 0.8 9 100/ 1043 = 0.077 which indicates that near the Curie temperature the embedding free energy is decreased by 2 Â ð0:077Þ Â ð1:6Þ % 0.25 kJ/mol. Clearly, magnetic disordering contributes a rather minor term to the embedding free energy.
The vibrational effects have been evaluated also. Figure 6a shows the vibrational free energy per atom computed for Fe 54 , Fe 53 Nb, Nb 54 , C 64 , (NbC) 32 , and Fe 54 C using Eq. 30. The very stiff NbC and especially diamond C lattices are evident from the much smaller vibrational free energies of these structures. The curves for Fe 54 ,Fe 53 Nb, and Fe 54 C largely overlap. The vibrational contribution to the Nb embedding energy is displayed in Fig. 6b. Nb expands the Fe lattice and thereby introduces lower frequency modes so that the vibrational contribution to the embedding free energy, particularly at higher temperatures is negative (exothermic). In the temperature range from 500 to 1100 K, the vibrational contribution is well described by DG a vib;emb ½Nb ¼ À0:1011 À 0:0042 T kJ/mol. At T = 1000 K, the vibrational contribution is seen to be about -4 kJ/mol. However, bear in mind that this result is based on a supercell calculation, which in the previous subsection was shown to be unreliable. The vibrational contribution has not yet been evaluated via the SQS structures.
Electronic contributions to the embedding free energy can be obtained from Eq. 28. Fe n Nb supercells give a nð F Þ Note thatK½NbC is defined with respect to weight percent Nb and C in ferrite  typically about 4 states/eV less than in the sum of nð F Þ of Fe n and Nb 1 . Analysis of nð F Þ of the SQS structures shows some scatter, but a least squares fit to the second-order polynomial in the composition gives a decrease of about 1.3 states/eV per Nb atom at infinite dilution. Again, a significant difference between supercell and SQS results, with the latter appearing more likely as pure Fe and pure Nb have nð F Þ values of only about 1 and 1.5 states/(eV atom). Taking the SQS results, we obtain thus DG elec ¼ 1:5 Â 10 À6 T 2 kJ/mol Nb dissolved in ferrite. At T = 1000 K this amounts to an endothermic contribution of about 1.5 kJ/mol. The magnetic disordering, vibrational, and electronic contributions are not small relative to the ab initio computed embedding enthalpy at T = 0 K. Figure 7 shows the total embedding free energy obtained Embedding free energy of C in ferrite The embedding free energy of C in ferrite at T = 0 K has been computed by two ab initio methods: the afore mentioned supercell method, and the CE method. At least two supercell results are available in the literature, h 0,Nb a = 0.74 eV = 71 kJ/mol [14] (Table III) and h 0,C a = 0.81 eV = 78 kJ/mol [44] (Table IV). Our own calculations are essentially identical to the latter. Both references [14,44] find their values confirmed by thermodynamic calculations based on experimental data [9,48], but some caution is required. When the solubility limit of C in ferrite c lmt is measured, as reported in Table I of Hasebe et al. [9], and the embedding free energy is extracted according to DG a emb ½C ¼ k B T lnðc lmt =3Þ; one obtains DH a emb ½C = 91 kJ/ mol and DS a emb ½C = 2.9 k B . Other solubility limit measurements give comparable values: DH a emb ½C = 101 kJ/mol and DS a emb ½C = 4.0 k B [45] and DH a emb ½C = 99 kJ/mol and DS a emb ½C = 3.8 k B [58]. As the embedding free energy computed with supercells for Nb was deemed unreliable for for T = 0 K, here also, a comparison with another method is called for. A CE can be carried out over the octahedral interstitial sublattice of bcc Fe. Of course, the interstitial sublattice can be vacant, or occupied by C. The highest possible C concentration occurs for FeC 3 when all interstices are occupied by C. We select a CE from a pool of clusters which satisfy two criteria: (a) up to four sites in a cluster, (b) no two sites within a cluster are farther apart than the fifth nearest neighbor. These criteria yield a pool of 26 clusters. Furthermore we compute the energies of all Fe n C m structures that can be generated by imposing that the primitive translations are 1.5 a bcc or shorter, this yields eight distinct types of periodic cells, and 389 distinct structures. After fully structurally relaxing these structures, we eliminate structures with Fe lattices that can no longer be classified as being bcc-based. By imposing that the CEs are complete [50], it is possible to examine all possible CEs that can be generated with our pool of clusters, we select the CE that has the smallest ''predictive error'' [52] or better known as leave-one-out cross-validation score [65]. We use this CE to compute the enthalpy of mixing, see Fig. 8, as described elsewhere [55]. Taking the slope of the mixing enthalpy at zero concentration C yields h 0,Nb a = 90 kJ/mol.
The magnetic contribution to the C embedding free energy can be estimated with Eq. 27 if the change in the Curie temperature per fraction C is known. Here we resort to the estimate by Hasebe et al. [9], DT C = -500 K, so that r C = -500/1043 = -0.48. It should be noted that the magnetization per atom decays more rapidly, q(M/M 0 )/ qc C = -1.37 in the limit of vanishing carbon concentration c C . Proceeding with r C = -0.48, we find that at T = 1000 K the embedding free energy is increased by 2 Â ðÀ0:48Þ Â ðÀ1:6Þ % 1.5 kJ/mol.
The vibrational contribution to the C embedding energy, as derived from the data shown in Fig. 6a, is displayed in Fig. 6b. C atoms leave the stiff diamond lattice when they are inserted in the Fe lattice and therefore clearly lower energy vibrational modes are created so that the vibrational contribution to the embedding free energy is strongly negative. In the temperature range of interest, from 500 to 1100 K, the vibrational free energy is approximately linear with T: DG a emb ½C = -3.068-0.0285 T kJ/mol. At T = 1000 K, the vibrational contribution is seen to be almost -32 kJ/mol. Electronic contributions to the embedding free energy are obtained from Eq. 28. In Fe 128 C supercells nð F Þ is about 10 states/eV greater than in the appropriate sum of Fe 128 and C diamond cells. An analysis of nð F Þ of the almost 400 structures for the CE gave an increase also in nð F Þ but just about 0.5 states/eV per dissolved C atom. If we accept the CE results we obtain thus DG elec ¼ À6 Â 10 À7 T 2 kJ/mol C dissolved in ferrite. The magnetic disordering and electronic contributions are small relative to the ab initio computed embedding enthalpy at T = 0 K, but the vibrational contribution is somewhat more significant. Figure 10 shows the total embedding free energy obtained from summing the T = 0 K enthalpy, the magnetic, vibrational, and electronic free energies (Fig. 9).

Free energy of NbC formation
The formation of NbC from the reference states gives a formation enthalpy of -104.3 kJ/mol NbC at T = 0 K, see Table 1. There is no magnetic free-energy contribution, so we focus on the vibrational and electronic terms. The vibrational free energy is favorable for NbC formation. The reason is that although NbC is hard and has high energy vibrational modes associated with the C species, it is still is not nearly as stiff as diamond so that the average of Nb-bcc and C-diamond has higher energy vibrational modes than NbC. In the temperature range of interest, from 500 to 1100 K, we may well represent the vibrational free energy by a linear function DG vib ½NbC = -1.706-0.0024 T kJ/mol NbC.
The electronic contribution is easily evaluated from the data in Table 1, where it is evident that NbC formation leads to the reduction of 0.86 states/(eV formula unit) at the Fermi level. Using Eq. 28 we find DG elec ½NbC ¼ 1 Â 10 À6 T 2 kJ/mol, which is a rather minor effect. The ab initio formation free energy can thus be given as DG f ½NbC = -106.0-0.0024 T ? 1 Â 10 À6 T 2 kJ/mol.

Ab initio solubility product
All parameters for computing the solubility product of stoichiometric NbC in ferrite have now been evaluated. Although the magnetic contribution to the free energy of NbC formation from the ferrite solid solution is highly nonlinear, the free energy of NbC formation from the pure elements and the Nb and C embedding free energies are all rather well-described by simple linear expressions in the temperature over the temperature range of interest for ferrite. In Fig. 11 the computed value of logðK a NbC Þ vs (1/T) is displayed. A linear fit with respect to 1/T is made both for the data below T C and for the data above T C , so that FM and PM A and B coefficients are found, which are listed in Table 3. The comparison with literature values is well within the noise of the experimental data [11,43,59,62]. The agreement is especially good with Hudd et al. [11], Turkdogan [62], and Taylor [59], with the ab initio solubility products generally being just a little lower, which in view of size effects and concomitant surface energies would be reasonable. The solubility product of Pichler et al. [43] clearly deviates from the other data for reasons that are not yet apparent. The good agreement is quite surprising in view of the many approximations that were made in the ab initio calculations: vacancies on the carbon sublattice in NbC were ignored, anharmonic effects in solid solution and precipitate phase were omitted; magnetic disordering was considered in a very crude manner as well; and supercell calculations for impurity properties were shown to be rather inconclusive. One now can wonder what the effect is of the various excitations on the computed solubility product: basically the electronic contribution could have been ignored, the magnetic contribution plays a minor role near and above the Curie temperature only, leaving the vibrational excitations as main effect. However, if all excitations are neglected, the A parameter increases by 1.7 %, while the B parameter decreases by 25 % over  Table 3. As could be expected the main contribution of the excitations is on the entropy and therefore on the B parameter. Ignoring the vibrational excitation particularly would have noticeably worsened the agreement with the experimental data.

Austenite solubility product
By using the experimentally determined values of DG cÀa emb ½Nb and DG cÀa emb ½C, we can now also obtainK c NbC with Eq. 19. Kaufman and Nesor [16] gave DG cÀa emb ½Nb = 210-3.556 T (J/mol) based on a Fe-Nb phase diagram assessment, while Sharma and Kirkaldy [46] as cited in Ref. [47] gave DG cÀa emb ½C = -64111.4-32.158 T (J/mol). Thus, the computedK c NbC is compared with various experimentally determined solubility products in Table 3. In general, a good agreement is seen, especially considering the scatter in the experimental data.

Conclusions
Solubility products are of great practical importance because they provide a quick insight on the likelihood on precipitation and dissolution. This is especially evident when alloys are dilute and when many alloying species are present because alternatives, such as a full multicomponent thermodynamic modeling, can be difficult and cumbersome. Nevertheless, ab initio modeling of solubility products is still in its infancy, especially so for steel. Here a first effort is made at predicting the solubility product of NbC in ferrite. A much simplified thermodynamic approach is chosen in order to get a transparent physical understanding of the physical effects that are essential for a reliable result. As a first step, ab initio T = 0 K formation enthalpies were used in conjunction with magnetic, vibrational, and electronic free energies. Many approximations were made in the ab initio calculations: vacancies on the carbon sublattice in NbC were ignored, anharmonic effects in solid solution and precipitate phase were omitted; magnetic disordering was considered in a very crude manner as well; and coupling between the various excitations was ignored. Thus the effect of lattice expansion was ignored, and also the effect of lattice expansion of the electronic density of states and the electronic free energy. Likewise, the effects of C deficiency in NbC on vibrational, electronic excitations, and on configurational entropy were ignored.
Moreover supercell calculations for impurity properties were shown to be rather inconclusive; giving an incorrect sign for the embedding enthalpy of Nb in ferrite, and quite probably a too low value for the embedding enthalpy of C in ferrite as well. It was shown that the n-atom supercell results exhibited n -1 scaling in the dilute limit only when similar supercell types were used. Supercells with sc, bcc, or fcc arrangements of impurities showed distinct convergence behavior toward the dilute limit. Constant volume supercells exhibited less scatter in the energetic and magnetic data than constant pressure supercells, presumably because supercell volume and shape optimization introduces non-systematic errors. Shortcomings of supercell calculations could be overcome by using special quasirandom structures or cluster expansions and it was surmised that this was due to the fact the SQSs results mimic the magnetic disorder on the Fe atoms that occurs under actual experimental conditions. This implies that the separation of the embedding free energy in terms of separate chemical and magnetic contributions is more complicated the common prescription followed here.
The fact that the ab initio computed solubility product agreed closely with the average of the experimental values gives confidence that more subtle effects, such as how other alloying elements affect the solubility product of NbC, in the spirit of Refs. [21,47], eventually can be modeled also.