Phase coexistence in multicomponent systems is ubiquitous in nature and technology. Examples range from diverse purification processes in the chemical industry [1][2][3] to the formation of membraneless organelles via liquid-liquid phase separation in biology. [4][5][6][7][8][9][10][11] Thorough characterization of multicomponent phase equilibria involves not only identifying the phases present but also determining their composition and abundance, as the distribution and composition of species across phases significantly affect system properties and functions. This information guides processing methods and underlies calculations in processsimulation software. Critically, these calculations can constitute a substantial fraction of the overall computational time dedicated to simulation. 12 Current multiphase flash calculation schemes require knowledge of the number of equilibrium phases, 12,13 are sensitive to initial guesses, 14 and can converge to spurious or trivial solutions if root-finding is not appropriately bounded and constrained. 15,16 Therefore, efficient and accurate methods for predicting equilibrium states are valuable for both industrial applications and fundamental research.
At equilibrium, species distribute across phases based on the extremization of an appropriate thermodynamic potential. For example, minimization of the Gibbs energy dictates equilibrium for a system at specified temperature T , pressure p, and global composition x i .
Equilibrium phase-coexistence arises when species partition into distinct phases with equal chemical potentials driven by the extremization, rather than forming a homogeneous mixture. For a system at fixed T and p, this yields
where µ π i is the chemical potential of species i ∈ {A, B, C, . . .} in phase π ∈ {α, β} with composition {x j } π = {x π A , x π B , x π C , . . . }. Eq. ( 1) constrains the equilibrium state of the system, as manifest in Gibbs' phase rule (i.e. F = N -P + 2 where F is the number of independent intensive relationships needed to specify a system of N species and P phases). Provided a thermodynamic model for describing the mixture behavior as a function of intensive variables, Eq. ( 1), or its equivalent at other conditions, functionally comprises N (P -1) equations to solve for the compositions of the various phases, with others fixed. The complexity of identifying equilibrium states can vary, even while the underlying thermodynamic framework is straightforward.
Determining the conditions, expected phases, and the chemical nature of species usually depends on appropriately parameterized equations-of-state or available free-energy models.
For condensed phases and binary mixtures, there are several simple free-energy models like the Margules equations, 17 the van Laar model, 18 or the Guggenheim-Scatchard/Redlich-Kister equation. 19 More complex models such as the Wilson models, non-random two-liquid (NRTL) models, 20 universal quasi-chemical theory (UNIQUAC), 21 UNIQUAC Functionalgroup Activity Coefficients (UNIFAC) models, 13,22,23 Flory-Huggins theory 24 can treat multicomponent systems. Although increasing complexity of the free-energy model or equationof-state can facilitate more accurate representation of physical systems, the underlying calculations and theoretical principles for phase behavior remain the same for simple and complex models alike.
Given a thermodynamic model, calculating phase stability and equilibrium compositions can be approached in various ways. Simple models and binary mixtures may yield algebraic relationships that can be handled analytically or resolved using simple numerical schemes, such as self-consistent iteration or Newton's method. However, characterizing multicomponent phase coexistence typically requires dedicated software and more sophisticated numerical algorithms. Many algorithms are designed to work for only a specific set or number of phases. 25 Direct solution methods based on Newton's root-finding algorithm can be effective but are computationally intensive and sensitive to the initial seed. Jindrova et al. refined Newton's algorithm and a successive substitution strategy to locate roots. Additionally, Nichita, [26][27][28] Jindrova, 14,16 and Castier 29 independently performed volume stability analysis to obtain better initial guesses for the substitution strategy. There has been significant development in generating phase diagrams using constrained backmapping search algorithms. 15,[30][31][32][33] Indirect solution methods, based on thermodynamic principles and geometric criteria established via stability analysis, offer alternative approaches. Examples include Korteweg's tangent construction 34 and Binous and Bellagi's arc extension method. 35 Michelsen's multiphase flash algorithm 36 minimizes the distance between the tangent plane and the free energy surface to identify coexisting phases. Homotopy methods have also been used to calculate critical and saturation properties of mixtures. 37,38 Additionally, Mao et al. 39 generalized phase-diagram construction to multicomponent systems using a convex-hull construction 40 applied to a discretized free-energy manifold, although accuracy and memory requirements depend on the mesh size. Overall, there is a need for simple, generalizable, and efficient methods for phase-coexistence calculations.
Machine learning (ML) techniques facilitate phase-coexistence calculations, offering prospective advantages relating to time-and memory-efficiency relative to more traditional optimization strategies. [41][42][43][44][45][46][47][48][49][50] However, many efforts only address the issue of phase stability and neglect consideration of phase composition. [42][43][44][45][46] Others have been restricted to binary systems with limited demonstration of more complex mixtures. [47][48][49][50] Recently, Flory-Huggins (FH) theory has been combined with ML to improve the interpretability and accuracy of mixture behavior predictions, but limitations exist in their ability to handle complex interactions and multicomponent systems beyond binary mixtures. 49,50 Nevertheless, such works highlight the potential of ML as part of a generalizable, accurate, efficient, and extensible framework for characterizing multicomponent phase behavior.
Here, we describe a data-driven workflow to characterize the phase behavior of multicomponent systems. Figure 1 illustrates the overall approach in the context of ternary systems described by Flory-Huggins (FH) theory. Using FH theory as a representative freeenergy model, we construct a series of phase diagrams across the model parameter space using labor-intensive methods. This data is then used to develop an ML surrogate model, based on neural network architectures, to predict the number, composition, and relative abundance of equilibrium phases from model parameters and total system composition.
Surrogate models optimized with and without physics-informed architectures and loss functions are compared. Errors are assessed for classification (number of equilibrium phases) and regression (composition and abundance of phases). Predictions from the surrogate model, which are computationally efficient and improvable, are then used to warm-start a simple optimization to precisely and accurately characterize the system's phase behavior. This procedure exemplifies an efficient, accurate, and extensible approach to phasecoexistence calculations.
Figure 1: Strategy for multi-component phase-coexistence prediction using machine learning. 1,036 ternary phase diagrams are generated using the algorithm arc continuation algorithm (Equations 6-8) and convex hull construction algorithm, 39 and are used as training data for a physics-informed machine learning (ML) model to classify phase regions and predict equilibrium phase compositions. The ML predictions serve as initial guesses for the Newton-CG method to obtain equilibrium composition predictions.
For demonstration, we consider the thermodynamics of ternary systems described by FH solution theory. Systems are comprised of species A, B, and C that occupy a lattice of n sites with volume V = nv 0 . The species can possess size disparities, reflected in their molar volumes v i . For a polymer comprised of N i monomers that each occupy a single lattice site,
Systems are incompressible such that V = i n i v i where n i is the mole number for species i. System composition is specified by the volume fractions
The dimensionless, intensive (per lattice site) Helmholtz energy of mixing follows as
where β = (k B T ) -1 is the inverse temperature with k B as Boltzmann's constant, and χ ij is the Flory-Huggins interaction parameter for species i and j with χ ii = 0; the summations are over all components (A, B, C). Altogether, the behavior of a system is determined by the composition ϕ = (ϕ A , ϕ B , ϕ C ), the molar volumes of the species v = (v A , v B , v C ), and the interaction parameters X = (χ AB , χ AC , χ BC ).
Up to a constant, chemical potentials are obtained by partial differentiation of the Helmholtz energy of mixing:
Using Eq. (2) in Eq. (3), this yields
where the summations exclude the species for which the chemical potential is being assessed (e.g., for µ A , the summation for j ̸ = i is equivalent to that for j ∈ {B, C}).
The thermodynamic stability of a mixture is assessed by considering the determinant of the Hessian matrix for the Helmholtz energy. For
the spinodal boundary of a ternary mixture is the locus of all compositions that solve
Critical points are identified by additionally considering constraints on third-order derivatives 51,52 given by
and
For fixed total particle density (ρ = n/V ) and constant temperature, Gibbs' phase rules indicate there can be at most three coexisting phases. To characterize three-phase coexistence, there are 12 variables. Nine correspond to the volume fractions in each phase: ϕ α , ϕ β , and ϕ γ , for which each ϕ π = (ϕ π A , ϕ π B , ϕ π C ). Three correspond to the fractional abundances of each phase-w α , w β , and w γ ). Criteria for chemical equilibrium applied to each species across each phase
provide six independent equations. For a system with a specified total composition, material balance constraints provide the remaining equations:
To characterize two-phase coexistence, the variable count is reduced to eight, with commensurate reduction by three equations from Eq. ( 9) and one equation from Eq. ( 10). Section 2.2 describes algorithmic approaches for determining equilibrium compositions.
Two different algorithms are used to characterize phase-coexistence based on the principles outlined in Section 2.1. For systems with at least one critical point and two-phase coexistence, an iterative and perturbative approach based on natural parameter continuation (NPC) is used to construct binodal curves originating from a critical point. Otherwise, the approach described by Mao et al. 39 based on convex hull construction (CHC) is used.
NPC is straightforward and computationally efficient but limited, while CHC is general but computationally intensive. Nevertheless, with this combination, the equilibrium composition of phases at coexistence can be reliably determined for models described by FH solution theory. The following algorithms are thus used to provide ground-truth results and requisite training data for the development of ML models (Section 2.3).
For two-phase coexistence, Eq. ( 9) is rearranged as
Provided a point on the coexistence curve ϕ * , a nearby point can be identified by solving a set of linear equations that enforce Eq. ( 12) following a small perturbation in the composi-tion:
j∈{A,B} π
where δϕ π j is the small perturbation in the composition of species j in phase π. Coexistence curves (i.e., a locus of equilibrium composition tuples) can then be constructed as follows:
1. Define tolerance parameters δ ∅ and δ 0 .
2. Identify and set the critical point to be ϕ * .
3. Generate a random, small perturbation on the composition δϕ, yielding two new compositions:
Use the compositions ϕ ′ and ϕ ′′ as initial guesses to solve Eq. ( 12), producing coexisting compositions ϕ α new and ϕ β new that are distinct from ϕ * .
5. Set ϕ π old ← ϕ π new and and use for Eq. ( 13). Set one of the δϕ π i (e.g., ϕ β B ) to a small perturbation and solve for the remaining δϕ π i to produce δϕ α ′ and δϕ β ′ .
6. Set ϕ ′ = ϕ α(0) + δϕ α ′ and ϕ ′′ = ϕ β(1) + δϕ β ′ and use as initial guesses to solve Eq. ( 12), producing new coexisting compositions ϕ α new and ϕ β new that are those distinct from those prior.
7. Repeat steps 5 and 6 until either ||ϕ α new -ϕ β new || < δ ∅ , which indicates a closure of the coexistence curves, or when any ϕ π i < δ 0 , which indicates termination at a composition boundary.
For the calculations described in this paper, δ ∅ = δ 0 = 10 -9 . Initial trials for random composition perturbations are set to have a magnitude of 10 -6 . Equations are solved numerically using fsolve from Python's SciPy module. Occasionally, the trial perturbations resulted in solutions that collapsed back to the critical point or other prior generated points, in which case new perturbations would be attempted with possibly different magnitudes.
For systems without critical points or valid coexistence curves extending from critical points, the utility of the NPC algorithm is limited. In such cases, we use the CHC to identify equilibrium compositions. On a free energy surface, compositions with equal chemical potential are cotangent, while stable compositions (D > 0) that are not cotangent with any other points exist as single phases. This information can be accurately reconstructed by creating a convex hull of the free energy surface and projecting it onto the composition space. We briefly remark on salient aspects of the algorithm as applied to a ternary system, but readers are referred to the work of Mao et al. 39 for a more complete description.
The composition space (ϕ A , ϕ B ) is discretized into a mesh of equilateral triangles, or twodimensional simplices. Using a finer mesh results in more accurate calculations but also increases computational cost and memory requirements; this work uses a simplex edge-length of 0.0002. After generation of the mesh, the free energy surface (FES) is also discretized into points defined by the tuple
the FES is calculated using the Quickhull algorithm. 53 The convex hull of a non-convex FES will necessarily deform the original simplices and facilitate the identification of cotangent points on the FES. If one of the projected simplices has three unstretched sides (maximum edge length within five times initial mesh size 39 ), the system is homogeneous (no phaseseparation). If two sides are stretched (side length greater than five times the initial mesh size), the two farthest vertices are cotangent, indicating two coexisting phases. If all three sides are stretched, the three vertices of the simplex are cotangent, indicating three coexisting phases. With graph theoretic techniques, the number of equilibrium phases and their compositions can be determined.
We explore machine learning algorithms as computationally expedient and generalizable alternatives to more traditional approaches for characterizing phase coexistence of multicomponent systems. Neural network architectures, with and without physics-informed loss functions, are optimized using data generated by the algorithms described in Section 2.2. The performance of the ML models is evaluated based on predicting the number of coexisting phases, their compositions, and relative abundance for FH models not featured in training data.
Dataset description. The dataset in this work is comprised of 1,036 phase diagrams:
107 diagrams (10%) with no phase separation (one phase), 538 diagrams (52%) with up to two-phase coexistence, and 391 diagrams (38%) with up to three-phase coexistence. Each phase diagram is produced using the methods of Section 2.2 with a distinct parameter set:
Parameters for the models are each selected from the range v i ∈ [1, 3] and χ ij ∈ [1, 3] where values for both ranges are discretized with a resolution of 0.1. Let s denote a parameter set and U denote the set of all possible parameter sets. With the given discretization, the total membership of U is then |U | = 21 6 . Initially, 750 possible parameter sets are randomly selected from U with uniform probability to form S ⊂ U ; care is taken to ensure that all parameter sets from this sampling are unique. From this initial sampling, only around 6.6% (≈ 50) of the selected parameter sets yielded three-phase coexistence. Using these parameter sets to define T ⊂ S, the representation of such rare systems is augmented by generating six additional parameter sets for each parameter set t ∈ T . Each new parameter set t ′ is generated from t by adding a Gaussian random vector X. In particular, we use t ′ = t + X with X ∼ N (0, σ 2 I) where σ = 0.005. All t ′ that yielded three-phase coexistence are collected and added to S, resulting in a final membership of |S| = 1, 036 parameter sets.
Input and output labels for the dataset are then generated as follows. First, the composition-space of the mixture is discretized into a uniform mesh with resolution 10 -4 .
For each parameter set, if there are more than 1,000 single-phase simplices, the centroid of 1,000 randomly chosen simplices is added to the database; otherwise, the centroid of all single-phase simplices is added. For double-phase simplices, if there are more than 1,000, a random point between the ends of 1,000 randomly chosen double-phase separations is generated; otherwise, a random point between each double-phase separation is added. For multiple three-phase separations, a uniform number of points is generated in each region, ensuring a total of 1,000 three-phase points in the database. Since the number of single and double-phase simplices are determined by the size of the discretization mesh, the number of data points per each parameter set can vary.
For each tuple (ϕ A , ϕ B ) the number of equilibrium phases, their compositions, and their abundances are recorded.
In this fashion, we define an input vector x =
that is linked to two outputs. The first output is a one-hot encoded classification vector y c ∈ R 3 , for which a nonzero entry indicates the presence of one, two, or three phases at equilibrium. The second output is a vector y r ≡
, which describes the composition and abundances of the equilibrium phases. The phases are ordered such that ϕ α A has the minimum value
If two phases have the same ϕ A , they are further ordered according to ϕ B . Such an ordering ensures a consistent representation of the equilibrium phases.
For systems with a single phase, ϕ α A , ϕ α B match the inputs ϕ A , ϕ B , and w α is set to unity; the abundance entries for phases β and γ are set to zero. However, the composition abundance entries for phases β and γ are assigned a value of 1/3. The value 1/3 is chosen to distribute errors uniformly across species. The absolute composition of these species in equilibrium will be determined by the abundance of the respective phases. For systems with two equilibrium phases, entries for the third phase compositions (i.e., ϕ γ A , ϕ γ B ) are set to 1/3, and the abundance w γ is set to zero. Model architectures. Figure 2 summarizes model architectures in this study; all models are implemented using PyTorch. 54 Every model takes as input x and predicts two outputs: ŷc ∈ R 3 and ŷr ∈ R 9 . The output ŷc is a classification vector containing the probabilities that x yields one, two, or three phases. The output ŷr is a regression vector with entries associated with predicted equilibrium compositions and weights (in the same order as for y r ).
The basic model architecture consists of three, fully-connected hidden layers, each with m (tunable hyperparameter) hidden units; this yields a hidden vector h ∈ R m . This hidden vector is then passed through a "classification layer" with softmax activation to yield ŷc . This vector is also fed into separate "regression layers" to predict the composition (ϕ α , ϕ β , ϕ γ ) and abundance (w). Each regression layer consists of three hidden units, representing the composition of A, B, and C for each phase, and the abundance of α, β, and γ phases. Sigmoid activation is applied to limit predicted values on compositions and abundances to be between zero and unity, which avoids obviously unphysical values; however, overall composition and abundance constraints are not enforced. Since the composition of C depends on A and B, only the predictions for A and B compositions are kept and combined with the abundance predictions to form ŷr .
We also consider a variation on the basic model architecture that enforces consistency between ŷr and the majority class featured in ŷc . This is achieved using a mask-layer that sets abundance entries in ŷr to zero based on the plurality class indicated in ŷc (see dashed box in Figure 2). For example, if one equilibrium phase is predicted, then abundance entries associated with the β and γ phase are set to zero. If two equilibrium phases are predicted, then entries associated with the γ phase are set to zero. If three equilibrium phases are predicted, then ŷr is preserved from the regression layer. Compositions of species for non-existent phases are set to 1/3, as described earlier. To enforce constraints on overall composition and abundance, the physics-informed (PI) model incorporates softmax activation functions, ensuring that predicted phase compositions and abundances sum to unity. As an alternative approach, following the masking, softmax normalization is applied to the abundances to en-sure their sum equals unity.
Loss functions. Models are optimized using loss functions that target raw numerical accuracy as well as physical sensibility. Both the simple baseline and PI models optimize a composite loss function
that combines losses for classification cross-entropy (CE), L CE , and regression mean absolute error (MAE) loss, L MAE . The weighting parameter (λ CE = 0.1) is determined empirically to balance loss magnitudes throughout training. While a perfect and physically meaningful model would necessarily minimize L base , with data limitations, simply minimizing the baseline loss function may not strictly satisfy all the criteria prescribed for thermodynamic systems at equilibrium. We therefore also consider augmented PI models (referred to as PI+) optimized with a composite loss function that includes additional regression targets
where λ split = 0.01, λ ∆µ = 0.01, and λ f = 0.001 (identical across all models) are weighting parameters chosen through empirical tuning. These values balance the influence of the PI losses to focus on minimizing L base while incorporating physical constraints. The specific functional forms for these PI losses are described next.
In Eq. ( 15), the additional loss terms aim to satisfy different constraints on the thermodynamics of physical systems. In particular, L split relates to constraints on the total composition of a given species distributed across phases:
Figure 2: Architectures of the machine learning models. In both the baseline and physicsinformed (PI) models, the parameter vector x is fed into the model to produce an intermediate hidden vector h. The hidden vector h produces two outputs: (1) a phase classification probability vector ŷc , trained with cross-entropy (CE) loss, and (2) an equilibrium composition and abundance vector ŷr , trained with mean absolute error (MAE) loss. Softmax activation is applied in PI models to ensure that the equilibrium composition and abundance vectors sum to unity. Optional functionalities (indicated by dashed lines and boxes) include a "mask", activated based on ŷc , which sets corresponding elements in ŷr to zero if an input is classified as one-or two-phase. For PI models, softmax renormalization is applied to the masked abundance to ensure the sum equals unity. Additionally, in PI models, a PI composite loss can be incorporated alongside the MAE loss during training for ŷr prediction.
The second term L ∆µ relates to the condition of equal chemical potentials for species across coexisting equilibrium phases. This loss is calculated as
where ∆µ ππ ′ i is as defined in Eq. ( 12), and the first two summations are over the groundtruth equilibrium coexisting phases (i.e., π, π ′ ∈ {α, β} for two-phase coexistence and π, π ′ ∈ {α, β, γ} for three-phase coexistence). The additional term L f promotes the minimization of the free energy of the equilibrium system:
We acknowledge there are various reasonable ways to constraint losses for physical constraints; the current work examines the overall strategy of incorporating physical information into the ML workflow rather than identifying optimal implementations.
Model training and assessment. To assess model generalizability and mitigate selection bias on test data, a nested five-fold cross-validation (CV) procedure is used. Stratified sampling is employed to evenly distribute diagrams featuring one, two, and three phases across the five folds. Then, five iterations are performed in a process referred to as the outer CV. Each iteration uses a unique fold as the test set and the remaining four folds as the overall training set to provide a more robust assessment of model performance.
The overall training set is further divided into training and validation sets, using a similar five-fold CV approach (inner CV) to the outer CV process. Each fold of the inner CV is trained with 10% of the training data for efficient hyperparameter optimization. Tunable hyperparameters include batch sizes of {5000, 10000, 20000}, learning rates of {0.001, 0.005, 0.01}, the presence (or absence) of a mask, and the number of neurons selected from {64, 128, 256} for each hidden layer. The optimal hyperparameter setting for each fold is identified by the highest average validation composite score across five sub-folds, calculated as the sum of the F 1 score for classification and the average R 2 score.
Each fold of the outer CV uses the optimal hyperparameter settings identified from its corresponding five-fold inner CV and retrains the model for up to 500 epochs.
The retraining involves selecting 80% of the overall training diagrams as the training set and 20% as the validation set, using the same stratified splitting. During the retraining process, the impact of training data sizes on model performance is assessed by using 1%, 5%, 10%, 20%, 30%, and up to 100% of the training data. Because each diagram contains a different number of data points, the number of training, validation, and test set data points ranges from 1,447,511 to 1,458,241, from 358,987 to 366,470, and from 452,308 to 460,074.
The nested CV strategy yields a mean and standard deviation of F 1 and R 2 scores as determined from the five-fold outer CV test sets. Given the imbalanced phase distribution in the dataset, the F 1 score evaluates classification performance, while the R 2 score assesses regression accuracy for the variables in y r .
We implement a post-inference optimization procedure to correct some deficiencies in ML model predictions. This procedure uses the predictions from the ML model as a warmstart on initial values for more traditional optimization algorithms (e.g., truncated Newton method). The objective function for minimization is
where 1 c is an indicator function equal to unity when the condition c is satisfied and zero otherwise and
For both Eq. ( 19) and ( 20), the first two summations are over predicted equilibrium coexisting phases (i.e., π, π ′ ∈ {α} for a single equilibrium phase, π, π ′ ∈ {α, β} for two-phase coexistence, and π, π ′ ∈ {α, β, γ} for three-phase coexistence). Eq. ( 20) is specifically relevant for two-phase coexistence and is minimized when the tie-line composition vectors are collinear and oriented in opposite directions. The indicator function 1 w α <1 excludes computation of L col when there is only a single predicted phase, as this term would otherwise diverge. In fact, this procedure has no effect when only a single equilibrium phase is predicted. Relatedly, we note that this algorithm is asymmetrically robust against erroneous misclassification of the system phase behavior. If the predicted number of phases exceeds the true number of phases, then converged solutions will "collapse" compositions onto those of the true equilibrium phases. However, this procedure will not identify the true solution if the predicted number of phases is fewer than the ground-truth number.
The final optimization employs the Newton-CG optimizer in scipy module in Python.
The Jacobian and Hessian matrix for the objective function are computed using the autograd package through automatic differentiation. The maximum number of iterations for optimization is limited to 10,000. If newly optimized compositions are within a tolerance of 10 -7 of the ideal value of the objective function, these values replace the predictions proffered by the ML model. Optimizations are only considered successful if they satisfy the stability criterion |H f | > 0 (see Eq. ( 6)).
With the standard loss functions (i.e. L CE and L MAE ) and a basic architecture, the ML model predicts phase separation and equilibrium compositions reasonably well. Figures 3a and 3b qualitatively depict performance in both classification and regression for some representative phase diagrams. Figures 3c and 3d quantitatively summarize results across all phase diagrams.
The ML model capably predicts the number of phases at equilibrium with an overall accuracy rate of about 97% (Figure 3c). The primary source of error (4.6%) stems from misclassifying three-phase points as two-phase, which is attributed to the relative paucity of three-phase splits (only 17% of the total data). Additionally, a significant number of twophase points (2.2%) are misclassified as one-phase. A closer inspection of predicted phase diagrams indicates that misclassifications mostly occur near binodals.
The model also performs well in predicting phase abundances and their compositions (Figure 3d). By inspection, there are vertical error regions in the predicted abundance for all phases at true abundances of 0 and 1. These errors stem from inaccurate predictions of equilibrium abundance for non-existent phases, such as phases β and γ in a one-phase region and phase γ in a two-phase region. This misprediction also leads to similar error regions for equilibrium compositions around 1/3.
Table 1 provides the baseline expectations for a standard ML model. It highlights nuances in regression performance across different phase regions. The single-phase region has the lowest average MAE (0.006), followed by the two-phase (0.023) and three-phase (0.037) regions. This trend suggests increasing difficulty in predicting compositions as the number of coexisting phases increases. Notably, all R 2 values remain high across all phases, with the three-phase region exhibiting a value above 0.88.
To build on the prior model, we evaluate the potential of incorporating additional physical information on prediction accuracy (Figure 4a). In particular, physical constraints on overall composition and abundance, along with several physics-informed losses (detailed in Section 2.3) are implemented, and classification masks are used to zero the abundances of non-existent phases. While the baseline, PI, and PI+ models without classification masks achieve comparable F 1 and R 2 scores, models with masks significantly underperform in equilibrium composition prediction (Figure 4b).
The baseline model exhibits lower accuracy compared to the PI and PI+ models in onephase and two-phase regions, while its performance is comparable to other models in threephase scenarios (Table 1). The PI and PI+ models show similar performance across all scenarios under these metrics. The coexistence curve predictions of both the PI and PI+ models are similar (Figures 5a, 6, and S6), producing smooth and physically sensible two-and three-phase coexistence curves. In contrast, the baseline model generates erratic two-phase coexistence curves that significantly deviate from the true curves. This is also evident from the distribution of the MAE loss in Figure 5b, where the baseline model (red) has a higher average MAE than the PI (blue) and PI+ (green) models, which perform similarly. The unphysical coexistence curves of the baseline model highlight the limitations of using broad performance metrics to assess improvements in predictive accuracy. Errors are better resolved by examining deviations in chemical equilibrium potential (∆µ) and split loss (L split ), where the baseline model shows significantly larger errors than both the PI and PI+ models.
To further examine the impact on composition and abundance constraints, we analyze two additional metrics: L unity , which relates to the volume fractions of each species within a given phase, and L weight , which measures overall material conservation. These metrics are defined as:
and
Since the PI and PI+ models enforce unity in composition and abundance through the softmax activation function, L unity and L weight remain zero for these models, whereas the baseline model violates these constraints (Figure 5b). The PI+ model, trained with additional constraints, demonstrates smaller deviations in chemical equilibrium potential (∆µ) and marginally improves composition prediction, split loss, and free energy minimization loss compared to the PI model. Overall, designing a physics-informed model architecture to enforce material constraints is essential; however, the addition of extra losses or masks complicates training without yielding significant improvements in phase classification or equilibrium composition prediction. Therefore, the PI architecture, without additional losses, emerges as the best practical choice for implementation. Post-ML optimization error analysis was performed on a random sample of 187 twophase and 76 three-phase equilibrium phase diagrams, using models trained on fold 1 data (Tables 2 and S1). The results indicate that Newton-CG optimization, initialized with predictions from ML models, achieves near-perfect success rates and significantly reduces deviations from true equilibrium compositions compared to individual ML model predictions.
After post-ML optimization, the PI model trained on the full dataset outperformed both the baseline and PI+ models in predicting two-phase and three-phase coexistence. The PI+ model trained on the full fold 1 dataset exhibited the lowest initial ML prediction accuracy, leading to the worst post-ML optimization performance. This underscores the importance of initial ML prediction accuracy in determining the effectiveness of the post-ML optimization process. Interestingly, all models achieved high success rates and minimal prediction errors in two-phase scenarios, even when trained with limited data. The PI+ model trained on 10% of the data demonstrated the highest accuracy in two-phase predictions. However, predicting three-phase equilibrium compositions remained challenging, with difficulties persisting even after post-ML optimization. This disparity likely stems from the relative scarcity of three-phase coexistence in the training set, which increases complexity and complicates precise prediction.
The post-ML optimization process is also efficient and parallelizable -taking less than 1 second to converge to the optimal solution (Table S1). The ML model training requires less than 1 MB for 140,000 parameters, a substantial reduction in memory usage compared to arc continuation or convex hull methods, which demand approximately 1 GB of storage and 50 GB of memory per run. Overall, the accurate ML predictions of equilibrium compositions enable rapid convergence to highly accurate solutions, offering significant advantages in both memory-and time-efficiency.
In this work, we presented an efficient and extensible machine learning-based approach for calculating phase coexistence in ternary systems. A neural network trained on phase coexistence data was able to predict the number and compositions of equilibrium phases for a solution prepared at a given composition under a specific mixing potential. Incorporating physical constraints into the neural network architecture enhanced prediction accuracy, while additional physics-informed losses offered no significant improvement. The physicsconstrained architecture produced higher-quality models with less data, offering advantages in scenarios where data acquisition is labor-or resource-intensive. However, the resulting models still exhibit errors that may be unacceptable for certain applications, such as process simulation software. To achieve precise results, a Newton conjugate gradient method was used, with machine-learning predictions serving as a warm start for optimization to determine final equilibrium phase compositions. This integration of neural networks with numerical refinement enabled rapid and accurate predictions of coexisting phases, their compositions, and abundances.
This work motivates several areas of future inquiry. Extensions to systems with more components would increase utility for complex industrial and biological processes. Expanding beyond the Flory-Huggins theory by incorporating other free energy models or data from molecular simulation, perhaps in a single framework, would further enhance its generalizability across diverse chemical systems. Additionally, exploring more advanced physics-informed learning strategies, incorporating uncertainty quantification, and refining neural network architectures could boost prediction efficiency and reliability. Collectively, these directions could enhance both the theoretical and practical impact of leveraging ML for phase coexistence calculations.