We develop a new optimized Metastable EXchange eQuilibrium ALgorithm (MEXQAL) which simplifies and extends the algorithm presented in Ghiorso (2013). In the original presentation of Ghiorso (2013), metastable exchange equilibrium was considered in the context of an omni-component phase (like melt or aqueous fluid), which can freely exchange atoms with all other coexisting (or precipitating) phases in the assemblage. Building upon this previous work, we seek to broaden the useful scope for this algorithm to include cases lacking an omni-component phase, where determining saturation affinities and compositions for an externally imposed chemical potential also has great general utility. Absent an omni-component phase, the stable assemblage consists of a collection of phases which together minimize the total Gibbs energy and match the total bulk composition of the system. Even in this case, a metastable exchange algorithm is highly useful, as briefly discussed in Ghiorso (2013) which introduced the concept of a hypothetical omni-component phase. Here we consider arbitrary phase assemblages which can generally lack an omni-component phase, where a metastable exchange algorithm can allow us to determine the most energetically favorable composition for each solution phase given a guess of the elemental chemical potentials of the system. With such a computational tool, we can then iteratively determine the elemental chemical potentials by minimizing the total Gibbs energy of the system.
There are numerous standard thermodynamic calculations that could make immediate use of a computationally efficient algorithm of this type. Within geological applications, calculating subsolidus phase assemblages is a typical application, used to determine in situ chemical and physical properties of a rock and extreme T/P conditions. Similarly, we can also use this algorithm to determine the equilibrium assemblage for a single solution phase that undergoes exsolution, where the bulk composition potentially lies inside of a miscibility gap. This method would be a critical part of the overall algorithm enabling the efficient and accurate calculation of tielines and 3-phase triangles for an immiscible bulk composition (as well as potential higher order coexisting assemblages). Metastable exchange equilibrium thus provides a useful reference frame for modeling equilibrium processes, even in the absence of an omni-component phase.
Metastable EXchange EQuilibrium ALgorithm [MEXQAL] overview
Our new Metastable EXchange EQuilibrium ALgorithm (MEXQAL) improves and generalizes the saturation state algorithm presented in Ghiorso (2013). It utilizes a similar iterative analytic approach, but simplifies the composition update calculation while utilizing additional derivative information. These changes improve numerical accuracy while cutting runtime by over an order of magnitude, reducing the number of Gibbs energy evaluations and remaining accurate for both major and trace components. The algorithm is outlined below.
We consider a solution phase in metastable exchange equilibrium with an externally imposed chemical potential (\(\hat{\mu}_i\)):
\[\Delta \mu_i = \hat{\mu}_i - \mu_i^0\]
For numerical stability, the imposed chemical potential is expressed as a deviation from the solution phase’s pure endmembers (\(\mu_i^0\)). Similarly, the chemical potential of the solution phase is expressed as a deviation from ideal behavior:
\[\mu_i^{\mathrm{xs}} \equiv RTm \log{\gamma_i} = \mu_i -\mu_i^0 - RTm \log X_i\]
where \(\mu_i^{\mathrm{xs}}\) is the excess non-ideal component of the chemical potential, \(RTm\) is the thermal energy scale modified by the atomic site multiplicity \(m\), and \(\gamma_i\) are the activity coefficients. The equilibrium affinity condition for the solution phase is then given by:
\[ \phi_i \equiv \Delta \hat{\mu}_i - \mu^\mathrm{xs}_i = A + RTm\log X_i \]
where \(A\) is the saturation affinity of the phase relative to the imposed chemical potential, and \(RTm \log X_i\) is the entropic ideal mixing contribution to the Gibbs energy. To obtain the equilibrium properties of the solution phase, we need to find the composition \(X_i\) that satisfies the above expression for all components, along with the closure condition that the composition is normalized (\(\sum_i X_i = 1\)). This will simultaneously also provide the phase saturation affinity, indicating whether the phase is over-saturated (\(A>0\)), under-saturated (\(A<0\)), or simply stable (\(A=0\)). Note that for an equilibrium assemblage calculation, we know that all phases present will be saturated (with affinities that are zero within tolerance), with all other possible phases having unfavorable negative affinities.
Iterative composition/affinity update
While the above equation cannot be analytically solved directly for arbitrarily complex solution phases, it can be manipulated to allow a step-wise analytic solution update; this enables an iterative approach to achieve a simultaneous self-consistent solution for composition and affinity. By rearranging and solving for the unknown phase composition, we obtain:
\[ X_i = \exp{\left(\frac{\phi_i -A}{RTm}\right)}\\ \hat{X}_i \equiv \exp{\left(\frac{\phi_i}{RTm}\right)} \]
where \(\hat{x}_i\) is the unnormalized compositional vector, with a constant of proportionality that depends on the unknown phase affinity. In this unnormalized expression, approximate values are known for every term since the current estimation of the chemical potential offset \(\phi_i\) is readily calculated for the current composition estimate. The actual composition is then trivially determined by normalizing the compositional vector, directly producing updated estimates for both the phase composition and the saturation affinity:
\[ X_i = \hat{X}_i/(\sum_i \hat{X}_i)\\ A = RTm \log\left(\sum_i \hat{X}_i\right) \]
This approach yields mathematically equivalent results to the nested fraction procedure described in Ghiorso (2013), while remaining simpler and improving numerical stability and accuracy. The self-consistent composition and affinity for the phase is determined using the above procedure to iteratively update the phase composition, and using that updated composition to determine an updated chemical potential offset \(\phi_i\). This is repeated until both the saturation affinity and the molar composition for each component converge to within the desired tolerance. We find that the compositional convergence shows optimal performance when monitored in log-space, ensuring that both major and minor or trace components converge to the same fractional degree (rather than absolute molar fraction).
Rapid convergence using Gibbs curvature
The standard iterative update scheme for determining phase saturation affinities and compositions can be further improved by making use of derivative information to dramatically speed up the calculation. The method presented above adopts a conservative approach, where the current composition determines a chemical potential update, which in turn determines the new compositional update. In this procedure, the chemical potential evaluation is the computationally-limited step, especially for complex solution phases involving multiple ordering parameters (as is common for geologically important phases like the pyroxenes or spinels). We can speed up the process by also evaluating the Gibbs energy curvature matrix, \(\partial^2_n G = d\mu_i/dn_j\), which expresses the (linear) compositional-dependence of the chemical potentials. This curvature matrix provides rapid approximate chemical potential updates, allowing us to get closer to the self-consistent answer without additional costly chemical potential updates. To achieve this, we adopt a nested convergence structure where approximate convergence is achieved within a fast inner loop, limiting the number of required iterations for the computationally intensive outer loop that updates chemical potentials and the curvature matrix. This proves to be a worthwhile optimization in typical applications because the additional computational cost of calculating the derivative matrix is outweighed by substantial iteration reduction.
A naive application of this information would simply impose a linear update of the chemical potentials after each compositional update. Though this does converge, the gains are rather modest since it approximates the Gibbs energy function with a quadratic surface, and completely misses the logarithmic nature of mixing entropy. A better approach is instead to transform the curvature matrix into log-composition space (e.g., \({d\mu_i}/{d\log{n_j}} = n_j \cdot {d\mu_i}/{dn_j}\)), since this better captures the combined effect of ideal and non-ideal mixing. Theoretical justification for this can be found in the logarithmic terms for composition and activity coefficients shown in the chemical potential expression above. Unsurprisingly, the physically-rooted nature of this approximation yields dramatic improvements in terms of effectiveness and convergence speed.
Example computational benchmarks
To test this improved algorithm, we can use a variety of thermodynamic solution phases. This will be explored in detail in a future post. But for now, we can see the immediate impact by testing it on the MELTS liquid phase, which has the advantage of being a rather difficult test requiring simultaneous convergence in 15 different components. In this case, we find that MEXQAL performs extremely well when prototyped using Python code. Whereas the original saturation algorithm of Ghiorso (2013) takes ~280 ms to converge using a “cold-start” (with no initial compositional guess), MEXQAL converges in only 8 ms. These times primarily reflect the number of iterations required as well as the sizable overhead occupied by the objective-c/python bridge software required to evaluate the phase properties. We therefore expect than an implementation of the new algorithm relying directly on cython-wrapped C code will be then again more than an order of magnitude faster.
To be completely transparent, it is thus useful to consider the number of iterations required for convergence. Here, we focus on the typical use-case, termed a “warm-start,” where a reasonable initial compositional guess is provided and we are simply refining that guess in response to changing conditions. Using MEXQAL without any derivative information, the algorithm converges for a warm-start in ~12-15 iterations. Utilizing the curvature approximation, this decreases to only two iterations of the expensive outer convergence loop.
Future Benchmarking
In a future post, we will extensively benchmark MEXQAL for a variety of important geological phases and under different conditions. This will provide invaluable information on how to best optimize the code implementation. Additionally, we will revisit the question of whether it is worthwhile to translate the code into objective-C to be able to directly utilize it without the bridge-code intermediary, which currently occupies the majority of the runtime.