A new algorithm to calculate binary phase diagrams

Computational Materials Science 159 (2019) 478–483

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A new algorithm to calculate binary phase diagrams a

a,⁎

b,⁎

Taibai Fu , Zhoushun Zheng , Yong Du Shuhong Liub a b

b

a

b

b

, Jiong Wang , Changfa Du , Bo Jin , Yuling Liu ,

T

School of Mathematics and Statistics, Central South University, Changsha, Hunan 410083, China State Key Laboratory of Powder Metallurgy, Central South University, Changsha, Hunan 410083, China

ARTICLE INFO

ABSTRACT

Keywords: Computational thermodynamics Equilibrium calculations Phase diagrams Linear system of equations Algorithm

A new algorithm to calculate binary phase diagrams is described. The condition numbers of phase matrices are analyzed, which can explain the stability during iterations. The linear systems of equations from thermodynamic calculation are simplified. Such simplifications significantly reduce the condition numbers of phase matrices and improve the computational efficiency. The Cu-Mg and Al-Ta phase diagrams are computed by using a homemade code to verify the efficiency of this algorithm over the method utilizing direct matrix inversion without simplifications. It is highly expected that the presently developed algorithm can be generalized to calculate multi-component phase diagrams efficiently.

1. Introduction Phase diagrams serve as the road map for materials design and process optimization by manipulating temperature, composition, pressure or other state variables to obtain desired microstructures including phase(s) and phase assembly [1,2]. Since experimental determinations of phase diagrams, in particular for multi-component systems, are costly and extremely time-consuming, CALPHAD (CALculation of PHAse Diagrams) approach was developed by Kaufman and Bernstein [3], who followed the pioneering work on phase diagram calculations by Van Laar in 1908 [4] and Meijering in 1950 [5]. Nowadays, a number of phase diagram calculation software packages are available, including Lukas program [6], ChemSage [7], MatCalc [8], ThermoCalc [9], Pandat [10], MTDATA [11], CaTCalc [12], FactSage [13], PyCalphad [14] and OpenCalphad [15]. Most of these software packages are for commercial purposes, and thus the detailed algorithms for phase diagram calculations are not presented in correlative papers [8–11,13]. More than one individual equilibrium calculation are needed in some softwares [6,7,12]. Most softwares cannot always guarantee that the computed equilibrium corresponds to the global minimum. The algorithms for equilibrium calculation are mostly based on the minimization of Gibbs energy, that was originary developed by White et al. [16]. The algorithm described in detail [17] was mainly derived by Eriksson and Rosen [18] as well as Hillert [19]. The calculation of the system of equations during the minimization was divided into two steps by Eriksson [20]. The implementation of the algorithm from Refs. [18–20] was explained in detail by Sundman et al. [17]. ⁎

All softwares calculate phase equilibria by using the iteration method. Different iterative schemes and the initial estimates can lead to problems such as the divergence, high computational cost and so on. Accordingly, some improvements have been developed by different groups [21–23]. For the well described algorithm [17] among the software packages [9,15], there exist linear systems derived from combining the Lagrangian multiplier method and Taylor expansion in the case of sublattice model, which is the most widely used thermodynamic model among various softwares. As the numbers of sublattices and constituents increase, the computational cost of inverting the coefficient matrices of the linear systems becomes larger and larger. Then the computational performance can be impeded significantly. Moreover, the condition numbers of coefficient matrices, which is used to express the sensitivity of matrix calculation to errors, become huge in general with the increase of numbers of sublattices and constituents. Numerical errors arise in initial estimates and solving equations, which may result in complex number when calculating the natural logarithms of constituent fractions in the iteration procedure. In the present work, we attempt to propose a new algorithm to calculate binary phase diagram for the sake of improving the existing algorithms and explain the divergence arising in iterative process from a condition number point of view. The new method greatly decreases the condition numbers of coefficient matrices and improves the efficiency of solving the linear systems of equations from equilibrium calculations. And we implement this algorithm by writing a home-made code and compare it with the widely used algorithm [9,17,19] in the field of computational thermodynamics.

Corresponding authors. E-mail addresses: [email protected] (Z. Zheng), [email protected] (Y. Du).

https://doi.org/10.1016/j.commatsci.2018.12.036 Received 1 August 2018; Received in revised form 16 November 2018; Accepted 17 December 2018 0927-0256/ © 2018 Elsevier B.V. All rights reserved.

Computational Materials Science 159 (2019) 478–483

T. Fu et al.

2. A widely used algorithm with inverting phase matrix

m1+ j=1

In thermodynamics, the Gibbs energy of each phase is modeled as a function of the temperature T, pressure P and constituent site fractions y or more fraction x. In multi-component system, for phase with K sublattices, let GM denote the Gibbs energy per mole of formula unit and n denote the number of moles of formula units. And let yik be the i-th constituent fraction in the sublattice k and y = (yik )i = 1, , mk ;k = 1, , K be all constituent fractions in phase . Using = {A, B, } be the set of components, then the number of moles of component A per mole forK m mula unit of the phase can be expressed as M A = k = 1 ak i =k1 bAi yik , where ak denotes the number of sites for sublattice k, and bAi is the stoichiometric coefficient of component A in the i-th constituent. For example, in a gas phase with H and O modeled (H2,O2,H2O), the stoichiometric coefficients of component H in 3 constituents are bH 1 = 2, bH 2 = 0 and bH3 = 2 , and the analogues of component O are bO1 = 0, bO2 = 2 and bO3 = 1, respectively [17]. In a system, the total n MA . number of moles of component A is formulated as NA = The equilibrium of a closed system under fixed T , P and overall composition will be defined by the following Gibbs energy minimization:

min G =

n GM (T , P , y )

yik T

k = 1,

, K;

NA =

n M A,

M A MB

k=1

2G M y12 2G M

k

(1

yik ) + i=1

µA (NA

n MA)

M A µA = 0

(5)

Thus, µA is the chemical potential. And for the partial derivatives of , we obtain L with respect to each constituent fraction of phase

µA A

MA yik

k

= 0, i = 1,

, mk ; k = 1,

,K (6)

derivative

GM yik

at

the

point

(T0, P0, y0 ) ,

GM,0 = GM (T0, P0, y0 ) , we have the first-order difference form for follows: 2G ,0 2G ,0 GM GM,0 M M = + T+ P+ yik yik yik T yik P

K

mt

t=1 j =1

2G ,0 M

yik yjt

yjt

letting GM yik

k

n

, mk ; k = 1,

,K

(8)

(9)

1

y1 y2

2G M y22

1

1

1

0

y1 y2 =

µA

MA y1

+ µB

MB y1

GM y1

µA

MA y2

+ µB

MB y2

GM y2

0

n

(10)

1

cond (M ) M 1 M

M + M

b b

(11)

0, b 0, M , b and x are the corrections of M , b , x , where |M| respectively. cond (M ) = M M 1 is the condition number of matrix M. The expression (11) means that when the cond (M ) is huge, the relative error of solution vector x caused by the corrections of coefficient matrix and right-hand-side vector may be a big number. The Cu-Mg binary system is used as an example to observe the condition number of the phase matrix. At T = 1000 K , the constitution fractions of phase Cu2Mg, which is described with the sublattice model (Cu,Mg)2(Cu,Mg)1, are given by yCu = 0.990332, yMg = 0.009667 and yCu = 0.000368, yMg = 0.999631, where and mean the first and second sublattices, respectively. And these data are sufficiently close to the equilibrium at current temperature.

as

(7)

The last term in Eq. (7), which is a summation over all constituents j on all sublattices t, is rewritten as just a summation over the constituents j. Inserting Eq. (7) into Eq. (6) and replacing the differentials with finite differences, we obtain the following linear equations about variables yj and

2G M y1 y2

x x

According to the first-order Taylor approximate expansion of the

partial

i = 1,

Then we get a set of new fraction constitutions with y = y + y . By repeating this process, we can obtain more reliable approximation results to equilibrium. The method of inverting the coefficient matrix directly is commonly used to solve the above linear system of equations, and the operations are obviously controlled by the size of the phase matrix. The larger the numbers of the sublattices and constituents are, the higher the size of the phase matrix is. As a result, the computational cost of inverting phase matrix will be increased greatly with the increase for the numbers of the sublattices and constituents. Besides, the accumulation of errors caused by iterative calculations can lead to the possibility of non-convergence or complex solutions. Now we analyze the phase matrix by using the concept of condition number. Taking into account the error estimation expression of the solution to the equations Mx = b , we get the following inequality:

A

A

n

yik

1

where µA and k are the Lagrange multipliers for all phases. Taking the partial derivatives of L with respect to the amount of moles formula unit of each stable phase , we have

GM yik

P,

yik P

GM,0

yik

It is clear that we can solve these two equations for the approximations of chemical potentials µA , µB directly. Thus, the fraction corrections yi can be calculated by inserting µA , µB into Eq. (8). In fact, omitting the phase indices, Eqs. (8) and (2) can be written for each phase:

(4)

L =n yik

2G ,0 M

M A,0

mk

n GM +

L = GM n

µA

A

GM µA µB = G M

MA MB

Here is the set of stable phases. With above constraints, we can obtain the Lagrange formula as follows [17]:

L=

=

For a binary system with fixed T , P and two stable phases , , for example, after providing the initial approximations of the fraction constitutions y , y , there are only two unknowns µA , µB and two equations in Eq. (5).

(3)

K

T

k

n

3. A new algorithm for binary phase diagram calculations

(2)

A

yj

The coefficient matrix of linear equations which are composed of Eqs. (2) and (8) is called the phase matrix [17]. We note that, with fixed T and P, the fraction corrections yj in Eq. (8) just depend on the chemical potentials µA . In [17], before calculating the chemical potentials µA , inverting phase matrices of stable phases must be done firstly. Indeed, under conditions as described below, we can obtain the chemical potentials directly without the inverted phase matrix.

mk

yik = 1,

yik y j 2G ,0 M

(1)

i=1

2G ,0 M

+ mK

: 479

Computational Materials Science 159 (2019) 478–483

T. Fu et al.

According to Eq. (10), the phase matrix becomes

5596 2199 39317 2199 573300 16872 39317 16872 7526272 58179 37278 4337 1 1 0 0 0 1

Mp =

58179 37278 4337 2772 0 1

1 1 0 0 0 0

reduced from 6 × 6 to 2 × 2 . Similarly, in the Cu-Mg binary system, for phase Cu2Mg, the simplified phase matrix is calculated as:

0 0 1 1 0 0

(

574497 1543 M¯p = 1543 7520370

The presently obtained condition number of M¯p is cond (M¯p ) = 13.090355, which is significantly less than the condition number cond (Mp) = 4.655703 × 1011 associated with the widely used algorithm [19]. That is to say that the relative error of the solution associated with the simplified linear equations is less influenced by the accumulation errors. Now, we calculate the Cu-Mg [24] and Al-Ta [25] binary phase diagrams by writing a home-made MATLAB code. Fig. 1 demonstrates the flow chart of the new algorithm to calculate binary phase diagrams. A detailed description for the procedure to calculate binary phase diagrams is as follows.

The condition number of Mp is calculated to be cond (MP ) = 4.655703 × 1011. This means that the phase matrix has a large condition number, which probably brings inaccurate fraction corrections yi . Accordingly, the accumulation errors can possibly result in unreasonable constituent (0, 1) during iterations. Obviously, the natural logafractions yi rithms of these constituent fractions are complex numbers. In this case, one cannot compute the whole phase diagram without some special handling, as this happens for many commercial software packages. In the binary system, notice that the elements of the last K rows and the last K columns are either 1 or 0. In view of the special structure of the phase matrix, Eq. (10) can be simplified by using the techniques for solving linear equations. For example, for a phase with the sublattice model (A, B)a1(A, B)a2 , we have the system of equations solved for fraction corrections yi as 2G M y12 2G M

y1 y2

2G M y1 y2

2G M y1 y3

2G M y1 y4

1 0

2G M y22 2G M

2G M y2 y3

2G M y2 y4

1 0

2G M y32 2G M

2G M y3 y4

0 1

y1 y2 y3 y4

2G M y1 y3

y2 y3

2G M y1 y4

2G M y2 y4

y3 y4

2G M y42

0 1

n

1 0

1 0

0 1

0 1

0 0 0 0

n

µA

MA y1

+ µB

MB y1

GM y1

µA

MA y2

+ µB

MB y2

GM y2

= µA

MA y3

+ µB

MB y3

GM y3

µA

MA y4

+ µB

MB y4

GM y4

1

Step 1. Provide a temperature T0 . Calculate the Gibbs energy values for all phases at different compositions, then get a set of discrete points and the convex hull. Find the first line segment from left to right and make sure that the ends of this line belong to two different phases. Set these two phases as the stable phases , for the starting of calculation. For each phase, based on the minimum of the distances from different discrete points to the line, select the initial values of constituent fractions. In general, the start temperature T0 is a relatively high temperature at which the equilibrium calculations have a low demand on the initial approximations of constituent fractions. If the initial values cannot be found at current temperature, provide an alternative temperature lower than T0 and repeat this step. Step 2. Calculate the equilibrium involving the chemical potentials, constituent fractions of stable phases as well as the constituent fractions and the driving forces of the metastable phases. If the driving force of some phase is positive, two new stable phases can be determined by the constituent fractions and the corresponding Gibbs energy values of phases , , then repeat this step. Step 3. Store these data involving start temperature T0 , stable and metastable phases and their constituent fractions as two “Start Data” with two temperature steps opposite in sign. Select one of “Start Data”. Step 4. Take T = T0 + T and calculate the equilibrium iteratively. If current temperature T is below 298.15, select the next one in “Start Data” and repeat this step. If two compositions of two stable phases are both close to 0 or 1, select the next one in “Start Data”, and repeat this step. If two compositions of two stable phases are both close to one neither 0 nor 1, go to step 5. If some driving force is positive, add the phase to the set of stable phases. Two new pair of phases are stored into “Start Data” with two temperature steps judged by the values of three compositions. Go to step 6. Otherwise, let T0 = T and repeat this step. For each T, it is necessary to obtain the initial values of constituent fractions before the iterative calculation. This process can be implemented by computing the discrete points and convex hull, which requires a high computational cost in the whole calculation of phase diagram. Hence in the calculation of two-phase region, the iteration results at previous temperature T0 can serve as the initial values at the current temperature T, which can improve the efficiency greatly. T , calculate the Gibbs energy values of these two Step 5. Let T0 = T phases at different compositions, find the other common tangent and calculate constituent fractions of all phases. Store T0 , T , stable and metastable phases and their temperature steps constituent fractions into “Start Data”. On the other hand, take

0 0

2

(12) Eliminating the terms , we then obtain the following equation. i

n

2G M y12 2G M

y1 y3

2G M

2G M

y1 y2

y1 y2

2G M

2G M

y1 y4

y2 y3

1 0

2G M y22 2G M

2G M y1 y3

2G M y2 y3

2G M y1 y4

2G M y2 y4

2G M y32

2G M

2G M

y3 y4

y3 y4

2G M y42

y2 y4

1 0

( = µ (

0 1

µA

MA y1

MA y2

A

MA y3

MA y4

)+µ ( )+µ (

y1 y2 y3 y4

0 1

B

MB y1

MB y2

B

MB y3

MB y4

) ( ) (

GM y1

GM y2

GM y3

GM y4

) )

0 0 (13) Note that the last two equations in Eq. (13) take the equivalent y1 , y4 = y3 . Hence, inserting these relations into the forms y2 = first two equations in Eq. (13) leads to the following system of equations:

Mp

µA y1 = y3 µA

( (

MA y1

MA y2

MA y3

MA y4

)+µ ( )+µ ( B

MB y1

MB y2

B

MB y3

MB y4

) ( ) (

GM y1

GM y2

GM y3

GM y4

) )

(14)

where Mp is defined by

M¯p =

2G M y12 2G M y1 y3

2 2G M y1 y4

2G M y1 y2

+

2G M y2 y3

2G M y22

+

2G M y1 y3 2G M y2 y4

2G M y2 y3 2G M y32

2

2G M y1 y4 2G M y3 y4

+

+

)

2G M y2 y4

2G M y42

(15) It is obvious that only addition and subtraction are performed during the process of this simplification. The size of the phase matrix is 480

Computational Materials Science 159 (2019) 478–483

T. Fu et al.

GM T, T

GM = M A µA + MB µB

={ , , }

(16)

Solve this equation for µA , µB and T , then yi can be obtained by inserting µA , µB and T in the next equation of each phase. m1+ j=1

+ mK

2G M

yik yj

yj 2G M

yik T

k

n

= T,

A

µA

i = 1,

MA

GM

yik

yik

, mk ; k = 1,

,K

(17)

And the above equations can be simplified as described in the present algorithm. In this way, the Cu-Mg and Al-Ta phase diagrams can be calculated automatically. The setting of the initial composition is not needed, and the numbers of moles of stable phases are not calculated. In our procedure, we focus on two-phase and three-phase equilibria, and single phase region comes automatically. Given a set of thermodynamic parameters for all the phases, this novel algorithm automatically calculates the stable phase diagram without either prior knowledge of the diagram or special user skills. In the calculation of Cu-Mg and Al-Ta phase diagrams, the step size of temperature T is 1. Additionally, we compare the running time of calculating binary phase diagrams by using the iteration Eq. (10) [9,17,19] for fraction corrections with that from presently proposed iteration Eq. (14). In particular, as the operating efficiency of different softwares varies, all computations for the Cu-Mg and Al-Ta phase diagrams are implemented by MATLAB. There is no sense to compare the running time of this MATLAB code with that of other softwares since different softwares use different computer language and are implemented in different operation systems. On the other hand, no detailed descriptions about almost all commercial softwares are presented. In this sense, the presently developed algorithm is described in detail for the sake of filling in this gap. By using the iteration Eq. (10) for fraction corrections y reported in the literature, the calculation of Fig. 2(a) runs generally more than 38 s, and the calculation of Fig. 3(a) spends about 240 s. In Fig. 2(a), the iterative results for phase Cu2Mg are non-real numbers at T1 (K )(400 < T1 < 500) and T2 (K )(300 < T2 < 400) . Consequently, the phase diagram is incomplete in the lower-left corner. And in Fig. 3(a), the constituent fractions for phase AlTa2 are also non-real at T3 (K )(400 < T3 < 500) . One then cannot get the phase diagram in the lower-right corner. On the contrary, we obtain the complete Cu-Mg and Al-Ta phase diagrams with the iteration Eq. (14) proposed in the present work. Moreover, the running time for Fig. 2(b) is less than 5 s usually, and the analogue of Fig. 3(b) is about 24 s. The unique features of the new algorithm are that the calculation of

Fig. 1. Flow chart for the procedure to calculate binary phase diagrams.

T = T0 + T , calculate the equilibrium and let T0 = T . Repeat this process until two compositions of two stable phases are close to each other. Select the next one in “Start Data” and go back to step 4. Step 6. With current T and constituent fractions of three stable phases as initial values, calculate three-phase equilibrium. Select the next one in “Start Data” and go back to step 4. The iterative equations for calculating three-phase equilibrium with variable T are as follows:

Fig. 2. Cu-Mg phase diagram computed by the present algorithm via MATLAB. (a) y are calculated by using Eq. (10). (b) y are calculated by using Eq. (14). 481

Computational Materials Science 159 (2019) 478–483

T. Fu et al.

Fig. 3. Al-Ta phase diagram computed by the present algorithm via MATLAB. (a) y are calculated by using Eq. (10). (b) y are calculated by using Eq. (14).

chemical potentials µA is independent of the inverse of the phase matrix and that the phase matrix can be simplified with only addition and subtraction, which decreases the computational expense of calculating fraction corrections greatly. Since the structure for the phase matrices is the same for binary and multi-component systems, only the dimension is changed when the present algorithm is extended from binary systems into multi-component ones. Generally, in m-component system, assuming m constituents in each sublattice, the size of the phase matrix can be reduced from (m + 1) K × (m + 1) K to (m 1) K × (m 1) K with K being the number of sublattices according to the present algorithm.

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.commatsci.2018.12.036. References [1] Y.A. Chang, S. Chen, F. Zhang, X. Yan, F. Xie, R. Schmid-Fetzer, W.A. Oates, Phase diagram calculation: past, present and future, Prog. Mater. Sci. 49 (3) (2004) 313–345. [2] R. Schmid-Fetzer, Phase diagrams: the beginning of wisdom, J. Phase Equilib. Diff. 35 (6) (2014) 735–760. [3] L. Kaufman, H. Bernstein, Computer Calculation of Phase Diagrams, Academic Press, New York, 1970. [4] J.J. Van Laar, Melting or solidification curves in binary system, Z. Phys. Chem. 63 (1908) 216–253. [5] J.L. Meijering, Segregation in regular ternary solutions: Part I, Philips Res. Rep. 5 (1950) 333–356. [6] H.L. Lukas, J. Weiss, E.T. Henig, Strategies for the calculation of phase diagrams, Calphad 6 (3) (1982) 229–251. [7] G. Eriksson, K. Hack, ChemSage–a computer program for the calculation of complex chemical equilibria, Metall. Mater. Trans. B 21 (6) (1990) 1013–1023. [8] E. Kozeschnik, B. Buchmayr, MatCalc – A Simulation Tool for Multicomponent Thermodynamics Diffusion and Phase Transformations, in: H. Cerjak, H.K.D.H. Bhadeshia (Eds.), Mathematical Modelling of Weld Phenomena 5 book 738, IOM Communications, London, 2001, pp. 349–361. [9] J.O. Andersson, T. Helander, L. Höglund, P. Shi, B. Sundman, Thermo-Calc & DICTRA, computational tools for materials science, Calphad 26 (2) (2002) 273–312. [10] S.L. Chen, S. Daniel, F. Zhang, Y.A. Chang, X.Y. Yan, F.Y. Xie, R. Schmid-Fetzer, W.A. Oates, The PANDAT software package and its applications, Calphad 26 (2) (2002) 175–188. [11] R.H. Davies, A.T. Dinsdale, J.A. Gisby, J.A.J. Robinson, S.M. Martin, MTDATAthermodynamic and phase equilibrium software from the national physical laboratory, Calphad 26 (2) (2002) 229–271. [12] K. Shobu, T. Tabaru, Development of new equilibrium calculation software: CaTCalc, Mater. Trans. 46 (6) (2005) 1175–1179. [13] C.W. Bale, E. Bélisle, P. Chartrand, S.A. Decterov, G. Eriksson, K. Hack, I.H. Jung, Y.B. Kang, J. Melançon, A.D. Pelton, C. Robelin, S. Petersen, FactSage thermochemical software and databases–recent developments, Calphad 33 (2) (2009) 295–311. [14] R. Otis, Z.K. Liu, pycalphad: CALPHAD-based Computational Thermodynamics in Python, J. Open Res. Softw. 5 (1) (2017) 1–11. [15] B. Sundman, U.R. Kattner, C. Sigli, M. Stratmann, R. Le Tellier, M. Palumbo, S.G. Fries, The OpenCalphad thermodynamic software interface, Comp. Mater. Sci. 125 (2016) 188–196. [16] W.B. White, S.M. Johnson, G.B. Dantzig, Chemical equilibrium in complex mixtures, J. Chem. Phys. 28 (5) (1958) 751–755. [17] B. Sundman, X.G. Lu, H. Ohtani, The implementation of an algorithm to calculate thermodynamic equilibria for multi-component systems with non-ideal phases in a free software, Comp. Mater. Sci. 101 (2015) 127–137. [18] G. Eriksson, E. Rosen, Thermodynamic studies of high temperature equilibrium. VIII. General equations for the calculation of equilibria in multi-phase systems, Chem. Scri. 19 (1973) 193–194. [19] M. Hillert, Some viewpoints on the use of a computer for calculating phase diagrams, Physica B+ C 103 (1) (1981) 31–40. [20] G. Eriksson, Thermodynamic studies of high temperature equilibria–XII:

4. Conclusions In this work, the condition numbers of the phase matrices are analyzed for the stability in the process of iteration. As the size and the condition number of the simplified phase matrix are decreased sharply, the accumulation of numerical errors arising in computing the constituent fractions is reduced accordingly. A home-made code has been written and used to compute the Cu-Mg and Al-Ta binary phase diagrams. It is demonstrated that the newly developed algorithm can establish the whole binary phase diagrams much faster and more stable than the method inverting the phase matrices directly without simplifications. Moreover, the present algorithm can be generalized to the calculation of multi-component phase diagrams highly efficiently. CRediT authorship contribution statement Taibai Fu: Conceptualization, Data curation, Formal analysis, Methodology, Software, Validation, Writing - original draft. Zhoushun Zheng: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation. Yong Du: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Writing - review & editing. Jiong Wang: Project administration, Supervision, Validation. Changfa Du: Data curation, Formal analysis, Validation. Bo Jin: Data curation, Investigation, Methodology, Formal analysis, Validation. Yuling Liu: Data curation, Validation. Shuhong Liu: Data curation, Validation. Acknowledgement This work was supported by the National Key Research and Development Program of China [Grant Nos. 2017YFB0701700, 2017YFB0305601].

482

Computational Materials Science 159 (2019) 478–483

T. Fu et al. SOLGASMIX, a computer program for calculation of equilibrium compositions in multiphase systems, Chem. Scri. 8 (1975) 100–103. [21] M. Emelianenko, Z.-K. Liu, Q. Du, A new algorithm for the automation of phase diagram calculation, Comp. Mater. Sci. 35 (1) (2006) 61–74. [22] M.H.A. Piro, S. Simunovic, T.M. Besmann, B.J. Lewis, W.T. Thompson, The thermochemistry library Thermochimica, Comp. Mater. Sci. 67 (2013) 266–272. [23] J. Snider, I. Griva, X. Sun, M. Emelianenko, Set based framework for Gibbs energy

minimization, Calphad 48 (2015) 18–26. [24] P. Liang, H.J. Seifert, H.L. Lukas, G. Ghosh, G. Effenberg, F. Aldinger, Thermodynamic modelling of the Cu-Mg-Zn ternary system, Calphad 22 (4) (1998) 527–544. [25] Y. Du, R. Schmid-Fetzer, Thermodynamic modeling of the Al-Ta system, J. Phase Equilib. 17 (4) (1996) 311–324.

483