# How to apply sublattice model for a ternary system?

I already asked a broader version of this question in the physics-SE, but after I got no answer, combined with further understanding of the problem, it was suggested to me to try asking here.

I'm calculating the phase diagram of a ternary system (Ag-Al-Cu, to be specific), which has some phases modelled as sublattices (using this publication for the description). I have managed to model all phases for the liquidus surface, except for the epsilon phase, (Ag,Al,Cu)(Ag,Cu).

I've used the CALPHAD book of Saunders (chapter 5.4, and specifically section 5.4.2 and its equation (5.37)) to try to understand the use of the sublattice model, but I still can't model the previously mentioned epsilon phase.

From the mentioned publication, I have the following parameters for the phase: $$\epsilon\text{-phase: } (\mathrm{Ag,Al,Cu})(\mathrm{Ag,Cu})\\ {}^0 G^\epsilon_\mathrm{Al:Ag} - H^\mathrm{SER}_\mathrm{Ag} - H^\mathrm{SER}_\mathrm{Al} = 160000 + GH^\mathrm{SER}_\mathrm{Ag} + GH^\mathrm{SER}_\mathrm{Al}\\ {}^0 G^\epsilon_\mathrm{Ag:Cu} - H^\mathrm{SER}_\mathrm{Ag} - H^\mathrm{SER}_\mathrm{Cu} = 85000 + GH^\mathrm{SER}_\mathrm{Ag} + GH^\mathrm{SER}_\mathrm{Cu}\\ {}^0 G^\epsilon_\mathrm{Cu:Ag} - H^\mathrm{SER}_\mathrm{Ag} - H^\mathrm{SER}_\mathrm{Cu} = 85000 + GH^\mathrm{SER}_\mathrm{Ag} + GH^\mathrm{SER}_\mathrm{Cu}$$

together with the following from COST507, which it is based on [from its "System Al-Cu" chapter]: $$\text{Phase AlCu-}\epsilon \quad [(\mathrm{Al,Cu})(\mathrm{Cu})]\\ G^\circ(T) - H^{\circ,\mathrm{fcc-}A1}_\mathrm{Al} (298.15\mathrm{ K}) - H^{\circ,\mathrm{fcc-}A1}_\mathrm{Cu} (298.15\mathrm{ K}) = \mathrm{G(Al:Cu)} = -36976 + 1.2 T + \mathrm{GHSER}_\mathrm{Al} + \mathrm{GHSER}_\mathrm{Cu}\\ \\ G^\circ(T) - 2.0 H^{\circ,\mathrm{fcc-}A1}_\mathrm{Cu} (298.15\mathrm{ K}) = \mathrm{G(Cu:Cu)} = 2.0 \mathrm{GBCC}_\mathrm{Cu}\\ L^{0,\mathrm{AlCu-}\epsilon}_\mathrm{Al,Cu:Cu} = 7600 - 24 T\\ L^{1,\mathrm{AlCu-}\epsilon}_\mathrm{Al,Cu:Cu} = -72000$$

Based on this, I calculate the Gibbs energy of the epsilon phase with the following python code:

# Phase epsilon, (Ag,Al,Cu):(Ag,Cu)
def AgAlCuGibbsGammaH (xAg, xAl, T):
# Thermodynamic parameters from litterature:
# Ternary (Witusiewicz2005, https://doi.org/10.1016/j.jallcom.2004.06.078 )
G0AlAg = 160000 + AgG0Fcc(T) + AlG0Fcc(T) # pure elements (AlG0Fcc, AgG0Fcc, etc) from Dinsdale1991, https://doi.org/10.1016/0364-5916(91)90030-N
G0AgCu = 85000 + AgG0Fcc(T) + CuG0Fcc(T)
G0CuAg = 85000 + AgG0Fcc(T) + CuG0Fcc(T)
# Binary (COST507, https://op.europa.eu/en/publication-detail/-/publication/35a95440-10ac-45f3-bd90-0fe0a0006aaf )
G0AlCu = -36976 + 1.2 * T + AlG0Fcc(T) + CuG0Fcc(T)
G0CuCu = 2 * CuG0Bcc(T)
L0AlCuCu = 7600 - 24 * T
L1AlCuCu = -72000

# Number on the sites
N1 = 1
N2 = 1

# Site fractions
yAg1 = 1 * xAg
yAl1 = 2 * xAl
yCu1 = 1 - yAg1 - yAl1

yAg2 = 1 * xAg
yCu2 = 1 - yAg2

# Pure, ideal and excess contributions to the Gibbs energy
pure = (
yAg1 * yCu2 * G0AgCu
+ yAl1 * yAg2 * G0AlAg
+ yAl1 * yCu2 * G0AlCu
+ yCu1 * yAg2 * G0CuAg
+ yCu1 * yCu2 * G0CuCu
)

ideal = R * T * (N1 * (yAg1 * np.log(yAg1) + yAl1 * np.log(yAl1) + yCu1 * np.log(yCu1)) +
N2 * (yAg2 * np.log(yAg2) + yCu2 * np.log(yCu2)))

excess = yAl1 * yCu1 * yCu2 * (L0AlCuCu + L1AlCuCu * (yAl1 - yCu1))

# Divide by total number of sites
dividend = N1 + N2
pure /= dividend
ideal /= dividend
excess /= dividend

# Gibbs energy is equal to the sum of the different parts
return pure + ideal + excess


But it does not give me the correct common tangent plane with the liquid phase, based on the peritectic points listed in the mentioned publication.

Am I calculating the site fractions incorrectly? Or am I doing a different mistake? Any help is greatly appreciated.

Calculating the Gibbs energy for the epsilon phase and the tangent plane to the liquid phase at this point, this close to $$x_\mathrm{Ag}=0$$ leads to the ideal contribution $$x_\mathrm{Ag} \log x_\mathrm{Ag}$$ getting unstable.
In the figure, the blue surface is the Gibbs energy of the liquid phase, orange is the epsilon three-component phase Gibbs energy, dashed black line the Al-Cu binary epsilon phase Gibbs energy, and solid black line with x's at the endpoints is the tangent. $$x$$ is silver mole fraction and $$y$$ is aluminium mole fraction.