diff --git a/tests/test_reaction.py b/tests/test_reaction.py
index 253f068b..07e549c1 100644
--- a/tests/test_reaction.py
+++ b/tests/test_reaction.py
@@ -223,7 +223,7 @@ def test_reactive_phase_equilibrium_no_kinetics():
     stream = tmo.Stream(
         H2O=1, Ethanol=5, LacticAcid=1
     )
-    stream.vle(T=360, P=101325, liquid_reaction=rxn)
+    stream.vle(T=360, P=101325, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['g'], 
         [0.0,
@@ -239,7 +239,7 @@ def test_reactive_phase_equilibrium_no_kinetics():
          1.8234109156732896]
     )
     
-    stream.vle(T=340, P=101325, liquid_reaction=rxn)
+    stream.vle(T=340, P=101325, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'], 
         stream.mol
@@ -248,7 +248,7 @@ def test_reactive_phase_equilibrium_no_kinetics():
         stream.imol['g'], 
         0
     )
-    stream.vle(T=450, P=101325, liquid_reaction=rxn)
+    stream.vle(T=450, P=101325, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['g'], 
         stream.mol
@@ -297,7 +297,7 @@ def rate(self, stream):
     )
     T = 360
     P = 101325
-    stream.vle(T=T, P=P, liquid_reaction=rxn)
+    stream.vle(T=T, P=P, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
         [0.026722998037919946,
@@ -321,7 +321,7 @@ def rate(self, stream):
     stream = tmo.Stream(
         H2O=2, Ethanol=5, LacticAcid=1, T=T,
     )
-    stream.vle(V=V, P=P, liquid_reaction=rxn)
+    stream.vle(V=V, P=P, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
         [0.026426291229780553,
@@ -343,7 +343,7 @@ def rate(self, stream):
     stream = tmo.Stream(
         H2O=2, Ethanol=5, LacticAcid=1, T=T,
     )
-    stream.vle(V=V, T=T, liquid_reaction=rxn)
+    stream.vle(V=V, T=T, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
         [0.026556271540719167,
@@ -362,7 +362,7 @@ def rate(self, stream):
     stream = tmo.Stream(
         H2O=2, Ethanol=5, LacticAcid=1, T=T,
     )
-    stream.vle(H=H, P=P, liquid_reaction=rxn)
+    stream.vle(H=H, P=P, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
         [0.026722998039327987,
diff --git a/thermosteam/_stream.py b/thermosteam/_stream.py
index b09ab135..e0fa2420 100644
--- a/thermosteam/_stream.py
+++ b/thermosteam/_stream.py
@@ -60,12 +60,39 @@ def __init__(self, stream, original, temporary):
         self.temporary = temporary
         
     def __enter__(self):
-        self.stream._phase._phase = self.temporary
+        stream = self.stream
+        stream._phase._phase = self.temporary
+        return stream
     
     def __exit__(self, type, exception, traceback):
         self.stream._phase._phase = self.original
         if exception: raise exception
    
+    
+class TemporaryStream:
+    __slots__ = ('stream', 'data', 'flow', 'T', 'P', 'phase')
+    
+    def __init__(self, stream, flow, T, P, phase):
+        self.stream = stream
+        self.data = stream.get_data()
+        self.flow = flow
+        self.T = T
+        self.P = P
+        self.phase = phase
+        
+    def __enter__(self):
+        stream = self.stream
+        self.data = stream.get_data()
+        if self.flow is not None: stream.imol.data[:] = self.flow
+        if self.T is not None: stream.T = self.T
+        if self.P is not None: stream.P = self.P
+        return stream
+    
+    def __exit__(self, type, exception, traceback):
+        self.stream.set_data(self.data)
+        if exception: raise exception
+   
+    
 class Equations:
     __slots__ = ('material', 'energy')
     
@@ -361,6 +388,9 @@ def __init__(self, ID: Optional[str]='',
             data = self._imol.data
             self.phases = [j for i, j in enumerate(self.phases) if data[i].any()]
 
+    def temporary(self, flow=None, T=None, P=None, phase=None):
+        return TemporaryStream(self, flow, T, P, phase)
+
     def temporary_phase(self, phase):
         return TemporaryPhase(self, self.phase, phase)
 
diff --git a/thermosteam/equilibrium/bubble_point.py b/thermosteam/equilibrium/bubble_point.py
index 1fb52a32..3dfb00a9 100644
--- a/thermosteam/equilibrium/bubble_point.py
+++ b/thermosteam/equilibrium/bubble_point.py
@@ -144,9 +144,9 @@ def _P_error(self, P, T, z_Psat_gamma, Psats, y):
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
         
-    def _T_error_reactive(self, T, P, z, dz, y, x, liquid_reaction):
+    def _T_error_reactive(self, T, P, z, dz, y, x, liquid_conversion):
         if T <= 0: raise InfeasibleRegion('negative temperature')
-        dz[:] = liquid_reaction.conversion(z, T, P, 'l')
+        dz[:] = liquid_conversion(z, T, P, 'l')
         x[:] = z + dz
         x /= x.sum()
         Psats = np.array([i(T) for i in self.Psats], dtype=float)
@@ -157,9 +157,9 @@ def _T_error_reactive(self, T, P, z, dz, y, x, liquid_reaction):
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
     
-    def _P_error_reactive(self, P, T, Psats, z, dz, y, x, liquid_reaction):
+    def _P_error_reactive(self, P, T, Psats, z, dz, y, x, liquid_conversion):
         if P <= 0: raise InfeasibleRegion('negative pressure')
-        dz[:] = liquid_reaction.conversion(z, T, P, 'l')
+        dz[:] = liquid_conversion(z, T, P, 'l')
         x[:] = z + dz
         x /= x.sum()
         z_Psat_gamma = x * Psats * self.gamma(x, T)
@@ -192,21 +192,21 @@ def _Py_ideal(self, z_Psat_gamma_pcf):
         y = z_Psat_gamma_pcf / P
         return P, y
     
-    def __call__(self, z, *, T=None, P=None, liquid_reaction=None):
+    def __call__(self, z, *, T=None, P=None, liquid_conversion=None):
         z = np.asarray(z, float)
         if T:
             if P: raise ValueError("may specify either T or P, not both")
-            P, *args = self.solve_Py(z, T, liquid_reaction)
+            P, *args = self.solve_Py(z, T, liquid_conversion)
         elif P:
-            T, *args = self.solve_Ty(z, P, liquid_reaction)
+            T, *args = self.solve_Ty(z, P, liquid_conversion)
         else:
             raise ValueError("must specify either T or P")
-        if liquid_reaction:
+        if liquid_conversion:
             return ReactiveBubblePointValues(T, P, self.IDs, z, *args)
         else:
             return BubblePointValues(T, P, self.IDs, z, *args)
     
-    def solve_Ty(self, z, P, liquid_reaction=None):
+    def solve_Ty(self, z, P, liquid_conversion=None):
         """
         Bubble point at given composition and pressure.
 
@@ -244,7 +244,7 @@ def solve_Ty(self, z, P, liquid_reaction=None):
             T = chemical.Tsat(P, check_validity=False) if P <= chemical.Pc else chemical.Tc
             y = z.copy()
             return T, fn.normalize(y)
-        elif liquid_reaction is None:
+        elif liquid_conversion is None:
             f = self._T_error
             z_norm = z / z.sum()
             z_over_P = z/P
@@ -268,7 +268,7 @@ def solve_Ty(self, z, P, liquid_reaction=None):
             dz = z_norm.copy()
             z_over_P = z / P
             T_guess, y = self._Ty_ideal(z_over_P)
-            args = (P, z_norm, dz, y, x, liquid_reaction)
+            args = (P, z_norm, dz, y, x, liquid_conversion)
             try:
                 T = flx.aitken_secant(f, T_guess, T_guess + 1e-3,
                                       self.T_tol, 5e-12, args,
@@ -281,7 +281,7 @@ def solve_Ty(self, z, P, liquid_reaction=None):
                                          checkiter=False, checkbounds=False)
             return T, dz, fn.normalize(y), x
     
-    def solve_Py(self, z, T, liquid_reaction=None):
+    def solve_Py(self, z, T, liquid_conversion=None):
         """
         Bubble point at given composition and temperature.
 
@@ -319,7 +319,7 @@ def solve_Py(self, z, T, liquid_reaction=None):
             P = chemical.Psat(T) if T <= chemical.Tc else chemical.Pc
             y = z.copy()
             return P, fn.normalize(y)
-        elif liquid_reaction is None:
+        elif liquid_conversion is None:
             if T > self.Tmax: T = self.Tmax
             elif T < self.Tmin: T = self.Tmin
             Psats = np.array([i(T) for i in self.Psats])
@@ -346,7 +346,7 @@ def solve_Py(self, z, T, liquid_reaction=None):
             dz = z_norm.copy()
             z_Psat_gamma = z * Psats * self.gamma(z_norm, T)
             P_guess, y = self._Py_ideal(z_Psat_gamma)
-            args = (T, Psats, z_norm, dz, y, x, liquid_reaction)
+            args = (T, Psats, z_norm, dz, y, x, liquid_conversion)
             try:
                 P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
                                       args, checkiter=False)
@@ -400,7 +400,7 @@ def __init__(self, chemicals=(), flasher=None):
     
     __call__ = BubblePoint.__call__
     
-    def solve_Ty(self, z, P, liquid_reaction=None):
+    def solve_Ty(self, z, P, liquid_conversion=None):
         """
         Bubble point at given composition and pressure.
 
@@ -442,7 +442,7 @@ def solve_Ty(self, z, P, liquid_reaction=None):
             T = results.T
         return T, fn.normalize(y)
     
-    def solve_Py(self, z, T, liquid_reaction=None):
+    def solve_Py(self, z, T, liquid_conversion=None):
         """
         Bubble point at given composition and temperature.
 
diff --git a/thermosteam/equilibrium/dew_point.py b/thermosteam/equilibrium/dew_point.py
index 5b9a20c5..40e53020 100644
--- a/thermosteam/equilibrium/dew_point.py
+++ b/thermosteam/equilibrium/dew_point.py
@@ -151,9 +151,9 @@ def _T_error(self, T, P, z_norm, zP, x):
         x[:] = self._solve_x(x_gamma, T, P, x)
         return 1 - x.sum()
     
-    def _T_error_reactive(self, T, P, z, dz, y, x, gas_reaction):
+    def _T_error_reactive(self, T, P, z, dz, y, x, gas_conversion):
         if T <= 0: raise InfeasibleRegion('negative temperature')
-        dz[:] = gas_reaction.conversion(z, T, P, 'g')
+        dz[:] = gas_conversion(z, T, P, 'g')
         y[:] = z + dz
         y /= y.sum()
         Psats = np.array([i(T) for i in self.Psats])
@@ -176,9 +176,9 @@ def _P_error(self, P, T, z_norm, z_over_Psats, Psats, x):
         x[:] = self._solve_x(x_gamma, T, P, x)
         return 1 - x.sum()
     
-    def _P_error_reactive(self, P, T, Psats, z, dz, y, x, gas_reaction):
+    def _P_error_reactive(self, P, T, Psats, z, dz, y, x, gas_conversion):
         if P <= 0: raise InfeasibleRegion('negative pressure')
-        dz[:] = gas_reaction.conversion(z, T, P, 'g')
+        dz[:] = gas_conversion(z, T, P, 'g')
         y[:] = z + dz
         y /= y.sum()
         x_gamma = y / Psats * P * self.phi(y, T, P) / self.pcf(T, P, Psats)
@@ -205,21 +205,21 @@ def _Px_ideal(self, z_over_Psats):
         x = z_over_Psats * P
         return P, x
     
-    def __call__(self, z, *, T=None, P=None, gas_reaction=None):
+    def __call__(self, z, *, T=None, P=None, gas_conversion=None):
         z = np.asarray(z, float)
         if T:
             if P: raise ValueError("may specify either T or P, not both")
-            P, *args = self.solve_Px(z, T, gas_reaction)
+            P, *args = self.solve_Px(z, T, gas_conversion)
         elif P:
-            T, *args = self.solve_Tx(z, P, gas_reaction)
+            T, *args = self.solve_Tx(z, P, gas_conversion)
         else:
             raise ValueError("must specify either T or P")
-        if gas_reaction:
+        if gas_conversion:
             return ReactiveDewPointValues(T, P, self.IDs, z, *args)
         else:
             return DewPointValues(T, P, self.IDs, z, *args)
     
-    def solve_Tx(self, z, P, gas_reaction=None):
+    def solve_Tx(self, z, P, gas_conversion=None):
         """
         Dew point given composition and pressure.
 
@@ -257,7 +257,7 @@ def solve_Tx(self, z, P, gas_reaction=None):
             T = chemical.Tsat(P, check_validity=False) if P <= chemical.Pc else chemical.Tc
             x = z.copy()
             return T, fn.normalize(x)
-        elif gas_reaction is None:
+        elif gas_conversion is None:
             f = self._T_error
             z_norm = z/z.sum()
             zP = z * P
@@ -282,7 +282,7 @@ def solve_Tx(self, z, P, gas_reaction=None):
             dz = z_norm.copy()
             zP = z * P
             T_guess, y = self._Tx_ideal(zP)
-            args = (P, z_norm, dz, y, x, gas_reaction)
+            args = (P, z_norm, dz, y, x, gas_conversion)
             try:
                 T = flx.aitken_secant(f, T_guess, T_guess + 1e-3,
                                       self.T_tol, 5e-12, args,
@@ -295,7 +295,7 @@ def solve_Tx(self, z, P, gas_reaction=None):
                                          checkiter=False, checkbounds=False)
             return T, dz, fn.normalize(y), x
     
-    def solve_Px(self, z, T, gas_reaction=None):
+    def solve_Px(self, z, T, gas_conversion=None):
         """
         Dew point given composition and temperature.
 
@@ -333,7 +333,7 @@ def solve_Px(self, z, T, gas_reaction=None):
             P = chemical.Psat(T) if T <= chemical.Tc else chemical.Pc
             x = z.copy()
             return P, fn.normalize(x)
-        elif gas_reaction is None:
+        elif gas_conversion is None:
             z_norm = z / z.sum()
             Psats = np.array([i(T) for i in self.Psats], dtype=float)
             z_over_Psats = z / Psats
@@ -359,7 +359,7 @@ def solve_Px(self, z, T, gas_reaction=None):
             Psats = np.array([i(T) for i in self.Psats], dtype=float)
             z_over_Psats = z / Psats
             P_guess, x = self._Px_ideal(z_over_Psats)
-            args = (T, Psats, z_norm, dz, y, x, gas_reaction)
+            args = (T, Psats, z_norm, dz, y, x, gas_conversion)
             try:
                 P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
                                       args, checkiter=False)
diff --git a/thermosteam/equilibrium/vle.py b/thermosteam/equilibrium/vle.py
index 2b56ea73..af1b1b5d 100644
--- a/thermosteam/equilibrium/vle.py
+++ b/thermosteam/equilibrium/vle.py
@@ -22,6 +22,7 @@
 from . import activity_coefficients as ac
 from .. import functional as fn
 from ..utils import Cache
+import thermosteam as tmo
 import numpy as np
 
 __all__ = ('VLE', 'VLECache')
@@ -41,16 +42,16 @@ def update_xV(xV, V, Ks, z):
     xV[:-1] = z/(1. + V * (Ks - 1.))
 
 def xV_iter_2n(xVlogK, pcf_Psat_over_P, T, P, z, f_gamma, gamma_args, f_phi,
-               n, gas_reaction, liquid_reaction, index):
+               n, gas_conversion, liquid_conversion, index):
     xVlogK = xVlogK.copy()
     xV = xVlogK[:-n]
     Ks = np.exp(xVlogK[-n:])
     x, y = xy(xV, Ks)
     Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P)
-    if gas_reaction: 
-        z = z + gas_reaction.conversion(material=y * xV[n], T=T, P=P, phase='g')
-    if liquid_reaction:
-        z = z + liquid_reaction.conversion(material=x * (1 - xV[n]), T=T, P=P, phase='l')
+    if gas_conversion: 
+        z = z + gas_conversion(material=y * xV[n], T=T, P=P, phase='g')
+    if liquid_conversion:
+        z = z + liquid_conversion(material=x * (1 - xV[n]), T=T, P=P, phase='l')
     V = binary.compute_phase_fraction_2N(z, Ks)
     update_xV(xV, V, Ks, z)
     xVlogK[-n:] = np.log(Ks)
@@ -59,17 +60,17 @@ def xV_iter_2n(xVlogK, pcf_Psat_over_P, T, P, z, f_gamma, gamma_args, f_phi,
 def xV_iter(
         xVlogK, pcf_Psat_over_P, T, P,
         z, z_light, z_heavy, f_gamma, gamma_args,
-        f_phi, n, gas_reaction, liquid_reaction, index
+        f_phi, n, gas_conversion, liquid_conversion, index
     ):
     xVlogK = xVlogK.copy()
     xV = xVlogK[:-n]
     V = xV[n]
     Ks = np.exp(xVlogK[-n:])
     x, y = xy(xV, Ks)
-    if gas_reaction: 
-        z = z + gas_reaction.conversion(material=y * V, T=T, P=P, phase='g')
-    if liquid_reaction:
-        z = z + liquid_reaction.conversion(material=x * (1 - V), T=T, P=P, phase='l')
+    if gas_conversion: 
+        z = z + gas_conversion(material=y * V, T=T, P=P, phase='g')
+    if liquid_conversion:
+        z = z + liquid_conversion(material=x * (1 - V), T=T, P=P, phase='l')
     Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P)
     V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, xV[-1], z_light, z_heavy)
     update_xV(xV, V, Ks, z)
@@ -304,7 +305,7 @@ def __init__(self, imol=None, thermal_condition=None,
         self._index = ()
     
     def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None,
-                 gas_reaction=None, liquid_reaction=None):
+                 gas_conversion=None, liquid_conversion=None):
         """
         Perform vapor-liquid equilibrium.
 
@@ -332,7 +333,17 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None,
         P or T (e.g., x and V is invalid).
         
         """
-        if gas_reaction or liquid_reaction:
+        if gas_conversion or liquid_conversion:
+            if gas_conversion:
+                if isinstance(gas_conversion, tmo.KineticReaction):
+                    gas_conversion = gas_conversion.conversion_handle(tmo.Stream(thermo=self.thermo))
+                elif isinstance(gas_conversion, (tmo.Reaction, tmo.ReactionItem)):
+                    gas_conversion = gas_conversion.conversion_handle(None)
+            if liquid_conversion:
+                if isinstance(liquid_conversion, tmo.KineticReaction):
+                    liquid_conversion = liquid_conversion.conversion_handle(tmo.Stream(thermo=self.thermo))
+                elif isinstance(liquid_conversion, (tmo.Reaction, tmo.ReactionItem)):
+                    liquid_conversion = liquid_conversion.conversion_handle(None)
             if x is not None or y is not None:
                 raise ValueError(
                     "can not pass either 'x' or 'y' arguments with reactions preset"
@@ -354,21 +365,21 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None,
         if T_spec:
             if P_spec:
                 try:
-                    self.set_thermal_condition(T, P, gas_reaction, liquid_reaction)
+                    self.set_thermal_condition(T, P, gas_conversion, liquid_conversion)
                 except NoEquilibrium:
                     thermal_condition = self._thermal_condition
                     thermal_condition.T = T
                     thermal_condition.P = P
             elif V_spec:
                 try:
-                    self.set_TV(T, V, gas_reaction, liquid_reaction)
+                    self.set_TV(T, V, gas_conversion, liquid_conversion)
                 except NoEquilibrium:
                     thermal_condition = self._thermal_condition
                     thermal_condition.T = T
             elif H_spec:
-                self.set_TH(T, H, gas_reaction, liquid_reaction)
+                self.set_TH(T, H, gas_conversion, liquid_conversion)
             elif S_spec:
-                self.set_TS(T, S, gas_reaction, liquid_reaction)
+                self.set_TS(T, S, gas_conversion, liquid_conversion)
             elif x_spec:
                 self.set_Tx(T, np.asarray(x))
             else: # y_spec
@@ -376,22 +387,22 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None,
         elif P_spec:
             if V_spec:
                 try:
-                    self.set_PV(P, V, gas_reaction, liquid_reaction)
+                    self.set_PV(P, V, gas_conversion, liquid_conversion)
                 except NoEquilibrium:
                     thermal_condition = self._thermal_condition
                     thermal_condition.P = P
             elif H_spec:
                 try:
-                    self.set_PH(P, H, gas_reaction, liquid_reaction)
+                    self.set_PH(P, H, gas_conversion, liquid_conversion)
                 except NoEquilibrium:
                     thermal_condition = self._thermal_condition
                     thermal_condition.P = P
             elif S_spec:
                 try:
-                    self.set_PS(P, S, gas_reaction, liquid_reaction)
+                    self.set_PS(P, S, gas_conversion, liquid_conversion)
                 except:
                     try:
-                        self.set_PS(P, S, gas_reaction, liquid_reaction)
+                        self.set_PS(P, S, gas_conversion, liquid_conversion)
                     except NoEquilibrium:
                         thermal_condition = self._thermal_condition
                         thermal_condition.P = P
@@ -631,16 +642,16 @@ def set_Py(self, P, y):
         self._thermal_condition.T, x = self._dew_point.solve_Tx(y, P) 
         self._lever_rule(x, y)
         
-    def set_thermal_condition(self, T, P, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_thermal_condition(self, T, P, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         thermal_condition = self._thermal_condition
         self._T = thermal_condition.T = T
         self._P = thermal_condition.P = P
         if self._N == 0: return
         if self._N == 1: return self._set_thermal_condition_chemical(T, P)
         # Check if there is equilibrium
-        if gas_reaction:
-            P_dew, dz_dew, x_dew, _ = self._dew_point.solve_Px(self._z, T, gas_reaction)
+        if gas_conversion:
+            P_dew, dz_dew, x_dew, _ = self._dew_point.solve_Px(self._z, T, gas_conversion)
             if P <= P_dew and not self._F_mol_heavy:
                 self._vapor_mol[self._index] = self._mol_vle + dz_dew * self._F_mol_vle
                 self._liquid_mol[self._index] = 0
@@ -652,8 +663,8 @@ def set_thermal_condition(self, T, P, gas_reaction=None, liquid_reaction=None):
                 self._liquid_mol[self._index] = 0
                 return
             dz_dew = None
-        if liquid_reaction:
-            P_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Py(self._z, T, liquid_reaction)
+        if liquid_conversion:
+            P_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Py(self._z, T, liquid_conversion)
             if P >= P_bubble and not self._F_mol_light:
                 self._vapor_mol[self._index] = 0
                 self._liquid_mol[self._index] = self._mol_vle + dz_bubble * self._F_mol_vle
@@ -670,12 +681,12 @@ def set_thermal_condition(self, T, P, gas_reaction=None, liquid_reaction=None):
         dP = (P_bubble - P_dew)
         V = (P - P_dew) / dP if dP > 1. else 0.5
         self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew)
-        v = self._solve_v(T, P, gas_reaction, liquid_reaction)
-        mol_vle = self._mol_vle + self._dmol_vle if (gas_reaction or liquid_reaction) else self._mol_vle 
+        v = self._solve_v(T, P, gas_conversion, liquid_conversion)
+        mol_vle = self._mol_vle + self._dmol_vle if (gas_conversion or liquid_conversion) else self._mol_vle 
         set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle)
         
-    def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_TV(self, T, V, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         thermal_condition = self._thermal_condition
         thermal_condition.T = self._T = T
         if self._N == 0: raise RuntimeError('no chemicals present to perform VLE')
@@ -683,8 +694,8 @@ def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
         if self._F_mol_heavy and V == 1.: V = 1. - 1e-3
         if self._F_mol_light and V == 0.: V = 1e-3
         if V == 1:
-            if gas_reaction:
-                P_dew, dz, y, x_dew = self._dew_point.solve_Px(self._z, T, gas_reaction)
+            if gas_conversion:
+                P_dew, dz, y, x_dew = self._dew_point.solve_Px(self._z, T, gas_conversion)
                 self._vapor_mol[self._index] = self._mol_vle + dz * self._F_mol_vle
             else:
                 P_dew, x_dew = self._dew_point.solve_Px(self._z, T)
@@ -692,7 +703,7 @@ def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
             self._liquid_mol[self._index] = 0
             thermal_condition.P = P_dew
         elif V == 0:
-            if liquid_reaction:
+            if liquid_conversion:
                 P_bubble, dz, y_bubble, x = self._bubble_point.solve_Py(self._z, T)
                 self._liquid_mol[self._index] = self._mol_vle + dz * self._F_mol_vle
             else:
@@ -701,25 +712,25 @@ def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
             self._vapor_mol[self._index] = 0
             thermal_condition.P = P_bubble
         else:
-            if liquid_reaction:
-                P_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Py(self._z, T, liquid_reaction)
+            if liquid_conversion:
+                P_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Py(self._z, T, liquid_conversion)
             else:
                 P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T)
                 dz_bubble = None
-            if gas_reaction:
-                P_dew, dz_dew, y, x_dew = self._dew_point.solve_Px(self._z, T, gas_reaction)
+            if gas_conversion:
+                P_dew, dz_dew, y, x_dew = self._dew_point.solve_Px(self._z, T, gas_conversion)
             else:
-                P_dew, x_dew = self._dew_point.solve_Px(self._z, T, gas_reaction)
+                P_dew, x_dew = self._dew_point.solve_Px(self._z, T, gas_conversion)
                 dz_dew = None
             self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew)
             if self._F_mol_light: P_bubble = 0.1 * self._bubble_point.Pmax + 0.9 * P_bubble
             if self._F_mol_heavy: P_dew = 0.1 * self._bubble_point.Pmin + 0.9 * P_dew
             
-            V_bubble = self._V_err_at_P(P_bubble, 0., gas_reaction, liquid_reaction)
+            V_bubble = self._V_err_at_P(P_bubble, 0., gas_conversion, liquid_conversion)
             if V_bubble > V:
                 F_mol = self._F_mol
                 mol = self._mol_vle
-                if liquid_reaction:
+                if liquid_conversion:
                     mol = mol + self._dmol_vle
                     F_mol = F_mol + self._dF_mol
                 F_mol_vapor = F_mol * V
@@ -728,11 +739,11 @@ def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
                 v[mask] = mol[mask]
                 P = P_bubble
             else:
-                V_dew = self._V_err_at_P(P_dew, 0., gas_reaction, liquid_reaction)
+                V_dew = self._V_err_at_P(P_dew, 0., gas_conversion, liquid_conversion)
                 if V_dew < V:
                     F_mol = self._F_mol
                     mol = self._mol_vle
-                    if gas_reaction:
+                    if gas_conversion:
                         mol = mol + self._dmol_vle
                         F_mol = F_mol + self._dF_mol
                     l = x_dew * F_mol * (1. - V)
@@ -745,18 +756,18 @@ def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
                         self._V_err_at_P,
                         P_bubble, P_dew, V_bubble - V, V_dew - V,
                         self._P, self.P_tol, self.V_tol,
-                        (V, gas_reaction, liquid_reaction),
+                        (V, gas_conversion, liquid_conversion),
                         checkiter=False, checkbounds=False,
                         maxiter=self.maxiter,
                     )
                     v = self._v
                     mol = self._mol_vle
-                    if gas_reaction or liquid_reaction:
+                    if gas_conversion or liquid_conversion:
                         mol = mol + self._dmol_vle
             
             self._P = thermal_condition.P = P
             set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol)
-            if liquid_reaction or gas_reaction:
+            if liquid_conversion or gas_conversion:
                 try: self._H_hat = (
                         self.mixture.xH(self._phase_data, T, P) + (self.chemicals.Hf * mol).sum()
                     ) / self._F_mass
@@ -765,8 +776,8 @@ def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None):
                 try: self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass
                 except: pass
 
-    def set_TH(self, T, H, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_TH(self, T, H, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         if self._N == 0: raise RuntimeError('no chemicals present to perform VLE')
         if self._N == 1: return self._set_TH_chemical(T, H)
         self._T = T
@@ -814,8 +825,8 @@ def set_TH(self, T, H, gas_reaction=None, liquid_reaction=None):
         self._P = self._thermal_condition.P = P   
         self._thermal_condition.T = T
     
-    def set_TS(self, T, S, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_TS(self, T, S, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         if self._N == 0: raise RuntimeError('no chemicals present to perform VLE')
         if self._N == 1: return self._set_TS_chemical(T, S)
         self._T = T
@@ -863,8 +874,8 @@ def set_TS(self, T, S, gas_reaction=None, liquid_reaction=None):
         self._P = self._thermal_condition.P = P   
         self._thermal_condition.T = T
     
-    def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_PV(self, P, V, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         self._thermal_condition.P = self._P = P
         if self._N == 0: raise RuntimeError('no chemicals present to perform VLE')
         if self._N == 1: return self._set_PV_chemical(P, V)
@@ -877,8 +888,8 @@ def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
         if self._F_mol_heavy and V == 1.: V = 1. - 1e-3
         if self._F_mol_light and V == 0.: V = 1e-3
         if V == 1 and not self._F_mol_heavy:
-            if gas_reaction:
-                T_dew, dz, y, x_dew = self._dew_point.solve_Tx(self._z, P, gas_reaction)
+            if gas_conversion:
+                T_dew, dz, y, x_dew = self._dew_point.solve_Tx(self._z, P, gas_conversion)
                 self._vapor_mol[self._index] = self._mol_vle + dz * self._F_mol_vle
             else:
                 T_dew, x_dew = self._dew_point.solve_Tx(self._z, P)
@@ -886,8 +897,8 @@ def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
             self._liquid_mol[self._index] = 0
             thermal_condition.T = T_dew
         elif V == 0 and not self._F_mol_light:
-            if liquid_reaction:
-                T_bubble, dz, y_bubble, x = self._bubble_point.solve_Ty(self._z, P, liquid_reaction)
+            if liquid_conversion:
+                T_bubble, dz, y_bubble, x = self._bubble_point.solve_Ty(self._z, P, liquid_conversion)
                 self._liquid_mol[self._index] = self._mol_vle + dz * self._F_mol_vle
             else:
                 T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P)
@@ -895,24 +906,24 @@ def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
             self._vapor_mol[self._index] = 0
             thermal_condition.T = T_bubble
         else:
-            if liquid_reaction:
-                T_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Ty(self._z, P, liquid_reaction)
+            if liquid_conversion:
+                T_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Ty(self._z, P, liquid_conversion)
             else:
                 T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P)
                 dz_bubble = None
-            if gas_reaction:
-                T_dew, dz_dew, y, x_dew = self._dew_point.solve_Tx(self._z, P, gas_reaction)
+            if gas_conversion:
+                T_dew, dz_dew, y, x_dew = self._dew_point.solve_Tx(self._z, P, gas_conversion)
             else:
                 T_dew, x_dew = self._dew_point.solve_Tx(self._z, P)
                 dz_dew = None
             self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew)
             if self._F_mol_light: T_bubble = 0.1 * self._bubble_point.Tmin + 0.9 * T_bubble
             if self._F_mol_heavy: T_dew = 0.1 * self._bubble_point.Tmax + 0.9 * T_dew
-            V_bubble = self._V_err_at_T(T_bubble, 0., gas_reaction, liquid_reaction)
+            V_bubble = self._V_err_at_T(T_bubble, 0., gas_conversion, liquid_conversion)
             if V_bubble > V:
                 F_mol = self._F_mol
                 mol = self._mol_vle
-                if liquid_reaction:
+                if liquid_conversion:
                     mol = mol + self._dmol_vle
                     F_mol = F_mol + self._dF_mol
                 F_mol_vapor = F_mol * V
@@ -921,11 +932,11 @@ def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
                 v[mask] = mol[mask]
                 T = T_bubble
             else:
-                V_dew = self._V_err_at_T(T_dew, 0., gas_reaction, liquid_reaction)
+                V_dew = self._V_err_at_T(T_dew, 0., gas_conversion, liquid_conversion)
                 if V_dew < V:
                     F_mol = self._F_mol
                     mol = self._mol_vle
-                    if gas_reaction:
+                    if gas_conversion:
                         mol = mol + self._dmol_vle
                         F_mol = F_mol + self._dF_mol
                     l = x_dew * F_mol * (1. - V)
@@ -938,17 +949,17 @@ def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
                         self._V_err_at_T,
                         T_bubble, T_dew, V_bubble - V, V_dew - V,
                         self._T, self.T_tol, self.V_tol,
-                        (V, gas_reaction, liquid_reaction),
+                        (V, gas_conversion, liquid_conversion),
                         checkiter=False, checkbounds=False,
                         maxiter=self.maxiter,
                     )
                     v = self._v
                     mol = self._mol_vle
-                    if gas_reaction or liquid_reaction:
+                    if gas_conversion or liquid_conversion:
                         mol = mol + self._dmol_vle
             self._T = thermal_condition.T = T
             set_flows(vapor_mol, liquid_mol, index, v, mol)
-            if liquid_reaction or gas_reaction:
+            if liquid_conversion or gas_conversion:
                 try: self._H_hat = (
                         self.mixture.xH(self._phase_data, T, P) + (self.chemicals.Hf * mol).sum(0)
                     ) / self._F_mass
@@ -957,8 +968,8 @@ def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None):
                 try: self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass
                 except: pass
     
-    def set_PS(self, P, S, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_PS(self, P, S, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         thermal_condition = self._thermal_condition
         thermal_condition.P = self._P = P
         if self._N == 0: 
@@ -1071,8 +1082,8 @@ def set_PS(self, P, S, gas_reaction=None, liquid_reaction=None):
             )
         self._S_hat = S_hat
     
-    def set_PH(self, P, H, gas_reaction=None, liquid_reaction=None):
-        self._setup(gas_reaction, liquid_reaction)
+    def set_PH(self, P, H, gas_conversion=None, liquid_conversion=None):
+        self._setup(gas_conversion, liquid_conversion)
         thermal_condition = self._thermal_condition
         thermal_condition.P = self._P = P
         if self._N == 0: 
@@ -1088,8 +1099,8 @@ def set_PH(self, P, H, gas_reaction=None, liquid_reaction=None):
         liquid_mol = self._liquid_mol
         
         # Check if subcooled liquid
-        if liquid_reaction:
-            T_bubble, dz_bubble, y_bubble, x = self._bubble_point.solve_Ty(self._z, P, liquid_reaction)
+        if liquid_conversion:
+            T_bubble, dz_bubble, y_bubble, x = self._bubble_point.solve_Ty(self._z, P, liquid_conversion)
             mol = self._mol_vle + dz_bubble * self._F_mol_vle
         else:
             T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P)
@@ -1099,18 +1110,18 @@ def set_PH(self, P, H, gas_reaction=None, liquid_reaction=None):
         vapor_mol[index] = 0
         liquid_mol[index] = mol
         H_bubble = self.mixture.xH(self._phase_data, T_bubble, P)
-        if liquid_reaction: 
+        if liquid_conversion: 
             Hf = (self.chemicals.Hf * mol).sum()
             H_bubble += Hf
         dH_bubble = H - H_bubble
         if dH_bubble <= 0:
-            if liquid_reaction: H -= Hf
+            if liquid_conversion: H -= Hf
             thermal_condition.T = self.mixture.xsolve_T_at_HP(self._phase_data, H, T_bubble, P)
             return
         
         # Check if super heated vapor
-        if gas_reaction:
-            T_dew, dz_dew, y, x_dew = self.dew_point.solve_Tx(self._z, P, gas_reaction)
+        if gas_conversion:
+            T_dew, dz_dew, y, x_dew = self.dew_point.solve_Tx(self._z, P, gas_conversion)
             mol = self._mol_vle + dz_dew
         else:
             T_dew, x_dew = self._dew_point.solve_Tx(self._z, P)
@@ -1124,12 +1135,12 @@ def set_PH(self, P, H, gas_reaction=None, liquid_reaction=None):
         vapor_mol[index] = mol
         liquid_mol[index] = 0
         H_dew = self.mixture.xH(self._phase_data, T_dew, P)
-        if gas_reaction: 
+        if gas_conversion: 
             Hf = (self.chemicals.Hf * mol).sum()
             H_dew += Hf
         dH_dew = H - H_dew
         if dH_dew >= 0:
-            if gas_reaction: H -= Hf
+            if gas_conversion: H -= Hf
             thermal_condition.T = self.mixture.xsolve_T_at_HP(self._phase_data, H, T_dew, P)
             return
         
@@ -1140,21 +1151,21 @@ def set_PH(self, P, H, gas_reaction=None, liquid_reaction=None):
         F_mass = self._F_mass
         H_hat = H/F_mass
         
-        if (H_hat_bubble:=self._H_hat_err_at_T(T_bubble, 0., gas_reaction, liquid_reaction)) > H_hat:
+        if (H_hat_bubble:=self._H_hat_err_at_T(T_bubble, 0., gas_conversion, liquid_conversion)) > H_hat:
             T = T_bubble
-        elif (H_hat_dew:=self._H_hat_err_at_T(T_dew, 0., gas_reaction, liquid_reaction)) < H_hat:
+        elif (H_hat_dew:=self._H_hat_err_at_T(T_dew, 0., gas_conversion, liquid_conversion)) < H_hat:
             T = T_dew
         else:
             T = flx.IQ_interpolation(
                 self._H_hat_err_at_T, T_bubble, T_dew, 
                 H_hat_bubble - H_hat, H_hat_dew - H_hat,
                 self._T, self.T_tol, self.H_hat_tol,
-                (H_hat, gas_reaction, liquid_reaction),
+                (H_hat, gas_conversion, liquid_conversion),
                 checkiter=False, checkbounds=False,
                 maxiter=self.maxiter
             )
          
-        if gas_reaction or liquid_reaction: 
+        if gas_conversion or liquid_conversion: 
             H -= (self.chemicals.Hf * (self._mol_vle + self._dmol_vle)).sum()
             
         # Make sure energy balance is correct by vaporizing a fraction
@@ -1235,9 +1246,9 @@ def _refresh_v(self, V, y_bubble, x_dew, dz_bubble=None, dz_dew=None): # TODO: u
             l[l < 0] = 1e-12
             self._x = fn.normalize(l, l.sum() + self._F_mol_heavy)
     
-    def _H_hat_err_at_T(self, T, H_hat, gas_reaction, liquid_reaction):
-        v = self._solve_v(T, self._P, gas_reaction, liquid_reaction)
-        if gas_reaction or liquid_reaction: 
+    def _H_hat_err_at_T(self, T, H_hat, gas_conversion, liquid_conversion):
+        v = self._solve_v(T, self._P, gas_conversion, liquid_conversion)
+        if gas_conversion or liquid_conversion: 
             mol_vle = self._mol_vle + self._dmol_vle
             set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle)
             H = self.mixture.xH(self._phase_data, T, self._P) + (self.chemicals.Hf * mol_vle).sum()
@@ -1248,9 +1259,9 @@ def _H_hat_err_at_T(self, T, H_hat, gas_reaction, liquid_reaction):
         self._H_hat = H / self._F_mass
         return self._H_hat - H_hat
     
-    def _H_hat_err_at_P(self, P, H_hat, gas_reaction=None, liquid_reaction=None):
-        v = self._solve_v(self._T, P, gas_reaction, liquid_reaction)
-        if gas_reaction or liquid_reaction: 
+    def _H_hat_err_at_P(self, P, H_hat, gas_conversion=None, liquid_conversion=None):
+        v = self._solve_v(self._T, P, gas_conversion, liquid_conversion)
+        if gas_conversion or liquid_conversion: 
             mol_vle = self._mol_vle + self._dmol_vle
             set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle)
             H = self.mixture.xH(self._phase_data, self._T, P) + (self.chemicals.Hf * mol_vle).sum()
@@ -1273,27 +1284,27 @@ def _S_hat_err_at_P(self, P, S_hat):
         self._S_hat = self.mixture.xS(self._phase_data, self._T, P)/self._F_mass
         return self._S_hat - S_hat
     
-    def _V_err_at_P(self, P, V, gas_reaction, liquid_reaction):
-        v = self._solve_v(self._T , P, gas_reaction, liquid_reaction).sum()
-        if gas_reaction or liquid_reaction:
+    def _V_err_at_P(self, P, V, gas_conversion, liquid_conversion):
+        v = self._solve_v(self._T , P, gas_conversion, liquid_conversion).sum()
+        if gas_conversion or liquid_conversion:
             F_mol_vle = self._F_mol_vle + self._dF_mol
         else:
             F_mol_vle = self._F_mol_vle
         return v / F_mol_vle - V
         
-    def _V_err_at_T(self, T, V, gas_reaction, liquid_reaction):
-        v = self._solve_v(T, self._P, gas_reaction, liquid_reaction).sum()
-        if gas_reaction or liquid_reaction:
+    def _V_err_at_T(self, T, V, gas_conversion, liquid_conversion):
+        v = self._solve_v(T, self._P, gas_conversion, liquid_conversion).sum()
+        if gas_conversion or liquid_conversion:
             F_mol_vle = self._F_mol_vle + self._dF_mol
         else:
             F_mol_vle = self._F_mol_vle
         return v / F_mol_vle - V
     
-    def _solve_v(self, T, P, gas_reaction=None, liquid_reaction=None):
+    def _solve_v(self, T, P, gas_conversion=None, liquid_conversion=None):
         """Solve for vapor mol"""
         method = self.method
         if method == 'shgo':
-            if gas_reaction or liquid_reaction:
+            if gas_conversion or liquid_conversion:
                 raise RuntimeError('shgo is not a valid method (yet) when reactions are present')
             gamma = self._gamma
             phi = self._phi
@@ -1312,8 +1323,8 @@ def _solve_v(self, T, P, gas_reaction=None, liquid_reaction=None):
                               self._bubble_point.Psats])
             pcf_Psats_over_P = self._pcf(T, P, Psats) * Psats / P
             self._T = T
-            y = self._solve_y(self._y, pcf_Psats_over_P, T, P, gas_reaction, liquid_reaction)
-            if gas_reaction or liquid_reaction:
+            y = self._solve_y(self._y, pcf_Psats_over_P, T, P, gas_conversion, liquid_conversion)
+            if gas_conversion or liquid_conversion:
                 self._v = v = (self._F_mol + self._dF_mol) * self._V * y
                 mol_vle = self._mol_vle + self._dmol_vle
             else:
@@ -1326,7 +1337,7 @@ def _solve_v(self, T, P, gas_reaction=None, liquid_reaction=None):
             raise RuntimeError(f"invalid method '{method}'")
         return v
     
-    def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_reaction):
+    def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_conversion, liquid_conversion):
         gamma = self._gamma
         x = self._x
         pcf_Psat_over_P = pcf_Psat_over_P
@@ -1339,12 +1350,12 @@ def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_reaction):
             args = (pcf_Psat_over_P, T, P, z, 
                     self._z_light, self._z_heavy, 
                     gamma.f, gamma.args, self._phi, n,
-                    gas_reaction, liquid_reaction, index)
+                    gas_conversion, liquid_conversion, index)
         else:
             f = xV_iter_2n
             args = (pcf_Psat_over_P, T, P, z, 
                     gamma.f, gamma.args, self._phi, n,
-                    gas_reaction, liquid_reaction, index)
+                    gas_conversion, liquid_conversion, index)
         xVlogK = np.zeros(2*x.size + 1)
         xVlogK[:n] = x
         xVlogK[n] = self._V
@@ -1363,12 +1374,12 @@ def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_reaction):
         else:
             Ks = (z / x - 1) / V + 1.
         self._z_last = z
-        if gas_reaction or liquid_reaction:
+        if gas_conversion or liquid_conversion:
             dmol_vle = 0
-            if gas_reaction: 
-                dmol_vle += gas_reaction.conversion(material=y * V, T=T, P=P, phase='g')
-            if liquid_reaction:
-                dmol_vle += liquid_reaction.conversion(material=x * (1 - V), T=T, P=P, phase='l')
+            if gas_conversion: 
+                dmol_vle += gas_conversion(material=y * V, T=T, P=P, phase='g')
+            if liquid_conversion:
+                dmol_vle += liquid_conversion(material=x * (1 - V), T=T, P=P, phase='l')
             dmol_vle *= self._F_mol
             self._dmol_vle = dmol_vle
             self._dF_mol = dmol_vle.sum()
@@ -1377,7 +1388,7 @@ def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_reaction):
             v = self._F_mol * V * x * Ks   
         return fn.normalize(v, v.sum() + self._F_mol_light)
     
-    def _setup(self, gas_reaction=None, liquid_reaction=None):
+    def _setup(self, gas_conversion=None, liquid_conversion=None):
         imol = self._imol
         self._phase_data = tuple(imol)
         self._liquid_mol = liquid_mol = imol['l']
@@ -1385,17 +1396,19 @@ def _setup(self, gas_reaction=None, liquid_reaction=None):
         mol = liquid_mol + vapor_mol
         if not mol.any(): raise NoEquilibrium('no chemicals to perform equilibrium')
         nonzero = set(mol.nonzero_keys())
-        if gas_reaction:
+        if gas_conversion:
             nonzero.update(
-                gas_reaction.stoichiometry.nonzero_keys()
+                gas_conversion.reaction.stoichiometry.nonzero_keys()
             )
-        if liquid_reaction:
+        if liquid_conversion:
             nonzero.update(
-                liquid_reaction.stoichiometry.nonzero_keys()
+                liquid_conversion.reaction.stoichiometry.nonzero_keys()
             )
         chemicals = self.chemicals
         if self._nonzero == nonzero:
             index = self._index
+            eq_chems = chemicals.tuple
+            eq_chems = [eq_chems[i] for i in index]
             reset = False
         else:
             # Set up indices for both equilibrium and non-equilibrium species
@@ -1405,10 +1418,10 @@ def _setup(self, gas_reaction=None, liquid_reaction=None):
             reset = True     
             self._nonzero = nonzero
             self._index = index
-            if gas_reaction and [i.ID for i in eq_chems] != [i.ID for i in gas_reaction.chemicals]:
-                gas_reaction.reset_thermo(gas_reaction.thermo.subset(eq_chems))
-            if liquid_reaction and [i.ID for i in eq_chems] != [i.ID for i in liquid_reaction.chemicals]:
-                liquid_reaction.reset_thermo(liquid_reaction.thermo.subset(eq_chems))
+        if gas_conversion and [i.ID for i in eq_chems] != [i.ID for i in gas_conversion.reaction.chemicals]:
+            gas_conversion.set_chemicals(eq_chems)
+        if liquid_conversion and [i.ID for i in eq_chems] != [i.ID for i in liquid_conversion.reaction.chemicals]:
+            liquid_conversion.set_chemicals(eq_chems)
         
         # Get overall composition
         self._F_mass = (chemicals.MW * mol).sum()
diff --git a/thermosteam/reaction/_reaction.py b/thermosteam/reaction/_reaction.py
index 346b9e14..dcf4ae4d 100644
--- a/thermosteam/reaction/_reaction.py
+++ b/thermosteam/reaction/_reaction.py
@@ -28,7 +28,7 @@
     'Reaction', 'ReactionItem', 'ReactionSet', 
     'ParallelReaction', 'SeriesReaction', 'ReactionSystem',
     'Rxn', 'RxnS', 'RxnI', 'PRxn', 'SRxn', 'RxnSys',
-    'KineticReaction',
+    'KineticReaction', 'KineticConversion',
     'KRxn',
 )
 
@@ -102,6 +102,25 @@ def as_material_array(material, basis, phases, chemicals):
 
 # %% Extent of reaction (stoichiometric)
 
+class Conversion:
+    __slots__ = ('reaction', 'index')
+    
+    def __init__(self, reaction):
+        self.reaction = reaction
+        self.index = None
+        
+    def set_chemicals(self, chemicals):
+        self.index = self.reaction.chemicals.indices([i.ID for i in chemicals])
+        
+    def __call__(self, material=None, T=None, P=None, phase=None):
+        if self.index is None:
+            return self.reaction.conversion(material)
+        else:
+            flow = SparseVector.from_size(self.reaction.chemicals.size)
+            flow[self.index] = material
+            return self.reaction.conversion(material)
+
+
 @chemicals_user
 class Reaction:
     """
@@ -121,7 +140,7 @@ class Reaction:
         Reactant conversion (fraction). Defaults to 1.
     chemicals : Chemicals, optional
         Chemicals corresponding to each entry in the stoichiometry array. 
-        Defaults to settings.chemicals
+        Defaults to settings.chemicals.
     basis: {'mol', 'wt'}, optional
         Basis of reaction. Defaults to 'mol'.
     
@@ -464,6 +483,9 @@ def conversion(self, material, T=None, P=None, phase=None, time=None):
         if config: material._imol.reset_chemicals(*config)
         return conversion
     
+    def conversion_handle(self, stream):
+        return Conversion(self)
+    
     def __call__(self, material, T=None, P=None, phase=None, time=None):
         values, config, original = as_material_array(
             material, self._basis, self._phases, self.chemicals
@@ -1724,7 +1746,29 @@ def _info(self):
 
 # %% Kinetic reactions (stoichiometric)
 
-@thermo_user
+class KineticConversion:
+    __slots__ = ('reaction', 'stream', 'index')
+    
+    def __init__(self, reaction, stream):
+        self.reaction = reaction
+        self.stream = stream
+        self.index = None
+        
+    def set_chemicals(self, chemicals):
+        self.index = self.stream.chemicals.indices([i.ID for i in chemicals])
+        
+    def __call__(self, material=None, T=None, P=None, phase=None):
+        if self.index is None:
+            with self.stream.temporary(material, T, P, phase) as stream:
+                return self.reaction.conversion(stream)
+        else:
+            flow = SparseVector.from_size(self.stream.chemicals.size)
+            flow[self.index] = material
+            with self.stream.temporary(flow, T, P, phase) as stream:
+                return self.reaction.conversion(stream)[self.index]
+
+
+@chemicals_user
 class KineticReaction:
     """
     Create an abstract KineticReaction object which defines a stoichiometric 
@@ -1741,8 +1785,9 @@ class KineticReaction:
                A dictionary of stoichiometric coefficients or a stoichiometric
                equation written as:
                i1 R1 + ... + in Rn -> j1 P1 + ... + jm Pm
-    thermo : Thermo, optional
-        Thermodynamic property package. Defaults to settings.thermo.
+    chemicals : Chemicals, optional
+        Chemicals corresponding to each entry in the stoichiometry array. 
+        Defaults to settings.chemicals.
     
     Other Parameters
     ----------------
@@ -1758,17 +1803,13 @@ class KineticReaction:
     
     Notes
     -----
-    A reaction object can react either a stream or an array. When a stream
-    is passed, it reacts either the mol or mass flow rate according to
-    the basis of the reaction object. When an array is passed, the array
-    elements are reacted regardless of what basis they are associated with.
+    A reaction object can react only a stream object (not an array).
     
     """
     rate = AbstractMethod
     volume = AbstractMethod
     phases = MaterialIndexer.phases
     __slots__ = (
-        '_stream',
         '_phases',
         '_thermo', 
         '_stoichiometry', 
@@ -1781,6 +1822,7 @@ class KineticReaction:
     istoichiometry = Reaction.istoichiometry
     stoichiometry = Reaction.stoichiometry
     reaction_chemicals = Reaction.reaction_chemicals
+    reset_chemicals = Reaction.reset_chemicals
     MWs = Reaction.MWs
     dH = Reaction.dH
     _rescale = Reaction._rescale
@@ -1794,14 +1836,12 @@ class KineticReaction:
     correct_mass_balance = Reaction.correct_mass_balance
     show = _ipython_display_ = Reaction.show
     
-    def __init__(self, reaction, thermo=None, reactant=None, *,
+    def __init__(self, reaction, chemicals=None, reactant=None, *,
                  phases=None,
                  check_mass_balance=False,
                  check_atomic_balance=False,
                  correct_atomic_balance=False):
-        thermo = self._load_thermo(thermo)
-        self._stream = tmo.Stream(thermo=thermo)
-        chemicals = thermo.chemicals
+        chemicals = self._load_chemicals(chemicals)
         if reaction:
             self._phases = phases = phase_tuple(phases) if phases else xprs.get_phases(reaction)
             if phases:
@@ -1830,61 +1870,7 @@ def __init__(self, reaction, thermo=None, reactant=None, *,
             self._stoichiometry = SparseVector.from_size(chemicals.size)
         self._rescale()
     
-    def reset_thermo(self, thermo):
-        if self._thermo is thermo: return
-        old_chemicals = self._thermo.chemicals
-        chemicals = thermo.chemicals
-        self._thermo = thermo
-        if old_chemicals is chemicals: return
-        phases = self.phases
-        stoichiometry = self._stoichiometry
-        reactant = self.reactant
-        if phases:
-            M, N = stoichiometry.shape
-            new_stoichiometry = SparseArray.from_shape([M, chemicals.size])
-            IDs = old_chemicals.IDs
-            for (i, j), k in stoichiometry.nonzero_items():
-                new_stoichiometry[i, chemicals.index(IDs[j])] = k
-            phase, reactant = reactant
-            reactant_index = phases.index(phase), chemicals.index(reactant)
-        else:
-            new_stoichiometry = SparseVector.from_size(chemicals.size)
-            IDs = old_chemicals.IDs
-            for i, j in stoichiometry.nonzero_items():
-                new_stoichiometry[chemicals.index(IDs[i])] = j
-            reactant_index = chemicals.index(reactant)
-        self._stoichiometry = new_stoichiometry
-        self._stream = tmo.Stream(thermo)
-        self._reactant_index = reactant_index
-    
-    def _get_stream(self, material, T=None, P=None, phase=None):
-        if isinstance(material, tmo.Stream):
-            try:
-                stream = self._stream
-            except:
-                if self._thermo is material._thermo:
-                    self._stream = stream = material.copy()
-                elif self.chemicals is material.chemicals:
-                    self._stream = stream = tmo.Stream(thermo=self._thermo)
-                else:
-                    raise ValueError('incompatible chemicals for kinetics')
-            stream.copy_like(material)
-            if T is not None: stream.T = T
-            if P is not None: stream.P = P
-            if phase is not None: stream.phase = phase
-        else:
-            try:
-                stream = self._stream
-            except:
-                self._stream = stream = tmo.Stream(thermo=self._thermo)
-            stream.imol.data[:] = material
-            if P is not None: stream.T = T
-            if P is not None: stream.P = P
-            stream.phase = phase
-        return stream
-    
     def conversion_bounds(self, material):
-        material, *_ = as_material_array(material, self._basis, self._phases, self.chemicals)
         stoichiometry = self._stoichiometry
         conversion = -material / stoichiometry
         bounds = (
@@ -1900,37 +1886,44 @@ def _equilibrium_objective(self, conversion, data):
         stream.set_data(data)
         return rate
     
-    def equilibrium(self, material, T, P, phase):
-        stream = self._get_stream(material, T, P, phase)
+    def equilibrium(self, stream):
         initial_condition = stream.get_data()
         f = self._equilibrium_objective
         args = (initial_condition,)
-        conversion = flx.IQ_interpolation(f, *self.conversion_bounds(stream), args=args)
+        conversion = flx.IQ_interpolation(f, *self.conversion_bounds(stream.imol.data), args=args)
         return conversion * self._stoichiometry
     
-    def conversion(self, material, T=None, P=None, phase=None, time=None):
-        values, config, original = as_material_array(
-            material, self._basis, self._phases, self.chemicals
-        )
+    def conversion_handle(self, stream):
+        return KineticConversion(self, stream)
+    
+    def _rate(self, stream):
+        conversion = self.volume(stream) * self.rate(stream)
+        mol = stream.mol
+        excess = mol + conversion
+        mask = excess < 0        
+        if excess[mask].any() and (mol[mask] > 0).all():
+            x = (-mol[mask] / conversion[mask]).min()
+            conversion *= x
+            assert 0 <= x <= 1
+        return conversion
+    
+    def conversion(self, stream, time=None):
+        if stream.chemicals is not self.chemicals: self.reset_chemicals(stream.chemicals)
         if time is None and self.volume: # Ideal continuous stirred tank CSTR
-            stream = self._get_stream(material, T, P, phase)
-            conversion = self.volume(stream) * self.rate(stream)
+            conversion = self._rate(stream)
         elif time is None: # Instantaneous reaction
-            conversion = self.equilibrium(material, T, P, phase)
+            conversion = self.equilibrium(stream)
         else: # Plug flow reaction (PFR)
             raise NotImplementedError('kinetic integration not yet implemented')
-        if config: material._imol.reset_chemicals(*config)
         return conversion
     
-    def __call__(self, material, T=None, P=None, phase=None, time=None):
-        values, config, original = as_material_array(
-            material, self._basis, self._phases, self.chemicals
-        )
+    def __call__(self, stream, time=None):
+        if stream.chemicals is not self.chemicals: self.reset_chemicals(stream.chemicals)
+        values = stream.imol.data
         if time is None and self.volume: # Ideal continuous stirred tank CSTR
-            stream = self._get_stream(material, T, P, phase)
-            values += self.volume(stream) * self.rate(stream)
+            values += self._rate(stream)
         elif time is None: # Instantaneous reaction
-            values += self.equilibrium(material, T, P, phase)
+            values += self.equilibrium(stream)
         else: # Plug flow reaction (PFR)
             raise NotImplementedError('kinetic integration not yet implemented')
         if tmo.reaction.CHECK_FEASIBILITY:
@@ -1948,8 +1941,6 @@ def __call__(self, material, T=None, P=None, phase=None, time=None):
                     values[negative_index] = 0.
         else:
             fn.remove_negligible_negative_values(values)
-        if original is not None: original[:] = values
-        if config: material._imol.reset_chemicals(*config)
     
     def to_df(self, index=None):
         columns = ['Stoichiometry (by mol)', 'Reactant']