Implementation of the UNIQUAC model in the OpenCalphad software

The UNIQUAC model is often used, for example in engineering, to obtain activity coefficients in multicomponent systems, while the CALPHAD method is known for its capability in phase stability assessment and equilibrium calculations. In this work, we combine them by representing the UNIQUAC model according to the CALPHAD method and implementing it in the OpenCalphad software. We explain the harmonization of nomenclature, the handling of the model parameters and the equations and partial derivatives needed for the implementation. The successful implementation is demonstrated with binary and multicomponent phase equilibrium calculations and comparisons with literature data. Additionally we show that the implementation of the UNIQUAC model in the OpenCalphad software allows for the calculation of various thermodynamic properties of the systems considered. The combination provides a convenient way to assess interaction parameters and calculate thermodynamic properties of phase equilibria. © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Generally, in chemical engineering two kinds of thermodynamic models are used: equations of state and models for the excess Gibbs free energy of mixing [1]. The models have in common that they contain parameters that have to be adjusted to represent experimental data. The models should be consistent with statistical mechanics [2] and thermodynamic constraints [3]. Typical earlier models are the Margules [4], Van Laar [5,6], Redlich-Kister [7], Scatchard-Hildebrand [8,9] and Flory-Huggins [10,11] equations. In 1964, G.M. Wilson introduced the local composition concept into the excess Gibbs energy model [12]. Some well-known models based on this concept are Wilson [12], Non-Random Two-Liquid (NRTL) [13], UNIversal QUAsiChemical (UNIQUAC) [14], and UNI-QUAC Functional group Activity Coefficients (UNIFAC) [15]. These are typically used for modeling Vapor-Liquid Equilibria (VLE), Liquid-Liquid-Equilibria (LLE) and Solid-Liquid-Equilibria (SLE) at relatively low pressures.
Although the mathematical expression of the UNIQUAC model is more complex than that of the NRTL model, UNIQUAC is used more than NRTL in the area of chemical engineering. The reason is that UNIQUAC has fewer adjustable parameters, two instead of three, which are less dependent on temperature, and can be applied to systems with larger size differences. A more detailed overview can be found, for example, in Polling et al. [16].
The UNIFAC model [15] was published simultaneously with UNIQUAC and is a group-contribution based equivalent of the UNIQUAC model. Importantly, UNIFAC is completely predictive, which makes it of great practical value in, for example, the chemical engineering community. The UNIFAC model is based on a growing databank of parameters that are obtained from experimental data. The development of UNIFAC over time is shown in Table 1, which shows that more and more phase equilibrium information and excess thermal properties are included, allowing the number of available functional groups to increase steadily.
In the 1970s, simultaneously with the development of the UNIQUAC and UNIFAC models, CALculation of PHAse Diagrams (CALPHAD)-type models [33,34] were also developed. Originally, they were mainly used in the fields of inorganic chemistry and metallurgy. Both of these model developments benefited from the availability of computers and have had a strong development during the past 50 years.
The CALPHAD technique is applied to alloys, ceramics and hightemperature processes involving binary, ternary and higher-order systems. Such a system typically includes 10e100 different crystalline phases in addition to gas and liquid phases. Each phase is described with a different Gibbs energy model. CALPHAD includes Long Range Ordering (LRO) for the crystalline phases and Short Range Ordering (SRO) in both solids and liquids. In addition, different kinds of excess model parameters are used to describe the binary, ternary and higher order interaction between the components of a phase. Dealing with so many different phases, a central part of the CALPHAD modeling is a unary database [35] where the Gibbs energy for each pure element in each phase is described as a function of temperature, T, relative to reference conditions, up to high temperatures of several thousand kelvin. These Gibbs energy functions are known as lattice stabilities and may include parameters for magnetic ordering. They are needed because an element may dissolve in many different crystalline phases for which the element itself is not stable. Extensive descriptions of CALPHAD models can be found, for example, in Lukas et al. [36], Hillert [3] and Saunders and Miodowski [37]. CALPHAD software is able to calculate multicomponent equilibria with hundreds of possible phases and multicomponent phase diagrams. The database contains, in addition to the unary data, model-dependent parameters that depend on the amounts of two, three or more components.
There are several commercial companies marketing software and databases for UNIQUAC/UNIFAC and CALPHAD. OpenCalphad is an open source software [38]. CALPHAD applications typically concern the temperatures up to 2000 K whereas UNIQUAC is normally applied at temperatures up to around 400 K. This is one of the main reasons for the different approaches used in the thermodynamic models. However, it is well possible to handle both UNIQUAC and CALPHAD-type models within the same software.
In this work, UNIQUAC is successfully incorporated into Open-Calphad. We show that the combination allows for an easy and accurate calculation of thermodynamic properties and phase diagrams, and the determination of UNIQUAC interaction parameters.

Exploring the differences between UNIQUAC model and CALPHAD method
There are several significant differences in the use of thermodynamics between UNIQUAC and CALPHAD: the main difference is that the UNIQUAC model has a combinatorial entropy based on the theory of Guggeneim [2], taking into account that the components can have very different sizes. In CALPHAD models, on the other hand, the components, normally atoms, have similar sizes and thus an ideal configurational entropy is used. Moreover, CALPHAD can model many crystalline phases simultaneously and takes longrange ordering into account. The use of ideal mixing on each sublattice describes the configurational entropy of complex solid phases with reasonable accuracy. For liquids, there are several models used in CALPHAD to describe short range ordering [39,40], but they do not account for polymerization or polarization.
Another difference is that UNIQUAC models are based on excess Gibbs energy or activity coefficients, whereas the CALPHAD models are always based on a molar Gibbs energy expression; however, as the activity coefficients are calculated from the molar Gibbs energy, this difference is not very important. For applications, a major difference is that CALPHAD software and databases usually consider more than 100 different solid crystalline phases with different molar Gibbs energy models, whereas applications of the UNIQUAC model usually involve few solid phases.
In the CALPHAD method, all models are described using a molar Gibbs energy for each phase and the models can be quite different for different phases. One reason to avoid activity coefficients is that their modeling may cause inconsistencies between the molar Gibbs energy and the chemical potentials.
The integration of both types of models in a single software allows the calculation of equilibria between complex solids described with CALPHAD models and liquids and polymers described with the UNIQUAC model. The free OpenCalphad software [38] was selected because it is publicly available, the source code is open for developing new models and it is written in the new Fortran standard. The OpenCalphad software can perform multicomponent equilibrium calculations and provide all thermodynamic properties of interest, such as enthalpies, entropies, chemical potentials and heat capacities as well as phase diagrams.

The molar Gibbs energy and the molar excess Gibbs energy
In CALPHAD, the molar Gibbs energy for a phase a is expressed as eq. (1): where the summation over + G a i represents the contribution from the lattice stabilities explained in section 1. The term multiplied with RT is the ideal configurational entropy and E G m is the excess Gibbs energy. In CALPHAD, it is preferred to use the E as a presuperscript because the normal superscript position is reserved for the phase label as CALPHAD normally deals with many different phases. The excess Gibbs energy is described by regular solution parameters, L, for binary and higher order interactions which are constants or linearly dependent on the temperature.
The chemical potential of a component i is calculated from the total Gibbs energy for a system: where G is the total Gibbs energy and N i the number of moles of component i. The chemical potentials are calculated from a molar Gibbs energy by combining the first derivatives of the molar Gibbs energy: where G m ¼ G=N is the molar Gibbs energy using mole fractions x i as composition variables. Temperature, pressure and the other mole fractions are constant when calculating the partial derivatives. This equation is one of the basics of the CALPHAD method and is derived, for example, in Refs. [36,41]. For phases with sublattices, a slightly more complicated equation is needed, as derived in Ref. [42].
The chemical potential m i is expressed using activity coefficients as in all thermodynamics books: where + m i is the reference chemical potential, a i is the activity and g i is the activity coefficient describing the deviation from ideality. For an ideal solution, g i ¼ 1.
When the chemical potentials, m i , are known, the molar Gibbs energy for a phase a is obtained by a simple summation: The CALPHAD models are typically used for high-temperature systems compared with UNIQUAC model as introduced before and can describe many different types of crystalline phases in addition to liquid and gas phases. In CALPHAD the ideal configurational entropy is normally used but it includes long range ordering in the crystalline phases, separate modeling of the ferromagnetic transition and complex excess models. There is interest to combine calculations using CALPHAD data for calculations at ambient temperatures to study, for example, corrosion with water and other fluids. Most CALPHAD calculations involve systems with 8-10 components and there is a particular unary database that provides data for the pure elements in different crystalline states as well as in liquid and gas species. This makes it possible to combine and extend descriptions of binary and ternary assessments to multicomponent alloys.
In detail CALPHAD always models the integral Gibbs energy, as in eq. (1), expressed as a function of temperature and composition using model parameters. The chemical potentials and activity coefficients are derived from this numerically.
Unlike CALPHAD, the UNIQUAC model, developed by Abrams and Prausnitz [14], is an expression of the molar excess Gibbs energy, g E . In addition, the activity coefficient of a component i in the liquid, g i , is derived from the partial derivative of n T , g E with respect to n i and expressed as (all the symbols are the same as those in the original paper): The derivation of expressions for the activity coefficients starting from eq. (11) is usually tedious and prone to errors. An alternative is to derive the activity coefficients from eqs (5)e(7). An important difference is that excess Gibbs energy in the UNIQUAC model includes a combinatorial entropy whereas CALPHAD, for substitutional solutions, uses an ideal configurational entropy. This difference will be discussed in detail in section 3.
In CALPHAD models, the ideal configurational entropy is modified using sublattices for LRO in crystalline phases. For liquids with strong SRO, such as molten salts or ionic components, there are special models like the 2-sublattice ionic liquid model [39] or the quasichemical model [40]. In addition, several excess parameters depending on binary and ternary interactions are used in the excess Gibbs energy in eq. (3). Chemical potentials and activity coefficients are calculated by the software using eq. (5).
Based on the reasoning so far, it is useful to represent the UNI-QUAC model in the CALPHAD method and implement it in software based on the CALPHAD method, such as, OpenCalphad. Analytical expressions of the first partial derivatives of molar Gibbs energy with respect to the components have been implemented in the software. Analytical expressions for the second derivatives are used to speed up convergence [41], but have no effect on the final calculated results. In this work, we did not implement the analytical expressions for the second derivatives.

The UNIQUAC model using CALPHAD nomenclature
The UNIQUAC model, derived by Adams and Prausnitz [14], presents the expression for the excess Gibbs energy as the sum of combinatorial, cmb G m , and residual, res G m , contributions. They are combined into an expression for the molar Gibbs energy, see eqs from 12 to 18: where q i is a surface-area parameter and r i a volume parameter, which are both component structural parameters for constituent i. z is the average number of nearest neighbors of a constituent, always assumed to be 10.
In accordance with the separate excess Gibbs energy terms for the configurational and residual contributions, we calculate the chemical potentials and activity coefficients using eq. (5) for each term separately.
where g c i and g r i represent the combinatorial and the residual contributions to the activity coefficients, respectively.

The configurational Gibbs energy in the UNIQUAC model
Abrams and Prausnitz [14] derived the configurational activity coefficient from Guggenheim [2] taking into account that the components usually have different interaction surfaces, q i , and volumes, r i . We start from the configurational molar Gibbs energy, cfg G m , eq. (18), based on the UNIQUAC model including the ideal term: The auxiliary function f is defined to simplify the calculations below: where, by definition, P

The first derivative of the configurational Gibbs energy
The derivative of cfg G m =RT with respect to x k is: From this we obtain (see Supplementary Material): The section "Derivatives of the configurational part" in the Supplementary Material shows that inserting eq (18) and (24) in eq. (5) results in a chemical potential which reproduces the configurational activity coefficient from Ref. [14]. The chemical potentials can also be used to recover the molar Gibbs energy using eq. (8). The derivatives of cfg G m =RT with respect to T and P are zero.

The second derivative of the configurational Gibbs energy
The numerical procedure in OpenCalphad requires also the second derivatives with respect to the mole fractions: as shown in "Derivatives of the configurational part" in the Supplementary Material.

The residual part of the UNIQUAC model
The residual integral Gibbs energy expression from Abrams and Prausnitz [14] is denoted res G m =RT in the CALPHAD method and written as: where Du ji ¼ u ji À u ii and u ii is a property of the pure component. This means that, in the general case, Du ij sDu ji and also t ij st ji . We introduce a new symbol, w ji with dimension K to describe t ji : The first order partial derivatives of residual Gibbs energy to component concentration and temperature are: The second order partial derivatives are: More detailed mathematical derivations are shown in "Derivatives of the residual part" in the Supplementary Material.
In order to confirm the consistency between chemical potential and Gibbs energy in the UNIQUAC model, eq. (5) is rearranged as: The residual contribution to the chemical potential or the activity coefficient has been calculated by summing the first derivatives of res G m =RT according to eq. (5) and it is shown to be the same as the equation in Abrams' paper [14]. For more detailed mathematical derivations, see "Confirmation the activity coefficients are the same" in the Supplementary Material.

Equilibrium calculations
In chemical engineering, especially in the separation of multicomponent mixtures, the objective of phase equilibrium calculations is to obtain equilibrium compositions of different phases. The compositions and thermodynamic properties in multicomponent multiphase systems can be calculated by various methods. With a model that provides the chemical potentials, or the activity coefficients, of the components, the equilibrium is calculated by finding the composition of the phases that gives the same chemical potentials of each component in all stable phases. This is usually a rapid and stable method if the set of stable phases is known beforehand, and has been used for the UNIQUAC model. In the CAL-PHAD method, the calculations sometimes involve hundreds of phases, and determining the set of stable phases is a major problem.
CALPHAD uses a Gibbs Energy Minimization (GEM) technique [41,43,44], taking into account all possible phases and then determining the set of phases and their compositions that give the lowest Gibbs energy. At this minimum, the components have the same chemical potential in each stable phase.
The OpenCalphad software uses an algorithm described in Ref. [41] and requires a model for the molar Gibbs energy. The algorithm also requires that the first and second derivatives of the Gibbs energy with respect to temperature, pressure and all constituents are implemented in the software in order to find the equilibrium in a fast and efficient way.
However, the situation for the UNIQUAC model is totally different. Previous attempts to use the Gibbs minimization technique to calculate equilibria with the UNIQUAC model [45e47] have been limited to binary and ternary systems. Some have used linear programming techniques which rely on the possibility to calculate the chemical potentials in each stable phase separately [46,47]. This is in fact different from to the idea of a Gibbs energy minimization, in which the set of stable phases is not known a priori because the method is more or less identical to equating the chemical potentials of the components in a preselected set of stable phases. The algebraic method developed by Iglesias-Silva et al. [45] can calculate phase equilibrium for any number of components and phases and their eq. (4) is similar to eq. (5) in this paper. However, their mathematical method to calculate the equilibrium is not generalized to multiphase systems and they have not introduced the entropy of mixing in their equations.
Talley et al. [48] compared UNIQUAC and CALPHAD modeling techniques. They claim significantly better fit to experimental data  Fig. 4 in Abrams and Prausnitz [14]. (b): Calculated miscibility gap vs. temperature for the Gibbs energy curve with q 1 ¼ q 2 ¼ 3.

Table 3
The residual parameters for the ternary and quaternary systems taken from Refs. [14,52].  using CALPHAD excess models, with an ideal mixing of the components with the FactSage [49] software. Therefore, combining the CALPHAD and UNIQUAC methods opens up possibilities for better equilibrium calculations and the exchange of ideas and experiences about models and methods between the CALPHAD and UNIQUAC communities. With the OpenCalphad software, this combination is now possible.
At present, only calculations for LLE are shown in this paper. The reason is that, in OpenCalphad, it is necessary to define a reference state for each element and provide a Gibbs energy function for each element or molecule relative to this reference state in each phase, gas, liquid and solid. It is possible to introduce either fugacities or non-ideal gas models, but, for calculations with several phases, these data must be introduced in a consistent way.

Results
The equations for the Gibbs energy calculation that include the contributions from UNIQUAC were translated into a Fortran script and added into the OpenCalphad software. This code is publicly accessible (http://www.opencalphad.com). This section is used to confirm the successful implementation of the UNIQUAC model in this software. First, the implementation of the UNIQUAC model in the OpenCalphad software is confirmed for artificial systems. Second, it is used to calculate thermodynamic properties and phase diagrams for binary, ternary, and quaternary systems. In the end, its ability to assess interaction parameters of a ternary system is tested.
All the systems in this section were chosen from the open literature.

Initial tests of the implementation
In order to verify the implementation of the UNIQUAC model in the OpenCalphad software, two steps are performed in this work. First, the implementation of the combinatorial excess Gibbs energy in the UNIQUAC model is tested by comparison of calculated results of thermodynamic properties obtained by Fortran code implemented in OpenCalphad, and that written independently based on the UNIQUAC model. In these calculations, structure parameters, qðAÞ ¼ 1:4, rðAÞ ¼ 0:92, qðBÞ ¼ 4:93, and rðBÞ ¼ 5:84, are used. The results are shown in Fig. 1. From this figure, we conclude that OpenCalphad calculates the combinatorial excess Gibbs energy according to the UNIQUAC model correctly.
Second, the residual part of the UNIQUAC model was implemented in OpenCalphad along with the combinatorial part. In order to verify the complete implementation of the UNIQUAC model, Fig. 4 of Abrams and Prausnitz [14] is reproduced here as it presents three typical situations in binary liquid-liquid systems: soluble (q 1 ¼ q 2 ¼ 2), critical state between soluble and insoluble (q 1 ¼ q 2 ¼ 2:5), and partially soluble with a miscibility gap (q 1 ¼ q 2 ¼ 3). The Gibbs energy curves as function of composition were calculated for three different sets of values of q using r 1 ¼ r 2 ¼ 3:3 and t 12 ¼ t 21 ¼ expð À 180=TÞ at 400 K. The results are shown in Fig. 2(a) and confirm the successful implementation. The miscibility gap for the system with q 1 ¼ q 2 ¼ 3 is shown in Fig. 2(b).

Calculation of binary systems
The implementation of the UNIQUAC model in the OpenCalphad software is tested for its capability to calculate the miscibility gap in the liquid phase using the system acetonitrile (A) -n-heptane (B). The calculations are based on the binary interaction parameters in Table 3 and compared to the experimental data of Palmer et al. [50]. The results are shown in Fig. 3. The miscibility gap is evident from the Gibbs energy curve, see Fig. 3(a) and the fact that the chemical potentials have a minimum and a maximum, see Fig. 3(b). The activity coefficients and the miscibility gap as a function of temperature have also been calculated, see Fig. 3(c) and (d). Again, good agreement is obtained between the experimental and calculated data, illustrating the ability of the implementation of the UNIQUAC model in OpenCalphad to calculate phase equilibria in binary liquid-liquid systems.

Calculation of ternary systems
In industrial processes, two types of ternary systems with partial miscibility are of interest. If there is one partially miscible binary, the ternary system is defined as type I; if there are two partially miscible binaries, the ternary system is defined as type II [53]. To illustrate these types of partial miscibility, we selected two ternary systems: acetonitrile-benzen-n-heptane (type I) and n-hexaneaniline-methylcyclopentane (type II). Their isothermal phase diagrams are calculated by the implemented UNIQUAC model in the OpenCalphad software and compared with experimental data [50,51], see Fig. 4. The interaction parameters for the calculations are taken from Tables 2 and 3.
Generally, the miscibility gap of a binary or multicomponent liquid system changes with temperature. Based on the parameters   Tables 2 and 3 Solid lines: calculated in this work; dash lines: calculated in literature [54].  Table 4 Binary interaction parameters for the system acetonitrile (A) -benzene (B) -nheptane (C), assessed by Ref. [52] and obtained in this work. in Tables 2 and 3, phase diagrams for the ternary system acetonitrile(A)-n-heptane(B)-benzen(C) at five different temperatures are calculated with the UNIQUAC model implemented in the OpenCalphad software. Fig. 5 shows that the miscibility gap closes at higher temperature. These results are in agreement with the reference results obtained by Y.C. Kim et al. [54].

Calculation of quaternary system
In order to confirm the possibility to calculate phase equilibria in multicomponent systems, an isothermal isopleth of the quaternary system 2,2,4-trimethylpentane(A) -furfural(B) -cyclohexane(C)benzene(D) is plotted in Fig. 6 at benzene mole fractions of 0.01, 0.10, 0.20 and 0.25. At low benzene mole fractions, there is a miscibility gap across the system from the 2,2,4-trimethylpentanefurfural binary to the furfural -2,2,4-trimethylpentane binary. At higher benzene contents, the miscibility gap closes from the 2,2,4trimethylpentane-furfural side. Note that there are no tie-lines in

Assessment of a ternary system
The previous subsections illustrate the successful implementation of the UNIQUAC model in OpenCalphad. The motivation of this work is to assess the interaction parameters and to calculate thermodynamic properties and phase equilibria in chemical engineering applications. In this subsection, we explore the possibility to assess the binary interaction parameters of a ternary system using the UNIQUAC model implemented in OpenCalphad, to predict the isothermal phase equilibria of the system, and then to compare the results with literature data. For this purpose, we use the ternary system acetonitrile (A) -benzene (B) -n-heptane (C). Table 4 compares three sets of parameters. The first set is taken from literature [52]. The other two sets are the parameter estimations of this work. The difference between these two sets of Table 5 Normalized Sum of Squared Errors (NSSE) of three sets of parameters listed in Table 4 Table 4; dash lines: calculated with parameters assessed by literature [52], see parameters is the different weight assigned to the experimental enthalpy data in the binary system, benzene (B) -n-heptane (C), because of the inconsistency between excess enthalpy and activity coefficient data (vide infra). The Normalized Sum of Squared Errors (NSSE) was used as a measure of the goodness of the assessment for the three sets of parameters. The parameter set obtained with assessment 2 resulted in the lowest NSSE, see Table 5. However, it is also important to consider the quality of the prediction of the ternary isothermal section, which will be discussed below.

Binary systems
The interaction parameters of acetonitrile (A) -benzene (B) are assessed based on experimental data [50,55e61]. The error bars of all the experimental data are estimated based on the error estimation made by D.A. Palmer et al. [50]. Figs. 7 and 8 show the comparisons between experimental and computational data at  Table 4; dash lines: calculated with parameters assessed by literature [52], see Table 4; symbols: experimental data observed in literature [50,58e61]. Fig. 9. Activities coefficients and excess enthalpy for the binary system acetonitrile (A) -n-heptane (C) at 318.15 K. Symbols: experimental data [50]; lines: calculated with interaction parameters from literature [52], see Table 4. different temperatures. Two types of computed results are shown based on parameters from literature [52] and parameters from this work.
From the comparison of the two figures, we conclude that the parameters in this work can predict the activity coefficients equally well as the literature parameters. However, excess enthalpies calculated with the parameters from this work are obviously more consistent with experimental data.
The interaction parameters for binary system, acetonitrile (A)n-heptane (C), are calculated from experimental data [50]. Literature values of the parameters [52] are used as a starting point. When used with the UNIQUAC model implemented in the Open-Calphad software, the parameters from the literature could not be improved upon. Therefore, they are accepted in this work. The comparison with experimental data is shown in Fig. 9.
In contrast to the previous binaries, no set of parameters could be found that adequately represented both activity coefficients and experimental excess enthalpy data [50,57,62] of the binary system benzene (B) -n-heptane (C). Therefore, two different weights, 0 and 0.5, are assigned to the experimental enthalpy data. The comparison between experimental and computational data for activity coefficients at 318 K and excess enthalpy at different temperatures is shown in Figs. 10 and 11, respectively.

Ternary system
The prediction of ternary phase equilibria is a good way to assess the capability of binary interaction parameter estimations. Therefore, isothermal phase diagrams for the ternary system, acetonitrile (A) -benzene (B) -n-heptane (C), are calculated with the binary interaction parameter sets obtained in this work, see Fig. 12. Moreover, an isothermal phase diagram for this system calculated with the parameter set reported in the literature [52] is shown in Fig. 4(a). When comparing the three figures, Fig. 12(a) shows the Fig. 10. Activity coefficients for the binary system benzene (B) -n-heptane (C). Symbols: experimental data [57]; dash lines: calculated with interaction parameters from literature [52], see Table 4; dash lines with dots and solid lines: calculated with parameters from this work, see Table 4. Fig. 11. Excess enthalpy for the binary system benzene (B) -n-heptane (C) at different temperatures. Symbols: experimental data [50,62]; solid lines: calculated by OpenCalphad software with parameters assessed in this work, see Table 4; dash lines: calculated with parameters assessed by literature [52], see Table 4.
best agreement between the calculated and experimental tie-line data. This implies that the interaction parameter set, assessment 1, is the best of the three parameter sets in Table 4. Calculated results based on the binary interaction sets in Table 4 and experimental LLE for the ternary system at 318 K are shown in Table 6.
It should be noted that the phase diagram from the literature, Fig. 4(a), was constructed using information from ternary tie-line data. In this work, on the other hand, the LLE ternary phase diagram is calculated using binary interaction parameters obtained from binary experimental data only.

Conclusion
The advantage of implementing the integral Gibbs energy expression and calculating chemical potentials and activity coefficients using a Gibbs energy minimizer like OpenCalphad is that thermodynamic consistency is guaranteed. For example, constraints such as Eq. (8) Table 4. Experimental data in both figures are from the work of Palmer et al. [50]. Table 6 Experimental and calculated LLE mole fraction for the ternary system acetonitrile (A) -benzene (B) -n-heptane (C) at 318 K. Interaction parameters from terms that have been added. Inside the framework of OpenCalphad, it would also be possible to add excess terms; see, for example, Eq.
(3). This work makes the UNIQUAC model available in OpenCalphad, which simplifies the use of the UNIQUAC model for multicomponent phase diagram calculations. It facilitates the use of thermodynamic properties like excess enthalpy and activity coefficients to assess the interaction parameters. Future applications include the possibility to combine calculations using the UNIQUAC model for the liquid phase together with the standard CALPHAD models for solid phases. In addition, this work could encourage applications of the UNIQUAC model to polymer systems by implementation of modified UNIQUAC models in the OpenCalphad software.