ChemistryLab.jl

Appearance

Solid solutions ​

Before this page

Thermochemistry and Standard states §2, where the standard state of an end-member is set.

A pure solid has  : it is its own standard state, its chemical potential does not depend on how much of it there is, and it either precipitates or it does not. A solid solution is a single crystalline phase whose composition varies continuously between end-members, and that changes the thermodynamics qualitatively: the activity of an end-member depends on the composition of the phase, so the phase can absorb a component gradually instead of appearing all at once.

This matters in cement more than anywhere else. C-S-H is not a compound with a formula but a phase of variable Ca/Si; the AFm phases exchange sulfate, carbonate and hydroxide on the same interlayer site; alkalis partition into hydrates rather than staying in solution. A model that only knows pure phases predicts sharp appearances and disappearances that are not observed.

1. Ideal mixing: the activity is the mole fraction ​

Put moles of each end-member into one phase, with mole fractions  . Even with no interaction at all, mixing lowers the Gibbs energy, because there are more ways to arrange a mixture than a pure solid. Per mole of phase,

which is entropy alone — the of mixing. Differentiating for one end-member gives its chemical potential and hence its activity:

That is IdealSolidSolutionModel, and it is the default. Note what it already buys: an end-member at   has   , so its saturation index is shifted by decades relative to the pure phase. A trace component is stabilized by being diluted in a host, which is why a solid solution can take up an ion that would never precipitate as its own phase.

The curve, its tangent, and the plane of the rest of the system ​

Per mole of phase, the Gibbs energy of a binary of end-members A and B is the straight line of the two end-members unmixed,     , plus , negative between the two: the curve lies below the line, and the same matter costs less mixed than apart (left panel). The tangent to the curve at a composition meets the edges   and   at the chemical potentials of the two end-members in the phase, and it is through these potentials that the phase exchanges matter with the rest of the system.

At equilibrium each element has one potential, shared by every phase, and those potentials give each end-member a potential in units of , which defines a plane over the compositions of the phase (right panel). A curve lying above that plane everywhere is the Gibbs energy of an absent phase, none of its compositions being as cheap as what the rest of the system already offers; a curve touching it is that of a present phase, at the composition of the contact; a curve dipping below it is that of a phase which would lower the Gibbs energy by forming, so that the state is not yet an equilibrium. With  , the phase is present exactly when    , which the certifying solver uses as an equation. For ideal mixing  , and   is then the saturation ratio end-member would have as a pure phase: an ideal solid solution forms as soon as reaches one, while every is still below one and no end-member would precipitate on its own. Solid solutions in a calculation reads these ratios on a C-S-H gel.

2. Non-ideal mixing: the excess Gibbs energy ​

Real end-members interact. Everything beyond the ideal term is collected into the excess Gibbs energy , and the activity coefficient is its partial molar derivative:

and then   . This is exactly what the code does: for these models _ss_log_activities! writes log(x[k] + ϵ) + _excess_ln_gamma(model, k, x, T) into the log-activity vector, so the two terms are visibly separate and a model supplies only the second. The model of section 7 is written differently, and has its own method.

3. RegularSolutionModel — one interaction energy per pair ​

The simplest non-ideal form keeps only pairwise contacts, with one energy per pair of end-members:

which is the symmetric multicomponent Margules expression, and it is _excess_ln_gamma(::RegularSolutionModel, …) line for line. For a binary it collapses to the form usually quoted:

In the picture of nearest neighbors behind the model, (J/mol) is the energy of an – contact less the mean of the – and – contacts, counted over the contacts of a mole. Fitted to data, it often has to depend on temperature or pressure, and it then mixes enthalpy and entropy, but its sign keeps the sign of the departure from ideal mixing (Richet, 2001; Sec. 11.4b, pp. 262–263), and that sign is physical:

  •   — unlike neighbors are unfavorable.  , the phase resists mixing, and above a threshold it unmixes;

  •   — unlike neighbors are favorable.  , mixing is stabilized beyond ideal, and the phase may order.

W must be square and symmetric, which a validating inner constructor enforces; the reason it is inner rather than outer is recorded in the source, and it is a mistake worth not repeating.

The threshold is  , and it is checkable ​

For a symmetric binary, the molar Gibbs energy of mixing is       . A phase is unstable to unmixing where that function is concave, and

which at   is  : negative as soon as  . So   — about 5 kJ/mol at 25 °C — is the critical point of a symmetric regular solution: below it one homogeneous phase, above it a miscibility gap that widens as grows. The threshold is evaluated, and the sign of the second derivative tabulated on either side of it, in Solid solution models, in numbers.

This matters in practice because nothing in the activity expression detects it. Past the threshold it goes on returning values, and those values describe a single phase that should have split. The check is therefore made at construction: SolidSolutionPhase evaluates the second derivative and refuses a model that unmixes at the temperature given, naming the interval, unless the phase is declared with two instances (section 6) or the check is switched off with check_convexity = false.

4. RedlichKisterModel — a binary that need not be symmetric ​

A regular solution is symmetric by construction: swapping the end-members leaves unchanged. Real binaries often are not, so the excess is expanded in a polynomial of the composition difference. As implemented, for two end-members,

with in J/mol. Read the roles off the formulas: is the symmetric term and is exactly a regular solution's ; and are the asymmetric corrections, so setting them to zero must reproduce the regular model. That is an identity between two independently written methods, and it is checked rather than trusted in Solid solution models, in numbers, where the two agree to the last digit.

RedlichKisterModel requires exactly two end-members, and SolidSolutionPhase refuses the combination at construction rather than letting a ternary reach an expression written for a binary. For three or more end-members the choices are the ideal model, RegularSolutionModel with a full matrix, or the asymmetric models of section 9: the subregular model, whose binaries are this model to the first order, the series of every order, and the van Laar model.

Published dimensionless parameters, and the temperature ​

Cemdata18 gives its non-ideal AFt and AFm binaries dimensionless Guggenheim parameters and (Lothenbach et al., 2019) (its Table 1, notes a, b, g and h, in data/literature/Lothenbach2019.json), and no temperature dependence for them. The package turns them into   with   K and keeps in J/mol at every temperature: the excess Gibbs energy is then all enthalpy, and the interaction it puts into at another temperature is , 7 % larger at 5 °C and 16 % smaller at 80 °C. Holding itself constant would make the excess all entropy and leave the same at every temperature. The two limits have classical names, a regular solution for the first and an athermal one for the second (Richet, 2001; Secs. 11.4b and 11.5e, pp. 262 and 270). Neither choice is in the source, and no data in the package decides between them; the edges of a gap measured at two temperatures would. What the choice does to the gap of the AFm SO₄/OH binary, monosulfate and C₄AH₁₃, whose printed edges are the 25 °C ones:

julia
using ChemistryLab, DynamicQuantities, Printf
p = literature_row("Lothenbach2019", "guggenheim_parameters", "AFm SO4/OH")
α0, α1 = ustrip(p.alpha0), ustrip(p.alpha1)
rk(T) = RedlichKisterModel(a0 = α0 * R_GAS * T, a1 = α1 * R_GAS * T)
@printf("printed gap: %.2f to %.2f\n", ustrip(p.gap_from), ustrip(p.gap_to))
println(" °C   a in J/mol held (the package)   α held")
for T_C in (5, 25, 50, 80)
    T = T_C + 273.15
    held_a = common_tangent(rk(298.15); T)
    held_α = common_tangent(rk(T); T)
    @printf("%3d       %.3f to %.3f            %.3f to %.3f\n", T_C, held_a..., held_α...)
end
printed gap: 0.50 to 0.97
 °C   a in J/mol held (the package)   α held
  5       0.478 to 0.977            0.500 to 0.970
 25       0.500 to 0.970            0.500 to 0.970
 50       0.528 to 0.959            0.500 to 0.970
 80       0.565 to 0.943            0.500 to 0.970

Both give the printed gap at 25 °C. The package's choice narrows it as the temperature rises, from 0.48–0.98 at 5 °C to 0.57–0.94 at 80 °C (the mole fraction of C₄AH₁₃, the OH end, as the edges are printed); the other keeps it where Cemdata18 prints it.

5. What the three models do to an activity ​

Evaluated over the whole composition range in Solid solution models, in numbers, with    kJ/mol — inside the stable range of §3, so that the numbers describe a phase that stays homogeneous. Two things to read off it.

A positive pushes the activity of a dilute end-member above its mole fraction: the host is rejecting it. A negative pulls it below: the host is stabilizing it. And at   all three models converge, because   as the phase becomes pure — the standard state of an end-member is the pure end-member, so they agree there by construction.

The ideal case is the one to keep in mind for cement. An end-member at   has  , so its saturation index sits three decades below the pure phase: a trace component is stabilized simply by being diluted in a host. That is how a solid solution takes up an ion which would never precipitate as a phase of its own, and it is the whole reason C-S-H and the AFm phases can absorb alkalis, sulfate and carbonate continuously instead of in jumps.

6. The spinodal, the common tangent, and why one amount is not enough ​

Sections 3 and 4 gave models whose excess term can be strong enough to make the molar Gibbs energy of mixing concave over an interval. This section is about what that means for the equilibrium, because it is not a detail of the model: it changes what the answer is.

Convexity is the hypothesis under everything else ​

For a binary write the molar Gibbs energy of mixing, in units of , as

The ideal part is convex everywhere — its second derivative is     — so unmixing is always the excess term's doing. Where

the phase is inside its spinodal, and there the straight line joining two compositions lies below the curve between them. A system at an intermediate overall composition therefore lowers its energy by separating into those two, and the minimum of is not a point but a pair.

For the symmetric regular solution of section 3 the criterion is exactly    at  , so unmixing begins at   — the threshold section 3 already quoted, here derived from the same inequality. spinodal_interval evaluates the second derivative on a grid rather than in closed form, which is why the same routine covers the asymmetric Redlich-Kister case without a separate derivation.

The answer inside a gap is the common tangent ​

Which pair? Write now for the molar Gibbs energy of the phase with its standard terms included,

the mole fraction of the end-member B and the energy of mixing above; the linear terms change neither the second derivative nor, as shown below, the pair. Let moles of the phase, of overall composition , be shared between two parts and holding and moles at compositions and . At fixed and the equilibrium minimizes

the conservation of the two end-members. With a multiplier for the first constraint and for the second, the Lagrangian is          , and its stationarity with respect to the compositions and to the amounts, both parts being present ( ), reads

The four equations say one thing: the straight line    passes through the curve at and at and has the slope of the curve at both, so that it is tangent to at the two compositions, the common tangent. Eliminating gives the two equations in two unknowns that common_tangent solves:

The multipliers are chemical potentials. The tangent to at a composition meets the vertical   at   , the chemical potential of A at that composition, and the vertical   at    , that of B, partial molar quantities being read as the intercepts of the tangent (Richet, 2001; Sec. 3.2c, Fig. 3.2, pp. 51–52). The stationarity therefore states that   and    in both parts: each end-member has the same chemical potential in the two compositions, which is the condition of equilibrium between any two phases (Richet, 2001; Sec. 7.3, Fig. 7.5, pp. 167–169).

Nothing requires the two points to be at the same height:    , which vanishes only for a horizontal tangent,  , as for a symmetric mixing energy between end-members of equal standard potentials. Nor are they the minima of , where   while here  . A term linear in , as the standard potentials add, tilts the curve and its tangent alike and moves neither point: the pair depends on the mixing model only.

What the split gains follows from the lever rule,     : the Gibbs energy of the two parts is   , the height of the tangent at , below the   of a single composition wherever the curve rises above its tangent. Between and the equilibrium is this mixture, and the energy follows the tangent rather than the curve. Stationarity alone does not select it, however: a single composition,   , satisfies the same equations with the tangent drawn at . What separates the two is a condition on the whole curve, that it lie nowhere below the line, which is the tangent-plane test of Michelsen (1982), the subject of a subsection below.

The common tangent construction is wider than the spinodal: the binodal contains the spinodal, and between the two the phase is metastable rather than unstable. Inside the spinodal any small fluctuation of composition lowers and grows by itself; between the spinodal and the binodal every small fluctuation raises it, and the phase splits only through a nucleus large enough to pay for its interface (Richet, 2001; Sec. 7.3c, p. 170). A minimization sees only the tangent.

The right panel states the difficulty the rest of this section resolves: inside the gap the equilibrium holds the phase at two compositions at once, in the proportions of the lever rule, which a formulation carrying one amount per end-member has no means to express.

Computing the pair from the model alone ​

The pair is not found by minimizing. It is computed from the model alone, by common_tangent: two equations in two unknowns, solved by Newton with ForwardDiff supplying the derivatives. Nothing about the rest of the chemical system enters — not the aqueous solution, not the other phases, not the element budget — which is why it costs microseconds and can be used as the starting point of a full equilibrium rather than as its result.

The construction is that of Glynn and Reardon (1990) for binary solid solutions, which PHREEQC (Parkhurst and Appelo, 2013) uses as well. It was validated here against a case with a closed form: for a symmetric model the pair must be symmetric about  , and it comes out at with a residual of   and the symmetry exact to the last bit. For a symmetric regular model the pair has a closed form,     (Richet, 2001; Sec. 10.4b, Eq. (10.29), p. 230), which these two compositions satisfy with   to the six figures given.

miscibility_split then applies the lever rule. Inside the gap the two compositions are fixed and only their proportions move with the overall composition — which is what makes a miscibility gap flat in a phase diagram, and what a single-composition answer cannot reproduce at any resolution.

Why a formulation with one amount per species cannot hold it ​

The composition vector of a ChemicalSystem has one entry per species, so a species belongs to exactly one phase and a phase has exactly one composition. A common-tangent pair is two compositions of one substance, and there is no way to write that down with one amount each.

So the formulation has to be given the substance twice. That is what SolidSolutionPhase(...; instances = 2) does: ChemicalSystem builds a second copy of each end-member under a derived symbol (monosulphate12#2) sharing the same thermodynamic record, and each copy belongs to its own phase. The groups stay disjoint, the duplicated columns of the conservation matrix are copies of columns already there — so the row rank, and with it conservation, is untouched — and the minimization is free to populate either lobe or both.

It is refused for a convex model, and the reason is worth stating because it is not conservatism. Two instances of a convex phase are degenerate: the energy is the same however the amount is split between them, so the minimum becomes a flat manifold and the optimizer is asked to pick a point on it arbitrarily. Inside a spinodal the common-tangent pair is unique and no such direction exists.

Michelsen's tangent-plane distance, which is what detects it ​

Representation is one problem; knowing that it is needed is another. A phase sitting inside its own spinodal satisfies every first-order condition — its end-members are stationary, the element balance closes, nothing absent is supersaturated — so stationarity alone certifies a non-minimum.

The test that separates them is Michelsen's (Michelsen, 1982): for a trial composition of the phase, the tangent-plane distance

measures how far the Gibbs surface at lies above the tangent plane drawn at the reported composition . If is negative anywhere, some other composition lies below that plane and the reported state is not a minimum. At   the distance is exactly zero and is a stationary point of the search, which is why the measure has to be started from the end-member corners — the other lobe is where the negative value lives.

OptimaSolver's phase_split_measure computes this for every present mole-fraction phase and folds the result into the certificate. Concretely: check_convexity = false stops being a silent loss of the proof — a concave declaration that does unmix now fails to certify, with the phase named.

A second composition, declared or added when needed ​

The certificate's measure of a phase is the sum over its members of their activities from the dual solution divided by their activity coefficients from the primal one, an estimate of the mole fraction the phase would take:

It is the phase stability index   of Kulik et al. (2013), who derive it from the same KKT conditions and call "a generalization of the saturation index". The package evaluates it in two places: as   (the ideal case,  ) when it repairs a start, and through the log-sum-exp of     in phase_split_measure.

Holding a gap takes the phase twice. Cemdata18 ships the AFm and AFt binaries under two names each, so that the database itself declares each binary twice; instances = 2 declares the second composition with a keyword. instances = :auto asks for less: the phase is solved with one composition, and a second instance is added only when the certificate finds it wanting to split. On a binary whose overall composition the element budget holds inside the gap, calcite and magnesite in equal amounts, that finds the common-tangent pair:

julia
using ChemistryLab, DynamicQuantities, Printf
sp = Dict(symbol(s) => s for s in build_species(datapath("slop98-inorganic-thermofun.json"); verbose = false))
gap = RedlichKisterModel(a0 = 14_000.0)
carbonate_ss = SolidSolutionPhase("carbonate", [sp["Cal"], sp["Mgs"]]; model = gap, instances = :auto)
cs = ChemicalSystem([sp[s] for s in split("H2O@ H+ OH- CO2@ HCO3- CO3-2 Ca+2 Mg+2 Cal Mgs")],
                    ["H2O@", "H+", "Ca+2", "Mg+2", "CO3-2", "Zz"]; solid_solutions = [carbonate_ss])
st = ChemicalState(cs)
set_quantity!(st, "H2O@", 1.0u"kg"); set_quantity!(st, "Cal", 0.025u"mol"); set_quantity!(st, "Mgs", 0.025u"mol")
eq, cert = equilibrate_certified(st; b = budget(st))
n = ustrip.(us"mol", eq.n)
@printf("certified %s, instances added to %s\n", cert.optimal, cert.instances_added)
for (g, ph) in zip(eq.system.ss_groups, eq.system.solid_solutions)
    @printf("  %-12s %.5f mol, x(calcite) = %.4f\n", name(ph), sum(n[g]), n[g[1]] / sum(n[g]))
end
@printf("common tangent of the model: %.4f and %.4f\n", common_tangent(gap)...)
certified true, instances added to ["carbonate"]
  carbonate    0.02495 mol, x(calcite) = 0.9963
  carbonate#2  0.02486 mol, x(calcite) = 0.0037
common tangent of the model: 0.0037 and 0.9963

It does not make a nonconvex problem easy: the certificate names a phase to split only once the one-composition solve has converged.

The executed counterpart of this section is the miscibility-gap page, which runs one cement three ways and reports the certificate each time.

7. Mixing on several sites: SublatticeModel ​

Sections 1 to 4 mix end-members. In a C-S-H the entities that actually mix are smaller: a few structural positions of the silicate chain and of the interlayer, each held by one species or another. An end-member is then one particular filling of those positions, and mixing is random on each position separately. That is the model of Kulik (2011) for CSH3T and of Myers et al. (2014) for the CNASH gel, and it is what SublatticeModel implements.

The model. A formula unit has sites , site counted times, and end-member puts the species on site . The fraction of site held by species is the total mole fraction of the end-members that put it there,    , and the configurational Gibbs energy is the ideal one of each site:

The implication is  , exact because every end-member fills every site once. One site counted once is section 1 again. Kulik (2011) and Myers et al. (2014) write the same model as a fictive activity coefficient  , which is what excess_ln_gamma_expression returns for this model.

CSH3T. Two bridging-tetrahedral sites hold either Si or Ca. TobH puts Si on both, T2C puts Ca on both, and the ordered T5C puts Ca on the first and Si on the second (Kulik, 2011; Eq. 19). The Cemdata18 records are half of its formula units, so on them each site counts one half:

julia
using ChemistryLab
subs = Dict(symbol(s) => s for s in build_species(datapath("cemdata18-thermofun.json"); verbose = false))
csh3t = sublattice_model("Kulik2011:csh3t", [subs[n] for n in ("CSH3T-TobH", "CSH3T-T5C", "CSH3T-T2C")])
csh3t.multiplicity
2-element Vector{Float64}:
 0.5
 0.5

sublattice_model finds that factor itself, by comparing each record's formula with the one Kulik (2011) prints, and refuses a record that is not such a multiple. At   the first site holds Si only through TobH, the second through TobH and T5C:

julia
site_fractions(csh3t, [0.2, 0.5, 0.3])
2-element Vector{Vector{Float64}}:
 [0.2, 0.8]
 [0.7, 0.3]

and the fictive activity coefficient of T5C is Eq. (20) of Kulik (2011), halved:

The block below builds that expression from the model and subtracts the formula above; the difference simplifies to zero:

julia
using Symbolics
x = [Symbolics.variable(:x, i) for i in 1:3]   # x₁ = TobH, x₂ = T5C, x₃ = T2C
λ = excess_ln_gamma_expression(csh3t, 2, 3)
simplify(expand(λ - (0.5log(x[1] + x[2]) + 0.5log(x[3] + x[2]) - log(x[2]))))

Two properties the solver relies on.

  • Convexity. is a sum of convex functions of the site fractions, which are linear in , so it is convex, and strictly so when the occupancy matrix (one row per site and species, one column per end-member) has full column rank, stored as rank. The activities are the gradient of that energy, so a certificate keeps its :global_minimum scope.

  • The dilute limit. As   only the site fractions of the species alone carries vanish with it:  , with the multiplicity of those sites, stored as exponents. T5C owns no species of its own,  : its activity stays finite as it disappears, and it can be absent from a CSH3T that is present.

julia
csh3t.exponents
3-element Vector{Float64}:
 0.5
 0.0
 0.5

What the solver does with it. The composition of a mixing phase is found by successive substitution, and on a sublattice model that iteration diverges: its rate is set by the sum of the multiplicities, nine for the CNASH gel. Such a phase is therefore inverted by Newton's method, and the members with   are declared to the certificate as able to leave the phase, which then tests them by the inequality a pure phase obeys (OptimaSolver's SolutionPhase, keywords newton and bounded_members). A system with a sublattice phase is first solved with that phase mixing its end-members ideally, and the answer is the start of the sublattice solve.

data/solid_solutions.toml ships CNASH_ss as the model of Myers et al. (2014) and CSH3T as Cemdata18 ships it, ideal between end-members; sublattice_model("Kulik2011:csh3t", members) gives the site form.

8. Sites with the energies of the compounds: CompoundEnergyModel ​

Section 7 mixes ideally on the sites and keeps the end-members' own standard energies, . That is exact only when the end-members are independent in the site fractions. The CASH+ model of C-S-H (Kulik et al., 2022) has two mixing sites: the bridging tetrahedron holds a silicate, a vacancy or a calcium (S, v, C), and the interlayer a vacancy or a calcium (v, C). Its six end-members are the six ways of filling them, TSvh to TCCh, but three site fractions fix the composition, and the six members are not independent. The reciprocal reaction

leaves every site fraction unchanged, yet its standard Gibbs energy is not zero. The compound energy formalism, in which Kulik et al. (2022) write the model, gives such a reaction its energy. CompoundEnergyModel implements it.

The Gibbs energy. The end-members must be every compound of the sites, each once. Compound puts species on site (the of section 7). Per formula unit,

The first term, the reference surface, equals at compound and is linear in the fractions of each site, so a reaction between compounds that changes no site fraction keeps its energy. The last term is a regular interaction between two species of one site, symmetric, with the parameters of Kulik et al. (2022). depends on the site fractions alone.

The activities. With  , the chemical potential of member is  . For a function of the site fractions,    , so

On the reference surface    for every site, because is linear in each site's fractions. With mixing sites, and writing   ,

The middle term vanishes at every compound. It also does not change when each is shifted by the energies of its elements, since those are a sum over the sites of the species' shares, which the reference surface interpolates exactly. So the convention of the database does not matter, but the energies of the members do: the solver hands them to the model at the temperature of the solve.

Three consequences.

  • The amounts are not unique. Every split of the same site fractions between the members has the same and the same element content: in CASH+, six amounts for three independent site fractions and one total. Left to choose, the solver returns one of them, which one depending on where the search started. It is therefore given one split, the product of the site fractions,  , by adding to the energy it minimizes  , with

    is the divergence of from the product of its own site fractions: it is never negative, and it vanishes only at that product, where its gradient   vanishes too. The minimum over the splits of one set of site fractions is therefore the product, with the energy and the chemical potentials of the model. The activity of a vanishing member now goes as its mole fraction, as in a mixture of end-members. Measured on twelve pastes of CASH+NK, the equilibria certify with or without , at the same cost; what adds is an answer that does not depend on the start.

  • Convexity must be decided. The reference surface is multilinear, and a reciprocal energy can make concave somewhere. In the site fractions, each site contributes a curvature of at least  , the configurational bound of section 9 plus the smallest tangent eigenvalue of its interactions , and with two sites the reference surface couples them through a constant matrix , the of the compounds projected on the two tangent spaces. So is convex when    . That proves the CASH+ core convex (34.3 against 13.2 at 25 °C), and a certificate on it keeps its :global_minimum scope as far as the gel goes; the bound does not decide CASH+NK, whose sampled Hessian shows no concave point, and its scope is then :local_minimum or :kkt_point, according to the curvature at the answer (mixing_convexity).

  • Euler's relation holds.    for any amounts with the site fractions , as it must for a Gibbs energy.

The block below checks the gradient and the last point on the core model, at a composition of the six members. The activities are compared with the gradient of   coded from the formulas above and differentiated by ForwardDiff; at the product split, and its gradient vanish:

julia
using ChemistryLab, DynamicQuantities, ForwardDiff, LinearAlgebra
subs = Dict(symbol(s) => s for s in build_species(datapath("cemdata18-cashplus.json"); verbose = false))
members = ["TSvh", "TSCh", "Tvvh", "TCvh", "TvCh", "TCCh"]
cash = compound_energy_model("Kulik2022:cashplus", [subs[m] for m in members])
T = 298.15
g = [ustrip(us"J/mol", subs[m][:ΔₐG⁰](T = T * u"K", P = 1.0e5u"Pa"; unit = true)) for m in members] ./ (ChemistryLab.R_GAS * T)
o = cash.lattice.occupancy
function nG(n; split = true)
    x = n ./ sum(n)
    y = site_fractions(cash, x)
    G = sum(g[j] * y[1][o[1, j]] * y[2][o[2, j]] for j in 1:6) + sum(v * log(v) for ys in y for v in ys)
    G += sum(W / (ChemistryLab.R_GAS * T) * y[s][i] * y[s][l] for (s, i, l, W) in cash.interactions)
    split && (G += sum(v * log(v) for v in x) - sum(v * log(v) for ys in y for v in ys))   # D(x)
    return sum(n) * G
end
n = [0.10, 0.25, 0.05, 0.20, 0.15, 0.25]
lna = ChemistryLab._ss_log_activities!(zeros(6), 1:6, n ./ sum(n), cash, T, 0.0, g)
μ = ForwardDiff.gradient(nG, n)
y = site_fractions(cash, n ./ sum(n))
xp = [y[1][o[1, j]] * y[2][o[2, j]] for j in 1:6]          # the product split
(gradient = maximum(abs, lna .+ g .- μ), euler = dot(n, μ) - nG(n),
 D_at_product = nG(xp) - nG(xp; split = false))
(gradient = 6.821210263296962e-13, euler = 0.0, D_at_product = 0.0)

All three are at the rounding of the arithmetic. The model and its data are in data/literature/Kulik2022.json, the sodium and potassium of Miron et al. (2022), Miron et al. (2022) in Miron2022a.json and Miron2022b.json, and the twelve end-members in the database cemdata18-cashplus.json. The CASH+ page computes the C-S-H in water and in alkali solutions with it.

9. More than two end-members ​

Asymmetric mixing of any number of end-members: SubregularSolutionModel ​

The regular model of section 3 is symmetric on every edge, and the Redlich–Kister model of section 4 is asymmetric but binary. Helffrich and Wood (1989) write the asymmetric Margules model for any number of end-members (their Eq. 5′):

On the edge – the bracket sums vanish and the excess is  : is at infinite dilution of in , and the converse. The ternary coefficients are a separate measurement. Helffrich and Wood (1989) show that the binaries do not determine them, even when every binary is symmetric, and that no coefficient of higher order appears.

Helffrich and Wood (1989) give the activity coefficients (their Eq. 6′) without the intermediate steps. They are short once the excess is written on the simplex, where    . The bracket of a pair is then

so each pair is a Redlich–Kister term of the first order,   and  , and two end-members are exactly the model of section 4 with those two coefficients. Write   in this form, as a function of . The excess of moles is ,  , and since   ,

The formula holds whatever is off the simplex, since reads only on it. For a pair,     has

and a triple, , has   equal to three times itself. Collecting, with every coefficient divided by ,

with   , which is _excess_ln_gamma(::SubregularSolutionModel, …) term by term. Setting   and every to zero leaves the regular model of section 3,  . The test suite checks this expression against (Helffrich and Wood, 1989; Eq. 6′), for each end-member of a quaternary with its twelve binary and four ternary coefficients, and Gibbs–Duhem and   as identities.

A series for each pair: MulticomponentRedlichKisterModel ​

Redlich and Kister (1948) write the excess of a binary as a series in the difference of the two mole fractions, which changes sign when the components are exchanged, and that of a ternary as the sum of its three binaries, each evaluated at the mole fractions of the ternary as they are, plus a term in the product of the three (their Eqs. 17 and 18):

The order of a pair matters for its odd terms, which change sign with it, and the order of a triple for and . Beyond three components Redlich and Kister (1948) add the further pairs and triples and a term in four mole fractions, which it does not write out and which is not implemented. Its coefficients are in units of  , the paper working with decimal logarithms.

The activity coefficients follow from Eq. (14) of Redlich and Kister (1948), the partial molar derivative written in mole fractions, which is the formula of the subregular model above. For a pair, with   ,   and  , the term   has

and for a triple, with      , the term   has

is the sum over every pair and triple of  , plus the derivatives with respect to of the terms that contain it, which is _excess_ln_gamma(::MulticomponentRedlichKisterModel, …). Two end-members give the model of section 4 with  . To the first order the series is the subregular model above, with    and   : Redlich and Kister (1948) call the first-order binary "the equation of Margules", and the subregular section shows why the ternary extensions coincide.

Redlich and Kister (1948) also give the ratio of two activity coefficients in a ternary (their Eq. 19). Its term in ,   , is the value on the binary: in the ternary, Eq. (14) gives   , which is that plus   . Every other term of Eq. (19) agrees with Eq. (14). The implementation follows Eq. (14), and the test suite checks it against the paper's worked ternary, heptane, methanol and toluene with the coefficients of its Eq. (21) and the association of methanol left out, through its Eqs. (22) and (23).

Unequal sizes: VanLaarModel ​

Holland and Powell (2003) make the regular model asymmetric by giving each end-member a size and weighting the mole fractions by it,  :

Only the ratios of the sizes matter. Equal sizes give back the regular model of section 3, which Holland and Powell (2003) call the symmetric formalism, and two end-members give the van Laar binary, asymmetric as soon as the sizes differ. The papers that fit and often make them depend on temperature and pressure; the model takes their values at the conditions of the calculation.

Holland and Powell (2003) write the activity coefficients (their Eq. 2) with interactions rescaled by the size of the end-member considered,   , their double sum running over  , the range the excess energy implies. A shorter route goes through the excess itself. Substituting ,

a quadratic form over a linear one. The excess of moles is then  , homogeneous of degree one in the amounts, and its derivative is

which is _excess_ln_gamma(::VanLaarModel, …). With equal sizes ,   and  , and the regular expression of section 3 comes back; for a binary it is Eqs. (4) and (5) of Holland and Powell (2003),    . The test suite checks it against those equations for the alkali feldspar of the paper, against its Eq. (13) for calcite, magnesite and dolomite, with the paper's parameters, and against .

Neither this model nor the series above has the Redlich-Kister form of section 4 with three coefficients, so their binary is scanned for a spinodal and its common tangent found from the excess itself, , in the mole fraction of the first end-member.

Convexity with more than two end-members ​

A one-dimensional scan decides a binary. With end-members the question is whether the Hessian of is positive on the tangent space of the simplex,   , at every composition, and mixing_convexity answers it by model.

For a regular model,    with  , and the Hessian is  . The first term is at least 2 on the tangent space: for a unit with  , the Cauchy–Schwarz inequality gives   , and the ℓ1 norm of such a is at least (its positive and negative parts have equal sums). Hence

an orthonormal basis of the tangent space. For two end-members this is  , the threshold of section 3; for more it is a sufficient condition. In the other direction, a pair with   is a witness: at the middle of that edge the second derivative along it is   . Between the two, the smallest eigenvalue of the projected Hessian is searched on a lattice of compositions, and a negative one is a witness too; finding none proves nothing, and the verdict is then :undecided. The subregular, the Redlich-Kister and the van Laar models are decided by that search alone beyond two end-members: their Hessian depends on the composition through the asymmetric and the ternary terms, and no bound of the kind above is written for them.

julia
using ChemistryLab, DynamicQuantities, Printf
RT = ChemistryLab.R_GAS * 298.15
W(a, b, c) = [0.0 a b; a 0.0 c; b c 0.0] .* RT
[mixing_convexity(RegularSolutionModel(W(1.9, 1.9, 1.9)), 3).verdict,
 mixing_convexity(RegularSolutionModel(W(3.0, 0.0, 0.0)), 3).verdict]
2-element Vector{Symbol}:
 :convex
 :nonconvex

A phase whose verdict is :nonconvex is refused at construction, as a concave binary is, and admitted with two instances. The certificate of an answer is not scoped :global_minimum unless every mixing phase is proved convex; it is then :local_minimum where the Hessian over the directions that conserve matter is positive definite at the answer, and :kkt_point otherwise.

A member that is a mixture of two others ​

Two different things can make one end-member the mixture of two others in composition. In an ordered member the difference of Gibbs energy is the point: the pentameric T5C of CSH3T is the average of TobH and T2C less an ordering energy (Kulik, 2011; Eq. 18), and the siliceous hydrogarnet C3AS0.41H5.18 lies between C3AH6 and C3AS0.84H4.32. When the Gibbs energy is the average as well, the phase holds one substance twice, once as a member and once as a mixture of two, and ideal mixing counts its configurations twice. SolidSolutionPhase warns when that holds within  ; the MgAl-OH-LDH ternary of Cemdata18 is the case, its M6A member the average of M4A and M8A:

julia
subs = Dict(symbol(s) => s for s in build_species(datapath("cemdata18-thermofun.json"); verbose = false))
g(s) = ustrip(us"kJ/mol", subs[s][:ΔₐG⁰](T = 298.15u"K", P = 1.0e5u"Pa"; unit = true))
@printf("MgAl-OH-LDH  M6A - (M4A + M8A)/2  = %+.4f kJ/mol\n",
        g("M6A-OH-LDH") - (g("M4A-OH-LDH") + g("M8A-OH-LDH")) / 2)
@printf("CSH3T        T5C - (TobH + T2C)/2 = %+.2f kJ/mol\n",
        g("CSH3T-T5C") - (g("CSH3T-TobH") + g("CSH3T-T2C")) / 2)
MgAl-OH-LDH  M6A - (M4A + M8A)/2  = +0.0003 kJ/mol
CSH3T        T5C - (TobH + T2C)/2 = -4.35 kJ/mol

The LDH difference, 0.3 J/mol, is about 10⁻⁴ RT, where the −4.35 kJ/mol of T5C is an ordering energy. The shipped LDH entry keeps the published model with acknowledge_degenerate = true.

10. How a solid solution is declared ​

A SolidSolutionPhase names its end-members and carries a model:

julia
phase = SolidSolutionPhase("CSHQ", [sp["CSHQ-JenD"], sp["CSHQ-JenH"],
                                    sp["CSHQ-TobD"], sp["CSHQ-TobH"]];
                           model = IdealSolidSolutionModel())

Every end-member must be AS_CRYSTAL, and the phase is checked against its model's arity and against the convexity of its mixing energy at construction; instances is how a declaration asks for more than one coexisting composition, which section 6 is about. Sets of end-members with their interaction parameters are read from a TOML through build_solid_solutions, which is also the only interaction-parameter file format in the package — data/solid_solutions.toml is the shipped example.

Two consequences to keep in mind when reading results:

  • a ChemicalSystem groups end-members into ss_groups, and every aqueous activity model must call the solid-solution branch or the end-members come back with  , i.e. treated as pure phases;

  • an end-member's saturation_indices entry is relative to its current mole fraction, since its activity is . For an end-member sitting at the solver's lower bound that is a statement about a vanishing phase, not about whether the solid solution would form.

See also ​