From ee76d72b83706210e3d85cff02e5749f45eeccdf Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 00:29:29 -0600
Subject: [PATCH 1/7] add EOS-mixture options

---
 thermosteam/_thermo.py         |   4 +-
 thermosteam/mixture/mixture.py | 517 ++++++++++++++++++++++-----------
 2 files changed, 356 insertions(+), 165 deletions(-)

diff --git a/thermosteam/_thermo.py b/thermosteam/_thermo.py
index d7c6786a..d19bd38f 100644
--- a/thermosteam/_thermo.py
+++ b/thermosteam/_thermo.py
@@ -12,7 +12,7 @@
 from . import equilibrium as eq
 from ._chemical import Chemical
 from ._chemicals import Chemicals
-from .mixture import Mixture
+from .mixture import Mixture, IdealMixture
 from .utils import read_only, cucumber
 
 __all__ = ('Thermo', 'IdealThermo')
@@ -138,7 +138,7 @@ def __init__(self, chemicals, mixture=None,
         if PCF is None: PCF = eq.MockPoyintingCorrectionFactors
         if not isinstance(chemicals, Chemicals): chemicals = Chemicals(chemicals, cache)
         if not mixture:
-            mixture = Mixture.from_chemicals(chemicals)
+            mixture = IdealMixture.from_chemicals(chemicals)
         elif not isinstance(mixture, Mixture): # pragma: no cover
             raise ValueError(f"mixture must be a '{Mixture.__name__}' object")
         chemicals.compile(skip_checks=skip_checks)
diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index bb09028d..f1263ee4 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -8,11 +8,14 @@
 """
 """
 import flexsolve as flx
+import thermosteam as tmo
 import numpy as np
 from math import exp
 from thermosteam import functional as fn
+from thermo.interaction_parameters import IPDB
+from thermo import eos_mix
 from .. import units_of_measure as thermo_units
-from ..base import PhaseHandle, MockPhaseTHandle, MockPhaseTPHandle, sparse
+from ..base import PhaseHandle, MockPhaseTHandle, MockPhaseTPHandle, SparseVector, sparse
 from .ideal_mixture_model import (
     SinglePhaseIdealTMixtureModel,
     IdealTMixtureModel, 
@@ -22,7 +25,7 @@
 )
 from .._chemicals import Chemical, CompiledChemicals, chemical_data_array
 
-__all__ = ('Mixture',)
+__all__ = ('Mixture', 'IdealMixture')
 
 # %% Functions for building mixture models
 
@@ -73,177 +76,45 @@ def xiter_T_at_SP(T, S, S_model, phase_mol, P, Cn_model, Cn_cache):
     return T * exp((S - S_model(phase_mol, T, P)) / Cn)
 
 
-# %% Ideal mixture
+# %% Abstract mixture
 
 class Mixture:
     """
-    Create an Mixture object for estimating mixture properties.
-    
-    Parameters
-    ----------
-    rule : str
-        Description of mixing rules used.
-    Cn : function(phase, mol, T)
-        Molar isobaric heat capacity mixture model [J/mol/K].
-    H : function(phase, mol, T)
-        Enthalpy mixture model [J/mol].
-    S : function(phase, mol, T, P)
-        Entropy mixture model [J/mol].
-    H_excess : function(phase, mol, T, P)
-        Excess enthalpy mixture model [J/mol].
-    S_excess : function(phase, mol, T, P)
-        Excess entropy mixture model [J/mol].
-    mu : function(phase, mol, T, P)
-        Dynamic viscosity mixture model [Pa*s].
-    V : function(phase, mol, T, P)
-        Molar volume mixture model [m^3/mol].
-    kappa : function(phase, mol, T, P)
-        Thermal conductivity mixture model [W/m/K].
-    Hvap : function(mol, T)
-        Heat of vaporization mixture model [J/mol]
-    sigma : function(mol, T, P)
-        Surface tension mixture model [N/m].
-    epsilon : function(mol, T, P)
-        Relative permitivity mixture model [-]
+    Abstract class for estimating mixture properties.
+    
+    Abstract methods
+    ----------------
+    Cn(phase, mol, T) : float
+        Molar isobaric heat capacity [J/mol/K].
+    H(phase, mol, T) : float
+        Enthalpy [J/mol].
+    S(phase, mol, T, P) : float
+        Entropy [J/mol].
+    mu(phase, mol, T, P) : float
+        Dynamic viscosity [Pa*s].
+    V(phase, mol, T, P) : float
+        Molar volume [m^3/mol].
+    kappa(phase, mol, T, P) : float
+        Thermal conductivity [W/m/K].
+    Hvap(mol, T) : float
+        Heat of vaporization [J/mol]
+    sigma(mol, T, P) : float
+        Surface tension [N/m].
+    epsilon(mol, T, P) : float
+        Relative permitivity [-]
+    
+    Abstract attributes
     MWs : 1d array[float]
         Component molecular weights [g/mol].
-    include_excess_energies=False : bool
-        Whether to include excess energies
-        in enthalpy and entropy calculations.
-    
-    Notes
-    -----
-    Although the mixture models are on a molar basis, this is only if the molar
-    data is normalized before the calculation (i.e. the `mol` parameter is 
-    normalized before being passed to the model).
-    
-    See also
-    --------
-    IdealTMixtureModel
-    IdealTPMixtureModel
-    
-    Attributes
-    ----------
-    rule : str
-        Description of mixing rules used.
-    include_excess_energies : bool
-        Whether to include excess energies
-        in enthalpy and entropy calculations.
-    Cn(phase, mol, T) : 
-        Mixture molar isobaric heat capacity [J/mol/K].
-    mu(phase, mol, T, P) : 
-        Mixture dynamic viscosity [Pa*s].
-    V(phase, mol, T, P) : 
-        Mixture molar volume [m^3/mol].
-    kappa(phase, mol, T, P) : 
-        Mixture thermal conductivity [W/m/K].
-    Hvap(mol, T, P) : 
-        Mixture heat of vaporization [J/mol]
-    sigma(mol, T, P) : 
-        Mixture surface tension [N/m].
-    epsilon(mol, T, P) : 
-        Mixture relative permitivity [-].
-    MWs : 1d-array[float]
-        Component molecular weights [g/mol].
     
     """
     maxiter = 20
     T_tol = 1e-6
-    __slots__ = ('rule',
-                 'rigorous_energy_balance',
-                 'include_excess_energies',
-                 'Cn', 'mu', 'V', 'kappa',
-                 'Hvap', 'sigma', 'epsilon',
-                 'MWs', '_H', '_H_excess', '_S', '_S_excess',
-    )
-    
-    def __init__(self, rule, Cn, H, S, H_excess, S_excess,
-                 mu, V, kappa, Hvap, sigma, epsilon,
-                 MWs, include_excess_energies=False):
-        self.rule = rule
-        self.include_excess_energies = include_excess_energies
-        self.Cn = Cn
-        self.mu = mu
-        self.V = V
-        self.kappa = kappa
-        self.Hvap = Hvap
-        self.sigma = sigma
-        self.epsilon = epsilon
-        self.MWs = MWs
-        self._H = H
-        self._S = S
-        self._H_excess = H_excess
-        self._S_excess = S_excess
-    
-    @classmethod
-    def from_chemicals(cls, chemicals, 
-                       include_excess_energies=False, 
-                       rule='ideal',
-                       cache=True):
-        """
-        Create a Mixture object from chemical objects.
-        
-        Parameters
-        ----------
-        chemicals : Iterable[Chemical]
-            For retrieving pure component chemical data.
-        include_excess_energies=False : bool
-            Whether to include excess energies in enthalpy and entropy calculations.
-        rule : str, optional
-            Mixing rule. Defaults to 'ideal'.
-        cache : optional
-            Whether or not to use cached chemicals and cache new chemicals. Defaults to True.
-    
-        See also
-        --------
-        :class:`~.mixture.Mixture`
-        :class:`~.IdealMixtureModel`
-    
-        Examples
-        --------
-        Calculate enthalpy of evaporation for a water and ethanol mixture:
-        
-        >>> from thermosteam import Mixture
-        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
-        >>> mixture.Hvap([0.2, 0.8], 350)
-        39750.62
-
-        Calculate density for a water and ethanol mixture in g/L:
-
-        >>> from thermosteam import Mixture
-        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
-        >>> mixture.get_property('rho', 'g/L', 'l', [0.2, 0.8], 350, 101325)
-        754.23
-        
-        """
-        if rule == 'ideal':
-            isa = isinstance
-            if isa(chemicals, CompiledChemicals):
-                MWs = chemicals.MW
-                chemicals = chemicals.tuple
-            else:
-                chemicals = [(i if isa(i, Chemical) else Chemical(i)) for i in chemicals]
-                MWs = chemical_data_array(chemicals, 'MW')
-            getfield = getattr
-            Cn =  create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
-            H =  create_mixture_model(chemicals, 'H', IdealTPMixtureModel)
-            S = create_mixture_model(chemicals, 'S', IdealEntropyModel)
-            H_excess = create_mixture_model(chemicals, 'H_excess', IdealTPMixtureModel)
-            S_excess = create_mixture_model(chemicals, 'S_excess', IdealTPMixtureModel)
-            mu = create_mixture_model(chemicals, 'mu', IdealTPMixtureModel)
-            V = create_mixture_model(chemicals, 'V', IdealTPMixtureModel)
-            kappa = create_mixture_model(chemicals, 'kappa', IdealTPMixtureModel)
-            Hvap = IdealHvapModel(chemicals)
-            sigma = SinglePhaseIdealTMixtureModel([getfield(i, 'sigma') for i in chemicals], 'sigma')
-            epsilon = SinglePhaseIdealTMixtureModel([getfield(i, 'epsilon') for i in chemicals], 'epsilon')
-            return cls(rule, Cn, H, S, H_excess, S_excess,
-                       mu, V, kappa, Hvap, sigma, epsilon, MWs, include_excess_energies)
-        else:
-            raise ValueError("rule '{rule}' is not available (yet)")
+    __slots__ = ()
     
     def MW(self, mol):
         """Return molecular weight [g/mol] given molar array [mol]."""
-        total_mol = sparse(mol).sum()
+        total_mol = mol.sum()
         return (mol * self.MWs).sum() / total_mol if total_mol else 0.
     
     def rho(self, phase, mol, T, P):
@@ -326,8 +197,8 @@ def H(self, phase, mol, T, P):
     
     def S(self, phase, mol, T, P):
         """Return entropy in [J/mol/K]."""
-        total_mol = mol.sum()
-        if total_mol == 0.: return 0.
+        if mol.__class__ is not SparseVector: mol = SparseVector(mol)
+        if not mol.dct: return 0.
         S = self._S(phase, mol, T, P)
         if self.include_excess_energies:
             S += self._S_excess(phase, mol, T, P)
@@ -435,4 +306,324 @@ def _info(self):
     def show(self):
         print(self._info())
     _ipython_display_ = show
+
+
+# %% Ideal mixture
+
+class IdealMixture(Mixture):
+    """
+    Create an Mixture object for estimating mixture properties.
+    
+    Parameters
+    ----------
+    Cn : function(phase, mol, T)
+        Molar isobaric heat capacity mixture model [J/mol/K].
+    H : function(phase, mol, T)
+        Enthalpy mixture model [J/mol].
+    S : function(phase, mol, T, P)
+        Entropy mixture model [J/mol].
+    H_excess : function(phase, mol, T, P)
+        Excess enthalpy mixture model [J/mol].
+    S_excess : function(phase, mol, T, P)
+        Excess entropy mixture model [J/mol].
+    mu : function(phase, mol, T, P)
+        Dynamic viscosity mixture model [Pa*s].
+    V : function(phase, mol, T, P)
+        Molar volume mixture model [m^3/mol].
+    kappa : function(phase, mol, T, P)
+        Thermal conductivity mixture model [W/m/K].
+    Hvap : function(mol, T)
+        Heat of vaporization mixture model [J/mol]
+    sigma : function(mol, T, P)
+        Surface tension mixture model [N/m].
+    epsilon : function(mol, T, P)
+        Relative permitivity mixture model [-]
+    MWs : 1d array[float]
+        Component molecular weights [g/mol].
+    include_excess_energies=False : bool
+        Whether to include excess energies
+        in enthalpy and entropy calculations.
+    
+    Notes
+    -----
+    Although the mixture models are on a molar basis, this is only if the molar
+    data is normalized before the calculation (i.e. the `mol` parameter is 
+    normalized before being passed to the model).
+    
+    See also
+    --------
+    IdealTMixtureModel
+    IdealTPMixtureModel
+    
+    Attributes
+    ----------
+    include_excess_energies : bool
+        Whether to include excess energies
+        in enthalpy and entropy calculations.
+    Cn(phase, mol, T) : 
+        Mixture molar isobaric heat capacity [J/mol/K].
+    mu(phase, mol, T, P) : 
+        Mixture dynamic viscosity [Pa*s].
+    V(phase, mol, T, P) : 
+        Mixture molar volume [m^3/mol].
+    kappa(phase, mol, T, P) : 
+        Mixture thermal conductivity [W/m/K].
+    Hvap(mol, T, P) : 
+        Mixture heat of vaporization [J/mol]
+    sigma(mol, T, P) : 
+        Mixture surface tension [N/m].
+    epsilon(mol, T, P) : 
+        Mixture relative permitivity [-].
+    MWs : 1d-array[float]
+        Component molecular weights [g/mol].
+    
+    """
+    __slots__ = (
+        'include_excess_energies',
+        'Cn', 'mu', 'V', 'kappa',
+        'Hvap', 'sigma', 'epsilon',
+        'MWs', '_H', '_H_excess', '_S', '_S_excess',
+    )
+    
+    def __init__(self, Cn, H, S, H_excess, S_excess,
+                 mu, V, kappa, Hvap, sigma, epsilon,
+                 MWs, include_excess_energies=False):
+        self.include_excess_energies = include_excess_energies
+        self.Cn = Cn
+        self.mu = mu
+        self.V = V
+        self.kappa = kappa
+        self.Hvap = Hvap
+        self.sigma = sigma
+        self.epsilon = epsilon
+        self.MWs = MWs
+        self._H = H
+        self._S = S
+        self._H_excess = H_excess
+        self._S_excess = S_excess
+    
+    @classmethod
+    def from_chemicals(cls, chemicals, 
+                       include_excess_energies=False, 
+                       cache=True):
+        """
+        Create a Mixture object from chemical objects.
+        
+        Parameters
+        ----------
+        chemicals : Iterable[Chemical]
+            For retrieving pure component chemical data.
+        include_excess_energies=False : bool
+            Whether to include excess energies in enthalpy and entropy calculations.
+        rule : str, optional
+            Mixing rule. Defaults to 'ideal'.
+        cache : optional
+            Whether or not to use cached chemicals and cache new chemicals. Defaults to True.
+    
+        See also
+        --------
+        :class:`~.mixture.Mixture`
+        :class:`~.IdealMixtureModel`
+    
+        Examples
+        --------
+        Calculate enthalpy of evaporation for a water and ethanol mixture:
+        
+        >>> from thermosteam import Mixture
+        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> mixture.Hvap([0.2, 0.8], 350)
+        39750.62
+
+        Calculate density for a water and ethanol mixture in g/L:
+
+        >>> from thermosteam import Mixture
+        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> mixture.get_property('rho', 'g/L', 'l', [0.2, 0.8], 350, 101325)
+        754.23
+        
+        """
+        isa = isinstance
+        if isa(chemicals, CompiledChemicals):
+            MWs = chemicals.MW
+            chemicals = chemicals.tuple
+        else:
+            chemicals = [(i if isa(i, Chemical) else Chemical(i)) for i in chemicals]
+            MWs = chemical_data_array(chemicals, 'MW')
+        getfield = getattr
+        Cn =  create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
+        H =  create_mixture_model(chemicals, 'H', IdealTPMixtureModel)
+        S = create_mixture_model(chemicals, 'S', IdealEntropyModel)
+        H_excess = create_mixture_model(chemicals, 'H_excess', IdealTPMixtureModel)
+        S_excess = create_mixture_model(chemicals, 'S_excess', IdealTPMixtureModel)
+        mu = create_mixture_model(chemicals, 'mu', IdealTPMixtureModel)
+        V = create_mixture_model(chemicals, 'V', IdealTPMixtureModel)
+        kappa = create_mixture_model(chemicals, 'kappa', IdealTPMixtureModel)
+        Hvap = IdealHvapModel(chemicals)
+        sigma = SinglePhaseIdealTMixtureModel([getfield(i, 'sigma') for i in chemicals], 'sigma')
+        epsilon = SinglePhaseIdealTMixtureModel([getfield(i, 'epsilon') for i in chemicals], 'epsilon')
+        return cls(Cn, H, S, H_excess, S_excess,
+                   mu, V, kappa, Hvap, sigma, epsilon, MWs, include_excess_energies)
+    
+    def __repr__(self):
+        return f"{type(self).__name__}(rule={repr(self.rule)}, ..., include_excess_energies={self.include_excess_energies})"
+    
+    def _info(self):
+        return (f"{type(self).__name__}(\n"
+                f"    include_excess_energies={self.include_excess_energies}\n"
+                 ")")
+    
+    def show(self):
+        print(self._info())
+    _ipython_display_ = show
+    
+# %% Thermo mixture
+
+class EOSMixture:
+    __slots__ = (
+        'chemicals', 'Cn', 'mu', 'V', 'kappa',
+        'Hvap', 'sigma', 'epsilon',
+        'MWs', 'chemicals', 
+    )
+    chemsep_db = None
+    
+    def __init__(self, chemicals, Cn, mu, V, kappa, Hvap, sigma, epsilon, MWs):
+        self.chemicals = chemicals
+        self.Cn = Cn
+        self.mu = mu
+        self.V = V
+        self.kappa = kappa
+        self.Hvap = Hvap
+        self.sigma = sigma
+        self.epsilon = epsilon
+        self.MWs = MWs
+    
+    def eos(self, phase, chemicals, zs, T, P):
+        key = (phase, chemicals)
+        cache = self.cache
+        if key in cache:
+            return cache[key]
+        else:
+            only_g = phase == 'g'
+            only_l = phase == 'l'
+            if only_g or only_l:
+                try:
+                    self._eos = eos = self._eos.to_TP_zs_fast(
+                        T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
+                        full_alphas=False
+                    )
+                except:
+                    data = tmo.ChemicalData(chemicals)
+                    if self.chemsep_db is None:
+                        kijs = None
+                    else:
+                        try:
+                            kijs = IPDB.get_ip_asymmetric_matrix(self.chemsep_db, data.CASs, 'kij')
+                        except:
+                            kijs = None
+                    self._eos = eos = self.EOS(
+                        zs=zs, T=T, P=P, Tcs=data.Tcs, Pcs=data.Pcs, omegas=data.omegas, kijs=kijs,
+                        only_g=only_g, only_l=only_l,
+                    )    
+        return eos
+    
+    def H(self, phase, mol, T, P):
+        """Return enthalpy [J/mol]."""
+        if mol.__class__ is not SparseVector: mol = SparseVector(mol)
+        chemicals = self.chemicals
+        H = sum([j * chemicals[i].H(phase, T, P) for i, j in mol.dct.items()])
+        if phase != 's':
+            subset_chemicals = tuple([chemicals[i] for i in mol.nonzero_keys()])
+            mol = [i for i in mol.nonzero_values()]
+            eos = self.eos(phase, subset_chemicals, mol, T, P)
+            if phase == 'l':
+                H += eos.H_dep_l
+            else:
+                H += eos.H_dep_g
+        return H
+    
+    def S(self, phase, mol, T, P):
+        """Return entropy [J/mol/K]."""
+        if mol.__class__ is not SparseVector: mol = SparseVector(mol)
+        chemicals = self.chemicals
+        S = sum([j * chemicals[i].S(phase, T, P) for i, j in mol.dct.items()])
+        if phase != 's':
+            subset_chemicals = tuple([chemicals[i] for i in mol.nonzero_keys()])
+            mol = [i for i in mol.nonzero_values()]
+            eos = self.eos(phase, subset_chemicals, mol, T, P)
+            if phase == 'l':
+                S += eos.S_dep_l
+            else:
+                S += eos.S_dep_g
+        return S
+    
+    @classmethod
+    def from_chemicals(cls, chemicals, cache=True):
+        """
+        Create a Mixture object from chemical objects.
         
+        Parameters
+        ----------
+        chemicals : Iterable[Chemical]
+            For retrieving pure component chemical data.
+        cache : optional
+            Whether or not to use cached chemicals and cache new chemicals. Defaults to True.
+    
+        See also
+        --------
+        :class:`~.mixture.Mixture`
+        :class:`~.IdealMixtureModel`
+    
+        Examples
+        --------
+        Calculate enthalpy of evaporation for a water and ethanol mixture:
+        
+        >>> from thermosteam import Mixture
+        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> mixture.Hvap([0.2, 0.8], 350)
+        39750.62
+
+        Calculate density for a water and ethanol mixture in g/L:
+
+        >>> from thermosteam import Mixture
+        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> mixture.get_property('rho', 'g/L', 'l', [0.2, 0.8], 350, 101325)
+        754.23
+        
+        """
+        isa = isinstance
+        if isa(chemicals, CompiledChemicals):
+            MWs = chemicals.MW
+            chemicals = chemicals.tuple
+        else:
+            chemicals = [(i if isa(i, Chemical) else Chemical(i)) for i in chemicals]
+            MWs = chemical_data_array(chemicals, 'MW')
+        getfield = getattr
+        Cn =  create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
+        mu = create_mixture_model(chemicals, 'mu', IdealTPMixtureModel)
+        V = create_mixture_model(chemicals, 'V', IdealTPMixtureModel)
+        kappa = create_mixture_model(chemicals, 'kappa', IdealTPMixtureModel)
+        Hvap = IdealHvapModel(chemicals)
+        sigma = SinglePhaseIdealTMixtureModel([getfield(i, 'sigma') for i in chemicals], 'sigma')
+        epsilon = SinglePhaseIdealTMixtureModel([getfield(i, 'epsilon') for i in chemicals], 'epsilon')
+        return cls(chemicals, Cn, mu, V, kappa, Hvap, sigma, epsilon, MWs)
+    
+    @classmethod
+    def subclass(cls, EOS, name=None):
+        if name is None: name = EOS.__name__.replace('MIX', '') + 'Mixture'
+        return type(name, (cls,), dict(EOS=EOS, cache={}))
+
+dct = globals()    
+clsnames = []
+for name in ('PRMIX', 'SRKMIX', 'PR78MIX', 'VDWMIX', 'PRSVMIX',
+             'PRSV2MIX', 'TWUPRMIX', 'TWUSRKMIX', 'APISRKMIX', 'IGMIX', 'RKMIX',
+             'PRMIXTranslatedConsistent', 'PRMIXTranslatedPPJP', 'PRMIXTranslated',
+             'SRKMIXTranslatedConsistent', 'PSRK', 'MSRKMIXTranslated',
+             'SRKMIXTranslated'):
+    cls = EOSMixture.subclass(getattr(eos_mix, name))
+    clsname = cls.__name__
+    clsnames.append(clsname)
+    dct[clsname] = cls
+
+dct['PRMixture'].chemsep_db = 'ChemSep PR'
+__all__ = (*__all__, *clsnames)
+del dct, clsnames
\ No newline at end of file

From c59d3e2b0e9ce813e3857221d0ff7035d010984f Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 08:20:55 -0600
Subject: [PATCH 2/7] get eos-mixture working

---
 thermosteam/mixture/mixture.py | 126 ++++++++++++++++++++++++---------
 1 file changed, 92 insertions(+), 34 deletions(-)

diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index f1263ee4..1ada29bc 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -50,28 +50,28 @@ def create_mixture_model(chemicals, var, Model):
 def iter_T_at_HP(T, H, H_model, phase, mol, P, Cn_model, Cn_cache):
     # Used to solve for temperature at given ethalpy 
     counter, Cn = Cn_cache
-    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase, mol, T)
+    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase, mol, T, P)
     Cn_cache[0] += 1
     return T + (H - H_model(phase, mol, T, P)) / Cn
 
 def xiter_T_at_HP(T, H, H_model, phase_mol, P, Cn_model, Cn_cache):
     # Used to solve for temperature at given ethalpy 
     counter, Cn = Cn_cache
-    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase_mol, T)
+    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase_mol, T, P)
     Cn_cache[0] += 1
     return T + (H - H_model(phase_mol, T, P)) / Cn
 
 def iter_T_at_SP(T, S, S_model, phase, mol, P, Cn_model, Cn_cache):
     # Used to solve for temperature at given entropy 
     counter, Cn = Cn_cache
-    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase, mol, T)
+    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase, mol, T, P)
     Cn_cache[0] += 1
     return T * exp((S - S_model(phase, mol, T, P)) / Cn)
 
 def xiter_T_at_SP(T, S, S_model, phase_mol, P, Cn_model, Cn_cache):
     # Used to solve for temperature at given entropy 
     counter, Cn = Cn_cache
-    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase_mol, T)
+    if not counter % 5: Cn_cache[1] = Cn = Cn_model(phase_mol, T, P)
     Cn_cache[0] += 1
     return T * exp((S - S_model(phase_mol, T, P)) / Cn)
 
@@ -295,7 +295,7 @@ def xkappa(self, phase_mol, T, P):
         return sum([self.kappa(phase, mol, T, P) for phase, mol in phase_mol])
     
     def __repr__(self):
-        return f"{type(self).__name__}(rule={repr(self.rule)}, ..., include_excess_energies={self.include_excess_energies})"
+        return f"{type(self).__name__}(..., include_excess_energies={self.include_excess_energies})"
     
     def _info(self):
         return (f"{type(self).__name__}(\n"
@@ -478,17 +478,19 @@ def show(self):
     
 # %% Thermo mixture
 
-class EOSMixture:
+class EOSMixture(Mixture):
     __slots__ = (
-        'chemicals', 'Cn', 'mu', 'V', 'kappa',
-        'Hvap', 'sigma', 'epsilon',
-        'MWs', 'chemicals', 
+        'chemicals', 'eos_chemicals', 'Cn_ideal', 'mu', 'V', 'kappa',
+        'Hvap', 'sigma', 'epsilon', 'H_ideal', 'S_ideal', 'MWs', 
     )
     chemsep_db = None
     
-    def __init__(self, chemicals, Cn, mu, V, kappa, Hvap, sigma, epsilon, MWs):
+    def __init__(self, chemicals, eos_chemicals, Cn_ideal, H_ideal, S_ideal, mu, V, kappa, Hvap, sigma, epsilon, MWs):
         self.chemicals = chemicals
-        self.Cn = Cn
+        self.eos_chemicals = eos_chemicals
+        self.Cn_ideal = Cn_ideal
+        self.H_ideal = H_ideal
+        self.S_ideal = S_ideal
         self.mu = mu
         self.V = V
         self.kappa = kappa
@@ -507,9 +509,9 @@ def eos(self, phase, chemicals, zs, T, P):
             only_l = phase == 'l'
             if only_g or only_l:
                 try:
-                    self._eos = eos = self._eos.to_TP_zs_fast(
+                    self._eos = eos = self._eos.to_TP_zs(
                         T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
-                        full_alphas=False
+                        fugacities=False
                     )
                 except:
                     data = tmo.ChemicalData(chemicals)
@@ -522,42 +524,92 @@ def eos(self, phase, chemicals, zs, T, P):
                             kijs = None
                     self._eos = eos = self.EOS(
                         zs=zs, T=T, P=P, Tcs=data.Tcs, Pcs=data.Pcs, omegas=data.omegas, kijs=kijs,
-                        only_g=only_g, only_l=only_l,
+                        only_g=only_g, only_l=only_l, fugacities=False, 
                     )    
         return eos
     
+    def Cn(self, phase, mol, T, P):
+        if mol.__class__ is not SparseVector: mol = SparseVector(mol)
+        chemicals = self.chemicals
+        dct = mol.dct
+        Cn = self.Cn_ideal(phase, mol, T, P)
+        if phase != 's':
+            eos_chemicals = self.eos_chemicals
+            chemical_subset = []
+            mol_subset = []
+            eos_mol = 0
+            for i, j in dct.items():
+                chemical = chemicals[i]
+                if chemical in eos_chemicals:
+                    chemical_subset.append(chemical)
+                    mol_subset.append(j)
+                    eos_mol += j
+            zs = [i / eos_mol for i in mol_subset]
+            eos = self.eos(
+                phase, tuple(chemical_subset), zs, T, P
+            )
+            if phase == 'l':
+                Cn += eos.Cp_dep_l * eos_mol
+            else:
+                Cn += eos.Cp_dep_g * eos_mol
+        return Cn
+    
     def H(self, phase, mol, T, P):
         """Return enthalpy [J/mol]."""
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
         chemicals = self.chemicals
-        H = sum([j * chemicals[i].H(phase, T, P) for i, j in mol.dct.items()])
+        dct = mol.dct
+        H = self.H_ideal(phase, mol, T, P)
         if phase != 's':
-            subset_chemicals = tuple([chemicals[i] for i in mol.nonzero_keys()])
-            mol = [i for i in mol.nonzero_values()]
-            eos = self.eos(phase, subset_chemicals, mol, T, P)
+            eos_chemicals = self.eos_chemicals
+            chemical_subset = []
+            mol_subset = []
+            eos_mol = 0
+            for i, j in dct.items():
+                chemical = chemicals[i]
+                if chemical in eos_chemicals:
+                    chemical_subset.append(chemical)
+                    mol_subset.append(j)
+                    eos_mol += j
+            zs = [i / eos_mol for i in mol_subset]
+            eos = self.eos(
+                phase, tuple(chemical_subset), zs, T, P
+            )
             if phase == 'l':
-                H += eos.H_dep_l
+                H += eos.H_dep_l * eos_mol
             else:
-                H += eos.H_dep_g
+                H += eos.H_dep_g * eos_mol
         return H
     
     def S(self, phase, mol, T, P):
         """Return entropy [J/mol/K]."""
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
         chemicals = self.chemicals
-        S = sum([j * chemicals[i].S(phase, T, P) for i, j in mol.dct.items()])
+        dct = mol.dct
+        S = self.S_ideal(phase, mol, T, P)
         if phase != 's':
-            subset_chemicals = tuple([chemicals[i] for i in mol.nonzero_keys()])
-            mol = [i for i in mol.nonzero_values()]
-            eos = self.eos(phase, subset_chemicals, mol, T, P)
+            eos_chemicals = self.eos_chemicals
+            chemical_subset = []
+            mol_subset = []
+            eos_mol = 0
+            for i, j in dct.items():
+                chemical = chemicals[i]
+                if chemical in eos_chemicals:
+                    chemical_subset.append(chemical)
+                    mol_subset.append(j)
+                    eos_mol += j
+            zs = [i / eos_mol for i in mol_subset]
+            eos = self.eos(
+                phase, tuple(chemical_subset), zs, T, P
+            )
             if phase == 'l':
-                S += eos.S_dep_l
+                S += eos.S_dep_l * eos_mol
             else:
-                S += eos.S_dep_g
+                S += eos.S_dep_g * eos_mol
         return S
     
     @classmethod
-    def from_chemicals(cls, chemicals, cache=True):
+    def from_chemicals(cls, chemicals, eos_chemicals, cache=True):
         """
         Create a Mixture object from chemical objects.
         
@@ -565,14 +617,11 @@ def from_chemicals(cls, chemicals, cache=True):
         ----------
         chemicals : Iterable[Chemical]
             For retrieving pure component chemical data.
+        eos_chemicals : Iterable[Chemical]
+            Chemicals with equations of state.
         cache : optional
             Whether or not to use cached chemicals and cache new chemicals. Defaults to True.
     
-        See also
-        --------
-        :class:`~.mixture.Mixture`
-        :class:`~.IdealMixtureModel`
-    
         Examples
         --------
         Calculate enthalpy of evaporation for a water and ethanol mixture:
@@ -598,20 +647,29 @@ def from_chemicals(cls, chemicals, cache=True):
             chemicals = [(i if isa(i, Chemical) else Chemical(i)) for i in chemicals]
             MWs = chemical_data_array(chemicals, 'MW')
         getfield = getattr
-        Cn =  create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
+        Cn = create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
+        H = create_mixture_model(chemicals, 'H', IdealTPMixtureModel)
+        S = create_mixture_model(chemicals, 'S', IdealEntropyModel)
         mu = create_mixture_model(chemicals, 'mu', IdealTPMixtureModel)
         V = create_mixture_model(chemicals, 'V', IdealTPMixtureModel)
         kappa = create_mixture_model(chemicals, 'kappa', IdealTPMixtureModel)
         Hvap = IdealHvapModel(chemicals)
         sigma = SinglePhaseIdealTMixtureModel([getfield(i, 'sigma') for i in chemicals], 'sigma')
         epsilon = SinglePhaseIdealTMixtureModel([getfield(i, 'epsilon') for i in chemicals], 'epsilon')
-        return cls(chemicals, Cn, mu, V, kappa, Hvap, sigma, epsilon, MWs)
+        return cls(chemicals, eos_chemicals, Cn, H, S, mu, V, kappa, Hvap, sigma, epsilon, MWs)
     
     @classmethod
     def subclass(cls, EOS, name=None):
         if name is None: name = EOS.__name__.replace('MIX', '') + 'Mixture'
         return type(name, (cls,), dict(EOS=EOS, cache={}))
 
+    def _info(self):
+        return (
+            f"{type(self).__name__}(\n"
+            f"    eos_chemicals=[{', '.join([i.ID for i in self.eos_chemicals])}]\n"
+             ")"
+        )
+
 dct = globals()    
 clsnames = []
 for name in ('PRMIX', 'SRKMIX', 'PR78MIX', 'VDWMIX', 'PRSVMIX',

From e4d5ccb6e8a61ccad2329c97a77be4519ba9dce4 Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 11:05:44 -0600
Subject: [PATCH 3/7] optimize

---
 thermosteam/mixture/mixture.py | 232 ++++++++++++++++++---------------
 1 file changed, 128 insertions(+), 104 deletions(-)

diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index 1ada29bc..f5ea06c3 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -204,59 +204,79 @@ def S(self, phase, mol, T, P):
             S += self._S_excess(phase, mol, T, P)
         return S
     
+    def _load_free_energy_args(self, *args): pass
+    def _load_xfree_energy_args(self, *args): pass
+    
     def solve_T_at_HP(self, phase, mol, H, T_guess, P):
         """Solve for temperature in Kelvin."""
-        args = (H, self.H, phase, mol, P, self.Cn, [0, None])
-        T_guess = flx.aitken(iter_T_at_HP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
-        T = iter_T_at_HP(T_guess, *args)
-        return (
-            flx.aitken_secant(
-                lambda T: self.H(phase, mol, T, P) - H,
-                x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+        self._load_free_energy_args(phase, mol, T_guess, P)
+        try:
+            args = (H, self.H, phase, mol, P, self.Cn, [0, None])
+            T_guess = flx.aitken(iter_T_at_HP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
+            T = iter_T_at_HP(T_guess, *args)
+            return (
+                flx.aitken_secant(
+                    lambda T: self.H(phase, mol, T, P) - H,
+                    x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+                )
+                if abs(T - T_guess) > self.T_tol else T
             )
-            if abs(T - T_guess) > self.T_tol else T
-        )
+        finally:
+            self._free_energy_args.clear()
         
     def xsolve_T_at_HP(self, phase_mol, H, T_guess, P):
         """Solve for temperature in Kelvin."""
         phase_mol = tuple(phase_mol)
-        args = (H, self.xH, phase_mol, P, self.xCn, [0, None])
-        T_guess = flx.aitken(xiter_T_at_HP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
-        T = xiter_T_at_HP(T_guess, *args)
-        return (
-            flx.aitken_secant(
-                lambda T: self.xH(phase_mol, T, P) - H,
-                x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+        self._load_xfree_energy_args(phase_mol, T_guess, P)
+        try:
+            free_energy_args = self._xfree_energy_args(phase_mol, T_guess, P)
+            args = (H, self.xH, phase_mol, P, self.xCn, [0, None], free_energy_args)
+            T_guess = flx.aitken(xiter_T_at_HP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
+            T = xiter_T_at_HP(T_guess, *args)
+            return (
+                flx.aitken_secant(
+                    lambda T: self.xH(phase_mol, T, P, free_energy_args) - H,
+                    x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+                )
+                if abs(T - T_guess) > self.T_tol else T
             )
-            if abs(T - T_guess) > self.T_tol else T
-        )
+        finally:
+            self._free_energy_args.clear()
     
     def solve_T_at_SP(self, phase, mol, S, T_guess, P):
         """Solve for temperature in Kelvin."""
-        args = (S, self.S, phase, mol, P, self.Cn, [0, None])
-        T_guess = flx.aitken(iter_T_at_SP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
-        T = iter_T_at_SP(T_guess, *args)
-        return (
-            flx.aitken_secant(
-                lambda T: self.S(phase, mol, T, P) - S,
-                x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+        self._load_free_energy_args(phase, mol, T_guess, P)
+        try:
+            args = (S, self.S, phase, mol, P, self.Cn, [0, None])
+            T_guess = flx.aitken(iter_T_at_SP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
+            T = iter_T_at_SP(T_guess, *args)
+            return (
+                flx.aitken_secant(
+                    lambda T: self.S(phase, mol, T, P) - S,
+                    x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+                )
+                if abs(T - T_guess) > self.T_tol else T
             )
-            if abs(T - T_guess) > self.T_tol else T
-        )
+        finally:
+            self._free_energy_args.clear()
         
     def xsolve_T_at_SP(self, phase_mol, S, T_guess, P):
         """Solve for temperature in Kelvin."""
         phase_mol = tuple(phase_mol)
-        args = (S, self.xS, phase_mol, P, self.xCn, [0, None])
-        T_guess = flx.aitken(xiter_T_at_SP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
-        T = xiter_T_at_SP(T_guess, *args)
-        return (
-            flx.aitken_secant(
-                lambda T: self.xS(phase_mol, T, P) - S,
-                x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+        self._load_xfree_energy_args(phase_mol, T_guess, P)
+        try:
+            args = (S, self.xS, phase_mol, P, self.xCn, [0, None])
+            T_guess = flx.aitken(xiter_T_at_SP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
+            T = xiter_T_at_SP(T_guess, *args)
+            return (
+                flx.aitken_secant(
+                    lambda T: self.xS(phase_mol, T, P) - S,
+                    x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
+                )
+                if abs(T - T_guess) > self.T_tol else T
             )
-            if abs(T - T_guess) > self.T_tol else T
-        )
+        finally:
+            self._free_energy_args.clear()
     
     def xCn(self, phase_mol, T, P=None):
         """Multi-phase mixture molar isobaric heat capacity [J/mol/K]."""
@@ -383,6 +403,7 @@ class IdealMixture(Mixture):
         'Cn', 'mu', 'V', 'kappa',
         'Hvap', 'sigma', 'epsilon',
         'MWs', '_H', '_H_excess', '_S', '_S_excess',
+        '_free_energy_args',
     )
     
     def __init__(self, Cn, H, S, H_excess, S_excess,
@@ -401,6 +422,7 @@ def __init__(self, Cn, H, S, H_excess, S_excess,
         self._S = S
         self._H_excess = H_excess
         self._S_excess = S_excess
+        self._free_energy_args = {}
     
     @classmethod
     def from_chemicals(cls, chemicals, 
@@ -482,6 +504,7 @@ class EOSMixture(Mixture):
     __slots__ = (
         'chemicals', 'eos_chemicals', 'Cn_ideal', 'mu', 'V', 'kappa',
         'Hvap', 'sigma', 'epsilon', 'H_ideal', 'S_ideal', 'MWs', 
+        '_free_energy_args',
     )
     chemsep_db = None
     
@@ -498,55 +521,70 @@ def __init__(self, chemicals, eos_chemicals, Cn_ideal, H_ideal, S_ideal, mu, V,
         self.sigma = sigma
         self.epsilon = epsilon
         self.MWs = MWs
+        self._free_energy_args = {}
     
-    def eos(self, phase, chemicals, zs, T, P):
-        key = (phase, chemicals)
+    def eos_args(self, phase, mol, T, P):
+        chemicals = self.chemicals
+        dct = mol.dct
+        eos_chemicals = self.eos_chemicals
+        chemical_subset = []
+        mol_subset = []
+        eos_mol = 0
+        for i, j in dct.items():
+            chemical = chemicals[i]
+            if chemical in eos_chemicals:
+                chemical_subset.append(chemical)
+                mol_subset.append(j)
+                eos_mol += j
+        zs = [i / eos_mol for i in mol_subset]
+        eos_chemicals = tuple(eos_chemicals)
+        key = (phase, eos_chemicals)
         cache = self.cache
+        only_g = phase == 'g'
+        only_l = phase == 'l'
         if key in cache:
-            return cache[key]
+            eos = cache[key].to_TP_zs(
+                T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
+                fugacities=False
+            )
         else:
-            only_g = phase == 'g'
-            only_l = phase == 'l'
-            if only_g or only_l:
+            data = tmo.ChemicalData(eos_chemicals)
+            if self.chemsep_db is None:
+                kijs = None
+            else:
                 try:
-                    self._eos = eos = self._eos.to_TP_zs(
-                        T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
-                        fugacities=False
-                    )
+                    kijs = IPDB.get_ip_asymmetric_matrix(self.chemsep_db, data.CASs, 'kij')
                 except:
-                    data = tmo.ChemicalData(chemicals)
-                    if self.chemsep_db is None:
-                        kijs = None
-                    else:
-                        try:
-                            kijs = IPDB.get_ip_asymmetric_matrix(self.chemsep_db, data.CASs, 'kij')
-                        except:
-                            kijs = None
-                    self._eos = eos = self.EOS(
-                        zs=zs, T=T, P=P, Tcs=data.Tcs, Pcs=data.Pcs, omegas=data.omegas, kijs=kijs,
-                        only_g=only_g, only_l=only_l, fugacities=False, 
-                    )    
-        return eos
+                    kijs = None
+            self.cache[key] = eos = self.EOS(
+                Tcs=data.Tcs, Pcs=data.Pcs, omegas=data.omegas, kijs=kijs,
+                T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
+                fugacities=False
+            )
+        return eos, eos_mol, dict(
+            P=P, zs=zs, only_g=only_g, only_l=only_l,
+            fugacities=False
+        )
+
+    def _load_free_energy_args(self, phase, mol, T, P):
+        self._free_energy_args[phase] = self.eos_args(phase, mol, T, P)
+    
+    def _load_xfree_energy_args(self, phase_mol, T, P):
+        fea = self._free_energy_args
+        for phase, mol in phase_mol: fea[phase] = self.eos_args(phase, mol, T, P)
     
     def Cn(self, phase, mol, T, P):
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
-        chemicals = self.chemicals
-        dct = mol.dct
         Cn = self.Cn_ideal(phase, mol, T, P)
         if phase != 's':
-            eos_chemicals = self.eos_chemicals
-            chemical_subset = []
-            mol_subset = []
-            eos_mol = 0
-            for i, j in dct.items():
-                chemical = chemicals[i]
-                if chemical in eos_chemicals:
-                    chemical_subset.append(chemical)
-                    mol_subset.append(j)
-                    eos_mol += j
-            zs = [i / eos_mol for i in mol_subset]
-            eos = self.eos(
-                phase, tuple(chemical_subset), zs, T, P
+            if phase in self._free_energy_args:
+                eos, eos_mol, eos_kwargs = self._free_energy_args[phase]
+            else:
+                eos, eos_mol, eos_kwargs = self.eos_args(
+                    phase, mol, T, P
+                )
+            eos = eos.to_TP_zs(
+                T=T, **eos_kwargs
             )
             if phase == 'l':
                 Cn += eos.Cp_dep_l * eos_mol
@@ -557,23 +595,16 @@ def Cn(self, phase, mol, T, P):
     def H(self, phase, mol, T, P):
         """Return enthalpy [J/mol]."""
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
-        chemicals = self.chemicals
-        dct = mol.dct
         H = self.H_ideal(phase, mol, T, P)
         if phase != 's':
-            eos_chemicals = self.eos_chemicals
-            chemical_subset = []
-            mol_subset = []
-            eos_mol = 0
-            for i, j in dct.items():
-                chemical = chemicals[i]
-                if chemical in eos_chemicals:
-                    chemical_subset.append(chemical)
-                    mol_subset.append(j)
-                    eos_mol += j
-            zs = [i / eos_mol for i in mol_subset]
-            eos = self.eos(
-                phase, tuple(chemical_subset), zs, T, P
+            if phase in self._free_energy_args:
+                eos, eos_mol, eos_kwargs = self._free_energy_args[phase]
+            else:
+                eos, eos_mol, eos_kwargs = self.eos_args(
+                    phase, mol, T, P
+                )
+            eos = eos.to_TP_zs(
+                T=T, **eos_kwargs
             )
             if phase == 'l':
                 H += eos.H_dep_l * eos_mol
@@ -584,23 +615,16 @@ def H(self, phase, mol, T, P):
     def S(self, phase, mol, T, P):
         """Return entropy [J/mol/K]."""
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
-        chemicals = self.chemicals
-        dct = mol.dct
         S = self.S_ideal(phase, mol, T, P)
         if phase != 's':
-            eos_chemicals = self.eos_chemicals
-            chemical_subset = []
-            mol_subset = []
-            eos_mol = 0
-            for i, j in dct.items():
-                chemical = chemicals[i]
-                if chemical in eos_chemicals:
-                    chemical_subset.append(chemical)
-                    mol_subset.append(j)
-                    eos_mol += j
-            zs = [i / eos_mol for i in mol_subset]
-            eos = self.eos(
-                phase, tuple(chemical_subset), zs, T, P
+            if phase in self._free_energy_args:
+                eos, eos_mol, eos_kwargs = self._free_energy_args[phase]
+            else:
+                eos, eos_mol, eos_kwargs = self.eos_args(
+                    phase, mol, T, P
+                )
+            eos = eos.to_TP_zs(
+                T=T, **eos_kwargs
             )
             if phase == 'l':
                 S += eos.S_dep_l * eos_mol

From 311cc767cdb263afcaf925f0ceab5b077a6f1bcf Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 14:39:49 -0600
Subject: [PATCH 4/7] more robust

---
 thermosteam/__init__.py                       |  8 +--
 thermosteam/_thermo.py                        | 15 ++---
 .../poyinting_correction_factors.py           | 19 ++++++-
 thermosteam/equilibrium/vle.py                |  6 +-
 thermosteam/mixture/mixture.py                | 55 ++++++++++---------
 5 files changed, 59 insertions(+), 44 deletions(-)

diff --git a/thermosteam/__init__.py b/thermosteam/__init__.py
index 303e95f1..c52db1ce 100644
--- a/thermosteam/__init__.py
+++ b/thermosteam/__init__.py
@@ -51,15 +51,14 @@
 from ._chemical import Chemical
 from ._chemicals import Chemicals, CompiledChemicals
 from ._thermal_condition import ThermalCondition
-from . import mixture
-from .mixture import Mixture
 from ._thermo import Thermo, IdealThermo
 from ._settings import settings, ProcessSettings
 from ._thermo_data import ThermoData
 from . import (
     indexer,
     reaction,
-    equilibrium
+    equilibrium,
+    mixture,
 )
 from ._stream import Stream
 from ._heat_and_power import Heat, Power
@@ -67,10 +66,11 @@
 from .base import functor
 from .reaction import *
 from .equilibrium import * 
+from .mixture import *
 
 __all__ = ('Chemical', 'ChemicalData', 'Chemicals', 'CompiledChemicals', 'Mixture', 'Thermo', 
            'IdealThermo', 'Stream', 'MultiStream', 'Heat', 'Power', 'ThermalCondition', 'ProcessSettings',
-           'mixture', 'ThermoData', *reaction.__all__, *equilibrium.__all__, 
+           'mixture', 'ThermoData', *reaction.__all__, *equilibrium.__all__,  *mixture.__all__,
            'indexer', 'settings', 'functor', 'functors', 'chemicals', 'base', 
            'equilibrium', 'units_of_measure', 'exceptions', 'functional', 
            'reaction', 'constants', 'utils', 'separations')
diff --git a/thermosteam/_thermo.py b/thermosteam/_thermo.py
index d19bd38f..aa4aa26f 100644
--- a/thermosteam/_thermo.py
+++ b/thermosteam/_thermo.py
@@ -49,8 +49,7 @@ class Thermo:
     >>> thermo.show() 
     Thermo(
         chemicals=CompiledChemicals([Ethanol, Water]),
-        mixture=Mixture(
-            rule='ideal', ...
+        mixture=IdealMixture(
             include_excess_energies=False
         ),
         Gamma=DortmundActivityCoefficients,
@@ -66,8 +65,7 @@ class Thermo:
     >>> ideal.show()
     IdealThermo(
         chemicals=CompiledChemicals([Ethanol, Water]),
-        mixture=Mixture(
-            rule='ideal', ...
+        mixture=IdealMixture(
             include_excess_energies=False
         ),
     )
@@ -77,7 +75,7 @@ class Thermo:
     >>> # Ideal
     >>> tmo.settings.set_thermo(ideal)
     >>> stream = tmo.Stream('stream', Water=100, Ethanol=100)
-    >>> stream.vle(T=361, P=101325)
+    >>> stream.vle(T=360, P=101325)
     >>> stream.show()
     MultiStream: stream
     phases: ('g', 'l'), T: 361 K, P: 101325 Pa
@@ -88,10 +86,10 @@ class Thermo:
     >>> # Modified Roult's law:                 
     >>> tmo.settings.set_thermo(thermo)
     >>> stream = tmo.Stream('stream', Water=100, Ethanol=100)
-    >>> stream.vle(T=360, P=101325) 
+    >>> stream.vle(T=361, P=101325) 
     >>> stream.show()
     MultiStream: stream
-    phases: ('g', 'l'), T: 360 K, P: 101325 Pa
+    phases: ('g', 'l'), T: 361 K, P: 101325 Pa
     flow (kmol/hr): (g) Ethanol  100
                         Water    100
     
@@ -102,8 +100,7 @@ class Thermo:
     >>> thermo.show()
     Thermo(
         chemicals=CompiledChemicals([Ethanol, Water]),
-        mixture=Mixture(
-            rule='ideal', ...
+        mixture=IdealMixture(
             include_excess_energies=False
         ),
         Gamma=DortmundActivityCoefficients,
diff --git a/thermosteam/equilibrium/poyinting_correction_factors.py b/thermosteam/equilibrium/poyinting_correction_factors.py
index 3cdab344..41e07a54 100644
--- a/thermosteam/equilibrium/poyinting_correction_factors.py
+++ b/thermosteam/equilibrium/poyinting_correction_factors.py
@@ -84,7 +84,20 @@ class IdealGasPoyintingCorrectionFactors(PoyintingCorrectionFactors):
     __slots__ = ('_chemicals',)
 
     def __call__(self, T, P, Psats=None):
-        vls = np.array([i.V('l', T, P) for i in self._chemicals], dtype=float)
+        vls = []
+        chemicals = []
+        index = []
+        for i, chemical in enumerate(self._chemicals):
+            try: vls.append(chemical.V('l', T, P))
+            except: continue
+            chemicals.append(chemical)
+            index.append(i)
+        vls = np.array(vls, dtype=float)
         if Psats is None:
-            Psats = np.array([i.Psat(T) for i in self._chemicals], dtype=float)
-        return ideal_gas_poyinting_correction_factors(T, P, vls, Psats)
\ No newline at end of file
+            Psats = np.array([i.Psat(T) for i in chemicals], dtype=float)
+        else:
+            Psats = Psats[index]
+        pcf = np.ones(len(self._chemicals))
+        pcf[index] = ideal_gas_poyinting_correction_factors(T, P, vls, Psats)
+        return pcf
+        
\ No newline at end of file
diff --git a/thermosteam/equilibrium/vle.py b/thermosteam/equilibrium/vle.py
index 00df23ef..7da6a820 100644
--- a/thermosteam/equilibrium/vle.py
+++ b/thermosteam/equilibrium/vle.py
@@ -757,7 +757,8 @@ def set_TV(self, T, V):
             
             self._P = thermal_condition.P = P
             set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol)
-            self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass
+            try: self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass
+            except: pass
 
     def set_TH(self, T, H):
         self._setup()
@@ -914,7 +915,8 @@ def set_PV(self, P, V):
                     v = self._v
             self._T = thermal_condition.T = T
             set_flows(vapor_mol, liquid_mol, index, v, mol)
-            self._H_hat = self.mixture.xH(self._phase_data, T, P)/self._F_mass
+            try: self._H_hat = self.mixture.xH(self._phase_data, T, P)/self._F_mass
+            except: pass
     
     def set_PS(self, P, S, stacklevel=0):
         self._setup()
diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index f5ea06c3..c7e04c05 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -114,6 +114,7 @@ class Mixture:
     
     def MW(self, mol):
         """Return molecular weight [g/mol] given molar array [mol]."""
+        if mol.__class__ is not SparseVector: mol = SparseVector(mol)
         total_mol = mol.sum()
         return (mol * self.MWs).sum() / total_mol if total_mol else 0.
     
@@ -284,22 +285,16 @@ def xCn(self, phase_mol, T, P=None):
     
     def xH(self, phase_mol, T, P):
         """Multi-phase mixture enthalpy [J/mol]."""
-        H = self._H
+        H = self.H
         phase_mol = tuple(phase_mol)
         H_total = sum([H(phase, mol, T, P) for phase, mol in phase_mol])
-        if self.include_excess_energies:
-            H_excess = self._H_excess
-            H_total += sum([H_excess(phase, mol, T, P) for phase, mol in phase_mol])
         return H_total
     
     def xS(self, phase_mol, T, P):
         """Multi-phase mixture entropy [J/mol/K]."""
-        S = self._S
+        S = self.S
         phase_mol = tuple(phase_mol)
         S_total = sum([S(phase, mol, T, P) for phase, mol in phase_mol])
-        if self.include_excess_energies:
-            S_excess = self._S_excess
-            S_total += sum([S_excess(phase, mol, T, P) for phase, mol in phase_mol])
         return S_total
     
     def xV(self, phase_mol, T, P):
@@ -469,7 +464,7 @@ def from_chemicals(cls, chemicals,
             MWs = chemicals.MW
             chemicals = chemicals.tuple
         else:
-            chemicals = [(i if isa(i, Chemical) else Chemical(i)) for i in chemicals]
+            chemicals = [(i if isa(i, Chemical) else Chemical(i, cache=cache)) for i in chemicals]
             MWs = chemical_data_array(chemicals, 'MW')
         getfield = getattr
         Cn =  create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
@@ -537,7 +532,7 @@ def eos_args(self, phase, mol, T, P):
                 mol_subset.append(j)
                 eos_mol += j
         zs = [i / eos_mol for i in mol_subset]
-        eos_chemicals = tuple(eos_chemicals)
+        eos_chemicals = tuple(chemical_subset)
         key = (phase, eos_chemicals)
         cache = self.cache
         only_g = phase == 'g'
@@ -576,7 +571,7 @@ def _load_xfree_energy_args(self, phase_mol, T, P):
     def Cn(self, phase, mol, T, P):
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
         Cn = self.Cn_ideal(phase, mol, T, P)
-        if phase != 's':
+        if phase != 's' and mol.dct:
             if phase in self._free_energy_args:
                 eos, eos_mol, eos_kwargs = self._free_energy_args[phase]
             else:
@@ -587,16 +582,18 @@ def Cn(self, phase, mol, T, P):
                 T=T, **eos_kwargs
             )
             if phase == 'l':
-                Cn += eos.Cp_dep_l * eos_mol
+                try: Cn += eos.Cn_dep_l * eos_mol
+                except: Cn += eos.Cn_dep_g * eos_mol
             else:
-                Cn += eos.Cp_dep_g * eos_mol
+                try: Cn += eos.Cn_dep_g * eos_mol
+                except: Cn += eos.Cn_dep_l * eos_mol
         return Cn
     
     def H(self, phase, mol, T, P):
         """Return enthalpy [J/mol]."""
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
         H = self.H_ideal(phase, mol, T, P)
-        if phase != 's':
+        if phase != 's' and mol.dct:
             if phase in self._free_energy_args:
                 eos, eos_mol, eos_kwargs = self._free_energy_args[phase]
             else:
@@ -607,16 +604,18 @@ def H(self, phase, mol, T, P):
                 T=T, **eos_kwargs
             )
             if phase == 'l':
-                H += eos.H_dep_l * eos_mol
+                try: H += eos.H_dep_l * eos_mol
+                except: H += eos.H_dep_g * eos_mol
             else:
-                H += eos.H_dep_g * eos_mol
+                try: H += eos.H_dep_g * eos_mol
+                except: H += eos.H_dep_l * eos_mol
         return H
     
     def S(self, phase, mol, T, P):
         """Return entropy [J/mol/K]."""
         if mol.__class__ is not SparseVector: mol = SparseVector(mol)
         S = self.S_ideal(phase, mol, T, P)
-        if phase != 's':
+        if phase != 's' and mol.dct:
             if phase in self._free_energy_args:
                 eos, eos_mol, eos_kwargs = self._free_energy_args[phase]
             else:
@@ -627,15 +626,17 @@ def S(self, phase, mol, T, P):
                 T=T, **eos_kwargs
             )
             if phase == 'l':
-                S += eos.S_dep_l * eos_mol
+                try: S += eos.S_dep_l * eos_mol
+                except: S += eos.S_dep_g * eos_mol
             else:
-                S += eos.S_dep_g * eos_mol
+                try: S += eos.S_dep_g * eos_mol
+                except: S += eos.S_dep_l * eos_mol
         return S
     
     @classmethod
-    def from_chemicals(cls, chemicals, eos_chemicals, cache=True):
+    def from_chemicals(cls, chemicals, eos_chemicals=None, cache=True):
         """
-        Create a Mixture object from chemical objects.
+        Create a EOSMixture object from chemical objects.
         
         Parameters
         ----------
@@ -650,15 +651,15 @@ def from_chemicals(cls, chemicals, eos_chemicals, cache=True):
         --------
         Calculate enthalpy of evaporation for a water and ethanol mixture:
         
-        >>> from thermosteam import Mixture
-        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> from thermosteam import PRMixture
+        >>> mixture = PRMixture.from_chemicals(['Water', 'Ethanol'])
         >>> mixture.Hvap([0.2, 0.8], 350)
         39750.62
 
         Calculate density for a water and ethanol mixture in g/L:
 
-        >>> from thermosteam import Mixture
-        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> from thermosteam import PRMixture
+        >>> mixture = PRMixture.from_chemicals(['Water', 'Ethanol'])
         >>> mixture.get_property('rho', 'g/L', 'l', [0.2, 0.8], 350, 101325)
         754.23
         
@@ -668,8 +669,10 @@ def from_chemicals(cls, chemicals, eos_chemicals, cache=True):
             MWs = chemicals.MW
             chemicals = chemicals.tuple
         else:
-            chemicals = [(i if isa(i, Chemical) else Chemical(i)) for i in chemicals]
+            chemicals = [(i if isa(i, Chemical) else Chemical(i, cache=cache)) for i in chemicals]
             MWs = chemical_data_array(chemicals, 'MW')
+        if eos_chemicals is None:
+            eos_chemicals = tuple(chemicals)
         getfield = getattr
         Cn = create_mixture_model(chemicals, 'Cn', IdealTMixtureModel)
         H = create_mixture_model(chemicals, 'H', IdealTPMixtureModel)

From e18feb873fd908fc19f84fbc942695b273e7366d Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 19:30:32 -0600
Subject: [PATCH 5/7] fix default rep

---
 thermosteam/mixture/mixture.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index c7e04c05..1c74471a 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -482,17 +482,16 @@ def from_chemicals(cls, chemicals,
                    mu, V, kappa, Hvap, sigma, epsilon, MWs, include_excess_energies)
     
     def __repr__(self):
-        return f"{type(self).__name__}(rule={repr(self.rule)}, ..., include_excess_energies={self.include_excess_energies})"
+        return f"{type(self).__name__}(...)"
     
     def _info(self):
-        return (f"{type(self).__name__}(\n"
-                f"    include_excess_energies={self.include_excess_energies}\n"
-                 ")")
+        return (f"{type(self).__name__}(...)")
     
     def show(self):
         print(self._info())
     _ipython_display_ = show
     
+    
 # %% Thermo mixture
 
 class EOSMixture(Mixture):

From 84add4c4066289387452079f5636f001e94bfe1e Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 19:36:10 -0600
Subject: [PATCH 6/7] fix recent bug

---
 thermosteam/mixture/mixture.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index 1c74471a..2e778ff9 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -230,13 +230,12 @@ def xsolve_T_at_HP(self, phase_mol, H, T_guess, P):
         phase_mol = tuple(phase_mol)
         self._load_xfree_energy_args(phase_mol, T_guess, P)
         try:
-            free_energy_args = self._xfree_energy_args(phase_mol, T_guess, P)
-            args = (H, self.xH, phase_mol, P, self.xCn, [0, None], free_energy_args)
+            args = (H, self.xH, phase_mol, P, self.xCn, [0, None])
             T_guess = flx.aitken(xiter_T_at_HP, T_guess, self.T_tol, args, self.maxiter, checkiter=False)
             T = xiter_T_at_HP(T_guess, *args)
             return (
                 flx.aitken_secant(
-                    lambda T: self.xH(phase_mol, T, P, free_energy_args) - H,
+                    lambda T: self.xH(phase_mol, T, P) - H,
                     x0=T_guess, x1=T, xtol=self.T_tol, ytol=0.
                 )
                 if abs(T - T_guess) > self.T_tol else T

From ce79be9b9324cd4a7bddebebf3fc4c2bb5ca46b0 Mon Sep 17 00:00:00 2001
From: Yoel Cortes-Pena <yoelcortes@gmail.com>
Date: Sat, 27 Jan 2024 20:14:40 -0600
Subject: [PATCH 7/7] finishing touches

---
 thermosteam/_settings.py       |  6 ++----
 thermosteam/_thermo.py         |  8 ++++----
 thermosteam/mixture/mixture.py | 24 ++++++++++--------------
 3 files changed, 16 insertions(+), 22 deletions(-)

diff --git a/thermosteam/_settings.py b/thermosteam/_settings.py
index 308c8813..fa6c9e95 100644
--- a/thermosteam/_settings.py
+++ b/thermosteam/_settings.py
@@ -63,8 +63,7 @@ class ProcessSettings:
     >>> settings.thermo.show()
     Thermo(
         chemicals=CompiledChemicals([Water]),
-        mixture=Mixture(
-            rule='ideal', ...
+        mixture=IdealMixture(...
             include_excess_energies=False
         ),
         Gamma=DortmundActivityCoefficients,
@@ -80,8 +79,7 @@ class ProcessSettings:
     Access defined mixture property algorithm:
         
     >>> settings.mixture.show()
-    Mixture(
-        rule='ideal', ...
+    IdealMixture(...
         include_excess_energies=False
     )
 
diff --git a/thermosteam/_thermo.py b/thermosteam/_thermo.py
index aa4aa26f..ca954c90 100644
--- a/thermosteam/_thermo.py
+++ b/thermosteam/_thermo.py
@@ -49,7 +49,7 @@ class Thermo:
     >>> thermo.show() 
     Thermo(
         chemicals=CompiledChemicals([Ethanol, Water]),
-        mixture=IdealMixture(
+        mixture=IdealMixture(...
             include_excess_energies=False
         ),
         Gamma=DortmundActivityCoefficients,
@@ -65,7 +65,7 @@ class Thermo:
     >>> ideal.show()
     IdealThermo(
         chemicals=CompiledChemicals([Ethanol, Water]),
-        mixture=IdealMixture(
+        mixture=IdealMixture(...
             include_excess_energies=False
         ),
     )
@@ -75,7 +75,7 @@ class Thermo:
     >>> # Ideal
     >>> tmo.settings.set_thermo(ideal)
     >>> stream = tmo.Stream('stream', Water=100, Ethanol=100)
-    >>> stream.vle(T=360, P=101325)
+    >>> stream.vle(T=361, P=101325)
     >>> stream.show()
     MultiStream: stream
     phases: ('g', 'l'), T: 361 K, P: 101325 Pa
@@ -100,7 +100,7 @@ class Thermo:
     >>> thermo.show()
     Thermo(
         chemicals=CompiledChemicals([Ethanol, Water]),
-        mixture=IdealMixture(
+        mixture=IdealMixture(...
             include_excess_energies=False
         ),
         Gamma=DortmundActivityCoefficients,
diff --git a/thermosteam/mixture/mixture.py b/thermosteam/mixture/mixture.py
index 2e778ff9..07976cf3 100644
--- a/thermosteam/mixture/mixture.py
+++ b/thermosteam/mixture/mixture.py
@@ -308,14 +308,12 @@ def xkappa(self, phase_mol, T, P):
         """Multi-phase mixture thermal conductivity [W/m/K]."""
         return sum([self.kappa(phase, mol, T, P) for phase, mol in phase_mol])
     
+        
     def __repr__(self):
-        return f"{type(self).__name__}(..., include_excess_energies={self.include_excess_energies})"
+        return f"{type(self).__name__}(...)"
     
     def _info(self):
-        return (f"{type(self).__name__}(\n"
-                f"    rule={repr(self.rule)}, ...\n"
-                f"    include_excess_energies={self.include_excess_energies}\n"
-                 ")")
+        return (f"{type(self).__name__}(...)")
     
     def show(self):
         print(self._info())
@@ -445,15 +443,15 @@ def from_chemicals(cls, chemicals,
         --------
         Calculate enthalpy of evaporation for a water and ethanol mixture:
         
-        >>> from thermosteam import Mixture
-        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> from thermosteam import IdealMixture
+        >>> mixture = IdealMixture.from_chemicals(['Water', 'Ethanol'])
         >>> mixture.Hvap([0.2, 0.8], 350)
         39750.62
 
         Calculate density for a water and ethanol mixture in g/L:
 
         >>> from thermosteam import Mixture
-        >>> mixture = Mixture.from_chemicals(['Water', 'Ethanol'])
+        >>> mixture = IdealMixture.from_chemicals(['Water', 'Ethanol'])
         >>> mixture.get_property('rho', 'g/L', 'l', [0.2, 0.8], 350, 101325)
         754.23
         
@@ -481,14 +479,12 @@ def from_chemicals(cls, chemicals,
                    mu, V, kappa, Hvap, sigma, epsilon, MWs, include_excess_energies)
     
     def __repr__(self):
-        return f"{type(self).__name__}(...)"
+        return f"{type(self).__name__}(..., include_excess_energies={self.include_excess_energies})"
     
     def _info(self):
-        return (f"{type(self).__name__}(...)")
-    
-    def show(self):
-        print(self._info())
-    _ipython_display_ = show
+        return (f"{type(self).__name__}(...\n"
+                f"    include_excess_energies={self.include_excess_energies}\n"
+                 ")")
     
     
 # %% Thermo mixture