Ternary Stability and LLE
*************************

This example demonstrates "STABILITY" and "LLE" calculations for the
ternary system of water/ethanol/benzene to identify stable and
unstable feed compositions and to trace the liquid-liquid phase
boundary and tie lines.

In AMS2026, the python package "ternary" is available in "amspython"
and can therefore be used directly for plotting ternary diagrams.


Python code for generating the input dataframes
===============================================

The helper script "TernaryStabilityLLE_common.py" is used to generate
the dataframes for the later calculations. In particular, it builds
the ternary composition grid used for the "STABILITY" calculations and
the composition scan used for the "LLE" calculations. The same script
also contains the shared plotting routine used afterward by both
examples.

   from dataclasses import dataclass
   from decimal import Decimal
   from typing import List, Tuple

   import matplotlib.pyplot as plt
   import numpy as np
   import pandas as pd


   @dataclass
   class Compound:
       name: str
       coskf: str


   def build_simplex_grid_points(
       n_components: int,
       min_frac: float = 0.1,
       step: float = 0.1,
   ) -> List[Tuple[float, ...]]:
       min_dec = Decimal(str(min_frac))
       step_dec = Decimal(str(step))

       if n_components < 2:
           raise ValueError("n_components must be >= 2.")
       if min_dec <= 0 or step_dec <= 0:
           raise ValueError("min_frac and step must be > 0.")
       if n_components * min_dec > 1:
           raise ValueError("n_components * min_frac must be <= 1.")

       precision = max(-min_dec.as_tuple().exponent, -step_dec.as_tuple().exponent)
       scale = 10**precision
       total_units = scale
       min_units = int(min_dec * scale)
       step_units = int(step_dec * scale)

       if Decimal(min_units) / scale != min_dec:
           raise ValueError("min_frac is incompatible with the chosen decimal precision.")
       if Decimal(step_units) / scale != step_dec:
           raise ValueError("step is incompatible with the chosen decimal precision.")

       def recurse(remaining_components: int, remaining_units: int) -> List[Tuple[int, ...]]:
           if remaining_components == 1:
               if remaining_units >= min_units and remaining_units % step_units == 0:
                   return [(remaining_units,)]
               return []

           points: List[Tuple[int, ...]] = []
           max_units = remaining_units - min_units * (remaining_components - 1)
           for units in range(min_units, max_units + 1, step_units):
               for tail in recurse(remaining_components - 1, remaining_units - units):
                   points.append((units,) + tail)
           return points

       integer_points = recurse(n_components, total_units)
       return [tuple(float(x / scale) for x in point) for point in integer_points]


   def build_ternary_grid_df(
       compounds: List[Compound],
       min_frac: float = 0.1,
       step: float = 0.1,
       temperature: float = 298.15,
       method: str = "COSMORS",
   ) -> pd.DataFrame:
       rows = []
       for mix_idx, point in enumerate(build_simplex_grid_points(3, min_frac=min_frac, step=step)):
           rows.append(
               {
                   "Method": method,
                   "Temperature": temperature,
                   "ncomp": 3,
                   "sys_idx": 0,
                   "mix_idx": mix_idx,
                   "s1": compounds[0].name,
                   "s2": compounds[1].name,
                   "s3": compounds[2].name,
                   "x1": point[0],
                   "x2": point[1],
                   "x3": point[2],
                   "unstable": pd.NA,
                   "converged": pd.NA,
                   "llle_detected": pd.NA,
                   "xI1": pd.NA,
                   "xI2": pd.NA,
                   "xI3": pd.NA,
                   "xII1": pd.NA,
                   "xII2": pd.NA,
                   "xII3": pd.NA,
               }
           )
       return pd.DataFrame(rows)


   def build_ternary_lle_df(
       compounds: List[Compound],
       s1_ratio: float,
       s2_ratio: float,
       s3_min: float,
       s3_max: float,
       num_points: int,
       temperature: float = 298.15,
       method: str = "COSMORS",
   ) -> pd.DataFrame:
       total = 1.0

       if s1_ratio <= 0 or s2_ratio <= 0:
           raise ValueError("s1_ratio and s2_ratio must be > 0.")
       if num_points < 2:
           raise ValueError("num_points must be >= 2.")
       if s3_min < 0 or s3_max > total or s3_min > s3_max:
           raise ValueError("Require 0 <= s3_min <= s3_max <= 1.")

       ratio_total = s1_ratio + s2_ratio
       s1_share = s1_ratio / ratio_total
       s2_share = s2_ratio / ratio_total

       rows = []
       for mix_idx, x3 in enumerate(np.linspace(s3_min, s3_max, num_points)):
           remaining = total - x3
           x1 = remaining * s1_share
           x2 = remaining * s2_share

           if any(value < -1e-12 or value > total + 1e-12 for value in (x1, x2, x3)):
               raise ValueError("Generated composition is outside the [0, 1] interval.")

           rows.append(
               {
                   "Method": method,
                   "Temperature": temperature,
                   "ncomp": 3,
                   "sys_idx": 0,
                   "mix_idx": mix_idx,
                   "s1": compounds[0].name,
                   "s2": compounds[1].name,
                   "s3": compounds[2].name,
                   "x1": float(x1),
                   "x2": float(x2),
                   "x3": float(x3),
                   "unstable": pd.NA,
                   "converged": pd.NA,
                   "llle_detected": pd.NA,
                   "xI1": pd.NA,
                   "xI2": pd.NA,
                   "xI3": pd.NA,
                   "xII1": pd.NA,
                   "xII2": pd.NA,
                   "xII3": pd.NA,
               }
           )

       return pd.DataFrame(rows)


   def _collect_feed_points(df: pd.DataFrame) -> List[Tuple[float, float, float]]:
       points: List[Tuple[float, float, float]] = []

       for row in df.itertuples():
           points.append((row.x1, row.x2, row.x3))

       return points


   def _collect_phase_points(df: pd.DataFrame, prefix: str) -> List[Tuple[float, float, float]]:
       points: List[Tuple[float, float, float]] = []

       for row in df.itertuples():
           points.append((getattr(row, f"{prefix}1"), getattr(row, f"{prefix}2"), getattr(row, f"{prefix}3")))

       return points


   def _collect_tie_lines(df: pd.DataFrame) -> List[Tuple[Tuple[float, float, float], Tuple[float, float, float]]]:
       tie_lines: List[Tuple[Tuple[float, float, float], Tuple[float, float, float]]] = []

       for row in df.itertuples():
           tie_lines.append(((row.xI1, row.xI2, row.xI3), (row.xII1, row.xII2, row.xII3)))

       return tie_lines


   def configure_tax(tax, labels: List[str], title: str) -> None:
       fontsize = 10
       tax.boundary(linewidth=2.0)
       tax.gridlines(color="black", multiple=0.05, linewidth=0.5, alpha=0.4)
       tax.bottom_axis_label(f"{labels[0]}(1)", fontsize=fontsize, offset=0.1)
       tax.right_axis_label(f"{labels[1]}(2)", fontsize=fontsize, offset=0.1)
       tax.left_axis_label(f"{labels[2]}(3)", fontsize=fontsize, offset=0.1)
       tax.ticks(
           axis="lbr",
           multiple=0.1,
           linewidth=1,
           tick_formats="%.1f",
           offset=0.015,
           fontsize=fontsize,
       )
       tax.clear_matplotlib_ticks()
       raw_ax = tax.get_axes()
       for spine in raw_ax.spines.values():
           spine.set_visible(False)
       raw_ax.set_xticks([])
       raw_ax.set_yticks([])
       tax.set_title(title)


   def plot_demo_figure(df_stability: pd.DataFrame, df_lle: pd.DataFrame):
       import ternary

       labels = [
           str(df_stability.iloc[0]["s1"]),
           str(df_stability.iloc[0]["s2"]),
           str(df_stability.iloc[0]["s3"]),
       ]
       fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

       _, tax_stab = ternary.figure(scale=1.0, ax=axes[0])
       configure_tax(tax_stab, labels, "STABILITY: stable vs unstable grid points")

       stable_df = df_stability.loc[~df_stability["unstable"].fillna(False).astype(bool)].copy()
       unstable_df = df_stability.loc[df_stability["unstable"].fillna(False).astype(bool)].copy()

       stable_points = _collect_feed_points(stable_df)
       unstable_points = _collect_feed_points(unstable_df)

       if stable_points:
           tax_stab.scatter(
               stable_points,
               s=42,
               marker="o",
               facecolors="none",
               edgecolors="#1b5e20",
               linewidths=1.2,
               label="Stable feed",
           )
       if unstable_points:
           tax_stab.scatter(
               unstable_points,
               s=42,
               marker="o",
               facecolors="none",
               edgecolors="#c04a00",
               linewidths=1.2,
               label="Unstable feed",
           )

       _, tax_lle = ternary.figure(scale=1.0, ax=axes[1])
       configure_tax(tax_lle, labels, "LLE: phase boundary")

       lle_plot_df = df_lle.loc[
           df_lle["converged"].fillna(False).astype(bool) & ~df_lle["llle_detected"].fillna(False).astype(bool)
       ].copy()
       if lle_plot_df.empty:
           raise ValueError("No converged non-LLLE LLE rows are available for this ternary system.")

       feed_points = _collect_feed_points(lle_plot_df)
       phase_i_points = _collect_phase_points(lle_plot_df, "xI")
       phase_ii_points = _collect_phase_points(lle_plot_df, "xII")
       tie_lines = _collect_tie_lines(lle_plot_df)

       tax_lle.scatter(
           feed_points,
           s=44,
           marker="^",
           facecolors="none",
           edgecolors="black",
           linewidths=1.0,
           label="Feed",
       )
       tax_lle.scatter(
           phase_i_points,
           s=34,
           marker="o",
           facecolors="none",
           edgecolors="tab:blue",
           linewidths=1.2,
           label="Phase I boundary",
       )
       tax_lle.scatter(
           phase_ii_points,
           s=34,
           marker="o",
           facecolors="none",
           edgecolors="tab:orange",
           linewidths=1.2,
           label="Phase II boundary",
       )

       for index, (point_i, point_ii) in enumerate(tie_lines):
           tax_lle.line(
               point_i,
               point_ii,
               color="0.65",
               linewidth=0.9,
               alpha=0.85,
               label="Tie line" if index == 0 else None,
           )

       for tax in (tax_stab, tax_lle):
           handles, _ = tax.get_axes().get_legend_handles_labels()
           if handles:
               tax.legend(loc="upper right", frameon=False)
           if hasattr(tax, "_redraw_labels"):
               tax._redraw_labels()

       fig.tight_layout()
       return fig


Python code using pyCRS
=======================

In the first script example, the pyCRS.Database and pyCRS.CRSManager
modules are used to build a small COSKF database and run the
"STABILITY" and "LLE" calculations in batch.

   import os
   from pathlib import Path

   import matplotlib.pyplot as plt
   from pyCRS.CRSManager import CRSSystem
   from pyCRS.Database import COSKFDatabase
   from scm.plams import CRSJob, JobRunner, config
   from TernaryStabilityLLE_common import Compound, build_ternary_grid_df, build_ternary_lle_df, plot_demo_figure

   TEMPERATURE = 298.15
   METHOD = "COSMORS"
   MAX_JOBS = 8
   SOLVENTS = [
       Compound("Water", "Water.coskf"),
       Compound("Ethanol", "Ethanol.coskf"),
       Compound("Benzene", "Benzene.coskf"),
   ]
   DB_NAME = "ternary_stability_lle_demo.db"
   DATABASE_PATH = CRSJob.database()
   IMAGE_NAME = "TernaryStabilityLLE.png"


   def prepare_database(db_name, compounds):
       db_path = Path(db_name)
       if db_path.exists():
           db_path.unlink()

       db = COSKFDatabase(db_name)
       for compound in compounds:
           coskf_path = os.path.join(DATABASE_PATH, compound.coskf)
           if not os.path.exists(coskf_path):
               raise FileNotFoundError(f"COSKF file not found: {coskf_path}")
           db.add_compound(coskf_file=compound.coskf, name=compound.name, coskf_path=DATABASE_PATH)

       return db


   def run_demo(df_stability, df_lle, db):
       cal_stability = CRSSystem()
       for _, row in df_stability.iterrows():
           mixture = {row["s1"]: row["x1"], row["s2"]: row["x2"], row["s3"]: row["x3"]}
           cal_stability.add_Mixture(
               mixture,
               temperature=row["Temperature"],
               problem_type="STABILITY",
               database=db,
               method=row["Method"],
               jobname=f"STABILITY_{row['mix_idx']}",
           )
       cal_stability.runCRSJob()

       for index, out in zip(df_stability.index, cal_stability.outputs):
           results = out.get_results()
           df_stability.loc[index, "unstable"] = bool(results.get("unstable", False))

       cal_lle = CRSSystem()
       for _, row in df_lle.iterrows():
           mixture = {row["s1"]: row["x1"], row["s2"]: row["x2"], row["s3"]: row["x3"]}
           cal_lle.add_Mixture(
               mixture,
               temperature=row["Temperature"],
               problem_type="LLE",
               database=db,
               method=row["Method"],
               jobname=f"LLE_{row['mix_idx']}",
           )
       cal_lle.runCRSJob()

       for index, out in zip(df_lle.index, cal_lle.outputs):
           results = out.get_results()
           converged = bool(results.get("converged", False))
           llle_detected = bool(results.get("llle_detected", False))
           df_lle.loc[index, "converged"] = converged
           df_lle.loc[index, "llle_detected"] = llle_detected
           if converged:
               x_i = results.get("xI")
               x_ii = results.get("xII")
               for component_index in range(3):
                   component_number = component_index + 1
                   if x_i is not None:
                       df_lle.loc[index, f"xI{component_number}"] = x_i[component_index][0]
                   if x_ii is not None:
                       df_lle.loc[index, f"xII{component_number}"] = x_ii[component_index][0]

       return df_stability, df_lle


   db = prepare_database(DB_NAME, SOLVENTS)

   config.default_jobrunner = JobRunner(parallel=True, maxjobs=MAX_JOBS)
   config.job.runscript.nproc = 1
   config.log.stdout = 1

   df_stability = build_ternary_grid_df(SOLVENTS, min_frac=0.1, step=0.1, temperature=TEMPERATURE, method=METHOD)
   df_lle = build_ternary_lle_df(
       SOLVENTS,
       s1_ratio=1.0,
       s2_ratio=1.0,
       s3_min=0.1,
       s3_max=0.9,
       num_points=9,
       temperature=TEMPERATURE,
       method=METHOD,
   )

   df_stability, df_lle = run_demo(df_stability, df_lle, db)
   fig = plot_demo_figure(df_stability, df_lle)

   print(df_stability[["mix_idx", "x1", "x2", "x3", "unstable"]])
   print(df_lle[["mix_idx", "x1", "x2", "x3", "converged", "llle_detected", "xI1", "xI2", "xI3", "xII1", "xII2", "xII3"]])
   fig.savefig(IMAGE_NAME, dpi=200, bbox_inches="tight")
   plt.show()

This code produces the following output:

   [image]Ternary stability grid and LLE phase-boundary plot for the
   water/ethanol/benzene system at 298.15 K.


Python code using PLAMS
=======================

In the second script example, the same workflow is set up directly
with PLAMS input setting blocks for each "CRSJob". This is useful when
the calculation setup needs to stay close to the underlying COSMO-RS
input.

   import os

   import matplotlib.pyplot as plt
   from scm.plams import CRSJob, JobRunner, Settings, config
   from TernaryStabilityLLE_common import Compound, build_ternary_grid_df, build_ternary_lle_df, plot_demo_figure

   TEMPERATURE = 298.15
   METHOD = "COSMORS"
   MAX_JOBS = 8
   SOLVENTS = [
       Compound("Water", "Water.coskf"),
       Compound("Ethanol", "Ethanol.coskf"),
       Compound("Benzene", "Benzene.coskf"),
   ]
   DATABASE_PATH = CRSJob.database()
   IMAGE_NAME = "TernaryStabilityLLE.png"


   def ensure_compound_files(compounds):
       for compound in compounds:
           coskf_path = os.path.join(DATABASE_PATH, compound.coskf)
           if not os.path.exists(coskf_path):
               raise FileNotFoundError(f"COSKF file not found: {coskf_path}")


   def run_demo(df_stability, df_lle):
       outputs = []
       for _, row in df_stability.iterrows():
           settings = Settings()
           settings.input.compound = []
           for component_index, compound in enumerate(SOLVENTS, start=1):
               block = Settings()
               block._h = os.path.join(DATABASE_PATH, compound.coskf)
               block.frac1 = float(row[f"x{component_index}"])
               block.name = compound.name.lower()
               settings.input.compound.append(block)

           settings.input.method = row["Method"]
           settings.input.property._h = "STABILITY"
           settings.input.temperature = row["Temperature"]

           job = CRSJob(settings=settings, name=f"STABILITY_{row['mix_idx']}")
           outputs.append(job.run())

       for index, out in zip(df_stability.index, outputs):
           results = out.get_results()
           df_stability.loc[index, "unstable"] = bool(results.get("unstable", False))

       outputs = []
       for _, row in df_lle.iterrows():
           settings = Settings()
           settings.input.compound = []
           for component_index, compound in enumerate(SOLVENTS, start=1):
               block = Settings()
               block._h = os.path.join(DATABASE_PATH, compound.coskf)
               block.frac1 = float(row[f"x{component_index}"])
               block.name = compound.name.lower()
               settings.input.compound.append(block)

           settings.input.method = row["Method"]
           settings.input.property._h = "LLE"
           settings.input.temperature = row["Temperature"]

           job = CRSJob(settings=settings, name=f"LLE_{row['mix_idx']}")
           outputs.append(job.run())

       for index, out in zip(df_lle.index, outputs):
           results = out.get_results()
           converged = bool(results.get("converged", False))
           llle_detected = bool(results.get("llle_detected", False))
           df_lle.loc[index, "converged"] = converged
           df_lle.loc[index, "llle_detected"] = llle_detected
           if converged:
               x_i = results.get("xI")
               x_ii = results.get("xII")
               for component_index in range(3):
                   component_number = component_index + 1
                   if x_i is not None:
                       df_lle.loc[index, f"xI{component_number}"] = x_i[component_index][0]
                   if x_ii is not None:
                       df_lle.loc[index, f"xII{component_number}"] = x_ii[component_index][0]

       return df_stability, df_lle


   ensure_compound_files(SOLVENTS)

   config.default_jobrunner = JobRunner(parallel=True, maxjobs=MAX_JOBS)
   config.job.runscript.nproc = 1
   config.log.stdout = 1

   df_stability = build_ternary_grid_df(SOLVENTS, min_frac=0.1, step=0.1, temperature=TEMPERATURE, method=METHOD)
   df_lle = build_ternary_lle_df(
       SOLVENTS,
       s1_ratio=1.0,
       s2_ratio=1.0,
       s3_min=0.1,
       s3_max=0.9,
       num_points=9,
       temperature=TEMPERATURE,
       method=METHOD,
   )

   df_stability, df_lle = run_demo(df_stability, df_lle)
   fig = plot_demo_figure(df_stability, df_lle)

   print(df_stability[["mix_idx", "x1", "x2", "x3", "unstable"]])
   print(df_lle[["mix_idx", "x1", "x2", "x3", "converged", "llle_detected", "xI1", "xI2", "xI3", "xII1", "xII2", "xII3"]])
   fig.savefig(IMAGE_NAME, dpi=200, bbox_inches="tight")
   plt.show()

This script generates the same ternary figure as shown above.


Download scripts
================

"Download TernaryStabilityLLE_common.py"

"Download TernaryStabilityLLE_pyCRS.py"

"Download TernaryStabilityLLE_PLAMS.py"
