diff --git a/thermosteam/_chemical_data.py b/thermosteam/_chemical_data.py
index 5fa1a4be..a5bd3bc9 100644
--- a/thermosteam/_chemical_data.py
+++ b/thermosteam/_chemical_data.py
@@ -7,12 +7,31 @@
 from collections.abc import Iterable
 import thermosteam as tmo
 
+class CachedProperty:
+    __slots__ = ('fget', 'key')
+    
+    def __new__(cls, fget):
+        self = super().__new__(cls)
+        self.fget = fget
+        self.key = fget.__name__
+        return self
+    
+    def __get__(self, data, _):
+        cache = data.cache
+        key = self.key
+        if key in cache:
+            value = cache[key]
+        else:
+            cache[key] = value = self.fget(data)
+        return value
+
+
 class ChemicalData:
     """
     Create a ChemicalData object that works as thermo ChemicalPropertiesPackage
     and ChemicalConstantsPackage objects.
     """
-    __slots__ = ('N', 'chemicals', 'index', 'cmps')
+    __slots__ = ('N', 'chemicals', 'index', 'cmps', 'cache')
     
     def __init__(self, chemicals):
         if isinstance(chemicals, tmo.CompiledChemicals):
@@ -23,6 +42,7 @@ def __init__(self, chemicals):
             self.chemicals = tuple(chemicals)
             self.N = len(chemicals)
         self.cmps = range(self.N)
+        self.cache = {}
     
     def __iter__(self):
         return self.chemicals.__iter__()
@@ -58,213 +78,227 @@ def __getitem__(self, IDs):
         else:
             raise KeyError("key must be a string or a tuple, not a '{type(IDs).__name__}' object")
     
-    @property
+    @CachedProperty
     def EnthalpyVaporizations(self):
         return [i.Hvap for i in self.chemicals]
-    @property
+    @CachedProperty
     def VaporPressures(self):
         return [i.Psat for i in self.chemicals]
-    @property
+    @CachedProperty
     def SublimationPressures(self):
         return [i.Psub for i in self.chemicals]
-    @property
+    @CachedProperty
     def EnthalpySublimations(self):
         return [i.Hsub for i in self.chemicals]
-    @property
+    @CachedProperty
     def SurfaceTensions(self):
         return [i.sigma for i in self.chemicals]
-    @property
+    @CachedProperty
     def PermittivityLiquids(self):
         return [i.epsilon for i in self.chemicals]
-    @property
+    @CachedProperty
     def CASs(self):
         return [i.CAS for i in self.chemicals]
-    @property
+    @CachedProperty
     def MWs(self):
         return [i.MW for i in self.chemicals]
-    @property
+    @CachedProperty
     def Tms(self):
         return [i.Tm for i in self.chemicals]
-    @property
+    @CachedProperty
     def Tbs(self):
         return [i.Tb for i in self.chemicals]
-    @property
+    @CachedProperty
     def Tcs(self):
         return [i.Tc for i in self.chemicals]
-    @property
+    @CachedProperty
     def Pcs(self):
         return [i.Pc for i in self.chemicals]
-    @property
+    @CachedProperty
     def Vcs(self):
         return [i.Vc for i in self.chemicals]
-    @property
+    @CachedProperty
     def omegas(self):
         return [i.omega for i in self.chemicals]
-    @property
+    @CachedProperty
     def Zcs(self):
         return [i.Zc for i in self.chemicals]
-    @property
+    @CachedProperty
     def Hfus_Tms(self):
         return [i.Hfus for i in self.chemicals]   
-    @property
+    @CachedProperty
     def Tts(self):
         return [i.Tt for i in self.chemicals]
-    @property
+    @CachedProperty
     def Pts(self):
         return [i.Pt for i in self.chemicals]
-    @property
+    @CachedProperty
     def dipoles(self):
         return [i.dipole for i in self.chemicals]
-    @property
+    @CachedProperty
     def charges(self):
         return [i.charge for i in self.chemicals]
-    @property
+    @CachedProperty
     def UNIFAC_groups(self):
         return [i.UNIFAC for i in self.chemicals]
-    @property
+    @CachedProperty
     def UNIFAC_Dortmund_groups(self):
         return [i.Dortmund for i in self.chemicals]
-    @property
+    @CachedProperty
     def PSRK_groups(self):
         return [i.PSRK for i in self.chemicals]
-    @property
+    @CachedProperty
     def similarity_variables(self):
         return [i.similarity_variable for i in self.chemicals]
-    @property
+    @CachedProperty
     def StielPolars(self):
         return [i.Stiel_Polar for i in self.chemicals]
-    @property
+    @CachedProperty
     def VolumeGases(self):
         return [i.V if i._locked_state else i.V.g for i in self.chemicals]
-    @property
+    @CachedProperty
     def VolumeLiquids(self):
         return [i.V if i._locked_state else i.V.l for i in self.chemicals]
-    @property
+    @CachedProperty
     def VolumeSolids(self):
         return [i.V if i._locked_state else i.V.s for i in self.chemicals]
-    @property
+    @CachedProperty
     def HeatCapacityGases(self):
         return [i.Cn if i._locked_state else i.Cn.g for i in self.chemicals]
-    @property
+    @CachedProperty
     def HeatCapacityLiquids(self):
         return [i.Cn if i._locked_state else i.Cn.l for i in self.chemicals]
-    @property
+    @CachedProperty
     def HeatCapacitySolids(self):
         return [i.Cn if i._locked_state else i.Cn.s for i in self.chemicals]
-    @property
+    @CachedProperty
     def ViscosityGases(self):
         return [i.mu if i._locked_state else i.mu.g for i in self.chemicals]
-    @property
+    @CachedProperty
     def ViscosityLiquids(self):
         return [i.mu if i._locked_state else i.mu.l for i in self.chemicals]
-    @property
+    @CachedProperty
     def ThermalConductivityGases(self):
         return [i.k if i._locked_state else i.k.g for i in self.chemicals]
-    @property
+    @CachedProperty
     def ThermalConductivityLiquids(self):
         return [i.k if i._locked_state else i.k.l for i in self.chemicals]
-    @property
+    @CachedProperty
     def rhocs(self):
         return [None if Vc is None else 1.0/Vc for Vc in self.chemicals.Vcs]
-    @property
+    @CachedProperty
     def rhocs_mass(self):
         return [None if None in (rhoc, MW) else rhoc * 1e-3 * MW
                 for rhoc, MW in zip(self.rhocs, self.MWs)]
-    @property
+    @CachedProperty
     def Hfus_Tms_mass(self):
         return [None if None in (Hfus, MW) else Hfus * 1000.0 / MW
                 for Hfus, MW in zip(self.Hfus_Tms, self.MWs)]
-    @property
+    @CachedProperty
     def Hvap_Tbs(self):
         return [Hvap(Tb) for Hvap, Tb in zip(self.EnthalpyVaporizations, self.Tbs)]
-    @property
+    @CachedProperty
     def Hvap_Tbs_mass(self):
         return [None if None in (Hvap, MW) else Hvap*1000.0/MW
                 for Hvap, MW in zip(self.Hvap_Tbs, self.Tbs)]
-    @property
+    @CachedProperty
     def Vmg_STPs(self):
         return [i.TP_dependent_property(298.15, 101325.0)
                 for i in self.chemicals.VolumeGases]
-    @property
+    @CachedProperty
     def rhog_STPs(self):
         return [1.0/V if V is not None else None for V in self.chemicals.Vmg_STPs]
-    @property
+    @CachedProperty
     def rhog_STPs_mass(self):
         return [1e-3*MW/V if V is not None else None for V, MW in zip(self.Vmg_STPs, self.MWs)]
-    @property
+    @CachedProperty
     def Vml_Tbs(self):
         return [None if Tb is None else i.T_dependent_property(Tb) for i, Tb in zip(self.VolumeLiquids, self.Tbs)]
-    @property
+    @CachedProperty
     def Vml_Tms(self):
         return [None if Tm is None else i.T_dependent_property(Tm)
                 for i, Tm in zip(self.VolumeLiquids, self.Tms)]
-    @property
+    @CachedProperty
     def Vml_STPs(self):
         return [i.T_dependent_property(298.15) for i in self.chemicals.VolumeLiquids]
-    @property
+    @CachedProperty
     def Vml_60Fs(self):
         return [i.T_dependent_property(288.7055555555555) for i in self.chemicals.VolumeLiquids]
-    @property
+    @CachedProperty
     def rhol_STPs(self):
         return [1.0/V if V is not None else None for V in self.chemicals.Vml_STPs]
-    @property
+    @CachedProperty
     def rhol_60Fs(self):
         return [1.0/V if V is not None else None for V in self.chemicals.Vml_60Fs]
-    @property
+    @CachedProperty
     def rhol_STPs_mass(self):
         return [1e-3*MW/V if V is not None else None for V, MW in zip(self.Vml_STPs, self.MWs)]
-    @property
+    @CachedProperty
     def rhol_60Fs_mass(self):
         return [1e-3*MW/V if V is not None else None for V, MW in zip(self.Vml_60Fs, self.MWs)]
-    @property
+    @CachedProperty
     def Hsub_Tts(self):
         return [None if Tt is None else i(Tt)
                 for i, Tt in zip(self.EnthalpySublimations, self.Tts)]
-    @property
+    @CachedProperty
     def Hsub_Tts_mass(self):
         return [Hsub*1000.0/MW if Hsub is not None else None for Hsub, MW in zip(self.Hsub_Tts, self.MWs)]
-    @property
+    @CachedProperty
     def Stockmayers(self):
-        return [Stockmayer(Tm=self.Tms[i], Tb=self.Tbs[i], Tc=self.Tcs[i], 
-                           Zc=self.Zcs[i], omega=self.omegas[i], CASRN=self.CASs[i]) 
+        Tms = self.Tms
+        Tbs = self.Tbs
+        Tcs = self.Tcs
+        Zcs = self.Zcs
+        omegas = self.omegas
+        CASs = self.CASs
+        return [Stockmayer(Tm=Tms[i], Tb=Tbs[i], Tc=Tcs[i], 
+                           Zc=Zcs[i], omega=omegas[i], CASRN=CASs[i]) 
                 for i in self.cmps]
-    @property
+    @CachedProperty
     def molecular_diameters(self):
-        return [molecular_diameter(Tc=self.Tcs[i], Pc=self.Pcs[i], Vc=self.Vcs[i],
-                                   Zc=self.Zcs[i], omega=self.omegas[i],
-                                   Vm=self.Vml_Tms[i], Vb=self.Vml_Tbs[i], 
-                                   CASRN=self.CASs[i])
+        Tcs = self.Tcs
+        Zcs = self.Zcs
+        omegas = self.omegas
+        CASs = self.CASs
+        Vml_Tms = self.Vml_Tms
+        Vml_Tb = self.Vml_Tb
+        Pcs = self.Pcs
+        Vcs = self.Vcs
+        return [molecular_diameter(Tc=Tcs[i], Pc=Pcs[i], Vc=Vcs[i],
+                                   Zc=Zcs[i], omega=omegas[i],
+                                   Vm=Vml_Tms[i], Vb=Vml_Tb[i], 
+                                   CASRN=CASs[i])
                 for i in self.cmps]
-    @property
+    @CachedProperty
     def Psat_298s(self):
         return [i(298.15) for i in self.chemicals.VaporPressures]
-    @property
+    @CachedProperty
     def Hvap_298s(self):
         return [o.T_dependent_property(298.15) for o in self.chemicals.EnthalpyVaporizations]
-    @property
+    @CachedProperty
     def Hvap_298s_mass(self):
         return [Hvap*1000.0/MW if Hvap is not None else None for Hvap, MW in zip(self.Hvap_298s, self.MWs)]
-    @property
+    @CachedProperty
     def Vms_Tms(self):
         return [None if Tm is None else i.T_dependent_property(Tm) for i, Tm in zip(self.VolumeSolids, self.Tms)]
-    @property
+    @CachedProperty
     def rhos_Tms(self):
         return [1.0/V if V is not None else None for V in self.chemicals.Vms_Tms]
-    @property
+    @CachedProperty
     def rhos_Tms_mass(self):
         return [1e-3*MW/V if V is not None else None for V, MW in zip(self.Vms_Tms, self.MWs)]
-    @property
+    @CachedProperty
     def sigma_STPs(self):
         return [i.T_dependent_property(298.15) for i in self.SurfaceTensions]
-    @property
+    @CachedProperty
     def sigma_Tbs(self):
         return [None if Tb is None else i.T_dependent_property(Tb)
                 for i, Tb in zip(self.SurfaceTensions, self.Tbs)]
-    @property
+    @CachedProperty
     def sigma_Tms(self):
         return [None if Tm is None else i.T_dependent_property(Tm)
                for i, Tm in zip(self.SurfaceTensions, self.Tms)]
-    @property
+    @CachedProperty
     def Van_der_Waals_volumes(self):
         N = self.N
         N_range = self.cmps
@@ -275,7 +309,7 @@ def Van_der_Waals_volumes(self):
             if groups is not None:
                 UNIFAC_Rs[i] = UNIFAC_RQ(groups)[0]
         return [Van_der_Waals_volume(UNIFAC_Rs[i]) if UNIFAC_Rs[i] is not None else None for i in N_range]
-    @property
+    @CachedProperty
     def Van_der_Waals_areas(self):
         N = self.N
         N_range = self.cmps
@@ -286,7 +320,7 @@ def Van_der_Waals_areas(self):
             if groups is not None:
                 UNIFAC_Qs[i] = UNIFAC_RQ(groups)[1]
         return [Van_der_Waals_area(UNIFAC_Qs[i]) if UNIFAC_Qs[i] is not None else None for i in N_range]
-    @property
+    @CachedProperty
     def Parachors(self):
         N = self.N
         Parachors = [None]*N
diff --git a/thermosteam/equilibrium/activity_coefficients.py b/thermosteam/equilibrium/activity_coefficients.py
index f0905cae..1e2b2a9a 100644
--- a/thermosteam/equilibrium/activity_coefficients.py
+++ b/thermosteam/equilibrium/activity_coefficients.py
@@ -20,6 +20,11 @@
 from .unifac import DOUFSG, DOUFIP2016, UFIP, UFSG, NISTUFSG, NISTUFIP
 from numba import njit, prange
 from .ideal import ideal
+import thermosteam as tmo
+from thermo.interaction_parameters import IPDB
+from fluids.constants import R_inv
+from thermo import eos_mix
+import numpy as np
 
 __all__ = ('ActivityCoefficients',
            'IdealActivityCoefficients',
@@ -466,3 +471,73 @@ def loggammacs(self):
     def psi(self):
         return psi_modified_UNIFAC
 
+
+class GCEOSActivityCoefficients(ActivityCoefficients):
+    """
+    Abstract class for estimating fugacity coefficients using a generic cubic equation of state 
+    when called with composition, temperature (K), and pressure (Pa).
+    
+    Parameters
+    ----------
+    chemicals : Iterable[:class:`~thermosteam.Chemical`]
+    
+    """
+    __slots__ = ('_chemicals', '_eos')
+    EOS = None # type[GCEOSMIX] Subclasses must implement this attribute.
+    
+    def __init__(self, chemicals):
+        self.chemicals = chemicals
+    
+    @classmethod
+    def subclass(cls, EOS, name=None):
+        if name is None: name = EOS.__name__[:-3] + 'ActivityCoefficients'
+        return type(name, (cls,), dict(EOS=EOS))
+    
+    @property
+    def chemicals(self):
+        """tuple[Chemical] All chemicals involved in the calculation of fugacity coefficients."""
+        return self._chemicals
+    @chemicals.setter
+    def chemicals(self, chemicals):
+        self._chemicals = tuple(chemicals)
+    
+    def eos(self, zs, T, P):
+        if zs.__class__ is np.ndarray: zs = [float(i) for i in zs]
+        try:
+            self._eos = eos = self._eos.to_TP_zs_fast(T=T, P=P, zs=zs, only_l=True, full_alphas=False)
+        except:
+            data = tmo.ChemicalData(self.chemicals)
+            try:
+                kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', 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_l=True
+            )
+        return eos
+    
+    def __call__(self, x, T, P=101325):
+        eos = self.eos(x, T, P)
+        try:
+            log_phis = np.array(eos.dlnphi_dns(eos.Z_l)) + eos.G_dep_l * R_inv / T
+            return np.exp(log_phis)
+        except:
+            return 1.
+    f = __call__
+    args = ()
+        
+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 = GCEOSActivityCoefficients.subclass(getattr(eos_mix, name))
+    clsname = cls.__name__
+    clsnames.append(clsname)
+    dct[clsname] = cls
+
+__all__ = (*__all__, *clsnames)
+del dct, clsnames
\ No newline at end of file
diff --git a/thermosteam/equilibrium/fugacity_coefficients.py b/thermosteam/equilibrium/fugacity_coefficients.py
index 78ff54fb..9b7ea739 100644
--- a/thermosteam/equilibrium/fugacity_coefficients.py
+++ b/thermosteam/equilibrium/fugacity_coefficients.py
@@ -8,11 +8,16 @@
 """
 """
 from .ideal import ideal
-
+import thermosteam as tmo
+from thermo.interaction_parameters import IPDB
+from thermo import eos_mix
+from fluids.constants import R_inv
+import numpy as np
 
 __all__ = ('FugacityCoefficients', 
            'IdealFugacityCoefficients')
 
+
 class FugacityCoefficients:
     """
     Abstract class for the estimation of fugacity coefficients. Non-abstract subclasses should implement the following methods:
@@ -57,6 +62,70 @@ def chemicals(self, chemicals):
     def __call__(self, y, T, P):
         return 1.
 
+
+class GCEOSFugacityCoefficients(FugacityCoefficients):
+    """
+    Abstract class for estimating fugacity coefficients using a generic cubic equation of state 
+    when called with composition, temperature (K), and pressure (Pa).
+    
+    Parameters
+    ----------
+    chemicals : Iterable[:class:`~thermosteam.Chemical`]
+    
+    """
+    __slots__ = ('_chemicals', '_eos')
+    EOS = None # type[GCEOSMIX] Subclasses must implement this attribute.
     
+    @classmethod
+    def subclass(cls, EOS, name=None):
+        if name is None: name = EOS.__name__[:-3] + 'FugacityCoefficients'
+        return type(name, (cls,), dict(EOS=EOS))
     
+    @property
+    def chemicals(self):
+        """tuple[Chemical] All chemicals involved in the calculation of fugacity coefficients."""
+        return self._chemicals
+    @chemicals.setter
+    def chemicals(self, chemicals):
+        self._chemicals = tuple(chemicals)
+    
+    def eos(self, zs, T, P):
+        if zs.__class__ is np.ndarray: zs = [float(i) for i in zs]
+        try:
+            self._eos = eos = self._eos.to_TP_zs_fast(T=T, P=P, zs=zs, only_g=True, full_alphas=False)
+        except:
+            data = tmo.ChemicalData(self.chemicals)
+            try:
+                kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', 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=True
+            )
+        return eos
     
+    def __call__(self, x, T, P=101325):
+        eos = self.eos(x, T, P)
+        try:
+            log_phis = np.array(eos.dlnphi_dns(eos.Z_g)) + eos.G_dep_g * R_inv / T
+            return np.exp(log_phis)
+        except:
+            return 1.
+    f = __call__
+    args = ()
+
+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 = GCEOSFugacityCoefficients.subclass(getattr(eos_mix, name))
+    clsname = cls.__name__
+    clsnames.append(clsname)
+    dct[clsname] = cls
+
+__all__ = (*__all__, *clsnames)
+del dct, clsnames
diff --git a/thermosteam/equilibrium/vle.py b/thermosteam/equilibrium/vle.py
index 79a9fa99..c9574866 100644
--- a/thermosteam/equilibrium/vle.py
+++ b/thermosteam/equilibrium/vle.py
@@ -286,8 +286,8 @@ class VLE(Equilibrium, phases='lg'):
     H_hat_tol = 1e-6
     S_hat_tol = 1e-6
     V_tol = 1e-6
-    x_tol = 1e-9
-    y_tol = 1e-9
+    x_tol = 1e-8
+    y_tol = 1e-8
     default_method = 'fixed-point'
     
     def __init__(self, imol=None, thermal_condition=None,